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

Script updates and minor setting changes

Browse Source
This commit is contained in:
Kenneth Jao 2022-01-29 20:00:53 -05:00
parent 583615ba4b
commit b32e0b1a8c
17 changed files with 515 additions and 319 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

101
scripts/defect_energy.py
View File

@ -1,101 +0,0 @@
import numpy as np, os
import matplotlib.pyplot as plt
from squish import Simulation, ordered
from squish.common import OUTPUT_DIR
from script_tools import RC_SETTINGS, get_args, get_data, format_data
NAME = "DefectEnergy"
def main():
sims_path, regen = get_args(
"Generates graph for Average Defects and Energy per Defect.",
"folder that contains simulations at various N",
)
def f():
data = {}
files = list(sims_path.iterdir())
for i, file in enumerate(files):
sim, frames = Simulation.load(file)
domain = sim.domain
e_hex = ordered.e_hex(domain)
defects, energy = [], []
for frame in frames:
if np.var(frame["stats"]["avg_radius"]) > 1e-8:
defects.append(
np.count_nonzero(frame["stats"]["site_edge_count"] != 6)
)
energy.append(100 * frame["energy"] / domain.n)
avg_defects = sum(defects) / (1 if len(defects) == 0 else len(defects))
m, b = np.polyfit(defects, energy, 1)
data[sim.domain.n] = [avg_defects, m]
hashes = int(21 * i / len(files))
print(
f'Processed N={sim.domain.n:03} |{"#"*hashes}{" "*(20-hashes)}|'
+ f" {i+1}/{len(files)} simulations processed.",
flush=True,
end="\r",
)
print(flush=True)
return format_data(
data, key_name="N", col_names=["Average Defects", "Energy Per Defect"]
)
plt.rcParams.update(RC_SETTINGS)
data = get_data(sims_path / (NAME + ".pkl"), f, regen=regen)
ns, defects, epds = data["N"], data["Average Defects"], data["Energy Per Defect"]
epds *= 100
fig = plt.figure(figsize=(18, 15))
gs = fig.add_gridspec(1, 1)
ax = fig.add_subplot(gs[0])
ax2 = ax.twinx()
ax3 = ax.twinx()
m0, b0 = np.polyfit(ns, defects, 1)
m1, b1 = np.polyfit(ns, defects * epds, 1)
ax.plot(ns, defects, color="C0", alpha=0.5)
ax.plot(ns, m0 * ns + b0, color="C0", linestyle="dashed")
ax2.plot(ns, epds, color="C1", alpha=0.5)
ax3.plot(ns, defects * epds / 10, color="C2", alpha=0.5)
ax3.plot(ns, (m1 * ns + b1) / 10, color="C2", linestyle="dashed")
ax.set_ylim(3, 37)
ax.set_yticks(np.arange(5, 40, 5))
ax2.set_ylim(-3 * 0.4, 18 + 3 * 0.4)
ax2.set_yticks(np.arange(0, 20, 3))
ax3.set_ylim(-3 * 0.5, 18 + 3 * 0.4)
ax3.set_yticks([])
ax3.spines["right"].set_visible(False)
ax3.spines.right.set_position(("axes", 1.11))
ax.grid(zorder=0)
ax.set_xlabel("N")
ax.set_ylabel("Average Defects", color="C0")
ax2.set_ylabel(r"Energy per Defect $\left[\times 10^{4} \right]$", color="C1")
ax3.set_ylabel(r"Defect Energy $\left[\times 10^{3} \right]$", color="C2")
fig.savefig(OUTPUT_DIR / (NAME + ".png"))
print(f"Wrote to {OUTPUT_DIR / (NAME + '.png')}")
if __name__ == "__main__":
os.environ["QT_log10GING_RULES"] = "*=false"
try:
main()
except KeyboardInterrupt:
print("Program terminated by user.")

68
scripts/defect_scatter.py Normal file
View File

@ -0,0 +1,68 @@
import argparse, numpy as np, os, math
import matplotlib.pyplot as plt
from multiprocessing import Pool, cpu_count
from squish import Simulation, ordered
from squish.common import OUTPUT_DIR
from script_tools import RC_SETTINGS
NAME = "DefectScatter"
def main():
parser = argparse.ArgumentParser(description="Graphs scatter of VEE and defects.")
parser.add_argument(
"sims_path", metavar="sim_dir", help="simulation to make scatter plot of"
)
args = parser.parse_args()
sim, frames = Simulation.load(args.sims_path)
defects, energy = [], []
e_hex = ordered.e_hex(sim.domain)
for frame in frames:
defects.append(np.count_nonzero(frame["stats"]["site_edge_count"] != 6))
energy.append(frame["energy"] / sim.domain.n - e_hex)
defects = np.array(defects)
energy = 100 * np.array(energy)
plt.rcParams.update(RC_SETTINGS)
fig = plt.figure(figsize=(18, 15))
gs = fig.add_gridspec(1, 1)
ax = fig.add_subplot(gs[0])
m, b = np.polyfit(defects, energy, 1)
ax.scatter(defects, energy, color="C0", marker="*", s=100)
ax.plot(np.arange(0, 64), np.arange(0, 64) * m + b, zorder=3, color="C1")
ax.scatter(0, b, color="C2", s=500, marker="*", zorder=50)
ax.annotate(
r"$\zeta_0$",
xy=(0, b),
xytext=(10, -50),
textcoords="offset points",
# arrowprops={"arrowstyle": "->", "color": "black"},
)
ax.set_xlim(0, 64)
ax.set_xticks(np.arange(0, 65, 8))
ax.set_ylim(0, 3.6)
ax.set_yticks(np.arange(0, 3.6, 0.5))
ax.set_xlabel("Number of Defects")
ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$")
ax.grid(zorder=0)
fig.savefig(OUTPUT_DIR / (NAME + ".png"))
print(f"Wrote to {OUTPUT_DIR / (NAME + '.png')}")
if __name__ == "__main__":
os.environ["QT_log10GING_RULES"] = "*=false"
try:
main()
except KeyboardInterrupt:
print("Program terminated by user.")

75
scripts/density_of_states.py
View File

@ -7,80 +7,47 @@ from pathlib import Path
from squish import Simulation, DomainParams
from squish.common import OUTPUT_DIR
from script_tools import RC_SETTINGS, get_args, get_data, get_simulation_data
def main():
# Loading arguments.
parser = argparse.ArgumentParser("Outputs width search data into diagrams")
parser.add_argument(
"sims_path",
metavar="path/to/data",
help="Path to simulation folders with generated data.pkl from aspect_diagrams.py",
sims_path, regen = get_args(
"Density of states of various N across alpha",
"folders that contains various N simulations to plot",
)
packages = []
args = parser.parse_args()
for fol in Path(args.sims_path).iterdir():
for fol in sims_path.iterdir():
if fol.is_file():
continue
store = fol / "data.pkl"
if store.is_file():
with open(store, "rb") as f:
packages.append(pickle.load(f))
else:
print(
f"{store} not found! Use aspect_diagrams.py to generate this file first."
)
data, n, r = get_data(
fol / "package.pkl", get_simulation_data, args=(fol,), regen=regen
)
domain = DomainParams(n, 1, 1, r)
if len(packages) == []:
print("No data.pkl files found, terminating")
return
packages.append([data, domain])
plt.rcParams.update(
{
"axes.titlesize": 45,
"axes.labelsize": 45,
"xtick.labelsize": 40,
"ytick.labelsize": 40,
"xtick.major.width": 2,
"ytick.major.width": 2,
"xtick.major.size": 5,
"ytick.major.size": 5,
"xtick.minor.width": 1,
"ytick.minor.width": 1,
"xtick.minor.size": 3,
"ytick.minor.size": 3,
"legend.fontsize": 40,
"lines.linewidth": 3.5,
"font.family": "cm",
"font.size": 40,
"text.usetex": True,
"text.latex.preamble": r"\usepackage{amsmath}",
"figure.constrained_layout.use": True,
}
)
packages.sort(key=lambda x: x[1].n)
plt.rcParams.update(RC_SETTINGS)
fig = plt.figure(figsize=(18, 15))
gs = fig.add_gridspec(1, 1)
ax = fig.add_subplot(gs[0])
packages.sort(key=lambda x: x[0].n)
for j, package in enumerate(packages):
domain, data, order_data, widths = package
data, domain = package
distinct_unordered = []
for width in widths:
equal_shape = list([c[1] for c in data["distinct"][width]])
distinct_unordered.append(equal_shape.count(False))
distinct_disordered = []
for ordered in data["distinct"]["Ordered"]:
distinct_disordered.append(np.count_nonzero(ordered == False))
ax.plot(widths, distinct_unordered, label=f"N={domain.n}")
ax.plot(data["all"]["alpha"], distinct_disordered, label=f"N={domain.n}")
widths = packages[0][3]
ax.set_xticks([round(w, 2) for w in widths[::10]])
ax.set_xticklabels([f"{round(w, 3):.2f}" for w in widths[::10]], rotation=90)
alphas = packages[0][0]["all"]["alpha"]
ax.set_xticks([round(w, 2) for w in alphas[::10]])
ax.set_xticklabels([f"{round(w, 3):.2f}" for w in alphas[::10]], rotation=90)
ax.set_xlim(0.3, 1.0)
ax.set_xlabel("Aspect Ratio")
ax.set_ylabel("Number of Distinct Equilibria")

173
scripts/frustration.py Normal file
View File

@ -0,0 +1,173 @@
import numpy as np, os, math
import matplotlib.pyplot as plt
from multiprocessing import Pool, cpu_count
from squish import Simulation, ordered
from squish.common import OUTPUT_DIR
from script_tools import RC_SETTINGS, get_args, get_data, format_data
NAME = "DefectEnergy"
NAME2 = "Frustration"
def proc(path):
sim, frames = Simulation.load(path)
domain = sim.domain
e_hex = ordered.e_hex(domain)
defects, energy = [], []
for frame in frames:
# if np.var(frame["stats"]["avg_radius"]) > 1e-8:
defects.append(np.count_nonzero(frame["stats"]["site_edge_count"] != 6))
energy.append(frame["energy"] / domain.n - e_hex)
configs = ordered.configurations(domain)
order_energies = []
for config in configs:
rbar = ordered.avg_radius(domain, config)
order_energies.append(
(
2 * domain.w * domain.h
+ 2 * math.pi * domain.n * (domain.r ** 2 - 2 * domain.r * rbar)
)
/ domain.n
- e_hex
)
avg_defects = sum(defects) / len(defects)
avg_vee = sum(energy) / len(energy)
m, b = np.polyfit(defects, energy, 1)
return [domain.n, avg_defects, avg_vee, m, b, min(order_energies)]
def main():
sims_path, regen = get_args(
"Generates graph for Average Defects and Energy per Defect for alpha=1",
"folder that contains simulations at various N",
)
def f():
data = {}
files = list(sims_path.iterdir())
with Pool(cpu_count()) as pool:
for i, res in enumerate(pool.imap_unordered(proc, files)):
data[res[0]] = res[1:]
hashes = int(21 * i / len(files))
print(
f'Processed N={res[0]} |{"#"*hashes}{" "*(20-hashes)}|'
+ f" {i+1}/{len(files)} simulations processed.",
flush=True,
end="\r",
)
print(flush=True)
return format_data(
data,
key_name="N",
col_names=[
"Average Defects",
"Average VEE",
"Energy Per Defect",
"Ground VEE",
"Minimum Ordered VEE",
],
)
plt.rcParams.update(RC_SETTINGS)
data = get_data(sims_path / (NAME + ".pkl"), f, regen=regen)
ns, avg_defects, avg_vees, epds, vee0s, min_order = (
data["N"],
data["Average Defects"],
100 * data["Average VEE"],
100 * data["Energy Per Defect"],
100 * data["Ground VEE"],
100 * data["Minimum Ordered VEE"],
)
fig = plt.figure(figsize=(18, 15))
gs = fig.add_gridspec(1, 1)
ax = fig.add_subplot(gs[0])
# ax2 = ax.twinx()
# ax3 = ax.twinx()
m0, b0 = np.polyfit(ns, avg_defects, 1)
# m1, b1 = np.polyfit(ns, avg_defects * epds, 1)
lns1 = ax.plot(
ns, avg_defects, color="C0", alpha=0.5, label=r"$\langle \mathrm{D} \rangle$"
)
ax.plot(ns, m0 * ns + b0, color="C0", linestyle="dashed")
lns2 = ax.plot(
ns, 100 * epds, color="C1", alpha=0.5, label=r"$\zeta \times 10^{4}$"
)
lns3 = ax.plot(
ns,
100 * avg_defects * epds / 10,
color="C2",
alpha=0.5,
label=r"$\mathcal{F} \times 10^{3}$",
)
# ax3.plot(ns, 100*(m1 * ns + b1) / 10, color="C2", linestyle="dashed")
ax.set_ylim(3, 37)
ax.set_yticks(np.arange(5, 40, 5))
# ax2.set_ylim(-3 * 0.4, 18 + 3 * 0.4)
# ax2.set_yticks(np.arange(0, 20, 3))
# ax3.set_ylim(-3 * 0.5, 18 + 3 * 0.4)
# ax3.set_yticks([])
# ax3.spines["right"].set_visible(False)
# ax3.spines.right.set_position(("axes", 1.11))
# lns = lns1 + lns2 + lns3
# labs = [l.get_label() for l in lns]
ax.grid(zorder=0)
# ax.legend(lns, labs, loc="lower right")
ax.legend()
ax.set_xlabel("N")
# ax.set_ylabel(r"$\langle \mathrm{D} \rangle$")
# ax2.set_ylabel("VEE")
# ax2.set_ylabel(r"$\zeta \left[\times 10^{4} \right]$", color="C1")
# ax3.set_ylabel(r"Defect Energy $\left[\times 10^{3} \right]$", color="C2")
fig.savefig(OUTPUT_DIR / (NAME + ".png"))
print(f"Wrote to {OUTPUT_DIR / (NAME + '.png')}")
return
fig = plt.figure(figsize=(30, 15))
gs = fig.add_gridspec(1, 1)
ax = fig.add_subplot(gs[0])
ax.plot(ns, vee0s, label="$\mathrm{VEE}_0$")
ax.plot(ns, avg_defects * epds, label=r"$\mathcal{F}$")
ax.plot(ns, avg_vees, label=r"$\langle \mathrm{VEE} \rangle$")
ax.plot(ns, min_order, label="Minimum Ordered", color="C5", alpha=0.5)
print(np.linalg.norm(avg_vees - (vee0s + avg_defects * epds)))
ax.grid(zorder=0)
ax.legend()
ax.set_ylim(0, 3)
ax.set_yticks(np.arange(0, 3.1, 0.5))
ax.set_xlabel("N")
ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$")
fig.savefig(OUTPUT_DIR / (NAME2 + ".png"))
print(f"Wrote to {OUTPUT_DIR / (NAME2 + '.png')}")
if __name__ == "__main__":
os.environ["QT_log10GING_RULES"] = "*=false"
try:
main()
except KeyboardInterrupt:
print("Program terminated by user.")

2
scripts/hexagon_alpha.py
View File

@ -68,7 +68,7 @@ def main():
ax.scatter(alphas, ns, s=5)
ax.set_xlim(-0.01, 1.1)
ax.set_xlim(-0.01, 1.01)
ax.set_xticks(np.arange(0, 1.01, 0.1))
ax.set_ylim(0, args.max_n)

70
scripts/ordered_rbar.py Normal file
View File

@ -0,0 +1,70 @@
import argparse, numpy as np, os
from multiprocessing import Pool, cpu_count
import matplotlib.pyplot as plt
from squish import Simulation, ordered
from squish.common import OUTPUT_DIR
from script_tools import RC_SETTINGS, get_data, format_data
NAME = "OrdereredRBar"
def main():
parser = argparse.ArgumentParser(
description="Ordered rbar vector for specified frames of a simulation"
)
parser.add_argument(
"sims_path",
metavar="sim_dir",
help="folder that contains of perturbations from an equilibrium",
)
parser.add_argument(
"frames", metavar="frames", type=int, help="frame numbers to select", nargs="+"
)
args = parser.parse_args()
rbars = [None] * len(args.frames)
sim, frames = Simulation.load(args.sims_path)
for i, frame in enumerate(frames):
if i not in args.frames:
continue
rbars[args.frames.index(i)] = sorted(frame["stats"]["avg_radius"])
plt.rcParams.update(RC_SETTINGS)
fig = plt.figure(figsize=(20, 10))
gs = fig.add_gridspec(1, 1)
ax = fig.add_subplot(gs[0])
alphabet = "abcdefhijklmnopqrstuvwxyz"
for i, rbar in enumerate(rbars):
ax.plot(list(range(len(rbar))), rbar, label=alphabet[i])
ax.set_xlim(-0.1, len(rbars[0]))
ax.set_xticks(np.arange(0, 401, 40))
# ax.set_ylim(0, args.max_n)
ax.plot(
list(range(len(rbars[0]))),
[ordered.R_HEX] * (len(rbars[0])),
linestyle="dotted",
color="black",
)
ax.set_ylabel(r"$\bar{r}$")
ax.grid(zorder=0)
ax.legend()
fig.savefig(OUTPUT_DIR / (NAME + ".png"))
print(f"Wrote to {OUTPUT_DIR / (NAME + '.png')}")
if __name__ == "__main__":
os.environ["QT_log10GING_RULES"] = "*=false"
try:
main()
except KeyboardInterrupt:
print("Program terminated by user.")

74
scripts/ordered_scatter.py Normal file
View File

@ -0,0 +1,74 @@
import argparse, numpy as np, os, math
import matplotlib.pyplot as plt
from squish import Simulation, ordered, DomainParams
from squish.common import OUTPUT_DIR
from script_tools import RC_SETTINGS, get_args, get_data, get_ordered_data
NAME = "OrderedScatter"
def main():
parser = argparse.ArgumentParser(
description="Graphs scatter of ordered energy at fixed N."
)
parser.add_argument("n", metavar="N", type=int, help="N to make scatter plot of")
parser.add_argument("r", metavar="R", type=float, help="natural radius of object")
parser.add_argument(
"--regenerate",
dest="regen",
action="store_true",
help="regenerates the cache file for processed data",
)
args = parser.parse_args()
ordered_data = get_data(
OUTPUT_DIR / "OrderedCache" / f"{args.n}.pkl",
get_ordered_data,
args=(DomainParams(args.n, 1, 1, args.r), np.linspace(0.3, 1, 141)),
regen=args.regen,
)
plt.rcParams.update(RC_SETTINGS)
fig = plt.figure(figsize=(15, 15))
gs = fig.add_gridspec(1, 1)
ax = fig.add_subplot(gs[0])
e_hex = ordered.e_hex(DomainParams(args.n, 1, 1, args.r))
for alpha, energies in zip(ordered_data["alpha"], ordered_data["Energy"]):
ax.scatter(
[alpha] * len(energies), 100 * (energies / args.n - e_hex), color="C0"
)
ax.set_xlim(0.3, 1)
ax.set_xticks(np.arange(0.3, 1.01, 0.1))
ax.set_ylim(-0.1, 10)
# ax.set_ylim(0, 3.6)
# ax.set_yticks(np.arange(0, 3.6, 0.5))
props = dict(boxstyle="round", facecolor="white", alpha=0.8, zorder=20)
ax.text(
0.873,
0.97,
f"N={args.n}",
transform=ax.transAxes,
verticalalignment="top",
bbox=props,
)
ax.set_xlabel("Aspect Ratio")
ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$")
ax.grid(zorder=0)
fig.savefig(OUTPUT_DIR / (NAME + ".png"))
print(f"Wrote to {OUTPUT_DIR / (NAME + '.png')}")
if __name__ == "__main__":
os.environ["QT_log10GING_RULES"] = "*=false"
try:
main()
except KeyboardInterrupt:
print("Program terminated by user.")

3
scripts/probability_of_disorder.py
View File

@ -66,7 +66,8 @@ def main():
ax2.plot(alphas, all_disorder_count, color="C4")
ax.set_xlim(0.3, 1)
ax.set_xticks(np.arange(0.3, 1.01, 0.1))
ax.set_xticks([round(w, 2) for w in alphas[::10]])
ax.set_xticklabels([f"{round(w, 3):.2f}" for w in alphas[::10]], rotation=90)
start, end = ax.get_ylim()
# space = np.linspace(0, end, 20)

19
scripts/script_tools.py
View File

@ -1,6 +1,6 @@
from __future__ import annotations
from typing import List
import argparse, pickle, numpy as np, math
import argparse, pickle, numpy as np, math, os
from pathlib import Path
from multiprocessing import Pool, cpu_count
import matplotlib.pyplot as plt
@ -52,7 +52,13 @@ def get_args(
def get_data(
path: Path, func: Callable[Any, Any], args: Tuple[Any] = (), regen: bool = False
) -> Any:
if path.is_file() and not regen:
if regen:
try:
os.remove(path)
except FileNotFoundError:
pass
if path.is_file():
with open(path, "rb") as f:
return pickle.load(f)
else:
@ -69,7 +75,8 @@ def format_data(
new_data = {}
new_data[key_name] = np.array([x[0] for x in data])
for i, col_name in enumerate(col_names):
if type(data[0][1][i]) is list:
col_value_type = type(data[0][1][i])
if col_value_type is list or col_value_type is tuple:
new_data[col_name] = [np.array(x[1][i]) for x in data]
else:
new_data[col_name] = np.array([x[1][i] for x in data])
@ -145,7 +152,8 @@ def ordered_data_proc(
else:
coercivities.append(eigs[0])
return (alpha, sorted(energies), sorted(coercivities)[::-1])
energies, coercivities = list(zip(*sorted(zip(energies, coercivities))))
return (alpha, energies, coercivities)
def get_simulation_data(filepath: Path) -> Tuple[Dict, numpy.ndarray, DomainParams]:
@ -189,6 +197,8 @@ def simulation_data_proc(file: Path) -> Tuple[float, List[float], List[float]]:
alls[1].append(np.var(frame_info["stats"]["avg_radius"]) <= 1e-8)
alls[2].append(np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6))
alls = list(zip(*sorted(zip(*alls))))
sim, frames = Simulation.load(file)
sim.frames = list(frames)
counts = sim.get_distinct()
@ -202,5 +212,6 @@ def simulation_data_proc(file: Path) -> Tuple[float, List[float], List[float]]:
)
distincts.append(counts)
distincts = list(zip(*sorted(zip(*distincts))))
return sim.domain.w / sim.domain.h, alls, distincts

14
scripts/vee.py
View File

@ -25,10 +25,6 @@ def main():
sims_path / "package.pkl", get_simulation_data, args=(sims_path,), regen=regen
)
domain, alphas = DomainParams(n, 1, 1, r), data["all"]["alpha"]
if regen:
os.remove(OUTPUT_DIR / "OrderedCache" / f"{n}".pkl)
ordered_data = get_data(
OUTPUT_DIR / "OrderedCache" / f"{n}.pkl",
get_ordered_data,
@ -45,10 +41,10 @@ def main():
min_disorder.append(min(disorder_energies))
max_disorder.append(max(disorder_energies))
min_order = []
min_order, min_order_coer = [], []
for i, energies in enumerate(ordered_data["Energy"]):
min_order.append(energies[0])
min_order_coer.append(ordered_data["Coercivity"][i][0])
e_hex = ordered.e_hex(domain)
min_disorder = np.array(min_disorder) / domain.n - e_hex
max_disorder = np.array(max_disorder) / domain.n - e_hex
@ -72,7 +68,13 @@ def main():
label="Max Disordered",
)
ax.scatter(
hex_ratios, [0] * len(hex_ratios), color="C2", s=120, marker="H", zorder=50
)
ax.set_xlim(0.3, 1)
ax.set_xticks([round(w, 2) for w in alphas[::10]])
ax.set_xticklabels([f"{round(w, 3):.2f}" for w in alphas[::10]], rotation=90)
# start, end = ax.get_ylim()
# space = np.linspace(0, end, 20)

2
scripts/vee_diff_pod.py
View File

@ -78,7 +78,7 @@ def main():
100 * (min_order - min_disorder), all_disorder_count, label=f"N={domain.n}"
)
ax.set_ylabel("Probability of Disordered Equilibria")
ax.set_ylabel("Perecent of Disordered Equilibria")
ax.set_xlabel(r"VEE Difference $\left[\times 10^{2}\right]$")
ax.yaxis.set_major_formatter(mtick.PercentFormatter())
ax.set_ylim(48, 102)

11
squish/diagram.py
View File

@ -13,6 +13,7 @@ from scipy.spatial import Voronoi, voronoi_plot_2d
from multiprocessing import Pool, cpu_count
from .common import DomainParams, OUTPUT_DIR
from squish import ordered
SYMM = np.array(
[[0, 0], [1, 0], [1, 1], [0, 1], [-1, 1], [-1, 0], [-1, -1], [0, -1], [1, -1]]
@ -187,11 +188,7 @@ class Diagram:
scale = 1.5
area = n <= 50
offset = (
2
- 2 * domain.r * (6 * 3 ** (-0.25) * math.sqrt(2) * math.atanh(0.5))
+ 2 * math.pi * domain.r ** 2
)
e_hex = ordered.e_hex(domain)
# Make color map axis.
divider = make_axes_locatable(ax)
@ -218,7 +215,7 @@ class Diagram:
polygon = [self.sim.voronois[i].vertices[k] for k in region]
ax.fill(
*zip(*polygon),
color=mapper.to_rgba(energies[j % n] - offset),
color=mapper.to_rgba(energies[j % n] - e_hex),
zorder=0,
)
@ -278,7 +275,7 @@ class Diagram:
ax.text(
0.065,
0.935,
f"VEE: {round(sum(energies)/n - offset, 8)}",
f"VEE: {round(sum(energies)/n - e_hex, 8)}",
transform=ax.transAxes,
verticalalignment="top",
bbox=props,

201
squish/energy.c
View File

@ -6969,7 +6969,7 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
PyObject *__pyx_t_5 = NULL;
__Pyx_memviewslice __pyx_t_6 = { 0, 0, { 0 }, { 0 }, { 0 } };
Py_ssize_t __pyx_t_7;
__pyx_t_6squish_4core_INT_T __pyx_t_8;
Py_ssize_t __pyx_t_8;
__pyx_t_6squish_4core_INT_T __pyx_t_9;
int __pyx_t_10;
__pyx_t_6squish_7voronoi_HalfEdge __pyx_t_11;
@ -6977,7 +6977,8 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
__pyx_t_6squish_4core_INT_T __pyx_t_13;
__pyx_t_6squish_4core_INT_T __pyx_t_14;
__pyx_t_6squish_4core_INT_T __pyx_t_15;
__Pyx_memviewslice __pyx_t_16 = { 0, 0, { 0 }, { 0 }, { 0 } };
__pyx_t_6squish_4core_INT_T __pyx_t_16;
__Pyx_memviewslice __pyx_t_17 = { 0, 0, { 0 }, { 0 }, { 0 } };
int __pyx_lineno = 0;
const char *__pyx_filename = NULL;
int __pyx_clineno = 0;
@ -7081,157 +7082,93 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
/* "squish/energy.pyx":288
*
* cdef INT_T i, j, k
* for i in prange(self.sites.shape[0], nogil=True): # <<<<<<<<<<<<<<
* for i in range(self.sites.shape[0]): # <<<<<<<<<<<<<<
* xi = _Site(i, &info)
* e = xi.edge(&xi)
*/
{
#ifdef WITH_THREAD
PyThreadState *_save;
Py_UNBLOCK_THREADS
__Pyx_FastGIL_Remember();
#endif
/*try:*/ {
if (unlikely(!__pyx_v_self->__pyx_base.sites.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 288, __pyx_L4_error)}
__pyx_t_7 = (__pyx_v_self->__pyx_base.sites.shape[0]);
if ((1 == 0)) abort();
{
#if ((defined(__APPLE__) || defined(__OSX__)) && (defined(__GNUC__) && (__GNUC__ > 2 || (__GNUC__ == 2 && (__GNUC_MINOR__ > 95)))))
#undef likely
#undef unlikely
#define likely(x) (x)
#define unlikely(x) (x)
#endif
__pyx_t_9 = (__pyx_t_7 - 0 + 1 - 1/abs(1)) / 1;
if (__pyx_t_9 > 0)
{
#ifdef _OPENMP
#pragma omp parallel private(__pyx_t_10, __pyx_t_11, __pyx_t_12)
#endif /* _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_xi)
#endif /* _OPENMP */
for (__pyx_t_8 = 0; __pyx_t_8 < __pyx_t_9; __pyx_t_8++){
{
__pyx_v_i = (__pyx_t_6squish_4core_INT_T)(0 + 1 * __pyx_t_8);
/* Initialize private variables to invalid values */
__pyx_v_j = ((__pyx_t_6squish_4core_INT_T)0xbad0bad0);
if (unlikely(!__pyx_v_self->__pyx_base.sites.memview)) {PyErr_SetString(PyExc_AttributeError,"Memoryview is not initialized");__PYX_ERR(0, 288, __pyx_L1_error)}
__pyx_t_7 = (__pyx_v_self->__pyx_base.sites.shape[0]);
__pyx_t_8 = __pyx_t_7;
for (__pyx_t_9 = 0; __pyx_t_9 < __pyx_t_8; __pyx_t_9+=1) {
__pyx_v_i = __pyx_t_9;
/* "squish/energy.pyx":289
/* "squish/energy.pyx":289
* cdef INT_T i, j, k
* for i in prange(self.sites.shape[0], nogil=True):
* for i in range(self.sites.shape[0]):
* xi = _Site(i, &info) # <<<<<<<<<<<<<<
* e = xi.edge(&xi)
*
*/
__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":290
* for i in prange(self.sites.shape[0], nogil=True):
/* "squish/energy.pyx":290
* for i in range(self.sites.shape[0]):
* xi = _Site(i, &info)
* e = xi.edge(&xi) # <<<<<<<<<<<<<<
*
* 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":292
/* "squish/energy.pyx":292
* e = xi.edge(&xi)
*
* j = 0 # <<<<<<<<<<<<<<
* while j < xi.edge_num(&xi):
* em, ep = e.prev(&e), e.next(&e)
*/
__pyx_v_j = 0;
__pyx_v_j = 0;
/* "squish/energy.pyx":293
/* "squish/energy.pyx":293
*
* 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))
*/
while (1) {
__pyx_t_10 = ((__pyx_v_j < __pyx_v_xi.edge_num((&__pyx_v_xi))) != 0);
if (!__pyx_t_10) break;
while (1) {
__pyx_t_10 = ((__pyx_v_j < __pyx_v_xi.edge_num((&__pyx_v_xi))) != 0);
if (!__pyx_t_10) break;
/* "squish/energy.pyx":294
/* "squish/energy.pyx":294
* 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))
*
*/
__pyx_t_11 = __pyx_v_e.prev((&__pyx_v_e));
__pyx_t_12 = __pyx_v_e.next((&__pyx_v_e));
__pyx_v_em = __pyx_t_11;
__pyx_v_ep = __pyx_t_12;
__pyx_t_11 = __pyx_v_e.prev((&__pyx_v_e));
__pyx_t_12 = __pyx_v_e.next((&__pyx_v_e));
__pyx_v_em = __pyx_t_11;
__pyx_v_ep = __pyx_t_12;
/* "squish/energy.pyx":295
/* "squish/energy.pyx":295
* while j < xi.edge_num(&xi):
* em, ep = e.prev(&e), e.next(&e)
* e.cache.H(&e, VoronoiContainer.calc_H(em, e)) # <<<<<<<<<<<<<<
*
* e = e.next(&e)
*/
(void)(__pyx_v_e.cache->H((&__pyx_v_e), __pyx_vtabptr_6squish_7voronoi_VoronoiContainer->calc_H(__pyx_v_em, __pyx_v_e)));
(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":297
/* "squish/energy.pyx":297
* e.cache.H(&e, VoronoiContainer.calc_H(em, e))
*
* e = e.next(&e) # <<<<<<<<<<<<<<
* j = j + 1
*
*/
__pyx_v_e = __pyx_v_e.next((&__pyx_v_e));
__pyx_v_e = __pyx_v_e.next((&__pyx_v_e));
/* "squish/energy.pyx":298
/* "squish/energy.pyx":298
*
* e = e.next(&e)
* j = j + 1 # <<<<<<<<<<<<<<
*
*
*/
__pyx_v_j = (__pyx_v_j + 1);
}
}
}
}
}
}
#if ((defined(__APPLE__) || defined(__OSX__)) && (defined(__GNUC__) && (__GNUC__ > 2 || (__GNUC__ == 2 && (__GNUC_MINOR__ > 95)))))
#undef likely
#undef unlikely
#define likely(x) __builtin_expect(!!(x), 1)
#define unlikely(x) __builtin_expect(!!(x), 0)
#endif
}
/* "squish/energy.pyx":288
*
* cdef INT_T i, j, k
* for i in prange(self.sites.shape[0], nogil=True): # <<<<<<<<<<<<<<
* xi = _Site(i, &info)
* e = xi.edge(&xi)
*/
/*finally:*/ {
/*normal exit:*/{
#ifdef WITH_THREAD
__Pyx_FastGIL_Forget();
Py_BLOCK_THREADS
#endif
goto __pyx_L5;
}
__pyx_L4_error: {
#ifdef WITH_THREAD
__Pyx_FastGIL_Forget();
Py_BLOCK_THREADS
#endif
goto __pyx_L1_error;
}
__pyx_L5:;
}
__pyx_v_j = (__pyx_v_j + 1);
}
}
/* "squish/energy.pyx":301
@ -7242,9 +7179,9 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
* e = xi.edge(&xi)
*/
__pyx_t_9 = __pyx_v_self->__pyx_base.n;
__pyx_t_8 = __pyx_t_9;
for (__pyx_t_13 = 0; __pyx_t_13 < __pyx_t_8; __pyx_t_13+=1) {
__pyx_v_i = __pyx_t_13;
__pyx_t_13 = __pyx_t_9;
for (__pyx_t_14 = 0; __pyx_t_14 < __pyx_t_13; __pyx_t_14+=1) {
__pyx_v_i = __pyx_t_14;
/* "squish/energy.pyx":302
*
@ -7411,9 +7348,9 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
* HE[2*i, 2*i+1] -= dsite.b
* HE[2*i+1, 2*i] -= dsite.c
*/
__pyx_t_14 = (2 * __pyx_v_i);
__pyx_t_15 = (2 * __pyx_v_i);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_14 * __pyx_v_HE.strides[0]) ) + __pyx_t_15 * __pyx_v_HE.strides[1]) )) -= __pyx_v_dsite.a;
__pyx_t_16 = (2 * __pyx_v_i);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_15 * __pyx_v_HE.strides[0]) ) + __pyx_t_16 * __pyx_v_HE.strides[1]) )) -= __pyx_v_dsite.a;
/* "squish/energy.pyx":333
*
@ -7422,9 +7359,9 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
* HE[2*i+1, 2*i] -= dsite.c
* HE[2*i+1, 2*i+1] -= dsite.d
*/
__pyx_t_15 = (2 * __pyx_v_i);
__pyx_t_14 = ((2 * __pyx_v_i) + 1);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_15 * __pyx_v_HE.strides[0]) ) + __pyx_t_14 * __pyx_v_HE.strides[1]) )) -= __pyx_v_dsite.b;
__pyx_t_16 = (2 * __pyx_v_i);
__pyx_t_15 = ((2 * __pyx_v_i) + 1);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_16 * __pyx_v_HE.strides[0]) ) + __pyx_t_15 * __pyx_v_HE.strides[1]) )) -= __pyx_v_dsite.b;
/* "squish/energy.pyx":334
* HE[2*i, 2*i] -= dsite.a
@ -7433,9 +7370,9 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
* HE[2*i+1, 2*i+1] -= dsite.d
*
*/
__pyx_t_14 = ((2 * __pyx_v_i) + 1);
__pyx_t_15 = (2 * __pyx_v_i);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_14 * __pyx_v_HE.strides[0]) ) + __pyx_t_15 * __pyx_v_HE.strides[1]) )) -= __pyx_v_dsite.c;
__pyx_t_15 = ((2 * __pyx_v_i) + 1);
__pyx_t_16 = (2 * __pyx_v_i);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_15 * __pyx_v_HE.strides[0]) ) + __pyx_t_16 * __pyx_v_HE.strides[1]) )) -= __pyx_v_dsite.c;
/* "squish/energy.pyx":335
* HE[2*i, 2*i+1] -= dsite.b
@ -7444,9 +7381,9 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
*
* # Calculating of q
*/
__pyx_t_16 = ((2 * __pyx_v_i) + 1);
__pyx_t_15 = ((2 * __pyx_v_i) + 1);
__pyx_t_14 = ((2 * __pyx_v_i) + 1);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_15 * __pyx_v_HE.strides[0]) ) + __pyx_t_14 * __pyx_v_HE.strides[1]) )) -= __pyx_v_dsite.d;
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_16 * __pyx_v_HE.strides[0]) ) + __pyx_t_15 * __pyx_v_HE.strides[1]) )) -= __pyx_v_dsite.d;
/* "squish/energy.pyx":338
*
@ -7772,9 +7709,9 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
* HE[2*i, 2*k+1] += tempm.b
* HE[2*i+1, 2*k] += tempm.c
*/
__pyx_t_14 = (2 * __pyx_v_i);
__pyx_t_15 = (2 * __pyx_v_k);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_14 * __pyx_v_HE.strides[0]) ) + __pyx_t_15 * __pyx_v_HE.strides[1]) )) += __pyx_v_tempm.a;
__pyx_t_15 = (2 * __pyx_v_i);
__pyx_t_16 = (2 * __pyx_v_k);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_15 * __pyx_v_HE.strides[0]) ) + __pyx_t_16 * __pyx_v_HE.strides[1]) )) += __pyx_v_tempm.a;
/* "squish/energy.pyx":405
*
@ -7783,9 +7720,9 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
* HE[2*i+1, 2*k] += tempm.c
* HE[2*i+1, 2*k+1] += tempm.d
*/
__pyx_t_15 = (2 * __pyx_v_i);
__pyx_t_14 = ((2 * __pyx_v_k) + 1);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_15 * __pyx_v_HE.strides[0]) ) + __pyx_t_14 * __pyx_v_HE.strides[1]) )) += __pyx_v_tempm.b;
__pyx_t_16 = (2 * __pyx_v_i);
__pyx_t_15 = ((2 * __pyx_v_k) + 1);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_16 * __pyx_v_HE.strides[0]) ) + __pyx_t_15 * __pyx_v_HE.strides[1]) )) += __pyx_v_tempm.b;
/* "squish/energy.pyx":406
* HE[2*i, 2*k] += tempm.a
@ -7794,9 +7731,9 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
* HE[2*i+1, 2*k+1] += tempm.d
*
*/
__pyx_t_14 = ((2 * __pyx_v_i) + 1);
__pyx_t_15 = (2 * __pyx_v_k);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_14 * __pyx_v_HE.strides[0]) ) + __pyx_t_15 * __pyx_v_HE.strides[1]) )) += __pyx_v_tempm.c;
__pyx_t_15 = ((2 * __pyx_v_i) + 1);
__pyx_t_16 = (2 * __pyx_v_k);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_15 * __pyx_v_HE.strides[0]) ) + __pyx_t_16 * __pyx_v_HE.strides[1]) )) += __pyx_v_tempm.c;
/* "squish/energy.pyx":407
* HE[2*i, 2*k+1] += tempm.b
@ -7805,9 +7742,9 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
*
* f = f.twin(&f)
*/
__pyx_t_15 = ((2 * __pyx_v_i) + 1);
__pyx_t_14 = ((2 * __pyx_v_k) + 1);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_15 * __pyx_v_HE.strides[0]) ) + __pyx_t_14 * __pyx_v_HE.strides[1]) )) += __pyx_v_tempm.d;
__pyx_t_16 = ((2 * __pyx_v_i) + 1);
__pyx_t_15 = ((2 * __pyx_v_k) + 1);
*((__pyx_t_6squish_4core_FLOAT_T *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_HE.data + __pyx_t_16 * __pyx_v_HE.strides[0]) ) + __pyx_t_15 * __pyx_v_HE.strides[1]) )) += __pyx_v_tempm.d;
/* "squish/energy.pyx":401
* k = k + self.n
@ -7853,7 +7790,7 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
*
* e = e.next(&e)
*/
goto __pyx_L19_break;
goto __pyx_L12_break;
/* "squish/energy.pyx":411
* f = f.twin(&f)
@ -7864,7 +7801,7 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
*/
}
}
__pyx_L19_break:;
__pyx_L12_break:;
/* "squish/energy.pyx":414
* break
@ -7931,12 +7868,12 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
__Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
__pyx_t_16 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_nn___pyx_t_6squish_4core_FLOAT_T(__pyx_t_3, PyBUF_WRITABLE); if (unlikely(!__pyx_t_16.memview)) __PYX_ERR(0, 417, __pyx_L1_error)
__pyx_t_17 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_nn___pyx_t_6squish_4core_FLOAT_T(__pyx_t_3, PyBUF_WRITABLE); if (unlikely(!__pyx_t_17.memview)) __PYX_ERR(0, 417, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
__PYX_XDEC_MEMVIEW(&__pyx_v_self->__pyx_base.hess, 0);
__pyx_v_self->__pyx_base.hess = __pyx_t_16;
__pyx_t_16.memview = NULL;
__pyx_t_16.data = NULL;
__pyx_v_self->__pyx_base.hess = __pyx_t_17;
__pyx_t_17.memview = NULL;
__pyx_t_17.data = NULL;
/* "squish/energy.pyx":419
* self.hess = -2*self.r*np.asarray(HE, dtype=FLOAT)
@ -8018,7 +7955,7 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
__pyx_t_3 = __Pyx_PyInt_TrueDivideObjC(__pyx_t_5, __pyx_int_2, 2, 0, 0); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 421, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_3);
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
__pyx_t_16 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_nn___pyx_t_6squish_4core_FLOAT_T(__pyx_t_3, PyBUF_WRITABLE); if (unlikely(!__pyx_t_16.memview)) __PYX_ERR(0, 421, __pyx_L1_error)
__pyx_t_17 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_nn___pyx_t_6squish_4core_FLOAT_T(__pyx_t_3, PyBUF_WRITABLE); if (unlikely(!__pyx_t_17.memview)) __PYX_ERR(0, 421, __pyx_L1_error)
__Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
/* "squish/energy.pyx":418
@ -8029,9 +7966,9 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
* + np.asarray(self.hess, dtype=FLOAT).T )
*/
__PYX_XDEC_MEMVIEW(&__pyx_v_self->__pyx_base.hess, 0);
__pyx_v_self->__pyx_base.hess = __pyx_t_16;
__pyx_t_16.memview = NULL;
__pyx_t_16.data = NULL;
__pyx_v_self->__pyx_base.hess = __pyx_t_17;
__pyx_t_17.memview = NULL;
__pyx_t_17.data = NULL;
/* "squish/energy.pyx":273
* self.grad = dedx
@ -8050,7 +7987,7 @@ static void __pyx_f_6squish_6energy_13RadialTEnergy_calc_hess(struct __pyx_obj_6
__Pyx_XDECREF(__pyx_t_4);
__Pyx_XDECREF(__pyx_t_5);
__PYX_XDEC_MEMVIEW(&__pyx_t_6, 1);
__PYX_XDEC_MEMVIEW(&__pyx_t_16, 1);
__PYX_XDEC_MEMVIEW(&__pyx_t_17, 1);
__Pyx_AddTraceback("squish.energy.RadialTEnergy.calc_hess", __pyx_clineno, __pyx_lineno, __pyx_filename);
__pyx_L0:;
__PYX_XDEC_MEMVIEW(&__pyx_v_HE, 1);

2
squish/energy.pyx
View File

@ -285,7 +285,7 @@ cdef class RadialTEnergy(VoronoiContainer):
cdef FLOAT_T [:,:] HE = np.zeros((2*self.n, 2*self.n))
cdef INT_T i, j, k
for i in prange(self.sites.shape[0], nogil=True):
for i in range(self.sites.shape[0]):
xi = _Site(i, &info)
e = xi.edge(&xi)

7
squish/ordered.py
View File

@ -7,14 +7,11 @@ import cmath
from fractions import Fraction
Config = Tuple[int, int]
R_HEX = (6 * 3 ** (-0.25) * sqrt(2) * atanh(0.5)) / (2 * pi)
def e_hex(domain: DomainParams) -> float:
return (
2
- 2 * domain.r * (6 * 3 ** (-0.25) * sqrt(2) * atanh(0.5))
+ 2 * pi * domain.r ** 2
)
return 2 - 2 * domain.r * (2 * pi) * R_HEX + 2 * pi * domain.r ** 2
def configurations(domain: DomainParams) -> List[Config]:

8
squish/simulation.py
View File

@ -217,6 +217,7 @@ class Flow(Simulation):
log: bool,
log_steps: int,
new_sites: Optional[numpy.ndarray] = None,
newton: bool = False,
) -> None:
if log:
print(f"Find - {self.domain}", flush=True)
@ -239,13 +240,14 @@ class Flow(Simulation):
change, grad = frame.iterate(self.step_size)
new_frame = self.energy.mode(*self.domain, frame.add_sites(change))
grad_norm = np.linalg.norm(grad)
end = timer()
grad_norm = np.linalg.norm(grad)
if self.adaptive:
error = change - grad * self.step_size
tol = 10 ** min(-4, -3 + log10(grad_norm))
# tol = 10 ** -7
tol = 10 ** min(-5, 2.3 * log10(grad_norm))
# tol = 10 ** -10
self.step_size *= (tol / np.linalg.norm(error)) ** 0.5
if not save:

4
squish/squish.py
View File

@ -175,9 +175,7 @@ def do_sim(args, file):
)
elif mode == "shrink":
check_params(
sim_params,
["width_change", "width_stop"],
{"width_change": "positive", "width_stop": "positive"},
sim_params, ["width_change", "width_stop"], {"width_stop": "positive"}
)
sim = Shrink(
domain,