crystal.polytope
1# https://mathsfromnothing.au/checking-if-a-point-is-inside-a-shape/ 2# Checking if a Point is Inside a Polytope 3 4# TODO mv BACKEND to config 5BACKEND='numpy' 6# BACKEND='torch' 7if BACKEND=='numpy': 8 from numpy import array, shape, ones, moveaxis, einsum, min, max, vstack, arange, meshgrid, linspace, matmul, concatenate 9 from numpy.linalg import det, inv 10 from numpy.random import rand 11elif BACKEND=='torch': 12 from torch import ones, moveaxis, einsum, vstack, det, rand, min, max, arange, meshgrid, concatenate 13 from torch import tensor as array 14import matplotlib.pyplot as plt 15from mpl_toolkits.mplot3d.art3d import Poly3DCollection 16from matplotlib.patches import Path, PathPatch 17from scipy.spatial import ConvexHull, convex_hull_plot_2d 18from scipy.spatial import Delaunay 19 20# def external_prod(A): 21# """Calculates the external product of ... x n-1 column vectors of length n. A generalization of the cross product. 22 23# Examples 24# ======== 25# >>> A=[[1,0,0],[0,1,0]] 26# >>> B=[[0,1,0],[1,0,0]] 27# >>> C=[[0,0,1],[1,0,0]] 28# >>> external_prod(A) 29# array([0., 0., 1.]) 30# >>> external_prod(B) 31# array([ 0., 0., -1.]) 32# >>> external_prod([A,B,C]) 33# array([[ 0., 0., 1.], 34# [ 0., 0., -1.], 35# [ 0., 1., 0.]]) 36# """ 37# A = array(A) 38# dim = shape(A) 39# (n,d) = dim[-2:] 40# if n!=d-1: 41# raise ValueError("A must be an (...,n-1,n) array!") 42# sgn = ones(d,int) 43# sgn[1::d] = -1 44# _det = det([A[...,[i for i in range(d) if i!=j]] for j in range(d)]) 45# _det = moveaxis(_det, -1, 0) 46# return einsum('i,...i->...i', sgn, _det) 47 48# def normal(face): 49# """ 50# TODO generalise to list of faces 51# """ 52# face = array(face) 53# edges = face[...,1:,:] - face[...,0,:] 54# return external_prod(edges) 55 56class Polytope(Delaunay): 57 """ 58 Polytope 59 60 Attributes 61 ---------- 62 points: boundary points in Cartesian coordinates 63 color: colour of each faces (optional) 64 n_points: number of points 65 n_faces: number of faces 66 dim: the dimension of the polytope 67 faces_cartesian_coordinates: faces in Caresian coordinates 68 plot: plot the Polytope 69 70 Examples 71 -------- 72 73 >>> z = 2**(1/2) 74 >>> points = [ [1,1, 1], [-1,1, 1], [-1,-1, 1], [1,-1, 1], 75 ... [z,0,-1], [ 0,z,-1], [-z, 0,-1], [0,-z,-1] ] # 8 points 76 >>> faces_a = [ [0,3,4], [1,0,5], [2,1,6], [3,2,7] ] 77 >>> faces_b = [ [4,5,0], [5,6,1], [6,7,2], [7,4,3] ] 78 >>> faces_c = [ [0,1,2], [4,5,6] ] 79 >>> faces_d = [ [2,3,0], [6,7,4] ] 80 >>> faces = faces_a + faces_b + faces_c + faces_d 81 >>> colors = ['teal']*4 + ['indianred']*4 + ['olivedrab']*2 + ['mediumpurple']*2 # 12 colors 82 >>> polytope = Polytope(points) 83 84 """ 85 def __init__(self, points:list, basis=None): 86 self.dim = len(points[0]) 87 self.color = "teal" 88 89 super().__init__(points) 90 91 # construct faces using Convex Hull (for visualisation purposes): 92 hull = ConvexHull(points) 93 indices = hull.simplices 94 self.faces = self.points[indices] 95 96 self.pts = [] 97 98 if type(basis)==type(None): 99 basis = array([[0 for _ in range(self.nd)] for _ in range(self.nd)]) 100 for d in range(self.nd): 101 basis[d,d] = 1 102 self.basis = array(basis) 103 self.inv_basis = inv(basis) 104 105 def __repr__(self): 106 """Returns object class name""" 107 return f'Polytope with {self.n_points} points' 108 109 @property 110 def n_points(self): 111 return len(self.points) 112 113 @property 114 def n_faces(self): 115 return len(self.faces) 116 117 @property 118 def nd(self): 119 return len(self.points[0]) 120 121 @property 122 def faces_cartesian_coordinates(self): 123 return [ [ self.points[f] for f in face] for face in self.faces] 124 125 def _get_boundary(self, points): 126 """d x 2 points describing d-dimensional minimum and maximum points 127 of a hypercube wrapping the points.""" 128 if BACKEND=='numpy': 129 mn = min(points, axis=(0)) 130 mx = max(points, axis=(0)) 131 return array([mn,mx]).T 132 elif BACKEND=='torch': 133 mn = (min(array(points), axis=(0)).values.unsqueeze(0)).tolist() 134 mx = (max(array(points), axis=(0)).values.unsqueeze(0)).tolist() 135 return array([mn,mx]).T 136 else: 137 raise ValueError(f'BACKEND={BACKEND} unsupported!') 138 139 @property 140 def corners(self): 141 """d x 2 points describing d-dimensional minimum and maximum points 142 of a hypercube wrapping the polytope.""" 143 return self._get_boundary(points=self.points) 144 145 @property 146 def frac_corners(self): 147 """d x 2 points describing d-dimensional minimum and maximum points 148 of a hypercube wrapping the polytope.""" 149 points = matmul(self.points, self.inv_basis.T) 150 return self._get_boundary(points=points) 151 152 def _scaled_lim(self): 153 try: 154 points = concatenate((self.points, self.pts), axis=0) 155 except: 156 points = self.points 157 158 lim = [] 159 for l in self._get_boundary(points=points): 160 scale = (l[1] - l[0])/50 161 l = [l[0]-scale, l[1]+scale] 162 lim.append(l) 163 return lim 164 165 def set_centre(self, centre): 166 points = [pt + centre for pt in self.points] 167 self.__init__(points=points, basis=self.basis) 168 self.centre = centre 169 170 def generate_lattice(self): 171 172 frac_corners = self.frac_corners 173 174 frac_points = [] 175 for i in range(self.dim): 176 frac_points.append(arange(frac_corners[i,0],frac_corners[i,1]+1,1)) 177 frac_points = vstack(meshgrid(*frac_points)).reshape(self.dim,-1).T 178 frac_points = array(frac_points).T 179 180 cartesian_points = matmul(self.basis, frac_points).T 181 182 return cartesian_points 183 184 def cut(self, points:list=[]): 185 """ 186 Divide the points into two subsets --inside and outside the boundary 187 defined by the polytope. 188 """ 189 if len(points) == 0: 190 points = self.generate_lattice() 191 points = array(points) 192 index = self.find_simplex(points) 193 inside = points[index > -1] 194 outside = points[index == -1] 195 self.pts.extend(inside) 196 self.pts.extend(outside) 197 return inside, outside 198 199 def plot(self, ax=None): 200 if self.dim==3: 201 if ax==None: 202 ax = plt.axes(projection='3d') 203 plt.tight_layout() 204 for f in self.faces: 205 face = Poly3DCollection([f]) 206 face.set_alpha(0.5) 207 face.set_edgecolor('k') 208 face.set_color(self.color) 209 ax.add_collection3d(face) 210 lim = self._scaled_lim() 211 ax.set(xlim=lim[0], ylim=lim[1], zlim=lim[2]) 212 ax.set_xlabel('x') 213 ax.set_ylabel('y') 214 ax.set_zlabel('z') 215 elif self.dim==2: 216 if ax==None: 217 ax = plt.axes() 218 plt.tight_layout() 219 points = array(self.points) 220 surface = vstack([points, points[0]]) 221 ax.add_patch(PathPatch(Path(surface), color=self.color, alpha=0.5)) 222 lim = self._scaled_lim() 223 ax.set(xlim=lim[0], ylim=lim[1]) 224 ax.set_xlabel('x') 225 ax.set_ylabel('y') 226 else: 227 raise ValueError("Dim!=2,3 not yet configured") 228 return ax 229 230def random(): 231 232 polytope = Polytope(points=[[-10, -10], [-10, 10], [10, 10], [10, -10]]) 233 234 points = 10*(rand(100, 2)*4 - 2) 235 inside, outside = polytope.cut(points) 236 237 ax = polytope.plot() 238 ax.scatter(*inside.T, c='g') 239 ax.scatter(*outside.T, c='r') 240 241 plt.show() 242 243def lines2d(): 244 245 polytope = Polytope(points=[[-10, -10], [-10, 10], [10, 10], [10, -10]]) 246 247 points = 10*(rand(100, 2)*4 - 2) 248 points = [[i*0.1,5] for i in range(-110,110)] 249 points.extend([[i*0.1,i*0.2] for i in range(-110,110)]) 250 inside, outside = polytope.cut(points) 251 252 ax = polytope.plot() 253 ax.scatter(*inside.T, c='g') 254 ax.scatter(*outside.T, c='r') 255 256 plt.show() 257 258# Examples 259 260def lattice(): 261 262 polytope = Polytope(points=[[5, 5], [10, 5], [10, 10]], basis = [[2,0.25],[-0.25,1]]) 263 264 inside, outside = polytope.cut() 265 266 ax = polytope.plot() 267 ax.scatter(*inside.T, c='g') 268 ax.scatter(*outside.T, c='r') 269 270 plt.show() 271 272def plot_3d(): 273 274 z = 5*2**(1/2) 275 points = [ [5,5, 5], [-5,5, 5], [-5,-5, 5], [5,-5, 5], 276 [z,0,-5], [ 0,z,-5], [-z, 0,-5], [0,-z,-5] ] # 8 points 277 faces_a = [ [0,3,4], [1,0,5], [2,1,6], [3,2,7] ] 278 faces_b = [ [4,5,0], [5,6,1], [6,7,2], [7,4,3] ] 279 faces_c = [ [0,1,2], [4,5,6] ] 280 faces_d = [ [2,3,0], [6,7,4] ] 281 faces = faces_a + faces_b + faces_c + faces_d 282 colors = ['teal']*4 + ['indianred']*4 + ['olivedrab']*2 + ['mediumpurple']*2 # 12 colors 283 polytope = Polytope(points) 284 inside, outside = polytope.cut() 285 ax = polytope.plot() 286 ax.scatter(*inside.T, c='g') 287 ax.scatter(*outside.T, c='r') 288 plt.show() 289 290# random() 291# lattice() 292# lines2d() 293# plot_3d()
BACKEND =
'numpy'
class
Polytope(scipy.spatial._qhull.Delaunay):
57class Polytope(Delaunay): 58 """ 59 Polytope 60 61 Attributes 62 ---------- 63 points: boundary points in Cartesian coordinates 64 color: colour of each faces (optional) 65 n_points: number of points 66 n_faces: number of faces 67 dim: the dimension of the polytope 68 faces_cartesian_coordinates: faces in Caresian coordinates 69 plot: plot the Polytope 70 71 Examples 72 -------- 73 74 >>> z = 2**(1/2) 75 >>> points = [ [1,1, 1], [-1,1, 1], [-1,-1, 1], [1,-1, 1], 76 ... [z,0,-1], [ 0,z,-1], [-z, 0,-1], [0,-z,-1] ] # 8 points 77 >>> faces_a = [ [0,3,4], [1,0,5], [2,1,6], [3,2,7] ] 78 >>> faces_b = [ [4,5,0], [5,6,1], [6,7,2], [7,4,3] ] 79 >>> faces_c = [ [0,1,2], [4,5,6] ] 80 >>> faces_d = [ [2,3,0], [6,7,4] ] 81 >>> faces = faces_a + faces_b + faces_c + faces_d 82 >>> colors = ['teal']*4 + ['indianred']*4 + ['olivedrab']*2 + ['mediumpurple']*2 # 12 colors 83 >>> polytope = Polytope(points) 84 85 """ 86 def __init__(self, points:list, basis=None): 87 self.dim = len(points[0]) 88 self.color = "teal" 89 90 super().__init__(points) 91 92 # construct faces using Convex Hull (for visualisation purposes): 93 hull = ConvexHull(points) 94 indices = hull.simplices 95 self.faces = self.points[indices] 96 97 self.pts = [] 98 99 if type(basis)==type(None): 100 basis = array([[0 for _ in range(self.nd)] for _ in range(self.nd)]) 101 for d in range(self.nd): 102 basis[d,d] = 1 103 self.basis = array(basis) 104 self.inv_basis = inv(basis) 105 106 def __repr__(self): 107 """Returns object class name""" 108 return f'Polytope with {self.n_points} points' 109 110 @property 111 def n_points(self): 112 return len(self.points) 113 114 @property 115 def n_faces(self): 116 return len(self.faces) 117 118 @property 119 def nd(self): 120 return len(self.points[0]) 121 122 @property 123 def faces_cartesian_coordinates(self): 124 return [ [ self.points[f] for f in face] for face in self.faces] 125 126 def _get_boundary(self, points): 127 """d x 2 points describing d-dimensional minimum and maximum points 128 of a hypercube wrapping the points.""" 129 if BACKEND=='numpy': 130 mn = min(points, axis=(0)) 131 mx = max(points, axis=(0)) 132 return array([mn,mx]).T 133 elif BACKEND=='torch': 134 mn = (min(array(points), axis=(0)).values.unsqueeze(0)).tolist() 135 mx = (max(array(points), axis=(0)).values.unsqueeze(0)).tolist() 136 return array([mn,mx]).T 137 else: 138 raise ValueError(f'BACKEND={BACKEND} unsupported!') 139 140 @property 141 def corners(self): 142 """d x 2 points describing d-dimensional minimum and maximum points 143 of a hypercube wrapping the polytope.""" 144 return self._get_boundary(points=self.points) 145 146 @property 147 def frac_corners(self): 148 """d x 2 points describing d-dimensional minimum and maximum points 149 of a hypercube wrapping the polytope.""" 150 points = matmul(self.points, self.inv_basis.T) 151 return self._get_boundary(points=points) 152 153 def _scaled_lim(self): 154 try: 155 points = concatenate((self.points, self.pts), axis=0) 156 except: 157 points = self.points 158 159 lim = [] 160 for l in self._get_boundary(points=points): 161 scale = (l[1] - l[0])/50 162 l = [l[0]-scale, l[1]+scale] 163 lim.append(l) 164 return lim 165 166 def set_centre(self, centre): 167 points = [pt + centre for pt in self.points] 168 self.__init__(points=points, basis=self.basis) 169 self.centre = centre 170 171 def generate_lattice(self): 172 173 frac_corners = self.frac_corners 174 175 frac_points = [] 176 for i in range(self.dim): 177 frac_points.append(arange(frac_corners[i,0],frac_corners[i,1]+1,1)) 178 frac_points = vstack(meshgrid(*frac_points)).reshape(self.dim,-1).T 179 frac_points = array(frac_points).T 180 181 cartesian_points = matmul(self.basis, frac_points).T 182 183 return cartesian_points 184 185 def cut(self, points:list=[]): 186 """ 187 Divide the points into two subsets --inside and outside the boundary 188 defined by the polytope. 189 """ 190 if len(points) == 0: 191 points = self.generate_lattice() 192 points = array(points) 193 index = self.find_simplex(points) 194 inside = points[index > -1] 195 outside = points[index == -1] 196 self.pts.extend(inside) 197 self.pts.extend(outside) 198 return inside, outside 199 200 def plot(self, ax=None): 201 if self.dim==3: 202 if ax==None: 203 ax = plt.axes(projection='3d') 204 plt.tight_layout() 205 for f in self.faces: 206 face = Poly3DCollection([f]) 207 face.set_alpha(0.5) 208 face.set_edgecolor('k') 209 face.set_color(self.color) 210 ax.add_collection3d(face) 211 lim = self._scaled_lim() 212 ax.set(xlim=lim[0], ylim=lim[1], zlim=lim[2]) 213 ax.set_xlabel('x') 214 ax.set_ylabel('y') 215 ax.set_zlabel('z') 216 elif self.dim==2: 217 if ax==None: 218 ax = plt.axes() 219 plt.tight_layout() 220 points = array(self.points) 221 surface = vstack([points, points[0]]) 222 ax.add_patch(PathPatch(Path(surface), color=self.color, alpha=0.5)) 223 lim = self._scaled_lim() 224 ax.set(xlim=lim[0], ylim=lim[1]) 225 ax.set_xlabel('x') 226 ax.set_ylabel('y') 227 else: 228 raise ValueError("Dim!=2,3 not yet configured") 229 return ax
Polytope
Attributes
points: boundary points in Cartesian coordinates color: colour of each faces (optional) n_points: number of points n_faces: number of faces dim: the dimension of the polytope faces_cartesian_coordinates: faces in Caresian coordinates plot: plot the Polytope
Examples
>>> z = 2**(1/2)
>>> points = [ [1,1, 1], [-1,1, 1], [-1,-1, 1], [1,-1, 1],
... [z,0,-1], [ 0,z,-1], [-z, 0,-1], [0,-z,-1] ] # 8 points
>>> faces_a = [ [0,3,4], [1,0,5], [2,1,6], [3,2,7] ]
>>> faces_b = [ [4,5,0], [5,6,1], [6,7,2], [7,4,3] ]
>>> faces_c = [ [0,1,2], [4,5,6] ]
>>> faces_d = [ [2,3,0], [6,7,4] ]
>>> faces = faces_a + faces_b + faces_c + faces_d
>>> colors = ['teal']*4 + ['indianred']*4 + ['olivedrab']*2 + ['mediumpurple']*2 # 12 colors
>>> polytope = Polytope(points)
Polytope(points: list, basis=None)
86 def __init__(self, points:list, basis=None): 87 self.dim = len(points[0]) 88 self.color = "teal" 89 90 super().__init__(points) 91 92 # construct faces using Convex Hull (for visualisation purposes): 93 hull = ConvexHull(points) 94 indices = hull.simplices 95 self.faces = self.points[indices] 96 97 self.pts = [] 98 99 if type(basis)==type(None): 100 basis = array([[0 for _ in range(self.nd)] for _ in range(self.nd)]) 101 for d in range(self.nd): 102 basis[d,d] = 1 103 self.basis = array(basis) 104 self.inv_basis = inv(basis)
corners
140 @property 141 def corners(self): 142 """d x 2 points describing d-dimensional minimum and maximum points 143 of a hypercube wrapping the polytope.""" 144 return self._get_boundary(points=self.points)
d x 2 points describing d-dimensional minimum and maximum points of a hypercube wrapping the polytope.
frac_corners
146 @property 147 def frac_corners(self): 148 """d x 2 points describing d-dimensional minimum and maximum points 149 of a hypercube wrapping the polytope.""" 150 points = matmul(self.points, self.inv_basis.T) 151 return self._get_boundary(points=points)
d x 2 points describing d-dimensional minimum and maximum points of a hypercube wrapping the polytope.
def
generate_lattice(self):
171 def generate_lattice(self): 172 173 frac_corners = self.frac_corners 174 175 frac_points = [] 176 for i in range(self.dim): 177 frac_points.append(arange(frac_corners[i,0],frac_corners[i,1]+1,1)) 178 frac_points = vstack(meshgrid(*frac_points)).reshape(self.dim,-1).T 179 frac_points = array(frac_points).T 180 181 cartesian_points = matmul(self.basis, frac_points).T 182 183 return cartesian_points
def
cut(self, points: list = []):
185 def cut(self, points:list=[]): 186 """ 187 Divide the points into two subsets --inside and outside the boundary 188 defined by the polytope. 189 """ 190 if len(points) == 0: 191 points = self.generate_lattice() 192 points = array(points) 193 index = self.find_simplex(points) 194 inside = points[index > -1] 195 outside = points[index == -1] 196 self.pts.extend(inside) 197 self.pts.extend(outside) 198 return inside, outside
Divide the points into two subsets --inside and outside the boundary defined by the polytope.
def
plot(self, ax=None):
200 def plot(self, ax=None): 201 if self.dim==3: 202 if ax==None: 203 ax = plt.axes(projection='3d') 204 plt.tight_layout() 205 for f in self.faces: 206 face = Poly3DCollection([f]) 207 face.set_alpha(0.5) 208 face.set_edgecolor('k') 209 face.set_color(self.color) 210 ax.add_collection3d(face) 211 lim = self._scaled_lim() 212 ax.set(xlim=lim[0], ylim=lim[1], zlim=lim[2]) 213 ax.set_xlabel('x') 214 ax.set_ylabel('y') 215 ax.set_zlabel('z') 216 elif self.dim==2: 217 if ax==None: 218 ax = plt.axes() 219 plt.tight_layout() 220 points = array(self.points) 221 surface = vstack([points, points[0]]) 222 ax.add_patch(PathPatch(Path(surface), color=self.color, alpha=0.5)) 223 lim = self._scaled_lim() 224 ax.set(xlim=lim[0], ylim=lim[1]) 225 ax.set_xlabel('x') 226 ax.set_ylabel('y') 227 else: 228 raise ValueError("Dim!=2,3 not yet configured") 229 return ax
def
random():
def
lines2d():
244def lines2d(): 245 246 polytope = Polytope(points=[[-10, -10], [-10, 10], [10, 10], [10, -10]]) 247 248 points = 10*(rand(100, 2)*4 - 2) 249 points = [[i*0.1,5] for i in range(-110,110)] 250 points.extend([[i*0.1,i*0.2] for i in range(-110,110)]) 251 inside, outside = polytope.cut(points) 252 253 ax = polytope.plot() 254 ax.scatter(*inside.T, c='g') 255 ax.scatter(*outside.T, c='r') 256 257 plt.show()
def
lattice():
def
plot_3d():
273def plot_3d(): 274 275 z = 5*2**(1/2) 276 points = [ [5,5, 5], [-5,5, 5], [-5,-5, 5], [5,-5, 5], 277 [z,0,-5], [ 0,z,-5], [-z, 0,-5], [0,-z,-5] ] # 8 points 278 faces_a = [ [0,3,4], [1,0,5], [2,1,6], [3,2,7] ] 279 faces_b = [ [4,5,0], [5,6,1], [6,7,2], [7,4,3] ] 280 faces_c = [ [0,1,2], [4,5,6] ] 281 faces_d = [ [2,3,0], [6,7,4] ] 282 faces = faces_a + faces_b + faces_c + faces_d 283 colors = ['teal']*4 + ['indianred']*4 + ['olivedrab']*2 + ['mediumpurple']*2 # 12 colors 284 polytope = Polytope(points) 285 inside, outside = polytope.cut() 286 ax = polytope.plot() 287 ax.scatter(*inside.T, c='g') 288 ax.scatter(*outside.T, c='r') 289 plt.show()