More scripts for post processing

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.")

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__':

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.")

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.")

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.")

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')