Commit b085b387 by Manuel Schrauth

### minor tweaks; major tidy-up; simplify a number of interfaces; easier import

parent 392b3785
This diff is collapsed.
This diff is collapsed.
 from .core import HyperbolicTiling __all__ = ['HyperbolicTiling'] \ No newline at end of file
 import numpy as np import copy # imports # relative imports from .hyperpolygon import HyperPolygon from .transformation import * from .util import fund_radius, remove_duplicates, find_num_of_pgons_73 from .distance import disk_distance from .old import remove_overflow from .util import fund_radius class VertexTilingClass: # the main object of this library # essentially represents a list of polygons which build the hyperbolic lattice class HyperbolicTiling: def __init__(self, p, q, nlayers): self.p = p # number of edges (and thus number of vertices) per polygon self.q = q # number of polygons that meet at each vertex ... ... @@ -103,194 +102,14 @@ class VertexTilingClass: # assign each polygon a unique number for num, poly in enumerate(self.polygons): poly.number = num+1 # based on my own ideas # uses Moebius transformations to compute adjacent polygons by effectively mirroring the existent polygon at its edges # works perfectly for {7,3}, but is not stable for non-{7,3} tessellations class HyperbolicTiling: def __init__(self, p, q, nlayers): self.p = p self.q = q self.s = fund_radius(self.p, self.q) / 2 # np.sin(2*np.pi / self.p) * fund_radius(self.p,self.q) self.phi = 2*np.pi/self.p # acute angle of a 1/p slice of the disk self.dgts = 10 # decimal place that is rounded -> precision self.cutoff = 1 # tiling condition that tessellations end after cutoff; no cutoff for ==1 avg_d = 2*self.s*(1/2+np.sin(self.phi/2)) self.nlayers = nlayers if self.cutoff == 1 else int(2 + np.floor(disk_distance(0, self.cutoff)/avg_d)) self.nsectors = 360 self.fund_poly = self.create_fundamental_polygon() # central polygon self.lpolygons = [[] for _ in range(self.nlayers)] # for each layer there is one subarray self.polygons = [] self.centerlist = [] # used to keep track of which polygons has already been drawn def __getitem__(self, idx): return self.polygons[idx] def __iter__(self): self.iterctr = 0 self.itervar = self.polygons[self.iterctr] return self def __next__(self): if self.iterctr < len(self.polygons): retval = self.polygons[self.iterctr] self.iterctr += 1 return retval else: raise StopIteration poly.idx = num+1 def __len__(self): return len(self.polygons) def create_fundamental_polygon(self): # constructs the vertices of the fundamental hyperbolic {p,q} polygon r = fund_radius(self.p, self.q) polygon = HyperPolygon(self.p, self.q) polygon.number = 1 polygon.layer = 1 for i in range(self.p): z = complex(r * np.cos(i*self.phi), r * np.sin(2 * np.pi * i / self.p)) polygon.verticesP[i] = z polygon.verticesW[:, i] = p2w(z) return polygon def generate(self): self.lpolygons[0].append(self.fund_poly) self.centerlist.append(np.round(self.fund_poly.centerP, self.dgts)) sectors = np.arange(0, self.nsectors+1) for l in range(1, self.nlayers): for pgon in self.lpolygons[l-1]: # finds all neighbors of layer l, only adds those that... for vert_ind in range(self.p): # finds the neighbor of each edge of pgon polycopy = copy.deepcopy(pgon) adj_pgon = self.generate_adj_poly(polycopy, vert_ind) # mirror pgon at edge center = np.round(adj_pgon.centerP, self.dgts) if center not in self.centerlist and abs(center) <= self.cutoff: # ...are no duplicates adj_pgon.find_angle(360) if adj_pgon.sector in sectors: # and in the same sector adj_pgon.layer = l+1 self.centerlist.append(center) self.lpolygons[l].append(adj_pgon) pnum = 1 for lst in self.lpolygons: # flattening the list, only including the right sector polygons for polygon in lst: if polygon.sector in sectors[:-1]: polygon.number = pnum # assigning each polygon a unique number for later purposes pnum += 1 self.polygons.append(polygon) # self.find_boundary_width() # used for finding the radial distance the outmost layer covers self.angular_replicate(copy.deepcopy(self.polygons)) # uses symmetry to fill the plane by rotating the slice def generate_adj_poly(self, polygon, vert_ind): z0 = polygon.verticesP[vert_ind] polygon.moeb_origin(z0) # map vertex to origin phi = np.arctan2(polygon.verticesP[vert_ind-1].imag, polygon.verticesP[vert_ind-1].real) phi += 2 * np.pi if phi < 0 else 0 # such that phi in [0, 2pi] polygon.moeb_rotate(phi) # rotate such that edge is purely real polygon.moeb_translate(self.s) # translate such that the edge is symmetric to the im-axis polygon.mirror() # reflect at edge polygon.moeb_inverse(z0, phi, self.s) # reverses all previous transformations return polygon def angular_replicate(self, polygons): # tessellates the disk by applying a rotation of 2pi/p to the pizza slice polygons.pop(0) # first pgon is in each 1/p sector for p in range(1, self.p): for ind, polygon in enumerate(polygons): pgon = copy.deepcopy(polygon) pgon.rotate(-p*self.phi) # - stems from definition of moeb_rotate() pgon.find_angle(360) self.polygons.append(pgon) for num, poly in enumerate(self.polygons): poly.number = num+1 # the following are experimental, untested, and currently not being used def generate_alt(self): p_ind = 0 pcounter = 2 self.draw_neighbors(self.fund_poly, self.nlayers-1) print(len(self.polygons)) print("soll = ", (find_num_of_pgons_73(self.nlayers) - 1) / 7 + 1) lim = 1+(find_num_of_pgons_73(self.nlayers)-1)/7 while pcounter < 1+(find_num_of_pgons_73(self.nlayers)-1)/7: for p_ind in range(0, int(lim)): # while pcounter < find_num_of_pgons73(self.nlayers): for vert_ind in range(self.p): polycopy = copy.deepcopy(self.polygons[p_ind]) polycopy = self.generate_adj_poly(polycopy, vert_ind) center = np.round(polycopy.centerP, self.dgts) if center not in self.centerlist: polycopy.find_angle(360) if polycopy.sector == 0: self.polygons.append(polycopy) self.centerlist.append(center) polycopy.number = pcounter pcounter += 1 p_ind += 1 self.polygons = remove_duplicates(self.polygons) print(len(self.polygons)) pgons_to_remove = (find_num_of_pgons_73(self.nlayers)-1)/7+1 print("to remove: ", len(self.polygons) - pgons_to_remove) self.polygons = remove_overflow(self.polygons, num=len(self.polygons) - pgons_to_remove) self.polygons.pop(0) # exclude first element to avoid duplicates polygons = copy.deepcopy(self.polygons) for p in range(0): self.angular_replicate(polygons) # self.polygons = remove_duplicates(self.polygons) self.polygons.insert(0, self.fund_poly) # reinsert first polygon def generate_recursive(self): if self.nlayers == 1: return self.draw_neighbors(self.fund_poly, self.nlayers - 1) def draw_neighbors(self, polygon, ltd): # finish n of 1 poly first, then move on polygon.number = ltd self.polygons.append(polygon) self.centerlist.append(np.round(polygon.centerP, self.dgts)) if ltd > 0: # inefficient neighbors = [] for v in range(self.p): pgon = copy.deepcopy(polygon) poly = self.generate_adj_poly(pgon, v) # print(abs(polygon.centerP)) center = np.round(poly.centerP, self.dgts) if center not in self.centerlist: poly.find_angle(360) if poly.sector == 0: poly.number = ltd-1 self.polygons.append(poly) self.centerlist.append(center) # neighbors.append(poly) self.draw_neighbors(poly, ltd-1) # for _ in neighbors: # self.draw_neighbors(_, ltd-1) def find_boundary_width(self): rmax, rmin = 0, disk_distance(0, self.polygons[1].centerP) for poly in self.polygons: if poly.layer == self.nlayers: d = disk_distance(0, poly.centerP) rmax = d if d>rmax else rmax rmin = d if d
 ... ... @@ -7,8 +7,8 @@ def weierstrass_distance(a, b): if arg < 1: return 0 else: return math.acosh(arg) # math.acosh is faster for scalars #return np.arccosh(arg) return math.acosh(arg) # for scalars math.acosh is usually faster than np.arccosh def disk_distance(z1, z2): ... ...
 import numpy as np import copy # relative imports from .hyperpolygon import HyperPolygon from .transformation import * from .util import fund_radius, remove_duplicates, find_num_of_pgons_73 from .distance import disk_distance from .old import remove_overflow # Tessellation algorithm written by F. Dusel # uses Moebius transformations to compute adjacent polygons by effectively mirroring the existent polygon at its edges # works perfectly for {7,3}, but is not stable for non-{7,3} tessellations class HyperbolicTiling73Experimental: def __init__(self, p, q, nlayers): self.p = p self.q = q self.s = fund_radius(self.p, self.q) / 2 # np.sin(2*np.pi / self.p) * fund_radius(self.p,self.q) self.phi = 2*np.pi/self.p # acute angle of a 1/p slice of the disk self.dgts = 10 # decimal place that is rounded -> precision self.cutoff = 1 # tiling condition that tessellations end after cutoff; no cutoff for ==1 avg_d = 2*self.s*(1/2+np.sin(self.phi/2)) self.nlayers = nlayers if self.cutoff == 1 else int(2 + np.floor(disk_distance(0, self.cutoff)/avg_d)) self.nsectors = 360 self.fund_poly = self.create_fundamental_polygon() # central polygon self.lpolygons = [[] for _ in range(self.nlayers)] # for each layer there is one subarray self.polygons = [] self.centerlist = [] # used to keep track of which polygons has already been drawn def __getitem__(self, idx): return self.polygons[idx] def __iter__(self): self.iterctr = 0 self.itervar = self.polygons[self.iterctr] return self def __next__(self): if self.iterctr < len(self.polygons): retval = self.polygons[self.iterctr] self.iterctr += 1 return retval else: raise StopIteration def __len__(self): return len(self.polygons) def create_fundamental_polygon(self): # constructs the vertices of the fundamental hyperbolic {p,q} polygon r = fund_radius(self.p, self.q) polygon = HyperPolygon(self.p, self.q) polygon.number = 1 polygon.layer = 1 for i in range(self.p): z = complex(r * np.cos(i*self.phi), r * np.sin(2 * np.pi * i / self.p)) polygon.verticesP[i] = z polygon.verticesW[:, i] = p2w(z) return polygon def generate(self): self.lpolygons[0].append(self.fund_poly) self.centerlist.append(np.round(self.fund_poly.centerP, self.dgts)) sectors = np.arange(0, self.nsectors+1) for l in range(1, self.nlayers): for pgon in self.lpolygons[l-1]: # finds all neighbors of layer l, only adds those that... for vert_ind in range(self.p): # finds the neighbor of each edge of pgon polycopy = copy.deepcopy(pgon) adj_pgon = self.generate_adj_poly(polycopy, vert_ind) # mirror pgon at edge center = np.round(adj_pgon.centerP, self.dgts) if center not in self.centerlist and abs(center) <= self.cutoff: # ...are no duplicates adj_pgon.find_angle(360) if adj_pgon.sector in sectors: # and in the same sector adj_pgon.layer = l+1 self.centerlist.append(center) self.lpolygons[l].append(adj_pgon) pnum = 1 for lst in self.lpolygons: # flattening the list, only including the right sector polygons for polygon in lst: if polygon.sector in sectors[:-1]: polygon.number = pnum # assigning each polygon a unique number for later purposes pnum += 1 self.polygons.append(polygon) # self.find_boundary_width() # used for finding the radial distance the outmost layer covers self.angular_replicate(copy.deepcopy(self.polygons)) # uses symmetry to fill the plane by rotating the slice def generate_adj_poly(self, polygon, vert_ind): z0 = polygon.verticesP[vert_ind] polygon.moeb_origin(z0) # map vertex to origin phi = np.arctan2(polygon.verticesP[vert_ind-1].imag, polygon.verticesP[vert_ind-1].real) phi += 2 * np.pi if phi < 0 else 0 # such that phi in [0, 2pi] polygon.moeb_rotate(phi) # rotate such that edge is purely real polygon.moeb_translate(self.s) # translate such that the edge is symmetric to the im-axis polygon.mirror() # reflect at edge polygon.moeb_inverse(z0, phi, self.s) # reverses all previous transformations return polygon def angular_replicate(self, polygons): # tessellates the disk by applying a rotation of 2pi/p to the pizza slice polygons.pop(0) # first pgon is in each 1/p sector for p in range(1, self.p): for ind, polygon in enumerate(polygons): pgon = copy.deepcopy(polygon) pgon.rotate(-p*self.phi) # - stems from definition of moeb_rotate() pgon.find_angle(360) self.polygons.append(pgon) for num, poly in enumerate(self.polygons): poly.number = num+1 # the following are experimental, untested, and currently not being used def generate_alt(self): p_ind = 0 pcounter = 2 self.draw_neighbors(self.fund_poly, self.nlayers-1) print(len(self.polygons)) print("soll = ", (find_num_of_pgons_73(self.nlayers) - 1) / 7 + 1) lim = 1+(find_num_of_pgons_73(self.nlayers)-1)/7 while pcounter < 1+(find_num_of_pgons_73(self.nlayers)-1)/7: for p_ind in range(0, int(lim)): # while pcounter < find_num_of_pgons73(self.nlayers): for vert_ind in range(self.p): polycopy = copy.deepcopy(self.polygons[p_ind]) polycopy = self.generate_adj_poly(polycopy, vert_ind) center = np.round(polycopy.centerP, self.dgts) if center not in self.centerlist: polycopy.find_angle(360) if polycopy.sector == 0: self.polygons.append(polycopy) self.centerlist.append(center) polycopy.number = pcounter pcounter += 1 p_ind += 1 self.polygons = remove_duplicates(self.polygons) print(len(self.polygons)) pgons_to_remove = (find_num_of_pgons_73(self.nlayers)-1)/7+1 print("to remove: ", len(self.polygons) - pgons_to_remove) self.polygons = remove_overflow(self.polygons, num=len(self.polygons) - pgons_to_remove) self.polygons.pop(0) # exclude first element to avoid duplicates polygons = copy.deepcopy(self.polygons) for p in range(0): self.angular_replicate(polygons) # self.polygons = remove_duplicates(self.polygons) self.polygons.insert(0, self.fund_poly) # reinsert first polygon def generate_recursive(self): if self.nlayers == 1: return self.draw_neighbors(self.fund_poly, self.nlayers - 1) def draw_neighbors(self, polygon, ltd): # finish n of 1 poly first, then move on polygon.number = ltd self.polygons.append(polygon) self.centerlist.append(np.round(polygon.centerP, self.dgts)) if ltd > 0: # inefficient neighbors = [] for v in range(self.p): pgon = copy.deepcopy(polygon) poly = self.generate_adj_poly(pgon, v) # print(abs(polygon.centerP)) center = np.round(poly.centerP, self.dgts) if center not in self.centerlist: poly.find_angle(360) if poly.sector == 0: poly.number = ltd-1 self.polygons.append(poly) self.centerlist.append(center) # neighbors.append(poly) self.draw_neighbors(poly, ltd-1) # for _ in neighbors: # self.draw_neighbors(_, ltd-1) def find_boundary_width(self): rmax, rmin = 0, disk_distance(0, self.polygons[1].centerP) for poly in self.polygons: if poly.layer == self.nlayers: d = disk_distance(0, poly.centerP) rmax = d if d>rmax else rmax rmin = d if d
 from math import floor from .transformation import * # defines a hyperbolic polygon class HyperPolygon: def __init__(self, p, q): self.p = p # no. of edges self.q = q # no. of adjacent polygons per vertex self.p = p # number of edges self.q = q # number of adjacent polygons per vertex # Poincare disk coordinates self.centerP = complex(0, 0) # center self.verticesP = np.zeros(shape=self.p, dtype=np.complex128) # vertices self.centerP = complex(0, 0) # poincare coordinates of the center self.centerW = np.array([0, 0, 1]) # Weierstrass (hyperboloid) coordinates self.centerW = np.array([0, 0, 1]) # center self.verticesW = np.zeros((3, self.p)) # vertices self.verticesP = np.zeros(shape=self.p, dtype=np.complex128) # POINCARE Disk 2D self.verticesW = np.zeros((3, self.p)) # WEIERSTRASS 3D self.idx = 1 # auxiliary scalar index; can be used, e.g, for easy identifaction inside a tessellation self.layer = 1 # encodes in which layer of a tessellation this polygons is located self.sector = 0 # index of the sector this polygons is located; can be used for finding neighbours more efficiently self.angle = 0 # angle between self.centerP and the positive x-axis self.val = 0 # assign a value (useful in any application) self.number = 1 # for counting them and coloring them successively self.layer = 1 self.angle = 0 # angle between self.centerP and the positive x-axis self.sector = 0 # sector for finding NN more efficiently self.val = 0 # field value for Ising/PDEQ solver def __eq__(self, other): # checks if two polygons are equivalent; currently not being used # checks whether two polygons are equal (within given numerical precision) # untested! def __eq__(self, other, digits=7): re1 = self.centerP.real im1 = self.centerP.imag re2 = other.centerP.real im2 = other.centerP.imag c1 = complex(round(re1, 5), round(im1, 5)) c2 = complex(round(re2, 5), round(im2, 5)) print("__ eq__ was used") # in case it happens by accident c1 = complex(round(re1, digits), round(im1, digits)) c2 = complex(round(re2, digits), round(im2, digits)) print("Warning: Equality operator of the HyperPolygon class is untested!") return c1 == c2 # transforms all points of the polygon by matrix tmat; only used by dunhams tiling class # transforms all points of the polygon by matrix tmat # only used by HyperbolicTilingDunham def transform(self, tmat): self.centerW = tmat @ self.centerW self.centerP = w2p(self.centerW) ... ... @@ -38,7 +47,9 @@ class HyperPolygon: self.verticesW[:, i] = tmat @ self.verticesW[:, i] self.verticesP[i] = w2p(self.verticesW[:, i]) def moeb_origin(self, z0): # transforms the whole polygon such that z0 is mapped to origin # transforms the entire polygon such that z0 is mapped to origin def moeb_origin(self, z0): self.centerP = moeb_origin_trafo(z0, self.centerP) self.centerW = p2w(self.centerP) # self.find_angle(360) # this might be superfluous ... ... @@ -47,33 +58,34 @@ class HyperPolygon: self.verticesP[i] = z self.verticesW[:, i] = p2w(self.verticesP[i]) def moeb_rotate(self, phi): # rotates each point of the polygon by phi self.centerP = moeb_rotate_trafo(self.centerP, -phi) # these two lines might be redundant self.centerW = p2w(self.centerP) # self.find_angle(360) # i think this is not needed anymore for i in range(self.p): z = moeb_rotate_trafo(self.verticesP[i], -phi) self.verticesP[i] = z self.verticesW[:, i] = p2w(self.verticesP[i]) def moeb_translate(self, s): self.centerP = moeb_translate_trafo(self.centerP, s) self.centerW = p2w(self.centerP) # self.find_angle(360) # i think this is not needed anymore for i in range(self.p): z = moeb_translate_trafo(self.verticesP[i], s) self.verticesP[i] = z self.verticesW[:, i] = p2w(self.verticesP[i]) def moeb_inverse(self, z0, phi=0, s=0): self.centerP = moeb_inverse_trafo(self.centerP, z0, -phi, s) self.centerW = p2w(self.centerP) # self.find_angle(360) # superfluous for i in range(self.p): z = moeb_inverse_trafo(self.verticesP[i], z0, -phi, s) self.verticesP[i] = z self.verticesW[:, i] = p2w(self.verticesP[i]) def rotate(self, phi): rotation = np.exp(complex(0, phi)) z = self.centerP ... ... @@ -86,6 +98,7 @@ class HyperPolygon: self.verticesP[i] = z self.verticesW[:, i] = p2w(self.verticesP[i]) def find_edges(self): # finds twice the amount of necessary edges! xedges = [] yedges = [] ... ... @@ -101,12 +114,14 @@ class HyperPolygon: return [xedges, yedges] def find_angle(self, k, offset=0): # has to be called after the inverse trafo... self.angle = np.arctan2(self.centerP.imag, self.centerP.real)*180/np.pi-offset # requires y first self.angle = np.round(self.angle, 10) self.angle += 360 if self.angle < 0 else 0 self.sector = floor(self.angle/(360/(k*self.p))) # k*p sectors; insert offset +0.1 here?! def mirror(self): self.centerP = complex(self.centerP.real, - self.centerP.imag) # mirror on axis Im(z)=0 self.centerW = p2w(self.centerP) ... ...
 import numpy as np from .distance import weierstrass_distance import time # wrapper to provide a nicer interface def find(tiling, nn_dist, which="optimized_slice"): if which == "optimized_slice": return find_nn_optimized_slice(tiling, nn_dist) # fastest elif which == "optimized": return find_nn_optimized(tiling, nn_dist) elif which == "brute_force":