kjao/squish
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Reorganized Cython modules and updated Radial[T] formula

Browse Source
This commit is contained in:
Kenneth Jao 2021-12-02 22:41:12 -05:00
parent 71fb42eb21
commit 85d5486141
17 changed files with 45154 additions and 16565 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

253
squish/_squish/_squish.pxd
View File

@ -1,253 +0,0 @@
cimport numpy as np
# Cython Types.
ctypedef np.int64_t INT_T
ctypedef np.float64_t FLOAT_T
# Stores initialization functions.
cdef struct Init:
IArray (*IArray)(INT_T*, (INT_T, INT_T)) nogil
FArray (*FArray)(FLOAT_T*, (INT_T, INT_T)) nogil
#IList (*IList)() nogil
BitSet (*BitSet)(INT_T) nogil
Vector2D (*Vector2D)(FLOAT_T, FLOAT_T) nogil
Matrix2x2 (*Matrix2x2)(FLOAT_T, FLOAT_T, FLOAT_T, FLOAT_T) nogil
SiteCacheMap (*SiteCacheMap)(INT_T, INT_T, INT_T, INT_T, INT_T) nogil
EdgeCacheMap (*EdgeCacheMap)(INT_T, INT_T, INT_T, INT_T, INT_T, INT_T, INT_T, INT_T,
INT_T, INT_T, INT_T, INT_T, INT_T, INT_T) nogil
VoronoiInfo (*VoronoiInfo)(INT_T [:, ::1], INT_T[:, ::1], FLOAT_T[:, ::1],
FLOAT_T[:, ::1], FLOAT_T[:, ::1], FLOAT_T[:, ::1],
EdgeCacheMap*) nogil
Site (*Site)(INT_T, VoronoiInfo*) nogil
HalfEdge (*HalfEdge)(INT_T, VoronoiInfo*) nogil
# Integer Array psuedo-class for continguous arrays.
cdef struct IArray:
INT_T* arr
(INT_T, INT_T) shape
INT_T (*get)(IArray*, (INT_T, INT_T)) nogil
void (*set)(IArray*, (INT_T, INT_T), INT_T) nogil
# Float Array psuedo-class for continguous arrays.
ctypedef struct FArray:
FLOAT_T* arr
(INT_T, INT_T) shape
FLOAT_T (*get)(FArray*, (INT_T, INT_T)) nogil
void (*set)(FArray*, (INT_T, INT_T), FLOAT_T) nogil
# Simple append-only dynamic integer array.
# ctypedef struct IList:
# INT_T* data
# INT_T size, length
# void (*append)(IList*, INT_T) nogil
# void (*free)(IList*) nogil
# Uses an array of bits to determine if value in set.
ctypedef struct BitSet:
INT_T* bits
bint (*add)(BitSet*, INT_T) nogil
void (*free)(BitSet*) nogil
# Psuedo-operator definitions.
ctypedef Vector2D* (*VectorSelfVecOp)(Vector2D*, Vector2D) nogil
ctypedef Vector2D (*VectorCopyVecOp)(Vector2D*, Vector2D) nogil
ctypedef Vector2D* (*VectorSelfSclOp)(Vector2D*, FLOAT_T) nogil
ctypedef Vector2D (*VectorCopySclOp)(Vector2D*, FLOAT_T) nogil
ctypedef Matrix2x2* (*MatrixSelfMatOp)(Matrix2x2*, Matrix2x2) nogil
ctypedef Matrix2x2 (*MatrixCopyMatOp)(Matrix2x2*, Matrix2x2) nogil
ctypedef Matrix2x2* (*MatrixSelfSclOp)(Matrix2x2*, FLOAT_T) nogil
ctypedef Matrix2x2 (*MatrixCopySclOp)(Matrix2x2*, FLOAT_T) nogil
ctypedef struct VectorSelfOps:
Vector2D* (*neg)(Vector2D*) nogil
VectorSelfVecOp vadd
VectorSelfVecOp vsub
VectorSelfVecOp vmul
VectorSelfVecOp vdiv
Vector2D* (*matmul)(Vector2D*, Matrix2x2) nogil
VectorSelfSclOp sadd
VectorSelfSclOp ssub
VectorSelfSclOp smul
VectorSelfSclOp sdiv
ctypedef struct VectorCopyOps:
Vector2D (*neg)(Vector2D*) nogil
VectorCopyVecOp vadd
VectorCopyVecOp vsub
VectorCopyVecOp vmul
VectorCopyVecOp vdiv
Vector2D (*matmul)(Vector2D*, Matrix2x2) nogil
VectorCopySclOp sadd
VectorCopySclOp ssub
VectorCopySclOp smul
VectorCopySclOp sdiv
ctypedef struct MatrixSelfOps:
Matrix2x2* (*neg)(Matrix2x2*) nogil
MatrixSelfMatOp madd
MatrixSelfMatOp msub
MatrixSelfMatOp mmul
MatrixSelfMatOp mdiv
MatrixSelfMatOp matmul
MatrixSelfSclOp sadd
MatrixSelfSclOp ssub
MatrixSelfSclOp smul
MatrixSelfSclOp sdiv
ctypedef struct MatrixCopyOps:
Matrix2x2 (*neg)(Matrix2x2*) nogil
MatrixCopyMatOp madd
MatrixCopyMatOp msub
MatrixCopyMatOp mmul
MatrixCopyMatOp mdiv
MatrixCopyMatOp matmul
MatrixCopySclOp sadd
MatrixCopySclOp ssub
MatrixCopySclOp smul
MatrixCopySclOp sdiv
# Psuedo-class for a 2-dimensional vector. No orientation.
ctypedef struct Vector2D:
FLOAT_T x, y
VectorSelfOps self
VectorCopyOps copy
bint (*equals)(Vector2D*, Vector2D) nogil
Vector2D (*rot)(Vector2D*) nogil
FLOAT_T (*dot)(Vector2D*, Vector2D) nogil
FLOAT_T (*mag)(Vector2D*) nogil
# Psuedo-class for a 2x2 matrix.
ctypedef struct Matrix2x2:
FLOAT_T a, b, c, d
MatrixSelfOps self
MatrixCopyOps copy
bint (*equals)(Matrix2x2*, Matrix2x2) nogil
Vector2D (*vecmul)(Matrix2x2*, Vector2D) nogil
# Psuedo-class that handles caching for sites.
ctypedef struct SiteCacheMap:
INT_T iarea, iperim, iisoparam, ienergy, iavg_radius
FLOAT_T (*area)(Site*, FLOAT_T) nogil
FLOAT_T (*perim)(Site*, FLOAT_T) nogil
FLOAT_T (*isoparam)(Site*, FLOAT_T) nogil
FLOAT_T (*energy)(Site*, FLOAT_T) nogil
FLOAT_T (*avg_radius)(Site*, FLOAT_T) nogil
# Psuedo-class that handles caching for edges.
ctypedef struct EdgeCacheMap:
INT_T iH, ila, ila_mag, ida, ida_mag, ixij, idVdv, iphi, iB, iF, ii2p,\
ilntan, icsc, size
Matrix2x2 (*H)(HalfEdge*, Matrix2x2) nogil
Vector2D (*la)(HalfEdge*, Vector2D) nogil
Vector2D (*da)(HalfEdge*, Vector2D) nogil
Vector2D (*xij)(HalfEdge*, Vector2D) nogil
Vector2D (*dVdv)(HalfEdge*, Vector2D) nogil
Vector2D (*i2p)(HalfEdge*, Vector2D) nogil
FLOAT_T (*la_mag)(HalfEdge*, FLOAT_T) nogil
FLOAT_T (*da_mag)(HalfEdge*, FLOAT_T) nogil
FLOAT_T (*phi)(HalfEdge*, FLOAT_T) nogil
FLOAT_T (*B)(HalfEdge*, FLOAT_T) nogil
FLOAT_T (*F)(HalfEdge*, FLOAT_T) nogil
FLOAT_T (*lntan)(HalfEdge*, FLOAT_T) nogil
FLOAT_T (*csc)(HalfEdge*, FLOAT_T) nogil
# Psuedo-class to just contain all pertaining info for sites and edges.
ctypedef struct VoronoiInfo:
IArray sites, edges
FArray points, vertices, site_cache, edge_cache
EdgeCacheMap* edge_cache_map
# Psuedo-class for a Site.
ctypedef struct Site:
INT_T arr_index
VoronoiInfo* info
SiteCacheMap* cache
INT_T (*index)(Site*) nogil
Vector2D (*vec)(Site*) nogil
HalfEdge (*edge)(Site*) nogil
INT_T (*edge_num)(Site*) nogil
# Psuedo-class for an HalfEdge.
ctypedef struct HalfEdge:
INT_T orig_arr_index, arr_index
VoronoiInfo* info
EdgeCacheMap* cache
INT_T (*origin_index)(HalfEdge*) nogil
Vector2D (*origin)(HalfEdge*) nogil
Site (*face)(HalfEdge*) nogil
HalfEdge (*next)(HalfEdge*) nogil
HalfEdge (*prev)(HalfEdge*) nogil
HalfEdge (*twin)(HalfEdge*) nogil
Matrix2x2 (*get_H)(HalfEdge*, Site) nogil
cdef class VoronoiContainer:
cdef readonly INT_T n
cdef readonly FLOAT_T w, h, r, energy
cdef FLOAT_T [2] dim
cdef FLOAT_T [:, ::1] points, vertices, site_cache, edge_cache, grad
cdef INT_T [:, ::1] sites, edges
cdef EdgeCacheMap* edge_cache_map
cdef dict __dict__
cdef void calculate_voronoi(VoronoiContainer self,
np.ndarray[FLOAT_T, ndim=2] site_arr) except *
cdef void generate_dcel(VoronoiContainer self) except *
cdef void common_cache(VoronoiContainer self) except *
cdef void precompute(self) except *
cdef void calc_grad(self) except *
cdef void get_statistics(VoronoiContainer self) except *
@staticmethod
cdef inline Matrix2x2 calc_H(HalfEdge, HalfEdge) nogil
@staticmethod
cdef inline bint sign(FLOAT_T [::1], FLOAT_T [::1], FLOAT_T [::1])
cdef class AreaEnergy(VoronoiContainer):
cdef readonly FLOAT_T minimum
cdef void precompute(self) except *
cdef void calc_grad(self) except *
cdef class RadialALEnergy(VoronoiContainer):
cdef void precompute(self) except *
cdef void calc_grad(self) except *
cdef class RadialTEnergy(VoronoiContainer):
cdef void precompute(self) except *
cdef void calc_grad(self) except *
cdef class Calc:
@staticmethod
cdef inline FLOAT_T phi(HalfEdge) nogil
@staticmethod
cdef inline Vector2D I2(HalfEdge, FLOAT_T, FLOAT_T) nogil
@staticmethod
cdef Vector2D radialt_edge_grad(HalfEdge, Site, FLOAT_T) nogil

3
squish/_squish/_squish.pyx
View File

@ -1,3 +0,0 @@
include "core.pyx"
include "voronoi_dcel.pyx"
include "energy.pyx"

385
squish/_squish/core.pyx
View File

@ -1,385 +0,0 @@
import array, scipy.spatial, numpy as np
from cython.parallel import parallel, prange
cimport numpy as np
from cpython cimport array
from libc.stdlib cimport malloc, realloc, calloc, free
from libc.math cimport isnan, NAN, pi as PI, M_PI_2 as PI_2, \
sqrt, log, sin, cos, tan, acos, fabs
from _squish cimport INT_T, FLOAT_T, Init, IArray, FArray, BitSet, Vector2D, Matrix2x2, \
VectorSelfOps, VectorCopyOps, MatrixSelfOps, MatrixCopyOps, \
SiteCacheMap, EdgeCacheMap, VoronoiInfo, Site, HalfEdge
#### Constants ####
INT = np.int64
FLOAT = np.float64
cdef FLOAT_T TAU = 2*PI
# In most cases, the amount of edges relevant to a gradient will
# not exceed this number. However, we assign a growth rate of 8 edges,
# when dynamically allocating.
cdef INT_T EDGE_ARR_SIZE = 32
cdef Init init
init.IArray, init.FArray, init.BitSet, init.Vector2D, init.Matrix2x2 = \
init_iarray, init_farray, init_bitset, init_vector2d, init_matrix2x2
cdef VectorSelfOps VSO
cdef VectorCopyOps VCO
cdef MatrixSelfOps MSO
cdef MatrixCopyOps MCO
VSO.neg, VSO.vadd, VSO.vsub, VSO.vmul, VSO.vdiv, VSO.sadd, VSO.ssub, VSO.smul, VSO.sdiv = \
v_neg_s, v_vadd_s, v_vsub_s, v_vmul_s, v_vdiv_s, v_sadd_s, v_ssub_s, v_smul_s, v_sdiv_s
VSO.matmul = v_matmul_s
VCO.neg, VCO.vadd, VCO.vsub, VCO.vmul, VCO.vdiv, VCO.sadd, VCO.ssub, VCO.smul, VCO.sdiv = \
v_neg_c, v_vadd_c, v_vsub_c, v_vmul_c, v_vdiv_c, v_sadd_c, v_ssub_c, v_smul_c, v_sdiv_c
VCO.matmul = v_matmul_c
MSO.neg, MSO.madd, MSO.msub, MSO.mmul, MSO.mdiv, MSO.sadd, MSO.ssub, MSO.smul, MSO.sdiv = \
m_neg_s, m_madd_s, m_msub_s, m_mmul_s, m_mdiv_s, m_sadd_s, m_ssub_s, m_smul_s, m_sdiv_s
MSO.matmul = m_matmul_s
MCO.neg, MCO.madd, MCO.msub, MCO.mmul, MCO.mdiv, MCO.sadd, MCO.ssub, MCO.smul, MCO.sdiv = \
m_neg_c, m_madd_c, m_msub_c, m_mmul_c, m_mdiv_c, m_sadd_c, m_ssub_c, m_smul_c, m_sdiv_c
MCO.matmul = m_matmul_c
cdef Vector2D NAN_VECTOR = init.Vector2D(NAN, NAN)
cdef Matrix2x2 NAN_MATRIX = init.Matrix2x2(NAN, NAN, NAN, NAN)
cdef FLOAT_T[18] SYMM = [0,0, 1,0, 1,1, 0,1, -1,1, -1,0, -1,-1, 0,-1, 1,-1]
cdef Matrix2x2 R = init.Matrix2x2(0, -1, 1, 0)
"""
If bound checking is desired, uncomment out ..._valid_indices functions.
"""
#### IArray Methods ####
cdef inline IArray init_iarray(INT_T* arr, (INT_T, INT_T) shape) nogil:
cdef IArray iarray
iarray.arr, iarray.shape = arr, shape
iarray.get = iarray_get
iarray.set = iarray_set
return iarray
cdef inline bint iarray_valid_indices(IArray* self, (INT_T, INT_T) index) nogil:
if index[0] > self.shape[0] or index[1] > self.shape[1]:
with gil:
raise IndexError(f"Index out of range for IArray with shape {self.shape}")
cdef inline INT_T iarray_get(IArray* self, (INT_T, INT_T) index) nogil:
#iarray_valid_indices(&self, index)
return self.arr[index[0]*self.shape[1] + index[1]]
cdef inline void iarray_set(IArray* self, (INT_T, INT_T) index, INT_T val) nogil:
#iarray_valid_indices(&self, index)
self.arr[index[0]*self.shape[1] + index[1]] = val
#### FArray Methods ####
cdef inline FArray init_farray(FLOAT_T* arr, (INT_T, INT_T) shape) nogil:
cdef FArray farray
farray.arr, farray.shape = arr, shape
farray.get = farray_get
farray.set = farray_set
return farray
cdef inline bint farray_valid_indices(FArray* self, (INT_T, INT_T) index) nogil:
if index[0] > self.shape[0] or index[1] > self.shape[1]:
with gil:
raise IndexError(f"Index out of range for FArray with shape {self.shape}")
cdef inline FLOAT_T farray_get(FArray* self, (INT_T, INT_T) index) nogil:
#iarray_valid_indices(&self, index)
return self.arr[index[0]*self.shape[1] + index[1]]
cdef inline void farray_set(FArray* self, (INT_T, INT_T) index, FLOAT_T val) nogil:
#iarray_valid_indices(&self, index)
self.arr[index[0]*self.shape[1] + index[1]] = val
#### IList Methods ####
# cdef inline IList init_ilist() nogil:
# cdef IList ilist
# ilist.size = EDGE_ARR_SIZE
# ilist.length = 0
# ilist.data = <INT_T*> malloc(self.size * sizeof(INT_T))
# ilist.append, ilist.free = ilist_append, ilist_free
# return ilist
# cdef inline void ilist_append(IList* self, INT_T) nogil:
# if self.size == self.length:
# ilist.data = <INT_T*> realloc((self.size+8) * sizeof(INT_T))
# self.size += 8
# self.data[self.length] == INT_T
# self.length += 1
# cdef inline void ilist_free(IList* self) nogil:
# free(self.data)
#### BitSet Methods ####
cdef inline BitSet init_bitset(INT_T elements) nogil:
cdef BitSet bitset
bitset.bits = <INT_T*> calloc(((elements/sizeof(INT_T))+1), sizeof(INT_T))
bitset.add, bitset.free = bitset_add, bitset_free
return bitset
cdef inline bint bitset_add(BitSet* self, INT_T val) nogil:
cdef INT_T index, rel_index, old
index = val/sizeof(INT_T)
old = self.bits[index]
rel_index = val - index*sizeof(INT_T)
self.bits[index] = (1 << rel_index) | old # New value.
return old == self.bits[index] # Means 1 was already there.
cdef inline void bitset_free(BitSet* self) nogil:
free(self.bits)
#### Vector2D Methods ####
"""
Prefix 'v' stands for vector, element by element operation.
Prefix 's' stands for scalar, broadcasted operation.
Suffix 'w' stands for write, overwriting current value.
Suffix 'n' stands for new, copying to a new location.
While it's possible to chain 'new' operations, when possible,
avoid this, so fewer objects are needed.
"""
cdef inline Vector2D init_vector2d(FLOAT_T x, FLOAT_T y) nogil:
cdef Vector2D vec
vec.x, vec.y = x, y
vec.self, vec.copy = VSO, VCO
vec.equals, vec.rot, vec.dot, vec.mag = v_equals, rot, dot, mag
return vec
cdef inline bint v_equals(Vector2D* self, Vector2D w) nogil:
return ((self.x == w.x) and (self.y == w.y))
cdef inline Vector2D* v_neg_s(Vector2D* self) nogil:
self.x = -self.x
self.y = -self.y
return self
cdef inline Vector2D* v_vadd_s(Vector2D* self, Vector2D w) nogil:
self.x += w.x
self.y += w.y
return self
cdef inline Vector2D* v_vsub_s(Vector2D* self, Vector2D w) nogil:
self.x -= w.x
self.y -= w.y
return self
cdef inline Vector2D* v_vmul_s(Vector2D* self, Vector2D w) nogil:
self.x *= w.x
self.y *= w.y
return self
cdef inline Vector2D* v_vdiv_s(Vector2D* self, Vector2D w) nogil:
self.x /= w.x
self.y /= w.y
return self
cdef inline Vector2D* v_sadd_s(Vector2D* self, FLOAT_T s) nogil:
self.x += s
self.y += s
return self
cdef inline Vector2D* v_ssub_s(Vector2D* self, FLOAT_T s) nogil:
self.x -= s
self.y -= s
return self
cdef inline Vector2D* v_smul_s(Vector2D* self, FLOAT_T s) nogil:
self.x *= s
self.y *= s
return self
cdef inline Vector2D* v_sdiv_s(Vector2D* self, FLOAT_T s) nogil:
self.x /= s
self.y /= s
return self
cdef inline Vector2D* v_matmul_s(Vector2D* self, Matrix2x2 m) nogil:
self.x, self.y = self.x*m.a + self.y*m.c, self.x*m.b + self.y*m.d
return self
cdef inline Vector2D v_neg_c(Vector2D* self) nogil:
return init.Vector2D(-self.x, -self.y)
cdef inline Vector2D v_vadd_c(Vector2D* self, Vector2D w) nogil:
return init.Vector2D(self.x + w.x, self.y + w.y)
cdef inline Vector2D v_vsub_c(Vector2D* self, Vector2D w) nogil:
return init.Vector2D(self.x - w.x, self.y - w.y)
cdef inline Vector2D v_vmul_c(Vector2D* self, Vector2D w) nogil:
return init.Vector2D(self.x * w.x, self.y * w.y)
cdef inline Vector2D v_vdiv_c(Vector2D* self, Vector2D w) nogil:
return init.Vector2D(self.x / w.x, self.y / w.y)
cdef inline Vector2D v_sadd_c(Vector2D* self, FLOAT_T s) nogil:
return init.Vector2D(self.x + s, self.y + s)
cdef inline Vector2D v_ssub_c(Vector2D* self, FLOAT_T s) nogil:
return init.Vector2D(self.x + s, self.y + s)
cdef inline Vector2D v_smul_c(Vector2D* self, FLOAT_T s) nogil:
return init.Vector2D(self.x * s, self.y * s)
cdef inline Vector2D v_sdiv_c(Vector2D* self, FLOAT_T s) nogil:
return init.Vector2D(self.x / s, self.y / s)
cdef inline Vector2D v_matmul_c(Vector2D* self, Matrix2x2 m) nogil:
return init.Vector2D(
self.x*m.a + self.y*m.c, self.x*m.b + self.y*m.d
)
cdef inline Vector2D rot(Vector2D* self) nogil:
return init.Vector2D(-self.y, self.x)
cdef inline FLOAT_T dot(Vector2D* self, Vector2D w) nogil:
return self.x*w.x + self.y*w.y
cdef inline FLOAT_T mag(Vector2D* self) nogil:
return <FLOAT_T>sqrt(<double>(self.x*self.x + self.y*self.y))
#### Matrix2x2 Methods ####
cdef inline Matrix2x2 init_matrix2x2(FLOAT_T a, FLOAT_T b, FLOAT_T c, FLOAT_T d) nogil:
cdef Matrix2x2 matrix
matrix.a, matrix.b, matrix.c, matrix.d = a, b, c, d
matrix.self, matrix.copy = MSO, MCO
matrix.equals, matrix.vecmul = m_equals, m_vecmul
return matrix
cdef inline bint m_equals(Matrix2x2* self, Matrix2x2 m) nogil:
return (
(self.a == m.a) and (self.b == m.b) and (self.c == m.c) and (self.d == m.d)
)
cdef inline Vector2D m_vecmul(Matrix2x2* self, Vector2D v) nogil:
return init.Vector2D(
self.a*v.x + self.b*v.y, self.c*v.x + self.d*v.y
)
cdef inline Matrix2x2* m_neg_s(Matrix2x2* self) nogil:
self.a, self.b, self.c, self.d = -self.a, -self.b, -self.c, -self.d
return self
cdef inline Matrix2x2* m_madd_s(Matrix2x2* self, Matrix2x2 m) nogil:
self.a += m.a
self.b += m.b
self.c += m.c
self.d += m.d
return self
cdef inline Matrix2x2* m_msub_s(Matrix2x2* self, Matrix2x2 m) nogil:
self.a -= m.a
self.b -= m.b
self.c -= m.c
self.d -= m.d
return self
cdef inline Matrix2x2* m_mmul_s(Matrix2x2* self, Matrix2x2 m) nogil:
self.a *= m.a
self.b *= m.b
self.c *= m.c
self.d *= m.d
return self
cdef inline Matrix2x2* m_mdiv_s(Matrix2x2* self, Matrix2x2 m) nogil:
self.a /= m.a
self.b /= m.b
self.c /= m.c
self.d /= m.d
return self
cdef inline Matrix2x2* m_sadd_s(Matrix2x2* self, FLOAT_T s) nogil:
self.a += s
self.b += s
self.c += s
self.d += s
return self
cdef inline Matrix2x2* m_ssub_s(Matrix2x2* self, FLOAT_T s) nogil:
self.a -= s
self.b -= s
self.c -= s
self.d -= s
return self
cdef inline Matrix2x2* m_smul_s(Matrix2x2* self, FLOAT_T s) nogil:
self.a *= s
self.b *= s
self.c *= s
self.d *= s
return self
cdef inline Matrix2x2* m_sdiv_s(Matrix2x2* self, FLOAT_T s) nogil:
self.a /= s
self.b /= s
self.c /= s
self.d /= s
return self
cdef inline Matrix2x2* m_matmul_s(Matrix2x2* self, Matrix2x2 m) nogil:
self.a, self.b, self.c, self.d = \
self.a*m.a + self.b*m.c, self.a*m.b + self.b*m.d, \
self.c*m.a + self.d*m.c, self.c*m.b + self.d*m.d
return self
cdef inline Matrix2x2 m_neg_c(Matrix2x2* self) nogil:
return init.Matrix2x2(-self.a, -self.b, -self.c, -self.d)
cdef inline Matrix2x2 m_madd_c(Matrix2x2* self, Matrix2x2 m) nogil:
return init.Matrix2x2(self.a+m.a, self.b+m.b, self.c+m.c, self.d+m.d)
cdef inline Matrix2x2 m_msub_c(Matrix2x2* self, Matrix2x2 m) nogil:
return init.Matrix2x2(self.a-m.a, self.b-m.b, self.c-m.c, self.d-m.d)
cdef inline Matrix2x2 m_mmul_c(Matrix2x2* self, Matrix2x2 m) nogil:
return init.Matrix2x2(self.a*m.a, self.b*m.b, self.c*m.c, self.d*m.d)
cdef inline Matrix2x2 m_mdiv_c(Matrix2x2* self, Matrix2x2 m) nogil:
return init.Matrix2x2(self.a/m.a, self.b/m.b, self.c/m.c, self.d/m.d)
cdef inline Matrix2x2 m_sadd_c(Matrix2x2* self, FLOAT_T s) nogil:
return init.Matrix2x2(self.a+s, self.b+s, self.c+s, self.d+s)
cdef inline Matrix2x2 m_ssub_c(Matrix2x2* self, FLOAT_T s) nogil:
return init.Matrix2x2(self.a-s, self.b-s, self.c-s, self.d-s)
cdef inline Matrix2x2 m_smul_c(Matrix2x2* self, FLOAT_T s) nogil:
return init.Matrix2x2(self.a*s, self.b*s, self.c*s, self.d*s)
cdef inline Matrix2x2 m_sdiv_c(Matrix2x2* self, FLOAT_T s) nogil:
return init.Matrix2x2(self.a/s, self.b/s, self.c/s, self.d/s)
cdef inline Matrix2x2 m_matmul_c(Matrix2x2* self, Matrix2x2 m) nogil:
return init.Matrix2x2(
self.a*m.a + self.b*m.c, self.a*m.b + self.b*m.d,
self.c*m.a + self.d*m.c, self.c*m.b + self.d*m.d
)

338
squish/_squish/energy.pyx
View File

@ -1,338 +0,0 @@
cdef class AreaEnergy(VoronoiContainer):
"""
Class for formulas relevant to the Area energy.
:param n: [int] how many sites to generate.
:param w: [float] width of the bounding domain.
:param h: [float] height of the bounding domain.
:param r: [float] radius of zero energy circle.
:param sites: [np.ndarray] collection of sites.
"""
attr_str = "area"
title_str = "Area"
def __init__(AreaEnergy self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r,
np.ndarray[FLOAT_T, ndim=2] site_arr):
self.edge_cache_map = &AREA_EDGE_CACHE_MAP
self.energy = 0.0
super().__init__(n, w, h, r, site_arr)
self.minimum = (<FLOAT_T>n)*(w*h/(<FLOAT_T>n)-PI*r**2)**2
cdef void precompute(self) except *:
cdef VoronoiInfo info = init.VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
cdef Site xi
cdef HalfEdge em, e, ep
cdef Vector2D vdiff
cdef FLOAT_T A = PI*self.r**2
cdef FLOAT_T energy = 0
cdef INT_T i, j
for i in prange(self.sites.shape[0], nogil=True):
xi = init.Site(i, &info)
e = xi.edge(&xi)
xi.cache.energy(&xi,
(xi.cache.area(&xi, NAN) - A)**2
)
if i < self.n:
energy += xi.cache.energy(&xi, NAN)
for j in range(xi.edge_num(&xi)):
em, ep = e.prev(&e), e.next(&e)
vdiff = em.origin(&em)
vdiff.self.vsub(&vdiff, ep.origin(&ep))
e.cache.dVdv(&e, R.vecmul(&R, vdiff))
e.cache.H(&e, VoronoiContainer.calc_H(em, e))
e = e.next(&e)
self.energy = energy
cdef void calc_grad(self) except *:
cdef VoronoiInfo info = init.VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
cdef Site xi, xf
cdef HalfEdge e, f
cdef Vector2D dedxi_p
cdef BitSet edge_set
cdef INT_T num_edges = self.edges.shape[0]
cdef FLOAT_T A = PI*self.r**2
cdef FLOAT_T [:, ::1] dedx = np.zeros((self.n, 2), dtype=FLOAT)
cdef INT_T i, j
for i in prange(self.n, nogil=True):
xi = init.Site(i, &info)
e = xi.edge(&xi)
edge_set = init.BitSet(num_edges)
for j in range(xi.edge_num(&xi)): # Looping through site edges.
f = e
while True: # Circling this vertex.
if not edge_set.add(&edge_set, f.arr_index):
xf = f.face(&f)
dedxi_p = f.cache.dVdv(&f, NAN_VECTOR) #dVdv
dedxi_p.self.smul(&dedxi_p, xf.cache.area(&xf, NAN) - A)
dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX))
dedx[i][0] -= dedxi_p.x
dedx[i][1] -= dedxi_p.y
f = f.twin(&f)
f = f.next(&f)
if f.arr_index == e.arr_index:
break
e = e.next(&e)
edge_set.free(&edge_set)
self.grad = dedx
cdef class RadialALEnergy(VoronoiContainer):
"""
Class for formulas relevant to the Area energy.
:param n: [int] how many sites to generate.
:param w: [float] width of the bounding domain.
:param h: [float] height of the bounding domain.
:param r: [float] radius of zero energy circle.
:param sites: [np.ndarray] collection of sites.
"""
attr_str = "radial-al"
title_str = "Radial[AL]"
def __init__(AreaEnergy self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r,
np.ndarray[FLOAT_T, ndim=2] site_arr):
#self.edge_cache_map = &AREA_EDGE_CACHE_MAP
self.energy = 0.0
super().__init__(n, w, h, r, site_arr)
cdef void precompute(self) except *:
cdef VoronoiInfo info = init.VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
pass
cdef void calc_grad(self) except *:
cdef VoronoiInfo info = init.VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
pass
cdef class RadialTEnergy(VoronoiContainer):
"""
Class for formulas relevant to the Area energy.
:param n: [int] how many sites to generate.
:param w: [float] width of the bounding domain.
:param h: [float] height of the bounding domain.
:param r: [float] radius of zero energy circle.
:param sites: [np.ndarray] collection of sites.
"""
attr_str = "radial-t"
title_str = "Radial[T]"
def __init__(AreaEnergy self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r,
np.ndarray[FLOAT_T, ndim=2] site_arr):
self.edge_cache_map = &RADIALT_EDGE_CACHE_MAP
self.energy = 0.0
super().__init__(n, w, h, r, site_arr)
cdef void precompute(self) except *:
cdef VoronoiInfo info = init.VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
cdef Site xi
cdef HalfEdge em, e
cdef Vector2D Rnla
# All energy has a 2pir_0 term.
cdef FLOAT_T [:] site_energy = np.full(self.sites.shape[0], TAU*self.r**2)
cdef FLOAT_T [:] avg_radii = np.zeros(self.sites.shape[0])
cdef FLOAT_T energy, r0, t, tp, B, lntan, csc
energy, r0 = 0, self.r
cdef INT_T i, j
for i in prange(self.sites.shape[0], nogil=True):
xi = init.Site(i, &info)
e = xi.edge(&xi)
for j in range(xi.edge_num(&xi)):
em = e.prev(&e)
e.cache.H(&e, VoronoiContainer.calc_H(em, e))
t = Calc.phi(e)
e.cache.phi(&e, t)
Rnla = e.cache.la(&e, NAN_VECTOR)
Rnla.self.neg(&Rnla)
Rnla = Rnla.rot(&Rnla)
if Rnla.x < 0:
e.cache.B(&e, -<FLOAT_T>acos(<double>(Rnla.y/e.cache.la_mag(&e, NAN))))
else:
e.cache.B(&e, <FLOAT_T>acos(<double>(Rnla.y/e.cache.la_mag(&e, NAN))))
e.cache.i2p(&e, Calc.I2(e, r0, t))
e = e.next(&e)
# For looping again to calculate integrals.
em = xi.edge(&xi)
for j in range(xi.edge_num(&xi)):
e = em.next(&em)
B = em.cache.B(&em, NAN)
t, tp = em.cache.phi(&em, NAN), e.cache.phi(&e, NAN)
lntan = <FLOAT_T>(log(fabs(tan(<double>((tp+B)/2))))) - \
<FLOAT_T>(log(fabs(tan(<double>((t+B)/2)))))
csc = 1/(<FLOAT_T>(sin(<double>(tp+B)))) - \
1/(<FLOAT_T>(sin(<double>(t+B))))
em.cache.lntan(&em, lntan)
em.cache.csc(&em, csc)
avg_radii[i] += (em.cache.F(&em, NAN)/em.cache.la_mag(&em, NAN))*lntan
em = em.next(&em)
site_energy[i] += 2*(xi.cache.area(&xi, NAN) - r0*avg_radii[i])
xi.cache.avg_radius(&xi, avg_radii[i]/TAU)
xi.cache.energy(&xi, site_energy[i])
if i < self.n:
energy += site_energy[i]
self.energy = energy
cdef void calc_grad(self) except *:
cdef VoronoiInfo info = init.VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
cdef Site xi
cdef HalfEdge e, fm, f
cdef Vector2D dedxi_p
cdef BitSet edge_set
cdef INT_T num_edges = self.edges.shape[0]
cdef FLOAT_T r0 = self.r
cdef FLOAT_T [:, ::1] dedx = np.zeros((self.n, 2), dtype=FLOAT)
cdef INT_T i, j
for i in prange(self.n, nogil=True):
xi = init.Site(i, &info)
e = xi.edge(&xi)
edge_set = init.BitSet(num_edges)
for j in range(xi.edge_num(&xi)): # Looping through site edges.
f = e
while True: # Circling this vertex.
fm = f.prev(&f)
if not edge_set.add(&edge_set, f.arr_index):
dedxi_p = Calc.radialt_edge_grad(f, xi, r0)
dedx[i][0] -= dedxi_p.x
dedx[i][1] -= dedxi_p.y
if not edge_set.add(&edge_set, fm.arr_index):
dedxi_p = Calc.radialt_edge_grad(fm, xi, r0)
dedx[i][0] -= dedxi_p.x
dedx[i][1] -= dedxi_p.y
f = f.twin(&f)
f = f.next(&f)
if f.arr_index == e.arr_index:
break
e = e.next(&e)
edge_set.free(&edge_set)
self.grad = dedx
cdef class Calc:
@staticmethod
cdef inline FLOAT_T phi(HalfEdge e) nogil:
cdef Vector2D da = e.cache.da(&e, NAN_VECTOR)
cdef FLOAT_T angle = <FLOAT_T>acos(<double>(da.x/e.cache.da_mag(&e, NAN)))
return angle if da.y >= 0 else TAU - angle
@staticmethod
cdef inline Vector2D I2(HalfEdge e, FLOAT_T r0, FLOAT_T t) nogil:
cdef Vector2D Rda = e.cache.da(&e, NAN_VECTOR)
Rda = Rda.rot(&Rda)
Rda.self.sdiv(&Rda, e.cache.da_mag(&e, NAN))
return Rda
@staticmethod
cdef Vector2D radialt_edge_grad(HalfEdge e, Site xi, FLOAT_T r0) nogil:
cdef Site xe
cdef HalfEdge ep
cdef Vector2D Rda, i2ps, fp, gterms, q
cdef Matrix2x2 ha, hap, hdiff
cdef FLOAT_T t1, t2, lntan, csc, sinB, cosB, sinBp, cosBp, F, A, B
xe = e.face(&e)
ep = e.next(&e)
F, A, B = e.cache.F(&e, NAN), e.cache.la_mag(&e, NAN), e.cache.B(&e, NAN)
t1, t2 = e.cache.phi(&e, NAN), ep.cache.phi(&ep, NAN)
lntan, csc = e.cache.lntan(&e, NAN), e.cache.csc(&e, NAN)
sinB, cosB = <FLOAT_T>(sin(<double>(B))), <FLOAT_T>(cos(<double>(B)))
sinBp, cosBp = <FLOAT_T>(sin(<double>(B-PI_2))), \
<FLOAT_T>(cos(<double>(B-PI_2)))
ha, hap = e.get_H(&e, xi), ep.get_H(&ep, xi)
hdiff = hap.copy.msub(&hap, ha)
# If edge is part of differentiated site.
if xe.index(&xe) == xi.index(&xi):
ha.self.msub(&ha, init.Matrix2x2(1.0, 0.0, 0.0, 1.0))
hap.self.msub(&hap, init.Matrix2x2(1.0, 0.0, 0.0, 1.0))
i2ps = ep.cache.i2p(&ep, NAN_VECTOR)
i2ps.self.matmul(&i2ps, hap)
q = e.cache.i2p(&e, NAN_VECTOR)
q.self.matmul(&q, ha)
i2ps.self.vsub(&i2ps, q)
Rda = e.cache.da(&e, NAN_VECTOR)
Rda = Rda.rot(&Rda)
fp = e.cache.la(&e, NAN_VECTOR)
fp.self.matmul(&fp, R.copy.matmul(&R, ha))
fp.self.vadd(&fp, Rda.copy.matmul(&Rda, hdiff))
fp.self.smul(&fp, lntan/A)
gterms = init.Vector2D(
cosBp*lntan + sinBp*csc,
cosB*lntan + sinB*csc
)
gterms.self.smul(&gterms, -F/A**2)
gterms = gterms.rot(&gterms)
gterms.self.matmul(&gterms, hdiff)
fp.self.vadd(&fp, gterms)
i2ps.self.vadd(&i2ps, fp)
i2ps.self.smul(&i2ps, -2*r0)
return i2ps

689
squish/_squish/voronoi_dcel.pyx
View File

@ -1,689 +0,0 @@
from _squish cimport SiteCacheMap, EdgeCacheMap, VoronoiInfo, Site, HalfEdge
#### Constants ####
init.SiteCacheMap, init.EdgeCacheMap, init.VoronoiInfo, init.Site, init.HalfEdge = \
init_sitecachemap, init_edgecachemap, init_voronoiinfo, init_site, init_halfedge
cdef SiteCacheMap SITE_CACHE_MAP = init.SiteCacheMap(0, 1, 2, 3, 4)
cdef EdgeCacheMap AREA_EDGE_CACHE_MAP = init.EdgeCacheMap(0, 4, 6, 8, 10, -1, 12, 13,
-1, -1, -1, -1, -1, 14)
cdef EdgeCacheMap RADIALT_EDGE_CACHE_MAP = init.EdgeCacheMap(0, 4, 6, 8, -1, 10, 12, 13,
14, 15, 16, 17, 18, 19)
#### SiteCacheMap Methods ####
cdef inline SiteCacheMap init_sitecachemap(INT_T iarea, INT_T iperim, INT_T iisoparam,
INT_T ienergy, INT_T iavg_radius) nogil:
cdef SiteCacheMap sc
sc.iarea, sc.iperim, sc.iisoparam, sc.ienergy, sc.iavg_radius = \
iarea, iperim, iisoparam, ienergy, iavg_radius
sc.area, sc.perim, sc.isoparam, sc.energy, sc.avg_radius = \
area, perim, isoparam, energy, avg_radius
return sc
cdef inline FLOAT_T area(Site* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.iarea)
)
else:
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.iarea), val)
return val
cdef inline FLOAT_T perim(Site* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.iperim)
)
else:
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.iperim), val)
return val
cdef inline FLOAT_T isoparam(Site* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.iisoparam)
)
else:
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.iisoparam), val)
return val
cdef inline FLOAT_T energy(Site* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.ienergy)
)
else:
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.ienergy), val)
return val
cdef inline FLOAT_T avg_radius(Site* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.iavg_radius)
)
else:
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.iavg_radius), val)
return val
#### EdgeCacheMap Methods ####
cdef inline EdgeCacheMap init_edgecachemap(INT_T iH, INT_T ila, INT_T ida, INT_T ixij,
INT_T idVdv, INT_T ii2p, INT_T ila_mag, INT_T ida_mag, INT_T iphi, INT_T iB,
INT_T iF, INT_T ilntan, INT_T icsc, INT_T size) nogil:
cdef EdgeCacheMap ec
ec.iH, ec.ila, ec.ida, ec.ixij, ec.idVdv, ec.ii2p, ec.ila_mag, ec.ida_mag, ec.iphi, \
ec.iB, ec.iF, ec.ilntan, ec.icsc = iH, ila, ida, ixij, idVdv, ii2p, \
ila_mag, ida_mag, iphi, iB, iF, ilntan, icsc
ec.size = size
ec.H, ec.la, ec.da, ec.xij, ec.dVdv, ec.i2p, ec.la_mag, ec.da_mag, ec.phi, ec.B, ec.F, \
ec.lntan, ec.csc = H, la, da, xij, dVdv, i2p, la_mag, da_mag, phi, B, F, lntan, csc
return ec
cdef inline Matrix2x2 H(HalfEdge* self, Matrix2x2 val) nogil:
if isnan(<double>val.a):
return init.Matrix2x2(
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.iH)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.iH+1)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.iH+2)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.iH+3)
),
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.iH), val.a)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.iH+1), val.b)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.iH+2), val.c)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.iH+3), val.d)
return val
cdef inline Vector2D la(HalfEdge* self, Vector2D val) nogil:
if isnan(<double>val.x):
return init.Vector2D(
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ila)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ila+1)
)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ila), val.x)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ila+1), val.y)
return val
cdef inline Vector2D da(HalfEdge* self, Vector2D val) nogil:
if isnan(<double>val.x):
return init.Vector2D(
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ida)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ida+1)
)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ida), val.x)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ida+1), val.y)
return val
cdef inline Vector2D xij(HalfEdge* self, Vector2D val) nogil:
if isnan(<double>val.x):
return init.Vector2D(
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ixij)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ixij+1)
)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ixij), val.x)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ixij+1), val.y)
return val
cdef inline Vector2D dVdv(HalfEdge* self, Vector2D val) nogil:
if isnan(<double>val.x):
return init.Vector2D(
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.idVdv)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.idVdv+1)
)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.idVdv), val.x)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.idVdv+1), val.y)
return val
cdef inline Vector2D i2p(HalfEdge* self, Vector2D val) nogil:
if isnan(<double>val.x):
return init.Vector2D(
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ii2p)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ii2p+1)
)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ii2p), val.x)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ii2p+1), val.y)
return val
cdef inline FLOAT_T la_mag(HalfEdge* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ila_mag)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ila_mag), val)
return val
cdef inline FLOAT_T da_mag(HalfEdge* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ida_mag)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ida_mag), val)
return val
cdef inline FLOAT_T phi(HalfEdge* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.iphi)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.iphi), val)
return val
cdef inline FLOAT_T B(HalfEdge* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.iB)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.iB), val)
return val
cdef inline FLOAT_T F(HalfEdge* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.iF)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.iF), val)
return val
cdef inline FLOAT_T lntan(HalfEdge* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ilntan)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ilntan), val)
return val
cdef inline FLOAT_T csc(HalfEdge* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.icsc)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.icsc), val)
return val
#### VoronoiInfo Methods ####
cdef inline VoronoiInfo init_voronoiinfo(INT_T [:, ::1] sites, INT_T [:, ::1] edges,
FLOAT_T [:, ::1] points, FLOAT_T [:, ::1] vertices,
FLOAT_T [:, ::1] site_cache, FLOAT_T [:, ::1] edge_cache,
EdgeCacheMap* edge_cache_map) nogil:
cdef VoronoiInfo info
info.sites = init_iarray(&sites[0,0], (<INT_T>sites.shape[0], <INT_T>sites.shape[1]))
info.edges = init_iarray(&edges[0,0], (<INT_T>edges.shape[0], <INT_T>edges.shape[1]))
info.points = init_farray(&points[0,0], (<INT_T>points.shape[0], <INT_T>points.shape[1]))
info.vertices = init_farray(&vertices[0,0],
(<INT_T>vertices.shape[0], <INT_T>vertices.shape[1])
)
info.site_cache = init_farray(&site_cache[0,0],
(<INT_T>site_cache.shape[0], <INT_T>site_cache.shape[1])
)
info.edge_cache = init_farray(&edge_cache[0,0],
(<INT_T>edge_cache.shape[0], <INT_T>edge_cache.shape[1])
)
info.edge_cache_map = edge_cache_map
return info
#### Site Methods ####
cdef inline Site init_site(INT_T arr_index, VoronoiInfo* info) nogil:
cdef Site site
site.arr_index, site.info, site.cache = arr_index, info, &SITE_CACHE_MAP
site.index, site.vec, site.edge, site.edge_num = index, vec, edge, edge_num
return site
cdef inline INT_T index(Site* self) nogil:
return self.info.sites.get(&self.info.sites, (self.arr_index, 0))
cdef inline Vector2D vec(Site* self) nogil:
return init.Vector2D(
self.info.points.get(&self.info.points, (self.index(self), 0)),
self.info.points.get(&self.info.points, (self.index(self), 1))
)
cdef inline HalfEdge edge(Site* self) nogil:
return init.HalfEdge(
self.info.sites.get(&self.info.sites, (self.arr_index, 1)), self.info
)
cdef inline INT_T edge_num(Site* self) nogil:
return self.info.sites.get(&self.info.sites, (self.arr_index, 2))
#### HalfEdge Methods ####
cdef inline HalfEdge init_halfedge(INT_T arr_index, VoronoiInfo* info) nogil:
cdef HalfEdge edge
edge.arr_index, edge.info, edge.cache = arr_index, info, info.edge_cache_map
edge.orig_arr_index = arr_index
edge.origin_index, edge.origin, edge.face, edge.next, edge.prev, edge.twin, edge.get_H = \
origin_index, origin, face, edge_next, prev, twin, get_H
return edge
cdef inline INT_T origin_index(HalfEdge* self) nogil:
return self.info.edges.get(&self.info.edges, (self.arr_index, 0))
cdef inline Vector2D origin(HalfEdge* self) nogil:
return init.Vector2D(
self.info.vertices.get(&self.info.vertices, (self.origin_index(self), 0)),
self.info.vertices.get(&self.info.vertices, (self.origin_index(self), 1))
)
cdef inline Site face(HalfEdge* self) nogil:
return init.Site(
self.info.edges.get(&self.info.edges, (self.arr_index, 1)), self.info
)
cdef inline HalfEdge edge_next(HalfEdge* self) nogil:
return init.HalfEdge(
self.info.edges.get(&self.info.edges, (self.arr_index, 2)), self.info
)
cdef inline HalfEdge prev(HalfEdge* self) nogil:
return init.HalfEdge(
self.info.edges.get(&self.info.edges, (self.arr_index, 3)), self.info
)
cdef inline HalfEdge twin(HalfEdge* self) nogil:
return init.HalfEdge(
self.info.edges.get(&self.info.edges, (self.arr_index, 4)), self.info
)
cdef inline Matrix2x2 get_H(HalfEdge* self, Site xi) nogil:
cdef INT_T this_e = self.origin_index(self)
cdef HalfEdge s_e = xi.edge(&xi)
cdef INT_T i
for i in range(xi.edge_num(&xi)):
if s_e.origin_index(&s_e) == this_e:
return s_e.cache.H(&s_e, NAN_MATRIX)
s_e = s_e.next(&s_e)
return init.Matrix2x2(0.0, 0.0, 0.0, 0.0)
cdef class VoronoiContainer:
"""
Class for Voronoi diagrams, stored in a modified DCEL.
:param n: [int] how many sites to generate.
:param w: [float] width of the bounding domain.
:param h: [float] height of the bounding domain.
:param r: [float] radius of zero energy circle.
:param sites: np.ndarray collection of sites.
"""
def __init__(VoronoiContainer self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r, object site_arr):
self.n, self.w, self.h, self.r = n, w, h, r
self.dim = [w, h]
self.calculate_voronoi(site_arr.astype(FLOAT))
self.generate_dcel()
self.common_cache()
self.precompute()
self.calc_grad()
self.get_statistics()
cdef void calculate_voronoi(VoronoiContainer self,
np.ndarray[FLOAT_T, ndim=2] site_arr) except *:
"""
Does all necessary computation and caching once points are set.
:param site_arr: initial points for this container.
"""
global SYMM
cdef np.ndarray[FLOAT_T, ndim=2] symm = np.asarray(SYMM).reshape(9,2)
cdef np.ndarray[FLOAT_T, ndim=1] dim = np.asarray(self.dim)
cdef np.ndarray[FLOAT_T, ndim=2] full_site_arr = np.empty((self.n*9+8, 2), dtype=FLOAT)
# Generate periodic sites and sites that bound periodic sites.
cdef INT_T i
for i in range(9):
full_site_arr[self.n*i:self.n*(i+1)] = site_arr + symm[i]*dim
if i > 0:
full_site_arr[9*self.n+i-1] = dim/2 + 2*dim*symm[i]
# Use SciPy to compute the Voronoi set.
self.scipy_vor = scipy.spatial.Voronoi(full_site_arr)
self.points = self.scipy_vor.points
self.vertices = self.scipy_vor.vertices
cdef void generate_dcel(VoronoiContainer self) except *:
cdef INT_T npoints = self.n*9+8
cdef array.array int_tmplt = array.array('q', [])
cdef np.ndarray[INT_T, ndim=1] offsets = np.zeros(self.n*9+1, dtype=INT)
cdef array.array vert_indices = array.clone(int_tmplt, 0, False)
# Flatten regions into array, so it can be used later.
cdef INT_T i
for i in range(self.n*9):
verts = self.scipy_vor.regions[self.scipy_vor.point_region[i]]
offsets[i+1] = offsets[i] + len(verts) # Build offsets.
vert_indices.extend(array.array('q', verts)) # Flatten
# Get vertices of original N sites.
cdef np.ndarray[INT_T, ndim=1] vert_indices_np = np.asarray(vert_indices)
cdef np.ndarray[INT_T, ndim=1] border_sites = np.unique(np.searchsorted(
np.asarray(offsets), # Check indices where below matches would be inserted
np.nonzero(np.isin( # Indices of other verts being in bound verts.
vert_indices_np[offsets[self.n]:], # Rest of the verts to check.
np.unique(vert_indices_np[:offsets[self.n]]) # Bound verts
))[0] + offsets[self.n],
side='right' # If on index == offset_number, should be part of the next site.
) - 1) # Subtract by one to get actual site number.
cdef INT_T border_num = len(border_sites)
# Build sites array.
# [Site Index, Edge Index/Offset, Edge Count]
self.sites = np.empty((self.n+border_num, 3), dtype=INT)
self.sites.base[:self.n, 0] = np.arange(self.n, dtype=INT)
self.sites.base[self.n:, 0] = border_sites
self.sites.base[:self.n+1, 1] = offsets[:self.n+1]
for i in range(self.n):
self.sites[i, 2] = self.sites[i+1, 1] - self.sites[i, 1]
cdef INT_T edge_count = offsets[self.n]
cdef INT_T diff
for i in range(border_num):
diff = offsets[border_sites[i]+1] - offsets[border_sites[i]]
edge_count += diff
self.sites[self.n+i, 2] = diff
if i < border_num - 1:
self.sites[self.n+i+1, 1] = self.sites[self.n+i, 1] + diff
# Build edges array
# [Origin Index, Site Index, Next Index, Prev Index, Twin Index]
self.edges = np.empty((edge_count, 5), dtype=INT)
cdef np.ndarray[INT_T, ndim=1] site_verts
cdef INT_T j, site_i, edge_i, edge_offset, vert_num, twin_index, prev_res
edge_indices = dict()
for i in range(self.n + border_num):
site_i = self.sites[i, 0]
edge_offset = self.sites[i, 1]
site_verts = vert_indices_np[offsets[site_i]:offsets[site_i+1]]
# Scipy outputs sorted vertices, but reverse if not counterclockwise.
if not VoronoiContainer.sign(self.points[site_i],
self.vertices[site_verts[0]], self.vertices[site_verts[1]]):
site_verts = np.flip(site_verts)
vert_num = offsets[site_i+1] - offsets[site_i]
for j in range(vert_num):
edge_i = edge_offset+j
self.edges[edge_i, 0] = site_verts[j]
self.edges[edge_i, 1] = i
# Add vert_num because of C modulo to get always positive.
self.edges[edge_i, 2] = (j+vert_num+1) % vert_num + edge_offset
self.edges[edge_i, 3] = (j+vert_num-1) % vert_num + edge_offset
# Get reversed tuple to theck for twin.
twin_index = edge_indices.get(
(site_verts[(j+1) % vert_num], site_verts[j]
), -1)
self.edges[edge_i, 4] = twin_index
if twin_index == -1:
edge_indices[(site_verts[j], site_verts[(j+1) % vert_num])] = \
j + edge_offset
else:
self.edges[twin_index, 4] = j + edge_offset
self.site_cache = np.empty((self.n + border_num, 5), dtype=FLOAT)
self.edge_cache = np.empty((edge_count, self.edge_cache_map.size), dtype=FLOAT)
cdef void common_cache(VoronoiContainer self) except *:
cdef VoronoiInfo info = init.VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
cdef Site xi
cdef HalfEdge em, ep
cdef Vector2D p, q, la, da, Rla
cdef FLOAT_T [:] area = np.zeros(self.sites.shape[0], dtype=FLOAT)
cdef FLOAT_T [:] perim = np.zeros(self.sites.shape[0], dtype=FLOAT)
cdef INT_T i, j
cdef FLOAT_T e_area, la_mag
for i in prange(self.sites.shape[0], nogil=True):
xi = init.Site(i, &info)
em = xi.edge(&xi)
for j in range(xi.edge_num(&xi)):
ep = em.next(&em)
p, q = em.origin(&em), ep.origin(&ep)
la, da = q.copy.vsub(&q, p), p.copy.vsub(&p, xi.vec(&xi)) # vp - vm, vm - xi
la_mag = la.mag(&la)
e_area = la.dot(&la, da.rot(&da))
Rla = la.rot(&la)
em.cache.la(&em, la)
em.cache.la_mag(&em, la_mag)
em.cache.da(&em, da)
em.cache.da_mag(&em, da.mag(&da))
em.cache.xij(&em, Rla.copy.smul(&Rla, -e_area/la.dot(&la, la)))
if info.edge_cache_map.iF != -1:
em.cache.F(&em, e_area)
area[i] += e_area
perim[i] += la_mag
em = em.next(&em)
xi.cache.area(&xi, area[i]/2)
xi.cache.perim(&xi, perim[i])
xi.cache.isoparam(&xi, 2*PI*area[i]/(perim[i]*perim[i]))
@staticmethod
cdef inline Matrix2x2 calc_H(HalfEdge em, HalfEdge ep) nogil:
cdef Vector2D xmv, xpv, im, mp, right, Rpm, Rim, f
cdef Matrix2x2 h
cdef FLOAT_T im2, mp2
# Vectors from xi to xm and xp.
xmv, xpv = em.cache.xij(&em, NAN_VECTOR), ep.cache.xij(&ep, NAN_VECTOR)
im, mp = xmv.copy.neg(&xmv), xmv.copy.vsub(&xmv, xpv) # -xmv, xmv - xpv
im2, mp2 = -(xmv.dot(&xmv, xmv)), xmv.dot(&xmv, xmv) - xpv.dot(&xpv, xpv)
# (-xmv*xmv, xmv*xmv - xpv*xpv)
right = init.Vector2D(im2, mp2)
Rpm, Rim = R.vecmul(&R, mp.copy.neg(&mp)), im.rot(&im) # R*-mp, R*im
h = init.Matrix2x2(Rpm.x, Rim.x, Rpm.y, Rim.y) # [Rpm | Rim], h is temporary.
f = h.vecmul(&h, right) # [Rpm | Rim]*right
h = R.copy.smul(&R, mp2*(2*mp.dot(&mp, Rim))) # fp*g, g is a scalar.
# (fp*g - f*gp)/(g**2). f is a column vector, gp = 2*Rpm is a row vector.
h.self.msub(&h, init.Matrix2x2(
f.x*2*Rpm.x, f.x*2*Rpm.y, f.y*2*Rpm.x, f.y*2*Rpm.y
))
h.self.sdiv(&h, (2*mp.dot(&mp, Rim))**2)
return h
@staticmethod
cdef inline bint sign(FLOAT_T [::1] ref, FLOAT_T [::1] p, FLOAT_T [::1] q):
"""
Outputs if p2 - self is counterclockwise of p1 - self.
:param p1: [List[float]] first vector
:param p2: [List[float]] second vector
:return: [bool] returns if counterclockwise.
"""
return ((q[0] - ref[0])*-(p[1] - ref[1]) + \
(q[1] - ref[1])*(p[0] - ref[0])) >= 0
# global ROT
# cdef np.ndarray[FLOAT_T, ndim=2] rot = np.asarray(ROT).reshape(2,2)
# return (q - ref).dot(rot.dot(p - ref)) >= 0
cdef void precompute(self) except *:
pass
cdef void calc_grad(self) except *:
pass
cdef void get_statistics(self) except *:
self.stats = {}
cache = self.site_cache[:self.n, :]
self.stats["site_areas"] = np.asarray(cache[:, SITE_CACHE_MAP.iarea])
self.stats["site_edge_count"] = np.asarray(self.sites[:self.n, 2])
self.stats["site_isos"] = np.asarray(cache[:, SITE_CACHE_MAP.iisoparam])
self.stats["site_energies"] = np.asarray(cache[:, SITE_CACHE_MAP.ienergy])
self.stats["avg_radius"] = np.asarray(cache[:, SITE_CACHE_MAP.iavg_radius])
self.stats["isoparam_avg"] = self.stats["site_areas"] / \
(PI*self.stats["avg_radius"]**2)
edges = np.asarray(self.edges)
mask = np.nonzero(edges[:, 0] != -1)[0]
all_edges = mask[(mask % 2 == 0)]
caches = edges[all_edges, 4]
edge_cache = np.asarray(self.edge_cache)
self.stats["edge_lengths"] = edge_cache[caches, self.edge_cache_map.ila_mag]
@property
def site_arr(self):
return np.asarray(self.points[:self.n], dtype=FLOAT)
@property
def vor_data(self):
return self.scipy_vor
@property
def gradient(self):
return np.asarray(self.grad, dtype=FLOAT)
def add_sites(self, add):
return (self.site_arr + add) % np.asarray(self.dim, dtype=FLOAT)
def iterate(self, FLOAT_T step):
k1 = self.gradient
k2 = self.__class__(self.n, self.w, self.h, self.r,
self.add_sites(step*k1)
).gradient
return (step/2)*(k1+k2), k1
def hessian(self, d: float) -> np.ndarray:
"""
Obtains the approximate Hessian.
:param d: [float] small d for approximation.
:return: 2Nx2N array that represents Hessian.
"""
HE = np.zeros((2*self.n, 2*self.n))
new_sites = np.copy(self.site_arr) # Maintain one copy for speed.
for i in range(self.n):
for j in range(2):
mod = self.w if j == 0 else self.h
new_sites[i][j] = (new_sites[i][j] + d) % mod
Ep = self.__class__(self.n, self.w, self.h, self.r, new_sites)
new_sites[i][j] = (new_sites[i][j] - 2*d) % mod
Em = self.__class__(self.n, self.w, self.h, self.r, new_sites)
new_sites[i][j] = (new_sites[i][j] + d) % mod
HE[:, 2*i+j] = ((Ep.gradient - Em.gradient)/(2*d)).flatten()
# Average out discrepencies, since it should be symmetric.
for i in range(2*self.n):
for j in range(i, 2*self.n):
HE[i][j] = (HE[i][j] + HE[j][i])/2
HE[j][i] = HE[i][j]
return HE

3
squish/common.py
View File

@ -3,7 +3,8 @@ from typing import List, Tuple, Union, Optional, Iterator, Generator
import pickle, numpy as np
from math import gcd
from pathlib import Path
from ._squish import VoronoiContainer, AreaEnergy, RadialALEnergy, RadialTEnergy
from squish.voronoi import VoronoiContainer
from squish.energy import AreaEnergy, RadialALEnergy, RadialTEnergy
OUTPUT_DIR = Path("squish_output")

8834
squish/core.c Normal file
View File

File diff suppressed because it is too large Load Diff

135
squish/core.pxd Normal file
View File

@ -0,0 +1,135 @@
cimport numpy as np
# Cython Types.
ctypedef np.int64_t INT_T
ctypedef np.float64_t FLOAT_T
ctypedef (INT_T, INT_T) Pair
# Integer Array psuedo-class for continguous arrays.
cdef struct IArray:
INT_T* arr
Pair shape
INT_T (*get)(IArray*, Pair) nogil
void (*set)(IArray*, Pair, INT_T) nogil
# Float Array psuedo-class for continguous arrays.
ctypedef struct FArray:
FLOAT_T* arr
Pair shape
FLOAT_T (*get)(FArray*, Pair) nogil
void (*set)(FArray*, Pair, FLOAT_T) nogil
# Simple append-only dynamic integer array.
# ctypedef struct IList:
# INT_T* data
# INT_T size, length
# void (*append)(IList*, INT_T) nogil
# void (*free)(IList*) nogil
# Uses an array of bits to determine if value in set.
ctypedef struct BitSet:
INT_T* bits
bint (*add)(BitSet*, INT_T) nogil
void (*free)(BitSet*) nogil
# Psuedo-operator definitions.
ctypedef Vector2D* (*VectorSelfVecOp)(Vector2D*, Vector2D) nogil
ctypedef Vector2D (*VectorCopyVecOp)(Vector2D*, Vector2D) nogil
ctypedef Vector2D* (*VectorSelfSclOp)(Vector2D*, FLOAT_T) nogil
ctypedef Vector2D (*VectorCopySclOp)(Vector2D*, FLOAT_T) nogil
ctypedef Matrix2x2* (*MatrixSelfMatOp)(Matrix2x2*, Matrix2x2) nogil
ctypedef Matrix2x2 (*MatrixCopyMatOp)(Matrix2x2*, Matrix2x2) nogil
ctypedef Matrix2x2* (*MatrixSelfSclOp)(Matrix2x2*, FLOAT_T) nogil
ctypedef Matrix2x2 (*MatrixCopySclOp)(Matrix2x2*, FLOAT_T) nogil
ctypedef struct VectorSelfOps:
Vector2D* (*neg)(Vector2D*) nogil
VectorSelfVecOp vadd
VectorSelfVecOp vsub
VectorSelfVecOp vmul
VectorSelfVecOp vdiv
Vector2D* (*matmul)(Vector2D*, Matrix2x2) nogil
VectorSelfSclOp sadd
VectorSelfSclOp ssub
VectorSelfSclOp smul
VectorSelfSclOp sdiv
ctypedef struct VectorCopyOps:
Vector2D (*neg)(Vector2D*) nogil
VectorCopyVecOp vadd
VectorCopyVecOp vsub
VectorCopyVecOp vmul
VectorCopyVecOp vdiv
Vector2D (*matmul)(Vector2D*, Matrix2x2) nogil
VectorCopySclOp sadd
VectorCopySclOp ssub
VectorCopySclOp smul
VectorCopySclOp sdiv
ctypedef struct MatrixSelfOps:
Matrix2x2* (*neg)(Matrix2x2*) nogil
MatrixSelfMatOp madd
MatrixSelfMatOp msub
MatrixSelfMatOp mmul
MatrixSelfMatOp mdiv
MatrixSelfMatOp matmul
MatrixSelfSclOp sadd
MatrixSelfSclOp ssub
MatrixSelfSclOp smul
MatrixSelfSclOp sdiv
ctypedef struct MatrixCopyOps:
Matrix2x2 (*neg)(Matrix2x2*) nogil
MatrixCopyMatOp madd
MatrixCopyMatOp msub
MatrixCopyMatOp mmul
MatrixCopyMatOp mdiv
MatrixCopyMatOp matmul
MatrixCopySclOp sadd
MatrixCopySclOp ssub
MatrixCopySclOp smul
MatrixCopySclOp sdiv
# Psuedo-class for a 2-dimensional vector. No orientation.
ctypedef struct Vector2D:
FLOAT_T x, y
VectorSelfOps self
VectorCopyOps copy
bint (*equals)(Vector2D*, Vector2D) nogil
Vector2D (*rot)(Vector2D*) nogil
FLOAT_T (*dot)(Vector2D*, Vector2D) nogil
FLOAT_T (*mag)(Vector2D*) nogil
# Psuedo-class for a 2x2 matrix.
ctypedef struct Matrix2x2:
FLOAT_T a, b, c, d
MatrixSelfOps self
MatrixCopyOps copy
bint (*equals)(Matrix2x2*, Matrix2x2) nogil
Vector2D (*vecmul)(Matrix2x2*, Vector2D) nogil
cdef IArray _IArray(INT_T*, Pair) nogil
cdef FArray _FArray(FLOAT_T*, Pair) nogil
cdef BitSet _BitSet(INT_T) nogil
cdef Vector2D _Vector2D(FLOAT_T, FLOAT_T) nogil
cdef Matrix2x2 _Matrix2x2(FLOAT_T, FLOAT_T, FLOAT_T, FLOAT_T) nogil

363
squish/core.pyx Normal file
View File

@ -0,0 +1,363 @@
from libc.stdlib cimport calloc, free
from libc.math cimport sqrt
from squish.core cimport INT_T, FLOAT_T, Pair, \
IArray, FArray, BitSet, Vector2D, Matrix2x2, \
VectorSelfOps, VectorCopyOps, MatrixSelfOps, MatrixCopyOps
#### Constants ####
# In most cases, the amount of edges relevant to a gradient will
# not exceed this number. However, we assign a growth rate of 8 edges,
# when dynamically allocating.
cdef INT_T EDGE_ARR_SIZE = 32
cdef VectorSelfOps VSO
cdef VectorCopyOps VCO
cdef MatrixSelfOps MSO
cdef MatrixCopyOps MCO
VSO.neg, VSO.vadd, VSO.vsub, VSO.vmul, VSO.vdiv, VSO.sadd, VSO.ssub, VSO.smul, VSO.sdiv = \
v_neg_s, v_vadd_s, v_vsub_s, v_vmul_s, v_vdiv_s, v_sadd_s, v_ssub_s, v_smul_s, v_sdiv_s
VSO.matmul = v_matmul_s
VCO.neg, VCO.vadd, VCO.vsub, VCO.vmul, VCO.vdiv, VCO.sadd, VCO.ssub, VCO.smul, VCO.sdiv = \
v_neg_c, v_vadd_c, v_vsub_c, v_vmul_c, v_vdiv_c, v_sadd_c, v_ssub_c, v_smul_c, v_sdiv_c
VCO.matmul = v_matmul_c
MSO.neg, MSO.madd, MSO.msub, MSO.mmul, MSO.mdiv, MSO.sadd, MSO.ssub, MSO.smul, MSO.sdiv = \
m_neg_s, m_madd_s, m_msub_s, m_mmul_s, m_mdiv_s, m_sadd_s, m_ssub_s, m_smul_s, m_sdiv_s
MSO.matmul = m_matmul_s
MCO.neg, MCO.madd, MCO.msub, MCO.mmul, MCO.mdiv, MCO.sadd, MCO.ssub, MCO.smul, MCO.sdiv = \
m_neg_c, m_madd_c, m_msub_c, m_mmul_c, m_mdiv_c, m_sadd_c, m_ssub_c, m_smul_c, m_sdiv_c
MCO.matmul = m_matmul_c
"""
If bound checking is desired, uncomment out ..._valid_indices functions.
"""
#### IArray Methods ####
cdef inline IArray _IArray(INT_T* arr, Pair shape) nogil:
cdef IArray iarray
iarray.arr, iarray.shape = arr, shape
iarray.get, iarray.set = &iarray_get, &iarray_set
return iarray
cdef inline bint iarray_valid_indices(IArray* self, Pair index) nogil:
if index[0] > self.shape[0] or index[1] > self.shape[1]:
with gil:
raise IndexError(f"Index out of range for IArray with shape {self.shape}")
cdef inline INT_T iarray_get(IArray* self, Pair index) nogil:
#iarray_valid_indices(&self, index)
return self.arr[index[0]*self.shape[1] + index[1]]
cdef inline void iarray_set(IArray* self, Pair index, INT_T val) nogil:
#iarray_valid_indices(&self, index)
self.arr[index[0]*self.shape[1] + index[1]] = val
#### FArray Methods ####
cdef inline FArray _FArray(FLOAT_T* arr, Pair shape) nogil:
cdef FArray farray
farray.arr, farray.shape = arr, shape
farray.get, farray.set = &farray_get, &farray_set
return farray
cdef inline bint farray_valid_indices(FArray* self, Pair index) nogil:
if index[0] > self.shape[0] or index[1] > self.shape[1]:
with gil:
raise IndexError(f"Index out of range for FArray with shape {self.shape}")
cdef inline FLOAT_T farray_get(FArray* self, Pair index) nogil:
#iarray_valid_indices(&self, index)
return self.arr[index[0]*self.shape[1] + index[1]]
cdef inline void farray_set(FArray* self, Pair index, FLOAT_T val) nogil:
#iarray_valid_indices(&self, index)
self.arr[index[0]*self.shape[1] + index[1]] = val
#### IList Methods ####
# cdef inline IList init_ilist() nogil:
# cdef IList ilist
# ilist.size = EDGE_ARR_SIZE
# ilist.length = 0
# ilist.data = <INT_T*> malloc(self.size * sizeof(INT_T))
# ilist.append, ilist.free = ilist_append, ilist_free
# return ilist
# cdef inline void ilist_append(IList* self, INT_T) nogil:
# if self.size == self.length:
# ilist.data = <INT_T*> realloc((self.size+8) * sizeof(INT_T))
# self.size += 8
# self.data[self.length] == INT_T
# self.length += 1
# cdef inline void ilist_free(IList* self) nogil:
# free(self.data)
#### BitSet Methods ####
cdef inline BitSet _BitSet(INT_T elements) nogil:
cdef BitSet bitset
bitset.bits = <INT_T*> calloc(((elements/sizeof(INT_T))+1), sizeof(INT_T))
bitset.add, bitset.free = &bitset_add, &bitset_free
return bitset
cdef inline bint bitset_add(BitSet* self, INT_T val) nogil:
cdef INT_T index, rel_index, old
index = val/sizeof(INT_T)
old = self.bits[index]
rel_index = val - index*sizeof(INT_T)
self.bits[index] = (1 << rel_index) | old # New value.
return old == self.bits[index] # Means 1 was already there.
cdef inline void bitset_free(BitSet* self) nogil:
free(self.bits)
#### Vector2D Methods ####
"""
Prefix 'v' stands for vector, element by element operation.
Prefix 's' stands for scalar, broadcasted operation.
Suffix 'w' stands for write, overwriting current value.
Suffix 'n' stands for new, copying to a new location.
While it's possible to chain 'new' operations, when possible,
avoid this, so fewer objects are needed.
"""
cdef inline Vector2D _Vector2D(FLOAT_T x, FLOAT_T y) nogil:
cdef Vector2D vec
vec.x, vec.y = x, y
vec.self, vec.copy = VSO, VCO
vec.equals, vec.rot, vec.dot, vec.mag = &v_equals, &rot, &dot, &mag
return vec
cdef inline bint v_equals(Vector2D* self, Vector2D w) nogil:
return ((self.x == w.x) and (self.y == w.y))
cdef inline Vector2D* v_neg_s(Vector2D* self) nogil:
self.x = -self.x
self.y = -self.y
return self
cdef inline Vector2D* v_vadd_s(Vector2D* self, Vector2D w) nogil:
self.x += w.x
self.y += w.y
return self
cdef inline Vector2D* v_vsub_s(Vector2D* self, Vector2D w) nogil:
self.x -= w.x
self.y -= w.y
return self
cdef inline Vector2D* v_vmul_s(Vector2D* self, Vector2D w) nogil:
self.x *= w.x
self.y *= w.y
return self
cdef inline Vector2D* v_vdiv_s(Vector2D* self, Vector2D w) nogil:
self.x /= w.x
self.y /= w.y
return self
cdef inline Vector2D* v_sadd_s(Vector2D* self, FLOAT_T s) nogil:
self.x += s
self.y += s
return self
cdef inline Vector2D* v_ssub_s(Vector2D* self, FLOAT_T s) nogil:
self.x -= s
self.y -= s
return self
cdef inline Vector2D* v_smul_s(Vector2D* self, FLOAT_T s) nogil:
self.x *= s
self.y *= s
return self
cdef inline Vector2D* v_sdiv_s(Vector2D* self, FLOAT_T s) nogil:
self.x /= s
self.y /= s
return self
cdef inline Vector2D* v_matmul_s(Vector2D* self, Matrix2x2 m) nogil:
self.x, self.y = self.x*m.a + self.y*m.c, self.x*m.b + self.y*m.d
return self
cdef inline Vector2D v_neg_c(Vector2D* self) nogil:
return _Vector2D(-self.x, -self.y)
cdef inline Vector2D v_vadd_c(Vector2D* self, Vector2D w) nogil:
return _Vector2D(self.x + w.x, self.y + w.y)
cdef inline Vector2D v_vsub_c(Vector2D* self, Vector2D w) nogil:
return _Vector2D(self.x - w.x, self.y - w.y)
cdef inline Vector2D v_vmul_c(Vector2D* self, Vector2D w) nogil:
return _Vector2D(self.x * w.x, self.y * w.y)
cdef inline Vector2D v_vdiv_c(Vector2D* self, Vector2D w) nogil:
return _Vector2D(self.x / w.x, self.y / w.y)
cdef inline Vector2D v_sadd_c(Vector2D* self, FLOAT_T s) nogil:
return _Vector2D(self.x + s, self.y + s)
cdef inline Vector2D v_ssub_c(Vector2D* self, FLOAT_T s) nogil:
return _Vector2D(self.x + s, self.y + s)
cdef inline Vector2D v_smul_c(Vector2D* self, FLOAT_T s) nogil:
return _Vector2D(self.x * s, self.y * s)
cdef inline Vector2D v_sdiv_c(Vector2D* self, FLOAT_T s) nogil:
return _Vector2D(self.x / s, self.y / s)
cdef inline Vector2D v_matmul_c(Vector2D* self, Matrix2x2 m) nogil:
return _Vector2D(
self.x*m.a + self.y*m.c, self.x*m.b + self.y*m.d
)
cdef inline Vector2D rot(Vector2D* self) nogil:
return _Vector2D(-self.y, self.x)
cdef inline FLOAT_T dot(Vector2D* self, Vector2D w) nogil:
return self.x*w.x + self.y*w.y
cdef inline FLOAT_T mag(Vector2D* self) nogil:
return <FLOAT_T>sqrt(<double>(self.x*self.x + self.y*self.y))
#### Matrix2x2 Methods ####
cdef inline Matrix2x2 _Matrix2x2(FLOAT_T a, FLOAT_T b, FLOAT_T c, FLOAT_T d) nogil:
cdef Matrix2x2 matrix
matrix.a, matrix.b, matrix.c, matrix.d = a, b, c, d
matrix.self, matrix.copy = MSO, MCO
matrix.equals, matrix.vecmul = &m_equals, &m_vecmul
return matrix
cdef inline bint m_equals(Matrix2x2* self, Matrix2x2 m) nogil:
return (
(self.a == m.a) and (self.b == m.b) and (self.c == m.c) and (self.d == m.d)
)
cdef inline Vector2D m_vecmul(Matrix2x2* self, Vector2D v) nogil:
return _Vector2D(
self.a*v.x + self.b*v.y, self.c*v.x + self.d*v.y
)
cdef inline Matrix2x2* m_neg_s(Matrix2x2* self) nogil:
self.a, self.b, self.c, self.d = -self.a, -self.b, -self.c, -self.d
return self
cdef inline Matrix2x2* m_madd_s(Matrix2x2* self, Matrix2x2 m) nogil:
self.a += m.a
self.b += m.b
self.c += m.c
self.d += m.d
return self
cdef inline Matrix2x2* m_msub_s(Matrix2x2* self, Matrix2x2 m) nogil:
self.a -= m.a
self.b -= m.b
self.c -= m.c
self.d -= m.d
return self
cdef inline Matrix2x2* m_mmul_s(Matrix2x2* self, Matrix2x2 m) nogil:
self.a *= m.a
self.b *= m.b
self.c *= m.c
self.d *= m.d
return self
cdef inline Matrix2x2* m_mdiv_s(Matrix2x2* self, Matrix2x2 m) nogil:
self.a /= m.a
self.b /= m.b
self.c /= m.c
self.d /= m.d
return self
cdef inline Matrix2x2* m_sadd_s(Matrix2x2* self, FLOAT_T s) nogil:
self.a += s
self.b += s
self.c += s
self.d += s
return self
cdef inline Matrix2x2* m_ssub_s(Matrix2x2* self, FLOAT_T s) nogil:
self.a -= s
self.b -= s
self.c -= s
self.d -= s
return self
cdef inline Matrix2x2* m_smul_s(Matrix2x2* self, FLOAT_T s) nogil:
self.a *= s
self.b *= s
self.c *= s
self.d *= s
return self
cdef inline Matrix2x2* m_sdiv_s(Matrix2x2* self, FLOAT_T s) nogil:
self.a /= s
self.b /= s
self.c /= s
self.d /= s
return self
cdef inline Matrix2x2* m_matmul_s(Matrix2x2* self, Matrix2x2 m) nogil:
self.a, self.b, self.c, self.d = \
self.a*m.a + self.b*m.c, self.a*m.b + self.b*m.d, \
self.c*m.a + self.d*m.c, self.c*m.b + self.d*m.d
return self
cdef inline Matrix2x2 m_neg_c(Matrix2x2* self) nogil:
return _Matrix2x2(-self.a, -self.b, -self.c, -self.d)
cdef inline Matrix2x2 m_madd_c(Matrix2x2* self, Matrix2x2 m) nogil:
return _Matrix2x2(self.a+m.a, self.b+m.b, self.c+m.c, self.d+m.d)
cdef inline Matrix2x2 m_msub_c(Matrix2x2* self, Matrix2x2 m) nogil:
return _Matrix2x2(self.a-m.a, self.b-m.b, self.c-m.c, self.d-m.d)
cdef inline Matrix2x2 m_mmul_c(Matrix2x2* self, Matrix2x2 m) nogil:
return _Matrix2x2(self.a*m.a, self.b*m.b, self.c*m.c, self.d*m.d)
cdef inline Matrix2x2 m_mdiv_c(Matrix2x2* self, Matrix2x2 m) nogil:
return _Matrix2x2(self.a/m.a, self.b/m.b, self.c/m.c, self.d/m.d)
cdef inline Matrix2x2 m_sadd_c(Matrix2x2* self, FLOAT_T s) nogil:
return _Matrix2x2(self.a+s, self.b+s, self.c+s, self.d+s)
cdef inline Matrix2x2 m_ssub_c(Matrix2x2* self, FLOAT_T s) nogil:
return _Matrix2x2(self.a-s, self.b-s, self.c-s, self.d-s)
cdef inline Matrix2x2 m_smul_c(Matrix2x2* self, FLOAT_T s) nogil:
return _Matrix2x2(self.a*s, self.b*s, self.c*s, self.d*s)
cdef inline Matrix2x2 m_sdiv_c(Matrix2x2* self, FLOAT_T s) nogil:
return _Matrix2x2(self.a/s, self.b/s, self.c/s, self.d/s)
cdef inline Matrix2x2 m_matmul_c(Matrix2x2* self, Matrix2x2 m) nogil:
return _Matrix2x2(
self.a*m.a + self.b*m.c, self.a*m.b + self.b*m.d,
self.c*m.a + self.d*m.c, self.c*m.b + self.d*m.d
)

4
squish/diagram.py
View File

@ -183,8 +183,8 @@ class Diagram:
defects[7]["x"].append(vec[0])
defects[7]["y"].append(vec[1])
ax.scatter(defects[5]["x"], defects[5]["y"], marker="p", color="red")
ax.scatter(defects[7]["x"], defects[7]["y"], marker="*", color="red")
ax.scatter(defects[5]["x"], defects[5]["y"], marker="p", color="C0")
ax.scatter(defects[7]["x"], defects[7]["y"], marker="*", color="C0")
ax.text(
0.05,

27737
squish/energy.c Normal file
View File

File diff suppressed because it is too large Load Diff

302
squish/energy.pyx Normal file
View File

@ -0,0 +1,302 @@
import numpy as np
from cython.parallel import parallel, prange
cimport numpy as np
from libc.math cimport NAN, pi as PI, atanh
from squish.core cimport INT_T, FLOAT_T, Vector2D, Matrix2x2, BitSet, \
_Vector2D, _Matrix2x2, _BitSet
from squish.voronoi cimport Site, HalfEdge, EdgeCacheMap, VoronoiInfo, \
_Site, _HalfEdge, _EdgeCacheMap, _VoronoiInfo, VoronoiContainer, \
R, NAN_MATRIX, NAN_VECTOR
#### Constants ####
INT = np.int64
FLOAT = np.float64
cdef EdgeCacheMap AREA_ECM = _EdgeCacheMap(0, 4, 6, 8, 10, 12, 13, -1, -1, 14)
cdef EdgeCacheMap RADIALT_ECM = _EdgeCacheMap(0, 4, 6, 8, -1, 10, 11, 12, 13, 14)
cdef class AreaEnergy(VoronoiContainer):
"""
Class for formulas relevant to the Area energy.
:param n: [int] how many sites to generate.
:param w: [float] width of the bounding domain.
:param h: [float] height of the bounding domain.
:param r: [float] radius of zero energy circle.
:param sites: [np.ndarray] collection of sites.
"""
attr_str = "area"
title_str = "Area"
def __init__(AreaEnergy self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r,
np.ndarray[FLOAT_T, ndim=2] site_arr):
self.edge_cache_map = &AREA_ECM
self.energy = 0.0
super().__init__(n, w, h, r, site_arr)
self.minimum = (<FLOAT_T>n)*(w*h/(<FLOAT_T>n)-PI*r**2)**2
cdef void precompute(self) except *:
cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
cdef Site xi
cdef HalfEdge em, e, ep
cdef Vector2D vdiff
cdef FLOAT_T [:] site_energy = np.full(self.sites.shape[0], PI*self.r**2)
cdef INT_T i, j
for i in prange(self.sites.shape[0], nogil=True):
xi = _Site(i, &info)
e = xi.edge(&xi)
site_energy[i] = (xi.cache.area(&xi, NAN) - site_energy[i])**2
xi.cache.energy(&xi, site_energy[i])
j = 0
while j < xi.edge_num(&xi):
em, ep = e.prev(&e), e.next(&e)
vdiff = em.origin(&em)
vdiff.self.vsub(&vdiff, ep.origin(&ep))
e.cache.dVdv(&e, R.vecmul(&R, vdiff))
e.cache.H(&e, VoronoiContainer.calc_H(em, e))
e = e.next(&e)
j = j + 1
self.energy = np.sum(site_energy[:self.n])
cdef void calc_grad(self) except *:
cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
cdef Site xi, xf
cdef HalfEdge e, f
cdef Vector2D dedxi_p
cdef BitSet edge_set
cdef INT_T num_edges = self.edges.shape[0]
cdef FLOAT_T A = PI*self.r**2
cdef FLOAT_T [:, ::1] dedx = np.zeros((self.n, 2), dtype=FLOAT)
cdef INT_T i, j
for i in prange(self.n, nogil=True):
xi = _Site(i, &info)
e = xi.edge(&xi)
edge_set = _BitSet(num_edges)
j = 0
while j < xi.edge_num(&xi): # Looping through site edges.
f = e
while True: # Circling this vertex.
if not edge_set.add(&edge_set, f.arr_index):
xf = f.face(&f)
dedxi_p = f.cache.dVdv(&f, NAN_VECTOR) #dVdv
dedxi_p.self.smul(&dedxi_p, xf.cache.area(&xf, NAN) - A)
dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX))
dedx[i][0] -= dedxi_p.x
dedx[i][1] -= dedxi_p.y
f = f.twin(&f)
f = f.next(&f)
if f.arr_index == e.arr_index:
break
e = e.next(&e)
j = j + 1
edge_set.free(&edge_set)
self.grad = dedx
cdef class RadialALEnergy(VoronoiContainer):
"""
Class for formulas relevant to the Area energy.
:param n: [int] how many sites to generate.
:param w: [float] width of the bounding domain.
:param h: [float] height of the bounding domain.
:param r: [float] radius of zero energy circle.
:param sites: [np.ndarray] collection of sites.
"""
attr_str = "radial-al"
title_str = "Radial[AL]"
def __init__(AreaEnergy self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r,
np.ndarray[FLOAT_T, ndim=2] site_arr):
#self.edge_cache_map = &AREA_EDGE_CACHE_MAP
self.energy = 0.0
super().__init__(n, w, h, r, site_arr)
cdef void precompute(self) except *:
cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
pass
cdef void calc_grad(self) except *:
cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
pass
cdef class RadialTEnergy(VoronoiContainer):
"""
Class for formulas relevant to the Area energy.
:param n: [int] how many sites to generate.
:param w: [float] width of the bounding domain.
:param h: [float] height of the bounding domain.
:param r: [float] radius of zero energy circle.
:param sites: [np.ndarray] collection of sites.
"""
attr_str = "radial-t"
title_str = "Radial[T]"
def __init__(AreaEnergy self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r,
np.ndarray[FLOAT_T, ndim=2] site_arr):
self.edge_cache_map = &RADIALT_ECM
self.energy = 0.0
super().__init__(n, w, h, r, site_arr)
cdef void precompute(self) except *:
cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
cdef Site xi
cdef HalfEdge e, em, ep
cdef Vector2D la
# All energy has a 2pir_0 term.
cdef FLOAT_T [:] site_energy = np.full(self.sites.shape[0], 2*PI*self.r**2)
cdef FLOAT_T [:] avg_radii = np.zeros(self.sites.shape[0])
cdef FLOAT_T sm, sp
cdef INT_T i, j
for i in prange(self.sites.shape[0], nogil=True):
xi = _Site(i, &info)
e = xi.edge(&xi)
j = 0
while j < xi.edge_num(&xi):
em, ep = e.prev(&e), e.next(&e)
e.cache.H(&e, VoronoiContainer.calc_H(em, e))
la = e.cache.la(&e, NAN_VECTOR)
sp = la.dot(&la, ep.cache.da(&ep, NAN_VECTOR)) # dap . la
sm = la.dot(&la, e.cache.da(&e, NAN_VECTOR)) # da . la
sp = sp / (ep.cache.da_mag(&ep, NAN) * e.cache.la_mag(&e, NAN))
sm = sm / (e.cache.da_mag(&e, NAN) * e.cache.la_mag(&e, NAN))
e.cache.calI(&e, <FLOAT_T>(atanh(<double>sp) - atanh(<double>sm)))
avg_radii[i] += (
e.cache.area_p(&e, NAN) * e.cache.calI(&e, NAN)
/ e.cache.la_mag(&e, NAN)
)
e = e.next(&e)
j = j + 1
site_energy[i] += 2*(xi.cache.area(&xi, NAN) - self.r*avg_radii[i])
xi.cache.avg_radius(&xi, avg_radii[i]/(2*PI))
xi.cache.energy(&xi, site_energy[i])
self.energy = np.sum(site_energy[:self.n])
cdef void calc_grad(self) except *:
cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
cdef Site xi
cdef HalfEdge e, fm, f, fp
cdef Vector2D temp1, temp2, temp3, dedxi_p
cdef BitSet edge_set
cdef INT_T num_edges = self.edges.shape[0]
cdef FLOAT_T [:, ::1] dedx = np.zeros((self.n, 2), dtype=FLOAT)
cdef INT_T i, j
for i in prange(self.n, nogil=True):
xi = _Site(i, &info)
e = xi.edge(&xi)
edge_set = _BitSet(num_edges)
j = 0
while j < xi.edge_num(&xi): # Looping through site edges.
f = e
dedxi_p = _Vector2D(0, 0)
# dedx (only x)
temp1 = f.cache.la(&f, NAN_VECTOR)
temp1 = temp1.rot(&temp1)
dedxi_p = temp1.copy.smul(
&temp1,
f.cache.calI(&f, NAN) / f.cache.la_mag(&f, NAN)
)
while True: # Circling this vertex.
fm, fp = f.prev(&f), f.next(&f)
if not edge_set.add(&edge_set, f.arr_index):
# ( -rot(dap) ) / |la|
temp1 = fp.cache.da(&fp, NAN_VECTOR)
temp1 = temp1.rot(&temp1)
temp1.self.sdiv(&temp1, -f.cache.la_mag(&f, NAN))
# la * area_p / |la|^2
temp3 = f.cache.la(&f, NAN_VECTOR)
temp3.self.smul(
&temp3,
f.cache.area_p(&f, NAN) / f.cache.la_mag(&f, NAN)**3
)
# Combine * calI
temp1.self.vadd(&temp1, temp3)
temp1.self.smul(&temp1, f.cache.calI(&f, NAN))
# rot(dam) / |lam|
temp2 = fm.cache.da(&fm, NAN_VECTOR)
temp2 = temp2.rot(&temp2)
temp2.self.sdiv(&temp2, fm.cache.la_mag(&fm, NAN))
# lam * area_pm / |lam|^2
temp3 = fm.cache.la(&fm, NAN_VECTOR)
temp3.self.smul(
&temp3,
fm.cache.area_p(&fm, NAN) / fm.cache.la_mag(&fm, NAN)**3
)
# Combine * calIm
temp2.self.vsub(&temp2, temp3)
temp2.self.smul(&temp2, fm.cache.calI(&fm, NAN))
temp1.self.vadd(&temp1, temp2)
temp1.self.matmul(&temp1, e.get_H(&e, xi))
dedxi_p.self.vadd(&dedxi_p, temp1)
f = f.twin(&f)
f = f.next(&f)
if f.arr_index == e.arr_index:
break
dedx[i][0] -= -2*self.r*dedxi_p.x
dedx[i][1] -= -2*self.r*dedxi_p.y
e = e.next(&e)
j = j + 1
edge_set.free(&edge_set)
self.grad = dedx

58
squish/ordered.py
View File

@ -8,28 +8,22 @@ Config = Tuple[int, int]
def configurations(domain: DomainParams) -> List[Config]:
n, w, h = domain.n, domain.w, domain.h
valid = []
mults = np.arange(n)
configs = np.dstack((np.repeat(mults, n).T, np.tile(mults, n).T))[0][1:]
for i in range(len(configs)):
eq_x = n if configs[i][0] == 0 else configs[i][0]
eq_y = n if configs[i][1] == 0 else configs[i][1]
n = domain.n
coprimes, valid = [], []
for i in range(n):
for j in range(i):
if gcd(i, j) == 1:
coprimes.append((j, i))
if gcd(eq_x, eq_y) != 1:
continue
vecs = (
configs[i]
* np.dstack((w * mults, h * mults)).swapaxes(0, 1)
/ n
% domain.dim
)
vmod2 = np.squeeze(np.matmul(vecs, vecs.transpose(0, 2, 1)))
vmodv = np.squeeze(vecs).dot(vecs[1].T).T.flatten()
if np.all(vmod2 >= vmodv):
valid.append(tuple(configs[i]))
coprimes = set(coprimes)
while len(coprimes) > 0:
first = coprimes.pop()
valid.append(first)
for i in range(2, n):
try:
coprimes.remove(((first[0] * i) % n, (first[1] * i) % n))
except KeyError:
pass
return valid
@ -38,11 +32,15 @@ def get_config_generators(
domain: DomainParams, config: Config
) -> Tuple[Config, Config]:
n, w, h = domain.n, domain.w, domain.h
q1 = sites(domain, config)
v = q1[1]
# Sites to check can ignore 0*v and v itself.
all_sites = np.concatenate((q1, q1 - [w, 0], q1 - [w, h], q1 - [0, h]))[2:]
q1 = sites(domain, config)[1:]
all_sites = np.concatenate(
(q1, q1 - np.array([w, 0]), q1 - np.array([0, h]), q1 - domain.dim)
)
# Sort sites by magnitude and smallest.
all_sites = sorted(list(all_sites), key=lambda x: np.linalg.norm(x))
v = all_sites[0] # Smallest vector set to v.
all_sites = np.array(all_sites[1:]) # Remove v from search set.
# Checking 0 < ax + by < v*v to make the sites are within the region.
tol = 1e-3
vdot = np.matmul(all_sites, v)
@ -51,11 +49,13 @@ def get_config_generators(
0, 1
) # Used for the next step, getting site*site
w = in_box[
v2 = in_box[
np.argmin(np.squeeze(np.matmul(in_box, in_box.transpose(0, 2, 1))))
].flatten()
return tuple(v), tuple(w)
if np.all(v == v2):
print(v, v2, n, w, h, config)
return tuple(v), tuple(v2)
def sites(domain: DomainParams, config: Config) -> numpy.ndarray:
@ -92,8 +92,10 @@ def avg_rp(d: float, l: float) -> float:
return (d / (4 * pi)) * log(tan(0.5 * (atan(l / d) + pi / 2)) ** 2)
def circumcenter(v: numpy.ndarray, w: numpy.ndarray) -> Config:
def circumcenter(v: numpy.ndarray, w: numpy.ndarray) -> numpy.ndarray:
det = 1 / (2 * rot(v).dot(w))
if rot(v).dot(w) == 0:
print(v, w)
v2, w2 = v.dot(v), w.dot(w)
c = np.empty((2,))
c[0], c[1] = w[1] * v2 - v[1] * w2, -w[0] * v2 + v[0] * w2

18
squish/simulation.py
View File

@ -2,7 +2,7 @@ from __future__ import annotations
from typing import Optional
import pickle, numpy as np
from math import log10
from math import log10, sqrt
from scipy.linalg import null_space
from timeit import default_timer as timer
from pathlib import Path
@ -57,6 +57,22 @@ class Simulation:
self.frames.append(self.energy.mode(*self.domain, points % self.domain.dim))
def normalize(self) -> None:
new_frames = []
for frame in self.frames:
aspect = frame.w / frame.h
new_domain = DomainParams(
frame.n, sqrt(frame.n * aspect), sqrt(frame.n / aspect), frame.r
)
new_points = frame.site_arr * np.array(
[new_domain.w / frame.w, new_domain.h / frame.h]
)
new_frames.append(self.energy.mode(*new_domain, new_points))
self.frames = new_frames
def get_distinct(self) -> List[int]:
"""Gets the distinct configurations based on the average radii of the sites.
and returns the number of configurations for each distinct configuration.

21762
squish/_squish/_squish.c → squish/voronoi.c
View File

File diff suppressed because it is too large Load Diff

99
squish/voronoi.pxd Normal file
View File

@ -0,0 +1,99 @@
cimport numpy as np
from squish.core cimport INT_T, FLOAT_T, IArray, FArray, Vector2D, Matrix2x2
# Psuedo-class that handles caching for sites.
ctypedef struct SiteCacheMap:
INT_T iarea, iperim, iisoparam, ienergy, iavg_radius, icentroid, imaxcenter
FLOAT_T (*area)(Site*, FLOAT_T) nogil
FLOAT_T (*perim)(Site*, FLOAT_T) nogil
FLOAT_T (*isoparam)(Site*, FLOAT_T) nogil
FLOAT_T (*energy)(Site*, FLOAT_T) nogil
FLOAT_T (*avg_radius)(Site*, FLOAT_T) nogil
Vector2D (*centroid)(Site*, Vector2D) nogil
Vector2D (*maxcenter)(Site*, Vector2D) nogil
# Psuedo-class that handles caching for edges.
ctypedef struct EdgeCacheMap:
INT_T iH, ila, ida, ixij, idVdv, ila_mag, ida_mag, iarea_p, icalI, size
Matrix2x2 (*H)(HalfEdge*, Matrix2x2) nogil
Vector2D (*la)(HalfEdge*, Vector2D) nogil
Vector2D (*da)(HalfEdge*, Vector2D) nogil
Vector2D (*xij)(HalfEdge*, Vector2D) nogil
Vector2D (*dVdv)(HalfEdge*, Vector2D) nogil
FLOAT_T (*la_mag)(HalfEdge*, FLOAT_T) nogil
FLOAT_T (*da_mag)(HalfEdge*, FLOAT_T) nogil
FLOAT_T (*area_p)(HalfEdge*, FLOAT_T) nogil
FLOAT_T (*calI)(HalfEdge*, FLOAT_T) nogil
# Psuedo-class to just contain all pertaining info for sites and edges.
ctypedef struct VoronoiInfo:
IArray sites, edges
FArray points, vertices, site_cache, edge_cache
EdgeCacheMap* edge_cache_map
# Psuedo-class for a Site.
ctypedef struct Site:
INT_T arr_index
VoronoiInfo* info
SiteCacheMap* cache
INT_T (*index)(Site*) nogil
Vector2D (*vec)(Site*) nogil
HalfEdge (*edge)(Site*) nogil
INT_T (*edge_num)(Site*) nogil
# Psuedo-class for an HalfEdge.
ctypedef struct HalfEdge:
INT_T orig_arr_index, arr_index
VoronoiInfo* info
EdgeCacheMap* cache
INT_T (*origin_index)(HalfEdge*) nogil
Vector2D (*origin)(HalfEdge*) nogil
Site (*face)(HalfEdge*) nogil
HalfEdge (*next)(HalfEdge*) nogil
HalfEdge (*prev)(HalfEdge*) nogil
HalfEdge (*twin)(HalfEdge*) nogil
Matrix2x2 (*get_H)(HalfEdge*, Site) nogil
cdef class VoronoiContainer:
cdef readonly INT_T n
cdef readonly FLOAT_T w, h, r, energy
cdef FLOAT_T [2] dim
cdef FLOAT_T [:, ::1] points, vertices, site_cache, edge_cache, grad
cdef INT_T [:, ::1] sites, edges
cdef EdgeCacheMap* edge_cache_map
cdef dict __dict__
cdef void calculate_voronoi(VoronoiContainer self,
np.ndarray[FLOAT_T, ndim=2] site_arr) except *
cdef void generate_dcel(VoronoiContainer self) except *
cdef void common_cache(VoronoiContainer self) except *
cdef void precompute(self) except *
cdef void calc_grad(self) except *
cdef void get_statistics(VoronoiContainer self) except *
@staticmethod
cdef inline Matrix2x2 calc_H(HalfEdge, HalfEdge) nogil
@staticmethod
cdef inline bint sign(FLOAT_T [::1], FLOAT_T [::1], FLOAT_T [::1])
cdef SiteCacheMap _SiteCacheMap(INT_T, INT_T, INT_T, INT_T, INT_T, INT_T, INT_T) nogil
cdef EdgeCacheMap _EdgeCacheMap(INT_T, INT_T, INT_T, INT_T, INT_T, INT_T, INT_T,
INT_T, INT_T, INT_T) nogil
cdef VoronoiInfo _VoronoiInfo(INT_T [:, ::1], INT_T[:, ::1], FLOAT_T[:, ::1],
FLOAT_T[:, ::1], FLOAT_T[:, ::1], FLOAT_T[:, ::1],
EdgeCacheMap*) nogil
cdef Site _Site(INT_T, VoronoiInfo*) nogil
cdef HalfEdge _HalfEdge(INT_T, VoronoiInfo*) nogil
cdef Vector2D NAN_VECTOR
cdef Matrix2x2 R, NAN_MATRIX

736
squish/voronoi.pyx Normal file
View File

@ -0,0 +1,736 @@
import array, scipy.spatial, numpy as np
from timeit import default_timer as timer
from cython.parallel import parallel, prange
cimport numpy as np
from cpython cimport array
from libc.math cimport isnan, NAN, pi as PI
from squish.core cimport INT_T, FLOAT_T, \
IArray, FArray, Vector2D, Matrix2x2, \
_IArray, _FArray, _Vector2D, _Matrix2x2
from squish.voronoi cimport SiteCacheMap, EdgeCacheMap, VoronoiInfo, Site, HalfEdge
#### Constants ####
INT = np.int64
FLOAT = np.float64
cdef Vector2D NAN_VECTOR = _Vector2D(NAN, NAN)
cdef Matrix2x2 NAN_MATRIX = _Matrix2x2(NAN, NAN, NAN, NAN)
cdef FLOAT_T[18] SYMM = [0,0, 1,0, 1,1, 0,1, -1,1, -1,0, -1,-1, 0,-1, 1,-1]
cdef Matrix2x2 R = _Matrix2x2(0, -1, 1, 0)
cdef SiteCacheMap SITE_CACHE_MAP = _SiteCacheMap(0, 1, 2, 3, 4, 5, -1)
#### SiteCacheMap Methods ####
cdef inline SiteCacheMap _SiteCacheMap(INT_T iarea, INT_T iperim, INT_T iisoparam,
INT_T ienergy, INT_T iavg_radius,
INT_T icentroid, INT_T imaxcenter) nogil:
cdef SiteCacheMap sc
sc.iarea, sc.iperim, sc.iisoparam, sc.ienergy, sc.iavg_radius = (
iarea, iperim, iisoparam, ienergy, iavg_radius
)
sc.icentroid, sc.imaxcenter = icentroid, imaxcenter
sc.area, sc.perim, sc.isoparam, sc.energy, sc.avg_radius = (
area, perim, isoparam, energy, avg_radius
)
sc.centroid, sc.maxcenter = centroid, maxcenter
return sc
cdef inline FLOAT_T area(Site* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.iarea)
)
else:
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.iarea), val)
return val
cdef inline FLOAT_T perim(Site* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.iperim)
)
else:
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.iperim), val)
return val
cdef inline FLOAT_T isoparam(Site* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.iisoparam)
)
else:
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.iisoparam), val)
return val
cdef inline FLOAT_T energy(Site* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.ienergy)
)
else:
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.ienergy), val)
return val
cdef inline FLOAT_T avg_radius(Site* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.iavg_radius)
)
else:
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.iavg_radius), val)
return val
cdef inline Vector2D centroid(Site* self, Vector2D val) nogil:
if isnan(<double>val.x):
return _Vector2D(
self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.icentroid)
),
self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.icentroid+1)
)
)
else:
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.icentroid), val.x)
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.icentroid+1), val.y)
return val
cdef inline Vector2D maxcenter(Site* self, Vector2D val) nogil:
if isnan(<double>val.x):
return _Vector2D(
self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.imaxcenter)
),
self.info.site_cache.get(&self.info.site_cache,
(self.arr_index, self.cache.imaxcenter+1)
)
)
else:
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.imaxcenter), val.x)
self.info.site_cache.set(&self.info.site_cache,
(self.arr_index, self.cache.imaxcenter+1), val.y)
return val
#### EdgeCacheMap Methods ####
cdef inline EdgeCacheMap _EdgeCacheMap(INT_T iH, INT_T ila, INT_T ida, INT_T ixij,
INT_T idVdv, INT_T ila_mag, INT_T ida_mag,
INT_T iarea_p, INT_T icalI, INT_T size) nogil:
cdef EdgeCacheMap ec
ec.iH, ec.ila, ec.ida, ec.ixij, ec.idVdv = iH, ila, ida, ixij, idVdv
ec.ila_mag, ec.ida_mag, ec.iarea_p, ec.icalI = ila_mag, ida_mag, iarea_p, icalI
ec.size = size
ec.H, ec.la, ec.da, ec.xij, ec.dVdv = H, la, da, xij, dVdv
ec.la_mag, ec.da_mag, ec.area_p, ec.calI = la_mag, da_mag, area_p, calI
return ec
cdef inline Matrix2x2 H(HalfEdge* self, Matrix2x2 val) nogil:
if isnan(<double>val.a):
return _Matrix2x2(
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.iH)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.iH+1)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.iH+2)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.iH+3)
),
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.iH), val.a)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.iH+1), val.b)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.iH+2), val.c)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.iH+3), val.d)
return val
cdef inline Vector2D la(HalfEdge* self, Vector2D val) nogil:
if isnan(<double>val.x):
return _Vector2D(
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ila)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ila+1)
)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ila), val.x)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ila+1), val.y)
return val
cdef inline Vector2D da(HalfEdge* self, Vector2D val) nogil:
if isnan(<double>val.x):
return _Vector2D(
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ida)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ida+1)
)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ida), val.x)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ida+1), val.y)
return val
cdef inline Vector2D xij(HalfEdge* self, Vector2D val) nogil:
if isnan(<double>val.x):
return _Vector2D(
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ixij)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ixij+1)
)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ixij), val.x)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ixij+1), val.y)
return val
cdef inline Vector2D dVdv(HalfEdge* self, Vector2D val) nogil:
if isnan(<double>val.x):
return _Vector2D(
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.idVdv)
),
self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.idVdv+1)
)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.idVdv), val.x)
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.idVdv+1), val.y)
return val
cdef inline FLOAT_T la_mag(HalfEdge* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ila_mag)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ila_mag), val)
return val
cdef inline FLOAT_T da_mag(HalfEdge* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ida_mag)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ida_mag), val)
return val
cdef inline FLOAT_T area_p(HalfEdge* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.iarea_p)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.iarea_p), val)
return val
cdef inline FLOAT_T calI(HalfEdge* self, FLOAT_T val) nogil:
if isnan(<double>val):
return self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.icalI)
)
else:
self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.icalI), val)
return val
#### VoronoiInfo Methods ####
cdef inline VoronoiInfo _VoronoiInfo(INT_T [:, ::1] sites, INT_T [:, ::1] edges,
FLOAT_T [:, ::1] points, FLOAT_T [:, ::1] vertices,
FLOAT_T [:, ::1] site_cache,
FLOAT_T [:, ::1] edge_cache,
EdgeCacheMap* edge_cache_map) nogil:
cdef VoronoiInfo info
info.sites = _IArray(&sites[0,0], (<INT_T>sites.shape[0], <INT_T>sites.shape[1]))
info.edges = _IArray(&edges[0,0], (<INT_T>edges.shape[0], <INT_T>edges.shape[1]))
info.points = _FArray(
&points[0,0],
(<INT_T>points.shape[0], <INT_T>points.shape[1])
)
info.vertices = _FArray(
&vertices[0,0],
(<INT_T>vertices.shape[0], <INT_T>vertices.shape[1])
)
info.site_cache = _FArray(
&site_cache[0,0],
(<INT_T>site_cache.shape[0], <INT_T>site_cache.shape[1])
)
info.edge_cache = _FArray(
&edge_cache[0,0],
(<INT_T>edge_cache.shape[0], <INT_T>edge_cache.shape[1])
)
info.edge_cache_map = edge_cache_map
return info
#### Site Methods ####
cdef inline Site _Site(INT_T arr_index, VoronoiInfo* info) nogil:
cdef Site site
site.arr_index, site.info, site.cache = arr_index, info, &SITE_CACHE_MAP
site.index, site.vec, site.edge, site.edge_num = index, vec, edge, edge_num
return site
cdef inline INT_T index(Site* self) nogil:
return self.info.sites.get(&self.info.sites, (self.arr_index, 0))
cdef inline Vector2D vec(Site* self) nogil:
return _Vector2D(
self.info.points.get(&self.info.points, (self.index(self), 0)),
self.info.points.get(&self.info.points, (self.index(self), 1))
)
cdef inline HalfEdge edge(Site* self) nogil:
return _HalfEdge(
self.info.sites.get(&self.info.sites, (self.arr_index, 1)), self.info
)
cdef inline INT_T edge_num(Site* self) nogil:
return self.info.sites.get(&self.info.sites, (self.arr_index, 2))
#### HalfEdge Methods ####
cdef inline HalfEdge _HalfEdge(INT_T arr_index, VoronoiInfo* info) nogil:
cdef HalfEdge e
e.arr_index, e.info, e.cache = arr_index, info, info.edge_cache_map
e.orig_arr_index = arr_index
e.origin_index, e.origin, e.face, e.next, e.prev, e.twin, e.get_H = (
origin_index, origin, face, edge_next, prev, twin, get_H
)
return e
cdef inline INT_T origin_index(HalfEdge* self) nogil:
return self.info.edges.get(&self.info.edges, (self.arr_index, 0))
cdef inline Vector2D origin(HalfEdge* self) nogil:
return _Vector2D(
self.info.vertices.get(&self.info.vertices, (self.origin_index(self), 0)),
self.info.vertices.get(&self.info.vertices, (self.origin_index(self), 1))
)
cdef inline Site face(HalfEdge* self) nogil:
return _Site(
self.info.edges.get(&self.info.edges, (self.arr_index, 1)), self.info
)
cdef inline HalfEdge edge_next(HalfEdge* self) nogil:
return _HalfEdge(
self.info.edges.get(&self.info.edges, (self.arr_index, 2)), self.info
)
cdef inline HalfEdge prev(HalfEdge* self) nogil:
return _HalfEdge(
self.info.edges.get(&self.info.edges, (self.arr_index, 3)), self.info
)
cdef inline HalfEdge twin(HalfEdge* self) nogil:
return _HalfEdge(
self.info.edges.get(&self.info.edges, (self.arr_index, 4)), self.info
)
cdef inline Matrix2x2 get_H(HalfEdge* self, Site xi) nogil:
cdef INT_T this_e = self.origin_index(self)
cdef HalfEdge s_e = xi.edge(&xi)
for _ in range(xi.edge_num(&xi)):
if s_e.origin_index(&s_e) == this_e:
return s_e.cache.H(&s_e, NAN_MATRIX)
s_e = s_e.next(&s_e)
return _Matrix2x2(0.0, 0.0, 0.0, 0.0)
cdef class VoronoiContainer:
"""
Class for Voronoi diagrams, stored in a modified DCEL.
:param n: [int] how many sites to generate.
:param w: [float] width of the bounding domain.
:param h: [float] height of the bounding domain.
:param r: [float] radius of zero energy circle.
:param sites: np.ndarray collection of sites.
"""
def __init__(VoronoiContainer self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r,
object site_arr):
self.n, self.w, self.h, self.r = n, w, h, r
self.dim = [w, h]
self.calculate_voronoi(site_arr.astype(FLOAT))
self.generate_dcel()
self.common_cache()
self.precompute()
self.calc_grad()
self.get_statistics()
cdef void calculate_voronoi(VoronoiContainer self,
np.ndarray[FLOAT_T, ndim=2] site_arr) except *:
"""
Does all necessary computation and caching once points are set.
:param site_arr: initial points for this container.
"""
global SYMM
cdef np.ndarray[FLOAT_T, ndim=2] symm = np.asarray(SYMM).reshape(9,2)
cdef np.ndarray[FLOAT_T, ndim=1] dim = np.asarray(self.dim)
cdef np.ndarray[FLOAT_T, ndim=2] full_site_arr = np.empty(
(self.n*9+8, 2),
dtype=FLOAT
)
# Generate periodic sites and sites that bound periodic sites.
cdef INT_T i
for i in range(9):
full_site_arr[self.n*i:self.n*(i+1)] = site_arr + symm[i]*dim
if i > 0:
full_site_arr[9*self.n+i-1] = dim/2 + 2*dim*symm[i]
# Use SciPy to compute the Voronoi set.
self.scipy_vor = scipy.spatial.Voronoi(full_site_arr)
self.points = self.scipy_vor.points
self.vertices = self.scipy_vor.vertices
cdef void generate_dcel(VoronoiContainer self) except *:
cdef array.array int_tmplt = array.array('q', [])
cdef np.ndarray[INT_T, ndim=1] offsets = np.zeros(self.n*9+1, dtype=INT)
cdef array.array vert_indices = array.clone(int_tmplt, 0, False)
# Flatten regions into array, so it can be used later.
cdef INT_T i
for i in range(self.n*9):
verts = self.scipy_vor.regions[self.scipy_vor.point_region[i]]
offsets[i+1] = offsets[i] + len(verts) # Build offsets.
vert_indices.extend(array.array('q', verts)) # Flatten
# Get vertices of original N sites.
cdef np.ndarray[INT_T, ndim=1] vert_indices_np = np.asarray(vert_indices)
cdef np.ndarray[INT_T, ndim=1] border_sites = np.unique(np.searchsorted(
np.asarray(offsets), # Check indices where below matches would be inserted
np.nonzero(np.isin( # Indices of other verts being in bound verts.
vert_indices_np[offsets[self.n]:], # Rest of the verts to check.
np.unique(vert_indices_np[:offsets[self.n]]) # Bound verts
))[0] + offsets[self.n],
side='right' # If on index == offset_number, should be part of the next site.
) - 1) # Subtract by one to get actual site number.
cdef INT_T border_num = len(border_sites)
# Build sites array.
# [Site Index, Edge Index/Offset, Edge Count]
self.sites = np.empty((self.n+border_num, 3), dtype=INT)
self.sites.base[:self.n, 0] = np.arange(self.n, dtype=INT)
self.sites.base[self.n:, 0] = border_sites
self.sites.base[:self.n+1, 1] = offsets[:self.n+1]
for i in range(self.n):
self.sites[i, 2] = self.sites[i+1, 1] - self.sites[i, 1]
cdef INT_T edge_count = offsets[self.n]
cdef INT_T diff
for i in range(border_num):
diff = offsets[border_sites[i]+1] - offsets[border_sites[i]]
edge_count += diff
self.sites[self.n+i, 2] = diff
if i < border_num - 1:
self.sites[self.n+i+1, 1] = self.sites[self.n+i, 1] + diff
# Build edges array
# [Origin Index, Site Index, Next Index, Prev Index, Twin Index]
self.edges = np.empty((edge_count, 5), dtype=INT)
cdef np.ndarray[INT_T, ndim=1] site_verts
cdef INT_T j, site_i, edge_i, edge_offset, vert_num, twin_index
edge_indices = dict()
for i in range(self.n + border_num):
site_i = self.sites[i, 0]
edge_offset = self.sites[i, 1]
site_verts = vert_indices_np[offsets[site_i]:offsets[site_i+1]]
# Scipy outputs sorted vertices, but reverse if not counterclockwise.
if not VoronoiContainer.sign(self.points[site_i],
self.vertices[site_verts[0]], self.vertices[site_verts[1]]):
site_verts = np.flip(site_verts)
vert_num = offsets[site_i+1] - offsets[site_i]
for j in range(vert_num):
edge_i = edge_offset+j
self.edges[edge_i, 0] = site_verts[j]
self.edges[edge_i, 1] = i
# Add vert_num because of C modulo to get always positive.
self.edges[edge_i, 2] = (j+vert_num+1) % vert_num + edge_offset
self.edges[edge_i, 3] = (j+vert_num-1) % vert_num + edge_offset
# Get reversed tuple to theck for twin.
twin_index = edge_indices.get(
(site_verts[(j+1) % vert_num], site_verts[j]
), -1)
self.edges[edge_i, 4] = twin_index
if twin_index == -1:
edge_indices[(site_verts[j], site_verts[(j+1) % vert_num])] = \
j + edge_offset
else:
self.edges[twin_index, 4] = j + edge_offset
self.site_cache = np.empty((self.n + border_num, 7), dtype=FLOAT)
self.edge_cache = np.empty((edge_count, self.edge_cache_map.size), dtype=FLOAT)
cdef void common_cache(VoronoiContainer self) except *:
cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
cdef Site xi
cdef HalfEdge em, ep
cdef Vector2D p, q, la, da, Rla, centroid, cent_part
cdef FLOAT_T [:] area = np.zeros(self.sites.shape[0], dtype=FLOAT)
cdef FLOAT_T [:] perim = np.zeros(self.sites.shape[0], dtype=FLOAT)
cdef INT_T i, j
cdef FLOAT_T area_p, la_mag
for i in prange(self.sites.shape[0], nogil=True):
xi = _Site(i, &info)
centroid = _Vector2D(0, 0)
em = xi.edge(&xi)
j = 0
while j < xi.edge_num(&xi):
ep = em.next(&em)
p, q = em.origin(&em), ep.origin(&ep)
# vp - vm, vm - xi
la, da = q.copy.vsub(&q, p), p.copy.vsub(&p, xi.vec(&xi))
la_mag = la.mag(&la)
area_p = la.dot(&la, da.rot(&da))
Rla = la.rot(&la)
# Calculating centroid.
cent_part = p.copy.vadd(&p, q)
cent_part.self.vadd(&cent_part, xi.vec(&xi))
centroid.self.vadd(&centroid, cent_part.copy.smul(&cent_part, area_p))
# Caching
em.cache.la(&em, la)
em.cache.la_mag(&em, la_mag)
em.cache.da(&em, da)
em.cache.da_mag(&em, da.mag(&da))
em.cache.area_p(&em, area_p)
em.cache.xij(&em, Rla.copy.smul(&Rla, -area_p/la.dot(&la, la)))
area[i] += area_p
perim[i] += la_mag
em = em.next(&em)
j = j + 1
xi.cache.area(&xi, area[i]/2)
xi.cache.perim(&xi, perim[i])
xi.cache.isoparam(&xi, 2*PI*area[i]/(perim[i]**2))
xi.cache.centroid(&xi, centroid.copy.sdiv(&centroid, 3*area[i]))
@staticmethod
cdef inline Matrix2x2 calc_H(HalfEdge em, HalfEdge ep) nogil:
cdef Vector2D xmv, xpv, im, mp, right, Rpm, Rim, f
cdef Matrix2x2 h
cdef FLOAT_T im2, mp2
# Vectors from xi to xm and xp.
xmv, xpv = em.cache.xij(&em, NAN_VECTOR), ep.cache.xij(&ep, NAN_VECTOR)
im, mp = xmv.copy.neg(&xmv), xmv.copy.vsub(&xmv, xpv) # -xmv, xmv - xpv
im2, mp2 = -(xmv.dot(&xmv, xmv)), xmv.dot(&xmv, xmv) - xpv.dot(&xpv, xpv)
# (-xmv*xmv, xmv*xmv - xpv*xpv)
right = _Vector2D(im2, mp2)
Rpm, Rim = R.vecmul(&R, mp.copy.neg(&mp)), im.rot(&im) # R*-mp, R*im
h = _Matrix2x2(Rpm.x, Rim.x, Rpm.y, Rim.y) # [Rpm | Rim], h is temporary.
f = h.vecmul(&h, right) # [Rpm | Rim]*right
h = R.copy.smul(&R, mp2*(2*mp.dot(&mp, Rim))) # fp*g, g is a scalar.
# (fp*g - f*gp)/(g**2). f is a column vector, gp = 2*Rpm is a row vector.
h.self.msub(&h, _Matrix2x2(
f.x*2*Rpm.x, f.x*2*Rpm.y, f.y*2*Rpm.x, f.y*2*Rpm.y
))
h.self.sdiv(&h, (2*mp.dot(&mp, Rim))**2)
return h
@staticmethod
cdef inline bint sign(FLOAT_T [::1] ref, FLOAT_T [::1] p, FLOAT_T [::1] q):
"""
Outputs if p2 - self is counterclockwise of p1 - self.
:param p1: [List[float]] first vector
:param p2: [List[float]] second vector
:return: [bool] returns if counterclockwise.
"""
return ((q[0] - ref[0])*-(p[1] - ref[1]) + \
(q[1] - ref[1])*(p[0] - ref[0])) >= 0
# global ROT
# cdef np.ndarray[FLOAT_T, ndim=2] rot = np.asarray(ROT).reshape(2,2)
# return (q - ref).dot(rot.dot(p - ref)) >= 0
cdef void precompute(self) except *:
pass
cdef void calc_grad(self) except *:
pass
cdef void get_statistics(self) except *:
self.stats = {}
cache = self.site_cache[:self.n, :]
self.stats["site_areas"] = np.asarray(cache[:, SITE_CACHE_MAP.iarea])
self.stats["site_edge_count"] = np.asarray(self.sites[:self.n, 2])
self.stats["site_isos"] = np.asarray(cache[:, SITE_CACHE_MAP.iisoparam])
self.stats["site_energies"] = np.asarray(cache[:, SITE_CACHE_MAP.ienergy])
self.stats["avg_radius"] = np.asarray(cache[:, SITE_CACHE_MAP.iavg_radius])
self.stats["centroids"] = np.asarray(
cache[:, SITE_CACHE_MAP.icentroid:SITE_CACHE_MAP.icentroid+2]
)
self.stats["isoparam_avg"] = self.stats["site_areas"] / \
(PI*self.stats["avg_radius"]**2)
edges = np.asarray(self.edges)
mask = np.nonzero(edges[:, 0] != -1)[0]
all_edges = mask[(mask % 2 == 0)]
caches = edges[all_edges, 4]
edge_cache = np.asarray(self.edge_cache)
self.stats["edge_lengths"] = edge_cache[caches, self.edge_cache_map.ila_mag]
@property
def site_arr(self):
return np.asarray(self.points[:self.n], dtype=FLOAT)
@property
def vor_data(self):
return self.scipy_vor
@property
def gradient(self):
return np.asarray(self.grad, dtype=FLOAT)
def add_sites(self, add):
return (self.site_arr + add) % np.asarray(self.dim, dtype=FLOAT)
def iterate(self, FLOAT_T step):
k1 = self.gradient
k2 = self.__class__(self.n, self.w, self.h, self.r,
self.add_sites(step*k1)
).gradient
return (step/2)*(k1+k2), k1
def hessian(self, d: float) -> np.ndarray:
"""
Obtains the approximate Hessian.
:param d: [float] small d for approximation.
:return: 2Nx2N array that represents Hessian.
"""
HE = np.zeros((2*self.n, 2*self.n))
new_sites = np.copy(self.site_arr) # Maintain one copy for speed.
for i in range(self.n):
for j in range(2):
mod = self.w if j == 0 else self.h
new_sites[i][j] = (new_sites[i][j] + d) % mod
Ep = self.__class__(self.n, self.w, self.h, self.r, new_sites)
new_sites[i][j] = (new_sites[i][j] - 2*d) % mod
Em = self.__class__(self.n, self.w, self.h, self.r, new_sites)
new_sites[i][j] = (new_sites[i][j] + d) % mod
HE[:, 2*i+j] = ((Ep.gradient - Em.gradient)/(2*d)).flatten()
# Average out discrepencies, since it should be symmetric.
for i in range(2*self.n):
for j in range(i, 2*self.n):
HE[i][j] = (HE[i][j] + HE[j][i])/2
HE[j][i] = HE[i][j]
return HE
def site_vert_arr(self): # -> List[np.ndarray]
cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
cdef INT_T i, j
cdef Site xi
cdef HalfEdge e
cdef Vector2D v
sites, site_verts = [], []
for i in range(self.n):
xi = _Site(i, &info)
v = xi.vec(&xi)
sites.append(np.array([v.x, v.y]))
verts = np.empty((xi.edge_num(&xi), 2))
e = xi.edge(&xi)
for j in range(xi.edge_num(&xi)):
v = e.origin(&e)
verts[j, 0], verts[j, 1] = v.x, v.y
e = e.next(&e)
site_verts.append(verts)
return sites, site_verts