diff --git a/scripts/aspect_diagrams.py b/scripts/aspect_diagrams.py new file mode 100644 index 0000000..910f9a1 --- /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 + domains.append( + DomainParams( + orig_domain.n, + math.sqrt(orig_domain.n * aspect), + math.sqrt(orig_domain.n / aspect), + orig_domain.r, + ) + ) + + # domains = [ + # DomainParams(orig_domain.n, w, orig_domain.h, orig_domain.r) for w in widths + # ] + + with Pool(cpu_count()) as pool: + energy_mins, energy_maxes, isoparam_mins, isoparam_maxes = {}, {}, {}, {} + for i, res in enumerate(pool.imap_unordered(order_process, domains)): + energy_mins[res[0]] = res[1] + energy_maxes[res[0]] = res[2] + isoparam_mins[res[0]] = res[3] + isoparam_maxes[res[0]] = res[4] + + hashes = int(21 * i / len(widths)) + print( + f'Generating at width {res[0]:.02f}... |{"#"*hashes}{" "*(20-hashes)}|' + + f" {i+1}/{len(widths)} completed.", + flush=True, + end="\r", + ) + + print(flush=True) + + data["energy_min"] = list([x[1] for x in sorted(energy_mins.items())]) + data["energy_max"] = list([x[1] for x in sorted(energy_maxes.items())]) + data["isoparam_min"] = list([x[1] for x in sorted(isoparam_mins.items())]) + data["isoparam_max"] = list([x[1] for x in sorted(isoparam_maxes.items())]) + + return data + + +def eq_file_process(file: Path) -> Tuple[float, List[float], List[float]]: + sim, frames = Simulation.load(file) + + alls = [] + for frame_info in frames: + alls.append( + [ + frame_info["energy"], + np.var(frame_info["stats"]["avg_radius"]) <= 1e-8, + np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6), + sum(frame_info["stats"]["site_energies"][: sim.domain.n]), + ] + ) + + sim, frames = Simulation.load(file) + sim.frames = list(frames) + counts = sim.get_distinct() + + distincts = [] + for j, frame_info in enumerate(sim.frames): + distincts.append( + [ + frame_info["energy"], + np.var(frame_info["stats"]["avg_radius"]) <= 1e-8, + np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6), + sum(frame_info["stats"]["site_energies"][: sim.domain.n]), + counts[j], + ] + ) + + return sim.domain.w / sim.domain.h, alls, distincts + + +def get_equilibria_data(filepath: Path) -> Tuple[Dict, numpy.ndarray, DomainParams]: + data = {"all": {}, "distinct": {}} + files = list(Path(filepath).iterdir()) + + with Pool(cpu_count()) as pool: + for i, res in enumerate(pool.imap_unordered(eq_file_process, files)): + data["all"][res[0]] = res[1] + data["distinct"][res[0]] = res[2] + + hashes = int(21 * i / len(files)) + print( + f'Loading simulations... |{"#"*hashes}{" "*(20-hashes)}|' + + f" {i+1}/{len(files)} simulations loaded.", + flush=True, + end="\r", + ) + print(flush=True) + + sim, frames = Simulation.load(files[0]) + widths = np.asarray(sorted(data["all"])) + domain = DomainParams(sim.domain.n, widths[-1], sim.domain.h, sim.domain.r) + return data, widths, domain + + +def axis_settings(ax, widths): + ax.grid(zorder=0) + ax.set_xticks([round(w, 2) for w in widths[::2]]) + ax.set_xticklabels([f"{round(w, 3):.2f}" for w in widths[::2]], rotation=90) + plt.subplots_adjust(0.07, 0.12, 0.97, 0.9) + + +def probability_of_disorder(data, widths, domain): + fig, ax = plt.subplots(figsize=(16, 8)) + all_disorder_count = [] + for width in widths: + equal_shape = list([c[1] for c in data["all"][width]]) + all_disorder_count.append( + 100 * equal_shape.count(False) / len(data["all"][width]) + ) + + ax.plot(widths, all_disorder_count) + axis_settings(ax, widths) + + ax.yaxis.set_major_formatter(mtick.PercentFormatter()) + ax.title.set_text(f"Probability of Disorder - N{domain.n}") + ax.set_xlabel("Aspect Ratio") + ax.set_ylabel("Disordered Equilibria") + boa_y_min = round(min(all_disorder_count) / 20) * 20 - 5 + ax.set_yticks(np.arange(boa_y_min, 100.01, 2.5)) + + return fig + + +def density_of_states(data, widths, domain): + fig, ax = plt.subplots(figsize=(16, 8)) + distinct_ordered, distinct_unordered = [], [] + for width in widths: + equal_shape = list([c[1] for c in data["distinct"][width]]) + distinct_ordered.append(equal_shape.count(True)) + distinct_unordered.append(equal_shape.count(False)) + + ax2 = ax.twinx() + ax.plot(widths, distinct_unordered, label="Unordered Equilibria", color="C0") + ax2.plot(widths, distinct_ordered, label="Ordered Equilibria", color="C1") + axis_settings(ax, widths) + plt.subplots_adjust(0.07, 0.12, 0.92, 0.9) + ax.title.set_text(f"Density of States - N{domain.n}") + ax.set_xlabel("Aspect Ratio") + ax.set_ylabel("Number of States (Disordered)", color="C0") + ax2.set_ylabel("Number of States (Ordered)", color="C1") + + dos_y_max_unorder = 1.05 * max(distinct_unordered) + dos_y_max_order = 1.05 * max(distinct_ordered) + ax.set_yticks(np.linspace(0, dos_y_max_unorder, 20).astype(int)) + # ax.set_yticks(np.arange(0, dos_y_max_unorder, round(dos_y_max_unorder/200, 1)*10)) + ax2.set_yticks(np.arange(0, dos_y_max_order)) + + return fig + + +def defect_density(data, widths, domain): + fig, ax = plt.subplots(figsize=(16, 8)) + + defects = [] + for width in widths: + defects.append( + sum([c[2] for c in data["all"][width] if not c[1]]) + / len(data["all"][width]) + ) + + ax.plot(widths, defects) + axis_settings(ax, widths) + ax.title.set_text(f"Average Defects - N{domain.n}") + ax.set_xlabel("Aspect Ratio") + ax.set_ylabel("Defects") + ax.set_yticks(np.arange(0, 1 + max(defects), 0.5)) + + return fig + + +def circle_isoparam(data, widths, order_data, domain): + fig, ax = plt.subplots(figsize=(16, 8)) + + ax2 = ax.twinx() + axis_settings(ax, widths) + plt.subplots_adjust(0.07, 0.12, 0.92, 0.9) + ax.title.set_text(f"Circular Isoparametric Ratio - N{domain.n}") + ax.set_xlabel("Aspect Ratio") + ax.set_ylabel("Maximum Ratio", color="C0") + ax2.set_ylabel("Minimum Ratio", color="C1") + + ax.plot(widths, order_data["isoparam_max"], label="Maximum", color="C0") + ax2.plot(widths, order_data["isoparam_min"], label="Minimum", color="C1") + + return fig + + +def reduced_energy(data, widths, order_data, domain): + fig, ax = plt.subplots(figsize=(16, 8)) + + ordered_energies, unordered_energies = [], [] + for width in widths: + ordered_energies.append([c[0] for c in data["distinct"][width] if c[1]]) + unordered_energies.append([c[0] for c in data["distinct"][width] if not c[1]]) + + for i in range(len(order_data["energy_min"])): + ordered_energies[i].append(order_data["energy_min"][i]) + ordered_energies[i].append(order_data["energy_max"][i]) + + min_order = np.asarray([min(width) for width in ordered_energies]) + max_order = np.asarray([max(width) for width in ordered_energies]) + min_unorder = np.asarray([min(width) for width in unordered_energies]) + max_unorder = np.asarray([max(width) for width in unordered_energies]) + + offset = np.array(min_order) + + min_unorder_off = min_unorder - offset + max_unorder_off = max_unorder - offset + ax.plot(widths, min_order - offset, color="C1") + # ax.plot(widths, max_order - offset, color='C1', linestyle='dotted') + ax.plot(widths, min_unorder_off, color="C0") + ax.plot(widths, max_unorder_off, color="C0", linestyle="dotted") + axis_settings(ax, widths) + + ax.title.set_text(f"Reduced Energy vs. Width - N{domain.n}") + ax.set_xlabel("Aspect Ratio") + ax.set_ylabel("Reduced Energy") + bif_y_max = np.max(np.abs(np.concatenate((min_unorder_off, max_unorder_off)))) + bif_top = np.arange( + 0, bif_y_max, round(bif_y_max / 20, -math.floor(math.log10(bif_y_max / 20))) + ) + ax.set_yticks(np.concatenate((-bif_top[1:][::-1], bif_top))) + + return fig + + +def defect_energy(data, widths, order_data, domain): + fig, ax = plt.subplots(figsize=(16, 8)) + + ordered_energies, unordered_energies = [], [] + for width in widths: + ordered_energies.append([c[0] for c in data["distinct"][width] if c[1]]) + unordered_energies.append([c[0] for c in data["distinct"][width] if not c[1]]) + + for i in range(len(order_data["energy_min"])): + ordered_energies[i].append(order_data["energy_min"][i]) + ordered_energies[i].append(order_data["energy_max"][i]) + + min_order = np.asarray([min(width) for width in ordered_energies]) + max_order = np.asarray([max(width) for width in ordered_energies]) + min_unorder = np.asarray([min(width) for width in unordered_energies]) + max_unorder = np.asarray([max(width) for width in unordered_energies]) + + offset = np.array(min_order) + + defect_a, defect_b = [], [] + for width in widths: + num_defects = [c[2] for c in data["all"][width]] + defect_energy = [c[3] for c in data["all"][width]] + m, b = np.polyfit(num_defects, defect_energy, 1) + + defect_a.append(m) + defect_b.append(b) + + ax2 = ax.twinx() + ax.plot(widths, defect_a, label="Energy per Defect", color="C0") + ax2.plot(widths, defect_b - offset, label="Relative Initial Energy", color="C1") + axis_settings(ax, widths) + plt.subplots_adjust(0.07, 0.12, 0.92, 0.9) + ax.title.set_text(f"Defect Energy - N{domain.n}") + ax.set_xlabel("Aspect Ratio") + ax.set_ylabel("Energy per Defect", color="C0") + ax2.set_ylabel("Relative Initial Energy", color="C1") + + return fig + + +def excess_energy(data, widths, order_data, domain): + fig, ax = plt.subplots(figsize=(16, 8)) + + ordered_energies, unordered_energies = [], [] + for width in widths: + ordered_energies.append([c[0] for c in data["distinct"][width] if c[1]]) + unordered_energies.append([c[0] for c in data["distinct"][width] if not c[1]]) + + for i in range(len(order_data["energy_min"])): + ordered_energies[i].append(order_data["energy_min"][i]) + ordered_energies[i].append(order_data["energy_max"][i]) + + min_order = np.asarray([min(width) for width in ordered_energies]) + max_order = np.asarray([max(width) for width in ordered_energies]) + min_unorder = np.asarray([min(width) for width in unordered_energies]) + max_unorder = np.asarray([max(width) for width in unordered_energies]) + + # Energy of regular hexagon with area 1 + offset = ( + 2 + - 2 * domain.r * (6 * 3 ** (-0.25) * math.sqrt(2) * math.atanh(0.5)) + + 2 * math.pi * domain.r ** 2 + ) + + min_order_off = min_order / domain.n - offset + min_unorder_off = min_unorder / domain.n - offset + max_unorder_off = max_unorder / domain.n - offset + + ax.plot(widths, min_order_off, color="C1", label="Minimum Ordered") + ax.plot(widths, min_unorder_off, color="C0", label="Minimum Disordered") + ax.plot( + widths, + max_unorder_off, + color="C0", + linestyle="dotted", + label="Maximum Disordered", + ) + # ax.plot( + # [min(widths), max(widths)], + # [offset, offset], + # color="C1", + # linestyle="dotted", + # label="Regular Energy", + # ) + + axis_settings(ax, widths) + ax.title.set_text(f"Energy at Aspect Ratios - N{domain.n}") + ax.set_xlabel("Aspect Ratio") + ax.set_ylabel("Excess Energy per Site") + ax.legend() + + start, end = ax.get_ylim() + ax.set_yticks(np.linspace(0, end, 20)) + ax.ticklabel_format(axis="y", style="sci") + + return fig + + +def main(): + # Loading arguments. + parser = argparse.ArgumentParser("Outputs width search data into diagrams") + parser.add_argument( + "sims_path", + metavar="path/to/data", + help="folder that contains simulation files, or cached data file.", + ) + parser.add_argument( + "-q", + "--quiet", + dest="quiet", + action="store_true", + default=False, + help="suppress all normal output", + ) + + args = parser.parse_args() + + # Obtain data from simulation files and generate single shape data. + data, widths, domain = get_equilibria_data(Path(args.sims_path)) + order_data = get_ordered_energies(domain, widths) + + fig_folder = OUTPUT_DIR / Path(f"AspectDiagrams - N{domain.n}") + fig_folder.mkdir(exist_ok=True) + + # Generating diagrams. + probability_of_disorder(data, widths, domain).savefig( + fig_folder / "Probability of Disorder.png" + ) + + density_of_states(data, widths, domain).savefig( + fig_folder / "Density Of States.png" + ) + + defect_density(data, widths, domain).savefig(fig_folder / "Defects.png") + + reduced_energy(data, widths, order_data, domain).savefig( + fig_folder / "Reduced Energy.png" + ) + + defect_energy(data, widths, order_data, domain).savefig( + fig_folder / "Defect Energy.png" + ) + + circle_isoparam(data, widths, order_data, domain).savefig( + fig_folder / "Circular Isoparametric Ratio.png" + ) + + excess_energy(data, widths, order_data, domain).savefig( + fig_folder / "Excess Energy.png" + ) + + print(f"Wrote to {fig_folder}.") + + +if __name__ == "__main__": + os.environ["QT_LOGGING_RULES"] = "*=false" + try: + main() + except KeyboardInterrupt: + print("Program terminated by user.") diff --git a/scripts/coercivity.py b/scripts/coercivity.py new file mode 100644 index 0000000..dc3ea41 --- /dev/null +++ b/scripts/coercivity.py @@ -0,0 +1,113 @@ +from __future__ import annotations +from typing import List, Tuple, Dict +import argparse, math, numpy as np, os, pickle +import matplotlib.pyplot as plt +import matplotlib.ticker as mtick +from multiprocessing import Pool, cpu_count +from pathlib import Path + +import squish.ordered as order +from squish import Simulation, DomainParams +from squish.common import Energy, OUTPUT_DIR + + +def axis_settings(ax, widths): + ax.invert_xaxis() + ax.grid(zorder=0) + ax.set_xticks([round(w, 2) for w in widths[::-2]]) + ax.set_xticklabels(ax.get_xticks(), rotation=90) + plt.subplots_adjust(0.07, 0.12, 0.97, 0.9) + + +def main(): + # Loading arguments. + parser = argparse.ArgumentParser("Outputs ordered equilibria lowest eigenvalues.") + parser.add_argument( + "n_objects", + metavar="N", + type=int, + help="folder that contains simulation files, or cached data file.", + ) + parser.add_argument( + "-q", + "--quiet", + dest="quiet", + action="store_true", + default=False, + help="suppress all normal output", + ) + + args = parser.parse_args() + + widths = np.round(np.linspace(3.0, 10.0, 141), 2) + + values = [] + with open("coercivity_eigs.pkl", "rb") as f: + store_data = pickle.load(f) + + for width in widths: + eig_vals = [] + for config, eigs in store_data[width].items(): + zero_ind = np.where(np.isclose(eigs, 0, atol=1e-4))[0][0] + if zero_ind == 0: + eig_vals.append(eigs[2]) + else: + eig_vals.append(eigs[0]) + + values.append(min(eig_vals)) + + # for i, width in enumerate(widths): + # domain = DomainParams(args.n_objects, width, 10, 4.0) + # eig_vals = [] + # store_data[width] = {} + # configs = order.configurations(domain) + # for j, config in enumerate(configs): + # if config == (1,0): + # continue + # points = order.sites(domain, config) + + # hess = Energy("radial-t").mode(*domain, points).hessian(10e-5) + # eigs = np.sort(np.linalg.eig(hess)[0])[::-1] + # store_data[width][config] = eigs + + # zero_ind = np.where(np.isclose(eigs, 0))[0][0] + # if zero_ind == 0: + # eig_vals.append(eigs[2]) + # else: + # eig_vals.append(eigs[0]) + + # hashes = int(21*j/len(widths)) + # print(f'Generating at {width}, {i+1}/{len(widths)}... |{"#"*hashes}{" "*(20-hashes)}|' + \ + # f' {j+1}/{len(configs)} configs done.', flush=True, end='\r') + + # print(flush=True) + + # with open("coercivity_eigs.pkl", "wb") as f: + # pickle.dump(store_data, f, pickle.HIGHEST_PROTOCOL) + + # return + fig, ax = plt.subplots(figsize=(12, 8)) + plt.subplots_adjust(0.07, 0.12, 0.97, 0.9) + + ax.plot(widths, values) + + ax.invert_xaxis() + ax.grid(zorder=0) + ax.set_xticks([round(w, 2) for w in widths[::-2]]) + ax.set_xticklabels(ax.get_xticks(), rotation=90) + + fig.suptitle("Coercivity") + # ax.set_xlim([0, 5]) + ax.set_xlabel("Width") + ax.set_ylabel("Eigenvalue") + + fig.savefig(OUTPUT_DIR / "Coercivity.png") + print(f"Wrote to {OUTPUT_DIR / 'Coercivity.png'}") + + +if __name__ == "__main__": + os.environ["QT_LOGGING_RULES"] = "*=false" + try: + main() + except KeyboardInterrupt: + print("Program terminated by user.") diff --git a/scripts/convergence.py b/scripts/convergence.py index 19667a0..f5cccea 100644 --- a/scripts/convergence.py +++ b/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 = {} diff --git a/scripts/cumulative_vee.py b/scripts/cumulative_vee.py index 16eb19a..242b86f 100644 --- a/scripts/cumulative_vee.py +++ b/scripts/cumulative_vee.py @@ -72,20 +72,20 @@ def main(): for i, count in enumerate(counts): my_cool_data[i + 1].append(count) - - ax.plot( - 100 * vees[: 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}", - ) + 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) diff --git a/scripts/cumulative_vee_2.py b/scripts/cumulative_vee_2.py index 7d59e87..e6e7314 100644 --- a/scripts/cumulative_vee_2.py +++ b/scripts/cumulative_vee_2.py @@ -167,8 +167,6 @@ def main(): pad = 100 alpha_prob_dens = {} for alpha, avg in alpha_avgs.items(): - if alpha == 1.0: - continue mov_avgs = np.array( [np.mean(avg[max(0, i - pad) : i + pad]) for i in range(len(avg))] ) @@ -183,18 +181,22 @@ def main(): d_avgs = np.gradient(mov_avgs, vees[1] - vees[0]) a1_prob_dens[n] = d_avgs / (np.sum(d_avgs) * (vees[1] - vees[0])) - + + print(list(alpha_prob_dens.keys())) for alpha, prob_dens in sorted(alpha_prob_dens.items()): if alpha in [0.3, 0.5, 0.7, 0.9, 1.0]: - continue + # continue ax.plot(100 * vees, prob_dens, label=fr"$\alpha={alpha:.1f}$") + ax.plot(100*vees, a1_prob_dens[395], zorder=50, label=r"$\alpha = 1.0$", color="black") + + # 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") + #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) 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/defects.py b/scripts/defects.py new file mode 100644 index 0000000..053948e --- /dev/null +++ b/scripts/defects.py @@ -0,0 +1,108 @@ +from __future__ import annotations +from typing import List +import argparse, pickle, numpy as np, os +from pathlib import Path +import matplotlib.pyplot as plt + +from squish import Simulation +from squish.common import OUTPUT_DIR + + +def main(): + parser = argparse.ArgumentParser("Graphs average defects at N.") + parser.add_argument( + "sims_path", + metavar="path/to/data", + help="folder that contains simulation files at various Ns.", + ) + parser.add_argument( + "-q", + "--quiet", + dest="quiet", + action="store_true", + default=False, + help="suppress all normal output", + ) + + args = parser.parse_args() + data = {} + + files = list(Path(args.sims_path).iterdir()) + for i, file in enumerate(files): + sim, frames = Simulation.load(file) + avg_defects = 0 + count = 0 + + for frame in frames: + if np.var(frame["stats"]["avg_radius"]) > 1e-8: + avg_defects += np.count_nonzero(frame["stats"]["site_edge_count"] != 6) + + count += 1 + + avg_defects /= 1 if count == 0 else count + data[sim.domain.n] = avg_defects + + hashes = int(21 * i / len(files)) + print( + f'Processed N={sim.domain.n:03} |{"#"*hashes}{" "*(20-hashes)}|' + + f" {i+1}/{len(files)} simulations processed.", + flush=True, + end="\r", + ) + + print(flush=True) + + data = sorted(data.items()) + ns, defects = np.array([x[0] for x in data]), np.array([x[1] for x in data]) + + corrected = [] + for i, x in enumerate(defects): + if x == 0: + corrected.append(defects[i + 1]) + else: + corrected.append(x) + + fig, ax = plt.subplots(1, 2, figsize=(16, 8)) + plt.subplots_adjust(0.07, 0.12, 0.97, 0.9) + + fig.suptitle("Defects vs. Number of Sites (N)") + + m0, b0 = np.polyfit(ns, defects, 1) + + ax[0].plot(ns, defects) + ax[0].plot(ns, m0 * ns + b0, label=f"Slope: {m0:.5f}") + ax[0].grid(zorder=0) + ax[0].legend() + ax[0].set_xlabel("N") + ax[0].set_ylabel("Average Defects") + + x, y = np.log(ns), np.log(corrected) + m, b = np.polyfit(x, y, 1) + + x2, y2 = x[40:], np.log(defects[40:]) + m2, b2 = np.polyfit(x2, y2, 1) + + ax[1].plot(x, y, linestyle="dotted", color="C0") + ax[1].plot(x, np.log(defects)) + + # ax[1].plot(x, m*x+b, label=f"All N: {m:.5f}") + ax[1].plot(x2, m2 * x2 + b2, label=f"N $\\geq$ 40: {m2:.5f}") + ax[1].grid(zorder=0) + ax[1].legend() + ax[1].set_xticks(np.arange(3.75, 6.25, 0.25)) + ax[1].set_yticks(np.arange(0, 4.5, 0.5)) + ax[1].set_xlim(3.75, 6) + ax[1].set_ylim(0, 4) + ax[1].set_xlabel("$\\ln$ N") + ax[1].set_ylabel("$\\ln$ Average Defects") + + fig.savefig(OUTPUT_DIR / "DefectsN.png") + print(f"Wrote to {OUTPUT_DIR / 'DefectsN.png'}") + + +if __name__ == "__main__": + os.environ["QT_log10GING_RULES"] = "*=false" + try: + main() + except KeyboardInterrupt: + print("Program terminated by user.") diff --git a/scripts/frustration_n.py b/scripts/frustration_n.py index c31f607..6903881 100644 --- a/scripts/frustration_n.py +++ b/scripts/frustration_n.py @@ -102,21 +102,21 @@ def main(): ) 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}$", - ) + #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, 37) + ax.set_ylim(3, 36) ax.set_yticks(np.arange(5, 40, 5)) # ax2.set_ylim(-3 * 0.4, 18 + 3 * 0.4) @@ -131,9 +131,9 @@ def main(): # labs = [l.get_label() for l in lns] ax.grid(zorder=0) # ax.legend(lns, labs, loc="lower right") - ax.legend() + #ax.legend() ax.set_xlabel("N") - # ax.set_ylabel(r"$\langle \mathrm{D} \rangle$") + ax.set_ylabel(r"$\langle \mathrm{D} \rangle$") # ax2.set_ylabel("VEE") # ax2.set_ylabel(r"$\zeta \left[\times 10^{4} \right]$", color="C1") # ax3.set_ylabel(r"Defect Energy $\left[\times 10^{3} \right]$", color="C2") diff --git a/scripts/heatmap.py b/scripts/heatmap.py new file mode 100644 index 0000000..2fd2b09 --- /dev/null +++ b/scripts/heatmap.py @@ -0,0 +1,302 @@ +from __future__ import annotations +from typing import List, Tuple, Dict +import argparse, math, numpy as np, os, pickle +import matplotlib.pyplot as plt +import matplotlib.ticker as mtick +from scipy.optimize import curve_fit +from multiprocessing import Pool, cpu_count +from pathlib import Path + +import squish.ordered as order +from squish import Simulation, DomainParams +from squish.common import OUTPUT_DIR + + +def order_process(domain: DomainParams) -> Tuple[float, float, float]: + energies, isoparams = [], [] + configs = order.configurations(domain) + for config in configs: + rbar = order.avg_radius(domain, config) + area = domain.w * domain.h / domain.n + + energies.append( + 2 * domain.w * domain.h + + 2 * math.pi * domain.n * (domain.r ** 2 - 2 * domain.r * rbar) + ) + + isoparams.append(math.pi * rbar ** 2 / area) + + return domain.w, min(energies), max(energies), min(isoparams), max(isoparams) + + +def get_ordered_energies(orig_domain: DomainParams, widths: np.ndarray) -> Dict: + data = {} + domains = [] + for w in widths: + aspect = w + domains.append( + DomainParams( + orig_domain.n, + math.sqrt(orig_domain.n * aspect), + math.sqrt(orig_domain.n / aspect), + orig_domain.r, + ) + ) + + # domains = [ + # DomainParams(orig_domain.n, w, orig_domain.h, orig_domain.r) for w in widths + # ] + + with Pool(cpu_count()) as pool: + energy_mins, energy_maxes, isoparam_mins, isoparam_maxes = {}, {}, {}, {} + for i, res in enumerate(pool.imap_unordered(order_process, domains)): + energy_mins[res[0]] = res[1] + energy_maxes[res[0]] = res[2] + isoparam_mins[res[0]] = res[3] + isoparam_maxes[res[0]] = res[4] + + hashes = int(21 * i / len(widths)) + print( + f'Generating at width {res[0]:.02f}... |{"#"*hashes}{" "*(20-hashes)}|' + + f" {i+1}/{len(widths)} completed.", + flush=True, + end="\r", + ) + + print(flush=True) + + data["energy_min"] = list([x[1] for x in sorted(energy_mins.items())]) + data["energy_max"] = list([x[1] for x in sorted(energy_maxes.items())]) + data["isoparam_min"] = list([x[1] for x in sorted(isoparam_mins.items())]) + data["isoparam_max"] = list([x[1] for x in sorted(isoparam_maxes.items())]) + + return data + + +def eq_file_process(file: Path) -> Tuple[float, List[float], List[float]]: + sim, frames = Simulation.load(file) + + alls = [] + for frame_info in frames: + alls.append( + [ + frame_info["energy"], + np.var(frame_info["stats"]["avg_radius"]) <= 1e-8, + np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6), + sum(frame_info["stats"]["site_energies"][: sim.domain.n]), + ] + ) + + sim, frames = Simulation.load(file) + sim.frames = list(frames) + counts = sim.get_distinct() + + distincts = [] + for j, frame_info in enumerate(sim.frames): + distincts.append( + [ + frame_info["energy"], + np.var(frame_info["stats"]["avg_radius"]) <= 1e-8, + np.count_nonzero(frame_info["stats"]["site_edge_count"] != 6), + sum(frame_info["stats"]["site_energies"][: sim.domain.n]), + counts[j], + ] + ) + + return sim.domain.w / sim.domain.h, alls, distincts + + +def get_equilibria_data(filepath: Path) -> Tuple[Dict, numpy.ndarray, DomainParams]: + data = {"all": {}, "distinct": {}} + files = list(Path(filepath).iterdir()) + + with Pool(cpu_count()) as pool: + for i, res in enumerate(pool.imap_unordered(eq_file_process, files)): + data["all"][res[0]] = res[1] + data["distinct"][res[0]] = res[2] + + hashes = int(21 * i / len(files)) + print( + f'Loading simulations... |{"#"*hashes}{" "*(20-hashes)}|' + + f" {i+1}/{len(files)} simulations loaded.", + flush=True, + end="\r", + ) + print(flush=True) + + sim, frames = Simulation.load(files[0]) + widths = np.asarray(sorted(data["all"])) + domain = DomainParams(sim.domain.n, widths[-1], sim.domain.h, sim.domain.r) + return data, widths, domain + + +def probability_of_disorder(data, widths, domain): + fig, ax = plt.subplots(figsize=(16, 8)) + all_disorder_count = [] + for width in widths: + equal_shape = list([c[1] for c in data["all"][width]]) + all_disorder_count.append( + 100 * equal_shape.count(False) / len(data["all"][width]) + ) + + return all_disorder_count + + +def excess_energy(data, widths, order_data, domain): + fig, ax = plt.subplots(figsize=(16, 8)) + + ordered_energies, unordered_energies = [], [] + for width in widths: + ordered_energies.append([c[0] for c in data["distinct"][width] if c[1]]) + unordered_energies.append([c[0] for c in data["distinct"][width] if not c[1]]) + + for i in range(len(order_data["energy_min"])): + ordered_energies[i].append(order_data["energy_min"][i]) + ordered_energies[i].append(order_data["energy_max"][i]) + + min_order = np.asarray([min(width) for width in ordered_energies]) + max_order = np.asarray([max(width) for width in ordered_energies]) + min_unorder = np.asarray([min(width) for width in unordered_energies]) + max_unorder = np.asarray([max(width) for width in unordered_energies]) + + return min_order - min_unorder + + +def sigmoid(x, x0, k): + return 100 / (1 + np.exp(-k * (x - x0))) + + +def main(): + # Loading arguments. + parser = argparse.ArgumentParser("Outputs width search data into diagrams") + parser.add_argument( + "sims_path", + metavar="path/to/data", + help="folder that contains simulation files of all searches for all N.", + ) + parser.add_argument( + "-q", + "--quiet", + dest="quiet", + action="store_true", + default=False, + help="suppress all normal output", + ) + + args = parser.parse_args() + + fig_folder = OUTPUT_DIR + fig_folder.mkdir(exist_ok=True) + + store = Path(args.sims_path) / "EEvsPoD.pkl" + + if store.is_file(): + with open(store, "rb") as f: + horiz, vert = pickle.load(f) + else: + horiz = [] + vert = [] + + for file in Path(args.sims_path).iterdir(): + # Obtain data from simulation files and generate single shape data. + data, widths, domain = get_equilibria_data(file) + order_data = get_ordered_energies(domain, widths) + + vert.append(probability_of_disorder(data, widths, domain)) + horiz.append(excess_energy(data, widths, order_data, domain)) + + horiz, vert = np.concatenate(horiz), np.concatenate(vert) + with open(store, "wb") as f: + pickle.dump((horiz, vert), f, pickle.HIGHEST_PROTOCOL) + + fig, ax = plt.subplots(figsize=(10, 10)) + + for i in range(2): + ax.scatter( + horiz[i * 141 : (i + 1) * 141], + vert[i * 141 : (i + 1) * 141], + alpha=0.5, + color=f"C{i}", + s=5, + ) + + start, end = ax.get_xlim() + + # popt, pcov = curve_fit(sigmoid, horiz, vert) + # x = np.linspace(start, end, 100) + # y = sigmoid(x, *popt) + # y = sigmoid(x, -1.35, 3) + # ax.plot(x, y, color="C1") + + plt.subplots_adjust(0.1, 0.1, 0.97, 0.93) + + ax.set_xticks(np.linspace(start, end, 10)) + ax.set_yticks(np.arange(0, 105, 5)) + ax.grid() + + ax.yaxis.set_major_formatter(mtick.PercentFormatter()) + ax.title.set_text("Excess Energy Difference vs. PoD") + ax.set_xlabel("Excess Energy Difference") + ax.set_ylabel("Probability of Disorder") + fig.savefig(OUTPUT_DIR / "Energy Diff and Probability") + return + + # with open("testing.pkl", "rb") as f: + # disorder_dict = pickle.load(f) + # widths = np.linspace(3.0, 10.0, 141) + # min_n, max_n = 60, 80 + + disorder_dict = {} + for file in Path(args.sims_path).iterdir(): + sim_data, widths, domain = get_equilibria_data(file) + + disorder_count = [] + for width in widths: + equal_shape = list([c[1] for c in sim_data["all"][width]]) + disorder_count.append( + 100 * equal_shape.count(False) / len(sim_data["all"][width]) + ) + + disorder_dict[domain.n] = disorder_count + + min_n, max_n = min(disorder_dict), max(disorder_dict) + filepath = f"Disorder Heatmap N{min_n}-{max_n}" + + # with open("testing.pkl", "wb") as f: + # pickle.dump(disorder_dict, f, pickle.HIGHEST_PROTOCOL) + + disorder_arr = np.zeros((max_n - min_n + 1, len(widths))) + for key, value in disorder_dict.items(): + disorder_arr[key - min_n] = np.asarray(value) + + fig, ax = plt.subplots(figsize=(12, 8)) + + extent = [min(widths), max(widths), min_n, max_n + 1] + ax.imshow( + disorder_arr, + cmap="plasma", + interpolation="nearest", + aspect="auto", + extent=extent, + ) + + ax.invert_xaxis() + ax.set_xticks([round(w, 2) for w in widths[::-2]]) + ax.set_xticklabels(ax.get_xticks(), rotation=90) + ax.set_yticks(list(range(min_n, max_n + 1))) + plt.subplots_adjust(0.07, 0.12, 0.97, 0.9) + + ax.title.set_text(filepath) + ax.set_xlabel("Width") + ax.set_ylabel("Number of Sites") + fig.savefig(OUTPUT_DIR / filepath) + + print(f"Wrote to {OUTPUT_DIR / filepath}.") + + +if __name__ == "__main__": + os.environ["QT_LOGGING_RULES"] = "*=false" + try: + main() + except KeyboardInterrupt: + print("Program terminated by user.") diff --git a/scripts/min_order_vee.py b/scripts/min_order_vee.py index e334cd7..229ba5a 100644 --- a/scripts/min_order_vee.py +++ b/scripts/min_order_vee.py @@ -1,5 +1,6 @@ import numpy as np, os, math import matplotlib.pyplot as plt +from scipy.optimize import curve_fit from multiprocessing import Pool, cpu_count from squish import Simulation, ordered, DomainParams @@ -64,17 +65,32 @@ def main(): data["Config"], ) - print(min_order, min_config) - fig = plt.figure(figsize=(20, 15)) gs = fig.add_gridspec(1, 1) ax = fig.add_subplot(gs[0]) - ax.scatter(ns, min_order) + pad = 5 + mov_avgs = np.array( + [np.mean(min_order[max(0, i - pad) : i + pad]) for i in range(len(min_order))] + ) + + ax.scatter(ns, mov_avgs, s=8) + + log_slope, _ = np.polyfit(np.log10(ns), np.log10(mov_avgs), 1) + print(log_slope) + + def g(x, k): + return k * (x ** log_slope) + + # return a * np.e ** (-b * (x - c)) + params, covs = curve_fit(g, ns, min_order, p0=[5000]) + ax.plot(np.arange(500, 5500), g(np.arange(500, 5500), *params), color="C1") + + print(params) ax.grid(zorder=0) - # ax.set_xlim(1000, 5000) + ax.set_xlim(800, 5200) ax.set_xlabel("N") ax.set_ylabel(r"VEE $\left[\times 10^{2}\right]$") diff --git a/scripts/near_neighbors.py b/scripts/near_neighbors.py new file mode 100644 index 0000000..7ca55e5 --- /dev/null +++ b/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() diff --git a/scripts/perturbations.py b/scripts/perturbations.py index ad83bb5..996dc3b 100644 --- a/scripts/perturbations.py +++ b/scripts/perturbations.py @@ -2,11 +2,14 @@ import argparse, numpy as np, os from pathlib import Path import matplotlib.pyplot as plt -from squish import DomainParams, Simulation, ordered -from squish.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(): @@ -14,9 +17,7 @@ def main(): description="Graphs perturbation graphs for a collection of simulations." ) parser.add_argument( - "sims_path", - metavar="sim_dir", - help="folder that contains of perturbations from an equilibrium", + "sim_store_path", metavar="sim_dir", help="folder to save simulations to" ) parser.add_argument( "end_path", @@ -32,41 +33,122 @@ def main(): ) args = parser.parse_args() - sims_path = Path(args.sims_path) + out_fol = Path(args.sim_store_path) - end_sim = Simulation.from_file(args.end_path) + (OUTPUT_DIR / out_fol).mkdir(exist_ok=True) + + end_eq = Simulation.from_file(args.end_path) + e_hex = ordered.e_hex(end_eq.domain) def f(): - data = {} - for file in sims_path.iterdir(): - if "k" not in file.name or file.is_file(): - continue - k = float(file.name.split("k")[-1]) - delta = 10 ** k + 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) - sim, frames = Simulation.load(file) - data[delta] = {"norm": [], "time": [], "k": k} + perturb = np.random.random_sample(end_eq.frames[0].site_arr.shape) + perturb /= np.linalg.norm(perturb) - for i, frame in enumerate(frames): - adjusted = frame["arr"] + ( - end_sim.frames[0].site_arr[0] - frame["arr"][0] - ) + 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) - data[delta]["norm"].append( - np.linalg.norm( - ordered.toroidal_distance( - end_sim.domain, adjusted, end_sim.frames[0].site_arr + 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]["time"].append(sim.step_size * i) + data[delta]["vee"].append(frame["energy"] / sim.domain.n - e_hex) - return data + k += 1 if k < 0 else 0.25 + m, _ = np.polyfit(data[delta]["time"], data[delta]["norm"], 1) + same_eq = m < 0 - data = get_data(sims_path / "PerturbData.pkl", f, regen=args.regen) + all_data.append({"vec": perturb, "data": data}) + + return all_data + + all_data = get_data(OUTPUT_DIR / out_fol / "PerturbData.pkl", f, regen=args.regen) + + end_vee = end_eq.frames[0].energy / end_eq.domain.n - e_hex + print(end_vee) + vees = [] + for dat in all_data: + for k, v in sorted(dat["data"].items()): + if v["norm"][-1] > 1e-3: + vees.append(dat["data"][k]["vee"][0]) + break + + eigs = np.sort(np.linalg.eigvalsh(end_eq.frames[0].hessian)) + + zero_ind = np.where(np.isclose(eigs, 0, atol=1e-8))[0] + if len(zero_ind) == 0: + coer = eigs[0] + elif zero_ind[0] == 0: + coer = eigs[2] + else: + coer = eigs[0] plt.rcParams.update(RC_SETTINGS) + fig = plt.figure(figsize=(15, 15)) + gs = fig.add_gridspec(1, 1) + ax = fig.add_subplot(gs[0]) + + ax.hist(vees, bins=np.linspace(0, 1, 30)) + ax.grid(zorder=0) + + props = dict(boxstyle="round", facecolor="white", alpha=0.8, zorder=20) + ax.text( + 0.60, + 0.96, + f"Min Escape = {min(vees):.6f}", + transform=ax.transAxes, + verticalalignment="top", + bbox=props, + ) + + ax.text( + 0.62, + 0.89, + f"Coercivity = {coer:.6f}", + transform=ax.transAxes, + verticalalignment="top", + bbox=props, + ) + + out = OUTPUT_DIR / f"{NAME2} - {args.sim_store_path}.png" + fig.savefig(out) + print(f"Wrote to {out}") + + return + fig = plt.figure(figsize=(30, 8)) gs = fig.add_gridspec(1, 1) ax = fig.add_subplot(gs[0]) diff --git a/scripts/torus_mesh.py b/scripts/torus_mesh.py index 6ca1e34..9830777 100644 --- a/scripts/torus_mesh.py +++ b/scripts/torus_mesh.py @@ -117,11 +117,6 @@ def torus_region(c, verts, height): v = np.hstack((v, height * np.ones((len(v), 1)))) return v, f, v_inds[sd - 1] - centers = np.array([[c[0], c[1], heights[i]]]) - v_top = np.hstack((verts, heights[i] * np.ones((len(verts), 1)))) - - face_verts = np.concatenate((centers, v_top)) - def main(): parser = argparse.ArgumentParser(description="Generates 3D model for equilibrium.") @@ -136,7 +131,7 @@ def main(): args = parser.parse_args() torus = args.torus - size = args.size + size = args.size * 10 # Get desired frame and load. sim, frames = Simulation.load(args.sims_path) diff --git a/scripts/vor_test.py b/scripts/vor_test.py new file mode 100644 index 0000000..0c0e1c1 --- /dev/null +++ b/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()