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

Modified scripts to use aspect ratios, Radial[T] formula update for gradient

Browse Source
This commit is contained in:
Kenneth Jao 2021-12-19 03:35:11 -05:00
parent 2e9ed2f072
commit 8cd48a1acd
7 changed files with 2069 additions and 1179 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
scripts/aspect_diagrams.py
View File

@ -25,14 +25,14 @@ def order_process(domain: DomainParams) -> Tuple[float, float, float]:
isoparams.append(math.pi * rbar ** 2 / area) isoparams.append(math.pi * rbar ** 2 / area)
return domain.w, min(energies), max(energies), min(isoparams), max(isoparams) return (domain.w, min(energies), max(energies), min(isoparams), max(isoparams))
def get_ordered_energies(orig_domain: DomainParams, widths: np.ndarray) -> Dict: def get_ordered_energies(orig_domain: DomainParams, widths: np.ndarray) -> Dict:
data = {} data = {}
domains = [] domains = []
for w in widths: for w in widths:
aspect = w / orig_domain.h aspect = w
domains.append( domains.append(
DomainParams( DomainParams(
orig_domain.n, orig_domain.n,
@ -102,7 +102,7 @@ def eq_file_process(file: Path) -> Tuple[float, List[float], List[float]]:
] ]
) )
return sim.domain.w, alls, distincts return sim.domain.w / sim.domain.h, alls, distincts
def get_equilibria_data(filepath: Path) -> Tuple[Dict, numpy.ndarray, DomainParams]: def get_equilibria_data(filepath: Path) -> Tuple[Dict, numpy.ndarray, DomainParams]:
@ -132,7 +132,7 @@ def get_equilibria_data(filepath: Path) -> Tuple[Dict, numpy.ndarray, DomainPara
def axis_settings(ax, widths): def axis_settings(ax, widths):
ax.grid(zorder=0) ax.grid(zorder=0)
ax.set_xticks([round(w, 2) for w in widths[::2]]) ax.set_xticks([round(w, 2) for w in widths[::2]])
ax.set_xticklabels([f"{round(w / 10, 3):.2f}" for w in widths[::2]], rotation=90) ax.set_xticklabels([f"{round(w, 3):.2f}" for w in widths[::2]], rotation=90)
plt.subplots_adjust(0.07, 0.12, 0.97, 0.9) plt.subplots_adjust(0.07, 0.12, 0.97, 0.9)

166
scripts/heatmap.py
View File

@ -1,14 +1,78 @@
from __future__ import annotations from __future__ import annotations
from typing import List, Tuple, Dict from typing import List, Tuple, Dict
import argparse, numpy as np, os, pickle import argparse, math, numpy as np, os, pickle
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib.ticker as mtick
from scipy.optimize import curve_fit
from multiprocessing import Pool, cpu_count from multiprocessing import Pool, cpu_count
from pathlib import Path from pathlib import Path
import squish.ordered as order
from squish import Simulation, DomainParams from squish import Simulation, DomainParams
from squish.common import OUTPUT_DIR from squish.common import OUTPUT_DIR
def order_process(domain: DomainParams) -> Tuple[float, float, float]:
energies, isoparams = [], []
configs = order.configurations(domain)
for config in configs:
rbar = order.avg_radius(domain, config)
area = domain.w * domain.h / domain.n
energies.append(
2 * domain.w * domain.h
+ 2 * math.pi * domain.n * (domain.r ** 2 - 2 * domain.r * rbar)
)
isoparams.append(math.pi * rbar ** 2 / area)
return domain.w, min(energies), max(energies), min(isoparams), max(isoparams)
def get_ordered_energies(orig_domain: DomainParams, widths: np.ndarray) -> Dict:
data = {}
domains = []
for w in widths:
aspect = w
domains.append(
DomainParams(
orig_domain.n,
math.sqrt(orig_domain.n * aspect),
math.sqrt(orig_domain.n / aspect),
orig_domain.r,
)
)
# domains = [
# DomainParams(orig_domain.n, w, orig_domain.h, orig_domain.r) for w in widths
# ]
with Pool(cpu_count()) as pool:
energy_mins, energy_maxes, isoparam_mins, isoparam_maxes = {}, {}, {}, {}
for i, res in enumerate(pool.imap_unordered(order_process, domains)):
energy_mins[res[0]] = res[1]
energy_maxes[res[0]] = res[2]
isoparam_mins[res[0]] = res[3]
isoparam_maxes[res[0]] = res[4]
hashes = int(21 * i / len(widths))
print(
f'Generating at width {res[0]:.02f}... |{"#"*hashes}{" "*(20-hashes)}|'
+ f" {i+1}/{len(widths)} completed.",
flush=True,
end="\r",
)
print(flush=True)
data["energy_min"] = list([x[1] for x in sorted(energy_mins.items())])
data["energy_max"] = list([x[1] for x in sorted(energy_maxes.items())])
data["isoparam_min"] = list([x[1] for x in sorted(isoparam_mins.items())])
data["isoparam_max"] = list([x[1] for x in sorted(isoparam_maxes.items())])
return data
def eq_file_process(file: Path) -> Tuple[float, List[float], List[float]]: def eq_file_process(file: Path) -> Tuple[float, List[float], List[float]]:
sim, frames = Simulation.load(file) sim, frames = Simulation.load(file)
@ -19,6 +83,7 @@ def eq_file_process(file: Path) -> Tuple[float, List[float], List[float]]:
frame_info["energy"], frame_info["energy"],
np.var(frame_info["stats"]["avg_radius"]) <= 1e-8, np.var(frame_info["stats"]["avg_radius"]) <= 1e-8,
np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6), np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6),
sum(frame_info["stats"]["site_energies"][: sim.domain.n]),
] ]
) )
@ -33,17 +98,17 @@ def eq_file_process(file: Path) -> Tuple[float, List[float], List[float]]:
frame_info["energy"], frame_info["energy"],
np.var(frame_info["stats"]["avg_radius"]) <= 1e-8, np.var(frame_info["stats"]["avg_radius"]) <= 1e-8,
np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6), np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6),
sum(frame_info["stats"]["site_energies"][: sim.domain.n]),
counts[j], counts[j],
] ]
) )
return sim.domain.w, alls, distincts return sim.domain.w / sim.domain.h, alls, distincts
def get_equilibria_data(filepath: Path) -> Tuple[Dict, numpy.ndarray, DomainParams]: def get_equilibria_data(filepath: Path) -> Tuple[Dict, numpy.ndarray, DomainParams]:
data = {"all": {}, "distinct": {}} data = {"all": {}, "distinct": {}}
files = list(Path(filepath).iterdir()) files = list(Path(filepath).iterdir())
sim, frames = Simulation.load(files[0])
with Pool(cpu_count()) as pool: with Pool(cpu_count()) as pool:
for i, res in enumerate(pool.imap_unordered(eq_file_process, files)): for i, res in enumerate(pool.imap_unordered(eq_file_process, files)):
@ -52,18 +117,55 @@ def get_equilibria_data(filepath: Path) -> Tuple[Dict, numpy.ndarray, DomainPara
hashes = int(21 * i / len(files)) hashes = int(21 * i / len(files))
print( print(
f'Loading simulations for N={sim.domain.n}... |{"#"*hashes}{" "*(20-hashes)}|' f'Loading simulations... |{"#"*hashes}{" "*(20-hashes)}|'
+ f" {i+1}/{len(files)} simulations loaded.", + f" {i+1}/{len(files)} simulations loaded.",
flush=True, flush=True,
end="\r", end="\r",
) )
print(flush=True) print(flush=True)
sim, frames = Simulation.load(files[0])
widths = np.asarray(sorted(data["all"])) widths = np.asarray(sorted(data["all"]))
domain = DomainParams(sim.domain.n, widths[-1], sim.domain.h, sim.domain.r) domain = DomainParams(sim.domain.n, widths[-1], sim.domain.h, sim.domain.r)
return data, widths, domain return data, widths, domain
def probability_of_disorder(data, widths, domain):
fig, ax = plt.subplots(figsize=(16, 8))
all_disorder_count = []
for width in widths:
equal_shape = list([c[1] for c in data["all"][width]])
all_disorder_count.append(
100 * equal_shape.count(False) / len(data["all"][width])
)
return all_disorder_count
def excess_energy(data, widths, order_data, domain):
fig, ax = plt.subplots(figsize=(16, 8))
ordered_energies, unordered_energies = [], []
for width in widths:
ordered_energies.append([c[0] for c in data["distinct"][width] if c[1]])
unordered_energies.append([c[0] for c in data["distinct"][width] if not c[1]])
for i in range(len(order_data["energy_min"])):
ordered_energies[i].append(order_data["energy_min"][i])
ordered_energies[i].append(order_data["energy_max"][i])
min_order = np.asarray([min(width) for width in ordered_energies])
max_order = np.asarray([max(width) for width in ordered_energies])
min_unorder = np.asarray([min(width) for width in unordered_energies])
max_unorder = np.asarray([max(width) for width in unordered_energies])
return min_order - min_unorder
def sigmoid(x, x0, k):
return 100 / (1 + np.exp(-k * (x - x0)))
def main(): def main():
# Loading arguments. # Loading arguments.
parser = argparse.ArgumentParser("Outputs width search data into diagrams") parser = argparse.ArgumentParser("Outputs width search data into diagrams")
@ -83,6 +185,62 @@ def main():
args = parser.parse_args() args = parser.parse_args()
fig_folder = OUTPUT_DIR
fig_folder.mkdir(exist_ok=True)
store = Path(args.sims_path) / "EEvsPoD.pkl"
if store.is_file():
with open(store, "rb") as f:
horiz, vert = pickle.load(f)
else:
horiz = []
vert = []
for file in Path(args.sims_path).iterdir():
# Obtain data from simulation files and generate single shape data.
data, widths, domain = get_equilibria_data(file)
order_data = get_ordered_energies(domain, widths)
vert.append(probability_of_disorder(data, widths, domain))
horiz.append(excess_energy(data, widths, order_data, domain))
horiz, vert = np.concatenate(horiz), np.concatenate(vert)
with open(store, "wb") as f:
pickle.dump((horiz, vert), f, pickle.HIGHEST_PROTOCOL)
fig, ax = plt.subplots(figsize=(10, 10))
for i in range(2):
ax.scatter(
horiz[i * 141 : (i + 1) * 141],
vert[i * 141 : (i + 1) * 141],
alpha=0.5,
color=f"C{i}",
s=5,
)
start, end = ax.get_xlim()
# popt, pcov = curve_fit(sigmoid, horiz, vert)
# x = np.linspace(start, end, 100)
# y = sigmoid(x, *popt)
# y = sigmoid(x, -1.35, 3)
# ax.plot(x, y, color="C1")
plt.subplots_adjust(0.1, 0.1, 0.97, 0.93)
ax.set_xticks(np.linspace(start, end, 10))
ax.set_yticks(np.arange(0, 105, 5))
ax.grid()
ax.yaxis.set_major_formatter(mtick.PercentFormatter())
ax.title.set_text("Excess Energy Difference vs. PoD")
ax.set_xlabel("Excess Energy Difference")
ax.set_ylabel("Probability of Disorder")
fig.savefig(OUTPUT_DIR / "Energy Diff and Probability")
return
# with open("testing.pkl", "rb") as f: # with open("testing.pkl", "rb") as f:
# disorder_dict = pickle.load(f) # disorder_dict = pickle.load(f)
# widths = np.linspace(3.0, 10.0, 141) # widths = np.linspace(3.0, 10.0, 141)

562
squish/energy.c
View File

@ -1574,28 +1574,28 @@ struct __pyx_t_6squish_7voronoi_SiteCacheMap {
* *
* # Psuedo-class that handles caching for edges. * # Psuedo-class that handles caching for edges.
* ctypedef struct EdgeCacheMap: # <<<<<<<<<<<<<< * ctypedef struct EdgeCacheMap: # <<<<<<<<<<<<<<
* INT_T iH, ila, ida, ixij, idVdv, ila_mag, ida_mag, iarea_p, icalI, size * INT_T iH, ila, ida, iya, idVdv, ila_mag, ida_mag, iya_mag, icalI, size
* *
*/ */
struct __pyx_t_6squish_7voronoi_EdgeCacheMap { struct __pyx_t_6squish_7voronoi_EdgeCacheMap {
__pyx_t_6squish_4core_INT_T iH; __pyx_t_6squish_4core_INT_T iH;
__pyx_t_6squish_4core_INT_T ila; __pyx_t_6squish_4core_INT_T ila;
__pyx_t_6squish_4core_INT_T ida; __pyx_t_6squish_4core_INT_T ida;
__pyx_t_6squish_4core_INT_T ixij; __pyx_t_6squish_4core_INT_T iya;
__pyx_t_6squish_4core_INT_T idVdv; __pyx_t_6squish_4core_INT_T idVdv;
__pyx_t_6squish_4core_INT_T ila_mag; __pyx_t_6squish_4core_INT_T ila_mag;
__pyx_t_6squish_4core_INT_T ida_mag; __pyx_t_6squish_4core_INT_T ida_mag;
__pyx_t_6squish_4core_INT_T iarea_p; __pyx_t_6squish_4core_INT_T iya_mag;
__pyx_t_6squish_4core_INT_T icalI; __pyx_t_6squish_4core_INT_T icalI;
__pyx_t_6squish_4core_INT_T size; __pyx_t_6squish_4core_INT_T size;
__pyx_t_6squish_4core_Matrix2x2 (*H)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_Matrix2x2); __pyx_t_6squish_4core_Matrix2x2 (*H)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_Matrix2x2);
__pyx_t_6squish_4core_Vector2D (*la)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_Vector2D); __pyx_t_6squish_4core_Vector2D (*la)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_Vector2D);
__pyx_t_6squish_4core_Vector2D (*da)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_Vector2D); __pyx_t_6squish_4core_Vector2D (*da)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_Vector2D);
__pyx_t_6squish_4core_Vector2D (*xij)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_Vector2D); __pyx_t_6squish_4core_Vector2D (*ya)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_Vector2D);
__pyx_t_6squish_4core_Vector2D (*dVdv)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_Vector2D); __pyx_t_6squish_4core_Vector2D (*dVdv)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_Vector2D);
__pyx_t_6squish_4core_FLOAT_T (*la_mag)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_FLOAT_T); __pyx_t_6squish_4core_FLOAT_T (*la_mag)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_FLOAT_T);
__pyx_t_6squish_4core_FLOAT_T (*da_mag)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_FLOAT_T); __pyx_t_6squish_4core_FLOAT_T (*da_mag)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_FLOAT_T);
__pyx_t_6squish_4core_FLOAT_T (*area_p)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_FLOAT_T); __pyx_t_6squish_4core_FLOAT_T (*ya_mag)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_FLOAT_T);
__pyx_t_6squish_4core_FLOAT_T (*calI)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_FLOAT_T); __pyx_t_6squish_4core_FLOAT_T (*calI)(__pyx_t_6squish_7voronoi_HalfEdge *, __pyx_t_6squish_4core_FLOAT_T);
}; };
@ -4142,7 +4142,7 @@ static void __pyx_f_6squish_6energy_10AreaEnergy_calc_grad(struct __pyx_obj_6squ
* dedxi_p = f.cache.dVdv(&f, NAN_VECTOR) #dVdv * dedxi_p = f.cache.dVdv(&f, NAN_VECTOR) #dVdv
* dedxi_p.self.smul(&dedxi_p, xf.cache.area(&xf, NAN) - A) # <<<<<<<<<<<<<< * dedxi_p.self.smul(&dedxi_p, xf.cache.area(&xf, NAN) - A) # <<<<<<<<<<<<<<
* dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX)) * dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX))
* dedx[i][0] -= dedxi_p.x * dedx[i][0] += dedxi_p.x
*/ */
(void)(__pyx_v_dedxi_p.self.smul((&__pyx_v_dedxi_p), (__pyx_v_xf.cache->area((&__pyx_v_xf), NAN) - __pyx_v_A))); (void)(__pyx_v_dedxi_p.self.smul((&__pyx_v_dedxi_p), (__pyx_v_xf.cache->area((&__pyx_v_xf), NAN) - __pyx_v_A)));
@ -4150,32 +4150,32 @@ static void __pyx_f_6squish_6energy_10AreaEnergy_calc_grad(struct __pyx_obj_6squ
* dedxi_p = f.cache.dVdv(&f, NAN_VECTOR) #dVdv * dedxi_p = f.cache.dVdv(&f, NAN_VECTOR) #dVdv
* dedxi_p.self.smul(&dedxi_p, xf.cache.area(&xf, NAN) - A) * dedxi_p.self.smul(&dedxi_p, xf.cache.area(&xf, NAN) - A)
* dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX)) # <<<<<<<<<<<<<< * dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX)) # <<<<<<<<<<<<<<
* dedx[i][0] -= dedxi_p.x * dedx[i][0] += dedxi_p.x
* dedx[i][1] -= dedxi_p.y * dedx[i][1] += dedxi_p.y
*/ */
(void)(__pyx_v_dedxi_p.self.matmul((&__pyx_v_dedxi_p), __pyx_v_e.cache->H((&__pyx_v_e), __pyx_v_6squish_7voronoi_NAN_MATRIX))); (void)(__pyx_v_dedxi_p.self.matmul((&__pyx_v_dedxi_p), __pyx_v_e.cache->H((&__pyx_v_e), __pyx_v_6squish_7voronoi_NAN_MATRIX)));
/* "squish/energy.pyx":102 /* "squish/energy.pyx":102
* dedxi_p.self.smul(&dedxi_p, xf.cache.area(&xf, NAN) - A) * dedxi_p.self.smul(&dedxi_p, xf.cache.area(&xf, NAN) - A)
* dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX)) * dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX))
* dedx[i][0] -= dedxi_p.x # <<<<<<<<<<<<<< * dedx[i][0] += dedxi_p.x # <<<<<<<<<<<<<<
* dedx[i][1] -= dedxi_p.y * dedx[i][1] += dedxi_p.y
* *
*/ */
__pyx_t_10 = __pyx_v_i; __pyx_t_10 = __pyx_v_i;
__pyx_t_11 = 0; __pyx_t_11 = 0;
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ ((char *) (((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_dedx.data + __pyx_t_10 * __pyx_v_dedx.strides[0]) )) + __pyx_t_11)) )) -= __pyx_v_dedxi_p.x; *((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ ((char *) (((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_dedx.data + __pyx_t_10 * __pyx_v_dedx.strides[0]) )) + __pyx_t_11)) )) += __pyx_v_dedxi_p.x;
/* "squish/energy.pyx":103 /* "squish/energy.pyx":103
* dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX)) * dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX))
* dedx[i][0] -= dedxi_p.x * dedx[i][0] += dedxi_p.x
* dedx[i][1] -= dedxi_p.y # <<<<<<<<<<<<<< * dedx[i][1] += dedxi_p.y # <<<<<<<<<<<<<<
* *
* f = f.twin(&f) * f = f.twin(&f)
*/ */
__pyx_t_10 = __pyx_v_i; __pyx_t_10 = __pyx_v_i;
__pyx_t_11 = 1; __pyx_t_11 = 1;
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ ((char *) (((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_dedx.data + __pyx_t_10 * __pyx_v_dedx.strides[0]) )) + __pyx_t_11)) )) -= __pyx_v_dedxi_p.y; *((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ ((char *) (((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_dedx.data + __pyx_t_10 * __pyx_v_dedx.strides[0]) )) + __pyx_t_11)) )) += __pyx_v_dedxi_p.y;
/* "squish/energy.pyx":97 /* "squish/energy.pyx":97
* f = e * f = e
@ -4187,7 +4187,7 @@ static void __pyx_f_6squish_6energy_10AreaEnergy_calc_grad(struct __pyx_obj_6squ
} }
/* "squish/energy.pyx":105 /* "squish/energy.pyx":105
* dedx[i][1] -= dedxi_p.y * dedx[i][1] += dedxi_p.y
* *
* f = f.twin(&f) # <<<<<<<<<<<<<< * f = f.twin(&f) # <<<<<<<<<<<<<<
* f = f.next(&f) * f = f.next(&f)
@ -5251,7 +5251,6 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_precompute(struct __pyx_obj_
__pyx_t_6squish_7voronoi_VoronoiInfo __pyx_v_info; __pyx_t_6squish_7voronoi_VoronoiInfo __pyx_v_info;
__pyx_t_6squish_7voronoi_Site __pyx_v_xi; __pyx_t_6squish_7voronoi_Site __pyx_v_xi;
__pyx_t_6squish_7voronoi_HalfEdge __pyx_v_e; __pyx_t_6squish_7voronoi_HalfEdge __pyx_v_e;
__pyx_t_6squish_7voronoi_HalfEdge __pyx_v_em;
__pyx_t_6squish_7voronoi_HalfEdge __pyx_v_ep; __pyx_t_6squish_7voronoi_HalfEdge __pyx_v_ep;
__pyx_t_6squish_4core_Vector2D __pyx_v_la; __pyx_t_6squish_4core_Vector2D __pyx_v_la;
__Pyx_memviewslice __pyx_v_site_energy = { 0, 0, { 0 }, { 0 }, { 0 } }; __Pyx_memviewslice __pyx_v_site_energy = { 0, 0, { 0 }, { 0 }, { 0 } };
@ -5273,11 +5272,9 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_precompute(struct __pyx_obj_
__pyx_t_6squish_4core_INT_T __pyx_t_10; __pyx_t_6squish_4core_INT_T __pyx_t_10;
__pyx_t_6squish_4core_INT_T __pyx_t_11; __pyx_t_6squish_4core_INT_T __pyx_t_11;
int __pyx_t_12; int __pyx_t_12;
__pyx_t_6squish_7voronoi_HalfEdge __pyx_t_13; __pyx_t_6squish_4core_INT_T __pyx_t_13;
__pyx_t_6squish_7voronoi_HalfEdge __pyx_t_14; __pyx_t_6squish_4core_INT_T __pyx_t_14;
__pyx_t_6squish_4core_INT_T __pyx_t_15; __pyx_t_6squish_4core_FLOAT_T __pyx_t_15;
__pyx_t_6squish_4core_INT_T __pyx_t_16;
__pyx_t_6squish_4core_FLOAT_T __pyx_t_17;
int __pyx_lineno = 0; int __pyx_lineno = 0;
const char *__pyx_filename = NULL; const char *__pyx_filename = NULL;
int __pyx_clineno = 0; int __pyx_clineno = 0;
@ -5451,11 +5448,11 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_precompute(struct __pyx_obj_
if (__pyx_t_11 > 0) if (__pyx_t_11 > 0)
{ {
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp parallel private(__pyx_t_12, __pyx_t_13, __pyx_t_14, __pyx_t_15, __pyx_t_16) #pragma omp parallel private(__pyx_t_12, __pyx_t_13, __pyx_t_14)
#endif /* _OPENMP */ #endif /* _OPENMP */
{ {
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp for lastprivate(__pyx_v_e) lastprivate(__pyx_v_em) lastprivate(__pyx_v_ep) firstprivate(__pyx_v_i) lastprivate(__pyx_v_i) lastprivate(__pyx_v_j) lastprivate(__pyx_v_la) lastprivate(__pyx_v_sm) lastprivate(__pyx_v_sp) lastprivate(__pyx_v_xi) #pragma omp for lastprivate(__pyx_v_e) lastprivate(__pyx_v_ep) firstprivate(__pyx_v_i) lastprivate(__pyx_v_i) lastprivate(__pyx_v_j) lastprivate(__pyx_v_la) lastprivate(__pyx_v_sm) lastprivate(__pyx_v_sp) lastprivate(__pyx_v_xi)
#endif /* _OPENMP */ #endif /* _OPENMP */
for (__pyx_t_10 = 0; __pyx_t_10 < __pyx_t_11; __pyx_t_10++){ for (__pyx_t_10 = 0; __pyx_t_10 < __pyx_t_11; __pyx_t_10++){
{ {
@ -5488,7 +5485,7 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_precompute(struct __pyx_obj_
* e = xi.edge(&xi) * e = xi.edge(&xi)
* j = 0 # <<<<<<<<<<<<<< * j = 0 # <<<<<<<<<<<<<<
* while j < xi.edge_num(&xi): * while j < xi.edge_num(&xi):
* em, ep = e.prev(&e), e.next(&e) * ep = e.next(&e)
*/ */
__pyx_v_j = 0; __pyx_v_j = 0;
@ -5496,8 +5493,8 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_precompute(struct __pyx_obj_
* e = xi.edge(&xi) * e = xi.edge(&xi)
* j = 0 * j = 0
* while j < xi.edge_num(&xi): # <<<<<<<<<<<<<< * while j < xi.edge_num(&xi): # <<<<<<<<<<<<<<
* em, ep = e.prev(&e), e.next(&e) * ep = e.next(&e)
* e.cache.H(&e, VoronoiContainer.calc_H(em, e)) * #e.cache.H(&e, VoronoiContainer.calc_H(em, e))
*/ */
while (1) { while (1) {
__pyx_t_12 = ((__pyx_v_j < __pyx_v_xi.edge_num((&__pyx_v_xi))) != 0); __pyx_t_12 = ((__pyx_v_j < __pyx_v_xi.edge_num((&__pyx_v_xi))) != 0);
@ -5506,26 +5503,14 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_precompute(struct __pyx_obj_
/* "squish/energy.pyx":191 /* "squish/energy.pyx":191
* j = 0 * j = 0
* while j < xi.edge_num(&xi): * while j < xi.edge_num(&xi):
* em, ep = e.prev(&e), e.next(&e) # <<<<<<<<<<<<<< * ep = e.next(&e) # <<<<<<<<<<<<<<
* e.cache.H(&e, VoronoiContainer.calc_H(em, e)) * #e.cache.H(&e, VoronoiContainer.calc_H(em, e))
* *
*/ */
__pyx_t_13 = __pyx_v_e.prev((&__pyx_v_e)); __pyx_v_ep = __pyx_v_e.next((&__pyx_v_e));
__pyx_t_14 = __pyx_v_e.next((&__pyx_v_e));
__pyx_v_em = __pyx_t_13;
__pyx_v_ep = __pyx_t_14;
/* "squish/energy.pyx":192
* 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)
*/
(void)(__pyx_v_e.cache->H((&__pyx_v_e), __pyx_vtabptr_6squish_7voronoi_VoronoiContainer->calc_H(__pyx_v_em, __pyx_v_e)));
/* "squish/energy.pyx":194 /* "squish/energy.pyx":194
* e.cache.H(&e, VoronoiContainer.calc_H(em, e)) * #e.cache.H(&e, VoronoiContainer.calc_H(em, e))
* *
* la = e.cache.la(&e, NAN_VECTOR) # <<<<<<<<<<<<<< * la = e.cache.la(&e, NAN_VECTOR) # <<<<<<<<<<<<<<
* sp = la.dot(&la, ep.cache.da(&ep, NAN_VECTOR)) # dap . la * sp = la.dot(&la, ep.cache.da(&ep, NAN_VECTOR)) # dap . la
@ -5574,22 +5559,22 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_precompute(struct __pyx_obj_
* *
* e.cache.calI(&e, <FLOAT_T>(atanh(<double>sp) - atanh(<double>sm))) # <<<<<<<<<<<<<< * e.cache.calI(&e, <FLOAT_T>(atanh(<double>sp) - atanh(<double>sm))) # <<<<<<<<<<<<<<
* *
* avg_radii[i] += ( * avg_radii[i] += e.cache.ya_mag(&e, NAN) * e.cache.calI(&e, NAN) / 2
*/ */
(void)(__pyx_v_e.cache->calI((&__pyx_v_e), ((__pyx_t_6squish_4core_FLOAT_T)(atanh(((double)__pyx_v_sp)) - atanh(((double)__pyx_v_sm)))))); (void)(__pyx_v_e.cache->calI((&__pyx_v_e), ((__pyx_t_6squish_4core_FLOAT_T)(atanh(((double)__pyx_v_sp)) - atanh(((double)__pyx_v_sm))))));
/* "squish/energy.pyx":203 /* "squish/energy.pyx":203
* e.cache.calI(&e, <FLOAT_T>(atanh(<double>sp) - atanh(<double>sm))) * e.cache.calI(&e, <FLOAT_T>(atanh(<double>sp) - atanh(<double>sm)))
* *
* avg_radii[i] += ( # <<<<<<<<<<<<<< * avg_radii[i] += e.cache.ya_mag(&e, NAN) * e.cache.calI(&e, NAN) / 2 # <<<<<<<<<<<<<<
* e.cache.area_p(&e, NAN) * e.cache.calI(&e, NAN) *
* / e.cache.la_mag(&e, NAN) * e = e.next(&e)
*/ */
__pyx_t_15 = __pyx_v_i; __pyx_t_13 = __pyx_v_i;
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_avg_radii.data + __pyx_t_15 * __pyx_v_avg_radii.strides[0]) )) += ((__pyx_v_e.cache->area_p((&__pyx_v_e), NAN) * __pyx_v_e.cache->calI((&__pyx_v_e), NAN)) / __pyx_v_e.cache->la_mag((&__pyx_v_e), NAN)); *((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_avg_radii.data + __pyx_t_13 * __pyx_v_avg_radii.strides[0]) )) += ((__pyx_v_e.cache->ya_mag((&__pyx_v_e), NAN) * __pyx_v_e.cache->calI((&__pyx_v_e), NAN)) / 2.0);
/* "squish/energy.pyx":208 /* "squish/energy.pyx":205
* ) * avg_radii[i] += e.cache.ya_mag(&e, NAN) * e.cache.calI(&e, NAN) / 2
* *
* e = e.next(&e) # <<<<<<<<<<<<<< * e = e.next(&e) # <<<<<<<<<<<<<<
* j = j + 1 * j = j + 1
@ -5597,7 +5582,7 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_precompute(struct __pyx_obj_
*/ */
__pyx_v_e = __pyx_v_e.next((&__pyx_v_e)); __pyx_v_e = __pyx_v_e.next((&__pyx_v_e));
/* "squish/energy.pyx":209 /* "squish/energy.pyx":206
* *
* e = e.next(&e) * e = e.next(&e)
* j = j + 1 # <<<<<<<<<<<<<< * j = j + 1 # <<<<<<<<<<<<<<
@ -5607,36 +5592,36 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_precompute(struct __pyx_obj_
__pyx_v_j = (__pyx_v_j + 1); __pyx_v_j = (__pyx_v_j + 1);
} }
/* "squish/energy.pyx":211 /* "squish/energy.pyx":208
* j = j + 1 * j = j + 1
* *
* site_energy[i] += 2*(xi.cache.area(&xi, NAN) - self.r*avg_radii[i]) # <<<<<<<<<<<<<< * 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.avg_radius(&xi, avg_radii[i]/(2*PI))
*/ */
__pyx_t_15 = __pyx_v_i; __pyx_t_13 = __pyx_v_i;
__pyx_t_16 = __pyx_v_i; __pyx_t_14 = __pyx_v_i;
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_site_energy.data + __pyx_t_16 * __pyx_v_site_energy.strides[0]) )) += (2.0 * (__pyx_v_xi.cache->area((&__pyx_v_xi), NAN) - (__pyx_v_self->__pyx_base.r * (*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_avg_radii.data + __pyx_t_15 * __pyx_v_avg_radii.strides[0]) )))))); *((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_site_energy.data + __pyx_t_14 * __pyx_v_site_energy.strides[0]) )) += (2.0 * (__pyx_v_xi.cache->area((&__pyx_v_xi), NAN) - (__pyx_v_self->__pyx_base.r * (*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_avg_radii.data + __pyx_t_13 * __pyx_v_avg_radii.strides[0]) ))))));
/* "squish/energy.pyx":213 /* "squish/energy.pyx":210
* site_energy[i] += 2*(xi.cache.area(&xi, NAN) - self.r*avg_radii[i]) * 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.avg_radius(&xi, avg_radii[i]/(2*PI)) # <<<<<<<<<<<<<<
* xi.cache.energy(&xi, site_energy[i]) * xi.cache.energy(&xi, site_energy[i])
* *
*/ */
__pyx_t_15 = __pyx_v_i; __pyx_t_13 = __pyx_v_i;
(void)(__pyx_v_xi.cache->avg_radius((&__pyx_v_xi), ((*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_avg_radii.data + __pyx_t_15 * __pyx_v_avg_radii.strides[0]) ))) / ((__pyx_t_6squish_4core_FLOAT_T)(2.0 * M_PI))))); (void)(__pyx_v_xi.cache->avg_radius((&__pyx_v_xi), ((*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_avg_radii.data + __pyx_t_13 * __pyx_v_avg_radii.strides[0]) ))) / ((__pyx_t_6squish_4core_FLOAT_T)(2.0 * M_PI)))));
/* "squish/energy.pyx":214 /* "squish/energy.pyx":211
* *
* xi.cache.avg_radius(&xi, avg_radii[i]/(2*PI)) * xi.cache.avg_radius(&xi, avg_radii[i]/(2*PI))
* xi.cache.energy(&xi, site_energy[i]) # <<<<<<<<<<<<<< * xi.cache.energy(&xi, site_energy[i]) # <<<<<<<<<<<<<<
* *
* self.energy = np.sum(site_energy[:self.n]) * self.energy = np.sum(site_energy[:self.n])
*/ */
__pyx_t_15 = __pyx_v_i; __pyx_t_13 = __pyx_v_i;
(void)(__pyx_v_xi.cache->energy((&__pyx_v_xi), (*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_site_energy.data + __pyx_t_15 * __pyx_v_site_energy.strides[0]) ))))); (void)(__pyx_v_xi.cache->energy((&__pyx_v_xi), (*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_site_energy.data + __pyx_t_13 * __pyx_v_site_energy.strides[0]) )))));
} }
} }
} }
@ -5676,16 +5661,16 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_precompute(struct __pyx_obj_
} }
} }
/* "squish/energy.pyx":216 /* "squish/energy.pyx":213
* xi.cache.energy(&xi, site_energy[i]) * xi.cache.energy(&xi, site_energy[i])
* *
* self.energy = np.sum(site_energy[:self.n]) # <<<<<<<<<<<<<< * self.energy = np.sum(site_energy[:self.n]) # <<<<<<<<<<<<<<
* *
* *
*/ */
__Pyx_GetModuleGlobalName(__pyx_t_7, __pyx_n_s_np); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 216, __pyx_L1_error) __Pyx_GetModuleGlobalName(__pyx_t_7, __pyx_n_s_np); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 213, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_7); __Pyx_GOTREF(__pyx_t_7);
__pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_7, __pyx_n_s_sum); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 216, __pyx_L1_error) __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_7, __pyx_n_s_sum); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 213, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3); __Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0; __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
__pyx_t_8.data = __pyx_v_site_energy.data; __pyx_t_8.data = __pyx_v_site_energy.data;
@ -5706,10 +5691,10 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_precompute(struct __pyx_obj_
0, 0,
1) < 0)) 1) < 0))
{ {
__PYX_ERR(0, 216, __pyx_L1_error) __PYX_ERR(0, 213, __pyx_L1_error)
} }
__pyx_t_7 = __pyx_memoryview_fromslice(__pyx_t_8, 1, (PyObject *(*)(char *)) __pyx_memview_get_nn___pyx_t_6squish_4core_FLOAT_T, (int (*)(char *, PyObject *)) __pyx_memview_set_nn___pyx_t_6squish_4core_FLOAT_T, 0);; if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 216, __pyx_L1_error) __pyx_t_7 = __pyx_memoryview_fromslice(__pyx_t_8, 1, (PyObject *(*)(char *)) __pyx_memview_get_nn___pyx_t_6squish_4core_FLOAT_T, (int (*)(char *, PyObject *)) __pyx_memview_set_nn___pyx_t_6squish_4core_FLOAT_T, 0);; if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 213, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_7); __Pyx_GOTREF(__pyx_t_7);
__PYX_XDEC_MEMVIEW(&__pyx_t_8, 1); __PYX_XDEC_MEMVIEW(&__pyx_t_8, 1);
__pyx_t_8.memview = NULL; __pyx_t_8.memview = NULL;
@ -5727,12 +5712,12 @@ __pyx_t_7 = __pyx_memoryview_fromslice(__pyx_t_8, 1, (PyObject *(*)(char *)) __p
__pyx_t_1 = (__pyx_t_4) ? __Pyx_PyObject_Call2Args(__pyx_t_3, __pyx_t_4, __pyx_t_7) : __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_7); __pyx_t_1 = (__pyx_t_4) ? __Pyx_PyObject_Call2Args(__pyx_t_3, __pyx_t_4, __pyx_t_7) : __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_7);
__Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0; __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0; __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 216, __pyx_L1_error) if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 213, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1); __Pyx_GOTREF(__pyx_t_1);
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0; __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_17 = __pyx_PyFloat_AsDouble(__pyx_t_1); if (unlikely((__pyx_t_17 == ((npy_float64)-1)) && PyErr_Occurred())) __PYX_ERR(0, 216, __pyx_L1_error) __pyx_t_15 = __pyx_PyFloat_AsDouble(__pyx_t_1); if (unlikely((__pyx_t_15 == ((npy_float64)-1)) && PyErr_Occurred())) __PYX_ERR(0, 213, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_v_self->__pyx_base.energy = __pyx_t_17; __pyx_v_self->__pyx_base.energy = __pyx_t_15;
/* "squish/energy.pyx":172 /* "squish/energy.pyx":172
* *
@ -5759,7 +5744,7 @@ __pyx_t_7 = __pyx_memoryview_fromslice(__pyx_t_8, 1, (PyObject *(*)(char *)) __p
__Pyx_RefNannyFinishContext(); __Pyx_RefNannyFinishContext();
} }
/* "squish/energy.pyx":219 /* "squish/energy.pyx":216
* *
* *
* cdef void calc_grad(self) except *: # <<<<<<<<<<<<<< * cdef void calc_grad(self) except *: # <<<<<<<<<<<<<<
@ -5771,15 +5756,7 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_grad(struct __pyx_obj_6
__pyx_t_6squish_7voronoi_VoronoiInfo __pyx_v_info; __pyx_t_6squish_7voronoi_VoronoiInfo __pyx_v_info;
__pyx_t_6squish_7voronoi_Site __pyx_v_xi; __pyx_t_6squish_7voronoi_Site __pyx_v_xi;
__pyx_t_6squish_7voronoi_HalfEdge __pyx_v_e; __pyx_t_6squish_7voronoi_HalfEdge __pyx_v_e;
__pyx_t_6squish_7voronoi_HalfEdge __pyx_v_fm;
__pyx_t_6squish_7voronoi_HalfEdge __pyx_v_f;
__pyx_t_6squish_7voronoi_HalfEdge __pyx_v_fp;
__pyx_t_6squish_4core_Vector2D __pyx_v_temp1;
__pyx_t_6squish_4core_Vector2D __pyx_v_temp2;
__pyx_t_6squish_4core_Vector2D __pyx_v_temp3;
__pyx_t_6squish_4core_Vector2D __pyx_v_dedxi_p; __pyx_t_6squish_4core_Vector2D __pyx_v_dedxi_p;
__pyx_t_6squish_4core_BitSet __pyx_v_edge_set;
__pyx_t_6squish_4core_INT_T __pyx_v_num_edges;
__Pyx_memviewslice __pyx_v_dedx = { 0, 0, { 0 }, { 0 }, { 0 } }; __Pyx_memviewslice __pyx_v_dedx = { 0, 0, { 0 }, { 0 }, { 0 } };
__pyx_t_6squish_4core_INT_T __pyx_v_i; __pyx_t_6squish_4core_INT_T __pyx_v_i;
__pyx_t_6squish_4core_INT_T __pyx_v_j; __pyx_t_6squish_4core_INT_T __pyx_v_j;
@ -5793,38 +5770,36 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_grad(struct __pyx_obj_6
__pyx_t_6squish_4core_INT_T __pyx_t_7; __pyx_t_6squish_4core_INT_T __pyx_t_7;
__pyx_t_6squish_4core_INT_T __pyx_t_8; __pyx_t_6squish_4core_INT_T __pyx_t_8;
int __pyx_t_9; int __pyx_t_9;
__pyx_t_6squish_7voronoi_HalfEdge __pyx_t_10; __pyx_t_6squish_4core_INT_T __pyx_t_10;
__pyx_t_6squish_7voronoi_HalfEdge __pyx_t_11; Py_ssize_t __pyx_t_11;
__pyx_t_6squish_4core_INT_T __pyx_t_12;
Py_ssize_t __pyx_t_13;
int __pyx_lineno = 0; int __pyx_lineno = 0;
const char *__pyx_filename = NULL; const char *__pyx_filename = NULL;
int __pyx_clineno = 0; int __pyx_clineno = 0;
__Pyx_RefNannySetupContext("calc_grad", 0); __Pyx_RefNannySetupContext("calc_grad", 0);
/* "squish/energy.pyx":220 /* "squish/energy.pyx":217
* *
* cdef void calc_grad(self) except *: * cdef void calc_grad(self) except *:
* cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points, # <<<<<<<<<<<<<< * cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points, # <<<<<<<<<<<<<<
* self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map) * self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
* *
*/ */
if (unlikely(!__pyx_v_self->__pyx_base.sites.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 220, __pyx_L1_error)} if (unlikely(!__pyx_v_self->__pyx_base.sites.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 217, __pyx_L1_error)}
if (unlikely(!__pyx_v_self->__pyx_base.edges.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 220, __pyx_L1_error)} if (unlikely(!__pyx_v_self->__pyx_base.edges.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 217, __pyx_L1_error)}
if (unlikely(!__pyx_v_self->__pyx_base.points.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 220, __pyx_L1_error)} if (unlikely(!__pyx_v_self->__pyx_base.points.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 217, __pyx_L1_error)}
/* "squish/energy.pyx":221 /* "squish/energy.pyx":218
* cdef void calc_grad(self) except *: * cdef void calc_grad(self) except *:
* cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points, * cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
* self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map) # <<<<<<<<<<<<<< * self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map) # <<<<<<<<<<<<<<
* *
* cdef Site xi * cdef Site xi
*/ */
if (unlikely(!__pyx_v_self->__pyx_base.vertices.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 221, __pyx_L1_error)} if (unlikely(!__pyx_v_self->__pyx_base.vertices.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 218, __pyx_L1_error)}
if (unlikely(!__pyx_v_self->__pyx_base.site_cache.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 221, __pyx_L1_error)} if (unlikely(!__pyx_v_self->__pyx_base.site_cache.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 218, __pyx_L1_error)}
if (unlikely(!__pyx_v_self->__pyx_base.edge_cache.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 221, __pyx_L1_error)} if (unlikely(!__pyx_v_self->__pyx_base.edge_cache.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 218, __pyx_L1_error)}
/* "squish/energy.pyx":220 /* "squish/energy.pyx":217
* *
* cdef void calc_grad(self) except *: * cdef void calc_grad(self) except *:
* cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points, # <<<<<<<<<<<<<< * cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points, # <<<<<<<<<<<<<<
@ -5833,31 +5808,21 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_grad(struct __pyx_obj_6
*/ */
__pyx_v_info = __pyx_f_6squish_7voronoi__VoronoiInfo(__pyx_v_self->__pyx_base.sites, __pyx_v_self->__pyx_base.edges, __pyx_v_self->__pyx_base.points, __pyx_v_self->__pyx_base.vertices, __pyx_v_self->__pyx_base.site_cache, __pyx_v_self->__pyx_base.edge_cache, __pyx_v_self->__pyx_base.edge_cache_map); __pyx_v_info = __pyx_f_6squish_7voronoi__VoronoiInfo(__pyx_v_self->__pyx_base.sites, __pyx_v_self->__pyx_base.edges, __pyx_v_self->__pyx_base.points, __pyx_v_self->__pyx_base.vertices, __pyx_v_self->__pyx_base.site_cache, __pyx_v_self->__pyx_base.edge_cache, __pyx_v_self->__pyx_base.edge_cache_map);
/* "squish/energy.pyx":228 /* "squish/energy.pyx":224
* cdef BitSet edge_set * cdef Vector2D dedxi_p
*
* cdef INT_T num_edges = self.edges.shape[0] # <<<<<<<<<<<<<<
*
* cdef FLOAT_T [:, ::1] dedx = np.zeros((self.n, 2), dtype=FLOAT)
*/
if (unlikely(!__pyx_v_self->__pyx_base.edges.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 228, __pyx_L1_error)}
__pyx_v_num_edges = (__pyx_v_self->__pyx_base.edges.shape[0]);
/* "squish/energy.pyx":230
* cdef INT_T num_edges = self.edges.shape[0]
* *
* cdef FLOAT_T [:, ::1] dedx = np.zeros((self.n, 2), dtype=FLOAT) # <<<<<<<<<<<<<< * cdef FLOAT_T [:, ::1] dedx = np.zeros((self.n, 2), dtype=FLOAT) # <<<<<<<<<<<<<<
* *
* cdef INT_T i, j * cdef INT_T i, j
*/ */
__Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 230, __pyx_L1_error) __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 224, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1); __Pyx_GOTREF(__pyx_t_1);
__pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 230, __pyx_L1_error) __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 224, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_2); __Pyx_GOTREF(__pyx_t_2);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__pyx_t_1 = __Pyx_PyInt_From_npy_int64(__pyx_v_self->__pyx_base.n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 230, __pyx_L1_error) __pyx_t_1 = __Pyx_PyInt_From_npy_int64(__pyx_v_self->__pyx_base.n); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 224, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1); __Pyx_GOTREF(__pyx_t_1);
__pyx_t_3 = PyTuple_New(2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 230, __pyx_L1_error) __pyx_t_3 = PyTuple_New(2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 224, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3); __Pyx_GOTREF(__pyx_t_3);
__Pyx_GIVEREF(__pyx_t_1); __Pyx_GIVEREF(__pyx_t_1);
PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_1); PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_1);
@ -5865,29 +5830,29 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_grad(struct __pyx_obj_6
__Pyx_GIVEREF(__pyx_int_2); __Pyx_GIVEREF(__pyx_int_2);
PyTuple_SET_ITEM(__pyx_t_3, 1, __pyx_int_2); PyTuple_SET_ITEM(__pyx_t_3, 1, __pyx_int_2);
__pyx_t_1 = 0; __pyx_t_1 = 0;
__pyx_t_1 = PyTuple_New(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 230, __pyx_L1_error) __pyx_t_1 = PyTuple_New(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 224, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1); __Pyx_GOTREF(__pyx_t_1);
__Pyx_GIVEREF(__pyx_t_3); __Pyx_GIVEREF(__pyx_t_3);
PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_t_3); PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_t_3);
__pyx_t_3 = 0; __pyx_t_3 = 0;
__pyx_t_3 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 230, __pyx_L1_error) __pyx_t_3 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 224, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3); __Pyx_GOTREF(__pyx_t_3);
__Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_n_s_FLOAT); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 230, __pyx_L1_error) __Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_n_s_FLOAT); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 224, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4); __Pyx_GOTREF(__pyx_t_4);
if (PyDict_SetItem(__pyx_t_3, __pyx_n_s_dtype, __pyx_t_4) < 0) __PYX_ERR(0, 230, __pyx_L1_error) if (PyDict_SetItem(__pyx_t_3, __pyx_n_s_dtype, __pyx_t_4) < 0) __PYX_ERR(0, 224, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_t_4 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_1, __pyx_t_3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 230, __pyx_L1_error) __pyx_t_4 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_1, __pyx_t_3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 224, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4); __Pyx_GOTREF(__pyx_t_4);
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0; __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0; __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__pyx_t_5 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_nn___pyx_t_6squish_4core_FLOAT_T(__pyx_t_4, PyBUF_WRITABLE); if (unlikely(!__pyx_t_5.memview)) __PYX_ERR(0, 230, __pyx_L1_error) __pyx_t_5 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_nn___pyx_t_6squish_4core_FLOAT_T(__pyx_t_4, PyBUF_WRITABLE); if (unlikely(!__pyx_t_5.memview)) __PYX_ERR(0, 224, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_v_dedx = __pyx_t_5; __pyx_v_dedx = __pyx_t_5;
__pyx_t_5.memview = NULL; __pyx_t_5.memview = NULL;
__pyx_t_5.data = NULL; __pyx_t_5.data = NULL;
/* "squish/energy.pyx":233 /* "squish/energy.pyx":227
* *
* cdef INT_T i, j * cdef INT_T i, j
* for i in prange(self.n, nogil=True): # <<<<<<<<<<<<<< * for i in prange(self.n, nogil=True): # <<<<<<<<<<<<<<
@ -5914,11 +5879,11 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_grad(struct __pyx_obj_6
if (__pyx_t_8 > 0) if (__pyx_t_8 > 0)
{ {
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp parallel private(__pyx_t_10, __pyx_t_11, __pyx_t_12, __pyx_t_13, __pyx_t_9) #pragma omp parallel private(__pyx_t_10, __pyx_t_11, __pyx_t_9)
#endif /* _OPENMP */ #endif /* _OPENMP */
{ {
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp for lastprivate(__pyx_v_dedxi_p) lastprivate(__pyx_v_e) lastprivate(__pyx_v_edge_set) lastprivate(__pyx_v_f) lastprivate(__pyx_v_fm) lastprivate(__pyx_v_fp) firstprivate(__pyx_v_i) lastprivate(__pyx_v_i) lastprivate(__pyx_v_j) lastprivate(__pyx_v_temp1) lastprivate(__pyx_v_temp2) lastprivate(__pyx_v_temp3) lastprivate(__pyx_v_xi) #pragma omp for lastprivate(__pyx_v_dedxi_p) lastprivate(__pyx_v_e) firstprivate(__pyx_v_i) lastprivate(__pyx_v_i) lastprivate(__pyx_v_j) lastprivate(__pyx_v_xi)
#endif /* _OPENMP */ #endif /* _OPENMP */
for (__pyx_t_7 = 0; __pyx_t_7 < __pyx_t_8; __pyx_t_7++){ for (__pyx_t_7 = 0; __pyx_t_7 < __pyx_t_8; __pyx_t_7++){
{ {
@ -5926,363 +5891,86 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_grad(struct __pyx_obj_6
/* Initialize private variables to invalid values */ /* Initialize private variables to invalid values */
__pyx_v_j = ((__pyx_t_6squish_4core_INT_T)0xbad0bad0); __pyx_v_j = ((__pyx_t_6squish_4core_INT_T)0xbad0bad0);
/* "squish/energy.pyx":234 /* "squish/energy.pyx":228
* cdef INT_T i, j * cdef INT_T i, j
* for i in prange(self.n, nogil=True): * for i in prange(self.n, nogil=True):
* xi = _Site(i, &info) # <<<<<<<<<<<<<< * xi = _Site(i, &info) # <<<<<<<<<<<<<<
* e = xi.edge(&xi) * e = xi.edge(&xi)
* edge_set = _BitSet(num_edges) *
*/ */
__pyx_v_xi = __pyx_f_6squish_7voronoi__Site(__pyx_v_i, (&__pyx_v_info)); __pyx_v_xi = __pyx_f_6squish_7voronoi__Site(__pyx_v_i, (&__pyx_v_info));
/* "squish/energy.pyx":235 /* "squish/energy.pyx":229
* for i in prange(self.n, nogil=True): * for i in prange(self.n, nogil=True):
* xi = _Site(i, &info) * xi = _Site(i, &info)
* e = xi.edge(&xi) # <<<<<<<<<<<<<< * e = xi.edge(&xi) # <<<<<<<<<<<<<<
* edge_set = _BitSet(num_edges)
* *
* j = 0
*/ */
__pyx_v_e = __pyx_v_xi.edge((&__pyx_v_xi)); __pyx_v_e = __pyx_v_xi.edge((&__pyx_v_xi));
/* "squish/energy.pyx":236 /* "squish/energy.pyx":231
* xi = _Site(i, &info)
* e = xi.edge(&xi) * e = xi.edge(&xi)
* edge_set = _BitSet(num_edges) # <<<<<<<<<<<<<<
*
* j = 0
*/
__pyx_v_edge_set = __pyx_f_6squish_4core__BitSet(__pyx_v_num_edges);
/* "squish/energy.pyx":238
* edge_set = _BitSet(num_edges)
* *
* j = 0 # <<<<<<<<<<<<<< * j = 0 # <<<<<<<<<<<<<<
* while j < xi.edge_num(&xi): # Looping through site edges. * while j < xi.edge_num(&xi): # Looping through site edges.
* f = e * dedxi_p = e.cache.ya(&e, NAN_VECTOR)
*/ */
__pyx_v_j = 0; __pyx_v_j = 0;
/* "squish/energy.pyx":239 /* "squish/energy.pyx":232
* *
* j = 0 * j = 0
* while j < xi.edge_num(&xi): # Looping through site edges. # <<<<<<<<<<<<<< * while j < xi.edge_num(&xi): # Looping through site edges. # <<<<<<<<<<<<<<
* f = e * dedxi_p = e.cache.ya(&e, NAN_VECTOR)
* dedxi_p = _Vector2D(0, 0) * dedxi_p.self.smul(
*/ */
while (1) { while (1) {
__pyx_t_9 = ((__pyx_v_j < __pyx_v_xi.edge_num((&__pyx_v_xi))) != 0); __pyx_t_9 = ((__pyx_v_j < __pyx_v_xi.edge_num((&__pyx_v_xi))) != 0);
if (!__pyx_t_9) break; if (!__pyx_t_9) break;
/* "squish/energy.pyx":240 /* "squish/energy.pyx":233
* j = 0 * j = 0
* while j < xi.edge_num(&xi): # Looping through site edges. * while j < xi.edge_num(&xi): # Looping through site edges.
* f = e # <<<<<<<<<<<<<< * dedxi_p = e.cache.ya(&e, NAN_VECTOR) # <<<<<<<<<<<<<<
* dedxi_p = _Vector2D(0, 0) * dedxi_p.self.smul(
* # dedx (only x) * &dedxi_p,
*/ */
__pyx_v_f = __pyx_v_e; __pyx_v_dedxi_p = __pyx_v_e.cache->ya((&__pyx_v_e), __pyx_v_6squish_7voronoi_NAN_VECTOR);
/* "squish/energy.pyx":241 /* "squish/energy.pyx":234
* while j < xi.edge_num(&xi): # Looping through site edges. * while j < xi.edge_num(&xi): # Looping through site edges.
* f = e * dedxi_p = e.cache.ya(&e, NAN_VECTOR)
* dedxi_p = _Vector2D(0, 0) # <<<<<<<<<<<<<< * dedxi_p.self.smul( # <<<<<<<<<<<<<<
* # dedx (only x) * &dedxi_p,
* temp1 = f.cache.la(&f, NAN_VECTOR) * e.cache.calI(&e, NAN) / e.cache.ya_mag(&e, NAN)
*/ */
__pyx_v_dedxi_p = __pyx_f_6squish_4core__Vector2D(0.0, 0.0); (void)(__pyx_v_dedxi_p.self.smul((&__pyx_v_dedxi_p), (__pyx_v_e.cache->calI((&__pyx_v_e), NAN) / __pyx_v_e.cache->ya_mag((&__pyx_v_e), NAN))));
/* "squish/energy.pyx":243 /* "squish/energy.pyx":239
* dedxi_p = _Vector2D(0, 0)
* # dedx (only x)
* temp1 = f.cache.la(&f, NAN_VECTOR) # <<<<<<<<<<<<<<
* temp1 = temp1.rot(&temp1)
* dedxi_p = temp1.copy.smul(
*/
__pyx_v_temp1 = __pyx_v_f.cache->la((&__pyx_v_f), __pyx_v_6squish_7voronoi_NAN_VECTOR);
/* "squish/energy.pyx":244
* # dedx (only x)
* temp1 = f.cache.la(&f, NAN_VECTOR)
* temp1 = temp1.rot(&temp1) # <<<<<<<<<<<<<<
* dedxi_p = temp1.copy.smul(
* &temp1,
*/
__pyx_v_temp1 = __pyx_v_temp1.rot((&__pyx_v_temp1));
/* "squish/energy.pyx":245
* 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)
*/
__pyx_v_dedxi_p = __pyx_v_temp1.copy.smul((&__pyx_v_temp1), (__pyx_v_f.cache->calI((&__pyx_v_f), NAN) / __pyx_v_f.cache->la_mag((&__pyx_v_f), NAN)));
/* "squish/energy.pyx":250
* ) * )
* *
* while True: # Circling this vertex. # <<<<<<<<<<<<<< * dedx[i][0] += 2*self.r*dedxi_p.x # <<<<<<<<<<<<<<
* fm, fp = f.prev(&f), f.next(&f) * dedx[i][1] += 2*self.r*dedxi_p.y
* *
*/ */
while (1) { __pyx_t_10 = __pyx_v_i;
__pyx_t_11 = 0;
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ ((char *) (((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_dedx.data + __pyx_t_10 * __pyx_v_dedx.strides[0]) )) + __pyx_t_11)) )) += ((2.0 * __pyx_v_self->__pyx_base.r) * __pyx_v_dedxi_p.x);
/* "squish/energy.pyx":251 /* "squish/energy.pyx":240
* *
* while True: # Circling this vertex. * dedx[i][0] += 2*self.r*dedxi_p.x
* fm, fp = f.prev(&f), f.next(&f) # <<<<<<<<<<<<<< * dedx[i][1] += 2*self.r*dedxi_p.y # <<<<<<<<<<<<<<
*
* if not edge_set.add(&edge_set, f.arr_index):
*/
__pyx_t_10 = __pyx_v_f.prev((&__pyx_v_f));
__pyx_t_11 = __pyx_v_f.next((&__pyx_v_f));
__pyx_v_fm = __pyx_t_10;
__pyx_v_fp = __pyx_t_11;
/* "squish/energy.pyx":253
* 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)
*/
__pyx_t_9 = ((!(__pyx_v_edge_set.add((&__pyx_v_edge_set), __pyx_v_f.arr_index) != 0)) != 0);
if (__pyx_t_9) {
/* "squish/energy.pyx":255
* 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))
*/
__pyx_v_temp1 = __pyx_v_fp.cache->da((&__pyx_v_fp), __pyx_v_6squish_7voronoi_NAN_VECTOR);
/* "squish/energy.pyx":256
* # ( -rot(dap) ) / |la|
* temp1 = fp.cache.da(&fp, NAN_VECTOR)
* temp1 = temp1.rot(&temp1) # <<<<<<<<<<<<<<
* temp1.self.sdiv(&temp1, -f.cache.la_mag(&f, NAN))
*
*/
__pyx_v_temp1 = __pyx_v_temp1.rot((&__pyx_v_temp1));
/* "squish/energy.pyx":257
* 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
*/
(void)(__pyx_v_temp1.self.sdiv((&__pyx_v_temp1), (-__pyx_v_f.cache->la_mag((&__pyx_v_f), NAN))));
/* "squish/energy.pyx":260
*
* # la * area_p / |la|^2
* temp3 = f.cache.la(&f, NAN_VECTOR) # <<<<<<<<<<<<<<
* temp3.self.smul(
* &temp3,
*/
__pyx_v_temp3 = __pyx_v_f.cache->la((&__pyx_v_f), __pyx_v_6squish_7voronoi_NAN_VECTOR);
/* "squish/energy.pyx":261
* # 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
*/
(void)(__pyx_v_temp3.self.smul((&__pyx_v_temp3), (__pyx_v_f.cache->area_p((&__pyx_v_f), NAN) / pow(__pyx_v_f.cache->la_mag((&__pyx_v_f), NAN), 3.0))));
/* "squish/energy.pyx":266
* )
* # Combine * calI
* temp1.self.vadd(&temp1, temp3) # <<<<<<<<<<<<<<
* temp1.self.smul(&temp1, f.cache.calI(&f, NAN))
*
*/
(void)(__pyx_v_temp1.self.vadd((&__pyx_v_temp1), __pyx_v_temp3));
/* "squish/energy.pyx":267
* # Combine * calI
* temp1.self.vadd(&temp1, temp3)
* temp1.self.smul(&temp1, f.cache.calI(&f, NAN)) # <<<<<<<<<<<<<<
*
* # rot(dam) / |lam|
*/
(void)(__pyx_v_temp1.self.smul((&__pyx_v_temp1), __pyx_v_f.cache->calI((&__pyx_v_f), NAN)));
/* "squish/energy.pyx":270
*
* # rot(dam) / |lam|
* temp2 = fm.cache.da(&fm, NAN_VECTOR) # <<<<<<<<<<<<<<
* temp2 = temp2.rot(&temp2)
* temp2.self.sdiv(&temp2, fm.cache.la_mag(&fm, NAN))
*/
__pyx_v_temp2 = __pyx_v_fm.cache->da((&__pyx_v_fm), __pyx_v_6squish_7voronoi_NAN_VECTOR);
/* "squish/energy.pyx":271
* # rot(dam) / |lam|
* temp2 = fm.cache.da(&fm, NAN_VECTOR)
* temp2 = temp2.rot(&temp2) # <<<<<<<<<<<<<<
* temp2.self.sdiv(&temp2, fm.cache.la_mag(&fm, NAN))
*
*/
__pyx_v_temp2 = __pyx_v_temp2.rot((&__pyx_v_temp2));
/* "squish/energy.pyx":272
* 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
*/
(void)(__pyx_v_temp2.self.sdiv((&__pyx_v_temp2), __pyx_v_fm.cache->la_mag((&__pyx_v_fm), NAN)));
/* "squish/energy.pyx":275
*
* # lam * area_pm / |lam|^2
* temp3 = fm.cache.la(&fm, NAN_VECTOR) # <<<<<<<<<<<<<<
* temp3.self.smul(
* &temp3,
*/
__pyx_v_temp3 = __pyx_v_fm.cache->la((&__pyx_v_fm), __pyx_v_6squish_7voronoi_NAN_VECTOR);
/* "squish/energy.pyx":276
* # 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
*/
(void)(__pyx_v_temp3.self.smul((&__pyx_v_temp3), (__pyx_v_fm.cache->area_p((&__pyx_v_fm), NAN) / pow(__pyx_v_fm.cache->la_mag((&__pyx_v_fm), NAN), 3.0))));
/* "squish/energy.pyx":281
* )
* # Combine * calIm
* temp2.self.vsub(&temp2, temp3) # <<<<<<<<<<<<<<
* temp2.self.smul(&temp2, fm.cache.calI(&fm, NAN))
*
*/
(void)(__pyx_v_temp2.self.vsub((&__pyx_v_temp2), __pyx_v_temp3));
/* "squish/energy.pyx":282
* # Combine * calIm
* temp2.self.vsub(&temp2, temp3)
* temp2.self.smul(&temp2, fm.cache.calI(&fm, NAN)) # <<<<<<<<<<<<<<
*
* temp1.self.vadd(&temp1, temp2)
*/
(void)(__pyx_v_temp2.self.smul((&__pyx_v_temp2), __pyx_v_fm.cache->calI((&__pyx_v_fm), NAN)));
/* "squish/energy.pyx":284
* temp2.self.smul(&temp2, fm.cache.calI(&fm, NAN))
*
* temp1.self.vadd(&temp1, temp2) # <<<<<<<<<<<<<<
* temp1.self.matmul(&temp1, e.get_H(&e, xi))
*
*/
(void)(__pyx_v_temp1.self.vadd((&__pyx_v_temp1), __pyx_v_temp2));
/* "squish/energy.pyx":285
*
* temp1.self.vadd(&temp1, temp2)
* temp1.self.matmul(&temp1, e.get_H(&e, xi)) # <<<<<<<<<<<<<<
*
* dedxi_p.self.vadd(&dedxi_p, temp1)
*/
(void)(__pyx_v_temp1.self.matmul((&__pyx_v_temp1), __pyx_v_e.get_H((&__pyx_v_e), __pyx_v_xi)));
/* "squish/energy.pyx":287
* temp1.self.matmul(&temp1, e.get_H(&e, xi))
*
* dedxi_p.self.vadd(&dedxi_p, temp1) # <<<<<<<<<<<<<<
*
* f = f.twin(&f)
*/
(void)(__pyx_v_dedxi_p.self.vadd((&__pyx_v_dedxi_p), __pyx_v_temp1));
/* "squish/energy.pyx":253
* 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)
*/
}
/* "squish/energy.pyx":289
* dedxi_p.self.vadd(&dedxi_p, temp1)
*
* f = f.twin(&f) # <<<<<<<<<<<<<<
* f = f.next(&f)
*
*/
__pyx_v_f = __pyx_v_f.twin((&__pyx_v_f));
/* "squish/energy.pyx":290
*
* f = f.twin(&f)
* f = f.next(&f) # <<<<<<<<<<<<<<
*
* if f.arr_index == e.arr_index:
*/
__pyx_v_f = __pyx_v_f.next((&__pyx_v_f));
/* "squish/energy.pyx":292
* f = f.next(&f)
*
* if f.arr_index == e.arr_index: # <<<<<<<<<<<<<<
* break
*
*/
__pyx_t_9 = ((__pyx_v_f.arr_index == __pyx_v_e.arr_index) != 0);
if (__pyx_t_9) {
/* "squish/energy.pyx":293
*
* if f.arr_index == e.arr_index:
* break # <<<<<<<<<<<<<<
*
* dedx[i][0] -= -2*self.r*dedxi_p.x
*/
goto __pyx_L13_break;
/* "squish/energy.pyx":292
* f = f.next(&f)
*
* if f.arr_index == e.arr_index: # <<<<<<<<<<<<<<
* break
*
*/
}
}
__pyx_L13_break:;
/* "squish/energy.pyx":295
* break
*
* dedx[i][0] -= -2*self.r*dedxi_p.x # <<<<<<<<<<<<<<
* dedx[i][1] -= -2*self.r*dedxi_p.y
*
*/
__pyx_t_12 = __pyx_v_i;
__pyx_t_13 = 0;
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ ((char *) (((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_dedx.data + __pyx_t_12 * __pyx_v_dedx.strides[0]) )) + __pyx_t_13)) )) -= ((-2.0 * __pyx_v_self->__pyx_base.r) * __pyx_v_dedxi_p.x);
/* "squish/energy.pyx":296
*
* dedx[i][0] -= -2*self.r*dedxi_p.x
* dedx[i][1] -= -2*self.r*dedxi_p.y # <<<<<<<<<<<<<<
* *
* e = e.next(&e) * e = e.next(&e)
*/ */
__pyx_t_12 = __pyx_v_i; __pyx_t_10 = __pyx_v_i;
__pyx_t_13 = 1; __pyx_t_11 = 1;
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ ((char *) (((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_dedx.data + __pyx_t_12 * __pyx_v_dedx.strides[0]) )) + __pyx_t_13)) )) -= ((-2.0 * __pyx_v_self->__pyx_base.r) * __pyx_v_dedxi_p.y); *((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ ((char *) (((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=0 */ (__pyx_v_dedx.data + __pyx_t_10 * __pyx_v_dedx.strides[0]) )) + __pyx_t_11)) )) += ((2.0 * __pyx_v_self->__pyx_base.r) * __pyx_v_dedxi_p.y);
/* "squish/energy.pyx":298 /* "squish/energy.pyx":242
* dedx[i][1] -= -2*self.r*dedxi_p.y * dedx[i][1] += 2*self.r*dedxi_p.y
* *
* e = e.next(&e) # <<<<<<<<<<<<<< * e = e.next(&e) # <<<<<<<<<<<<<<
* j = j + 1 * j = j + 1
@ -6290,23 +5978,15 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_grad(struct __pyx_obj_6
*/ */
__pyx_v_e = __pyx_v_e.next((&__pyx_v_e)); __pyx_v_e = __pyx_v_e.next((&__pyx_v_e));
/* "squish/energy.pyx":299 /* "squish/energy.pyx":243
* *
* e = e.next(&e) * e = e.next(&e)
* j = j + 1 # <<<<<<<<<<<<<< * j = j + 1 # <<<<<<<<<<<<<<
* *
* edge_set.free(&edge_set) * self.grad = dedx
*/ */
__pyx_v_j = (__pyx_v_j + 1); __pyx_v_j = (__pyx_v_j + 1);
} }
/* "squish/energy.pyx":301
* j = j + 1
*
* edge_set.free(&edge_set) # <<<<<<<<<<<<<<
* self.grad = dedx
*/
__pyx_v_edge_set.free((&__pyx_v_edge_set));
} }
} }
} }
@ -6320,7 +6000,7 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_grad(struct __pyx_obj_6
#endif #endif
} }
/* "squish/energy.pyx":233 /* "squish/energy.pyx":227
* *
* cdef INT_T i, j * cdef INT_T i, j
* for i in prange(self.n, nogil=True): # <<<<<<<<<<<<<< * for i in prange(self.n, nogil=True): # <<<<<<<<<<<<<<
@ -6339,16 +6019,16 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_grad(struct __pyx_obj_6
} }
} }
/* "squish/energy.pyx":302 /* "squish/energy.pyx":245
* j = j + 1
* *
* edge_set.free(&edge_set)
* self.grad = dedx # <<<<<<<<<<<<<< * self.grad = dedx # <<<<<<<<<<<<<<
*/ */
__PYX_XDEC_MEMVIEW(&__pyx_v_self->__pyx_base.grad, 0); __PYX_XDEC_MEMVIEW(&__pyx_v_self->__pyx_base.grad, 0);
__PYX_INC_MEMVIEW(&__pyx_v_dedx, 0); __PYX_INC_MEMVIEW(&__pyx_v_dedx, 0);
__pyx_v_self->__pyx_base.grad = __pyx_v_dedx; __pyx_v_self->__pyx_base.grad = __pyx_v_dedx;
/* "squish/energy.pyx":219 /* "squish/energy.pyx":216
* *
* *
* cdef void calc_grad(self) except *: # <<<<<<<<<<<<<< * cdef void calc_grad(self) except *: # <<<<<<<<<<<<<<

85
squish/energy.pyx
View File

@ -99,8 +99,8 @@ cdef class AreaEnergy(VoronoiContainer):
dedxi_p = f.cache.dVdv(&f, NAN_VECTOR) #dVdv dedxi_p = f.cache.dVdv(&f, NAN_VECTOR) #dVdv
dedxi_p.self.smul(&dedxi_p, xf.cache.area(&xf, NAN) - A) dedxi_p.self.smul(&dedxi_p, xf.cache.area(&xf, NAN) - A)
dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX)) dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX))
dedx[i][0] -= dedxi_p.x dedx[i][0] += dedxi_p.x
dedx[i][1] -= dedxi_p.y dedx[i][1] += dedxi_p.y
f = f.twin(&f) f = f.twin(&f)
f = f.next(&f) f = f.next(&f)
@ -174,7 +174,7 @@ cdef class RadialTEnergy(VoronoiContainer):
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map) self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
cdef Site xi cdef Site xi
cdef HalfEdge e, em, ep cdef HalfEdge e, ep
cdef Vector2D la cdef Vector2D la
# All energy has a 2pir_0 term. # All energy has a 2pir_0 term.
@ -188,8 +188,8 @@ cdef class RadialTEnergy(VoronoiContainer):
e = xi.edge(&xi) e = xi.edge(&xi)
j = 0 j = 0
while j < xi.edge_num(&xi): while j < xi.edge_num(&xi):
em, ep = e.prev(&e), e.next(&e) ep = e.next(&e)
e.cache.H(&e, VoronoiContainer.calc_H(em, e)) #e.cache.H(&e, VoronoiContainer.calc_H(em, e))
la = e.cache.la(&e, NAN_VECTOR) la = e.cache.la(&e, NAN_VECTOR)
sp = la.dot(&la, ep.cache.da(&ep, NAN_VECTOR)) # dap . la sp = la.dot(&la, ep.cache.da(&ep, NAN_VECTOR)) # dap . la
@ -200,10 +200,7 @@ cdef class RadialTEnergy(VoronoiContainer):
e.cache.calI(&e, <FLOAT_T>(atanh(<double>sp) - atanh(<double>sm))) e.cache.calI(&e, <FLOAT_T>(atanh(<double>sp) - atanh(<double>sm)))
avg_radii[i] += ( avg_radii[i] += e.cache.ya_mag(&e, NAN) * e.cache.calI(&e, NAN) / 2
e.cache.area_p(&e, NAN) * e.cache.calI(&e, NAN)
/ e.cache.la_mag(&e, NAN)
)
e = e.next(&e) e = e.next(&e)
j = j + 1 j = j + 1
@ -221,11 +218,8 @@ cdef class RadialTEnergy(VoronoiContainer):
self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map) self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
cdef Site xi cdef Site xi
cdef HalfEdge e, fm, f, fp cdef HalfEdge e
cdef Vector2D temp1, temp2, temp3, dedxi_p cdef Vector2D 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 FLOAT_T [:, ::1] dedx = np.zeros((self.n, 2), dtype=FLOAT)
@ -233,70 +227,19 @@ cdef class RadialTEnergy(VoronoiContainer):
for i in prange(self.n, nogil=True): for i in prange(self.n, nogil=True):
xi = _Site(i, &info) xi = _Site(i, &info)
e = xi.edge(&xi) e = xi.edge(&xi)
edge_set = _BitSet(num_edges)
j = 0 j = 0
while j < xi.edge_num(&xi): # Looping through site edges. while j < xi.edge_num(&xi): # Looping through site edges.
f = e dedxi_p = e.cache.ya(&e, NAN_VECTOR)
dedxi_p = _Vector2D(0, 0) dedxi_p.self.smul(
# dedx (only x) &dedxi_p,
temp1 = f.cache.la(&f, NAN_VECTOR) e.cache.calI(&e, NAN) / e.cache.ya_mag(&e, NAN)
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. dedx[i][0] += 2*self.r*dedxi_p.x
fm, fp = f.prev(&f), f.next(&f) dedx[i][1] += 2*self.r*dedxi_p.y
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) e = e.next(&e)
j = j + 1 j = j + 1
edge_set.free(&edge_set)
self.grad = dedx self.grad = dedx

2240
squish/voronoi.c
View File

File diff suppressed because it is too large Load Diff

6
squish/voronoi.pxd
View File

@ -16,18 +16,18 @@ ctypedef struct SiteCacheMap:
# Psuedo-class that handles caching for edges. # Psuedo-class that handles caching for edges.
ctypedef struct EdgeCacheMap: ctypedef struct EdgeCacheMap:
INT_T iH, ila, ida, ixij, idVdv, ila_mag, ida_mag, iarea_p, icalI, size INT_T iH, ila, ida, iya, idVdv, ila_mag, ida_mag, iya_mag, icalI, size
Matrix2x2 (*H)(HalfEdge*, Matrix2x2) nogil Matrix2x2 (*H)(HalfEdge*, Matrix2x2) nogil
Vector2D (*la)(HalfEdge*, Vector2D) nogil Vector2D (*la)(HalfEdge*, Vector2D) nogil
Vector2D (*da)(HalfEdge*, Vector2D) nogil Vector2D (*da)(HalfEdge*, Vector2D) nogil
Vector2D (*xij)(HalfEdge*, Vector2D) nogil Vector2D (*ya)(HalfEdge*, Vector2D) nogil
Vector2D (*dVdv)(HalfEdge*, Vector2D) nogil Vector2D (*dVdv)(HalfEdge*, Vector2D) nogil
FLOAT_T (*la_mag)(HalfEdge*, FLOAT_T) nogil FLOAT_T (*la_mag)(HalfEdge*, FLOAT_T) nogil
FLOAT_T (*da_mag)(HalfEdge*, FLOAT_T) nogil FLOAT_T (*da_mag)(HalfEdge*, FLOAT_T) nogil
FLOAT_T (*area_p)(HalfEdge*, FLOAT_T) nogil FLOAT_T (*ya_mag)(HalfEdge*, FLOAT_T) nogil
FLOAT_T (*calI)(HalfEdge*, FLOAT_T) nogil FLOAT_T (*calI)(HalfEdge*, FLOAT_T) nogil
# Psuedo-class to just contain all pertaining info for sites and edges. # Psuedo-class to just contain all pertaining info for sites and edges.

181
squish/voronoi.pyx
View File

@ -129,16 +129,16 @@ cdef inline Vector2D maxcenter(Site* self, Vector2D val) nogil:
#### EdgeCacheMap Methods #### #### EdgeCacheMap Methods ####
cdef inline EdgeCacheMap _EdgeCacheMap(INT_T iH, INT_T ila, INT_T ida, INT_T ixij, cdef inline EdgeCacheMap _EdgeCacheMap(INT_T iH, INT_T ila, INT_T ida, INT_T iya,
INT_T idVdv, INT_T ila_mag, INT_T ida_mag, INT_T idVdv, INT_T ila_mag, INT_T ida_mag,
INT_T iarea_p, INT_T icalI, INT_T size) nogil: INT_T iya_mag, INT_T icalI, INT_T size) nogil:
cdef EdgeCacheMap ec cdef EdgeCacheMap ec
ec.iH, ec.ila, ec.ida, ec.ixij, ec.idVdv = iH, ila, ida, ixij, idVdv ec.iH, ec.ila, ec.ida, ec.iya, ec.idVdv = iH, ila, ida, iya, idVdv
ec.ila_mag, ec.ida_mag, ec.iarea_p, ec.icalI = ila_mag, ida_mag, iarea_p, icalI ec.ila_mag, ec.ida_mag, ec.iya_mag, ec.icalI = ila_mag, ida_mag, iya_mag, icalI
ec.size = size ec.size = size
ec.H, ec.la, ec.da, ec.xij, ec.dVdv = H, la, da, xij, dVdv ec.H, ec.la, ec.da, ec.ya, ec.dVdv = H, la, da, ya, dVdv
ec.la_mag, ec.da_mag, ec.area_p, ec.calI = la_mag, da_mag, area_p, calI ec.la_mag, ec.da_mag, ec.ya_mag, ec.calI = la_mag, da_mag, ya_mag, calI
return ec return ec
@ -203,21 +203,21 @@ cdef inline Vector2D da(HalfEdge* self, Vector2D val) nogil:
(self.arr_index, self.cache.ida+1), val.y) (self.arr_index, self.cache.ida+1), val.y)
return val return val
cdef inline Vector2D xij(HalfEdge* self, Vector2D val) nogil: cdef inline Vector2D ya(HalfEdge* self, Vector2D val) nogil:
if isnan(<double>val.x): if isnan(<double>val.x):
return _Vector2D( return _Vector2D(
self.info.edge_cache.get(&self.info.edge_cache, self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ixij) (self.arr_index, self.cache.iya)
), ),
self.info.edge_cache.get(&self.info.edge_cache, self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.ixij+1) (self.arr_index, self.cache.iya+1)
) )
) )
else: else:
self.info.edge_cache.set(&self.info.edge_cache, self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ixij), val.x) (self.arr_index, self.cache.iya), val.x)
self.info.edge_cache.set(&self.info.edge_cache, self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.ixij+1), val.y) (self.arr_index, self.cache.iya+1), val.y)
return val return val
cdef inline Vector2D dVdv(HalfEdge* self, Vector2D val) nogil: cdef inline Vector2D dVdv(HalfEdge* self, Vector2D val) nogil:
@ -257,14 +257,14 @@ cdef inline FLOAT_T da_mag(HalfEdge* self, FLOAT_T val) nogil:
(self.arr_index, self.cache.ida_mag), val) (self.arr_index, self.cache.ida_mag), val)
return val return val
cdef inline FLOAT_T area_p(HalfEdge* self, FLOAT_T val) nogil: cdef inline FLOAT_T ya_mag(HalfEdge* self, FLOAT_T val) nogil:
if isnan(<double>val): if isnan(<double>val):
return self.info.edge_cache.get(&self.info.edge_cache, return self.info.edge_cache.get(&self.info.edge_cache,
(self.arr_index, self.cache.iarea_p) (self.arr_index, self.cache.iya_mag)
) )
else: else:
self.info.edge_cache.set(&self.info.edge_cache, self.info.edge_cache.set(&self.info.edge_cache,
(self.arr_index, self.cache.iarea_p), val) (self.arr_index, self.cache.iya_mag), val)
return val return val
cdef inline FLOAT_T calI(HalfEdge* self, FLOAT_T val) nogil: cdef inline FLOAT_T calI(HalfEdge* self, FLOAT_T val) nogil:
@ -536,7 +536,7 @@ cdef class VoronoiContainer:
cdef Site xi cdef Site xi
cdef HalfEdge em, ep cdef HalfEdge em, ep
cdef Vector2D p, q, la, da, Rla, centroid, cent_part cdef Vector2D p, q, la, da, ya, Rla, centroid, cent_part
cdef FLOAT_T [:] area = np.zeros(self.sites.shape[0], dtype=FLOAT) 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 FLOAT_T [:] perim = np.zeros(self.sites.shape[0], dtype=FLOAT)
@ -556,6 +556,7 @@ cdef class VoronoiContainer:
la_mag = la.mag(&la) la_mag = la.mag(&la)
area_p = la.dot(&la, da.rot(&da)) area_p = la.dot(&la, da.rot(&da))
Rla = la.rot(&la) Rla = la.rot(&la)
ya = Rla.copy.smul(&Rla, -2*area_p/la.dot(&la, la))
# Calculating centroid. # Calculating centroid.
cent_part = p.copy.vadd(&p, q) cent_part = p.copy.vadd(&p, q)
@ -567,8 +568,8 @@ cdef class VoronoiContainer:
em.cache.la_mag(&em, la_mag) em.cache.la_mag(&em, la_mag)
em.cache.da(&em, da) em.cache.da(&em, da)
em.cache.da_mag(&em, da.mag(&da)) em.cache.da_mag(&em, da.mag(&da))
em.cache.area_p(&em, area_p) em.cache.ya(&em, ya)
em.cache.xij(&em, Rla.copy.smul(&Rla, -area_p/la.dot(&la, la))) em.cache.ya_mag(&em, ya.mag(&ya))
area[i] += area_p area[i] += area_p
perim[i] += la_mag perim[i] += la_mag
@ -584,29 +585,20 @@ cdef class VoronoiContainer:
@staticmethod @staticmethod
cdef inline Matrix2x2 calc_H(HalfEdge em, HalfEdge ep) nogil: cdef inline Matrix2x2 calc_H(HalfEdge em, HalfEdge ep) nogil:
cdef Vector2D xmv, xpv, im, mp, right, Rpm, Rim, f cdef Vector2D xj, xk, Rxjk, v
cdef Matrix2x2 h cdef Matrix2x2 top
cdef FLOAT_T im2, mp2
# Vectors from xi to xm and xp. xj, xk = em.cache.ya(&em, NAN_VECTOR), ep.cache.ya(&ep, NAN_VECTOR)
xmv, xpv = em.cache.xij(&em, NAN_VECTOR), ep.cache.xij(&ep, NAN_VECTOR) Rxjk = xk.copy.vsub(&xk, xj)
im, mp = xmv.copy.neg(&xmv), xmv.copy.vsub(&xmv, xpv) # -xmv, xmv - xpv Rxjk.self.smul(&Rxjk, 2)
im2, mp2 = -(xmv.dot(&xmv, xmv)), xmv.dot(&xmv, xmv) - xpv.dot(&xpv, xpv) Rxjk = Rxjk.rot(&Rxjk)
# (-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. v = ep.cache.da(&ep, NAN_VECTOR)
f = h.vecmul(&h, right) # [Rpm | Rim]*right top = R.copy.smul(&R, xj.dot(&xj, xj) - xk.dot(&xk, xk))
h = R.copy.smul(&R, mp2*(2*mp.dot(&mp, Rim))) # fp*g, g is a scalar. top.self.msub(&top, _Matrix2x2(v.x*Rxjk.x, v.x*Rxjk.y, v.y*Rxjk.x, v.y*Rxjk.y))
# (fp*g - f*gp)/(g**2). f is a column vector, gp = 2*Rpm is a row vector. top.self.sdiv(&top, -Rxjk.dot(&Rxjk, xj))
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
return _Matrix2x2(top.a, top.c, top.b, top.d)
@staticmethod @staticmethod
cdef inline bint sign(FLOAT_T [::1] ref, FLOAT_T [::1] p, FLOAT_T [::1] q): cdef inline bint sign(FLOAT_T [::1] ref, FLOAT_T [::1] p, FLOAT_T [::1] q):
@ -677,10 +669,10 @@ cdef class VoronoiContainer:
self.add_sites(step*k1) self.add_sites(step*k1)
).gradient ).gradient
return (step/2)*(k1+k2), k1 return -(step/2)*(k1+k2), -k1
def hessian(self, d: float) -> np.ndarray: def approx_hessian(self, d: float) -> np.ndarray:
""" """
Obtains the approximate Hessian. Obtains the approximate Hessian.
:param d: [float] small d for approximation. :param d: [float] small d for approximation.
@ -707,6 +699,117 @@ cdef class VoronoiContainer:
return HE return HE
def radialt_hessian(self) -> np.ndarray:
HE = np.zeros((2*self.n, 2*self.n))
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, xk
cdef HalfEdge e, em, ep, fj, fk
cdef Vector2D temp1, temp2, dau, dapu
cdef Matrix2x2 sigI, sigJ, sigK, p, q, tempm, toI, toJ, toK
cdef INT_T i, j, k, z
for i in range(self.n):
xi = _Site(i, &info)
ep = xi.edge(&xi)
e = ep.prev(&ep)
j = 0
while j < xi.edge_num(&xi):
e.cache.H(&e, VoronoiContainer.calc_H(e, ep))
e = e.next(&e)
j = j + 1
for i in range(self.n):
xi = _Site(i, &info)
e = xi.edge(&xi)
z = 0
while z < xi.edge_num(&xi):
em, ep = e.prev(&e), e.next(&e)
fj, fk = em.twin(&em), e.twin(&e)
fj, fk = fj.next(&fj), fk.next(&fk)
xj, xk = fj.face(&fj), fk.face(&fk)
j, k = xj.index(&xj) % self.n, xk.index(&xk) % self.n
if k < 0:
k = k + self.n
if j < 0:
j = j + self.n
sigI = e.cache.H(&e, NAN_MATRIX)
sigJ = fj.cache.H(&fj, NAN_MATRIX)
sigK = fk.cache.H(&fk, NAN_MATRIX)
### Calculating of p
temp1 = e.cache.la(&e, NAN_VECTOR)
temp1 = temp1.rot(&temp1)
temp1.self.sdiv(
&temp1,
e.cache.la_mag(&e, NAN) * e.cache.ya_mag(&e, NAN) / 2
)
dau = e.cache.da(&e, NAN_VECTOR)
dau.self.sdiv(&dau, e.cache.da_mag(&e, NAN))
dapu = ep.cache.da(&ep, NAN_VECTOR)
dapu.self.sdiv(&dapu, ep.cache.da_mag(&ep, NAN))
temp2 = dapu.copy.vsub(&dapu, dau)
temp2 = temp2.rot(&temp2)
p = _Matrix2x2(
temp1.x*temp2.x, temp1.x*temp2.y,
temp1.y*temp2.x, temp1.y*temp2.y
)
### Calculating of q
temp2 = e.cache.la(&e, NAN_VECTOR)
temp1 = temp2.rot(&temp2)
q = _Matrix2x2(
temp1.x*temp2.x, temp1.x*temp2.y,
temp1.y*temp2.x, temp1.y*temp2.y
)
q.self.sdiv(&q, e.cache.la_mag(&e, NAN)**2)
q.self.msub(&q, R)
q.self.smul(
&q,
e.cache.calI(&e, NAN) / e.cache.la_mag(&e, NAN)
)
temp2 = em.cache.la(&em, NAN_VECTOR)
temp1 = temp2.rot(&temp2)
tempm = _Matrix2x2(
temp1.x*temp2.x, temp1.x*temp2.y,
temp1.y*temp2.x, temp1.y*temp2.y
)
tempm.self.sdiv(&tempm, em.cache.la_mag(&em, NAN)**2)
tempm = R.copy.msub(&R, tempm)
tempm.self.smul(
&tempm,
em.cache.calI(&em, NAN) / em.cache.la_mag(&em, NAN)
)
q.self.madd(&q, tempm)
# Calculating components that go to the respective sites
toI = sigI.copy.mmul(&sigI, p.copy.madd(&p, q))
toI.self.msub(&toI, p)
toJ = sigJ.copy.mmul(&sigJ, p.copy.madd(&p, q))
toK = sigK.copy.mmul(&sigK, p.copy.madd(&p, q))
HE[2*i: 2*(i+1), 2*i: 2*(i+1)] += np.array([[toI.a, toI.b], [toI.c, toI.d]])
HE[2*i: 2*(i+1), 2*j: 2*(j+1)] += np.array([[toJ.a, toJ.b], [toJ.c, toJ.d]])
HE[2*i: 2*(i+1), 2*k: 2*(k+1)] += np.array([[toK.a, toK.b], [toK.c, toK.d]])
e = e.next(&e)
z = z + 1
return -2*self.r*HE
def site_vert_arr(self): # -> List[np.ndarray] def site_vert_arr(self): # -> List[np.ndarray]
cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points, cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,