diff --git a/scripts/check_width_exists.py b/scripts/check_width_exists.py
index c108839..88acdd7 100644
--- a/scripts/check_width_exists.py
+++ b/scripts/check_width_exists.py
@@ -1,16 +1,22 @@
 from pathlib import Path
-import sys, numpy as np
+import numpy as np, pickle, sys
+
+from squish import Simulation
 
 def main():
 	n = int(sys.argv[1])
 	all_widths = set(np.round(np.arange(3, 10.05, 0.05), 2))
-	for file in Path(f"simulations/Radial[T]T - N{n}R4.0").iterdir():
-		i = file.name.index("x")
-		all_widths.remove(float(file.name[i-4:i]))
+	for file in Path(f"squish_output/Radial[T]Search - N{n} - 500").iterdir():
+		sim, frames = Simulation.load(file / 'data.squish')
+
+		if sim.domain.n == n:
+			try:
+				all_widths.remove(next(frames)["domain"][1])
+			except StopIteration:
+				pass
 
 	remain_widths = sorted(list(all_widths))[::-1]
-	print(remain_widths)
-	print([int(round((10-w)/.05)) for w in remain_widths])
+	print("Remaining:", remain_widths)
 
 
 if __name__ == "__main__":
diff --git a/scripts/convergence.py b/scripts/convergence.py
index e33e542..08972fe 100644
--- a/scripts/convergence.py
+++ b/scripts/convergence.py
@@ -4,6 +4,7 @@ import argparse, pickle, numpy as np, os
 from pathlib import Path
 import matplotlib.pyplot as plt
 
+from squish import Simulation
 
 def main():
 	parser = argparse.ArgumentParser("Graphs convergence graphs for a collection of simulations.")
@@ -16,12 +17,14 @@ def main():
 	data = {}
 
 	for file in Path(args.sims_path).iterdir():
-		with open(file, "rb") as f:
-			all_info, _ = pickle.load(f)
+		sim = Simulation.load(file / "data.squish")
+		sim_info = next(sim)
+
+		step = sim_info["step_size"]
+		for frame in sim:
 
-		step = float(file.name[:-4].split("-")[1])
 		data[step] = {"times": [], "values": [], "diffs": []}
-		for i, frame_info in enumerate(all_info):
+		for i, frame_info in enumerate(sim):
 			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"]))
diff --git a/scripts/width_diagrams.py b/scripts/width_diagrams.py
index 3a1f884..5ca8298 100644
--- a/scripts/width_diagrams.py
+++ b/scripts/width_diagrams.py
@@ -1,101 +1,93 @@
 from __future__ import annotations
-from typing import List
-import os, math, argparse, numpy as np, pickle
+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
 
-from ..simulation import Diagram, Simulation
-from .._squish import AreaEnergy, RadialALEnergy, RadialTEnergy
+import squish.ordered as order
+from squish import Simulation, DomainParams
+from squish.common import OUTPUT_DIR
 
-ENERGY_R_STR = {AreaEnergy: "Area", RadialALEnergy: "Radial[AL]", RadialTEnergy: "Radial[T]"}
-ENERGY_I_STR = {AreaEnergy: "area", RadialALEnergy: "radial-al", RadialTEnergy: "radial-t"}
-I_TO_R = {"area": "Area","radial-t": "Radial[AL]", "radial-t": "Radial[T]"}
+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)))
 
-def get_torus_config_energies(n: int, widths: np.ndarray, h: float, r: float,
-								energy: str) -> Tuple:
-	sim_file = SIM_FOLDER / f"{I_TO_R[energy]} - TorusConfigEnergy - N{n}.data"
-	if sim_file.is_file():
-		with open(sim_file, "rb") as data:
-			return pickle.load(data)
-
-	torus_min_energies, torus_max_energies = np.empty(widths.shape), np.empty(widths.shape)
-	torus_min_configs, torus_max_configs = [None]*len(widths), [None]*len(widths)
-
-	for i, w in enumerate(widths):
-		sim = Simulation(n, w, h, r, energy)
-		configs = []
-
-		for j in range(2):
-			for c in range(1,n):	# Ignore 0, tends to error.
-				config = (1,c) if j == 0 else (c,1)
-				sim.add_frame(torus=config)
-				configs.append(config)
-
-				# eigs = np.sort(np.linalg.eig(sim.frames[-1].hessian(10e-5))[0])
-				# if eigs[0] > 1e-4:
-				# 	del sim.frames[-1]
-				# 	del config[-1]
-
-				hashes = int(21*i/len(widths))
-				print(f'Generating at width {w:.02f}... ' + \
-						f'|{"#"*hashes}{" "*(20-hashes)}| {i+1}/{len(widths)}, ' + \
-						f'{c + (n-1)*j}/{2*(n-1)} completed.', flush=True, end='\r')
-
-		pair = list(zip(configs,[frame.energy for frame in sim.frames]))
-		torus_min_configs[i], torus_min_energies[i] = min(pair, key=lambda x: x[1])
-		torus_max_configs[i], torus_max_energies[i] = max(pair, key=lambda x: x[1])
-
-	print(flush=True)
-
-	out_tup = (torus_min_energies, torus_max_energies, torus_min_configs, torus_max_configs)
-	with open(sim_file, "wb") as output:
-		pickle.dump(out_tup, output)
-
-	return out_tup
+	return domain.w, min(energies), max(energies)
 
 
-def get_equilibria_data(filepath: Path):
-	if filepath.is_file():
-		with open(filepath, "rb") as data:
-			return pickle.load(data)
+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 / '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())
 
-	for i, file in enumerate(files):
-		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,
-				np.count_nonzero(frame.stats["site_edge_count"] != 6)
-			])
+	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]
 
-		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,
-				np.count_nonzero(frame.stats["site_edge_count"] != 6)
-			])
-
-		hashes = int(21*i/len(files))
-		print(f'Loading simulations... |{"#"*hashes}{" "*(20-hashes)}|' + \
+			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)
+		print(flush=True)
 
+	sim, frames = Simulation.load(files[0] / 'data.squish')
 	widths = np.asarray(sorted(data["all"]))
-	n, h, r, energy = sim.n, sim.h, sim.r, sim.energy
+	domain = DomainParams(sim.domain.n, widths[-1], sim.domain.h, sim.domain.r)
+	return data, widths, domain
 
-	sim_file = SIM_FOLDER / f"{ENERGY_R_STR[energy]} - EquilibriaData - N{n}.data"
-	out_tup = (widths, data, n, h, r, ENERGY_I_STR[energy])
-	with open(sim_file, "wb") as output:
-		pickle.dump(out_tup, output)
-
-	return out_tup
 
 def axis_settings(ax, widths):
 	ax.invert_xaxis()
@@ -107,7 +99,7 @@ def axis_settings(ax, widths):
 
 def main():
 	# Loading arguments.
-	parser = argparse.ArgumentParser("Compiles the equilibriums for each width into a diagram.")
+	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,
@@ -115,13 +107,10 @@ def main():
 
 	args = parser.parse_args()
 
-	widths, data, n, h, r, energy = get_equilibria_data(Path(args.sims_path))
+	data, widths, domain = get_equilibria_data(Path(args.sims_path))
+	order_data = get_ordered_energies(domain, widths)
 
-	torus_min_energies, torus_max_energies, _, _ = get_torus_config_energies(
-		n, widths, h, r, energy
-	)
-
-	fig_folder = Path(f"figures/ShrinkEnergyComparison - N{n}")
+	fig_folder = OUTPUT_DIR / Path(f"ShrinkEnergyComparison - N{domain.n}")
 	fig_folder.mkdir(exist_ok=True)
 
 	# Torus minimum energies used as reference.
@@ -136,7 +125,7 @@ def main():
 	ax.plot(widths, all_disorder_count)
 	axis_settings(ax, widths)
 	ax.yaxis.set_major_formatter(mtick.PercentFormatter())
-	ax.title.set_text('Probability of Disorder')
+	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
@@ -152,15 +141,22 @@ def main():
 		distinct_ordered.append(equal_shape.count(True))
 		distinct_unordered.append(equal_shape.count(False))
 
-	ax.plot(widths, distinct_unordered, label="Unordered Equilibria")
-	ax.plot(widths, distinct_ordered, label="Ordered Equilibria")
-	ax.legend()
+	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('Density of States')
+	ax.title.set_text(f"Density of States - N{domain.n}")
 	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))
+	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
@@ -172,7 +168,7 @@ def main():
 
 	ax.plot(widths, defects)
 	axis_settings(ax, widths)
-	ax.title.set_text('Average Defects')
+	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))
@@ -187,53 +183,37 @@ def main():
 		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(torus_min_energies)):
-		ordered_energies[i].append(torus_min_energies[i])
-		ordered_energies[i].append(torus_max_energies[i])
-
+	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(torus_min_energies[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])
 
-	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')
+
+	min_unorder_off = min_unorder - min_order
+	max_unorder_off = max_unorder - min_order
+	ax.plot(widths, min_order - min_order, 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)
 
 	for i in null_unorder:
-		ax.scatter(widths[i], min_unorder[i] - torus_min_energies[i],
+		ax.scatter(widths[i], 0,
 			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] - offset[i],
 		# 	marker='X', edgecolors="blue", facecolors='none', s=100, zorder=4)
 
-	# for i, marker in enumerate(min_markers):
-	# 	if marker:
-	# 		ax.scatter(widths[i], min_energies[i]-torus_min_energies[i],
-	# 			marker='H', color="orange", s=20, zorder=4)
-	# 	else:
-	# 		ax.scatter(widths[i], min_energies[i]-torus_min_energies[i],
-	# 			marker='d', color="blue", s=20, zorder=4)
-
-	# for i, marker in enumerate(max_markers):
-	# 	if marker:
-	# 		ax.scatter(widths[i], max_energies[i]-torus_min_energies[i],
-	# 			marker='H', edgecolors="orange", s=20, facecolors='none', zorder=4)
-	# 	else:
-	# 		ax.scatter(widths[i], max_energies[i]-torus_min_energies[i],
-	# 			marker='d', edgecolors="blue", s=20, facecolors='none', zorder=4)
-
-	ax.title.set_text('Reduced Energy vs. Width')
+	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))))
@@ -245,8 +225,6 @@ def main():
 
 if __name__ == '__main__':
 	os.environ["QT_LOGGING_RULES"] = "*=false"
-	SIM_FOLDER = Path(f"simulations/ShrinkEnergyComparison")
-	SIM_FOLDER.mkdir(exist_ok=True)
 	try:
 		main()
 	except KeyboardInterrupt:
diff --git a/squish/common.py b/squish/common.py
index 5594deb..5aa5ff8 100644
--- a/squish/common.py
+++ b/squish/common.py
@@ -3,7 +3,7 @@ from typing import List, Tuple, Union, Optional, Iterator, Generator
 import pickle, numpy as np
 from math import gcd
 from pathlib import Path
-from ._squish import AreaEnergy, RadialALEnergy, RadialTEnergy
+from ._squish import VoronoiContainer, AreaEnergy, RadialALEnergy, RadialTEnergy
 
 OUTPUT_DIR = Path("squish_output")
 OUTPUT_DIR.mkdir(exist_ok=True)
@@ -84,7 +84,7 @@ class Energy:
 			except KeyError:
 				raise ValueError(f"\'{mode}\' is not a valid energy!")
 		else:
-			if mode is not VoronoiContainer and issubclaass(mode, VoronoiContainer):
+			if mode is not VoronoiContainer and issubclass(mode, VoronoiContainer):
 				raise ValueError("Provided class is not a valid energy!")
 			self.mode = mode
 
diff --git a/squish/ordered.py b/squish/ordered.py
index cdd9c0d..962cb91 100644
--- a/squish/ordered.py
+++ b/squish/ordered.py
@@ -36,8 +36,9 @@ def get_config_generators(domain: DomainParams, config: Config) -> Tuple[Config,
 	all_sites = np.concatenate((q1, q1-[w,0], q1-[w,h], q1-[0,h]))[2:]
 
 	# Checking 0 < ax + by < v*v to make the sites are within the region.
+	tol = 1e-3
 	vdot = np.matmul(all_sites, v)
-	in_box = all_sites[np.where((0 <= vdot) & (vdot <= v.dot(v)))[0]]
+	in_box = all_sites[np.where((-tol <= vdot) & (vdot <= (v.dot(v)+tol)))[0]]
 	in_box = np.expand_dims(in_box, 0).swapaxes(0,1)	# Used for the next step, getting site*site
 
 	w = in_box[np.argmin(np.squeeze(np.matmul(in_box, in_box.transpose(0,2,1))))].flatten()
diff --git a/squish/simulation.py b/squish/simulation.py
index 03dde75..8d2f8fe 100644
--- a/squish/simulation.py
+++ b/squish/simulation.py
@@ -62,8 +62,14 @@ class Simulation:
 
 		distinct_avg_radii, distinct_count, new_frames = [], [], []
 
+
 		for frame in self.frames:
-			avg_radii = np.sort(frame.stats["avg_radius"])
+			try:
+				stats = frame.stats
+			except AttributeError:
+				stats = frame["stats"]	# When we have a loaded simulations.
+
+			avg_radii = np.sort(stats["avg_radius"])
 			is_in = False
 			for i, dist_radii in enumerate(distinct_avg_radii):
 				if np.allclose(avg_radii, dist_radii, atol=1e-5):
@@ -73,6 +79,7 @@ class Simulation:
 
 			if not is_in:
 				distinct_avg_radii.append(avg_radii)
+				distinct_count.append(1)
 				new_frames.append(frame)
 
 		self.frames = new_frames
@@ -97,7 +104,7 @@ class Simulation:
 		f = self[index]
 		info = {
 			"arr": f.site_arr,
-			"domain": (f.n, f.h, f.w, f.r),
+			"domain": (f.n, f.w, f.h, f.r),
 			"energy": f.energy,
 			"stats": f.stats
 		}
@@ -123,8 +130,13 @@ class Simulation:
 	def load(path: str) -> Tuple[Simulation, Generator]:
 		def frames() -> Dict:
 			with open(path, 'rb') as infile:
+				first = True
 				while True:
 					try:
+						if first:
+							pickle.load(infile)
+							first = False
+							continue
 						yield pickle.load(infile)
 					except EOFError:
 						break
@@ -139,31 +151,6 @@ class Simulation:
 			return sim, frames()
 
 
-	@staticmethod
-	def load_old(filename: str) -> Simulation:
-		"""
-		Loads the points at every point into a file.
-		:param filename: [str] name of the file
-		"""
-		frames = []
-		with open(filename, 'rb') as data:
-			all_info, sim_class = pickle.load(data)
-			if type(sim_class) == str:
-				sim_class = {"flow": Flow, "search": Search, "shrink": Shrink}[sim_class]
-
-			if all_info[0]["params"][0] in [53, 59, 61, 83, 131]:
-				thres = 1e-5
-			else:
-				thres = 1e-6
-			sim = sim_class(DomainParams(*all_info[0]["params"]), Energy("radial-t"), 0.05, thres, True, 0.1, 500)
-			for frame_info in all_info:
-				frames.append(sim.energy.mode(*frame_info["params"], frame_info["arr"]))
-				#frames[-1].stats = frame_info["stats"]
-
-			sim.frames = frames
-		return sim
-
-
 class Flow(Simulation):
 	"""Finds an equilibrium from initial sites.
 
@@ -408,4 +395,6 @@ STR_TO_SIM = {
 	"flow": Flow,
 	"search": Search,
 	"shrink": Shrink
-}
\ No newline at end of file
+}
+
+simulation = Simulation
\ No newline at end of file