Changed version number

requirements.txt

@@ -19,7 +19,7 @@ keyring==23.2.1
 kiwisolver==1.3.2
 MarkupSafe==2.0.1
 matplotlib==3.4.3
-numpy==1.21.2
+numpy==1.22.1
 packaging==21.0
 pep517==0.11.0
 Pillow==8.3.2

scripts/continuation.py

@@ -0,0 +1,230 @@
+import argparse, numpy as np, os
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
+from multiprocessing import Pool, cpu_count
+from pathlib import Path
+
+from squish import Simulation, DomainParams, Energy, ordered, Diagram
+from squish.common import OUTPUT_DIR
+from script_tools import (
+    RC_SETTINGS,
+    get_args,
+    get_data,
+    get_simulation_data,
+    get_ordered_data,
+)
+
+NAME = "Continuation"
+
+
+def proc(combo):
+    dia, min_order, min_disorder, hex_ratios, alphas, sim_alphas, sim_energies, offset, length, num_frames = (
+        combo
+    )
+    sim = dia.sim
+    out_folder = sim.path / "continuation"
+
+    plt.rcParams.update(RC_SETTINGS)
+    plt.rcParams.update({"figure.constrained_layout.use": False})
+    fig = plt.figure(figsize=(30, 15))
+
+    e_hex = ordered.e_hex(sim.domains[0])
+
+    for i in range(length):
+        gs = fig.add_gridspec(1, 2)
+        gs.update(left=0, right=0.98, top=0.93, bottom=0.12, wspace=0.08)
+        ax = fig.add_subplot(gs[1])
+
+        curr_alpha = sim.domains[i].w / sim.domains[i].h
+        curr_vee = 100 * (sim.energies[i] / sim.domains[i].n - e_hex)
+
+        ax.plot(
+            alphas,
+            100 * min_order,
+            color="C1",
+            label="Min Ordered",
+            zorder=10,
+            linestyle="dotted",
+        )
+        ax.plot(
+            alphas,
+            100 * min_disorder,
+            color="C0",
+            label="Min Disordered",
+            zorder=11,
+            linestyle="dotted",
+        )
+        ax.plot(
+            sim_alphas[:175],
+            100 * sim_energies[:175],
+            color="C2",
+            label="_nolegend_",
+            linestyle="dashed",
+        )
+        ax.plot(
+            sim_alphas[174:], 100 * sim_energies[174:], color="C2", label="Continuation"
+        )
+
+        ax.scatter(
+            hex_ratios, [0] * len(hex_ratios), color="C2", s=120, marker="H", zorder=50
+        )
+
+        ax.scatter(
+            curr_alpha,
+            curr_vee,
+            s=250,
+            facecolors="none",
+            edgecolors="C6",
+            linewidth=4,
+            zorder=100,
+        )
+
+        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)
+
+        space = np.arange(0, 10, 0.5)
+        ax.set_ylim(-0.15, 9.5)
+        ax.set_yticks(space[::-1])
+        ax.ticklabel_format(axis="y", style="sci")
+
+        ax.set_xlabel("Aspect Ratio")
+        ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$")
+
+        ax.legend(loc="upper center")
+        # ax.legend()
+        ax.grid(zorder=0)
+
+        props = dict(boxstyle="round", facecolor="white", alpha=0.8, zorder=20)
+        ax.text(
+            0.873,
+            0.97,
+            f"N={sim.domains[i].n}",
+            transform=ax.transAxes,
+            verticalalignment="top",
+            bbox=props,
+        )
+
+        ax1 = fig.add_subplot(gs[0])
+        dia.voronoi_plot(i, ax1)
+        fig.savefig(out_folder / f"img{i+offset:05}.png")
+        fig.clear()
+
+        j = len(list((sim.path / "continuation").iterdir()))
+        hashes = int(21 * j / num_frames)
+        print(
+            f'Generating frames... |{"#"*hashes}{" "*(20-hashes)}|'
+            + f" {j}/{num_frames} frames rendered.",
+            flush=True,
+            end="\r",
+        )
+
+    print(flush=True)
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Makes video of analytic continuation."
+    )
+    parser.add_argument(
+        "sims_path",
+        metavar="sim_dir",
+        help="folder that contains simulation data at various aspect ratios for some N",
+    )
+    parser.add_argument(
+        "shrink_sim",
+        metavar="continuation_sim",
+        help="simulation that contains the continuation data",
+    )
+    args = parser.parse_args()
+    sims_path = Path(args.sims_path)
+
+    data, n, r = get_data(
+        sims_path / "package.pkl", get_simulation_data, args=(sims_path,)
+    )
+    domain, alphas = DomainParams(n, 1, 1, r), data["all"]["alpha"]
+    ordered_data = get_data(
+        OUTPUT_DIR / "OrderedCache" / f"{n}.pkl",
+        get_ordered_data,
+        args=(domain, alphas),
+    )
+
+    min_disorder, max_disorder = [], []
+    for i, energies in enumerate(data["distinct"]["Energy"]):
+        disorder_energies = []
+        for j, energy in enumerate(energies):
+            if not data["distinct"]["Ordered"][i][j]:
+                disorder_energies.append(energy)
+        min_disorder.append(min(disorder_energies))
+        max_disorder.append(max(disorder_energies))
+
+    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
+    min_order = np.array(min_order) / domain.n - e_hex
+
+    hex_ratios = ordered.hexagon_alpha(domain.n)
+
+    sim = Simulation.from_file(Path(args.shrink_sim))
+    sim_alphas = np.array([x.w / x.h for x in sim.frames])
+    sim_energies = np.array([x.energy for x in sim.frames]) / sim.domain.n - e_hex
+    dia = Diagram(sim, ["voronoi"])
+
+    out_folder = sim.path / "continuation"
+    out_folder.mkdir(exist_ok=True)
+
+    combo_list = []
+    for i in range(cpu_count()):
+        start, end = (
+            int(i * len(sim) / cpu_count()),
+            int((i + 1) * len(sim) / cpu_count()),
+        )
+        new_dia = Diagram(None, ["voronoi"])
+        new_dia.sim = dia.sim.slice(list(range(start, end)))
+        combo_list.append(
+            (
+                new_dia,
+                min_order,
+                min_disorder,
+                hex_ratios,
+                alphas,
+                sim_alphas,
+                sim_energies,
+                start,
+                len(sim.frames[start:end]),
+                len(sim),
+            )
+        )
+
+    print("Starting figure generation...", flush=True)
+    with Pool(cpu_count()) as pool:
+        for _ in pool.imap_unordered(proc, combo_list):
+            pass
+
+    video_path = sim.path / (NAME + ".mp4")
+
+    fps = 30
+    print("Assembling MP4...", flush=True)
+    os.system(
+        f"ffmpeg -hide_banner -loglevel error -r {fps} -i"
+        + f' "{sim.path}/continuation/img%05d.png"'
+        + f" -c:v libx264 -crf 18 -preset slow -pix_fmt yuv420p -vf"
+        + f' "scale=trunc(iw/2)*2:trunc(ih/2)*2" -f mp4 "{video_path}"'
+    )
+
+    print(f"Wrote to {video_path}.")
+
+
+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_cell_vee.py

@@ -0,0 +1,78 @@
+import numpy as np, os, csv
+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,
+    get_simulation_data,
+    get_ordered_data,
+)
+
+NAME = "Cumulative-CellVEE"
+ALPHA = 1.0
+
+
+def main():
+    sims_path, regen = get_args(
+        "Anti-cumulative distribution of VEE and percent of equilibria for fixed alpha",
+        "folders that contains various N simulations to plot",
+    )
+
+    sim, frames = Simulation.load(sims_path)
+    vees = np.linspace(-3, 3, 10000)
+    e_hex = ordered.e_hex(sim.domain)
+    energies = {"all": np.empty(0, dtype=float)}
+    counts = {}
+    for frame in frames:
+        energy = frame["stats"]["site_energies"] - e_hex
+        defects = np.count_nonzero(frame["stats"]["site_edge_count"] != 6)
+        if defects not in energies:
+            energies[defects] = np.empty(0, dtype=float)
+        energies[defects] = np.append(energies[defects], energy)
+        energies["all"] = np.append(energies["all"], energy)
+
+    all_count = None
+    for defect, energy in energies.items():
+        count = np.empty(vees.shape, dtype=float)
+        for i, vee in enumerate(vees):
+            count[i] = np.count_nonzero(energy >= vee)
+        count = 100 * count / len(energy)
+        if defect == "all":
+            all_count = count
+        else:
+            counts[defect] = count
+
+    plt.rcParams.update(RC_SETTINGS)
+
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    for defect, count in sorted(counts.items()):
+        if defect == min(counts.keys()):
+            ax.plot(vees, count, label=fr"$\mathrm{{D}}={defect}$", color="C0")
+        elif defect == max(counts.keys()):
+            ax.plot(vees, count, label=fr"$\mathrm{{D}}={defect}$", color="C1")
+        else:
+            ax.plot(vees, count, linestyle="dotted", alpha=0.3)
+    ax.plot(vees, all_count, label="All", color="black")
+
+    ax.set_xlim(-2.5, 2)
+    ax.yaxis.set_major_formatter(mtick.PercentFormatter())
+
+    ax.set_xlabel(r"VEE")
+    ax.set_ylabel("Percent of Voronoi Regions")
+
+    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/cumulative_vee.py

@@ -0,0 +1,108 @@
+import numpy as np, os, csv
+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,
+    get_simulation_data,
+    get_ordered_data,
+)
+
+NAME = "Cumulative-VEE"
+ALPHA = 1.0
+
+
+def main():
+    sims_path, regen = get_args(
+        "Anti-cumulative distribution of VEE and percent of equilibria for fixed alpha",
+        "folders that contains various N simulations to plot",
+    )
+
+    packages = []
+    for fol in sims_path.iterdir():
+        if fol.is_file():
+            continue
+
+        data, n, r = get_data(
+            fol / "package.pkl", get_simulation_data, args=(fol,), regen=regen
+        )
+        domain, alphas = DomainParams(n, 1, 1, r), data["all"]["alpha"]
+        ordered_data = get_data(
+            OUTPUT_DIR / "OrderedCache" / f"{n}.pkl",
+            get_ordered_data,
+            args=(domain, alphas),
+            regen=regen,
+        )
+
+        packages.append([data, ordered_data, domain])
+
+    packages.sort(key=lambda x: x[2].n)
+
+    plt.rcParams.update(RC_SETTINGS)
+
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    my_cool_data = [["VEE", 61, 67, 73, 81, 84, 100]]
+    for vee in np.linspace(0, 0.06, 10000):
+        my_cool_data.append([vee])
+
+    for j, package in enumerate(packages):
+        data, ordered_data, domain = package
+        e_hex = ordered.e_hex(domain)
+
+        alpha_index = np.where(data["all"]["alpha"] == ALPHA)[0][0]
+
+        energies = data["all"]["Energy"][alpha_index] / domain.n - e_hex
+
+        min_order = ordered_data["Energy"][alpha_index][0] / domain.n - e_hex
+
+        vees = np.linspace(0, 0.06, 10000)
+        index = np.argmin(np.abs(vees - min_order))
+
+        counts = np.empty(vees.shape, dtype=float)
+        for i, vee in enumerate(vees):
+            counts[i] = np.count_nonzero(energies >= vee)
+        counts = 100 * counts / len(energies)
+
+        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(
+        #    100 * vees[: index + 1],
+        #    counts[: index + 1],
+        #    label=f"N={domain.n}",
+        #    color=f"C{j}",
+        #)
+        #ax.plot(
+        #    100 * vees[index:],
+        #    counts[index:],
+        #    label=f"_nolegend_",
+        #    linestyle="dotted",
+        #    color=f"C{j}",
+        #)
+
+    with open("cumulative-vee.csv", "w") as csvfile:
+        writer = csv.writer(csvfile)
+        writer.writerows(my_cool_data)
+
+    ax.set_xlim(0, 6.3)
+    ax.yaxis.set_major_formatter(mtick.PercentFormatter())
+
+    ax.set_xlabel(r"VEE $\left[\times 10^{2}\right]$")
+    ax.set_ylabel("Percent of Equilibria")
+
+    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/cumulative_vee_2.py

@@ -0,0 +1,215 @@
+import numpy as np, os, csv
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
+from scipy.optimize import curve_fit
+
+from squish import Simulation, DomainParams, ordered
+from squish.common import OUTPUT_DIR
+from script_tools import (
+    RC_SETTINGS,
+    get_args,
+    get_data,
+    get_simulation_data,
+    get_ordered_data,
+)
+
+NAME = "Cumulative-VEE-AcrossN"
+NAME2 = "Cumulative-VEE-RelError"
+NAME3 = "Cumulative-VEE-Alphas"
+
+
+def main():
+    sims_path, regen = get_args(
+        "Anti-cumulative distribution of VEE and percent of equilibria for fixed alpha",
+        "folders that contains various N simulations to plot",
+    )
+
+    vees = np.linspace(0, 0.06, 10000)
+
+    def f(path):
+        files = list(path.iterdir())
+        files = [f for f in files if f.is_dir()]
+        data = {}
+
+        for i, fol in enumerate(files):
+            sim, frames = Simulation.load(fol)
+            e_hex = ordered.e_hex(sim.domain)
+            energies = []
+            for frame in frames:
+                energies.append(frame["energy"] / sim.domain.n - e_hex)
+            energies = np.array(energies)
+
+            counts = np.empty(vees.shape, dtype=float)
+            for j, vee in enumerate(vees):
+                counts[j] = np.count_nonzero(energies >= vee)
+            counts = 100 * counts / len(energies)
+
+            data[sim.domain.n] = counts
+
+            hashes = int(20 * i / len(files))
+            print(
+                f'Loading simulations... |{"#"*hashes}{" "*(20-hashes)}|'
+                + f" {i+1}/{len(files)} simulations loaded.",
+                flush=True,
+                end="\r",
+            )
+
+        print(flush=True)
+        return data
+
+    all_data = {}
+    alpha_files = list(sims_path.iterdir())
+    alpha_files = [f for f in alpha_files if f.is_dir()]
+    for fol in alpha_files:
+        all_data[float(fol.name)] = get_data(
+            fol / "Cumulative2.pkl", f, args=(fol,), regen=regen
+        )
+
+    data = all_data[1.0]
+    points = []
+    avgs = {}
+    # print((data[106] + data[104]) / 2 - data[105])
+    for n in range(105, 396):
+        arr = np.zeros(data[100].shape)
+        for i in range(1, 6):
+            arr += data[n - i] + data[n + i]
+        arr += data[n]
+        arr /= 11
+        avgs[n] = arr
+        points.append(np.linalg.norm(data[n] - arr) / np.linalg.norm(arr))
+
+    plt.rcParams.update(RC_SETTINGS)
+
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    ax.plot(100 * vees, avgs[395], color="black", label=r"Average", zorder=50)
+
+    for n, counts in sorted(data.items()):
+        if n in range(390, 401):
+            if n == 392:
+                ax.plot(100 * vees, counts, label="N=392", color="C0")
+            elif n == 397:
+                ax.plot(100 * vees, counts, label="N=397", color="C1")
+            elif n == 395:
+                ax.plot(100 * vees, counts, color="black", linestyle="dashed")
+            else:
+                ax.plot(
+                    100 * vees,
+                    counts,
+                    label="_nolegend_",
+                    linestyle="dashed",
+                    alpha=0.5,
+                )
+
+    ax.set_xlim(0, 6.3)
+    ax.yaxis.set_major_formatter(mtick.PercentFormatter())
+
+    ax.set_xlabel(r"VEE $\left[\times 10^{2}\right]$")
+    ax.set_ylabel("Percent of Equilibria")
+
+    ax.grid(zorder=0)
+    ax.legend()
+
+    fig.savefig(OUTPUT_DIR / (NAME + ".png"))
+    print(f"Wrote to {OUTPUT_DIR / (NAME + '.png')}")
+
+    fig = plt.figure(figsize=(30, 8))
+    gs = fig.add_gridspec(1, 1)
+
+    ax = fig.add_subplot(gs[0])
+    ns = np.array(range(105, 396))
+    ax.scatter(ns, points)
+
+    log_slope, _ = np.polyfit(np.log10(ns), np.log10(points), 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, points, p0=[5000])
+    ax.plot(np.arange(50, 500), g(np.arange(50, 500), *params), color="C1")
+
+    # with open("cumul.csv", "w") as csvfile:
+    #    csvwriter = csv.writer(csvfile)
+    #    csvwriter.writerows(zip(ns, points))
+    print(params)
+    ax.set_xlim(95, 405)
+    ax.set_ylim(0, 0.24)
+    ax.set_xlabel("N")
+
+    ax.grid(zorder=0)
+    # ax.legend()
+
+    fig.savefig(OUTPUT_DIR / (NAME2 + ".png"))
+    print(f"Wrote to {OUTPUT_DIR / (NAME2 + '.png')}")
+
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    alpha_avgs = {}
+    for alpha, data in all_data.items():
+        if len(data.values()) == 0 or alpha == 1.0:
+            continue
+        alpha_avgs[alpha] = sum(data.values()) / 11
+
+    for alpha, data in all_data.items():
+        if alpha == 1.0:
+            continue
+        for n, counts in data.items():
+            # ax.plot(100 * vees, counts, linestyle="dotted", alpha=0.4)
+            continue
+
+    pad = 100
+    alpha_prob_dens = {}
+    for alpha, avg in alpha_avgs.items():
+        mov_avgs = np.array(
+            [np.mean(avg[max(0, i - pad) : i + pad]) for i in range(len(avg))]
+        )
+        d_avgs = np.gradient(mov_avgs, vees[1] - vees[0])
+        alpha_prob_dens[alpha] = d_avgs / (np.sum(d_avgs) * (vees[1] - vees[0]))
+
+    a1_prob_dens = {}
+    for n, avg in avgs.items():
+        mov_avgs = np.array(
+            [np.mean(avg[max(0, i - pad) : i + pad]) for i in range(len(avg))]
+        )
+
+        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
+            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")
+    
+    # Probability densities of alpha=1 across N
+
+    #for n, prob_dens in sorted(a1_prob_dens.items()):
+    #    if n in [255, 315, 335, 355, 375]:
+    #        ax.plot(100 * vees, prob_dens, label=f"N={n}")
+    #    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)
+    # ax.yaxis.set_major_formatter(mtick.PercentFormatter())
+
+    ax.set_xlabel(r"VEE $\left[\times 10^{2}\right]$")
+    # ax.set_ylabel("Percent of Equilibria")
+
+    ax.grid(zorder=0)
+    ax.legend()
+    fig.savefig(OUTPUT_DIR / (NAME3 + ".png"))
+    print(f"Wrote to {OUTPUT_DIR / (NAME3 + '.png')}")
+
+
+if __name__ == "__main__":
+    main()

scripts/defect_scatter.py

@@ -0,0 +1,76 @@
+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)
+
+    counts = {}
+    for i in range(70):
+        counts[i] = np.count_nonzero(defects == i)
+    print(counts)
+    for k, v in counts.items():
+        if k % 2 == 1 and v > 0:
+            print(k, v)
+
+    plt.rcParams.update(RC_SETTINGS)
+    fig = plt.figure(figsize=(15, 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=70, zorder=2)
+    ax.plot(np.arange(0, 65), np.arange(0, 65) * 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.7)
+    ax.set_yticks(np.arange(0, 3.6, 0.5))
+
+    ax.set_xlabel("$\mathrm{D}$")
+    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.")

scripts/density_of_states.py

@@ -0,0 +1,63 @@
+from __future__ import annotations
+from typing import List, Tuple, Dict
+import argparse, math, numpy as np, os, pickle, itertools
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
+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():
+    sims_path, regen = get_args(
+        "Density of states of various N across alpha",
+        "folders that contains various N simulations to plot",
+    )
+
+    packages = []
+    for fol in sims_path.iterdir():
+        if fol.is_file():
+            continue
+
+        data, n, r = get_data(
+            fol / "package.pkl", get_simulation_data, args=(fol,), regen=regen
+        )
+        domain = DomainParams(n, 1, 1, r)
+
+        packages.append([data, domain])
+
+    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])
+
+    for j, package in enumerate(packages):
+        data, domain = package
+
+        distinct_disordered = []
+        for ordered in data["distinct"]["Ordered"]:
+            distinct_disordered.append(np.count_nonzero(ordered == False))
+
+        ax.plot(data["all"]["alpha"], distinct_disordered, label=f"N={domain.n}")
+
+    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")
+
+    ax.grid(zorder=0)
+    ax.legend(loc="center right", fancybox=True, bbox_to_anchor=(1.34, 0.5))
+    fig.savefig(OUTPUT_DIR / "DoS.png")
+
+    print(f"Wrote to {OUTPUT_DIR / 'DoS.png'}")
+
+
+if __name__ == "__main__":
+    main()

scripts/frustration_alpha.py

@@ -0,0 +1,185 @@
+import numpy as np, os, math
+import matplotlib.pyplot as plt
+from multiprocessing import Pool, cpu_count
+
+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, 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)
+
+    avg_defects = sum(defects) / len(defects)
+    avg_vee = sum(energy) / len(energy)
+    m, b = np.polyfit(defects, energy, 1)
+
+    alpha = domain.w / domain.h
+
+    return [alpha, tuple(domain), avg_defects, avg_vee, m, b]
+
+
+def main():
+    sims_path, regen = get_args(
+        "Generates graph for Frustration at fixed N across alpha",
+        "folder that contains simulations at various alpha for fixed N",
+    )
+
+    def f():
+        data = {}
+
+        files = list(sims_path.iterdir())
+        files = [x for x in files if x.is_dir()]
+        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 alpha={res[0]:.3f} |{"#"*hashes}{" "*(20-hashes)}|'
+                    + f" {i+1}/{len(files)} simulations processed.",
+                    flush=True,
+                    end="\r",
+                )
+
+            print(flush=True)
+
+        return format_data(
+            data,
+            key_name="alpha",
+            col_names=[
+                "Domain",
+                "Average Defects",
+                "Average VEE",
+                "Energy Per Defect",
+                "Ground VEE",
+            ],
+        )
+
+    plt.rcParams.update(RC_SETTINGS)
+    data = get_data(sims_path / (NAME + ".pkl"), f, regen=regen)
+    domain = DomainParams(*data["Domain"][0])
+    # ordered_data = get_data(
+    #     OUTPUT_DIR / "OrderedCache" / f"{domain.n}.pkl",
+    #     get_ordered_data,
+    #     args=(domain, data["alpha"]),
+    #     regen=regen,
+    # )
+    # min_order = []
+    # for i, energies in enumerate(ordered_data["Energy"]):
+    # min_order.append(energies[0])
+
+    alphas, avg_defects, avg_vees, epds, vee0s = (
+        data["alpha"],
+        data["Average Defects"],
+        100 * data["Average VEE"],
+        100 * data["Energy Per Defect"],
+        100 * data["Ground VEE"],
+        # 100 * (np.array(min_order) / domain.n - ordered.e_hex(domain)),
+    )
+
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    m0, b0 = np.polyfit(alphas, avg_defects, 1)
+    # m1, b1 = np.polyfit(ns, avg_defects * epds, 1)
+
+    lns1 = ax.plot(
+        alphas,
+        avg_defects,
+        color="C0",
+        alpha=0.5,
+        label=r"$\langle \mathrm{D} \rangle$",
+    )
+    ax.plot(alphas, m0 * alphas + b0, color="C0", linestyle="dashed")
+
+    lns2 = ax.plot(
+        alphas, 100 * epds, color="C1", alpha=0.5, label=r"$\zeta_1 \times 10^{4}$"
+    )
+
+    lns3 = ax.plot(
+        alphas,
+        10 * avg_defects * epds,
+        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_xlim(0.3, 1.0)
+
+    # 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")
+
+    out = OUTPUT_DIR / (NAME + f" - N{domain.n}.png")
+    fig.savefig(out)
+    print(f"Wrote to {out}")
+
+    # return
+    fig = plt.figure(figsize=(20, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    ax.plot(alphas, vee0s, label="$\zeta_0$", color="C5")
+    ax.plot(alphas, avg_defects * epds, label=r"$\mathcal{F}$", color="C2")
+    ax.plot(alphas, avg_vees, label=r"$\langle \mathrm{VEE} \rangle$", color="C4")
+    # ax.plot(alphas, min_order, label="Minimum Ordered", color="C7", alpha=0.5)
+
+    print(np.linalg.norm(avg_vees - (vee0s + avg_defects * epds)))
+
+    ax.grid(zorder=0)
+    ax.legend()
+
+    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)
+
+    # ax.set_ylim(-0.1, 9.2)
+    # ax.set_yticks(np.arange(0, 9.6, 1))
+
+    ax.set_xlabel("Aspect Ratio")
+    ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$")
+
+    out = OUTPUT_DIR / (NAME2 + f" - N{domain.n}.png")
+    fig.savefig(out)
+    print(f"Wrote to {out}")
+
+
+if __name__ == "__main__":
+    os.environ["QT_log10GING_RULES"] = "*=false"
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Program terminated by user.")

scripts/frustration_n.py

@@ -0,0 +1,176 @@
+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=(15, 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_1 \times 10^{4}$"
+    #)
+    #
+    #lns3 = ax.plot(
+    #    ns,
+    #    10 * avg_defects * epds,
+    #    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_xlim(93, 407)
+
+    ax.set_ylim(3, 36)
+    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=(20, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    ax.plot(ns, vee0s, label="$\zeta_0$", color="C5")
+    ax.plot(ns, avg_defects * epds, label=r"$\mathcal{F}$", color="C2")
+    ax.plot(ns, avg_vees, label=r"$\langle \mathrm{VEE} \rangle$", color="C4")
+    ax.scatter(ns, min_order, label="Minimum Ordered", color="C7", alpha=0.8)
+
+    print(np.linalg.norm(avg_vees - (vee0s + avg_defects * epds)))
+
+    ax.grid(zorder=0)
+    ax.legend()
+
+    ax.set_xlim(95, 405)
+
+    ax.set_ylim(-0.1, 9.2)
+    ax.set_yticks(np.arange(0, 9.6, 1))
+
+    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.")

scripts/hexagon_alpha.py

@@ -0,0 +1,89 @@
+import argparse, numpy as np, os
+from multiprocessing import Pool, cpu_count
+import matplotlib.pyplot as plt
+
+from squish import ordered
+from squish.common import OUTPUT_DIR
+from script_tools import RC_SETTINGS, get_data, format_data
+
+NAME = "RegHexTessRatios"
+
+
+def get_ratios(n: int):
+    return (n, ordered.hexagon_alpha(n))
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Generates graph for regular hexagonal tessellation ratios"
+    )
+    parser.add_argument(
+        "max_n", metavar="N", type=int, help="maximum N  of which to calculate"
+    )
+
+    parser.add_argument(
+        "--regenerate",
+        dest="regen",
+        action="store_true",
+        help="regenerates the cache file for processed data",
+    )
+
+    args = parser.parse_args()
+
+    def f():
+        data = {"alpha": [], "N": []}
+
+        with Pool(cpu_count()) as pool:
+            for i, res in enumerate(
+                pool.imap_unordered(get_ratios, range(2, args.max_n + 1))
+            ):
+                for ratio in res[1]:
+                    data["N"].append(res[0])
+                    data["alpha"].append(ratio)
+
+                hashes = int(21 * i / (args.max_n - 1))
+                print(
+                    f'Processed N={res[0]} |{"#"*hashes}{" "*(20-hashes)}|'
+                    + f" {i+1}/{args.max_n-1}",
+                    flush=True,
+                    end="\r",
+                )
+
+            print(flush=True)
+        return data
+
+    data = get_data(OUTPUT_DIR / "OrderedCache" / (NAME + ".pkl"), f, 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])
+
+    alphas, ns = [], []
+    for alpha, n in zip(data["alpha"], data["N"]):
+        if n <= args.max_n:
+            alphas.append(alpha)
+            ns.append(n)
+
+    ax.scatter(alphas, ns, s=5)
+
+    ax.set_xlim(-0.01, 1.01)
+    ax.set_xticks(np.arange(0, 1.01, 0.1))
+    ax.set_ylim(0, args.max_n)
+
+    ax.set_ylabel(r"$N$")
+    ax.set_xlabel("Aspect Ratio")
+    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.")

scripts/min_order_vee.py

@@ -0,0 +1,107 @@
+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
+from squish.common import OUTPUT_DIR
+from script_tools import RC_SETTINGS, get_args, get_data, format_data
+
+NAME = "MinOrderVEE"
+
+
+def get_min_vee(n):
+    domain = DomainParams(n, math.sqrt(n), math.sqrt(n), 4)
+    e_hex = ordered.e_hex(domain)
+
+    configs = ordered.configurations(domain)
+    min_vee, min_config = 1e10, None
+
+    for config in configs:
+        try:
+            rbar = ordered.avg_radius(domain, config)
+        except:
+            print(tuple(domain), config)
+
+        vee = (
+            2 * domain.w * domain.h
+            + 2 * math.pi * domain.n * (domain.r ** 2 - 2 * domain.r * rbar)
+        ) / domain.n - e_hex
+
+        if vee < min_vee and vee >= 0:
+            min_vee = vee
+            min_config = config
+
+    return (domain.n, min_vee, min_config)
+
+
+def main():
+    def f():
+        data = {}
+        ns = list(range(1000, 5001))
+        with Pool(cpu_count()) as pool:
+            for i, res in enumerate(pool.imap_unordered(get_min_vee, ns)):
+                data[res[0]] = res[1:]
+
+                hashes = int(21 * i / len(ns))
+                print(
+                    f'Processed N={res[0]} |{"#"*hashes}{" "*(20-hashes)}|'
+                    + f" {i+1}/{len(ns)} simulations processed.",
+                    flush=True,
+                    end="\r",
+                )
+
+            print(flush=True)
+
+        return format_data(
+            data, key_name="N", col_names=["Minimum Ordered VEE", "Config"]
+        )
+
+    plt.rcParams.update(RC_SETTINGS)
+    data = get_data(OUTPUT_DIR / "OrderedCache" / (NAME + ".pkl"), f)
+    ns, min_order, min_config = (
+        data["N"],
+        100 * data["Minimum Ordered VEE"],
+        data["Config"],
+    )
+
+    fig = plt.figure(figsize=(20, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    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(800, 5200)
+
+    ax.set_xlabel("N")
+    ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$")
+
+    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.")

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/ordered_rbar.py

@@ -0,0 +1,72 @@
+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 = "OrderedRBar"
+
+
+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=(18, 14.4))
+    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_ylim(0.51, 0.615)
+    ax.set_yticks(np.arange(0.52, 0.62, 0.03))
+    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.")

scripts/ordered_scatter.py

@@ -0,0 +1,79 @@
+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", s=25
+        )
+
+    hex_ratios = ordered.hexagon_alpha(args.n)
+    ax.scatter(
+        hex_ratios, [0] * len(hex_ratios), color="C2", s=200, marker="H", zorder=10
+    )
+
+    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.")

scripts/perturbations.py

@@ -1,73 +1,179 @@
-from __future__ import annotations
-from typing import List
-import argparse, pickle, numpy as np, os
+import argparse, numpy as np, os
 from pathlib import Path
 import matplotlib.pyplot as plt
 
-from squish import Simulation
-from squish.common import OUTPUT_DIR
+from squish import ordered
+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():
     parser = argparse.ArgumentParser(
-        "Graphs perturbation graphs for a collection of simulations."
+        description="Graphs perturbation graphs for a collection of simulations."
     )
     parser.add_argument(
-        "sims_path",
-        metavar="path/to/data",
-        help="folder that contains simulations of perturbations from an equilibrium.",
+        "sim_store_path", metavar="sim_dir", help="folder to save simulations to"
     )
     parser.add_argument(
         "end_path",
-        metavar="path/to/equilbrium",
-        help="NumPy binary (.npy) file that contains the equilibrium to compare to.",
+        metavar="eq_path",
+        help="simulation that contains the equilibrium to compare to.",
     )
+
     parser.add_argument(
-        "-q",
-        "--quiet",
-        dest="quiet",
+        "--regenerate",
+        dest="regen",
         action="store_true",
-        default=False,
-        help="suppress all normal output",
+        help="regenerates the cache file for processed data",
     )
 
     args = parser.parse_args()
+    out_fol = Path(args.sim_store_path)
 
-    end = np.load(args.end_path)
+    (OUTPUT_DIR / out_fol).mkdir(exist_ok=True)
 
-    data = {}
-    for file in Path(args.sims_path).iterdir():
-        k = float(file.name.split("k")[-1])
-        delta = 10 ** k
+    end_eq = Simulation.from_file(args.end_path)
+    e_hex = ordered.e_hex(end_eq.domain)
 
-        sim, frames = Simulation.load(file)
-        data[delta] = {"norm": [], "time": [], "k": k}
+    def f():
+        all_data = []
+        for j in range(50):
+            pert_out_fol = out_fol / f"Vector{j:03}"
+            (OUTPUT_DIR / pert_out_fol).mkdir(exist_ok=True)
 
-        for i, frame in enumerate(frames):
-            adjusted = frame["arr"] + (end[0] - frame["arr"][0])
+            perturb = np.random.random_sample(end_eq.frames[0].site_arr.shape)
+            perturb /= np.linalg.norm(perturb)
 
-            data[delta]["norm"].append(np.linalg.norm(adjusted - end))
-            data[delta]["time"].append(sim.step_size * i)
+            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)
 
-    fig, ax = plt.subplots(figsize=(12, 8))
-    plt.subplots_adjust(0.07, 0.12, 0.97, 0.9)
+                sim, frames = Simulation.load(OUTPUT_DIR / k_out_fol)
+
+                for i, frame in enumerate(frames):
+                    adjusted = frame["arr"] + (
+                        end_eq.frames[0].site_arr[0] - frame["arr"][0]
+                    )
+
+                    data[delta]["norm"].append(
+                        np.linalg.norm(
+                            ordered.toroidal_distance(
+                                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)
+
+                k += 1 if k < 0 else 0.25
+                m, _ = np.polyfit(data[delta]["time"], data[delta]["norm"], 1)
+                same_eq = m < 0
+
+            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])
 
     for delta in sorted(data):
         ax.plot(
-            np.log10(np.array(data[delta]["time"]) + 1),
-            np.log10(data[delta]["norm"]),
+            np.array(data[delta]["time"]),
+            np.array(data[delta]["norm"]),
             label=f"k = {data[delta]['k']}",
         )
 
-    fig.suptitle("Equilibrium Perturbations")
-    ax.grid(zorder=0)
-    # ax.set_xlim([0, 5])
-    ax.legend()
-    ax.set_xlabel("$\\log_{10}$ Time")
-    ax.set_ylabel("$\\log_{10}$ L2 Norm of Difference")
+    # ax.set_title(r"Relaxation of Perturbations")
 
-    fig.savefig(OUTPUT_DIR / "Equilibrium Perturbations.png")
-    print(f"Wrote to {OUTPUT_DIR / 'Equilibrium Perturbations.png'}")
+    ax.set_xlim([0, 60])
+    ax.set_yscale("log")
+
+    ax.set_xlabel(r"Time")
+    ax.set_ylabel(r"$\|\mathbf{x}-\mathbf{x_e}\|_2$")
+
+    h, l = ax.get_legend_handles_labels()
+    ax.legend(h[::-1], l[::-1])
+    ax.grid(zorder=0)
+
+    fig.savefig(OUTPUT_DIR / (NAME + ".png"))
+    print(f"Wrote to {OUTPUT_DIR / (NAME + '.png')}")
 
 
 if __name__ == "__main__":

scripts/probability_of_disorder.py

@@ -0,0 +1,118 @@
+import numpy as np, os
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
+
+from squish import Simulation, DomainParams, Energy, ordered
+from squish.common import OUTPUT_DIR
+from script_tools import (
+    RC_SETTINGS,
+    get_args,
+    get_data,
+    get_simulation_data,
+    get_ordered_data,
+)
+
+NAME = "PoD"
+
+
+def main():
+    sims_path, regen = get_args(
+        "Various graphs for single N data.",
+        "folder that contains simulation data at various aspect ratios for some N",
+    )
+
+    data, n, r = get_data(
+        sims_path / "package.pkl", get_simulation_data, args=(sims_path,), regen=regen
+    )
+    domain, alphas = DomainParams(n, 1, 1, r), data["all"]["alpha"]
+    ordered_data = get_data(
+        OUTPUT_DIR / "OrderedCache" / f"{n}.pkl",
+        get_ordered_data,
+        args=(domain, alphas),
+        regen=regen,
+    )
+
+    min_disorder, max_disorder = [], []
+    for i, energies in enumerate(data["distinct"]["Energy"]):
+        disorder_energies = []
+        for j, energy in enumerate(energies):
+            if not data["distinct"]["Ordered"][i][j]:
+                disorder_energies.append(energy)
+        min_disorder.append(min(disorder_energies))
+        max_disorder.append(max(disorder_energies))
+
+    min_order = []
+    for i, energies in enumerate(ordered_data["Energy"]):
+        min_order.append(energies[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
+    min_order = np.array(min_order) / domain.n - e_hex
+
+    all_disorder_count = []
+    for disorders in data["all"]["Ordered"]:
+        all_disorder_count.append(
+            100 * np.count_nonzero(disorders == False) / len(disorders)
+        )
+
+    plt.rcParams.update(RC_SETTINGS)
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+    ax2 = ax.twinx()
+
+    ax.plot(alphas, 100 * (min_order - min_disorder), color="C2")
+    ax2.plot(alphas, all_disorder_count, color="C4")
+
+    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)
+    space = np.arange(-2.3, 3.3, 0.5)
+    ax.set_ylim(-2.5, 3.4)
+    ax.set_yticks(space)
+    ax.ticklabel_format(axis="y", style="sci")
+
+    # start, end = ax2.get_ylim()
+    # space = np.linspace(start, end, 20)
+    # space += 100 - space[-2]
+    space = np.linspace(50, 100, len(space))
+    ax2.set_ylim(
+        space[0] - (space[1] - space[0]) * 0.4, space[-1] + (space[1] - space[0]) * 0.6
+    )
+    ax2.set_yticks(space)
+    # ax2.set_ylim(start, space[-1])
+    # ax2.set_yticks(space[1:-1])
+    ax2.yaxis.set_major_formatter(mtick.PercentFormatter())
+
+    ax.set_xlabel("Aspect Ratio")
+    ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$", color="C2")
+    ax2.set_ylabel("POD", color="C4")
+
+    # ax.legend(loc=(0.23, 0.5))
+    ax.grid(zorder=0)
+
+    props = dict(boxstyle="round", facecolor="white", alpha=0.8, zorder=20)
+    ax.text(
+        0.85,
+        0.92,
+        f"N={domain.n}",
+        transform=ax.transAxes,
+        verticalalignment="top",
+        bbox=props,
+    )
+
+    out = OUTPUT_DIR / f"{NAME}-N{domain.n}.png"
+    fig.savefig(out)
+    print(f"Wrote to {out}.")
+
+
+if __name__ == "__main__":
+    os.environ["QT_LOGGING_RULES"] = "*=false"
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Program terminated by user.")

scripts/rbar.py

@@ -0,0 +1,85 @@
+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_args, get_data, format_data
+
+NAME = "RBar"
+
+
+def main():
+    sims_path, regen = get_args(
+        "Variance of average radius for as N gets larger",
+        "folders that contains various N simulations to plot",
+    )
+
+    def f():
+        files = list(sims_path.iterdir())
+        files = [f for f in files if f.is_dir()]
+        data = {}
+
+        for i, fol in enumerate(files):
+            sim, frames = Simulation.load(fol)
+
+            variances = []
+            for frame in frames:
+                rbars = np.sort(frame["stats"]["avg_radius"])
+                variances.append(np.var(rbars))
+            avg_variance = sum(variances) / len(variances)
+
+            data[sim.domain.n] = [avg_variance]
+
+            hashes = int(20 * i / len(files))
+            print(
+                f'Loading simulations... |{"#"*hashes}{" "*(20-hashes)}|'
+                + f" {i+1}/{len(files)} simulations loaded.",
+                flush=True,
+                end="\r",
+            )
+
+        print(flush=True)
+        return format_data(data, key_name="N", col_names=["Average Variance"])
+
+    data = get_data(sims_path / (NAME + ".pkl"), f, regen=regen)
+
+    plt.rcParams.update(RC_SETTINGS)
+
+    fig = plt.figure(figsize=(18, 14.4))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    primes = [
+        i
+        for i in range(min(data["N"]), max(data["N"]))
+        if len(ordered.divisors(i)) == 2
+    ]
+
+    prime_points = []
+    for i, n in enumerate(data["N"]):
+        if n in primes:
+            prime_points.append(data["Average Variance"][i])
+
+    ax.plot(data["N"], 1e5 * data["Average Variance"])
+    ax.scatter(primes, 1e5 * np.array(prime_points), marker="o", s=60)
+    ax.set_xlim(95, 405)
+    ax.set_ylim(0, 11)
+    # ax.set_ylim(0, args.max_n)
+
+    ax.set_xlabel("N")
+    ax.set_ylabel(r"$\sigma \left[\times 10^5 \right]$")
+    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.")

scripts/script_tools.py

@@ -0,0 +1,217 @@
+from __future__ import annotations
+from typing import List
+import argparse, pickle, numpy as np, math, os
+from pathlib import Path
+from multiprocessing import Pool, cpu_count
+import matplotlib.pyplot as plt
+
+from squish import DomainParams, Simulation, Energy, ordered
+from squish.common import OUTPUT_DIR
+
+RC_SETTINGS = {
+    "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,
+    "font.family": "cm",
+    "font.size": 40,
+    "text.usetex": True,
+    "text.latex.preamble": r"\usepackage{amsmath, amsfonts}",
+    "figure.constrained_layout.use": True,
+}
+
+
+def get_args(
+    script_desc: str, path_desc: str, cache_file: bool = True
+) -> Tuple[Path, bool]:
+    parser = argparse.ArgumentParser(description=script_desc)
+    parser.add_argument("sims_path", metavar="sim_dir", help=path_desc)
+
+    parser.add_argument(
+        "--regenerate",
+        dest="regen",
+        action="store_true",
+        help="regenerates the cache file for processed data",
+    )
+
+    args = parser.parse_args()
+
+    return (Path(args.sims_path), args.regen)
+
+
+def get_data(
+    path: Path, func: Callable[Any, Any], args: Tuple[Any] = (), regen: bool = False
+) -> Any:
+    if regen:
+        try:
+            os.remove(path)
+        except FileNotFoundError:
+            pass
+
+    if path.is_file():
+        with open(path, "rb") as f:
+            return pickle.load(f)
+    else:
+        data = func(*args)
+        with open(path, "wb") as f:
+            pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)
+        return data
+
+
+def format_data(
+    data: Dict[Any, Any], key_name: str, col_names=List[str]
+) -> Dict[Any, Any]:
+    data = sorted(data.items())
+    new_data = {}
+    new_data[key_name] = np.array([x[0] for x in data])
+    for i, col_name in enumerate(col_names):
+        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])
+    return new_data
+
+
+def get_ordered_data(orig_domain: DomainParams, asps: np.ndarray) -> Dict:
+    data = {}
+    domains = []
+    for alpha in asps:
+        domains.append(
+            [
+                DomainParams(
+                    orig_domain.n,
+                    math.sqrt(orig_domain.n * alpha),
+                    math.sqrt(orig_domain.n / alpha),
+                    orig_domain.r,
+                ),
+                alpha,
+            ]
+        )
+
+    with Pool(cpu_count()) as pool:
+        for i, res in enumerate(pool.imap_unordered(ordered_data_proc, domains)):
+            data[res[0]] = res[1:]
+
+            hashes = int(21 * i / len(asps))
+            print(
+                f'Generating at width {res[0]:.02f}... |{"#"*hashes}{" "*(20-hashes)}|'
+                + f" {i+1}/{len(asps)} completed.",
+                flush=True,
+                end="\r",
+            )
+
+        print(flush=True)
+
+    return format_data(data, key_name="alpha", col_names=["Energy", "Coercivity"])
+
+
+def ordered_data_proc(
+    dom_tup: Tuple[DomainParams, float]
+) -> Tuple[float, float, float, float]:
+    domain, alpha = dom_tup
+    energies, coercivities = [], []
+    e_hex = ordered.e_hex(domain)
+
+    configs = ordered.configurations(domain)
+    for config in configs:
+        # Causes errors, so ignore.
+        if config[0] == 0 or config[1] == 0:
+            continue
+        # rbar = ordered.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)
+
+        sites = ordered.sites(domain, config)
+        frame = Energy("radial-t").mode(*domain, sites)
+        energies.append(frame.energy)
+
+        eigs = np.sort(np.linalg.eigvalsh(frame.hessian))
+
+        zero_ind = np.where(np.isclose(eigs, 0, atol=1e-8))[0]
+        if len(zero_ind) == 0:
+            coercivities.append(eigs[0])
+        elif zero_ind[0] == 0:
+            coercivities.append(eigs[2])
+        else:
+            coercivities.append(eigs[0])
+
+    energies, coercivities = list(zip(*sorted(zip(energies, coercivities))))
+    return (alpha, energies, coercivities)
+
+
+def get_simulation_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(simulation_data_proc, 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)
+
+    data["all"] = format_data(
+        data["all"], key_name="alpha", col_names=["Energy", "Ordered", "Defects"]
+    )
+
+    data["distinct"] = format_data(
+        data["distinct"],
+        key_name="alpha",
+        col_names=["Energy", "Ordered", "Defects", "Hits"],
+    )
+
+    sim, frames = Simulation.load(files[0])
+    return data, sim.domain.n, sim.domain.r
+
+
+def simulation_data_proc(file: Path) -> Tuple[float, List[float], List[float]]:
+    sim, frames = Simulation.load(file)
+
+    alls = [[], [], []]
+    for frame_info in frames:
+        alls[0].append(frame_info["energy"])
+        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()
+
+    distincts = [[], [], []]
+    for j, frame_info in enumerate(sim.frames):
+        distincts[0].append(frame_info["energy"])
+        distincts[1].append(np.var(frame_info["stats"]["avg_radius"]) <= 1e-8)
+        distincts[2].append(
+            np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6)
+        )
+
+    distincts.append(counts)
+    distincts = list(zip(*sorted(zip(*distincts))))
+
+    return sim.domain.w / sim.domain.h, alls, distincts

scripts/torus_img.py

@@ -0,0 +1,110 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def torus_mesh(n, r, R):
+    theta = np.linspace(0, 2 * np.pi, n)
+    phi = np.linspace(0, 2 * np.pi, n)
+    theta, phi = np.meshgrid(theta, phi)
+    x = (R + r * np.cos(theta)) * np.cos(phi)
+    y = (R + r * np.cos(theta)) * np.sin(phi)
+    z = r * np.sin(theta)
+    return x, y, z
+
+
+def torus_line(n, r, R, u, v):
+    phi = np.linspace(0, 2 * np.pi * u, n)
+    theta = (v / u) * phi
+    theta, phi = theta % (2 * np.pi), phi % (2 * np.pi)
+
+    x = (R + r * np.cos(theta)) * np.cos(phi)
+    y = (R + r * np.cos(theta)) * np.sin(phi)
+    z = r * np.sin(theta)
+
+    line = np.column_stack([x, y, z])
+    x, y, z = R * np.cos(phi), R * np.sin(phi), 0 * phi
+    norm = line - np.column_stack([x, y, z])
+
+    return line, norm
+
+
+def torus_points(r, R, u, v, k):
+    phi = np.linspace(0, 2 * np.pi * u, k + 1)
+    theta = (v / u) * phi
+    theta, phi = theta % (2 * np.pi), phi % (2 * np.pi)
+
+    x = (R + r * np.cos(theta)) * np.cos(phi)
+    y = (R + r * np.cos(theta)) * np.sin(phi)
+    z = r * np.sin(theta)
+
+    return x, y, z
+
+
+def get_obscure_segs(arr):
+    bounds = [0, 0]
+    seg = 1
+    for i in range(1, len(arr)):
+        if arr[i] != arr[i - 1]:
+            bounds[seg] = i
+            bounds.append(i)
+            seg += 1
+    return bounds + [len(arr) - 1]
+
+
+def main():
+    fig = plt.figure(figsize=(10, 10), constrained_layout=True)
+    ax = fig.add_subplot(projection="3d")
+
+    ax.set_box_aspect(aspect=(1, 1, 1))
+
+    alpha = np.sqrt(3) * 2 / 5
+    size = 5
+
+    azim, elev = 0, 40
+    R, r = size, size * (1 - alpha) / alpha
+    u, v = 5, 2
+    k = 20
+
+    print(alpha, r / (r + R))
+
+    ax.set_zlim(-(r + R), r + R)
+    ax.view_init(elev, azim)
+
+    ax.plot_surface(
+        *torus_mesh(100, r, R), rcount=200, ccount=200, color="white", alpha=0.4
+    )
+
+    # Calculate when surface is facing or away and apply dotted vs solid lines.
+    az, el = azim * np.pi / 180 - np.pi, elev * np.pi / 180 - np.pi / 2
+    camera_vec = np.array(
+        [np.sin(el) * np.cos(az), np.sin(el) * np.sin(az), np.cos(el)]
+    )
+    line, norm = torus_line(500, r, R, u, v)
+    cond = np.dot(norm, camera_vec) >= 0
+    bounds = get_obscure_segs(cond)
+
+    # Did not account for normal vectors intersecting with main body, so
+    # we have manual overrides...
+    print(bounds)
+    bounds[3] = 220
+    print(bounds)
+
+    for i in range(len(bounds) - 1):
+        seg = line[bounds[i] : bounds[i + 1] + 1, :]
+        lx, ly, lz = seg[:, 0], seg[:, 1], seg[:, 2]
+        style = "solid" if cond[bounds[i]] else "dashed"
+        if i == 3:
+            style = "dashed"
+        if i == 4:
+            style = "dashed"
+        ax.plot(lx, ly, lz, color="blue", linewidth=2, linestyle=style)
+
+    # Display points.
+    ax.scatter(*torus_points(r, R, u, v, k), color="purple", s=50)
+    # ax.plot(*torus_line(200, 2, 1, 2, 1), color="orange", linewidth=3)
+    plt.show()
+    # plt.savefig(f"Torus - N{k}({u},{v}).png", bbox_inches="tight")
+
+
+if __name__ == "__main__":
+    main()

scripts/torus_mesh.py

@@ -0,0 +1,325 @@
+import argparse, numpy as np, os
+from stl.stl import ASCII
+from stl import mesh
+from squish import DomainParams, Simulation, ordered
+from squish.common import OUTPUT_DIR
+from scipy.spatial import Voronoi
+
+SUBDIVISION_AMOUNT = 8
+
+
+def not_in_order(i, j):
+    if i < j:
+        return j - i != 1
+    else:
+        return i - j == 1
+
+
+def centroid(site, verts):
+    area, c = 0, 0
+    v = verts - site
+    # print(v)
+    for i in range(len(v)):
+        x, y = v[i], v[(i + 1) % len(v)]
+        a = np.cross(x, y)
+        area += a
+        c += a * (x + y)
+
+    return c / (3 * area) + site
+
+
+def flat_sheet(k):
+    x = np.linspace(-np.pi, np.pi, 101)
+    x = np.transpose([np.tile(a, len(a)), np.repeat(a, len(a))])
+    all_verts = np.hstack((a, np.zeros((len(a), 1))))
+
+    xx = np.array([[x for x in range(101 * 100) if (x - 100) % 101 > 0]]).T
+
+    b = np.hstack((xx, xx + 1, xx + 101))
+    c = np.hstack((xx + 1, xx + 102, xx + 101))
+    all_faces = np.concatenate((b, c))
+
+    return all_verts, all_faces
+
+
+def torus_transform(v, R, r):
+    new_verts = np.empty(v.shape)
+    # Rotate by 90 degrees, and then shift to range [0, 2\pi]
+    v[:, :2] = np.matmul(np.array([[0, -1], [1, 0]]), v[:, :2].T).T
+    v[:, :2] = v[:, :2] % (2 * np.pi)
+    # v[:, 0], v[:, 1] = v[:, 1], v[:, 0]
+
+    new_verts[:, 0] = (R + r * np.cos(v[:, 1])) * np.cos(v[:, 0])
+    new_verts[:, 1] = (R + r * np.cos(v[:, 1])) * np.sin(v[:, 0])
+    new_verts[:, 2] = r * np.sin(v[:, 1])
+    return new_verts
+
+
+def flat_region(c, verts, height):
+    m = len(verts)
+
+    centers = np.hstack((np.array([c, c]), np.array([[height], [0]])))
+    v_top = np.hstack((verts, height * np.ones((m, 1))))
+    v_bot = np.hstack((verts, np.zeros((m, 1))))
+
+    v = np.concatenate((centers, v_top, v_bot))
+
+    vi = np.atleast_2d(np.arange(m, dtype=int)).T
+    cent_top, cent_bot = np.zeros((m, 1), dtype=int), np.ones((m, 1), dtype=int)
+    vert_top_m, vert_top_p = vi + 2, (vi + 1) % m + 2
+    vert_bot_m, vert_bot_p = vi + m + 2, (vi + 1) % m + m + 2
+
+    top_face = np.hstack((cent_top, vert_top_m, vert_top_p))
+    bot_face = np.hstack((cent_bot, vert_bot_p, vert_bot_m))
+
+    f = np.concatenate((top_face, bot_face))
+    return v, f, 2
+
+
+def torus_region(c, verts, height):
+    sd = SUBDIVISION_AMOUNT
+    d = verts - c
+    m = len(verts)
+    # 1 + 6 + 12 ... verts, 1 + 3 + 5 ... faces.
+    v, f = (
+        np.empty((1 + m * sd * (sd - 1) // 2, 2)),
+        np.empty((m * (sd - 1) ** 2, 3), dtype=int),
+    )
+    v[0] *= 0
+
+    # Vertices are ordered by rings going out from the center.
+    v_inds = np.arange(-1, sd)
+    v_inds = 1 + m * v_inds * (v_inds + 1) // 2
+    v_inds[0] = 0
+
+    for i in range(1, sd):
+        fj, fk = m * (i - 1) ** 2, m * i ** 2
+        region_verts = [
+            np.linspace(i * d[j] / (sd - 1), i * d[(j + 1) % m] / (sd - 1), i + 1)[:-1]
+            for j in range(m)
+        ]
+        v[v_inds[i] : v_inds[i + 1]] = np.concatenate(region_verts)
+
+        faces = []
+        for j in range(m):
+            # v_off_inds = v_inds + j * np.arange(sd + 1)
+            r1, r2 = (
+                (j * (i - 1) + np.arange(i + 1)) % (m * (i - 1)) + v_inds[i - 1],
+                (j * i + np.arange(i + 2)) % (m * i) + v_inds[i],
+            )
+            r1, r2 = np.atleast_2d(r1).T, np.atleast_2d(r2).T
+            faces.append(np.hstack((r1[:-1], r2[:-2], r2[1:-1])))
+            faces.append(np.hstack((r1[:-2], r2[1:-2], r1[1:-1])))
+
+        f[fj:fk] = np.concatenate(faces)
+
+    v += c
+    v = np.hstack((v, height * np.ones((len(v), 1))))
+    return v, f, v_inds[sd - 1]
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Generates 3D model for equilibrium.")
+    parser.add_argument(
+        "sims_path", metavar="sim_dir", help="simulation to obtain equilibrium from"
+    )
+    parser.add_argument("frame_num", metavar="FRAME_NUM", type=int, help="frame number")
+    parser.add_argument("size", metavar="SIZE", type=int, help="width in mm")
+    parser.add_argument(
+        "--torus", dest="torus", action="store_true", help="generates torus model"
+    )
+    args = parser.parse_args()
+
+    torus = args.torus
+    size = args.size * 10
+
+    # Get desired frame and load.
+    sim, frames = Simulation.load(args.sims_path)
+    frames = list(frames)
+    # num = 100
+    # nums = [838. 348. 14. 664, 725, 974]
+    # frames.sort(key=lambda x: x["energy"])
+    # frame = frames[-32]
+    # for i, f in enumerate(frames):
+    #    if i == 858:
+    frame = frames[args.frame_num]
+    n, w, h = frame["domain"][0], frame["domain"][1], frame["domain"][2]
+    frame = sim.energy.mode(*frame["domain"], frame["arr"])
+
+    # Set up size and scaling variables.
+    dim = np.array([w, h])
+    vor = frame.vor_data
+
+    alpha = w / h
+    if torus:
+        r = alpha * size / 2
+        R = size / 2 - r
+        height_bounds = (0.5 * r, 1.25 * r)
+        scale = np.array([2 * np.pi, 2 * np.pi])
+    else:
+        height_bounds = (size / 10, 0.75 * size)
+        scale = np.array([alpha * size, size])
+
+    # Rescale and translate domain.
+    points = (vor.points - dim / 2) * scale / dim
+    vertices = (vor.vertices - dim / 2) * scale / dim
+
+    # Obtain height scaling.
+    heights = frame.stats["site_energies"] - ordered.e_hex(sim.domain)
+    heights -= np.min(heights)
+    heights *= (height_bounds[1] - height_bounds[0]) / np.max(heights)
+    heights += height_bounds[0]
+
+    # Prepare oriented regions and site_vert list
+    sites = points[:n]
+    regions = []
+    for i, x in enumerate(vor.point_region[:n]):
+        region = vor.regions[x]
+        x = sites[i]
+        p, q = vertices[region[0]] - x, vertices[region[1]] - x
+        regions.append(region if np.cross(p, q) >= 0 else region[::-1])
+
+    site_verts = []
+    for region in regions:
+        site_verts.append(vertices[region])
+
+    all_verts, all_faces = [], []
+    offsets = np.empty((n,), dtype=int)
+    for i, verts in enumerate(site_verts):
+        c = centroid(sites[i], verts)
+        if torus:
+            face_verts, faces, offset = torus_region(c, verts, heights[i])
+        else:
+            face_verts, faces, offset = flat_region(c, verts, heights[i])
+
+        old_len = len(all_verts)
+        if i == 0:
+            all_verts, all_faces = face_verts, faces
+        else:
+            all_verts = np.concatenate((all_verts, face_verts))
+            all_faces = np.concatenate((all_faces, faces + old_len))
+
+        offsets[i] = old_len + offset
+
+    # Merge regions
+    for i, x in enumerate(vor.ridge_points):
+        x, y = x[0], x[1]  # Site indices.
+        no_x, no_y = x >= sim.domain.n, y >= sim.domain.n
+        if no_x and no_y:
+            continue
+
+        adj_verts = vor.ridge_vertices[i]
+        v1, v2 = vertices[adj_verts[0]], vertices[adj_verts[1]]
+
+        bot_face, top_face = None, None
+        if torus:
+            sd = SUBDIVISION_AMOUNT - 1
+            if no_x:
+                x, y = y, x
+            y = y % n
+
+            v1_x_ind = regions[x].index(adj_verts[0])
+            v2_x_ind = regions[x].index(adj_verts[1])
+            # Need to find matching vertices
+            v1, v2 = v1 % (2 * np.pi), v2 % (2 * np.pi)
+            y_verts = vertices[regions[y]] % (2 * np.pi)
+            v1_y_ind = np.argmin(np.linalg.norm(y_verts - v1, axis=1))
+            v2_y_ind = np.argmin(np.linalg.norm(y_verts - v2, axis=1))
+
+            # We want x lower than y.
+            off_x, off_y = offsets[x], offsets[y]
+            if heights[x] > heights[y]:
+                v1_x_ind, v1_y_ind = v1_y_ind, v1_x_ind
+                v2_x_ind, v2_y_ind = v2_y_ind, v2_x_ind
+                off_x, off_y = off_y, off_x
+
+            if not_in_order(v1_x_ind, v2_x_ind):
+                v1_x_ind, v2_x_ind = v2_x_ind, v1_x_ind
+
+            if not_in_order(v1_y_ind, v2_y_ind):
+                v1_y_ind, v2_y_ind = v2_y_ind, v1_y_ind
+
+            sds = np.atleast_2d(np.arange(sd)).T
+            xm, ym = sd * v1_x_ind + sds, sd * v1_y_ind + sds
+            xp, yp = xm + 1, ym + 1
+            xp[-1], yp[-1] = sd * v2_x_ind, sd * v2_y_ind
+
+            xm += off_x
+            xp += off_x
+            ym += off_y
+            yp += off_y
+
+            bot_face = np.hstack((xm, ym[::-1], xp))
+            top_face = np.hstack((xm, yp[::-1], ym[::-1]))
+
+        else:
+            if no_x ^ no_y:
+                if no_x:
+                    x = y
+                m = len(regions[x])
+                v1_ind = regions[x].index(adj_verts[0])
+                v2_ind = regions[x].index(adj_verts[1])
+                if not_in_order(v1_ind, v2_ind):
+                    v1_ind, v2_ind = v2_ind, v1_ind
+
+                v1_ind += offsets[x]
+                v2_ind += offsets[x]
+                bot_face = np.array([[v1_ind, m + v1_ind, m + v2_ind]], dtype=int)
+                top_face = np.array([[v1_ind, m + v2_ind, v2_ind]], dtype=int)
+            else:
+                v1_x_ind = regions[x].index(adj_verts[0])
+                v2_x_ind = regions[x].index(adj_verts[1])
+                v1_y_ind = regions[y].index(adj_verts[0])
+                v2_y_ind = regions[y].index(adj_verts[1])
+
+                # We want x lower than y.
+                off_x, off_y = offsets[x], offsets[y]
+                if heights[x] > heights[y]:
+                    v1_x_ind, v1_y_ind = v1_y_ind, v1_x_ind
+                    v2_x_ind, v2_y_ind = v2_y_ind, v2_x_ind
+                    off_x, off_y = off_y, off_x
+
+                if not_in_order(v1_x_ind, v2_x_ind):
+                    v1_x_ind, v2_x_ind = v2_x_ind, v1_x_ind
+
+                if not_in_order(v1_y_ind, v2_y_ind):
+                    v1_y_ind, v2_y_ind = v2_y_ind, v1_y_ind
+
+                v1_x_ind += off_x
+                v2_x_ind += off_x
+                v1_y_ind += off_y
+                v2_y_ind += off_y
+
+                bot_face = np.array([[v1_x_ind, v1_y_ind, v2_x_ind]], dtype=int)
+                top_face = np.array([[v1_x_ind, v2_y_ind, v1_y_ind]], dtype=int)
+
+        all_faces = np.concatenate((all_faces, bot_face, top_face))
+
+    if torus:
+        all_verts = torus_transform(all_verts, R, all_verts[:, 2])
+
+    # Output into mesh
+    surf = mesh.Mesh(np.zeros(all_faces.shape[0], dtype=mesh.Mesh.dtype))
+    for i, f in enumerate(all_faces):
+        for j in range(3):
+            surf.vectors[i][j] = all_verts[f[j], :]
+
+    m_type = "Torus" if torus else "Flat"
+    out = OUTPUT_DIR / f"N{n} - {m_type} - {args.frame_num}.stl"
+    # temp = f"n{n}{m_type.lower()}{num}.stl"
+    surf.save(out)
+    # os.rename(temp, out)
+    print(f"Wrote to {out}.")
+    # n, alpha = 20, np.sqrt(3) * 2 / 5
+    # u, v = 5, 2
+    # domain = DomainParams(n, alpha, 1, 1)
+    # size = 1
+
+    # R, r = size, size * (1 - alpha) / alpha
+
+    # sites = ordered.sites(domain, (u, v))
+
+
+if __name__ == "__main__":
+    np.seterr(divide="ignore")
+    main()

scripts/vee.py

@@ -0,0 +1,114 @@
+import numpy as np, os
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
+
+from squish import Simulation, DomainParams, Energy, ordered
+from squish.common import OUTPUT_DIR
+from script_tools import (
+    RC_SETTINGS,
+    get_args,
+    get_data,
+    get_simulation_data,
+    get_ordered_data,
+)
+
+NAME = "VEE"
+
+
+def main():
+    sims_path, regen = get_args(
+        "Various graphs for single N data.",
+        "folder that contains simulation data at various aspect ratios for some N",
+    )
+
+    data, n, r = get_data(
+        sims_path / "package.pkl", get_simulation_data, args=(sims_path,), regen=regen
+    )
+    domain, alphas = DomainParams(n, 1, 1, r), data["all"]["alpha"]
+    ordered_data = get_data(
+        OUTPUT_DIR / "OrderedCache" / f"{n}.pkl",
+        get_ordered_data,
+        args=(domain, alphas),
+        regen=regen,
+    )
+
+    min_disorder, max_disorder = [], []
+    for i, energies in enumerate(data["distinct"]["Energy"]):
+        disorder_energies = []
+        for j, energy in enumerate(energies):
+            if not data["distinct"]["Ordered"][i][j]:
+                disorder_energies.append(energy)
+        min_disorder.append(min(disorder_energies))
+        max_disorder.append(max(disorder_energies))
+
+    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
+    min_order = np.array(min_order) / domain.n - e_hex
+
+    hex_ratios = ordered.hexagon_alpha(domain.n)
+
+    plt.rcParams.update(RC_SETTINGS)
+    fig = plt.figure(figsize=(15, 15))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    ax.plot(alphas, 100 * min_order, color="C1", label="Min Ordered", zorder=10)
+    ax.plot(alphas, 100 * min_disorder, color="C0", label="Min Disordered")
+
+    ax.plot(
+        alphas,
+        100 * max_disorder,
+        color="C0",
+        linestyle="dotted",
+        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)
+
+    space = np.arange(0, 6.6, 0.3)
+    ax.set_ylim(-0.15, 6.6)
+    ax.set_yticks(space[:-1])
+    ax.ticklabel_format(axis="y", style="sci")
+
+    ax.set_xlabel("Aspect Ratio")
+    ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$")
+
+    ax.legend(loc=(0.23, 0.5))
+    # ax.legend()
+    ax.grid(zorder=0)
+
+    props = dict(boxstyle="round", facecolor="white", alpha=0.8, zorder=20)
+    ax.text(
+        0.873,
+        0.97,
+        f"N={domain.n}",
+        transform=ax.transAxes,
+        verticalalignment="top",
+        bbox=props,
+    )
+
+    out = OUTPUT_DIR / f"{NAME}-N{domain.n}.png"
+    fig.savefig(out)
+    print(f"Wrote to {out}.")
+
+
+if __name__ == "__main__":
+    os.environ["QT_LOGGING_RULES"] = "*=false"
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Program terminated by user.")

scripts/vee_diff_pod.py

@@ -0,0 +1,95 @@
+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,
+    get_simulation_data,
+    get_ordered_data,
+)
+
+NAME = "VEEDiff-PoD"
+
+
+def main():
+    sims_path, regen = get_args(
+        "Scatter plot of VEE difference and disordered equilibria.",
+        "folders that contains various N simulations to plot",
+    )
+
+    packages = []
+    for fol in sims_path.iterdir():
+        if fol.is_file():
+            continue
+
+        data, n, r = get_data(
+            fol / "package.pkl", get_simulation_data, args=(fol,), regen=regen
+        )
+        domain, alphas = DomainParams(n, 1, 1, r), data["all"]["alpha"]
+        ordered_data = get_data(
+            OUTPUT_DIR / "OrderedCache" / f"{n}.pkl",
+            get_ordered_data,
+            args=(domain, alphas),
+            regen=regen,
+        )
+
+        packages.append([data, ordered_data, domain])
+
+    packages.sort(key=lambda x: x[2].n)
+
+    plt.rcParams.update(RC_SETTINGS)
+
+    fig = plt.figure(figsize=(20, 9))
+    gs = fig.add_gridspec(1, 1)
+    ax = fig.add_subplot(gs[0])
+
+    for j, package in enumerate(packages):
+        data, ordered_data, domain = package
+
+        min_disorder, max_disorder = [], []
+        for i, energies in enumerate(data["distinct"]["Energy"]):
+            disorder_energies = []
+            for j, energy in enumerate(energies):
+                if not data["distinct"]["Ordered"][i][j]:
+                    disorder_energies.append(energy)
+            min_disorder.append(min(disorder_energies))
+            max_disorder.append(max(disorder_energies))
+
+        min_order = []
+        for i, energies in enumerate(ordered_data["Energy"]):
+            min_order.append(energies[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
+        min_order = np.array(min_order) / domain.n - e_hex
+
+        all_disorder_count = []
+        for disorders in data["all"]["Ordered"]:
+            all_disorder_count.append(
+                100 * np.count_nonzero(disorders == False) / len(disorders)
+            )
+
+        ax.scatter(
+            100 * (min_order - min_disorder), all_disorder_count, label=f"N={domain.n}"
+        )
+
+    ax.set_ylabel("POD")
+    ax.set_xlabel(r"VEE Difference $\left[\times 10^{2}\right]$")
+    ax.yaxis.set_major_formatter(mtick.PercentFormatter())
+    ax.set_ylim(48, 102)
+    ax.set_yticks(np.arange(50, 101, 10))
+
+    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/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()

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = squish
-version = 0.1.5.post2
+version = 0.1.7
 author = Kenneth Jao
 author_email = ksjdragon@gmail.com
 description = squish is Python program which perform simulations for the flow of 'soft' or 'compressible' objects under some energy in a periodic domain.
@@ -21,9 +21,9 @@ zip_safe = False
 packages = squish
 python_requires = >= 3.8
 install_requires =
-	numpy == 1.21.2
-	scipy == 1.7.1
-	matplotlib == 3.4.3
+	numpy >= 1.21.2
+	scipy >= 1.7.1
+	matplotlib >= 3.4.3
 
 [options.entry_points]
 console_scripts =

squish/core.pxd

@@ -51,6 +51,7 @@ ctypedef Matrix2x2 (*MatrixCopySclOp)(Matrix2x2*, FLOAT_T) nogil
 
 ctypedef struct VectorSelfOps:
     Vector2D* (*neg)(Vector2D*) nogil
+    Vector2D* (*rot)(Vector2D*) nogil
 
     VectorSelfVecOp vadd
     VectorSelfVecOp vsub
@@ -66,6 +67,7 @@ ctypedef struct VectorSelfOps:
 
 ctypedef struct VectorCopyOps:
     Vector2D (*neg)(Vector2D*) nogil
+    Vector2D (*rot)(Vector2D*) nogil
 
     VectorCopyVecOp vadd
     VectorCopyVecOp vsub
@@ -81,6 +83,7 @@ ctypedef struct VectorCopyOps:
 
 ctypedef struct MatrixSelfOps:
     Matrix2x2* (*neg)(Matrix2x2*) nogil
+    Matrix2x2* (*T)(Matrix2x2*) nogil
 
     MatrixSelfMatOp madd
     MatrixSelfMatOp msub
@@ -96,6 +99,7 @@ ctypedef struct MatrixSelfOps:
 
 ctypedef struct MatrixCopyOps:
     Matrix2x2 (*neg)(Matrix2x2*) nogil
+    Matrix2x2 (*T)(Matrix2x2*) nogil
 
     MatrixCopyMatOp madd
     MatrixCopyMatOp msub
@@ -115,7 +119,7 @@ ctypedef struct Vector2D:
     VectorCopyOps copy
 
     bint (*equals)(Vector2D*, Vector2D) nogil
-    Vector2D (*rot)(Vector2D*) nogil
+    Matrix2x2 (*vecmul)(Vector2D*, Vector2D) nogil
     FLOAT_T (*dot)(Vector2D*, Vector2D) nogil
     FLOAT_T (*mag)(Vector2D*) nogil
 

squish/core.pyx

@@ -18,19 +18,19 @@ cdef MatrixCopyOps MCO
 
 VSO.neg, VSO.vadd, VSO.vsub, VSO.vmul, VSO.vdiv, VSO.sadd, VSO.ssub, VSO.smul, VSO.sdiv = \
     v_neg_s, v_vadd_s, v_vsub_s, v_vmul_s, v_vdiv_s, v_sadd_s, v_ssub_s, v_smul_s, v_sdiv_s
-VSO.matmul = v_matmul_s
+VSO.matmul, VSO.rot = v_matmul_s, rot_s
 
 VCO.neg, VCO.vadd, VCO.vsub, VCO.vmul, VCO.vdiv, VCO.sadd, VCO.ssub, VCO.smul, VCO.sdiv = \
     v_neg_c, v_vadd_c, v_vsub_c, v_vmul_c, v_vdiv_c, v_sadd_c, v_ssub_c, v_smul_c, v_sdiv_c
-VCO.matmul = v_matmul_c
+VCO.matmul, VCO.rot = v_matmul_c, rot_c
 
 MSO.neg, MSO.madd, MSO.msub, MSO.mmul, MSO.mdiv, MSO.sadd, MSO.ssub, MSO.smul, MSO.sdiv = \
     m_neg_s, m_madd_s, m_msub_s, m_mmul_s, m_mdiv_s, m_sadd_s, m_ssub_s, m_smul_s, m_sdiv_s
-MSO.matmul = m_matmul_s
+MSO.matmul, MSO.T = m_matmul_s, m_transpose_s
 
 MCO.neg, MCO.madd, MCO.msub, MCO.mmul, MCO.mdiv, MCO.sadd, MCO.ssub, MCO.smul, MCO.sdiv = \
     m_neg_c, m_madd_c, m_msub_c, m_mmul_c, m_mdiv_c, m_sadd_c, m_ssub_c, m_smul_c, m_sdiv_c
-MCO.matmul = m_matmul_c
+MCO.matmul, MCO.T = m_matmul_c, m_transpose_c
 
 """
 If bound checking is desired, uncomment out ..._valid_indices functions.
@@ -143,7 +143,7 @@ cdef inline Vector2D _Vector2D(FLOAT_T x, FLOAT_T y) nogil:
     vec.x, vec.y = x, y
     vec.self, vec.copy = VSO, VCO
 
-    vec.equals, vec.rot, vec.dot, vec.mag = &v_equals, &rot, &dot, &mag
+    vec.equals, vec.vecmul, vec.dot, vec.mag = &v_equals, &v_vecmul, &dot, &mag
 
     return vec
 
@@ -200,6 +200,10 @@ cdef inline Vector2D* v_matmul_s(Vector2D* self, Matrix2x2 m) nogil:
     self.x, self.y = self.x*m.a + self.y*m.c, self.x*m.b + self.y*m.d
     return self
 
+cdef inline Vector2D* rot_s(Vector2D* self) nogil:
+    self.x, self.y = -self.y, self.x
+    return self
+
 cdef inline Vector2D v_neg_c(Vector2D* self) nogil:
     return _Vector2D(-self.x, -self.y)
 
@@ -232,7 +236,7 @@ cdef inline Vector2D v_matmul_c(Vector2D* self, Matrix2x2 m) nogil:
         self.x*m.a + self.y*m.c, self.x*m.b + self.y*m.d
     )
 
-cdef inline Vector2D rot(Vector2D* self) nogil:
+cdef inline Vector2D rot_c(Vector2D* self) nogil:
     return _Vector2D(-self.y, self.x)
 
 cdef inline FLOAT_T dot(Vector2D* self, Vector2D w) nogil:
@@ -241,6 +245,9 @@ cdef inline FLOAT_T dot(Vector2D* self, Vector2D w) nogil:
 cdef inline FLOAT_T mag(Vector2D* self) nogil:
     return <FLOAT_T>sqrt(<double>(self.x*self.x + self.y*self.y))
 
+cdef inline Matrix2x2 v_vecmul(Vector2D* self, Vector2D v) nogil:
+    return _Matrix2x2(self.x*v.x, self.x*v.y, self.y*v.x, self.y*v.y)
+
 
 #### Matrix2x2 Methods ####
 
@@ -329,6 +336,10 @@ cdef inline Matrix2x2* m_matmul_s(Matrix2x2* self, Matrix2x2 m) nogil:
         self.c*m.a + self.d*m.c, self.c*m.b + self.d*m.d
     return self
 
+cdef inline Matrix2x2* m_transpose_s(Matrix2x2* self) nogil:
+    self.b, self.c = self.c, self.b
+    return self
+
 cdef inline Matrix2x2 m_neg_c(Matrix2x2* self) nogil:
     return _Matrix2x2(-self.a, -self.b, -self.c, -self.d)
 
@@ -361,3 +372,6 @@ cdef inline Matrix2x2 m_matmul_c(Matrix2x2* self, Matrix2x2 m) nogil:
         self.a*m.a + self.b*m.c, self.a*m.b + self.b*m.d,
         self.c*m.a + self.d*m.c, self.c*m.b + self.d*m.d
     )
+
+cdef inline Matrix2x2 m_transpose_c(Matrix2x2* self) nogil:
+    return _Matrix2x2(self.a, self.c, self.b, self.d)

squish/diagram.py

@@ -1,12 +1,19 @@
 from __future__ import annotations
 from typing import Tuple, List, Optional
 
-import matplotlib.pyplot as plt, numpy as np, os
+import matplotlib as mpl
+import matplotlib.cm as cm
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable
+
+import numpy as np, os, math
 from matplotlib.ticker import MaxNLocator, FormatStrFormatter
-from scipy.spatial import Voronoi, voronoi_plot_2d
+from scipy.spatial import Voronoi
 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]]
@@ -118,96 +125,164 @@ class Diagram:
         self.diagrams = np.atleast_2d(diagrams)
         self.cumulative = cumulative
 
-    def generate_frame(self, frame: int, mode: str, fol: str, name: str = None) -> None:
+    def generate_frame(
+        self, fig: Figure, frame: int, mode: str, fol: str, name: str = None
+    ) -> None:
         if mode not in ["save", "open"]:
             raise ValueError("Not a valid mode for diagrams!")
 
         shape = self.diagrams.shape
-        fig, axes = plt.subplots(*shape, figsize=(shape[1] * 8, shape[0] * 8))
+        gs = fig.add_gridspec(shape[0], shape[1])
+        # fig, axes = plt.subplots(*shape, figsize=(shape[1] * 15, shape[0] * 15))
+
         if self.diagrams.shape == (1, 1):
-            getattr(self, str(self.diagrams[0][0]) + "_plot")(frame, axes)
+            ax = fig.add_subplot(gs[0])
+            getattr(self, str(self.diagrams[0][0]) + "_plot")(frame, ax)
         else:
-            axes = np.atleast_2d(axes)
+            # axes = np.atleast_2d(axes)
             it = np.nditer(self.diagrams, flags=["multi_index"])
             for diagram in it:
                 if diagram == "":
                     continue
-                getattr(self, str(diagram) + "_plot")(frame, axes[it.multi_index])
 
-        plt.tight_layout()
+                ax = fig.add_subplot(gs[it.multi_index])
+                getattr(self, str(diagram) + "_plot")(frame, ax)
+
+        # plt.tight_layout()
         if name is None:
             name = f"img{frame:05}.png"
 
         if mode == "save":
             plt.savefig(self.sim.path / fol / name)
-            plt.close(fig)
+            fig.clear()
         elif mode == "show":
             plt.show()
 
     def voronoi_plot(self, i: int, ax: AxesSubplot) -> None:
         domain = self.sim.domains[i]
         n, w, h = domain.n, domain.w, domain.h
+        flip = w <= h
         scale = 1.5
-        area = n <= 60
 
-        voronoi_plot_2d(
-            self.sim.voronois[i], ax, show_vertices=False, point_size=7 - n / 100
-        )
-        ax.plot([-w, 2 * w], [0, 0], "r")
-        ax.plot([-w, 2 * w], [h, h], "r")
-        ax.plot([0, 0], [-h, 2 * h], "r")
-        ax.plot([w, w], [-h, 2 * h], "r")
-        ax.axis("equal")
-        ax.set_xlim([(1 - scale) * w / 2, (1 + scale) * w / 2])
-        ax.set_ylim([(1 - scale) * h / 2, (1 + scale) * h / 2])
-        ax.title.set_text("Voronoi Visualization")
+        e_hex = ordered.e_hex(domain)
+        line_color, point_color = "white", "lightseagreen"
 
-        props = dict(boxstyle="round", facecolor="wheat", alpha=0.8)
+        # Make color map axis.
+        divider = make_axes_locatable(ax)
+        cax = divider.append_axes("right", size="5%", pad=0.05)
 
-        defects = {5: {"x": [], "y": []}, 7: {"x": [], "y": []}}
+        vor = self.sim.voronois[i]
 
-        for j in range(n):
+        # Mark location axis with 3 ticks.
+        ax.set_xticks(np.round([0, w / 2, w], 2))
+        ax.set_yticks(np.round([0, h / 2, h], 2))
+
+        # Obtain site energies, and make color map.
+        energies = self.sim.stats[i]["site_energies"][:n]
+        norm = mpl.colors.Normalize(vmin=-2, vmax=2, clip=True)
+        mapper = cm.ScalarMappable(norm=norm, cmap=cm.magma)
+        cbar = plt.colorbar(mapper, cax=cax, ticks=np.linspace(-2, 2, 5))
+
+        SYMM = []
+        if flip:
+            dup = int((h // w + 1) // 2 + 1)
+            for i in range(-dup, dup + 1):
+                SYMM += [[i, -1], [i, 0], [i, 1]]
+        else:
+            dup = int((w // h + 1) // 2 + 1)
+            for i in range(-dup, dup + 1):
+                SYMM += [[-1, i], [0, i], [1, i]]
+
+        SYMM = np.array(SYMM) * np.array([w, h])
+
+        regions = [vor.regions[x] for x in vor.point_region[:n]]
+        edges = np.array([len(r) for r in regions])
+
+        # Fill polygons with energy colormap.
+        sizes = [40 - n / 100, 200 - n / 100, 40 - n / 100, 200 - n / 100, 40 - n / 100]
+        markers = [".", "p", ".", "*", "."]
+        for j in range(5):
+            sites = vor.points[np.argwhere(edges == j + 4)]
+            if len(sites) == 0:
+                continue
             for s in SYMM:
-                vec = self.sim.voronois[i].points[j] + s * self.sim.domains[i].dim
+                ax.scatter(
+                    *(sites + s).T,
+                    marker=markers[j],
+                    s=sizes[j],
+                    color=point_color,
+                    zorder=1,
+                )
 
-                if area:
-                    txt = ax.text(
-                        *vec, str(round(self.sim.stats[i]["site_areas"][j], 3))
-                    )
-                    txt.set_clip_on(True)
+        polygons = []
+        for region in regions:
+            polygon = np.empty((2 * len(SYMM), len(region)), dtype=float)
+            for i, s in enumerate(SYMM):
+                polygon[i * 2 : (i + 1) * 2, :] = (vor.vertices[region] + s).T
+            polygons.append(polygon)
 
-                if self.sim.stats[i]["site_edge_count"][j] == 5:
-                    defects[5]["x"].append(vec[0])
-                    defects[5]["y"].append(vec[1])
-                elif self.sim.stats[i]["site_edge_count"][j] == 7:
-                    defects[7]["x"].append(vec[0])
-                    defects[7]["y"].append(vec[1])
+        colors = [mapper.to_rgba(energies[i] - e_hex) for i in range(n)]
+        for col, poly_set in zip(colors, polygons):
+            ax.fill(*poly_set, color=col, edgecolor="black", zorder=0)
 
-        ax.scatter(defects[5]["x"], defects[5]["y"], marker="p", color="C0")
-        ax.scatter(defects[7]["x"], defects[7]["y"], marker="*", color="C0")
+        ax.axis("equal")
 
+        # Periodic border
+        if flip:
+            bot, top = (1 - scale) * h / 2, (1 + scale) * h / 2
+            left, right = w / 2 - 0.88 * (top - bot) / 2, w / 2 + 0.88 * (top - bot) / 2
+            ax.plot([left, right], [0, 0], line_color, linewidth="4")
+            ax.plot([left, right], [h, h], line_color, linewidth="4")
+            for i in range(-dup, dup + 2):
+                ax.plot([i * w, i * w], [-h, 2 * h], line_color, linewidth="4")
+        else:
+            left, right = (1 - 0.85 * scale) * w / 2, (1 + 0.85 * scale) * w / 2
+            bot, top = (h / 2 - (scale * w) / 2, h / 2 + (scale * w) / 2)
+
+            ax.plot([0, 0], [bot, top], line_color, linewidth="4")
+            ax.plot([w, w], [bot, top], line_color, linewidth="4")
+            for i in range(-dup, dup + 2):
+                ax.plot([-w, 2 * w], [i * h, i * h], line_color, linewidth="4")
+
+        ax.set_xlim(left, right)
+        ax.set_ylim(bot, top)
+        ax.set_title("Voronoi Tessellation")
+
+        # Make VEE text in top left.
+        props = dict(boxstyle="round", facecolor="white", alpha=0.8, zorder=20)
         ax.text(
-            0.05,
-            0.95,
-            f"Energy: {self.sim.energies[i]}",
+            0.065,
+            0.935,
+            f"VEE: {sum(energies)/n - e_hex: .8f}",
             transform=ax.transAxes,
-            fontsize=14,
             verticalalignment="top",
             bbox=props,
         )
 
     def energy_plot(self, i: int, ax: AxesSubplot) -> None:
         ax.set_xlim([0, len(self.sim)])
+        e_hex = ordered.e_hex(self.sim.domains[i])
 
-        energies = self.sim.energies[: (i + 1)]
-        ax.plot(list(range(i + 1)), energies)
-        ax.title.set_text("Energy vs. Time")
+        vees = np.array(self.sim.energies) / self.sim.domains[i].n - e_hex
+
+        ax.plot(list(range(len(vees))), vees)
+        ax.scatter(
+            i,
+            vees[i],
+            s=250,
+            facecolors="none",
+            edgecolors="C6",
+            linewidth=4,
+            zorder=100,
+        )
+
+        # ax.title.set_text("VEE")
         # max_value = round(self.sim[0].energy)
         # min_value = round(self.sim[-1].energy)
         # diff = max_value-min_value
         # ax.set_yticks(np.arange(int(min_value-diff/5), int(max_value+diff/5), diff/25))
-        ax.set_xlabel("Iterations")
-        ax.set_ylabel("Energy")
+        ax.set_xlabel("Frame")
+        ax.set_ylabel("VEE")
         ax.grid()
 
     def site_areas_plot(self, i: int, ax: AxesSubplot) -> None:
@@ -305,14 +380,32 @@ class Diagram:
     def eigs_plot(self, i: int, ax: AxesSubplot) -> None:
         try:
             eigs = self.sim.stats[i]["eigs"]
+            for i, eig in enumerate(eigs):
+                if eig > 1e-4:
+                    ax.annotate(
+                        f"Coercivity: {eig:.5f}",
+                        xy=(i, eig),
+                        xytext=(50, -50),
+                        textcoords="offset points",
+                        arrowprops={"arrowstyle": "->", "color": "black"},
+                    )
+                    break
+
             ax.plot(
-                list(range(len(eigs))), eigs, marker="o", linestyle="dashed", color="C0"
+                list(range(40)),
+                eigs[:40],
+                marker="o",
+                linestyle="dashed",
+                color="C0",
+                markersize=8,
             )
-            ax.plot([0, len(eigs)], [0, 0], color="red")
+
+            ax.plot([0, 40], [0, 0], color="red")
         except KeyError:
             ax.text(0.5, 0.5, "Not Computed")
 
-        ax.title.set_text("Hessian Eigenvalues")
+        ax.grid()
+        ax.set_title("Sorted Hessian Eigenvalues")
         ax.set_xlabel("")
         ax.set_ylabel("Value")
 
@@ -321,11 +414,8 @@ class Diagram:
         self.sim.path.mkdir(exist_ok=True)
         (self.sim.path / fol).mkdir(exist_ok=True)
         combo_list = []
-        for i in range(cpu_count()):
-            start, end = (
-                int(i * len(frames) / cpu_count()),
-                int((i + 1) * len(frames) / cpu_count()),
-            )
+        for i in range(1):
+            start, end = (int(i * len(frames) / 1), int((i + 1) * len(frames) / 1))
             new_dia = Diagram(None, self.diagrams, self.cumulative)
             new_dia.sim = self.sim.slice(frames[start:end])
             combo_list.append(
@@ -345,10 +435,10 @@ class Diagram:
         if mode not in ["use_all", "sample"]:
             raise ValueError("Not a valid mode for videos!")
 
+        fps = 30
         if mode == "use_all":
             frames = list(range(len(self.sim)))
         elif mode == "sample":
-            fps = 30
             if len(self.sim) < fps * time:
                 frames = list(range(len(self.sim)))
                 fps = len(self.sim) / time
@@ -378,8 +468,36 @@ class Diagram:
 
 def render_frame_range(combo: Tuple[Diagram, str, int, int, int]) -> None:
     self, fol, offset, length, num_frames = combo
+
+    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,
+            "font.family": "cm",
+            "font.size": 40,
+            "text.usetex": True,
+            "text.latex.preamble": r"\usepackage{amsmath}",
+            "figure.constrained_layout.use": True,
+        }
+    )
+
+    # Generate figure here and pass in to prevent matplotlib memory leak.
+    shape = self.diagrams.shape
+    fig = plt.figure(figsize=(shape[1] * 15, shape[0] * 15))
     for i in range(length):
-        self.generate_frame(i, "save", fol, f"img{i+offset:05}.png")
+        self.generate_frame(fig, i, "save", fol, f"img{i+offset:05}.png")
         i = len(list((self.sim.path / fol).iterdir()))
         hashes = int(21 * i / num_frames)
         print(
@@ -388,3 +506,4 @@ def render_frame_range(combo: Tuple[Diagram, str, int, int, int]) -> None:
             flush=True,
             end="\r",
         )
+    plt.close(fig)

squish/energy.pyx

@@ -4,10 +4,10 @@ from cython.parallel import parallel, prange
 cimport numpy as np
 from libc.math cimport NAN, pi as PI, atanh
 from squish.core cimport INT_T, FLOAT_T, Vector2D, Matrix2x2, BitSet, \
-	_Vector2D, _Matrix2x2, _BitSet
+    _Vector2D, _Matrix2x2, _BitSet
 from squish.voronoi cimport Site, HalfEdge, EdgeCacheMap, VoronoiInfo, \
-	_Site, _HalfEdge, _EdgeCacheMap, _VoronoiInfo, VoronoiContainer, \
-	R, NAN_MATRIX, NAN_VECTOR
+    _Site, _HalfEdge, _EdgeCacheMap, _VoronoiInfo, VoronoiContainer, \
+    R, NAN_MATRIX, NAN_VECTOR
 
 #### Constants ####
 
@@ -18,228 +18,405 @@ cdef EdgeCacheMap AREA_ECM = _EdgeCacheMap(0, 4, 6, 8, 10, 12, 13, -1, -1, 14)
 cdef EdgeCacheMap RADIALT_ECM = _EdgeCacheMap(0, 4, 6, 8, -1, 10, 11, 12, 13, 14)
 
 cdef class AreaEnergy(VoronoiContainer):
-	"""
-	Class for formulas relevant to the Area energy.
-	:param n: [int] how many sites to generate.
-	:param w: [float] width of the bounding domain.
-	:param h: [float] height of the bounding domain.
-	:param r: [float] radius of zero energy circle.
-	:param sites: [np.ndarray] collection of sites.
-	"""
+    """
+    Class for formulas relevant to the Area energy.
+    :param n: [int] how many sites to generate.
+    :param w: [float] width of the bounding domain.
+    :param h: [float] height of the bounding domain.
+    :param r: [float] radius of zero energy circle.
+    :param sites: [np.ndarray] collection of sites.
+    """
 
-	attr_str = "area"
-	title_str = "Area"
+    attr_str = "area"
+    title_str = "Area"
 
-	def __init__(AreaEnergy self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r,
-					np.ndarray[FLOAT_T, ndim=2] site_arr):
-		self.edge_cache_map = &AREA_ECM
-		self.energy = 0.0
+    def __init__(AreaEnergy self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r,
+                    np.ndarray[FLOAT_T, ndim=2] site_arr):
+        self.edge_cache_map = &AREA_ECM
+        self.energy = 0.0
 
-		super().__init__(n, w, h, r, site_arr)
-		self.minimum = (<FLOAT_T>n)*(w*h/(<FLOAT_T>n)-PI*r**2)**2
+        super().__init__(n, w, h, r, site_arr)
+        self.minimum = (<FLOAT_T>n)*(w*h/(<FLOAT_T>n)-PI*r**2)**2
 
 
-	cdef void precompute(self) except *:
-		cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
-			self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
+    cdef void precompute(self) except *:
+        cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
+            self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
 
-		cdef Site xi
-		cdef HalfEdge em, e, ep
-		cdef Vector2D vdiff
+        cdef Site xi
+        cdef HalfEdge em, e, ep
+        cdef Vector2D vdiff
 
-		cdef FLOAT_T [:] site_energy = np.full(self.sites.shape[0], PI*self.r**2)
+        cdef FLOAT_T [:] site_energy = np.full(self.sites.shape[0], PI*self.r**2)
 
-		cdef INT_T i, j
-		for i in prange(self.sites.shape[0], nogil=True):
-			xi = _Site(i, &info)
-			e = xi.edge(&xi)
+        cdef INT_T i, j
+        for i in prange(self.sites.shape[0], nogil=True):
+            xi = _Site(i, &info)
+            e = xi.edge(&xi)
 
-			site_energy[i] = (xi.cache.area(&xi, NAN) - site_energy[i])**2
-			xi.cache.energy(&xi, site_energy[i])
+            site_energy[i] = (xi.cache.area(&xi, NAN) - site_energy[i])**2
+            xi.cache.energy(&xi, site_energy[i])
 
-			j = 0
-			while j < xi.edge_num(&xi):
-				em, ep = e.prev(&e), e.next(&e)
-				vdiff = em.origin(&em)
-				vdiff.self.vsub(&vdiff, ep.origin(&ep))
-				e.cache.dVdv(&e, R.vecmul(&R, vdiff))
-				e.cache.H(&e, VoronoiContainer.calc_H(em, e))
+            j = 0
+            while j < xi.edge_num(&xi):
+                em, ep = e.prev(&e), e.next(&e)
+                vdiff = em.origin(&em)
+                vdiff.self.vsub(&vdiff, ep.origin(&ep))
+                e.cache.dVdv(&e, R.vecmul(&R, vdiff))
+                e.cache.H(&e, VoronoiContainer.calc_H(em, e))
 
-				e = e.next(&e)
-				j = j + 1
+                e = e.next(&e)
+                j = j + 1
 
-		self.energy = np.sum(site_energy[:self.n])
+        self.energy = np.sum(site_energy[:self.n])
 
 
-	cdef void calc_grad(self) except *:
-		cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
-				self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
+    cdef void calc_grad(self) except *:
+        cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
+                self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
 
-		cdef Site xi, xf
-		cdef HalfEdge e, f
-		cdef Vector2D dedxi_p
-		cdef BitSet edge_set
+        cdef Site xi, xf
+        cdef HalfEdge e, f
+        cdef Vector2D dedxi_p
+        cdef BitSet edge_set
 
-		cdef INT_T num_edges = self.edges.shape[0]
-		cdef FLOAT_T A = PI*self.r**2
+        cdef INT_T num_edges = self.edges.shape[0]
+        cdef FLOAT_T A = PI*self.r**2
 
-		cdef FLOAT_T [:, ::1] dedx = np.zeros((self.n, 2), dtype=FLOAT)
+        cdef FLOAT_T [:, ::1] dedx = np.zeros((self.n, 2), dtype=FLOAT)
 
-		cdef INT_T i, j
-		for i in prange(self.n, nogil=True):
-			xi = _Site(i, &info)
-			e = xi.edge(&xi)
-			edge_set = _BitSet(num_edges)
-			j = 0
-			while j < xi.edge_num(&xi):	# Looping through site edges.
-				f = e
-				while True:	# Circling this vertex.
-					if not edge_set.add(&edge_set, f.arr_index):
-						xf = f.face(&f)
-						dedxi_p = f.cache.dVdv(&f, NAN_VECTOR)	#dVdv
-						dedxi_p.self.smul(&dedxi_p, xf.cache.area(&xf, NAN) - A)
-						dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX))
-						dedx[i][0] += dedxi_p.x
-						dedx[i][1] += dedxi_p.y
+        cdef INT_T i, j
+        for i in prange(self.n, nogil=True):
+            xi = _Site(i, &info)
+            e = xi.edge(&xi)
+            edge_set = _BitSet(num_edges)
+            j = 0
+            while j < xi.edge_num(&xi): # Looping through site edges.
+                f = e
+                while True: # Circling this vertex.
+                    if not edge_set.add(&edge_set, f.arr_index):
+                        xf = f.face(&f)
+                        dedxi_p = f.cache.dVdv(&f, NAN_VECTOR)  #dVdv
+                        dedxi_p.self.smul(&dedxi_p, xf.cache.area(&xf, NAN) - A)
+                        dedxi_p.self.matmul(&dedxi_p, e.cache.H(&e, NAN_MATRIX))
+                        dedx[i][0] += dedxi_p.x
+                        dedx[i][1] += dedxi_p.y
 
-					f = f.twin(&f)
-					f = f.next(&f)
-					if f.arr_index == e.arr_index:
-						break
+                    f = f.twin(&f)
+                    f = f.next(&f)
+                    if f.arr_index == e.arr_index:
+                        break
 
-				e = e.next(&e)
-				j = j + 1
-			edge_set.free(&edge_set)
-		self.grad = dedx
+                e = e.next(&e)
+                j = j + 1
+            edge_set.free(&edge_set)
+        self.grad = dedx
+
+
+    cdef void calc_hess(self) except *:
+        d = 10e-5
+        HE = np.zeros((2*self.n, 2*self.n))
+        new_sites = np.copy(self.site_arr)    # Maintain one copy for speed.
+        for i in range(self.n):
+            for j in range(2):
+                mod = self.w if j == 0 else self.h
+                new_sites[i][j] = (new_sites[i][j] + d) % mod
+                Ep = self.__class__(self.n, self.w, self.h, self.r, new_sites)
+                new_sites[i][j] = (new_sites[i][j] - 2*d) % mod
+                Em = self.__class__(self.n, self.w, self.h, self.r, new_sites)
+                new_sites[i][j] = (new_sites[i][j] + d) % mod
+
+                HE[:, 2*i+j] = ((Ep.gradient - Em.gradient)/(2*d)).flatten()
+
+        # Average out discrepencies, since it should be symmetric.
+        for i in range(2*self.n):
+            for j in range(i, 2*self.n):
+                HE[i][j] = (HE[i][j] + HE[j][i])/2
+                HE[j][i] = HE[i][j]
+        self.hess = HE
 
 
 cdef class RadialALEnergy(VoronoiContainer):
-	"""
-	Class for formulas relevant to the Area energy.
-	:param n: [int] how many sites to generate.
-	:param w: [float] width of the bounding domain.
-	:param h: [float] height of the bounding domain.
-	:param r: [float] radius of zero energy circle.
-	:param sites: [np.ndarray] collection of sites.
-	"""
+    """
+    Class for formulas relevant to the Area energy.
+    :param n: [int] how many sites to generate.
+    :param w: [float] width of the bounding domain.
+    :param h: [float] height of the bounding domain.
+    :param r: [float] radius of zero energy circle.
+    :param sites: [np.ndarray] collection of sites.
+    """
 
-	attr_str = "radial-al"
-	title_str = "Radial[AL]"
+    attr_str = "radial-al"
+    title_str = "Radial[AL]"
 
 
-	def __init__(AreaEnergy self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r,
-					np.ndarray[FLOAT_T, ndim=2] site_arr):
-		#self.edge_cache_map = &AREA_EDGE_CACHE_MAP
-		self.energy = 0.0
+    def __init__(AreaEnergy self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r,
+                    np.ndarray[FLOAT_T, ndim=2] site_arr):
+        #self.edge_cache_map = &AREA_EDGE_CACHE_MAP
+        self.energy = 0.0
 
-		super().__init__(n, w, h, r, site_arr)
+        super().__init__(n, w, h, r, site_arr)
 
 
-	cdef void precompute(self) except *:
-		cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
-			self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
+    cdef void precompute(self) except *:
+        cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
+            self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
 
-		pass
+        pass
 
 
-	cdef void calc_grad(self) except *:
-		cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
-				self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
+    cdef void calc_grad(self) except *:
+        cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
+                self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
 
-		pass
+        pass
+
+    cdef void calc_hess(self) except *:
+        pass
 
 
 cdef class RadialTEnergy(VoronoiContainer):
-	"""
-	Class for formulas relevant to the Area energy.
-	:param n: [int] how many sites to generate.
-	:param w: [float] width of the bounding domain.
-	:param h: [float] height of the bounding domain.
-	:param r: [float] radius of zero energy circle.
-	:param sites: [np.ndarray] collection of sites.
-	"""
+    """
+    Class for formulas relevant to the Area energy.
+    :param n: [int] how many sites to generate.
+    :param w: [float] width of the bounding domain.
+    :param h: [float] height of the bounding domain.
+    :param r: [float] radius of zero energy circle.
+    :param sites: [np.ndarray] collection of sites.
+    """
 
-	attr_str = "radial-t"
-	title_str = "Radial[T]"
-	def __init__(AreaEnergy self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r,
-					np.ndarray[FLOAT_T, ndim=2] site_arr):
-		self.edge_cache_map = &RADIALT_ECM
-		self.energy = 0.0
+    attr_str = "radial-t"
+    title_str = "Radial[T]"
+    def __init__(AreaEnergy self, INT_T n, FLOAT_T w, FLOAT_T h, FLOAT_T r,
+                    np.ndarray[FLOAT_T, ndim=2] site_arr):
+        self.edge_cache_map = &RADIALT_ECM
+        self.energy = 0.0
 
-		super().__init__(n, w, h, r, site_arr)
+        super().__init__(n, w, h, r, site_arr)
 
 
-	cdef void precompute(self) except *:
-		cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
-			self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
+    cdef void precompute(self) except *:
+        cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
+            self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
 
-		cdef Site xi
-		cdef HalfEdge e, ep
-		cdef Vector2D la
+        cdef Site xi
+        cdef HalfEdge e, ep
+        cdef Vector2D la
 
-		# All energy has a 2pir_0 term.
-		cdef FLOAT_T [:] site_energy = np.full(self.sites.shape[0], 2*PI*self.r**2)
-		cdef FLOAT_T [:] avg_radii = np.zeros(self.sites.shape[0])
-		cdef FLOAT_T sm, sp
+        # All energy has a 2pir_0 term.
+        cdef FLOAT_T [:] site_energy = np.full(self.sites.shape[0], 2*PI*self.r**2)
+        cdef FLOAT_T [:] avg_radii = np.zeros(self.sites.shape[0])
+        cdef FLOAT_T sm, sp
 
-		cdef INT_T i, j
-		for i in prange(self.sites.shape[0], nogil=True):
-			xi = _Site(i, &info)
-			e = xi.edge(&xi)
-			j = 0
-			while j < xi.edge_num(&xi):
-				ep = e.next(&e)
-				#e.cache.H(&e, VoronoiContainer.calc_H(em, e))
+        cdef INT_T i, j
+        for i in prange(self.sites.shape[0], nogil=True):
+            xi = _Site(i, &info)
+            e = xi.edge(&xi)
+            j = 0
+            while j < xi.edge_num(&xi):
+                ep = e.next(&e)
+                #e.cache.H(&e, VoronoiContainer.calc_H(em, e))
 
-				la = e.cache.la(&e, NAN_VECTOR)
-				sp = la.dot(&la, ep.cache.da(&ep, NAN_VECTOR)) # dap . la
-				sm = la.dot(&la, e.cache.da(&e, NAN_VECTOR)) # da . la
+                la = e.cache.la(&e, NAN_VECTOR)
+                sp = la.dot(&la, ep.cache.da(&ep, NAN_VECTOR)) # dap . la
+                sm = la.dot(&la, e.cache.da(&e, NAN_VECTOR)) # da . la
 
-				sp = sp / (ep.cache.da_mag(&ep, NAN) * e.cache.la_mag(&e, NAN))
-				sm = sm / (e.cache.da_mag(&e, NAN) * e.cache.la_mag(&e, NAN))
+                sp = sp / (ep.cache.da_mag(&ep, NAN) * e.cache.la_mag(&e, NAN))
+                sm = sm / (e.cache.da_mag(&e, NAN) * e.cache.la_mag(&e, NAN))
 
-				e.cache.calI(&e, <FLOAT_T>(atanh(<double>sp) - atanh(<double>sm)))
+                e.cache.calI(&e, <FLOAT_T>(atanh(<double>sp) - atanh(<double>sm)))
 
-				avg_radii[i] += e.cache.ya_mag(&e, NAN) * e.cache.calI(&e, NAN) / 2
+                avg_radii[i] += e.cache.ya_mag(&e, NAN) * e.cache.calI(&e, NAN) / 2
 
-				e = e.next(&e)
-				j = j + 1
+                e = e.next(&e)
+                j = j + 1
 
-			site_energy[i] += 2*(xi.cache.area(&xi, NAN) - self.r*avg_radii[i])
+            site_energy[i] += 2*(xi.cache.area(&xi, NAN) - self.r*avg_radii[i])
 
-			xi.cache.avg_radius(&xi, avg_radii[i]/(2*PI))
-			xi.cache.energy(&xi, site_energy[i])
+            xi.cache.avg_radius(&xi, avg_radii[i]/(2*PI))
+            xi.cache.energy(&xi, site_energy[i])
 
-		self.energy = np.sum(site_energy[:self.n])
+        self.energy = np.sum(site_energy[:self.n])
 
 
-	cdef void calc_grad(self) except *:
-		cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
-				self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
+    cdef void calc_grad(self) except *:
+        cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
+                self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
 
-		cdef Site xi
-		cdef HalfEdge e
-		cdef Vector2D dedxi_p
+        cdef Site xi
+        cdef HalfEdge e
+        cdef Vector2D dedxi_p
 
-		cdef FLOAT_T [:, ::1] dedx = np.zeros((self.n, 2), dtype=FLOAT)
+        cdef FLOAT_T [:, ::1] dedx = np.zeros((self.n, 2), dtype=FLOAT)
 
-		cdef INT_T i, j
-		for i in prange(self.n, nogil=True):
-			xi = _Site(i, &info)
-			e = xi.edge(&xi)
+        cdef INT_T i, j
+        for i in prange(self.n, nogil=True):
+            xi = _Site(i, &info)
+            e = xi.edge(&xi)
 
-			j = 0
-			while j < xi.edge_num(&xi):	# Looping through site edges.
-				dedxi_p = e.cache.ya(&e, NAN_VECTOR)
-				dedxi_p.self.smul(
-					&dedxi_p,
-					e.cache.calI(&e, NAN) / e.cache.ya_mag(&e, NAN)
-				)
+            j = 0
+            while j < xi.edge_num(&xi): # Looping through site edges.
+                dedxi_p = e.cache.ya(&e, NAN_VECTOR)
+                dedxi_p.self.smul(
+                    &dedxi_p,
+                    e.cache.calI(&e, NAN) / e.cache.ya_mag(&e, NAN)
+                )
 
-				dedx[i][0] += 2*self.r*dedxi_p.x
-				dedx[i][1] += 2*self.r*dedxi_p.y
+                dedx[i][0] += 2*self.r*dedxi_p.x
+                dedx[i][1] += 2*self.r*dedxi_p.y
 
-				e = e.next(&e)
-				j = j + 1
+                e = e.next(&e)
+                j = j + 1
 
-		self.grad = dedx
+        self.grad = dedx
+
+    cdef void calc_hess(self) except *:
+        cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
+                                             self.vertices, self.site_cache,
+                                             self.edge_cache, self.edge_cache_map)
+        cdef Site xi, xk
+        cdef HalfEdge em, e, ep, f
+        cdef Vector2D temp1, temp2
+        cdef Matrix2x2 dsite, q, tempm
+        cdef BitSet edge_set
+
+        cdef INT_T num_edges = self.edges.shape[0]
+
+        cdef FLOAT_T [:,:] HE = np.zeros((2*self.n, 2*self.n))
+
+        cdef INT_T i, j, k
+        for i in range(self.sites.shape[0]):
+            xi = _Site(i, &info)
+            e = xi.edge(&xi)
+
+            j = 0
+            while j < xi.edge_num(&xi):
+                em, ep = e.prev(&e), e.next(&e)
+                e.cache.H(&e, VoronoiContainer.calc_H(em, e))
+
+                e = e.next(&e)
+                j = j + 1
+
+
+        for i in range(self.n):
+            xi = _Site(i, &info)
+            e = xi.edge(&xi)
+            edge_set = _BitSet(num_edges)
+
+            j = 0
+            while j < xi.edge_num(&xi):
+                em, ep = e.prev(&e), e.next(&e)
+
+                # Calculating of p
+                temp1 = em.cache.la(&em, NAN_VECTOR)
+                temp1.self.rot(&temp1)
+                temp1.self.sdiv(
+                    &temp1,
+                    em.cache.la_mag(&em, NAN) * em.cache.ya_mag(&em, NAN) / 2
+                )
+                temp2 = e.cache.la(&e, NAN_VECTOR)
+                temp2.self.rot(&temp2)
+                temp2.self.sdiv(
+                    &temp2,
+                    e.cache.la_mag(&e, NAN) * e.cache.ya_mag(&e, NAN) / 2
+                )
+
+                temp1.self.vsub(&temp1, temp2)  # (lm/Am - l/A)
+
+                temp2 = e.cache.da(&e, NAN_VECTOR) # rot(d) / |d|
+                temp2.self.sdiv(&temp2, e.cache.da_mag(&e, NAN))
+                temp2.self.rot(&temp2)
+
+                dsite = temp1.vecmul(&temp1, temp2)
+
+                HE[2*i, 2*i] -= dsite.a
+                HE[2*i, 2*i+1] -= dsite.b
+                HE[2*i+1, 2*i] -= dsite.c
+                HE[2*i+1, 2*i+1] -= dsite.d
+
+                # Calculating of q
+                temp2 = e.cache.la(&e, NAN_VECTOR)
+                temp1 = temp2.copy.rot(&temp2)
+                q = temp1.vecmul(&temp1, temp2)
+                q.self.sdiv(&q, e.cache.la_mag(&e, NAN)**2)
+
+                q.self.msub(&q, R)
+                q.self.smul(
+                    &q,
+                    e.cache.calI(&e, NAN) / e.cache.la_mag(&e, NAN)
+                )
+
+                temp1 = e.cache.la(&e, NAN_VECTOR)
+                temp1.self.rot(&temp1)
+                tempm = temp1.vecmul(&temp1, temp1)
+                tempm.self.smul(
+                    &tempm,
+                    ep.cache.da_mag(&ep, NAN) - e.cache.da_mag(&e, NAN)
+                )
+                tempm.self.sdiv(
+                    &tempm,
+                    e.cache.la_mag(&e, NAN)**3 * e.cache.ya_mag(&e, NAN) / 2
+                )
+
+                q.self.madd(&q, tempm)
+
+                # Minus portion
+                temp2 = em.cache.la(&em, NAN_VECTOR)
+                temp1 = temp2.copy.rot(&temp2)
+                tempm = temp1.vecmul(&temp1, temp2)
+                tempm.self.sdiv(&tempm, em.cache.la_mag(&em, NAN)**2)
+
+                tempm = R.copy.msub(&R, tempm)
+                tempm.self.smul(
+                    &tempm,
+                    em.cache.calI(&em, NAN) / em.cache.la_mag(&em, NAN)
+                )
+
+                q.self.madd(&q, tempm)
+
+                temp1 = em.cache.la(&em, NAN_VECTOR)
+                temp1.self.rot(&temp1)
+                tempm = temp1.vecmul(&temp1, temp1)
+                tempm.self.smul(
+                    &tempm,
+                    e.cache.da_mag(&e, NAN) - em.cache.da_mag(&em, NAN)
+                )
+                tempm.self.sdiv(
+                    &tempm,
+                    em.cache.la_mag(&em, NAN)**3 * em.cache.ya_mag(&em, NAN) / 2
+                )
+
+                q.self.msub(&q, tempm)
+
+                # Add p to q, so p = p, q = p+q
+                q.self.madd(&q, dsite)
+
+                f = e
+                while True:
+                    xk = f.face(&f)
+                    k = xk.index(&xk) % self.n
+                    if k < 0:
+                        k = k + self.n
+
+                    if not edge_set.add(&edge_set, f.arr_index):
+                        tempm = q.copy.matmul(&q, f.cache.H(&f, NAN_MATRIX))
+
+                        HE[2*i, 2*k] += tempm.a
+                        HE[2*i, 2*k+1] += tempm.b
+                        HE[2*i+1, 2*k] += tempm.c
+                        HE[2*i+1, 2*k+1] += tempm.d
+
+                    f = f.twin(&f)
+                    f = f.next(&f)
+                    if f.arr_index == e.arr_index:
+                        break
+
+                e = e.next(&e)
+                j = j + 1
+            edge_set.free(&edge_set)
+        self.hess = -2*self.r*np.asarray(HE, dtype=FLOAT)
+        self.hess = (
+            ( np.asarray(self.hess, dtype=FLOAT)
+            + np.asarray(self.hess, dtype=FLOAT).T )
+            / 2
+        )

squish/ordered.py

@@ -2,9 +2,23 @@ from __future__ import annotations
 from typing import List, Tuple
 import numpy as np
 from numpy.linalg import norm as mag
-from math import gcd, sqrt, log, tan, atan, pi
+from math import gcd, sqrt, log, tan, atan, atanh, pi
+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 * (2 * pi) * R_HEX + 2 * pi * domain.r ** 2
+
+
+def toroidal_distance(domain: DomainParams, x: np.ndarray, y: np.ndarray) -> float:
+    dim = np.array([domain.w, domain.h])
+    absdist = np.abs(y - x)
+    wrap = dim * (absdist >= dim / 2)
+    return np.linalg.norm(wrap - absdist, axis=1)
 
 
 def configurations(domain: DomainParams) -> List[Config]:
@@ -25,36 +39,46 @@ def configurations(domain: DomainParams) -> List[Config]:
             except KeyError:
                 pass
 
+    for i in range(len(valid)):
+        valid.append((valid[i][1], valid[i][0]))
+
     return valid
 
 
 def get_config_generators(
     domain: DomainParams, config: Config
 ) -> Tuple[Config, Config]:
+    # x = sites(domain, config)
+    # x = x[x[:, 1] <= (domain.h / 2)]
+    # mag_x = toroidal_distance(domain, x, np.zeros(x.shape))
+    # x = x[np.argsort(mag_x)]
+    # print(x)
+
+    # return [tuple(x[1]), tuple(x[2])]
+
     n, w, h = domain.n, domain.w, domain.h
-    q1 = sites(domain, config)[1:]
-    all_sites = np.concatenate(
-        (q1, q1 - np.array([w, 0]), q1 - np.array([0, h]), q1 - domain.dim)
+    x = sites(domain, config)[1:]  # Remove 0
+    # Concatenate quadrants 1, 2 and 4.
+    x = np.concatenate(
+        (x, x - np.array([w, 0]), x - np.array([0, h]), x - np.array([w, h]))
     )
-    # Sort sites by magnitude and smallest.
-    all_sites = sorted(list(all_sites), key=lambda x: np.linalg.norm(x))
-    v = all_sites[0]  # Smallest vector set to v.
 
-    all_sites = np.array(all_sites[1:])  # Remove v from search set.
+    # Reduce search space
+    # x = x[(np.abs(x[:, 0]) <= (domain.w / 2)) * (np.abs(x[:, 1]) <= (domain.h / 2))]
+    mag_x = np.linalg.norm(x, axis=1)
+    x = x[np.argsort(mag_x)]  # Sort by magnitude
+
+    v = x[0]
+    x = x[1:]  # Remove v from search set.
     # Checking 0 < ax + by < v*v to make the sites are within the region.
-    tol = 1e-3
-    vdot = np.matmul(all_sites, v)
-    in_box = all_sites[np.where((-tol <= vdot) & (vdot <= (v.dot(v) + tol)))[0]]
-    in_box = np.expand_dims(in_box, 0).swapaxes(
-        0, 1
-    )  # Used for the next step, getting site*site
+    tol = 1e-8
+    vdot = np.matmul(x, v)
+    x = x[(-tol <= vdot) * (vdot <= v.dot(v) + tol)]
+    mag_x = np.linalg.norm(x, axis=1)
+    x = x[np.argsort(mag_x)]
 
-    v2 = in_box[
-        np.argmin(np.squeeze(np.matmul(in_box, in_box.transpose(0, 2, 1))))
-    ].flatten()
+    v2 = x[0]
 
-    if np.all(v == v2):
-        print(v, v2, n, w, h, config)
     return tuple(v), tuple(v2)
 
 
@@ -95,7 +119,8 @@ def avg_rp(d: float, l: float) -> float:
 def circumcenter(v: numpy.ndarray, w: numpy.ndarray) -> numpy.ndarray:
     det = 1 / (2 * rot(v).dot(w))
     if rot(v).dot(w) == 0:
-        print(v, w)
+        pass
+        # print(v, w)
     v2, w2 = v.dot(v), w.dot(w)
     c = np.empty((2,))
     c[0], c[1] = w[1] * v2 - v[1] * w2, -w[0] * v2 + v[0] * w2
@@ -106,3 +131,71 @@ def rot(v: numpy.ndarray) -> numpy.ndarray:
     w = np.copy(v)
     w[0], w[1] = -w[1], w[0]
     return w
+
+
+def divisors(n: int) -> List[int]:
+    divs = [[i, n // i] for i in range(1, int(sqrt(n)) + 1) if n % i == 0]
+    return sorted(set(list(sum(divs, []))))
+
+
+def factorize(n: int) -> Dict[int, int]:
+    primes = [i for i in range(1, n + 1) if len(divisors(i)) == 2]
+    prime_fac = {}
+    for prime in primes:
+        i = 0
+        while n % prime ** (i + 1) == 0:
+            i += 1
+        if i > 0:
+            prime_fac[prime] = i
+
+    return prime_fac
+
+
+def hexagon_alpha(n: int, fraction: bool = False) -> List[int]:
+    if n % 2 == 1:
+        return []
+
+    fac = factorize(n)
+    q = 1
+    for prime in fac:
+        if prime % 3 != 2:
+            q *= prime ** fac[prime]
+
+    divq = divisors(q)
+    ratios, thres = [], 1 / sqrt(3)
+    for g in divq:
+        us = n // (2 * g)
+        divu = divisors(us)
+        for u in divu:
+            d = 2 * g * u * u
+            f = Fraction(n, d)
+            if f <= thres:
+                ratios.append(f)
+            else:
+                ratios.append(Fraction(d, 3 * n))
+
+    ratios = sorted(set(ratios))
+    if fraction:
+        return ratios
+    else:
+        return [float(x) * sqrt(3) for x in ratios]
+
+
+def hexagon_alpha_brute(n: int):
+    w = cmath.rect(1, 2 * pi / 3)
+    divs = divisors(n / 2)
+    ratios, thres = [], 1 / sqrt(3)
+    for a in range(n):
+        for b in range(1, a + 1):
+            if a == 0 and b == 0:
+                continue
+            z2, g = a * a - a * b + b * b, gcd(a - 2 * b, 2 * a - b)
+            if z2 // g in divs:
+                f = Fraction(n, 2 * z2)
+                if f <= thres:
+                    ratios.append(f)
+                else:
+                    ratios.append(Fraction(2 * z2, 3 * n))
+
+    ratios = sorted(set([float(x) * sqrt(3) for x in ratios]))
+    return ratios

squish/simulation.py

@@ -59,11 +59,16 @@ class Simulation:
 
     def normalize(self) -> None:
         new_frames = []
+        first = True
         for frame in self.frames:
             aspect = frame.w / frame.h
+
             new_domain = DomainParams(
                 frame.n, sqrt(frame.n * aspect), sqrt(frame.n / aspect), frame.r
             )
+            if first:
+                self.domain = new_domain
+                first = False
 
             new_points = frame.site_arr * np.array(
                 [new_domain.w / frame.w, new_domain.h / frame.h]
@@ -159,6 +164,9 @@ class Simulation:
         sim, frames = Simulation.load(path)
         for frame in frames:
             sim.frames.append(sim.energy.mode(*frame["domain"], frame["arr"]))
+            for k, v in frame["stats"].items():
+                if k not in sim.frames[-1].stats:
+                    sim.frames[-1].stats[k] = v
 
         return sim
 
@@ -233,13 +241,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) / (self.domain.n ** 0.5)
+
             if self.adaptive:
                 error = change - grad * self.step_size
                 tol = 10 ** min(-3, -2 + log10(grad_norm))
-
+                # tol = 10 ** -10
                 self.step_size *= (tol / np.linalg.norm(error)) ** 0.5
 
             if not save:
@@ -315,7 +324,7 @@ class Search(Simulation):
         new_sites: Optional[numpy.ndarray] = None,
     ) -> None:
         if log:
-            print(f"Travel - {self.domain}", flush=True)
+            print(f"Search - {self.domain} - {len(self)}", flush=True)
         if save and len(self) == 0:
             self.save(self.initial_data, True)
 
@@ -328,16 +337,15 @@ class Search(Simulation):
             sim.run(False, log, log_steps)
 
             self.frames.append(sim[-1])
+            # Get Hessian,and check nullity. If > 2, perturb.
+            hess = self.frames[i].hessian
+            eigs = np.sort(np.linalg.eig(hess)[0])
+            self.frames[i].stats["eigs"] = eigs
+
             if save:
                 self.save(self.frame_data(i))
             if log:
                 print(f"Equilibrium: {i:04}\n", flush=True)
-
-            # Get Hessian,and check nullity. If > 2, perturb.
-            hess = self.frames[i].hessian(10e-5)
-            eigs = np.sort(np.linalg.eig(hess)[0])
-            self.frames[i].stats["eigs"] = eigs
-
             zero_eigs = np.count_nonzero(
                 np.isclose(eigs, np.zeros((len(eigs),)), atol=1e-4)
             )
@@ -386,10 +394,7 @@ class Shrink(Simulation):
     ) -> None:
         super().__init__(domain, energy, name=name)
         self.step_size, self.thres, self.adaptive = step_size, thres, adaptive
-        self.delta, self.stop_width = (
-            self.domain.w * delta / 100,
-            self.domain.w * stop_width,
-        )
+        self.delta, self.stop_width = delta, stop_width
 
     @property
     def initial_data(self) -> Dict:
@@ -418,8 +423,19 @@ class Shrink(Simulation):
             self.save(self.initial_data, True)
 
         width = self.domain.w if len(self.frames) == 0 else self.frames[-1].w
+        if len(self.frames) != 0:
+            new_sites = self.frames[-1].site_arr
+            width = self.frames[-1].w - self.delta
+        else:
+            width = self.domain.w
         i = 0
-        while width >= self.stop_width:
+
+        if delta < 0:
+            cond = width <= self.stop_width
+        if delta > 0:
+            cond = width <= self.stop_width
+
+        while cond:
             # Get to equilibrium.
             new_domain = DomainParams(
                 self.domain.n, width, self.domain.h, self.domain.r

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,

squish/voronoi.pxd

@@ -66,7 +66,7 @@ cdef class VoronoiContainer:
     cdef readonly INT_T n
     cdef readonly FLOAT_T w, h, r, energy
     cdef FLOAT_T [2] dim
-    cdef FLOAT_T [:, ::1] points, vertices, site_cache, edge_cache, grad
+    cdef FLOAT_T [:, ::1] points, vertices, site_cache, edge_cache, grad, hess
     cdef INT_T [:, ::1] sites, edges
     cdef EdgeCacheMap* edge_cache_map
     cdef dict __dict__
@@ -77,6 +77,7 @@ cdef class VoronoiContainer:
     cdef void common_cache(VoronoiContainer self) except *
     cdef void precompute(self) except *
     cdef void calc_grad(self) except *
+    cdef void calc_hess(self) except *
     cdef void get_statistics(VoronoiContainer self) except *
 
     @staticmethod

squish/voronoi.pyx

@@ -1,4 +1,4 @@
-import array, scipy.spatial, numpy as np
+import array, scipy.spatial, numpy as np, math
 from cython.parallel import parallel, prange
 
 cimport numpy as np
@@ -6,8 +6,8 @@ from cpython cimport array
 from libc.math cimport isnan, NAN, pi as PI
 
 from squish.core cimport INT_T, FLOAT_T, \
-    IArray, FArray, Vector2D, Matrix2x2, \
-    _IArray, _FArray, _Vector2D, _Matrix2x2
+    IArray, FArray, Vector2D, Matrix2x2, BitSet, \
+    _IArray, _FArray, _Vector2D, _Matrix2x2, _BitSet
 
 from squish.voronoi cimport SiteCacheMap, EdgeCacheMap, VoronoiInfo, Site, HalfEdge
 
@@ -554,8 +554,8 @@ cdef class VoronoiContainer:
                 # vp - vm, vm - xi
                 la, da = q.copy.vsub(&q, p), p.copy.vsub(&p, xi.vec(&xi))
                 la_mag = la.mag(&la)
-                area_p = la.dot(&la, da.rot(&da))
-                Rla = la.rot(&la)
+                area_p = la.dot(&la, da.copy.rot(&da))
+                Rla = la.copy.rot(&la)
                 ya = Rla.copy.smul(&Rla, -2*area_p/la.dot(&la, la))
 
                 # Calculating centroid.
@@ -591,14 +591,14 @@ cdef class VoronoiContainer:
         xj, xk = em.cache.ya(&em, NAN_VECTOR), ep.cache.ya(&ep, NAN_VECTOR)
         Rxjk = xk.copy.vsub(&xk, xj)
         Rxjk.self.smul(&Rxjk, 2)
-        Rxjk = Rxjk.rot(&Rxjk)
+        Rxjk.self.rot(&Rxjk)
 
         v = ep.cache.da(&ep, NAN_VECTOR)
         top = R.copy.smul(&R, xj.dot(&xj, xj) - xk.dot(&xk, xk))
         top.self.msub(&top, _Matrix2x2(v.x*Rxjk.x, v.x*Rxjk.y, v.y*Rxjk.x, v.y*Rxjk.y))
         top.self.sdiv(&top, -Rxjk.dot(&Rxjk, xj))
 
-        return _Matrix2x2(top.a, top.c, top.b, top.d)
+        return top
 
     @staticmethod
     cdef inline bint sign(FLOAT_T [::1] ref, FLOAT_T [::1] p, FLOAT_T [::1] q):
@@ -621,6 +621,9 @@ cdef class VoronoiContainer:
     cdef void calc_grad(self) except *:
         pass
 
+    cdef void calc_hess(self) except *:
+        pass
+
     cdef void get_statistics(self) except *:
         self.stats = {}
         cache = self.site_cache[:self.n, :]
@@ -647,6 +650,7 @@ cdef class VoronoiContainer:
         edge_cache = np.asarray(self.edge_cache)
 
         self.stats["edge_lengths"] = edge_cache[caches, self.edge_cache_map.ila_mag]
+        self.stats["site_distances"] = edge_cache[caches, self.edge_cache_map.iya_mag]
 
     @property
     def site_arr(self):
@@ -660,6 +664,11 @@ cdef class VoronoiContainer:
     def gradient(self):
         return np.asarray(self.grad, dtype=FLOAT)
 
+    @property
+    def hessian(self):
+        self.calc_hess()
+        return np.asarray(self.hess, dtype=FLOAT)
+
     def add_sites(self, add):
         return (self.site_arr + add) % np.asarray(self.dim, dtype=FLOAT)
 
@@ -671,146 +680,6 @@ cdef class VoronoiContainer:
 
         return -(step/2)*(k1+k2), -k1
 
-
-    def approx_hessian(self, d: float) -> np.ndarray:
-        """
-        Obtains the approximate Hessian.
-        :param d: [float] small d for approximation.
-        :return: 2Nx2N array that represents Hessian.
-        """
-        HE = np.zeros((2*self.n, 2*self.n))
-        new_sites = np.copy(self.site_arr)    # Maintain one copy for speed.
-        for i in range(self.n):
-            for j in range(2):
-                mod = self.w if j == 0 else self.h
-                new_sites[i][j] = (new_sites[i][j] + d) % mod
-                Ep = self.__class__(self.n, self.w, self.h, self.r, new_sites)
-                new_sites[i][j] = (new_sites[i][j] - 2*d) % mod
-                Em = self.__class__(self.n, self.w, self.h, self.r, new_sites)
-                new_sites[i][j] = (new_sites[i][j] + d) % mod
-
-                HE[:, 2*i+j] = ((Ep.gradient - Em.gradient)/(2*d)).flatten()
-
-        # Average out discrepencies, since it should be symmetric.
-        for i in range(2*self.n):
-            for j in range(i, 2*self.n):
-                HE[i][j] = (HE[i][j] + HE[j][i])/2
-                HE[j][i] = HE[i][j]
-
-        return HE
-
-    def radialt_hessian(self) -> np.ndarray:
-        HE = np.zeros((2*self.n, 2*self.n))
-        cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
-                                             self.vertices, self.site_cache,
-                                             self.edge_cache, self.edge_cache_map)
-        cdef Site xi, xk
-        cdef HalfEdge e, em, ep, fj, fk
-        cdef Vector2D temp1, temp2, dau, dapu
-        cdef Matrix2x2 sigI, sigJ, sigK, p, q, tempm, toI, toJ, toK
-
-        cdef INT_T i, j, k, z
-        for i in range(self.n):
-            xi = _Site(i, &info)
-            ep = xi.edge(&xi)
-            e = ep.prev(&ep)
-            j = 0
-            while j < xi.edge_num(&xi):
-                e.cache.H(&e, VoronoiContainer.calc_H(e, ep))
-                e = e.next(&e)
-                j = j + 1
-
-
-        for i in range(self.n):
-            xi = _Site(i, &info)
-            e = xi.edge(&xi)
-
-            z = 0
-            while z < xi.edge_num(&xi):
-                em, ep = e.prev(&e), e.next(&e)
-                fj, fk = em.twin(&em), e.twin(&e)
-                fj, fk = fj.next(&fj), fk.next(&fk)
-
-                xj, xk = fj.face(&fj), fk.face(&fk)
-                j, k = xj.index(&xj) % self.n, xk.index(&xk) % self.n
-                if k < 0:
-                    k = k + self.n
-                if j < 0:
-                    j = j + self.n
-
-                sigI = e.cache.H(&e, NAN_MATRIX)
-                sigJ = fj.cache.H(&fj, NAN_MATRIX)
-                sigK = fk.cache.H(&fk, NAN_MATRIX)
-
-                ### Calculating of p
-                temp1 = e.cache.la(&e, NAN_VECTOR)
-                temp1 = temp1.rot(&temp1)
-                temp1.self.sdiv(
-                    &temp1,
-                    e.cache.la_mag(&e, NAN) * e.cache.ya_mag(&e, NAN) / 2
-                )
-
-                dau = e.cache.da(&e, NAN_VECTOR)
-                dau.self.sdiv(&dau, e.cache.da_mag(&e, NAN))
-                dapu = ep.cache.da(&ep, NAN_VECTOR)
-                dapu.self.sdiv(&dapu, ep.cache.da_mag(&ep, NAN))
-
-                temp2 = dapu.copy.vsub(&dapu, dau)
-                temp2 = temp2.rot(&temp2)
-
-                p = _Matrix2x2(
-                    temp1.x*temp2.x, temp1.x*temp2.y,
-                    temp1.y*temp2.x, temp1.y*temp2.y
-                )
-
-                ### Calculating of q
-                temp2 = e.cache.la(&e, NAN_VECTOR)
-                temp1 = temp2.rot(&temp2)
-                q = _Matrix2x2(
-                    temp1.x*temp2.x, temp1.x*temp2.y,
-                    temp1.y*temp2.x, temp1.y*temp2.y
-                )
-                q.self.sdiv(&q, e.cache.la_mag(&e, NAN)**2)
-
-                q.self.msub(&q, R)
-                q.self.smul(
-                    &q,
-                    e.cache.calI(&e, NAN) / e.cache.la_mag(&e, NAN)
-                )
-
-                temp2 = em.cache.la(&em, NAN_VECTOR)
-                temp1 = temp2.rot(&temp2)
-                tempm = _Matrix2x2(
-                    temp1.x*temp2.x, temp1.x*temp2.y,
-                    temp1.y*temp2.x, temp1.y*temp2.y
-                )
-                tempm.self.sdiv(&tempm, em.cache.la_mag(&em, NAN)**2)
-
-                tempm = R.copy.msub(&R, tempm)
-                tempm.self.smul(
-                    &tempm,
-                    em.cache.calI(&em, NAN) / em.cache.la_mag(&em, NAN)
-                )
-
-                q.self.madd(&q, tempm)
-
-                # Calculating components that go to the respective sites
-                toI = sigI.copy.mmul(&sigI, p.copy.madd(&p, q))
-                toI.self.msub(&toI, p)
-
-                toJ = sigJ.copy.mmul(&sigJ, p.copy.madd(&p, q))
-                toK = sigK.copy.mmul(&sigK, p.copy.madd(&p, q))
-
-                HE[2*i: 2*(i+1), 2*i: 2*(i+1)] += np.array([[toI.a, toI.b], [toI.c, toI.d]])
-                HE[2*i: 2*(i+1), 2*j: 2*(j+1)] += np.array([[toJ.a, toJ.b], [toJ.c, toJ.d]])
-                HE[2*i: 2*(i+1), 2*k: 2*(k+1)] += np.array([[toK.a, toK.b], [toK.c, toK.d]])
-
-                e = e.next(&e)
-                z = z + 1
-
-
-        return -2*self.r*HE
-
     def site_vert_arr(self): # -> List[np.ndarray]
         cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
             self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
@@ -836,3 +705,43 @@ cdef class VoronoiContainer:
             site_verts.append(verts)
 
         return sites, site_verts
+
+    def second_near_neighbor(self):
+        cdef VoronoiInfo info = _VoronoiInfo(self.sites, self.edges, self.points,
+            self.vertices, self.site_cache, self.edge_cache, self.edge_cache_map)
+        
+        cdef INT_T i, j, k    
+        cdef Site xi, xj, xk
+        cdef HalfEdge e, f
+        cdef Vector2D v
+        
+        snn = []
+        for i in range(self.n):
+            xi = _Site(i, &info)
+            e = xi.edge(&xi)
+            first_neighbors = np.empty((xi.edge_num(&xi),), dtype=INT)
+            for j in range(xi.edge_num(&xi)):
+                f = e.twin(&e)
+                xj = f.face(&f)
+                first_neighbors[j] = xj.index(&xj) % self.n
+                e = e.next(&e)
+
+            second_neighbors = set()
+            for j in range(len(first_neighbors)):
+                xj = _Site(first_neighbors[j], &info)
+                e = xj.edge(&xj)
+                for k in range(xj.edge_num(&xj)):
+                    f = e.twin(&e)
+                    xk = f.face(&f)
+                    second_neighbors.add(xk.index(&xk) % self.n)
+                    e = e.next(&e)
+            
+            for j in range(len(first_neighbors)):
+                second_neighbors.remove(first_neighbors[j])
+            second_neighbors.remove(i)
+            
+            snn.append(list(second_neighbors))
+        
+        return snn
+            
+