diff --git a/scripts/coercion.py b/scripts/coercion.py new file mode 100644 index 0000000..4220327 --- /dev/null +++ b/scripts/coercion.py @@ -0,0 +1,92 @@ +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(.07, .12, .97, .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 = [] + store_data = {} + + 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[zero_ind-1]) + + 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(.07, .12, .97, .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.legend() + ax.set_xlabel("Width") + ax.set_ylabel("Smallest positive 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.") \ No newline at end of file diff --git a/scripts/convergence.py b/scripts/convergence.py index 08972fe..5b64521 100644 --- a/scripts/convergence.py +++ b/scripts/convergence.py @@ -5,6 +5,8 @@ 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 convergence graphs for a collection of simulations.") @@ -17,14 +19,11 @@ def main(): data = {} for file in Path(args.sims_path).iterdir(): - sim = Simulation.load(file / "data.squish") - sim_info = next(sim) - - step = sim_info["step_size"] - for frame in sim: + sim, frames = Simulation.load(file / "data.squish") + step = sim.step_size data[step] = {"times": [], "values": [], "diffs": []} - for i, frame_info in enumerate(sim): + for i, frame_info in enumerate(frames): data[step]["times"].append(step*i) data[step]["values"].append(np.linalg.norm(frame_info["arr"])) data[step]["diffs"].append(np.linalg.norm(all_info[-1]["arr"] - frame_info["arr"])) @@ -47,7 +46,8 @@ def main(): ax[1].set_xlabel("Time") ax[1].set_ylabel("L2 Norm of Difference") - fig.savefig("figures/Equilibrium Convergence.png") + fig.savefig(OUTPUT_DIR / "Equilibrium Convergence.png") + print(f"Wrote to {OUTPUT_DIR / 'Equilibrium Convergence.png'}") if __name__ == '__main__': diff --git a/scripts/defects.py b/scripts/defects.py new file mode 100644 index 0000000..a179f31 --- /dev/null +++ b/scripts/defects.py @@ -0,0 +1,88 @@ +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 = {} + + + for file in Path(args.sims_path).iterdir(): + sim, frames = Simulation.load(file / "data.squish") + 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 + + 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(.07, .12, .97, .9) + + fig.suptitle("Defects at 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.log10(ns), np.log10(corrected) + m, b = np.polyfit(x, y, 1) + + x2, y2 = x[14:], np.log10(defects[14:]) + m2, b2 = np.polyfit(x2, y2, 1) + + ax[1].plot(x, y, linestyle='dotted', color='C0') + ax[1].plot(x, np.log10(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$ 25: {m2:.5f}") + ax[1].grid(zorder=0) + ax[1].legend() + ax[1].set_xlabel("log10 N") + ax[1].set_ylabel("log10 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.") \ No newline at end of file diff --git a/scripts/heatmap.py b/scripts/heatmap.py new file mode 100644 index 0000000..22962eb --- /dev/null +++ b/scripts/heatmap.py @@ -0,0 +1,119 @@ +from __future__ import annotations +from typing import List, Tuple, Dict +import argparse, numpy as np, os, pickle +import matplotlib.pyplot as plt +from multiprocessing import Pool, cpu_count +from pathlib import Path + +from squish import Simulation, DomainParams +from squish.common import OUTPUT_DIR + + +def eq_file_process(file: Path) -> Tuple[float, List[float], List[float]]: + sim, frames = Simulation.load(file / 'data.squish') + + 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 / 'data.squish') + 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()) + sim, frames = Simulation.load(files[0] / 'data.squish') + + 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 for N={sim.domain.n}... |{"#"*hashes}{" "*(20-hashes)}|' + \ + f' {i+1}/{len(files)} simulations loaded.', flush=True, end='\r') + print(flush=True) + + widths = np.asarray(sorted(data["all"])) + domain = DomainParams(sim.domain.n, widths[-1], sim.domain.h, sim.domain.r) + return data, widths, domain + + +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() + + # 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(.07, .12, .97, .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.") \ No newline at end of file diff --git a/scripts/perturbations.py b/scripts/perturbations.py new file mode 100644 index 0000000..ec40ec6 --- /dev/null +++ b/scripts/perturbations.py @@ -0,0 +1,61 @@ +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 perturbation graphs for a collection of simulations.") + parser.add_argument('sims_path', metavar='path/to/data', + help="folder that contains simulations of perturbations from an equilibrium.") + parser.add_argument('end_path', metavar='path/to/equilbrium', + help="NumPy binary (.npy) file that contains the equilibrium to compare to.") + parser.add_argument('-q', '--quiet', dest='quiet', action='store_true', default=False, + help="suppress all normal output") + + args = parser.parse_args() + + end = np.load(args.end_path) + + data = {} + for file in Path(args.sims_path).iterdir(): + k = float(file.name.split('k')[-1]) + delta = 10**k + + sim, frames = Simulation.load(file / 'data.squish') + data[delta] = {"norm": [], "time": [], "k": k} + + for i, frame in enumerate(frames): + adjusted = frame["arr"] + (end[0] - frame["arr"][0]) + + data[delta]["norm"].append(np.linalg.norm(adjusted - end)) + data[delta]["time"].append(sim.step_size * i) + + fig, ax = plt.subplots(figsize=(12, 8)) + plt.subplots_adjust(.07, .12, .97, .9) + + for delta in sorted(data): + ax.plot(np.log10(np.array(data[delta]["time"])+1), np.log10(data[delta]["norm"]), + label=f"k = {data[delta]['k']}") + + fig.suptitle("Equilibrium Perturbations") + ax.grid(zorder=0) + #ax.set_xlim([0, 5]) + ax.legend() + ax.set_xlabel("Log Time") + ax.set_ylabel("Log L2 Norm of Difference") + + fig.savefig(OUTPUT_DIR / "Equilibrium Perturbations.png") + print(f"Wrote to {OUTPUT_DIR / 'Equilibrium Perturbations.png'}") + + +if __name__ == '__main__': + os.environ["QT_LOGGING_RULES"] = "*=false" + try: + main() + except KeyboardInterrupt: + print("Program terminated by user.") \ No newline at end of file diff --git a/scripts/width_diagrams.py b/scripts/width_diagrams.py index 4aa0114..61398bd 100644 --- a/scripts/width_diagrams.py +++ b/scripts/width_diagrams.py @@ -199,10 +199,12 @@ def main(): 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 - min_order - max_unorder_off = max_unorder - min_order - ax.plot(widths, min_order - min_order, color='C1') + 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')