Packaged into working program with arguments

.gitignore

@@ -7,3 +7,4 @@ src/packsim.c
 
 figures
 simulations
+new_simulations

README.md

@@ -1 +1,145 @@
-# packsim
+# PackSim
+
+PackSim is Python program which perform simulations for the flow of 'soft' or 'compressible' objects under some energy in a periodic domain.
+
+
+## Installing
+Currently, this isn't packaged as a standalone program. Python 3.8+ is needed to use this.
+
+## Building
+
+Before running, first install the necessary requirements. A virtual environment is recommended.
+
+```bash
+pip install -r requirements.txt
+```
+It is also necessary to build the Cython components. So, simply run:
+
+```bash
+foo@bar:/path/to/packsim: ./build.sh
+```
+If video output of the simulation is desired, `ffmpeg` **must** be installed on system. Otherwise, it's all ready to go!
+
+## Usage
+
+The primary usage is to run `packsim.py` to output simulation results. There are many parameters, so it's recommended that you use a JSON configuration file. There are two required sections, `calculation` and `simulation`.
+
+### Domain
+
+The 'domain' section defines the basic setup. The `n_objects` parameter represents the number of objects in the domain of `width` by  `height`. The `natural_radius` parameter is the radius of the circle that represents the natural 'zero energy' state of the object. Lastly, the `energy` parameter is just the type of energy to simulate with. There are currently three: `'area', 'radial-al', 'radial-t'`. These parameters can also be overridden by a command line argument.
+
+
+
+```json
+{
+  "domain": {
+    "n_objects": 15,
+    "width": 10.0,
+    "height": 10.0,
+    "natural_radius": 4.0,
+    "energy": "radial-t",
+    "points": [  // Optional
+      [1,1], [2,2], [3,3], [4,4], [5,5],
+      [1,2], [2,3], [3,4], [4,5], [5,6],
+      [1,3], [2,4], [3,5], [3,6], [5,7]
+    ]
+  },
+  ...
+}
+```
+
+### Simulation
+
+There are currently three simulation modes, and the configuration changes slightly for each one.
+
+
+#### Flow
+This mode simulates the relaxing of the objects to its equilibrium. The threshold is the sufficient condition where the gradient is sufficiently close to zero. Specifically, the simulation stops when the L1 norm of the gradient divided by the number of objects is less than the threshold. The `step-size` parameter only represents the *initial* step size. This is because adaptive step size is employed.
+
+```json
+{
+  ...
+  "simulation": {
+    "mode": "flow",
+    "step_size": 0.05,
+    "threshold": 0.0001
+  },
+  ...
+}
+```
+
+#### Search
+This mode searches for equilibrium until `eq_stop_count` equilibria are found. Additionally, if the nullity of the Hessian at the equilibrium is greater than 2. (This is due to periodicity, as any translation is another equilibrium.) In this case, the `manifold_step_size` parameter is used to traverse along it.
+
+```json
+{
+  ...
+  "simulation": {
+    "mode": "search",
+    "step_size": 0.05,
+    "threshold": 0.0001,
+    "eq_stop_count": 100,
+    "manifold_step_size": 0.1
+  },
+  ...
+}
+```
+
+#### Shrink
+This mode simulates the the change in the equilibrium as the width is decreased. Both the `width_change` and `width_stop` parameters should be set as a percentage of the width.
+```json
+{
+  ...
+  "simulation": {
+    "mode": "shrink",
+    "step_size": 0.05,
+    "threshold": 0.0001,
+    "width_change": 1,
+    "width_stop": 0.3
+  },
+  ...
+}
+```
+
+Additionally, a filename is automatically generated by default, but it's also possible to provide it with the `name` parameter. The `save_sim` parameter determines whether or not the `.sim` file is saved. This allows simulations to be loaded again at a later time, if other processing is desired. Setting points is also possible, with the `points` parameter.
+
+```json
+{
+  ...
+  "simulation": {
+    ...
+    "name": "my_awesome_simulation", // Optional
+    "save_sim": true
+    ...
+  }
+  ...
+}
+```
+
+### Diagram
+It's also possible to visualize the data, but this section is optional. It's possible to output the frames as individual images [`img`], or as a video [`mp4`]. (Note: rendering as an `mp4` requires `ffmpeg` to be installed on your system.)
+
+While it's possible to customize the figures that are drawn and outputted, there are already a few preset modes: `animate`, `energy`, `stats`, and `eigs`, which outputs the simulation steps as a video, a video with the graph of energy, compiled statistics, and the eigenvalues, respectively.
+
+```json
+{
+  ...
+  "diagram": {
+    "filetype": "mp4",
+    "figures": "animate",
+  }
+}
+```
+
+Now with, all that configuration nonsense out of the way, to run, simply type
+```bash
+foo@bar:/path/to/packsim: python3 packsim my_sim.json
+```
+where `my_sim.json` is your configuration file.
+
+## Contributions
+
+This project welcomes contributors, so feel free to make a pull request!
+
+## License
+[GNU GPLv3.0](https://choosealicense.com/licenses/gpl-3.0/)

build.sh

@@ -1,3 +1,4 @@
+rm -f packsim_core*
 cd src
 python3 setup.py build_ext --inplace --quiet
 mv *.so ../

check_width_exists.py

@@ -1,3 +1,5 @@
+#!/usr/bin/env python3
+
 from pathlib import Path
 import sys, numpy as np
 

convert_old.py

@@ -0,0 +1,22 @@
+from __future__ import annotations
+from typing import List
+from simulation import Diagram, Simulation
+import argparse, numpy as np
+import matplotlib.pyplot as plt
+from pathlib import Path
+import os
+
+def main():
+	nums = [13, 23, 57, 59, 83, 131, 179]
+
+	for n in nums:
+		Path(f"new_simulations/Radial[T]T - N{n}R4.0").mkdir(exist_ok=True)
+		for file in Path(f"simulations/Radial[T]T - N{n}R4.0").iterdir():
+			sim = Simulation.load(file)
+			sim.get_distinct()
+			sim.save(file.name[:-4].replace("x10.0", "x10.00"))
+	#sim = Simulation.load("simulations/AreaT - N30R4 - 10x10.sim")
+
+
+if __name__ == '__main__':
+	main()

packsim.py

@@ -1,68 +1,141 @@
-#!/usr/bin/env python3
-
 from __future__ import annotations
-import argparse, json
+from typing import Dict
+import argparse, json, numpy as np, os
+from shutil import which
 from simulation import Diagram, Flow, Search, Shrink
 
+from packsim_core import RadialTEnergy
 
-def get_diagram(sim, t):
+dia_presets = {
-	if t == "flow":
+	"animate": [["voronoi"]],
-		diagram = Diagram(sim, np.array([["voronoi", "energy"]]))
+	"energy": [["voronoi", "energy"]],
-	elif t == "stats":
+	"stats": [
-		diagram = Diagram(sim, np.array([
 		["voronoi", "eigs", "site_edge_count"],
 		["site_isos", "site_energies", "edge_lengths"]
-				]), cumulative=False)
+	],
-	elif t == "eigs":
+	"eigs": [["voronoi", "eigs"]]
-		diagram = Diagram(sim, np.array([["voronoi", "eigs"]]))
+}
-	elif t == "shrink":
-		diagram = Diagram(sim, np.array([["voronoi", "avg_radius", "isoparam_avg"]]), cumulative=False)
 
-	return diagram
+def check_params(container: Dict, needed: List[str], valid: Dict):
+	"""
+	Checks container for the necessary items, and raises
+	an error if the parameter is not found.
+	:param container: [Dict] contains the submitted parameters.
+	:param needed: [List[str]] contains the needed parameters.
+	:param valid: [Dict] if there are specific valid parameters,
+	will also check for those.
+	"""
+	for need in needed:
+		if need not in container:
+			raise ValueError(f"Parameter \'{need}\' is required.")
+
+		if need in valid:
+			if type(valid[need]) is list:
+				if container[need] not in valid[need]:
+					raise ValueError(f"Parameter \'{need}\' must be one of these values: " + \
+										 f"{str(valid[need])[1:-1]}.")
+			elif valid[need] == "positive":
+				if container[need] < 0:
+					raise ValueError(f"Parameter \'{need}\' must be positive.")
 
 
 def main():
 	# Loading configuration and settings.
-	parser = argparse.ArgumentParser("Processes packing simulations.")
+	parser = argparse.ArgumentParser("PackSim")
-	parser.add_argument('sim_conf', metavar='path/to/config',
+	parser.add_argument('sim_conf', metavar='/path/to/config.json',
 						help="configuration file for a simulation")
 	parser.add_argument('-q', '--quiet', dest='quiet', action='store_true', default=False,
 						help="suppress all normal output")
-	parser.add_argument('-l', '--log', dest='log_steps', action='store_true', default=50,
+	parser.add_argument('-l', '--log', dest='log_steps', default=50, type=int,
 						help="number of iterations before logging")
-	parser.add_argument('-i', '--input', dest='input_file')
+
-	parser.add_argument('-o', '--output', dest='output_file')
+	parser.add_argument('--n_objects', dest='n', type=int, help="objects in domain")
+	parser.add_argument('--width', dest='w', type=float, help="width of domain")
+	parser.add_argument('--height', dest='h', type=float, help="height of domain")
+	parser.add_argument('--natural_radius', dest='r', type=float, help="natural radius of object")
+	parser.add_argument('--energy', dest='energy', help="energy type of system")
 
 	args = parser.parse_args()
-
-	if args.input_file is None:
 	config_sim(args)
-	else:
+	# if args.input_file is None:
-		loaded_sim(args)
+	# 	config_sim(args)
+	# else:
+	# 	loaded_sim(args)
 
 
 def config_sim(args):
 	with open(args.sim_conf) as f:
 		params = json.load(f)
 
-	calc_params, sim_params = params["calc"], params["sim"]
+	check_params(params, ["domain", "simulation"], {})
-	n, w, h, r, energy = calc_params["n_objects"], calc_params["width"], calc_params["height"], \
+	dmn_params, sim_params = params["domain"], params["simulation"]
-				calc_params["natural_radius"], calc_params["energy"]
 
-	mode, thres, step = sim_params["mode"], sim_params["threshold"], sim_params["step_size"]
+	overrides = {args.n: "n_objects", args.w: "width", args.h: "height", args.r: "natural_radius",
+				 args.energy: "energy"}
+	for arg, arg_name in overrides.items():
+		if arg is not None:
+			dmn_params[arg_name] = arg
+
+	check_params(dmn_params, ["n_objects", "width", "height", "natural_radius", "energy"], {
+		"n_objects": "positive", "width": "positive", "height": "positive",
+		"natural_radius": "positive", "energy": ["area", "radial-al", "radial-t"]
+	})
+	n, w, h, r, energy = dmn_params["n_objects"], dmn_params["width"], dmn_params["height"], \
+				dmn_params["natural_radius"], dmn_params["energy"]
+
+
+
+	points = None
+	if "points" in dmn_params:
+		points = np.asarray(dmn_params["points"])
+
+	check_params(sim_params, ["mode", "step_size", "threshold", "save_sim"], {
+		"mode": ["flow", "search", "shrink"], "step_size": "positive", "threshold": "positive"
+	})
+	mode, step, thres, save_sim = sim_params["mode"], sim_params["step_size"], \
+									sim_params["threshold"], sim_params["save_sim"]
+
+	name = sim_params.get("name")
 
-	# Running simulation
 	if mode == "flow":
 		sim = Flow(n, w, h, r, energy, thres, step)
 	elif mode == "search":
+		check_params(sim_params, ["manifold_step_size"], {"manifold_step_size": "positive"})
 		sim = Search(n, w, h, r, energy, thres, step, sim_params["manifold_step"], 
 						sim_params["count"])
 	elif mode == "shrink":
-		sim = Shrink(n, w, h, r, energy, thres, step, sim_params["delta_width"], 
+		check_params(sim_params, ["width_change", "width_stop"], {
-						sim_params["stop_width"])
+			"width_change": "positive", "width_stop": "positive"
+		})
+		sim = Shrink(n, w, h, r, energy, thres, step, sim_params["width_change"],
+						sim_params["width_stop"])
 
-	sim.initialize()
+	save_diagram = False
+	if "diagram" in params:
+		save_diagram = True
+		dia_params = params["diagram"]
+		check_params(dia_params, ["filetype", "figures"], {
+			"filetype": ["img", "mp4"]
+		})
+		if dia_params["filetype"] == "mp4":
+			if which("ffmpeg") is None:
+				raise ValueError("The program 'ffmpeg' needs to be installed on your system.")
+
+		if type(dia_params["figures"]) is str:
+			dia_params["figures"] = np.asarray(dia_presets[dia_params["figures"]])
+		else:
+			dia_params["figures"] = np.asarray(dia_params["figures"])
+
+	sim.initialize(points)
 	sim.run(not args.quiet, args.log_steps)
+	if save_sim: sim.save(filename=name)
+
+	if save_diagram:
+		diagram = Diagram(sim, dia_params["figures"])
+		if dia_params["filetype"] == "img":
+			diagram.render_static(0, filename=name)
+		elif dia_params["filetype"] == "mp4":
+			diagram.render_video(filename=name)
 
 
 def loaded_sim(args):
@@ -70,6 +143,7 @@ def loaded_sim(args):
 
 
 if __name__ == '__main__':
+	os.environ["QT_LOGGING_RULES"] = "*=false"
 	try:
 		main()
 	except KeyboardInterrupt:

shrink_energy_comparison.py

@@ -1,9 +1,7 @@
-#!/usr/bin/env python3
-
 from __future__ import annotations
 from typing import List
 from simulation import Diagram, Simulation
-import argparse, numpy as np
+import os, argparse, numpy as np
 import matplotlib.pyplot as plt
 from pathlib import Path
 
@@ -13,10 +11,19 @@ def get_torus_config_energies(n: int, widths: np.ndarray, h: float, r: float,
 	torus_min_energies, torus_max_energies = np.empty(widths.shape), np.empty(widths.shape)
 	for i, w in enumerate(widths):
 		sim = Simulation(n, w, h, r, energy)
+		configs = []
 
+		for j in range(1):
 			for c in range(1,n):	# Ignore 0, tends to error.
-			sim.add_frame(torus=(1,c))
+				config = (1,c) if j == 0 else (c,1)
-			sim.add_frame(torus=(c,1))
+				sim.add_frame(torus=config)
+				configs.append(config)
+
+				eigs = np.sort(np.linalg.eig(sim.frames[-1].hessian(10e-5))[0])
+				zero_eigs = np.count_nonzero(np.isclose(eigs, np.zeros((len(eigs),)), atol=1e-4))
+				if zero_eigs != 2:
+					del sim.frames[-1]
+					del config[-1]
 
 				hashes = int(21*i/len(widths))
 				print(f'Generating at width {w:.02f}... ' + \
@@ -30,47 +37,6 @@ def get_torus_config_energies(n: int, widths: np.ndarray, h: float, r: float,
 	return torus_min_energies, torus_max_energies
 
 
-# def equal_shape_eigs(n, widths, h, r):
-# 	n,w,h,r = 57, 10, 10, 4	# Domain settings
-# 	thres, step_size = 10e-5, 5e-2	# Simulation settings
-# 	log_steps = 50
-# 	energy = "radial-t"
-
-# 	sims = [None]*n*2
-# 	energies = {}
-
-# 	for x in range(1,n):
-# 		sim = TravelEQ(n, w, h, r, energy, thres, step_size, log_steps)
-# 		sim2 = TravelEQ(n, w, h, r, energy, thres, step_size, log_steps)
-# 		#frame = FindEQ(n, w, h, r, "radial-t", POOL, thres, step_size, log_steps)
-# 		for j in range(141):
-# 			sim.w = 10-j*.05
-# 			sim2.w = 10-j*.05
-# 			sim.add_frame(None, (1,x), 0)
-# 			sim2.add_frame(None, (x, 1), 0)
-# 		#sim.initialize(torus=(1,x))
-
-# 		energies[(1,x)] = sim[0].energy
-# 		energies[(x,1)] = sim2[0].energy
-
-# 		sims[x] = list([y.energy for y in sim.frames])
-# 		sims[x+n] = list([y.energy for y in sim2.frames])
-# 		#k1 = np.concatenate(sim.frames[0].process(sim.frames[0].grad, sim.frames[0].get_ranges()))
-# 		#print(np.linalg.norm(k1))
-
-# 		# hess = sim.frames[0].hessian(10e-5)
-# 		# eigs = np.sort(np.linalg.eig(hess)[0])
-# 		# sim.frames[0].stats["eigs"] = eigs
-		
-# 		# diagram = Diagram(sim, np.array([["voronoi", "eigs"]]))
-# 		#diagram = Diagram(sim, np.array([["voronoi"]]))
-# 		#diagram.render_static(0, filename=f'EqualShape/EqualShapeN{n}/{str((1, x))}')
-
-# 	print(min(energies, key=energies.get))
-	
-# 	return sims
-
-
 def main():
 	# Loading arguments.
 	parser = argparse.ArgumentParser("Compiles the equilibriums for each width into a diagram.")
@@ -94,6 +60,11 @@ def main():
 	print(flush=True)
 	sims.sort(key=lambda x: x.w)
 
+	sim_folder = Path(f"simulations/ShrinkEnergyComparison")
+	fig_folder = Path(f"figures/ShrinkEnergyComparison - N{sims[0].n}")
+	sim_folder.mkdir(exist_ok=True)
+	fig_folder.mkdir(exist_ok=True)
+
 	widths = np.asarray([sim.w for sim in sims])
 
 	min_frames = [min(sim.frames, key=lambda x: x.energy) for sim in sims]
@@ -106,22 +77,17 @@ def main():
 		sims[0].n, widths, sims[0].h, sims[0].r, sims[0].energy
 	)
 
-	min_markers = [np.var(frame.stats["site_areas"]) <= 1e-8 for frame in min_frames]
+	min_markers = [np.var(frame.stats["avg_radius"]) <= 1e-8 for frame in min_frames]
-	max_markers = [np.var(frame.stats["site_areas"]) <= 1e-8 for frame in max_frames]
+	max_markers = [np.var(frame.stats["avg_radius"]) <= 1e-8 for frame in max_frames]
 
 	# Torus minimum energies used as reference.
 
 	fig, ax = plt.subplots(figsize=(16, 8))
 	#ax.plot(widths, nums)
-	# for i, equal_sim in enumerate(equal_sims):
+	ax.plot(widths, [len(sim.frames) for sim in sims])
-	# 	if i in [0, n]:
+	fig.savefig(folder / "")
-	# 		continue
-
-	# 	ax.plot(widths, 
-	# 		np.asarray(equal_sims[i]) - reference,
-	# 		color="orange", alpha=0.5, linewidth=0.5, zorder=3
-	# 	)
 
+	fig, ax = plt.subplots(figsize=(16, 8))
 	ax.plot(widths, torus_min_energies - torus_min_energies, color='C1')
 	ax.plot(widths, min_energies - torus_min_energies, color='C0')
 	ax.plot(widths, max_energies - torus_min_energies, color='C0', linestyle='dotted')
@@ -150,15 +116,16 @@ def main():
 	ax.set_ylabel("Reduced Energy")
 	ax.grid(zorder=0)
 
-	#ax.set_xticks([round(w,2) for w in widths[::-2]])
+	ax.set_xticks([round(w,2) for w in widths[::-2]])
 	#ax.set_yticks(np.arange(-920, 1120, 40))
-	#ax.set_xticklabels(ax.get_xticks(), rotation = 90)
+	ax.set_xticklabels(ax.get_xticks(), rotation = 90)
 
 	plt.tight_layout()
 
 	fig.savefig(f"figures/WidthsEnergyComparison - N{sims[0].n}.png")
 
 if __name__ == "__main__":
+	os.environ["QT_LOGGING_RULES"] = "*=false"
 	try:
 		main()
 	except KeyboardInterrupt:

simulation.py

@@ -5,7 +5,7 @@ import matplotlib.pyplot as plt
 from matplotlib.ticker import MaxNLocator, FormatStrFormatter
 import os, math, random, time, pickle, scipy, numpy as np
 
-from packsim import VoronoiContainer, AreaEnergy, RadialALEnergy, RadialTEnergy
+from packsim_core import VoronoiContainer, AreaEnergy, RadialALEnergy, RadialTEnergy
 from timeit import default_timer as timer
 
 
@@ -25,7 +25,7 @@ def gen_filepath(sim: Simulation, ext: str, parent_dir='figures') -> str:
 				RadialTEnergy: "Radial[T]"}[sim.energy]
 	mode = {Flow: "F", Search: "T", Shrink: "S"}[type(sim)]
 
-	base_filename = f'{energy}{mode} - N{sim[0].n}R{sim[0].r} - {round(sim[0].w, 2):.2f}x{sim[0].h}'
+	base_filename = f'{energy}{mode} - N{sim[0].n}R{round(sim[0].r, 1) :.1f} - {round(sim[0].w, 2):.2f}x{round(sim[0].h, 2):.2f}'
 	base_path = f'{parent_dir}/{base_filename}'
 	
 	i = 1
@@ -275,12 +275,16 @@ class Diagram():
 			os.mkdir(path)
 			for frame in range(i, j+1):
 				self.generate_frame(frame)
-				if frame % 20 == 0:
+
-					print(f'Rendered frame {frame}/{length} : {100*frame/length:.2f}%')
+				hashes = int(21*i/(j+1))
+				print(f'Generating frames... |{"#"*hashes}{" "*(20-hashes)}|' + \
+					f' {i+1}/{j+1} frames rendered.', flush=True, end='\r')
+
 				plt.savefig(f'{path}/img{frame:03}.png')
 				plt.close()
 
-			print(f'Wrote to folder \"{path}\"')
+			print(flush=True)
+			print(f'Wrote to folder \"{path}\"', flush=True)
 			
 
 	def render_video(self, time = 30, fps = None, filename = None):
@@ -305,15 +309,18 @@ class Diagram():
 		except FileExistsError: 
 			pass
 
-		print("Generating frames...")
 		frames = min(len(self.sim), int(fps * time))
 		for j in range(frames):
 			self.generate_frame(int(j*step))
-			if j % 20 == 0:
+			hashes = int(21*j/frames)
-				print(f'Rendered frame {j}/{frames} : {100*j/frames:.2f}%')
+			print(f'Generating frames... |{"#"*hashes}{" "*(20-hashes)}|' + \
+					f' {j+1}/{frames} frames rendered.', flush=True, end='\r')
+
 			plt.savefig(f'figures/temp/img{j:03}.png')
 			plt.close()
 
+		print(flush=True)
+
 
 		if filename is None:
 			path = gen_filepath(self.sim, "mp4")
@@ -321,7 +328,7 @@ class Diagram():
 			path = f'figures/{filename}.mp4'
 
 		# Convert to gif.
-		print("Assembling MP4...")
+		print("Assembling MP4...", flush=True)
 		os.system(f'ffmpeg -hide_banner -loglevel error -r {fps} -i figures/temp/img%03d.png' + \
 				  f' -c:v libx264 -crf 18 -preset slow -pix_fmt yuv420p -vf' + \
 				  f' "scale=trunc(iw/2)*2:trunc(ih/2)*2" -f mp4 "{path}"')
@@ -331,7 +338,7 @@ class Diagram():
 			os.remove(f'figures/temp/img{j:03}.png')
 
 		os.rmdir("figures/temp")
-		print(f'Wrote to \"{path}\".')
+		print(f'Wrote to \"{path}\".', flush=True)
 
 
 class Simulation:
@@ -428,7 +435,7 @@ class Simulation:
 		else:
 			values = self.frames[i].stats[stat]
 
-		bins = 9
+		#bins = 9
 		if np.var(values) <= 1e-8:
 			hist = np.zeros((bins,))
 			val = np.average(values)
@@ -452,26 +459,24 @@ class Simulation:
 		# return (labels, count, colors)
 
 
-	def get_distinct(self) -> Simulation:
+	def get_distinct(self):
-		distinct_eigs = []
+		distinct_hists = []
 		new_frames = []
 
 		for frame in self.frames:
-			new_eigs = frame.stats["eigs"]
+			new_hist = np.histogram(frame.stats["avg_radius"], bins=10)[0]
 
 			is_in = False
-			for eigs in distinct_eigs:
+			for hist in distinct_hists:
-				if np.allclose(new_eigs, eigs, atol=1e-4):
+				if np.allclose(new_hist, hist, atol=1e-8):
 					is_in = True
+					break
 
 			if not is_in:
-				distinct_eigs.append(new_eigs)
+				distinct_hists.append(new_hist)
 				new_frames.append(frame)
-				continue
 
-		new_sim = self.__class__(self.n, self.w, self.h, self.r, self.energy)
+		self.frames = new_frames
-		new_sim.frames = new_frames
-		return new_sim
 
 
 	def save(self, filename: str = None):
@@ -482,25 +487,22 @@ class Simulation:
 		if filename is None:
 			path = gen_filepath(self, "sim", "simulations")
 		else:
-			path = f'simulations/{filename}.sim'
+			path = f'new_simulations/{filename}.sim'
-
-		# Convert sites to NumPy array.
-		arr = np.zeros((len(self.frames), self.n, 2))
-		arr = np.stack([frame.site_arr for frame in self.frames])
-		#stats = [frame.stats for frame in self.frames]
 
 		all_info = []
 		for frame in self.frames:
 			frame_info = dict()
 			frame_info["arr"] = frame.site_arr
-			frame_info["energy"] = {AreaEnergy: "Area", RadialALEnergy: "Radial[AL]", 
+			frame_info["energy"] = {AreaEnergy: "area", RadialALEnergy: "radial-al",
-									RadialTEnergy: "Radial[T]"}[sim.energy]
+									RadialTEnergy: "radial-t"}[self.energy]
 			frame_info["params"] = (frame.n, frame.w, frame.h, frame.r)
 			all_info.append(frame_info)
 
+		class_name = {Flow: "flow", Search: "search", Shrink: "shrink"}[self.__class__]
+
 		with open(path, 'wb') as output:
-			pickle.dump((all_info, self.__class__), output, pickle.HIGHEST_PROTOCOL)
+			pickle.dump((all_info, class_name), output, pickle.HIGHEST_PROTOCOL)
-		print("Wrote to " + path)
+		print("Wrote to " + path, flush=True)
 
 
 	@staticmethod
@@ -512,10 +514,11 @@ class Simulation:
 		frames = []
 		with open(filename, 'rb') as data:
 			all_info, sim_class = pickle.load(data)
+
 			sim = sim_class(*all_info[0]["params"], all_info[0]["energy"], 0,0,0,0)
 			for frame_info in all_info:
-				frames.append(frame_info["energy"](*frame_info["params"], frame_info["arr"]))
+				frames.append(sim.energy(*frame_info["params"], frame_info["arr"]))
-				frames[-1].stats = frame_info["stats"]
+				#frames[-1].stats = frame_info["stats"]
 
 			sim.frames = frames
 		return sim	
@@ -563,50 +566,46 @@ class Flow(Simulation):
 		"""
 		if log:
 			print(f'Find - N = {self.n}, R = {self.r}, {self.w} X {self.h}', flush=True)
-		i, grad_mag = 0, float('inf')
+		i, grad_norm = 0, float('inf')
-
-		## Replace with adaptive step size eventually!!!!!!
 
 		trial = 2
-		while grad_mag > self.thres:	# Get to threshold.
+		while grad_norm > self.thres:	# Get to threshold.
-			# Iterate and generate next frame using Euler method. 
+			# Iterate and generate next frame using RK-3
 			start = timer()
-			new_sites, DE = self.frames[i].iterate(self.step_size)
+			change, grad = self.frames[i].iterate(self.step_size)
-			orig_step = self.energy(self.n, self.w, self.h, self.r, new_sites)
+			new_frame = self.energy(self.n, self.w, self.h, self.r,
-			grad_mag = np.linalg.norm(DE)
+							self.frames[i].add_sites(change))
+			grad_norm = np.sum(np.absolute(grad))/self.n
 			end = timer()
 
-
+			if new_frame.energy < self.frames[i].energy:	# If energy decreases.
-			if orig_step.energy < self.frames[i].energy:	# If energy decreases.
 				if trial < 20:	# Try increasing step size for 10 times.
 					factor = 1 + .1**trial
 
-					test_step = self.energy(self.n, self.w, self.h, self.r, 
+					test_frame = self.energy(self.n, self.w, self.h, self.r,
-											self.frames[i].add_sites(self.step_size*factor*DE))
+											self.frames[i].add_sites(change*factor))
 					# If increased step has less energy than original step.
-					if test_step.energy < orig_step.energy:
+					if test_frame.energy < new_frame.energy:
 						self.step_size *= factor
 						trial = max(2, trial-1)
-						grad_mag = np.linalg.norm(DE)
+						new_frame = test_frame
-						new_step = test_step
 					else:	# Otherwise, increases trials, and use original.
 						trial += 1
-						new_step = orig_step
-				else:
-					new_step = orig_step
 			else:	# Step size too large, decrease and reset trial counter.
-				self.step_size /= (1 + .1**(trial-1))
+				trial = 2
-				trial = 1
+				shrink_factor = 1.5
-				new_sites, DE = self.frames[i].iterate(self.step_size)
+				new_frame = self.energy(self.n, self.w, self.h, self.r,
-				new_step = self.energy(self.n, self.w, self.h, self.r, new_sites)
+								self.frames[i].add_sites(change/shrink_factor))
-			self.frames.append(new_step)
+				self.step_size /= shrink_factor
+
+			#self.step_size *= abs(.01/np.linalg.norm(error))**(1/3)
+			#self.step_size = max(10e-4, self.step_size)
+			self.frames.append(new_frame)
 
-			self.step_size = max(10e-4, self.step_size)
 			i += 1
-
 			if(log and i % log_steps == 0): 
 				print(f'Iteration: {i:05} | Energy: {self.frames[i].energy: .5f}' + \
-				 	  f' | Gradient: {grad_mag:.8f} | Step: {self.step_size: .5f} | ' + \
+				 	  f' | Gradient: {grad_norm:.8f} | Step: {self.step_size: .5f} | ' + \
 				 	  f'Time: {end-start: .3f}', flush=True)
 
 
@@ -656,7 +655,7 @@ class Search(Simulation):
 
 			self.frames[i] = sim[-1] # Replace frame with equilibrium frame.
 			if log:
-				print(f'Equilibrium: {i:04}\n')
+				print(f'Equilibrium: {i:04}\n', flush=True)
 			# Calculate kernel, and travel in some direction.
 
 			hess = self.frames[i].hessian(10e-5)
@@ -731,4 +730,3 @@ class Shrink(Simulation):
 		del self.frames[0]
 
 TravelEQ = Search
-

src/core.pyx

@@ -6,7 +6,7 @@ from cpython cimport array
 from libc.stdlib cimport malloc, realloc, calloc, free
 from libc.math cimport isnan, NAN, pi as PI, M_PI_2 as PI_2, \
 				sqrt, log, sin, cos, tan, acos, fabs
-from packsim cimport INT_T, FLOAT_T, Init, IArray, FArray, BitSet, Vector2D, Matrix2x2, \
+from packsim_core cimport INT_T, FLOAT_T, Init, IArray, FArray, BitSet, Vector2D, Matrix2x2, \
 					VectorSelfOps, VectorCopyOps, MatrixSelfOps, MatrixCopyOps, \
 					SiteCacheMap, EdgeCacheMap, VoronoiInfo, Site, HalfEdge
 

src/setup.py

@@ -2,7 +2,7 @@ from setuptools import Extension, setup
 from Cython.Build import cythonize
 import numpy
 
-MODULE_NAME = "packsim"
+MODULE_NAME = "packsim_core"
 
 ext_modules = [
     Extension(

src/voronoi_dcel.pyx

@@ -1,4 +1,4 @@
-from packsim cimport SiteCacheMap, EdgeCacheMap, VoronoiInfo, Site, HalfEdge
+from packsim_core cimport SiteCacheMap, EdgeCacheMap, VoronoiInfo, Site, HalfEdge
 
 #### Constants ####
 
@@ -687,11 +687,16 @@ cdef class VoronoiContainer:
 				self.add_sites(step*k1/2)
 		).gradient
 
+		lower = step*(-k1+ 2*k2)
 		k3 = self.__class__(self.n, self.w, self.h, self.r, 
-				self.add_sites(step*(-k1+ 2*k2))
+				self.add_sites(lower)
 		).gradient
 
-		return self.add_sites((step/6)*(k1+2*k2+k3)), k1
+		higher = (step/6)*(k1+2*k2+k3)
+		#new_sites = self.add_sites(higher)
+		#error = higher - lower
+		
+		return higher, k1
 
 	def hessian(self, d: float) -> np.ndarray:
 		"""

test_sim.json

@@ -1,33 +1,19 @@
 {
-  "calc": {
+  "domain": {
-    "n_objects": 10,
+    "n_objects": 100,
     "width": 10.0,
     "height": 10.0,
     "natural_radius": 4.0,
     "energy": "radial-t"
   },
-  "sim": {
+  "simulation": {
     "mode": "flow",
     "step_size": 0.05,
-    "threshold": 0.0001
+    "threshold": 0.00001,
+    "save_sim": true
+  },
+  "diagram": {
+    "filetype": "mp4",
+    "figures": "energy"
   }
 }
-
-/*
-# ARR = np.array([
-#   [1, 1], [3, 1], [5, 1],
-#   [1, 3], [3, 3], [5, 3],
-#   [1, 5], [3, 5], [5, 5],
-#   [1, 7], [3, 7], [5, 7],
-#   [1, 9], [3, 9], [5, 9],
-#   [7, 1], [8, 1], [9, 1],
-#   [7, 2], [8, 2], [9, 2],
-#   [7, 3], [8, 3], [9, 3],
-#   [7, 4], [8, 4], [9, 4],
-#   [7, 5], [8, 5], [9, 5],
-#   [7, 6], [8, 6], [9, 6],
-#   [7, 7], [8, 7], [9, 7],
-#   [7, 8], [8, 8], [9, 8],
-#   [7, 9], [8, 9], [9, 9],
-# ], dtype=float)
-*/