From 71fb42eb21092bb84d6e007ce9da63418d99025d Mon Sep 17 00:00:00 2001 From: Kenneth Jao Date: Thu, 2 Dec 2021 22:40:49 -0500 Subject: [PATCH] Updated scripts --- scripts/aspect_diagrams.py | 424 +++++++++++++++++++++++++++++++++++++ scripts/defectenergy.py | 39 ++++ scripts/maxcenter.py | 202 ++++++++++++++++++ scripts/width_diagrams.py | 276 ------------------------ 4 files changed, 665 insertions(+), 276 deletions(-) create mode 100644 scripts/aspect_diagrams.py create mode 100644 scripts/defectenergy.py create mode 100644 scripts/maxcenter.py delete mode 100644 scripts/width_diagrams.py diff --git a/scripts/aspect_diagrams.py b/scripts/aspect_diagrams.py new file mode 100644 index 0000000..2ba5d59 --- /dev/null +++ b/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.") diff --git a/scripts/defectenergy.py b/scripts/defectenergy.py new file mode 100644 index 0000000..9cfdb9a --- /dev/null +++ b/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.") diff --git a/scripts/maxcenter.py b/scripts/maxcenter.py new file mode 100644 index 0000000..bf9a5da --- /dev/null +++ b/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.") diff --git a/scripts/width_diagrams.py b/scripts/width_diagrams.py deleted file mode 100644 index 40a7f09..0000000 --- a/scripts/width_diagrams.py +++ /dev/null @@ -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.")