Implemented more efficient radial-t gradient calculation

.gitignore

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

convert_old.py

@@ -1,22 +0,0 @@
-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()

shrink_energy_comparison.py

@@ -1,14 +1,26 @@
 from __future__ import annotations
 from typing import List
-from simulation import Diagram, Simulation
+import os, argparse, numpy as np, pickle
-import os, argparse, numpy as np
 import matplotlib.pyplot as plt
 from pathlib import Path
 
+from simulation import Diagram, Simulation
+from packsim_core import AreaEnergy, RadialALEnergy, RadialTEnergy
+
+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 get_torus_config_energies(n: int, widths: np.ndarray, h: float, r: float,
-								energy: str) -> Tuple[np.ndarray, np.ndarray]:
+								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 = []
@@ -19,37 +31,36 @@ def get_torus_config_energies(n: int, widths: np.ndarray, h: float, r: float,
 				sim.add_frame(torus=config)
 				configs.append(config)
 
-				eigs = np.sort(np.linalg.eig(sim.frames[-1].hessian(10e-5))[0])
+				# 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 eigs[0] > 1e-4:
-				if zero_eigs != 2:
+				# 	del sim.frames[-1]
-					del sim.frames[-1]
+				# 	del config[-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)}, {2*c}/{2*(n-1)}' + \
+						f'|{"#"*hashes}{" "*(20-hashes)}| {i+1}/{len(widths)}, ' + \
-						f' completed.', flush=True, end='\r')
+						f'{c + (n-1)*j}/{2*(n-1)} completed.', flush=True, end='\r')
 
-		torus_min_energies[i] = min([frame.energy for frame in sim.frames])
+		pair = list(zip(configs,[frame.energy for frame in sim.frames]))
-		torus_max_energies[i] = max([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)
-	return torus_min_energies, torus_max_energies
+
+	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
 
 
-def main():
+def get_equilibria_data(filepath: Path):
-	# Loading arguments.
+	if filepath.is_file():
-	parser = argparse.ArgumentParser("Compiles the equilibriums for each width into a diagram.")
+		with open(filepath, "rb") as data:
-	parser.add_argument('sims_path', metavar='path/to/folder', 
+			return pickle.load(data)
-						help="folder that contains simulation files.")
-	parser.add_argument('-q', '--quiet', dest='quiet', action='store_true', default=False,
-						help="suppress all normal output")
-	parser.add_argument('-o', '--output', dest='output_file')
-
-	args = parser.parse_args()
 
 	sims = []
-	files = list(Path(args.sims_path).iterdir())
+	files = list(Path(filepath).iterdir())
 
 	for i, file in enumerate(files):
 		sims.append(Simulation.load(file))
@@ -60,72 +71,149 @@ 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]
+	pairs_at_widths = []
-	max_frames = [max(sim.frames, key=lambda x: x.energy) for sim in sims]
+	for sim in sims:
+		pairs_at_widths.append([(frame.energy, np.var(frame.stats["avg_radius"]) <= 1e-8) \
+										 for frame in sim.frames])
 
-	min_energies = np.asarray([frame.energy for frame in min_frames])
+	# min_frames = [min(sim.frames, key=lambda x: x.energy) for sim in sims]
-	max_energies = np.asarray([frame.energy for frame in max_frames])
+	# max_frames = [max(sim.frames, key=lambda x: x.energy) for sim in sims]
 
-	torus_min_energies, torus_max_energies = get_torus_config_energies(
+	# min_energies = np.asarray([frame.energy for frame in min_frames])
-		sims[0].n, widths, sims[0].h, sims[0].r, sims[0].energy
+	# 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])
+	with open(sim_file, "wb") as output:
+		pickle.dump(out_tup, output)
+
+	return out_tup
+
+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)
+
+
+def main():
+	# Loading arguments.
+	parser = argparse.ArgumentParser("Compiles the equilibriums for each width into a diagram.")
+	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,
+						help="suppress all normal output")
+
+	args = parser.parse_args()
+
+	widths, energy_shape_tups, 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
 	)
 
-	min_markers = [np.var(frame.stats["avg_radius"]) <= 1e-8 for frame in min_frames]
+	fig_folder = Path(f"figures/ShrinkEnergyComparison - N{n}")
-	max_markers = [np.var(frame.stats["avg_radius"]) <= 1e-8 for frame in max_frames]
+	fig_folder.mkdir(exist_ok=True)
 
 	# Torus minimum energies used as reference.
-
+	plt.tight_layout()
+	# Density of States diagram.
 	fig, ax = plt.subplots(figsize=(16, 8))
-	#ax.plot(widths, nums)
-	ax.plot(widths, [len(sim.frames) for sim in sims])
-	fig.savefig(folder / "")
 
+	distinct_ordered, distinct_unordered = [], []
+	for energy_shapes in energy_shape_tups:
+		equal_shape = list([tup[1] for tup in energy_shapes])
+		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()
+	axis_settings(ax, widths)
+	ax.title.set_text('Density of States')
+	ax.set_xlabel("Width")
+	ax.set_ylabel("Number of States")
+	fig.savefig(fig_folder / "Density Of States.png")
+
+	# Bifurcation diagram
 	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')
-	#ax.plot(widths, torus_max_energies - torus_min_energies, color='C1', linestyle='dotted')
 
-	for i, marker in enumerate(min_markers):
+	ordered_energies, unordered_energies = [], []
-		if marker:
+	for energy_shapes in energy_shape_tups:
-			ax.scatter(widths[i], min_energies[i]-torus_min_energies[i],
+		order_width, unorder_width = [], []
-				marker='H', color="orange", s=20, zorder=4)
+		for pair in energy_shapes:
-		else:
+			if pair[1]:
-			ax.scatter(widths[i], min_energies[i]-torus_min_energies[i],
+				order_width.append(pair[0])
-				marker='d', color="blue", s=20, zorder=4)
+			else:
-
+				unorder_width.append(pair[0])
-	for i, marker in enumerate(max_markers):
+		ordered_energies.append(order_width)
-		if marker:
+		unordered_energies.append(unorder_width)
-			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.invert_xaxis()
+	for i in range(len(torus_min_energies)):
+		ordered_energies[i].append(torus_min_energies[i])
+		ordered_energies[i].append(torus_max_energies[i])
+
+
+	null_unorder = []
+	for i, width in enumerate(unordered_energies):
+		if len(width) == 0:
+			null_unorder.append(i)
+			for x in unordered_energies[i-1]:
+				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])
+
+	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, max_unorder - torus_min_energies, 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],
+			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.set_xlabel("Width")
 	ax.set_ylabel("Reduced Energy")
-	ax.grid(zorder=0)
 
-	ax.set_xticks([round(w,2) for w in widths[::-2]])
+	fig.savefig(fig_folder / "Bifurcation.png")
-	#ax.set_yticks(np.arange(-920, 1120, 40))
-	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"
+	SIM_FOLDER = Path(f"simulations/ShrinkEnergyComparison")
+	SIM_FOLDER.mkdir(exist_ok=True)
 	try:
 		main()
 	except KeyboardInterrupt:

simulation.py

@@ -514,7 +514,7 @@ class Simulation:
 		frames = []
 		with open(filename, 'rb') as data:
 			all_info, sim_class = pickle.load(data)
-
+			sim_class = {"flow": Flow, "search": Search, "shrink": Shrink}[sim_class]
 			sim = sim_class(*all_info[0]["params"], all_info[0]["energy"], 0,0,0,0)
 			for frame_info in all_info:
 				frames.append(sim.energy(*frame_info["params"], frame_info["arr"]))

src/energy.pyx

@@ -142,12 +142,12 @@ cdef class RadialTEnergy(VoronoiContainer):
 
 		cdef Site xi
 		cdef HalfEdge em, e
-		cdef Vector2D Rnla, i2p
+		cdef Vector2D Rnla
 
 		# All energy has a 2pir_0 term.
 		cdef FLOAT_T [:] site_energy = np.full(self.sites.shape[0], TAU*self.r**2)
 		cdef FLOAT_T [:] avg_radii = np.zeros(self.sites.shape[0])
-		cdef FLOAT_T energy, r0, t, tp, B, lntan, cot, cscm, cscp, FA, int_r2d, int_rd
+		cdef FLOAT_T energy, r0, t, tp, B, lntan, csc
 		energy, r0 = 0, self.r
 
 		cdef INT_T i, j
@@ -169,8 +169,7 @@ cdef class RadialTEnergy(VoronoiContainer):
 				else:
 					e.cache.B(&e, <FLOAT_T>acos(<double>(Rnla.y/e.cache.la_mag(&e, NAN))))		
 				
-				i2p = Calc.I2(e, r0, t)
+				e.cache.i2p(&e, Calc.I2(e, r0, t))
-				e.cache.i2p(&e, i2p)
 				e = e.next(&e)
 
 			# For looping again to calculate integrals.
@@ -183,25 +182,18 @@ cdef class RadialTEnergy(VoronoiContainer):
 				lntan = <FLOAT_T>(log(fabs(tan(<double>((tp+B)/2))))) - \
 							<FLOAT_T>(log(fabs(tan(<double>((t+B)/2)))))
 
-				cot = -1/(<FLOAT_T>(tan(<double>(tp+B)))) + \
+				csc = 1/(<FLOAT_T>(sin(<double>(tp+B)))) - \
-						1/(<FLOAT_T>(tan(<double>(t+B))))	
+						1/(<FLOAT_T>(sin(<double>(t+B))))
-
-				cscm, cscp = 1/(<FLOAT_T>(sin(<double>(t+B)))), \
-								1/(<FLOAT_T>(sin(<double>(tp+B))))
 				
 				em.cache.lntan(&em, lntan)
-				em.cache.cot(&em, cot)
+				em.cache.csc(&em, csc)
-				em.cache.csc(&em, cscp-cscm)
-				em.cache.csc2(&em, cscp**2 - cscm**2)
-				FA = (em.cache.F(&em, NAN)/em.cache.la_mag(&em, NAN))
 
-				int_r2d, int_rd = FA**2*cot, FA*lntan
+				avg_radii[i] += (em.cache.F(&em, NAN)/em.cache.la_mag(&em, NAN))*lntan
-
-				avg_radii[i] += int_rd
-				site_energy[i] += int_r2d - 2*r0*int_rd
 
 				em = em.next(&em)
 
+			site_energy[i] += 2*(xi.cache.area(&xi, NAN) - r0*avg_radii[i])
+
 			xi.cache.avg_radius(&xi, avg_radii[i]/TAU)
 			xi.cache.energy(&xi, site_energy[i])
 			if i < self.n:
@@ -266,13 +258,7 @@ cdef class Calc:
 	cdef inline Vector2D I2(HalfEdge e, FLOAT_T r0, FLOAT_T t) nogil:
 		cdef Vector2D Rda = e.cache.da(&e, NAN_VECTOR)
 		Rda = Rda.rot(&Rda)
-
+		Rda.self.sdiv(&Rda, e.cache.da_mag(&e, NAN))
-		cdef Vector2D Rcircle = init.Vector2D(
-			-<FLOAT_T>sin(<double>t), <FLOAT_T>cos(<double>t)
-		)
-		cdef FLOAT_T p = e.cache.F(&e, NAN) / Rcircle.dot(&Rcircle, e.cache.la(&e, NAN_VECTOR))
-		p = ((p - r0)**2)/(Rda.dot(&Rda, Rda))
-		Rda.self.smul(&Rda, p)
 
 		return Rda
 
@@ -283,21 +269,20 @@ cdef class Calc:
 		cdef Vector2D Rda, i2ps, fp, gterms, q
 		cdef Matrix2x2 ha, hap, hdiff
 
-		cdef FLOAT_T t1, t2, lntan, cot, csc, csc2, sinB, cosB, sinBp, cosBp, F, A, B
+		cdef FLOAT_T t1, t2, lntan, csc, sinB, cosB, sinBp, cosBp, F, A, B
 
 		xe = e.face(&e)
 		ep = e.next(&e)
 		F, A, B = e.cache.F(&e, NAN), e.cache.la_mag(&e, NAN), e.cache.B(&e, NAN)
 		t1, t2 = e.cache.phi(&e, NAN), ep.cache.phi(&ep, NAN)
-		lntan, cot, csc, csc2 = e.cache.lntan(&e, NAN), e.cache.cot(&e, NAN), \
+
-								e.cache.csc(&e, NAN), e.cache.csc2(&e, NAN)
+		lntan, csc = e.cache.lntan(&e, NAN), e.cache.csc(&e, NAN)
 
 		sinB, cosB = <FLOAT_T>(sin(<double>(B))), <FLOAT_T>(cos(<double>(B)))
 		sinBp, cosBp = <FLOAT_T>(sin(<double>(B-PI_2))), \
 						<FLOAT_T>(cos(<double>(B-PI_2)))
 
 
-
 		ha, hap = e.get_H(&e, xi), ep.get_H(&ep, xi)
 		hdiff = hap.copy.msub(&hap, ha)
 		# If edge is part of differentiated site.
@@ -319,25 +304,20 @@ cdef class Calc:
 		fp = e.cache.la(&e, NAN_VECTOR)
 		fp.self.matmul(&fp, R.copy.matmul(&R, ha))
 		fp.self.vadd(&fp, Rda.copy.matmul(&Rda, hdiff))
-		fp.self.smul(&fp, (F/A**2)*cot - (r0/A)*lntan)
+		fp.self.smul(&fp, lntan/A)
 		
 		gterms = init.Vector2D(
 			cosBp*lntan + sinBp*csc,
 			cosB*lntan + sinB*csc
 		)
-		gterms.self.smul(&gterms, r0*F/A**2)
+		gterms.self.smul(&gterms, -F/A**2)
 
-		q = init.Vector2D(
-			0.5*sinBp*csc2 + cosBp*cot,
-			0.5*sinB*csc2 + cosB*cot
-		)
-		q.self.smul(&q, -F**2/A**3)
-
-		gterms.self.vadd(&gterms, q)
 		gterms = gterms.rot(&gterms)
 		gterms.self.matmul(&gterms, hdiff)
 
 		fp.self.vadd(&fp, gterms)
-		fp.self.smul(&fp, 2)
 
-		return i2ps.copy.vadd(&i2ps, fp)
+		i2ps.self.vadd(&i2ps, fp)
+		i2ps.self.smul(&i2ps, -2*r0)
+
+		return i2ps

src/packsim_core.pxd

@@ -14,7 +14,7 @@ cdef struct Init:
 	Matrix2x2 (*Matrix2x2)(FLOAT_T, FLOAT_T, FLOAT_T, FLOAT_T) nogil
 	SiteCacheMap (*SiteCacheMap)(INT_T, INT_T, INT_T, INT_T, INT_T) nogil
 	EdgeCacheMap (*EdgeCacheMap)(INT_T, INT_T, INT_T, INT_T, INT_T, INT_T, INT_T, INT_T,
-					INT_T, INT_T, INT_T, INT_T, INT_T, INT_T, INT_T, INT_T) nogil
+					INT_T, INT_T, INT_T, INT_T, INT_T, INT_T) nogil
 	VoronoiInfo (*VoronoiInfo)(INT_T [:, ::1], INT_T[:, ::1], FLOAT_T[:, ::1], 
 								FLOAT_T[:, ::1], FLOAT_T[:, ::1], FLOAT_T[:, ::1],
 								EdgeCacheMap*) nogil
@@ -156,7 +156,7 @@ ctypedef struct SiteCacheMap:
 # Psuedo-class that handles caching for edges.
 ctypedef struct EdgeCacheMap:
 	INT_T iH, ila, ila_mag, ida, ida_mag, ixij, idVdv, iphi, iB, iF, ii2p,\
-			ilntan, icot, icsc, icsc2, size
+			ilntan, icsc, size
 
 	Matrix2x2 (*H)(HalfEdge*, Matrix2x2) nogil
 
@@ -172,9 +172,7 @@ ctypedef struct EdgeCacheMap:
 	FLOAT_T (*B)(HalfEdge*, FLOAT_T) nogil
 	FLOAT_T (*F)(HalfEdge*, FLOAT_T) nogil	
 	FLOAT_T (*lntan)(HalfEdge*, FLOAT_T) nogil
-	FLOAT_T (*cot)(HalfEdge*, FLOAT_T) nogil
 	FLOAT_T (*csc)(HalfEdge*, FLOAT_T) nogil
-	FLOAT_T (*csc2)(HalfEdge*, FLOAT_T) nogil
 
 # Psuedo-class to just contain all pertaining info for sites and edges.
 ctypedef struct VoronoiInfo:

src/voronoi_dcel.pyx

@@ -8,9 +8,9 @@ init.SiteCacheMap, init.EdgeCacheMap, init.VoronoiInfo, init.Site, init.HalfEdge
 cdef SiteCacheMap SITE_CACHE_MAP = init.SiteCacheMap(0, 1, 2, 3, 4)
 
 cdef EdgeCacheMap AREA_EDGE_CACHE_MAP = init.EdgeCacheMap(0, 4, 6, 8, 10, -1, 12, 13,
-															-1, -1, -1, -1, -1, -1, -1, 14)
+															-1, -1, -1, -1, -1, 14)
 cdef EdgeCacheMap RADIALT_EDGE_CACHE_MAP = init.EdgeCacheMap(0, 4, 6, 8, -1, 10, 12, 13,
-																14, 15, 16, 17, 18, 19, 20, 21)
+																14, 15, 16, 17, 18, 19)
 
 #### SiteCacheMap Methods ####
 
@@ -81,16 +81,15 @@ cdef inline FLOAT_T avg_radius(Site* self, FLOAT_T val) nogil:
 
 cdef inline EdgeCacheMap init_edgecachemap(INT_T iH, INT_T ila, INT_T ida, INT_T ixij, 
 				INT_T idVdv, INT_T ii2p, INT_T ila_mag, INT_T ida_mag, INT_T iphi, INT_T iB,
-				INT_T iF, INT_T ilntan, INT_T icot, INT_T icsc, INT_T icsc2, INT_T size) nogil:
+				INT_T iF, INT_T ilntan, INT_T icsc, INT_T size) nogil:
 	cdef EdgeCacheMap ec
 	ec.iH, ec.ila, ec.ida, ec.ixij, ec.idVdv, ec.ii2p, ec.ila_mag, ec.ida_mag, ec.iphi, \
-		ec.iB, ec.iF, ec.ilntan, ec.icot, ec.icsc, ec.icsc2 = iH, ila, ida, ixij, idVdv, ii2p, \
+		ec.iB, ec.iF, ec.ilntan, ec.icsc = iH, ila, ida, ixij, idVdv, ii2p, \
-					 ila_mag, ida_mag, iphi, iB, iF, ilntan, icot, icsc, icsc2
+					 ila_mag, ida_mag, iphi, iB, iF, ilntan, icsc
 	ec.size = size
 
 	ec.H, ec.la, ec.da, ec.xij, ec.dVdv, ec.i2p, ec.la_mag, ec.da_mag, ec.phi, ec.B, ec.F, \
-		 ec.lntan, ec.cot, ec.csc, ec.csc2 = H, la, da, xij, dVdv, i2p, la_mag, da_mag, phi, \
+		 ec.lntan, ec.csc = H, la, da, xij, dVdv, i2p, la_mag, da_mag, phi, B, F, lntan, csc
-		 										B, F, lntan, cot, csc, csc2
 
 	return ec
 
@@ -267,16 +266,6 @@ cdef inline FLOAT_T lntan(HalfEdge* self, FLOAT_T val) nogil:
 			(self.arr_index, self.cache.ilntan), val)
 		return val
 
-cdef inline FLOAT_T cot(HalfEdge* self, FLOAT_T val) nogil:
-	if isnan(<double>val):
-		return self.info.edge_cache.get(&self.info.edge_cache,
-			(self.arr_index, self.cache.icot)
-		)
-	else:
-		self.info.edge_cache.set(&self.info.edge_cache,
-			(self.arr_index, self.cache.icot), val)
-		return val
-
 cdef inline FLOAT_T csc(HalfEdge* self, FLOAT_T val) nogil:
 	if isnan(<double>val):
 		return self.info.edge_cache.get(&self.info.edge_cache,
@@ -287,16 +276,6 @@ cdef inline FLOAT_T csc(HalfEdge* self, FLOAT_T val) nogil:
 			(self.arr_index, self.cache.icsc), val)
 		return val
 
-cdef inline FLOAT_T csc2(HalfEdge* self, FLOAT_T val) nogil:
-	if isnan(<double>val):
-		return self.info.edge_cache.get(&self.info.edge_cache,
-			(self.arr_index, self.cache.icsc2)
-		)
-	else:
-		self.info.edge_cache.set(&self.info.edge_cache,
-			(self.arr_index, self.cache.icsc2), val)
-		return val
-
 
 #### VoronoiInfo Methods ####