Added more detail to shrink diagrams and fixed unique equilibria detection for ordered equilibria.

shrink_energy_comparison.py

@@ -1,7 +1,8 @@
 from __future__ import annotations
 from typing import List
-import os, argparse, numpy as np, pickle
+import os, math, argparse, numpy as np, pickle
 import matplotlib.pyplot as plt
+import matplotlib.ticker as mtick
 from pathlib import Path
 
 from simulation import Diagram, Simulation
@@ -59,38 +60,31 @@ def get_equilibria_data(filepath: Path):
 		with open(filepath, "rb") as data:
 			return pickle.load(data)
 
-	sims = []
+	data = {"all": {}, "distinct": {}}
 	files = list(Path(filepath).iterdir())
 
 	for i, file in enumerate(files):
-		sims.append(Simulation.load(file))
+		sim = Simulation.load(file)
+		data["all"][sim.w] = []
+		for frame in sim.frames:
+			data["all"][sim.w].append([frame.energy, np.var(frame.stats["avg_radius"]) <= 1e-8])
+
+		sim.get_distinct()
+		data["distinct"][sim.w] = []
+		for frame in sim.frames:
+			data["distinct"][sim.w].append([frame.energy,
+							np.var(frame.stats["avg_radius"]) <= 1e-8])
 
 		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)
-	sims.sort(key=lambda x: x.w)
 
-	widths = np.asarray([sim.w for sim in sims])
+	widths = np.asarray(sorted(data["all"]))
-
+	n, h, r, energy = sim.n, sim.h, sim.r, sim.energy
-	pairs_at_widths = []
-	for sim in sims:
-		pairs_at_widths.append([(frame.energy, np.var(frame.stats["avg_radius"]) <= 1e-8) \
-										 for frame in sim.frames])
-
-	# min_frames = [min(sim.frames, key=lambda x: x.energy) for sim in sims]
-	# max_frames = [max(sim.frames, key=lambda x: x.energy) for sim in sims]
-
-	# min_energies = np.asarray([frame.energy for frame in min_frames])
-	# max_energies = np.asarray([frame.energy for frame in max_frames])
-
-	# min_markers = [np.var(frame.stats["avg_radius"]) <= 1e-8 for frame in min_frames]
-	# max_markers = [np.var(frame.stats["avg_radius"]) <= 1e-8 for frame in max_frames]
-
-	n, h, r, energy = sims[0].n, sims[0].h, sims[0].r, sims[0].energy
 
 	sim_file = SIM_FOLDER / f"{ENERGY_R_STR[energy]} - EquilibriaData - N{n}.data"
-	out_tup = (widths, pairs_at_widths, n, h, r, ENERGY_I_STR[energy])
+	out_tup = (widths, data, n, h, r, ENERGY_I_STR[energy])
 	with open(sim_file, "wb") as output:
 		pickle.dump(out_tup, output)
 
@@ -101,6 +95,7 @@ def axis_settings(ax, widths):
 	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(.05, .12, .97, .9)
 
 
 def main():
@@ -113,7 +108,7 @@ def main():
 
 	args = parser.parse_args()
 
-	widths, energy_shape_tups, n, h, r, energy = get_equilibria_data(Path(args.sims_path))
+	widths, data, n, h, r, energy = get_equilibria_data(Path(args.sims_path))
 
 	torus_min_energies, torus_max_energies, _, _ = get_torus_config_energies(
 		n, widths, h, r, energy
@@ -123,13 +118,29 @@ def main():
 	fig_folder.mkdir(exist_ok=True)
 
 	# Torus minimum energies used as reference.
-	plt.tight_layout()
+
+	# Basin of attraction 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)
+	ax.yaxis.set_major_formatter(mtick.PercentFormatter())
+	ax.title.set_text('Basin of Attraction')
+	ax.set_xlabel("Width")
+	ax.set_ylabel("Disordered Equilibria")
+	ax.set_yticks(np.arange(0,105, 5))
+	fig.savefig(fig_folder / "Basin of Attraction.png")
+
+
 	# Density of States diagram.
 	fig, ax = plt.subplots(figsize=(16, 8))
-
 	distinct_ordered, distinct_unordered = [], []
-	for energy_shapes in energy_shape_tups:
+	for width in widths:
-		equal_shape = list([tup[1] for tup in energy_shapes])
+		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))
 
@@ -140,22 +151,18 @@ def main():
 	ax.title.set_text('Density of States')
 	ax.set_xlabel("Width")
 	ax.set_ylabel("Number of States")
+	dos_y_max = 1.05*max(distinct_ordered + distinct_unordered)
+	ax.set_yticks(np.arange(0, dos_y_max, round(dos_y_max/200, 1)*10))
 	fig.savefig(fig_folder / "Density Of States.png")
 
+
 	# Bifurcation diagram
 	fig, ax = plt.subplots(figsize=(16, 8))
 
 	ordered_energies, unordered_energies = [], []
-	for energy_shapes in energy_shape_tups:
+	for width in widths:
-		order_width, unorder_width = [], []
+		ordered_energies.append([c[0] for c in data["distinct"][width] if c[1]])
-		for pair in energy_shapes:
+		unordered_energies.append([c[0] for c in data["distinct"][width] if not c[1]])
-			if pair[1]:
-				order_width.append(pair[0])
-			else:
-				unorder_width.append(pair[0])
-		ordered_energies.append(order_width)
-		unordered_energies.append(unorder_width)
-
 
 	for i in range(len(torus_min_energies)):
 		ordered_energies[i].append(torus_min_energies[i])
@@ -163,30 +170,29 @@ def main():
 
 
 	null_unorder = []
-	for i, width in enumerate(unordered_energies):
+	for i, energies in enumerate(unordered_energies):
-		if len(width) == 0:
+		if len(energies) == 0:
 			null_unorder.append(i)
-			for x in unordered_energies[i-1]:
+			energies.append(torus_min_energies[i])
-				width.append(x)
-			#width = unordered_energies[i-1]
-			#width.append(max(unordered_energies[i-1]))
 
 	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])
 
+	min_unorder_off = min_unorder - torus_min_energies
+	max_unorder_off = max_unorder - torus_min_energies
 	ax.plot(widths, min_order - torus_min_energies, color='C1')
 	#ax.plot(widths, max_order - torus_min_energies, color='C1', linestyle='dotted')
-	ax.plot(widths, min_unorder - torus_min_energies, color='C0')
+	ax.plot(widths, min_unorder_off, color='C0')
-	ax.plot(widths, max_unorder - torus_min_energies, color='C0', linestyle='dotted')
+	ax.plot(widths, max_unorder_off, color='C0', linestyle='dotted')
 	axis_settings(ax, widths)
 
 	for i in null_unorder:
 		ax.scatter(widths[i], min_unorder[i] - torus_min_energies[i],
 			marker='X', color="blue", s=50, zorder=4)
-		ax.scatter(widths[i], max_unorder[i] - torus_min_energies[i],
+		# ax.scatter(widths[i], max_unorder[i] - torus_min_energies[i],
-			marker='X', edgecolors="blue", facecolors='none', s=100, zorder=4)
+		# 	marker='X', edgecolors="blue", facecolors='none', s=100, zorder=4)
 
 	# for i, marker in enumerate(min_markers):
 	# 	if marker:
@@ -207,9 +213,13 @@ def main():
 	ax.title.set_text('Reduced Energy vs. Width')
 	ax.set_xlabel("Width")
 	ax.set_ylabel("Reduced Energy")
-
+	bif_y_max = np.max(np.abs(np.concatenate((min_unorder_off, max_unorder_off))))
+	ax.set_yticks(np.arange(-bif_y_max, bif_y_max, round(bif_y_max/20, \
+											-math.floor(math.log10(bif_y_max/20)))))
 	fig.savefig(fig_folder / "Bifurcation.png")
 
+	print(f"Wrote to {fig_folder}.")
+
 if __name__ == "__main__":
 	os.environ["QT_LOGGING_RULES"] = "*=false"
 	SIM_FOLDER = Path(f"simulations/ShrinkEnergyComparison")

simulation.py

@@ -461,20 +461,27 @@ class Simulation:
 
 	def get_distinct(self):
 		distinct_hists = []
+		distinct_avg_radii = []
 		new_frames = []
 
 		for frame in self.frames:
-			new_hist = np.histogram(frame.stats["avg_radius"], bins=10)[0]
+			if np.var(frame.stats["avg_radius"]) <= 1e-8:
+				avg_radii = np.average(frame.stats["avg_radius"])
+				if not np.any(np.isclose(avg_radii, distinct_avg_radii, atol=1e-5)):
+					distinct_avg_radii.append(avg_radii)
+					new_frames.append(frame)
+			else:
+				new_hist = np.histogram(frame.stats["avg_radius"], bins=8)[0]
 
-			is_in = False
+				is_in = False
-			for hist in distinct_hists:
+				for hist in distinct_hists:
-				if np.allclose(new_hist, hist, atol=1e-8):
+					if np.all(hist == new_hist):
-					is_in = True
+						is_in = True
-					break
+						break
 
-			if not is_in:
+				if not is_in:
-				distinct_hists.append(new_hist)
+					distinct_hists.append(new_hist)
-				new_frames.append(frame)
+					new_frames.append(frame)
 
 		self.frames = new_frames