Added calculations for ordered sites

squish/common.py

@@ -1,6 +1,7 @@
 from __future__ import annotations
-from typing import List, Union, Optional, Iterator, Generator
+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
 
@@ -29,12 +30,6 @@ def generate_filepath(sim: SimulationMode, fol: Union[str, Path]) -> Path:
 	return real_path
 
 
-def torus_sites(n: int, w: float, h: float, L: Tuple[int, int]) -> numpy.ndarray:
-	dim = np.array([[w, h]])
-	L = np.array(L)
-	return (np.array([1,1])/2 + np.concatenate([(i*dim*L/n) for i in range(n)])) % dim
-
-
 class DomainParams:
 	"""Container for basic domain parameters
 

squish/ordered.py

@@ -0,0 +1,86 @@
+from __future__ import annotations
+from typing import List, Tuple
+import numpy as np
+from numpy.linalg import norm as mag
+from math import gcd, sqrt, log, tan, atan, pi
+
+Config = Tuple[int, int]
+
+def configurations(domain: DomainParams) -> List[Config]:
+	n, w, h = domain.n, domain.w, domain.h
+	valid = []
+	mults = np.arange(n)
+	configs = np.dstack((np.repeat(mults,n).T, np.tile(mults, n).T))[0][1:]
+	for i in range(len(configs)):
+		eq_x = n if configs[i][0] == 0 else configs[i][0]
+		eq_y = n if configs[i][1] == 0 else configs[i][1]
+
+		if gcd(eq_x, eq_y) != 1:
+			continue
+
+		vecs = configs[i]*np.dstack((w*mults, h*mults)).swapaxes(0,1)/n % domain.dim
+		vmod2 = np.squeeze(np.matmul(vecs, vecs.transpose(0,2,1)))
+		vmodv = np.squeeze(vecs).dot(vecs[1].T).T.flatten()
+
+		if np.all(vmod2 >= vmodv):
+			valid.append(tuple(configs[i]))
+
+	return valid
+
+
+def get_config_generators(domain: DomainParams, config: Config) -> Tuple[Config, Config]:
+	n, w, h = domain.n, domain.w, domain.h
+	q1 = sites(domain, config)
+	v = q1[1]
+	# Sites to check can ignore 0*v and v itself.
+	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.
+	vdot = np.matmul(all_sites, v)
+	in_box = all_sites[np.where((0 <= vdot) & (vdot <= v.dot(v)))[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()
+
+	return tuple(v), tuple(w)
+
+
+def sites(domain: DomainParams, config: Config) -> numpy.ndarray:
+	n, w, h = domain.n, domain.w, domain.h
+	config, mults = np.array(config), np.arange(domain.n)
+	return (config*np.dstack((w*mults, h*mults))[0]/n) % domain.dim
+
+
+def area(domain: DomainParams, config: Config) -> float:
+	v, w = get_config_generators(domain, config)
+	v, w = np.array(v), np.array(w)
+	c = circumcenter(v, w)
+
+	return mag(v)*mag(v/2 - c) + mag(w)*mag(w/2-c) + mag(v-w)*mag((v+w)/2-c)
+
+
+def avg_radius(domain: DomainParams, config: Config) -> float:
+	v, w = get_config_generators(domain, config)
+	v, w = np.array(v), np.array(w)
+	c = circumcenter(v, w)
+
+	return 2*(avg_rp(mag(v), 2*mag(v/2 - c)) + avg_rp(mag(w), 2*mag(w/2-c)) + \
+			 avg_rp(mag(v-w),2*mag((v+w)/2-c)))
+
+
+def avg_rp(d: float, l: float) -> float:
+	return (d/(4*pi))*log(tan(.5*(atan(l/d)+pi/2))**2)
+
+
+def circumcenter(v: numpy.ndarray, w: numpy.ndarray) -> Config:
+	det = 1/(2*rot(v).dot(w))
+	v2, w2 = v.dot(v), w.dot(w)
+	c = np.empty((2,))
+	c[0], c[1] = w[1]*v2 - v[1]*w2,-w[0]*v2 + v[0]*w2
+	return det*c
+
+
+def rot(v: numpy.ndarray) -> numpy.ndarray:
+	w = np.copy(v)
+	w[0], w[1] = -w[1], w[0]
+	return w

squish/simulation.py

@@ -87,6 +87,12 @@ class Simulation:
 			pickle.dump(info, out, pickle.HIGHEST_PROTOCOL)
 
 
+	def save_all(self) -> None:
+		self.save(self.initial_data)
+		for i in range(len(self.frames)):
+			self.save(self.frame_data(i))
+
+
 	def frame_data(self, index: int) -> None:
 		f = self[index]
 		info = {
@@ -145,10 +151,13 @@ class Simulation:
 			if type(sim_class) == str:
 				sim_class = {"flow": Flow, "search": Search, "shrink": Shrink}[sim_class]
 
-
-			sim = sim_class(*all_info[0]["params"], "radial-t", 0,0)
+			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(*frame_info["params"], frame_info["arr"]))
+				frames.append(sim.energy.mode(*frame_info["params"], frame_info["arr"]))
 				#frames[-1].stats = frame_info["stats"]
 
 			sim.frames = frames