Updated scripts

scripts/aspect_diagrams.py

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

scripts/defectenergy.py

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

scripts/maxcenter.py

@@ -0,0 +1,202 @@
+from __future__ import annotations
+import matplotlib.pyplot as plt
+import os, numpy as np
+import cmath, math, pickle
+
+
+def main():
+
+    with open("site_verts.pkl", "rb") as f:
+        sites, site_verts = pickle.load(f)
+
+    for i in range(400):
+        verts = [
+            as_complex(site_verts[i][j] - sites[i]) for j in range(len(site_verts[i]))
+        ]
+        plot_2d(verts, f"squish_output/maxcenters_sim/{i:03}.png")
+    return
+    v = [
+        0.266 + 0.87j,
+        -0.626 + 0.747j,
+        -0.976 - 0.046j,
+        -0.283 - 0.873j,
+        0.676 - 0.447j,
+        0.875 + 0.414j,
+    ]
+
+    # v = [v[0], v[1], v[3]]
+
+    line = np.linspace(-1, 1, 120)
+    line2 = np.linspace(-1, 1, 120)
+    X, Y = np.meshgrid(line, line2)
+    Z = np.empty(X.shape)
+    DZ = np.empty(X.shape, dtype="complex")
+    HZ = np.empty((X.shape[0], X.shape[1], 2, 2))
+    for i, x in enumerate(line):
+        for j, y in enumerate(line2):
+            rad, deriv, hess = average_radius(x + 1j * y, v, l)
+            Z[j][i] = rad
+            DZ[j][i] = deriv
+            HZ[j][i] = hess
+
+    max_indices = np.unravel_index(np.argmax(Z), Z.shape)
+
+    fig = plt.figure()
+    ax1 = fig.add_subplot(111, projection="3d")
+    ax1.contour(
+        X,
+        Y,
+        Z,
+        np.linspace(4, 5.7, 15),
+        rstride=1,
+        cstride=1,
+        cmap="viridis",
+        edgecolor="none",
+    )
+    ax1.scatter(X[max_indices], Y[max_indices], Z[max_indices])
+
+    cent = centroid(v, l)
+    maxcent = maxcenter(v, l)
+    ax1.scatter(cent.real, cent.imag, 3)
+    ax1.scatter(maxcent.real, maxcent.imag, 3)
+
+    print(maxcent)
+    print(abs(maxcent - cent))
+
+    ax1.view_init(elev=90, azim=270)
+    plt.show()
+
+    ax2 = fig.add_subplot(111, projection="3d")
+    ax2.contour(
+        X,
+        Y,
+        DZ.real,
+        np.linspace(-3, 3, 9),
+        rstride=1,
+        cstride=1,
+        cmap="viridis",
+        edgecolor="none",
+    )
+    ax2.contour(
+        X,
+        Y,
+        DZ.imag,
+        np.linspace(-3, 3, 9),
+        rstride=1,
+        cstride=1,
+        cmap="viridis",
+        edgecolor="none",
+    )
+
+    ax2.view_init(elev=90, azim=270)
+    ax2.scatter(X[max_indices], Y[max_indices], Z[max_indices])
+    # ax2.plot_surface(X, Y, DZy, rstride=1, cstride=1, cmap="viridis", edgecolor="none")
+    # for vert in v:
+    #    ax.scatter(vert.real, vert.imag, 5)
+
+    plt.savefig("TestPolygonAverageRadius.png")
+    plt.show()
+
+
+def plot_2d(v: List[complex], name: str):
+    l = get_l(v)
+    cent, maxcent = centroid(v, l), maxcenter(v, l)
+
+    fig, ax = plt.subplots(1, figsize=(8, 8))
+
+    for i in range(len(v)):
+        va, vap = v[i], v[(i + 1) % len(v)]
+        ax.plot([va.real, vap.real], [va.imag, vap.imag], color="black")
+
+    ax.scatter(cent.real, cent.imag, label="Centroid")
+    ax.scatter(maxcent.real, maxcent.imag, label="Maxcenter")
+    ax.legend()
+    ax.grid()
+
+    plt.savefig(name)
+    plt.close()
+
+
+def get_l(v):
+    l = []
+    for i in range(len(v)):
+        l.append(v[(i + 1) % len(v)] - v[i])
+    return l
+
+
+def generate_hexagon():
+    angles = np.sort(np.random.random_sample((6,)))
+    while np.any(np.diff(angles) >= 0.5):
+        angles = np.sort(np.random.random_sample((6,)))
+    angles *= 2 * math.pi
+
+    mags = np.random.random_sample((3,))
+    mags = np.array([mags[0], 1, mags[1], 1, mags[2], 1])
+
+    v = []
+    for mag, angle in zip(mags, angles):
+        v.append(cmath.rect(1, angle))
+
+    return v
+
+
+def centroid(v, l):
+    area, cent = 0, 0
+    for i in range(len(v)):
+        jdi = v[i] * 1j
+        A = (jdi.conjugate() * l[i] + jdi * l[i].conjugate()) / 4
+        area += A
+        cent += (2 * v[i] + l[i]) * A
+
+    return (1 / (3 * area)) * cent
+
+
+def average_radius(x, v, l):
+    radius, deriv, hess = [], [], []
+    for i in range(len(v)):
+        jdi = (v[i] - x) * 1j
+        A = (jdi.conjugate() * l[i] + jdi * l[i].conjugate()) / 2
+        k = -1j * l[i]
+
+        da, dap = v[i] - x, v[i] - x + l[i]
+        dau, dapu = da / abs(da), dap / abs(dap)
+        kcu = k.conjugate() / abs(k)
+        z, zp = kcu * dau, kcu * dapu
+
+        int_rad = 2 * (cmath.atan(zp) - cmath.atan(z)) / (1j * abs(k))
+
+        radius.append(A * int_rad)
+        deriv.append(-k * int_rad)
+        hess.append(np.dot(as_vector(k).T, as_vector(1j * (dapu - dau) / A)))
+
+    if True in [x.real < 0 for x in radius]:
+        return 0, 0, 0
+    else:
+        return sum(radius).real, sum(deriv), sum(hess)
+
+
+def as_vector(c):
+    return np.atleast_2d(np.array([c.real, c.imag]))
+
+
+def as_complex(v):
+    v = v.flatten()
+    return v[0] + 1j * v[1]
+
+
+def maxcenter(v, l, delta=1e-8):
+    above_thres = True
+    x = centroid(v, l)
+    while above_thres:
+        rad, deriv, hess = average_radius(x, v, l)
+        above_thres = np.linalg.norm(deriv) > delta
+        x -= as_complex(as_vector(deriv).dot(np.linalg.inv(hess)))
+    return x
+
+
+if __name__ == "__main__":
+    os.environ["QT_LOGGING_RULES"] = "*=false"
+    try:
+        main()
+    except KeyboardInterrupt:
+        print("Program terminated by user.")

scripts/width_diagrams.py

@@ -1,276 +0,0 @@
-from __future__ import annotations
-from typing import List, Tuple, Dict
-import argparse, math, numpy as np, os
-import matplotlib.pyplot as plt
-import matplotlib.ticker as mtick
-from multiprocessing import Pool, cpu_count
-from pathlib import Path
-
-import squish.ordered as order
-from squish import Simulation, DomainParams
-from squish.common import OUTPUT_DIR
-
-
-def order_process(domain: DomainParams) -> Tuple[float, float, float]:
-    energies = []
-    configs = order.configurations(domain)
-    for config in configs:
-        energies.append(
-            2 * domain.w * domain.h
-            + 2
-            * math.pi
-            * domain.n
-            * (domain.r ** 2 - 2 * domain.r * order.avg_radius(domain, config))
-        )
-
-    return domain.w, min(energies), max(energies)
-
-
-def get_ordered_energies(orig_domain: DomainParams, widths: np.ndarray) -> Dict:
-    data = {}
-    domains = [
-        DomainParams(orig_domain.n, w, orig_domain.h, orig_domain.r) for w in widths
-    ]
-
-    with Pool(cpu_count()) as pool:
-        mins, maxes = {}, {}
-        for i, res in enumerate(pool.imap_unordered(order_process, domains)):
-            mins[res[0]] = res[1]
-            maxes[res[0]] = res[2]
-
-            hashes = int(21 * i / len(widths))
-            print(
-                f'Generating at width {res[0]:.02f}... |{"#"*hashes}{" "*(20-hashes)}|'
-                + f" {i+1}/{len(widths)} completed.",
-                flush=True,
-                end="\r",
-            )
-
-        print(flush=True)
-
-        data["min"] = list([x[1] for x in sorted(mins.items())])
-        data["max"] = list([x[1] for x in sorted(maxes.items())])
-
-    return data
-
-
-def eq_file_process(file: Path) -> Tuple[float, List[float], List[float]]:
-    sim, frames = Simulation.load(file)
-
-    alls = []
-    for frame_info in frames:
-        alls.append(
-            [
-                frame_info["energy"],
-                np.var(frame_info["stats"]["avg_radius"]) <= 1e-8,
-                np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6),
-            ]
-        )
-
-    sim, frames = Simulation.load(file)
-    sim.frames = list(frames)
-    counts = sim.get_distinct()
-
-    distincts = []
-    for j, frame_info in enumerate(sim.frames):
-        distincts.append(
-            [
-                frame_info["energy"],
-                np.var(frame_info["stats"]["avg_radius"]) <= 1e-8,
-                np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6),
-                counts[j],
-            ]
-        )
-
-    return sim.domain.w, alls, distincts
-
-
-def get_equilibria_data(filepath: Path) -> Tuple[Dict, numpy.ndarray, DomainParams]:
-    data = {"all": {}, "distinct": {}}
-    files = list(Path(filepath).iterdir())
-
-    with Pool(cpu_count()) as pool:
-        for i, res in enumerate(pool.imap_unordered(eq_file_process, files)):
-            data["all"][res[0]] = res[1]
-            data["distinct"][res[0]] = res[2]
-
-            hashes = int(21 * i / len(files))
-            print(
-                f'Loading simulations... |{"#"*hashes}{" "*(20-hashes)}|'
-                + f" {i+1}/{len(files)} simulations loaded.",
-                flush=True,
-                end="\r",
-            )
-        print(flush=True)
-
-    sim, frames = Simulation.load(files[0])
-    widths = np.asarray(sorted(data["all"]))
-    domain = DomainParams(sim.domain.n, widths[-1], sim.domain.h, sim.domain.r)
-    return data, widths, domain
-
-
-def axis_settings(ax, widths):
-    ax.invert_xaxis()
-    ax.grid(zorder=0)
-    ax.set_xticks([round(w, 2) for w in widths[::-2]])
-    ax.set_xticklabels(ax.get_xticks(), rotation=90)
-    plt.subplots_adjust(0.07, 0.12, 0.97, 0.9)
-
-
-def main():
-    # Loading arguments.
-    parser = argparse.ArgumentParser("Outputs width search data into diagrams")
-    parser.add_argument(
-        "sims_path",
-        metavar="path/to/data",
-        help="folder that contains simulation files, or cached data file.",
-    )
-    parser.add_argument(
-        "-q",
-        "--quiet",
-        dest="quiet",
-        action="store_true",
-        default=False,
-        help="suppress all normal output",
-    )
-
-    args = parser.parse_args()
-
-    data, widths, domain = get_equilibria_data(Path(args.sims_path))
-    order_data = get_ordered_energies(domain, widths)
-
-    fig_folder = OUTPUT_DIR / Path(f"ShrinkEnergyComparison - N{domain.n}")
-    fig_folder.mkdir(exist_ok=True)
-
-    # Torus minimum energies used as reference.
-
-    # Probability of disorder diagram.
-    fig, ax = plt.subplots(figsize=(16, 8))
-    all_disorder_count = []
-    for width in widths:
-        equal_shape = list([c[1] for c in data["all"][width]])
-        all_disorder_count.append(
-            100 * equal_shape.count(False) / len(data["all"][width])
-        )
-
-    ax.plot(widths, all_disorder_count)
-    axis_settings(ax, widths)
-
-    with open("N83-prob.txt", "w") as f:
-        f.write(", ".join([str(x) for x in widths]) + "\n")
-        f.write(", ".join([str(x) for x in all_disorder_count]))
-
-    ax.yaxis.set_major_formatter(mtick.PercentFormatter())
-    ax.title.set_text(f"Probability of Disorder - N{domain.n}")
-    ax.set_xlabel("Width")
-    ax.set_ylabel("Disordered Equilibria")
-    boa_y_min = round(min(all_disorder_count) / 20) * 20 - 5
-    ax.set_yticks(np.arange(boa_y_min, 100.01, 2.5))
-    fig.savefig(fig_folder / "Probability of Disorder.png")
-
-    # Density of States diagram.
-    fig, ax = plt.subplots(figsize=(16, 8))
-    distinct_ordered, distinct_unordered = [], []
-    for width in widths:
-        equal_shape = list([c[1] for c in data["distinct"][width]])
-        distinct_ordered.append(equal_shape.count(True))
-        distinct_unordered.append(equal_shape.count(False))
-
-    ax2 = ax.twinx()
-    ax.plot(widths, distinct_unordered, label="Unordered Equilibria", color="C0")
-    ax2.plot(widths, distinct_ordered, label="Ordered Equilibria", color="C1")
-    axis_settings(ax, widths)
-    ax.title.set_text(f"Density of States - N{domain.n}")
-    ax.set_xlabel("Width")
-    ax.set_ylabel("Number of States (Disordered)", color="C0")
-    ax2.set_ylabel("Number of States (Ordered)", color="C1")
-
-    dos_y_max_unorder = 1.05 * max(distinct_unordered)
-    dos_y_max_order = 1.05 * max(distinct_ordered)
-    ax.set_yticks(np.linspace(0, dos_y_max_unorder, 20).astype(int))
-    # ax.set_yticks(np.arange(0, dos_y_max_unorder, round(dos_y_max_unorder/200, 1)*10))
-    ax2.set_yticks(np.arange(0, dos_y_max_order))
-
-    fig.savefig(fig_folder / "Density Of States.png")
-
-    # Defect density diagram
-    fig, ax = plt.subplots(figsize=(16, 8))
-
-    defects = []
-    for width in widths:
-        defects.append(
-            sum([c[2] for c in data["all"][width] if not c[1]])
-            / len(data["all"][width])
-        )
-
-    ax.plot(widths, defects)
-    axis_settings(ax, widths)
-    ax.title.set_text(f"Average Defects - N{domain.n}")
-    ax.set_xlabel("Width")
-    ax.set_ylabel("Defects")
-    ax.set_yticks(np.arange(0, 1 + max(defects), 0.5))
-    fig.savefig(fig_folder / "Defects.png")
-
-    # Bifurcation diagram
-    fig, ax = plt.subplots(figsize=(16, 8))
-
-    ordered_energies, unordered_energies = [], []
-    for width in widths:
-        ordered_energies.append([c[0] for c in data["distinct"][width] if c[1]])
-        unordered_energies.append([c[0] for c in data["distinct"][width] if not c[1]])
-
-    for i in range(len(order_data["min"])):
-        ordered_energies[i].append(order_data["min"][i])
-        ordered_energies[i].append(order_data["max"][i])
-
-    null_unorder = []
-    for i, energies in enumerate(unordered_energies):
-        if len(energies) == 0:
-            null_unorder.append(i)
-            energies.append(order_data["min"][i])
-
-    min_order = np.asarray([min(width) for width in ordered_energies])
-    max_order = np.asarray([max(width) for width in ordered_energies])
-    min_unorder = np.asarray([min(width) for width in unordered_energies])
-    max_unorder = np.asarray([max(width) for width in unordered_energies])
-
-    offset = np.array(order_data["min"])
-    # offset = np.array(min_order)
-
-    min_unorder_off = min_unorder - offset
-    max_unorder_off = max_unorder - offset
-    ax.plot(widths, min_order - offset, color="C1")
-    # ax.plot(widths, max_order - offset, color='C1', linestyle='dotted')
-    ax.plot(widths, min_unorder_off, color="C0")
-    ax.plot(widths, max_unorder_off, color="C0", linestyle="dotted")
-    axis_settings(ax, widths)
-
-    with open("N83-od.txt", "w") as f:
-        f.write(", ".join([str(x) for x in widths]) + "\n")
-        f.write(", ".join([str(x) for x in min_unorder_off]) + "\n")
-        f.write(", ".join([str(x) for x in max_unorder_off]))
-
-    for i in null_unorder:
-        ax.scatter(widths[i], 0, marker="X", color="blue", s=50, zorder=4)
-        # ax.scatter(widths[i], max_unorder[i] - offset[i],
-        #     marker='X', edgecolors="blue", facecolors='none', s=100, zorder=4)
-
-    ax.title.set_text(f"Reduced Energy vs. Width - N{domain.n}")
-    ax.set_xlabel("Width")
-    ax.set_ylabel("Reduced Energy")
-    bif_y_max = np.max(np.abs(np.concatenate((min_unorder_off, max_unorder_off))))
-    bif_top = np.arange(
-        0, bif_y_max, round(bif_y_max / 20, -math.floor(math.log10(bif_y_max / 20)))
-    )
-    ax.set_yticks(np.concatenate((-bif_top[1:][::-1], bif_top)))
-    fig.savefig(fig_folder / "Bifurcation.png")
-
-    print(f"Wrote to {fig_folder}.")
-
-
-if __name__ == "__main__":
-    os.environ["QT_LOGGING_RULES"] = "*=false"
-    try:
-        main()
-    except KeyboardInterrupt:
-        print("Program terminated by user.")