From b32e0b1a8c0ce49be61f9984e88720517760133b Mon Sep 17 00:00:00 2001 From: Kenneth Jao Date: Sat, 29 Jan 2022 20:00:53 -0500 Subject: [PATCH] Script updates and minor setting changes --- scripts/defect_energy.py | 101 --------------- scripts/defect_scatter.py | 68 ++++++++++ scripts/density_of_states.py | 75 +++-------- scripts/frustration.py | 173 +++++++++++++++++++++++++ scripts/hexagon_alpha.py | 2 +- scripts/ordered_rbar.py | 70 ++++++++++ scripts/ordered_scatter.py | 74 +++++++++++ scripts/probability_of_disorder.py | 3 +- scripts/script_tools.py | 19 ++- scripts/vee.py | 14 +- scripts/vee_diff_pod.py | 2 +- squish/diagram.py | 11 +- squish/energy.c | 201 ++++++++++------------------- squish/energy.pyx | 2 +- squish/ordered.py | 7 +- squish/simulation.py | 8 +- squish/squish.py | 4 +- 17 files changed, 515 insertions(+), 319 deletions(-) delete mode 100644 scripts/defect_energy.py create mode 100644 scripts/defect_scatter.py create mode 100644 scripts/frustration.py create mode 100644 scripts/ordered_rbar.py create mode 100644 scripts/ordered_scatter.py diff --git a/scripts/defect_energy.py b/scripts/defect_energy.py deleted file mode 100644 index 4dc6d70..0000000 --- a/scripts/defect_energy.py +++ /dev/null @@ -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.") diff --git a/scripts/defect_scatter.py b/scripts/defect_scatter.py new file mode 100644 index 0000000..fc7e8da --- /dev/null +++ b/scripts/defect_scatter.py @@ -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.") diff --git a/scripts/density_of_states.py b/scripts/density_of_states.py index eb27405..60297e0 100644 --- a/scripts/density_of_states.py +++ b/scripts/density_of_states.py @@ -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") diff --git a/scripts/frustration.py b/scripts/frustration.py new file mode 100644 index 0000000..fe030e0 --- /dev/null +++ b/scripts/frustration.py @@ -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.") diff --git a/scripts/hexagon_alpha.py b/scripts/hexagon_alpha.py index b4ca6b6..9e748c3 100644 --- a/scripts/hexagon_alpha.py +++ b/scripts/hexagon_alpha.py @@ -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) diff --git a/scripts/ordered_rbar.py b/scripts/ordered_rbar.py new file mode 100644 index 0000000..08edad8 --- /dev/null +++ b/scripts/ordered_rbar.py @@ -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.") diff --git a/scripts/ordered_scatter.py b/scripts/ordered_scatter.py new file mode 100644 index 0000000..e8b153d --- /dev/null +++ b/scripts/ordered_scatter.py @@ -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.") diff --git a/scripts/probability_of_disorder.py b/scripts/probability_of_disorder.py index 11e20aa..e6efa01 100644 --- a/scripts/probability_of_disorder.py +++ b/scripts/probability_of_disorder.py @@ -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) diff --git a/scripts/script_tools.py b/scripts/script_tools.py index fc847b4..f04042f 100644 --- a/scripts/script_tools.py +++ b/scripts/script_tools.py @@ -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 diff --git a/scripts/vee.py b/scripts/vee.py index 4dec74a..7ce2b23 100644 --- a/scripts/vee.py +++ b/scripts/vee.py @@ -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) diff --git a/scripts/vee_diff_pod.py b/scripts/vee_diff_pod.py index 15aa91a..d32048c 100644 --- a/scripts/vee_diff_pod.py +++ b/scripts/vee_diff_pod.py @@ -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) diff --git a/squish/diagram.py b/squish/diagram.py index 932907d..c320e6e 100644 --- a/squish/diagram.py +++ b/squish/diagram.py @@ -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, diff --git a/squish/energy.c b/squish/energy.c index 2e7bad2..2eb2f96 100644 --- a/squish/energy.c +++ b/squish/energy.c @@ -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); diff --git a/squish/energy.pyx b/squish/energy.pyx index d395640..00c000f 100644 --- a/squish/energy.pyx +++ b/squish/energy.pyx @@ -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) diff --git a/squish/ordered.py b/squish/ordered.py index 0fa581f..52f6900 100644 --- a/squish/ordered.py +++ b/squish/ordered.py @@ -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]: diff --git a/squish/simulation.py b/squish/simulation.py index 65359ef..c01eb30 100644 --- a/squish/simulation.py +++ b/squish/simulation.py @@ -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: diff --git a/squish/squish.py b/squish/squish.py index 7702bf7..9ef39a4 100644 --- a/squish/squish.py +++ b/squish/squish.py @@ -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,