Updated scripts for post-processing data.

scripts/aspect_diagrams.py

@@ -0,0 +1,424 @@
+from __future__ import annotations
+from typing import List, Tuple, Dict
+import argparse, math, numpy as np, os
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
+from multiprocessing import Pool, cpu_count
+from pathlib import Path
+
+import squish.ordered as order
+from squish import Simulation, DomainParams
+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]]:
+    sim, frames = Simulation.load(file)
+
+    alls = []
+    for frame_info in frames:
+        alls.append(
+            [
+                frame_info["energy"],
+                np.var(frame_info["stats"]["avg_radius"]) <= 1e-8,
+                np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6),
+                sum(frame_info["stats"]["site_energies"][: sim.domain.n]),
+            ]
+        )
+
+    sim, frames = Simulation.load(file)
+    sim.frames = list(frames)
+    counts = sim.get_distinct()
+
+    distincts = []
+    for j, frame_info in enumerate(sim.frames):
+        distincts.append(
+            [
+                frame_info["energy"],
+                np.var(frame_info["stats"]["avg_radius"]) <= 1e-8,
+                np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6),
+                sum(frame_info["stats"]["site_energies"][: sim.domain.n]),
+                counts[j],
+            ]
+        )
+
+    return sim.domain.w / sim.domain.h, alls, distincts
+
+
+def get_equilibria_data(filepath: Path) -> Tuple[Dict, numpy.ndarray, DomainParams]:
+    data = {"all": {}, "distinct": {}}
+    files = list(Path(filepath).iterdir())
+
+    with Pool(cpu_count()) as pool:
+        for i, res in enumerate(pool.imap_unordered(eq_file_process, files)):
+            data["all"][res[0]] = res[1]
+            data["distinct"][res[0]] = res[2]
+
+            hashes = int(21 * i / len(files))
+            print(
+                f'Loading simulations... |{"#"*hashes}{" "*(20-hashes)}|'
+                + f" {i+1}/{len(files)} simulations loaded.",
+                flush=True,
+                end="\r",
+            )
+        print(flush=True)
+
+    sim, frames = Simulation.load(files[0])
+    widths = np.asarray(sorted(data["all"]))
+    domain = DomainParams(sim.domain.n, widths[-1], sim.domain.h, sim.domain.r)
+    return data, widths, domain
+
+
+def axis_settings(ax, widths):
+    ax.grid(zorder=0)
+    ax.set_xticks([round(w, 2) for w in widths[::2]])
+    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)
+
+
+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])
+        )
+
+    ax.plot(widths, all_disorder_count)
+    axis_settings(ax, widths)
+
+    ax.yaxis.set_major_formatter(mtick.PercentFormatter())
+    ax.title.set_text(f"Probability of Disorder - N{domain.n}")
+    ax.set_xlabel("Aspect Ratio")
+    ax.set_ylabel("Disordered Equilibria")
+    boa_y_min = round(min(all_disorder_count) / 20) * 20 - 5
+    ax.set_yticks(np.arange(boa_y_min, 100.01, 2.5))
+
+    return fig
+
+
+def density_of_states(data, widths, domain):
+    fig, ax = plt.subplots(figsize=(16, 8))
+    distinct_ordered, distinct_unordered = [], []
+    for width in widths:
+        equal_shape = list([c[1] for c in data["distinct"][width]])
+        distinct_ordered.append(equal_shape.count(True))
+        distinct_unordered.append(equal_shape.count(False))
+
+    ax2 = ax.twinx()
+    ax.plot(widths, distinct_unordered, label="Unordered Equilibria", color="C0")
+    ax2.plot(widths, distinct_ordered, label="Ordered Equilibria", color="C1")
+    axis_settings(ax, widths)
+    plt.subplots_adjust(0.07, 0.12, 0.92, 0.9)
+    ax.title.set_text(f"Density of States - N{domain.n}")
+    ax.set_xlabel("Aspect Ratio")
+    ax.set_ylabel("Number of States (Disordered)", color="C0")
+    ax2.set_ylabel("Number of States (Ordered)", color="C1")
+
+    dos_y_max_unorder = 1.05 * max(distinct_unordered)
+    dos_y_max_order = 1.05 * max(distinct_ordered)
+    ax.set_yticks(np.linspace(0, dos_y_max_unorder, 20).astype(int))
+    # ax.set_yticks(np.arange(0, dos_y_max_unorder, round(dos_y_max_unorder/200, 1)*10))
+    ax2.set_yticks(np.arange(0, dos_y_max_order))
+
+    return fig
+
+
+def defect_density(data, widths, domain):
+    fig, ax = plt.subplots(figsize=(16, 8))
+
+    defects = []
+    for width in widths:
+        defects.append(
+            sum([c[2] for c in data["all"][width] if not c[1]])
+            / len(data["all"][width])
+        )
+
+    ax.plot(widths, defects)
+    axis_settings(ax, widths)
+    ax.title.set_text(f"Average Defects - N{domain.n}")
+    ax.set_xlabel("Aspect Ratio")
+    ax.set_ylabel("Defects")
+    ax.set_yticks(np.arange(0, 1 + max(defects), 0.5))
+
+    return fig
+
+
+def circle_isoparam(data, widths, order_data, domain):
+    fig, ax = plt.subplots(figsize=(16, 8))
+
+    ax2 = ax.twinx()
+    axis_settings(ax, widths)
+    plt.subplots_adjust(0.07, 0.12, 0.92, 0.9)
+    ax.title.set_text(f"Circular Isoparametric Ratio - N{domain.n}")
+    ax.set_xlabel("Aspect Ratio")
+    ax.set_ylabel("Maximum Ratio", color="C0")
+    ax2.set_ylabel("Minimum Ratio", color="C1")
+
+    ax.plot(widths, order_data["isoparam_max"], label="Maximum", color="C0")
+    ax2.plot(widths, order_data["isoparam_min"], label="Minimum", color="C1")
+
+    return fig
+
+
+def reduced_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])
+
+    offset = np.array(min_order)
+
+    min_unorder_off = min_unorder - offset
+    max_unorder_off = max_unorder - offset
+    ax.plot(widths, min_order - offset, color="C1")
+    # ax.plot(widths, max_order - offset, color='C1', linestyle='dotted')
+    ax.plot(widths, min_unorder_off, color="C0")
+    ax.plot(widths, max_unorder_off, color="C0", linestyle="dotted")
+    axis_settings(ax, widths)
+
+    ax.title.set_text(f"Reduced Energy vs. Width - N{domain.n}")
+    ax.set_xlabel("Aspect Ratio")
+    ax.set_ylabel("Reduced Energy")
+    bif_y_max = np.max(np.abs(np.concatenate((min_unorder_off, max_unorder_off))))
+    bif_top = np.arange(
+        0, bif_y_max, round(bif_y_max / 20, -math.floor(math.log10(bif_y_max / 20)))
+    )
+    ax.set_yticks(np.concatenate((-bif_top[1:][::-1], bif_top)))
+
+    return fig
+
+
+def defect_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])
+
+    offset = np.array(min_order)
+
+    defect_a, defect_b = [], []
+    for width in widths:
+        num_defects = [c[2] for c in data["all"][width]]
+        defect_energy = [c[3] for c in data["all"][width]]
+        m, b = np.polyfit(num_defects, defect_energy, 1)
+
+        defect_a.append(m)
+        defect_b.append(b)
+
+    ax2 = ax.twinx()
+    ax.plot(widths, defect_a, label="Energy per Defect", color="C0")
+    ax2.plot(widths, defect_b - offset, label="Relative Initial Energy", color="C1")
+    axis_settings(ax, widths)
+    plt.subplots_adjust(0.07, 0.12, 0.92, 0.9)
+    ax.title.set_text(f"Defect Energy - N{domain.n}")
+    ax.set_xlabel("Aspect Ratio")
+    ax.set_ylabel("Energy per Defect", color="C0")
+    ax2.set_ylabel("Relative Initial Energy", color="C1")
+
+    return fig
+
+
+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])
+
+    # Energy of regular hexagon with area 1
+    offset = (
+        2
+        - 2 * domain.r * (6 * 3 ** (-0.25) * math.sqrt(2) * math.atanh(0.5))
+        + 2 * math.pi * domain.r ** 2
+    )
+
+    min_order_off = min_order / domain.n - offset
+    min_unorder_off = min_unorder / domain.n - offset
+    max_unorder_off = max_unorder / domain.n - offset
+
+    ax.plot(widths, min_order_off, color="C1", label="Minimum Ordered")
+    ax.plot(widths, min_unorder_off, color="C0", label="Minimum Disordered")
+    ax.plot(
+        widths,
+        max_unorder_off,
+        color="C0",
+        linestyle="dotted",
+        label="Maximum Disordered",
+    )
+    # ax.plot(
+    #     [min(widths), max(widths)],
+    #     [offset, offset],
+    #     color="C1",
+    #     linestyle="dotted",
+    #     label="Regular Energy",
+    # )
+
+    axis_settings(ax, widths)
+    ax.title.set_text(f"Energy at Aspect Ratios - N{domain.n}")
+    ax.set_xlabel("Aspect Ratio")
+    ax.set_ylabel("Excess Energy per Site")
+    ax.legend()
+
+    start, end = ax.get_ylim()
+    ax.set_yticks(np.linspace(0, end, 20))
+    ax.ticklabel_format(axis="y", style="sci")
+
+    return fig
+
+
+def main():
+    # Loading arguments.
+    parser = argparse.ArgumentParser("Outputs width search data into diagrams")
+    parser.add_argument(
+        "sims_path",
+        metavar="path/to/data",
+        help="folder that contains simulation files, or cached data file.",
+    )
+    parser.add_argument(
+        "-q",
+        "--quiet",
+        dest="quiet",
+        action="store_true",
+        default=False,
+        help="suppress all normal output",
+    )
+
+    args = parser.parse_args()
+
+    # Obtain data from simulation files and generate single shape data.
+    data, widths, domain = get_equilibria_data(Path(args.sims_path))
+    order_data = get_ordered_energies(domain, widths)
+
+    fig_folder = OUTPUT_DIR / Path(f"AspectDiagrams - N{domain.n}")
+    fig_folder.mkdir(exist_ok=True)
+
+    # Generating diagrams.
+    probability_of_disorder(data, widths, domain).savefig(
+        fig_folder / "Probability of Disorder.png"
+    )
+
+    density_of_states(data, widths, domain).savefig(
+        fig_folder / "Density Of States.png"
+    )
+
+    defect_density(data, widths, domain).savefig(fig_folder / "Defects.png")
+
+    reduced_energy(data, widths, order_data, domain).savefig(
+        fig_folder / "Reduced Energy.png"
+    )
+
+    defect_energy(data, widths, order_data, domain).savefig(
+        fig_folder / "Defect Energy.png"
+    )
+
+    circle_isoparam(data, widths, order_data, domain).savefig(
+        fig_folder / "Circular Isoparametric Ratio.png"
+    )
+
+    excess_energy(data, widths, order_data, domain).savefig(
+        fig_folder / "Excess Energy.png"
+    )
+
+    print(f"Wrote to {fig_folder}.")
+
+
+if __name__ == "__main__":
+    os.environ["QT_LOGGING_RULES"] = "*=false"
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Program terminated by user.")

scripts/coercivity.py

@@ -0,0 +1,113 @@
+from __future__ import annotations
+from typing import List, Tuple, Dict
+import argparse, math, numpy as np, os, pickle
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
+from multiprocessing import Pool, cpu_count
+from pathlib import Path
+
+import squish.ordered as order
+from squish import Simulation, DomainParams
+from squish.common import Energy, OUTPUT_DIR
+
+
+def axis_settings(ax, widths):
+    ax.invert_xaxis()
+    ax.grid(zorder=0)
+    ax.set_xticks([round(w, 2) for w in widths[::-2]])
+    ax.set_xticklabels(ax.get_xticks(), rotation=90)
+    plt.subplots_adjust(0.07, 0.12, 0.97, 0.9)
+
+
+def main():
+    # Loading arguments.
+    parser = argparse.ArgumentParser("Outputs ordered equilibria lowest eigenvalues.")
+    parser.add_argument(
+        "n_objects",
+        metavar="N",
+        type=int,
+        help="folder that contains simulation files, or cached data file.",
+    )
+    parser.add_argument(
+        "-q",
+        "--quiet",
+        dest="quiet",
+        action="store_true",
+        default=False,
+        help="suppress all normal output",
+    )
+
+    args = parser.parse_args()
+
+    widths = np.round(np.linspace(3.0, 10.0, 141), 2)
+
+    values = []
+    with open("coercivity_eigs.pkl", "rb") as f:
+        store_data = pickle.load(f)
+
+    for width in widths:
+        eig_vals = []
+        for config, eigs in store_data[width].items():
+            zero_ind = np.where(np.isclose(eigs, 0, atol=1e-4))[0][0]
+            if zero_ind == 0:
+                eig_vals.append(eigs[2])
+            else:
+                eig_vals.append(eigs[0])
+
+        values.append(min(eig_vals))
+
+    # for i, width in enumerate(widths):
+    # 	domain = DomainParams(args.n_objects, width, 10, 4.0)
+    # 	eig_vals = []
+    # 	store_data[width] = {}
+    # 	configs = order.configurations(domain)
+    # 	for j, config in enumerate(configs):
+    # 		if config == (1,0):
+    # 			continue
+    # 		points = order.sites(domain, config)
+
+    # 		hess = Energy("radial-t").mode(*domain, points).hessian(10e-5)
+    # 		eigs = np.sort(np.linalg.eig(hess)[0])[::-1]
+    # 		store_data[width][config] = eigs
+
+    # 		zero_ind = np.where(np.isclose(eigs, 0))[0][0]
+    # 		if zero_ind == 0:
+    # 			eig_vals.append(eigs[2])
+    # 		else:
+    # 			eig_vals.append(eigs[0])
+
+    # 		hashes = int(21*j/len(widths))
+    # 		print(f'Generating at {width}, {i+1}/{len(widths)}... |{"#"*hashes}{" "*(20-hashes)}|' + \
+    # 				f' {j+1}/{len(configs)} configs done.', flush=True, end='\r')
+
+    # print(flush=True)
+
+    # with open("coercivity_eigs.pkl", "wb") as f:
+    # 	pickle.dump(store_data, f, pickle.HIGHEST_PROTOCOL)
+
+    # return
+    fig, ax = plt.subplots(figsize=(12, 8))
+    plt.subplots_adjust(0.07, 0.12, 0.97, 0.9)
+
+    ax.plot(widths, values)
+
+    ax.invert_xaxis()
+    ax.grid(zorder=0)
+    ax.set_xticks([round(w, 2) for w in widths[::-2]])
+    ax.set_xticklabels(ax.get_xticks(), rotation=90)
+
+    fig.suptitle("Coercivity")
+    # ax.set_xlim([0, 5])
+    ax.set_xlabel("Width")
+    ax.set_ylabel("Eigenvalue")
+
+    fig.savefig(OUTPUT_DIR / "Coercivity.png")
+    print(f"Wrote to {OUTPUT_DIR / 'Coercivity.png'}")
+
+
+if __name__ == "__main__":
+    os.environ["QT_LOGGING_RULES"] = "*=false"
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Program terminated by user.")

scripts/convergence.py

@@ -17,14 +17,6 @@ def main():
         metavar="path/to/data",
         help="folder that contains simulation files at various step sizes.",
     )
-    parser.add_argument(
-        "-q",
-        "--quiet",
-        dest="quiet",
-        action="store_true",
-        default=False,
-        help="suppress all normal output",
-    )
 
     args = parser.parse_args()
     data = {}

scripts/cumulative_vee.py

@@ -72,20 +72,20 @@ def main():
 
         for i, count in enumerate(counts):
             my_cool_data[i + 1].append(count)
-
+        ax.plot(100*vees, counts, label=f"N={domain.n}")
-        ax.plot(
+        #ax.plot(
-            100 * vees[: index + 1],
+        #    100 * vees[: index + 1],
-            counts[: index + 1],
+        #    counts[: index + 1],
-            label=f"N={domain.n}",
+        #    label=f"N={domain.n}",
-            color=f"C{j}",
+        #    color=f"C{j}",
-        )
+        #)
-        ax.plot(
+        #ax.plot(
-            100 * vees[index:],
+        #    100 * vees[index:],
-            counts[index:],
+        #    counts[index:],
-            label=f"_nolegend_",
+        #    label=f"_nolegend_",
-            linestyle="dotted",
+        #    linestyle="dotted",
-            color=f"C{j}",
+        #    color=f"C{j}",
-        )
+        #)
 
     with open("cumulative-vee.csv", "w") as csvfile:
         writer = csv.writer(csvfile)

scripts/cumulative_vee_2.py

@@ -167,8 +167,6 @@ def main():
     pad = 100
     alpha_prob_dens = {}
     for alpha, avg in alpha_avgs.items():
-        if alpha == 1.0:
-            continue
         mov_avgs = np.array(
             [np.mean(avg[max(0, i - pad) : i + pad]) for i in range(len(avg))]
         )
@@ -184,17 +182,21 @@ def main():
         d_avgs = np.gradient(mov_avgs, vees[1] - vees[0])
         a1_prob_dens[n] = d_avgs / (np.sum(d_avgs) * (vees[1] - vees[0]))
     
+    print(list(alpha_prob_dens.keys()))
     for alpha, prob_dens in sorted(alpha_prob_dens.items()):
         if alpha in [0.3, 0.5, 0.7, 0.9, 1.0]:
-            continue
+            # continue
             ax.plot(100 * vees, prob_dens, label=fr"$\alpha={alpha:.1f}$")
+    ax.plot(100*vees, a1_prob_dens[395], zorder=50, label=r"$\alpha = 1.0$", color="black")
     
-    for n, prob_dens in sorted(a1_prob_dens.items()):
+    # Probability densities of alpha=1 across N
-        if n in [255, 315, 335, 355, 375]:
+
-            ax.plot(100 * vees, prob_dens, label=f"N={n}")
+    #for n, prob_dens in sorted(a1_prob_dens.items()):
-        if n == 395:
+    #    if n in [255, 315, 335, 355, 375]:
-            ax.plot(100 * vees, prob_dens, label=f"N={n}", color="black")
+    #        ax.plot(100 * vees, prob_dens, label=f"N={n}")
-            # ax.plot(100 * vees, prob_dens, label=r"$\alpha = 1.0$", zorder=50, color="black")
+    #    if n == 395:
+    #        ax.plot(100 * vees, prob_dens, label=f"N={n}", color="black")
+    #        # ax.plot(100 * vees, prob_dens, label=r"$\alpha = 1.0$", zorder=50, color="black")
     # ax.plot(100 * vees, avgs[395], zorder=50, color="black")
 
     ax.set_xlim(0, 6.3)

scripts/defectenergy.py

@@ -0,0 +1,39 @@
+from squish import Simulation
+import matplotlib.pyplot as plt
+import os, numpy as np
+
+
+def main():
+    sim, frames = Simulation.load(
+        "squish_output/Radial[T]Search - N11-400 - 10.00x10.00 - 500/Radial[T]Search - N397 - 10.00x10.00"
+    )
+
+    defect, energy = [], []
+    for frame_info in frames:
+        defect.append(np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6))
+        energy.append(sum(frame_info["stats"]["site_energies"][:400]))
+
+    fig, ax = plt.subplots(1, figsize=(8, 8))
+    plt.subplots_adjust(0.1, 0.12, 0.97, 0.9)
+
+    fig.suptitle("Defects vs. Energy")
+    ax.set_xlabel("Defects")
+    ax.set_ylabel("Energy")
+    ax.grid(zorder=0)
+    ax.set_xticks(np.arange(0, 64, 2))
+    ax.scatter(defect, energy, zorder=3, color="C0", marker="*")
+
+    m, b = np.polyfit(defect, energy, 1)
+    ax.plot(
+        defect, np.array(defect) * m + b, zorder=3, color="C1", label=f"Slope: {m:.4f}"
+    )
+    ax.legend()
+    fig.savefig("DefectEnergyN397-10.00.png")
+
+
+if __name__ == "__main__":
+    os.environ["QT_LOGGING_RULES"] = "*=false"
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Program terminated by user.")

scripts/defects.py

@@ -0,0 +1,108 @@
+from __future__ import annotations
+from typing import List
+import argparse, pickle, numpy as np, os
+from pathlib import Path
+import matplotlib.pyplot as plt
+
+from squish import Simulation
+from squish.common import OUTPUT_DIR
+
+
+def main():
+    parser = argparse.ArgumentParser("Graphs average defects at N.")
+    parser.add_argument(
+        "sims_path",
+        metavar="path/to/data",
+        help="folder that contains simulation files at various Ns.",
+    )
+    parser.add_argument(
+        "-q",
+        "--quiet",
+        dest="quiet",
+        action="store_true",
+        default=False,
+        help="suppress all normal output",
+    )
+
+    args = parser.parse_args()
+    data = {}
+
+    files = list(Path(args.sims_path).iterdir())
+    for i, file in enumerate(files):
+        sim, frames = Simulation.load(file)
+        avg_defects = 0
+        count = 0
+
+        for frame in frames:
+            if np.var(frame["stats"]["avg_radius"]) > 1e-8:
+                avg_defects += np.count_nonzero(frame["stats"]["site_edge_count"] != 6)
+
+            count += 1
+
+        avg_defects /= 1 if count == 0 else count
+        data[sim.domain.n] = avg_defects
+
+        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)
+
+    data = sorted(data.items())
+    ns, defects = np.array([x[0] for x in data]), np.array([x[1] for x in data])
+
+    corrected = []
+    for i, x in enumerate(defects):
+        if x == 0:
+            corrected.append(defects[i + 1])
+        else:
+            corrected.append(x)
+
+    fig, ax = plt.subplots(1, 2, figsize=(16, 8))
+    plt.subplots_adjust(0.07, 0.12, 0.97, 0.9)
+
+    fig.suptitle("Defects vs. Number of Sites (N)")
+
+    m0, b0 = np.polyfit(ns, defects, 1)
+
+    ax[0].plot(ns, defects)
+    ax[0].plot(ns, m0 * ns + b0, label=f"Slope: {m0:.5f}")
+    ax[0].grid(zorder=0)
+    ax[0].legend()
+    ax[0].set_xlabel("N")
+    ax[0].set_ylabel("Average Defects")
+
+    x, y = np.log(ns), np.log(corrected)
+    m, b = np.polyfit(x, y, 1)
+
+    x2, y2 = x[40:], np.log(defects[40:])
+    m2, b2 = np.polyfit(x2, y2, 1)
+
+    ax[1].plot(x, y, linestyle="dotted", color="C0")
+    ax[1].plot(x, np.log(defects))
+
+    # ax[1].plot(x, m*x+b, label=f"All N: {m:.5f}")
+    ax[1].plot(x2, m2 * x2 + b2, label=f"N $\\geq$ 40: {m2:.5f}")
+    ax[1].grid(zorder=0)
+    ax[1].legend()
+    ax[1].set_xticks(np.arange(3.75, 6.25, 0.25))
+    ax[1].set_yticks(np.arange(0, 4.5, 0.5))
+    ax[1].set_xlim(3.75, 6)
+    ax[1].set_ylim(0, 4)
+    ax[1].set_xlabel("$\\ln$ N")
+    ax[1].set_ylabel("$\\ln$ Average Defects")
+
+    fig.savefig(OUTPUT_DIR / "DefectsN.png")
+    print(f"Wrote to {OUTPUT_DIR / 'DefectsN.png'}")
+
+
+if __name__ == "__main__":
+    os.environ["QT_log10GING_RULES"] = "*=false"
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Program terminated by user.")

scripts/frustration_n.py

@@ -102,21 +102,21 @@ def main():
     )
     ax.plot(ns, m0 * ns + b0, color="C0", linestyle="dashed")
 
-    lns2 = ax.plot(
+    #lns2 = ax.plot(
-        ns, 100 * epds, color="C1", alpha=0.5, label=r"$\zeta_1 \times 10^{4}$"
+    #    ns, 100 * epds, color="C1", alpha=0.5, label=r"$\zeta_1 \times 10^{4}$"
-    )
+    #)
-
+    #
-    lns3 = ax.plot(
+    #lns3 = ax.plot(
-        ns,
+    #    ns,
-        10 * avg_defects * epds,
+    #    10 * avg_defects * epds,
-        color="C2",
+    #    color="C2",
-        alpha=0.5,
+    #    alpha=0.5,
-        label=r"$\mathcal{F} \times 10^{3}$",
+    #    label=r"$\mathcal{F} \times 10^{3}$",
-    )
+    #)
     # ax3.plot(ns, 100*(m1 * ns + b1) / 10, color="C2", linestyle="dashed")
     ax.set_xlim(93, 407)
 
-    ax.set_ylim(3, 37)
+    ax.set_ylim(3, 36)
     ax.set_yticks(np.arange(5, 40, 5))
 
     # ax2.set_ylim(-3 * 0.4, 18 + 3 * 0.4)
@@ -131,9 +131,9 @@ def main():
     # labs = [l.get_label() for l in lns]
     ax.grid(zorder=0)
     # ax.legend(lns, labs, loc="lower right")
-    ax.legend()
+    #ax.legend()
     ax.set_xlabel("N")
-    # ax.set_ylabel(r"$\langle \mathrm{D} \rangle$")
+    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")

scripts/heatmap.py

@@ -0,0 +1,302 @@
+from __future__ import annotations
+from typing import List, Tuple, Dict
+import argparse, math, numpy as np, os, pickle
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
+from scipy.optimize import curve_fit
+from multiprocessing import Pool, cpu_count
+from pathlib import Path
+
+import squish.ordered as order
+from squish import Simulation, DomainParams
+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]]:
+    sim, frames = Simulation.load(file)
+
+    alls = []
+    for frame_info in frames:
+        alls.append(
+            [
+                frame_info["energy"],
+                np.var(frame_info["stats"]["avg_radius"]) <= 1e-8,
+                np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6),
+                sum(frame_info["stats"]["site_energies"][: sim.domain.n]),
+            ]
+        )
+
+    sim, frames = Simulation.load(file)
+    sim.frames = list(frames)
+    counts = sim.get_distinct()
+
+    distincts = []
+    for j, frame_info in enumerate(sim.frames):
+        distincts.append(
+            [
+                frame_info["energy"],
+                np.var(frame_info["stats"]["avg_radius"]) <= 1e-8,
+                np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6),
+                sum(frame_info["stats"]["site_energies"][: sim.domain.n]),
+                counts[j],
+            ]
+        )
+
+    return sim.domain.w / sim.domain.h, alls, distincts
+
+
+def get_equilibria_data(filepath: Path) -> Tuple[Dict, numpy.ndarray, DomainParams]:
+    data = {"all": {}, "distinct": {}}
+    files = list(Path(filepath).iterdir())
+
+    with Pool(cpu_count()) as pool:
+        for i, res in enumerate(pool.imap_unordered(eq_file_process, files)):
+            data["all"][res[0]] = res[1]
+            data["distinct"][res[0]] = res[2]
+
+            hashes = int(21 * i / len(files))
+            print(
+                f'Loading simulations... |{"#"*hashes}{" "*(20-hashes)}|'
+                + f" {i+1}/{len(files)} simulations loaded.",
+                flush=True,
+                end="\r",
+            )
+        print(flush=True)
+
+    sim, frames = Simulation.load(files[0])
+    widths = np.asarray(sorted(data["all"]))
+    domain = DomainParams(sim.domain.n, widths[-1], sim.domain.h, sim.domain.r)
+    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():
+    # Loading arguments.
+    parser = argparse.ArgumentParser("Outputs width search data into diagrams")
+    parser.add_argument(
+        "sims_path",
+        metavar="path/to/data",
+        help="folder that contains simulation files of all searches for all N.",
+    )
+    parser.add_argument(
+        "-q",
+        "--quiet",
+        dest="quiet",
+        action="store_true",
+        default=False,
+        help="suppress all normal output",
+    )
+
+    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:
+    # 	disorder_dict = pickle.load(f)
+    # 	widths = np.linspace(3.0, 10.0, 141)
+    # 	min_n, max_n = 60, 80
+
+    disorder_dict = {}
+    for file in Path(args.sims_path).iterdir():
+        sim_data, widths, domain = get_equilibria_data(file)
+
+        disorder_count = []
+        for width in widths:
+            equal_shape = list([c[1] for c in sim_data["all"][width]])
+            disorder_count.append(
+                100 * equal_shape.count(False) / len(sim_data["all"][width])
+            )
+
+        disorder_dict[domain.n] = disorder_count
+
+    min_n, max_n = min(disorder_dict), max(disorder_dict)
+    filepath = f"Disorder Heatmap N{min_n}-{max_n}"
+
+    # with open("testing.pkl", "wb") as f:
+    # 	pickle.dump(disorder_dict, f, pickle.HIGHEST_PROTOCOL)
+
+    disorder_arr = np.zeros((max_n - min_n + 1, len(widths)))
+    for key, value in disorder_dict.items():
+        disorder_arr[key - min_n] = np.asarray(value)
+
+    fig, ax = plt.subplots(figsize=(12, 8))
+
+    extent = [min(widths), max(widths), min_n, max_n + 1]
+    ax.imshow(
+        disorder_arr,
+        cmap="plasma",
+        interpolation="nearest",
+        aspect="auto",
+        extent=extent,
+    )
+
+    ax.invert_xaxis()
+    ax.set_xticks([round(w, 2) for w in widths[::-2]])
+    ax.set_xticklabels(ax.get_xticks(), rotation=90)
+    ax.set_yticks(list(range(min_n, max_n + 1)))
+    plt.subplots_adjust(0.07, 0.12, 0.97, 0.9)
+
+    ax.title.set_text(filepath)
+    ax.set_xlabel("Width")
+    ax.set_ylabel("Number of Sites")
+    fig.savefig(OUTPUT_DIR / filepath)
+
+    print(f"Wrote to {OUTPUT_DIR / filepath}.")
+
+
+if __name__ == "__main__":
+    os.environ["QT_LOGGING_RULES"] = "*=false"
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Program terminated by user.")

scripts/min_order_vee.py

@@ -1,5 +1,6 @@
 import numpy as np, os, math
 import matplotlib.pyplot as plt
+from scipy.optimize import curve_fit
 from multiprocessing import Pool, cpu_count
 
 from squish import Simulation, ordered, DomainParams
@@ -64,17 +65,32 @@ def main():
         data["Config"],
     )
 
-    print(min_order, min_config)
-
     fig = plt.figure(figsize=(20, 15))
     gs = fig.add_gridspec(1, 1)
     ax = fig.add_subplot(gs[0])
 
-    ax.scatter(ns, min_order)
+    pad = 5
+    mov_avgs = np.array(
+        [np.mean(min_order[max(0, i - pad) : i + pad]) for i in range(len(min_order))]
+    )
+
+    ax.scatter(ns, mov_avgs, s=8)
+
+    log_slope, _ = np.polyfit(np.log10(ns), np.log10(mov_avgs), 1)
+    print(log_slope)
+
+    def g(x, k):
+        return k * (x ** log_slope)
+
+    # return a * np.e ** (-b * (x - c))
+    params, covs = curve_fit(g, ns, min_order, p0=[5000])
+    ax.plot(np.arange(500, 5500), g(np.arange(500, 5500), *params), color="C1")
+
+    print(params)
 
     ax.grid(zorder=0)
 
-    # ax.set_xlim(1000, 5000)
+    ax.set_xlim(800, 5200)
 
     ax.set_xlabel("N")
     ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$")

scripts/near_neighbors.py

@@ -0,0 +1,115 @@
+import numpy as np, os
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
+
+from squish import Simulation, DomainParams, ordered
+from squish.common import OUTPUT_DIR
+from script_tools import RC_SETTINGS, get_args, get_data
+
+NAME = "NearNeighbors"
+second_neigh = True
+
+
+def main():
+    sims_path, regen = get_args(
+        "Probability distribution of near neighbrs at distances",
+        "simulation folder of equilibriums to plot",
+    )
+
+    dists = np.linspace(0, 5, 25000)
+
+    def f():
+        all_data = {}
+        for n in [400]:
+            if second_neigh:
+                sim = Simulation.from_file(sims_path / f"Radial[T]Search - N{n} - 2500")
+                print("Loaded simulation.")
+            else:
+                sim, frames = Simulation.load(
+                    sims_path / f"Radial[T]Search - N{n} - 2500"
+                )
+
+            counts = np.empty(dists.shape, dtype=float)
+            total_dists = np.array([], dtype=float)
+            if second_neigh:
+                for i, frame in enumerate(sim.frames):
+                    if np.count_nonzero(frame.stats["site_edge_count"] != 6) != 6:
+                        continue
+                    total_dists = np.append(total_dists, frame.stats["site_distances"])
+                    snn = frame.second_near_neighbor()
+
+                    snn_dists = np.empty((sum([len(x) for x in snn])), dtype=float)
+                    ind = 0
+                    for j, site_snn in enumerate(snn):
+                        site_snn_dists = ordered.toroidal_distance(
+                            sim.domain,
+                            frame.site_arr[site_snn],
+                            np.ones((len(site_snn), 2)) * frame.site_arr[j],
+                        )
+                        snn_dists[ind : ind + len(site_snn_dists)] = site_snn_dists
+                        ind += len(site_snn_dists)
+
+                    total_dists = np.append(total_dists, snn_dists)
+
+                    hashes = int(21 * i / len(sim))
+                    print(
+                        f'Processing N={n} at frame {i}... |{"#"*hashes}{" "*(20-hashes)}|'
+                        + f" {i+1}/{len(sim)} completed.",
+                        flush=True,
+                        end="\r",
+                    )
+                print(flush=True)
+            else:
+                for frame in frames:
+                    if (
+                        np.count_nonzero(frame["stats"]["site_edge_count"] != 6) / n
+                        < 0.08
+                    ):
+                        continue
+                    total_dists = np.append(
+                        total_dists, frame["stats"]["site_distances"]
+                    )
+
+            for i, dist in enumerate(dists):
+                counts[i] = np.count_nonzero(total_dists <= dist)
+            counts = 100 * counts / len(total_dists)
+
+            pad = 20
+            counts = np.array(
+                [np.mean(counts[max(0, i - pad) : i + pad]) for i in range(len(counts))]
+            )
+
+            grad = np.gradient(counts, dists[1] - dists[0])
+            all_data[n] = grad / (np.sum(grad) * (dists[1] - dists[0]))
+
+        return all_data
+
+    all_data = get_data(sims_path / "NearNeighborsD6.pkl", f, regen=regen)
+    data_2 = get_data(sims_path / "NearNeighborsD60.pkl", f, regen=regen)
+
+    plt.rcParams.update(RC_SETTINGS)
+
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    for n, counts in sorted(all_data.items()):
+        ax.plot(dists, counts, label=f"D=6")
+
+    for n, counts in sorted(data_2.items()):
+        ax.plot(dists, counts, label=f"D=60")
+
+        ax.set_xlim(0, 3)
+
+    ax.set_xlabel(r"Distance")
+    # ax.set_ylabel("Percent of Distances")
+
+    ax.grid(zorder=0)
+    ax.legend()
+
+    fig.savefig(OUTPUT_DIR / (NAME + ".png"))
+    print(f"Wrote to {OUTPUT_DIR / (NAME + '.png')}")
+
+
+if __name__ == "__main__":
+    main()

scripts/perturbations.py

@@ -2,11 +2,14 @@ import argparse, numpy as np, os
 from pathlib import Path
 import matplotlib.pyplot as plt
 
-from squish import DomainParams, Simulation, ordered
+from squish import ordered
-from squish.common import OUTPUT_DIR
+from squish.common import OUTPUT_DIR, DomainParams, Energy
+from squish.simulation import Simulation, Flow
+
 from script_tools import RC_SETTINGS, get_data, format_data
 
 NAME = "Perturbations"
+NAME2 = "MinimumEscapeVEE"
 
 
 def main():
@@ -14,9 +17,7 @@ def main():
         description="Graphs perturbation graphs for a collection of simulations."
     )
     parser.add_argument(
-        "sims_path",
+        "sim_store_path", metavar="sim_dir", help="folder to save simulations to"
-        metavar="sim_dir",
-        help="folder that contains of perturbations from an equilibrium",
     )
     parser.add_argument(
         "end_path",
@@ -32,41 +33,122 @@ def main():
     )
 
     args = parser.parse_args()
-    sims_path = Path(args.sims_path)
+    out_fol = Path(args.sim_store_path)
 
-    end_sim = Simulation.from_file(args.end_path)
+    (OUTPUT_DIR / out_fol).mkdir(exist_ok=True)
+
+    end_eq = Simulation.from_file(args.end_path)
+    e_hex = ordered.e_hex(end_eq.domain)
 
     def f():
-        data = {}
+        all_data = []
-        for file in sims_path.iterdir():
+        for j in range(50):
-            if "k" not in file.name or file.is_file():
+            pert_out_fol = out_fol / f"Vector{j:03}"
-                continue
+            (OUTPUT_DIR / pert_out_fol).mkdir(exist_ok=True)
-            k = float(file.name.split("k")[-1])
-            delta = 10 ** k
 
-            sim, frames = Simulation.load(file)
+            perturb = np.random.random_sample(end_eq.frames[0].site_arr.shape)
-            data[delta] = {"norm": [], "time": [], "k": k}
+            perturb /= np.linalg.norm(perturb)
+
+            data = {}
+            k = -4
+            same_eq = True
+            while same_eq:
+                print(f"Testing k={k} on vector {j:03}")
+                k_out_fol = pert_out_fol / f"EQk{k}"
+                delta = 10 ** k
+                data[delta] = {"norm": [], "time": [], "vee": [], "k": k}
+                if not pert_out_fol.is_dir():
+                    this_pert = perturb * delta
+                    sim = Flow(
+                        end_eq.domain,
+                        end_eq.energy,
+                        end_eq.step_size,
+                        end_eq.thres,
+                        True,
+                        name=k_out_fol,
+                    )
+                    sim.run(True, True, 50, end_eq.frames[0].site_arr + this_pert)
+
+                sim, frames = Simulation.load(OUTPUT_DIR / k_out_fol)
 
                 for i, frame in enumerate(frames):
                     adjusted = frame["arr"] + (
-                    end_sim.frames[0].site_arr[0] - frame["arr"][0]
+                        end_eq.frames[0].site_arr[0] - frame["arr"][0]
                     )
 
                     data[delta]["norm"].append(
                         np.linalg.norm(
                             ordered.toroidal_distance(
-                            end_sim.domain, adjusted, end_sim.frames[0].site_arr
+                                end_eq.domain, adjusted, end_eq.frames[0].site_arr
                             )
                         )
                     )
                     data[delta]["time"].append(sim.step_size * i)
+                    data[delta]["vee"].append(frame["energy"] / sim.domain.n - e_hex)
 
-        return data
+                k += 1 if k < 0 else 0.25
+                m, _ = np.polyfit(data[delta]["time"], data[delta]["norm"], 1)
+                same_eq = m < 0
 
-    data = get_data(sims_path / "PerturbData.pkl", f, regen=args.regen)
+            all_data.append({"vec": perturb, "data": data})
+
+        return all_data
+
+    all_data = get_data(OUTPUT_DIR / out_fol / "PerturbData.pkl", f, regen=args.regen)
+
+    end_vee = end_eq.frames[0].energy / end_eq.domain.n - e_hex
+    print(end_vee)
+    vees = []
+    for dat in all_data:
+        for k, v in sorted(dat["data"].items()):
+            if v["norm"][-1] > 1e-3:
+                vees.append(dat["data"][k]["vee"][0])
+                break
+
+    eigs = np.sort(np.linalg.eigvalsh(end_eq.frames[0].hessian))
+
+    zero_ind = np.where(np.isclose(eigs, 0, atol=1e-8))[0]
+    if len(zero_ind) == 0:
+        coer = eigs[0]
+    elif zero_ind[0] == 0:
+        coer = eigs[2]
+    else:
+        coer = eigs[0]
 
     plt.rcParams.update(RC_SETTINGS)
 
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    ax.hist(vees, bins=np.linspace(0, 1, 30))
+    ax.grid(zorder=0)
+
+    props = dict(boxstyle="round", facecolor="white", alpha=0.8, zorder=20)
+    ax.text(
+        0.60,
+        0.96,
+        f"Min Escape = {min(vees):.6f}",
+        transform=ax.transAxes,
+        verticalalignment="top",
+        bbox=props,
+    )
+
+    ax.text(
+        0.62,
+        0.89,
+        f"Coercivity = {coer:.6f}",
+        transform=ax.transAxes,
+        verticalalignment="top",
+        bbox=props,
+    )
+
+    out = OUTPUT_DIR / f"{NAME2} - {args.sim_store_path}.png"
+    fig.savefig(out)
+    print(f"Wrote to {out}")
+
+    return
+
     fig = plt.figure(figsize=(30, 8))
     gs = fig.add_gridspec(1, 1)
     ax = fig.add_subplot(gs[0])

scripts/torus_mesh.py

@@ -117,11 +117,6 @@ def torus_region(c, verts, height):
     v = np.hstack((v, height * np.ones((len(v), 1))))
     return v, f, v_inds[sd - 1]
 
-    centers = np.array([[c[0], c[1], heights[i]]])
-    v_top = np.hstack((verts, heights[i] * np.ones((len(verts), 1))))
-
-    face_verts = np.concatenate((centers, v_top))
-
 
 def main():
     parser = argparse.ArgumentParser(description="Generates 3D model for equilibrium.")
@@ -136,7 +131,7 @@ def main():
     args = parser.parse_args()
 
     torus = args.torus
-    size = args.size
+    size = args.size * 10
 
     # Get desired frame and load.
     sim, frames = Simulation.load(args.sims_path)

scripts/vor_test.py

@@ -0,0 +1,43 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.spatial import Voronoi, voronoi_plot_2d
+from squish import ordered, DomainParams
+from script_tools import RC_SETTINGS
+from pathlib import Path
+
+N = 64
+C = (0, 0)
+
+out_fol = Path(f"N{N}C{C}")
+out_fol.mkdir(exist_ok=True)
+
+
+def render(i, dom, vor):
+    plt.rcParams.update(RC_SETTINGS)
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+    voronoi_plot_2d(vor, ax=ax, show_points=True, show_vertices=False)
+    ax.set_xlim(0, dom.w)
+    ax.set_ylim(0, dom.h)
+    fig.savefig(out_fol / f"{i:03}.png")
+
+
+def get_full_points(dom, points):
+    SYMM = np.array(
+        [[0, 0], [1, 0], [1, 1], [0, 1], [-1, 1], [-1, 0], [-1, -1], [0, -1], [1, -1]]
+    )
+    return np.concatenate([points + dom.dim * s for s in SYMM])
+
+
+def main():
+    d = DomainParams(N, 8, 8, 4)
+    # init = get_full_points(d, ordered.sites(d, C))
+    init = get_full_points(d, np.random.random_sample((64, 2)) * d.dim)
+    for i in range(10):
+        vor = Voronoi(init if i == 0 else vor.vertices)
+        render(i, d, vor)
+
+
+if __name__ == "__main__":
+    main()