Removes gradient descent and begins refactor

pathfinder/gradient_descent.py

@@ -1,29 +0,0 @@
-import numpy
-
-
-def find_sols(cost_array, jacobian, step_size=0.001, max_iterations=5, initial=(0, 0), desired_cost=1e-6, step_size_tries=30):
-	desired_cost_squared = desired_cost**2
-	current = numpy.array(initial)
-	iterations = 0
-	curr_cost = numpy.array(cost_array(current))
-
-	def total_cost(x):
-		cost = numpy.array(cost_array(x))
-		return cost.dot(cost)
-
-	while iterations < max_iterations:
-		curr_cost = numpy.array(cost_array(current))
-		if curr_cost.dot(curr_cost) < desired_cost_squared:
-			return ("Finished early", iterations, current)
-		gradient = .5 * numpy.matmul(numpy.transpose(jacobian(current)), curr_cost)
-
-		tries = step_size_tries
-		current_step = step_size
-		while total_cost(current - current_step * gradient) > (total_cost(current) - 0.5 * current_step * gradient.dot(gradient)):
-			current_step = current_step * .8
-			tries -= 1
-			if tries == 0:
-				return ("hit minimum step size", iterations, current)
-		current = current - current_step * gradient
-		iterations += 1
-	return ("Ran out of iterations", iterations, current)

pathfinder/model/__init__.py

@@ -1,4 +1,5 @@
-from pathfinder.model.dot import Dot
+from pathfinder.model.dot import DotMeasurement
 from pathfinder.model.model import DotDipoleModel
+from pathfinder.model.staticdipole import StaticDipole
 
-__all__ = ['Dot', 'DotDipoleModel', ]
+__all__ = ['DotMeasurement', 'DotDipoleModel', 'StaticDipole']

pathfinder/model/dot.py

@@ -2,9 +2,11 @@ from dataclasses import dataclass
 import numpy
 import numpy.typing
 
+from pathfinder.model.staticdipole import StaticDipole
+
 
 @dataclass
-class Dot():
+class DotMeasurement():
 	'''
 	Representation of a dot measuring static dipoles.
 
@@ -13,20 +15,32 @@ class Dot():
 	v : float
 		The voltage measured at the dot.
 	r : numpy.ndarray
-		The number of dots used to sample the potential.
+		The position of the dot.
 	'''
 	v: float
-	r: numpy.typing.ArrayLike
+	r: numpy.ndarray
 
 	def __post_init__(self) -> None:
 		self.r = numpy.array(self.r)
 
 	def v_for_point(self, pt: numpy.ndarray) -> float:
+		'''
+		Returns the voltage for a static dipole represented as a phase space point.
+
+		Parameters
+		----------
+		pt: A length 6 numpy array representing a dipole, as (px, py, pz, sx, sy, sz).
+
+		Returns
+		----------
+		Returns the voltage from that dipole at the point.
+		'''
 		p = pt[0:3]  # hardcoded here because chances
 		s = pt[3:6]  # are we'll only ever work in 3d.
 
-		diff = self.r - s
+		dipole = StaticDipole(p, s)
-		return p.dot(diff) / (numpy.linalg.norm(diff)**3)
+
+		return dipole.v_at_position(self.r)
 
 	def cost(self, pts: numpy.ndarray) -> float:
 		# 6 because dipole in 3d has 6 degrees of freedom.

pathfinder/model/model.py

@@ -1,7 +1,7 @@
 from typing import Callable, Sequence
 import numpy
 
-from pathfinder.model.dot import Dot
+from pathfinder.model.dot import DotMeasurement
 
 
 class DotDipoleModel():
@@ -16,7 +16,7 @@ class DotDipoleModel():
 	n: int
 		The number of dipoles to assume.
 	'''
-	def __init__(self, dots: Sequence[Dot], n: int) -> None:
+	def __init__(self, dots: Sequence[DotMeasurement], n: int) -> None:
 		self.dots = dots
 		self.m = len(dots)
 		self.n = n
@@ -29,3 +29,9 @@ class DotDipoleModel():
 			return numpy.array([dot.cost(pt) for dot in self.dots])
 
 		return costs_to_return
+
+	def jac(self) -> Callable[[numpy.ndarray], numpy.ndarray]:
+		def jac_to_return(pts: numpy.ndarray) -> numpy.ndarray:
+			return numpy.array([dot.jac(pts) for dot in self.dots])
+
+		return jac_to_return

pathfinder/model/staticdipole.py

@@ -0,0 +1,31 @@
+from dataclasses import dataclass
+import numpy
+import numpy.typing
+
+
+@dataclass
+class StaticDipole():
+	'''
+	Representation of a static dipoles, either known or guessed.
+
+	Parameters
+	----------
+	p : numpy.ndarray
+		The static dipole moment.
+	s : numpy.ndarray
+		The position of the dipole.
+	'''
+	p: numpy.ndarray
+	s: numpy.ndarray
+
+	def __post_init__(self) -> None:
+		'''
+			Coerce the inputs into numpy arrays.
+		'''
+		self.p = numpy.array(self.p)
+		self.s = numpy.array(self.s)
+
+	def v_at_position(self, r: numpy.typing.ArrayLike) -> float:
+
+		diff = r - self.s
+		return self.p.dot(diff) / (numpy.linalg.norm(diff)**3)

tests/model/test_dot.py

@@ -5,14 +5,14 @@ import pathfinder.model as model
 
 
 def test_dot():
-	dot = model.Dot(0.235, (1, 2, 3))
+	dot = model.DotMeasurement(0.235, (1, 2, 3))
 	assert dot.v == 0.235
 	numpy.testing.assert_array_equal(dot.r, (1, 2, 3), "These arrays should have been equal!")
 
 
 def test_dot_v_from_dipole():
 	# for a dot located at (1, 2, 3)
-	dot = model.Dot(50, (1, 2, 3))
+	dot = model.DotMeasurement(50, (1, 2, 3))
 
 	# and dipole located at (4, 7, 11) with p=(8, 9, 10)
 	pt = numpy.array((8, 9, 10, 4, 7, 11))
@@ -27,7 +27,7 @@ def test_dot_v_from_dipole():
 
 def test_dot_jac():
 	# for a dot located at (1, 2, 3)
-	dot = model.Dot(50, (1, 2, 3))
+	dot = model.DotMeasurement(50, (1, 2, 3))
 
 	# and dipole located at (4, 7, 11) with p=(8, 9, 10)
 	pt = numpy.array((8, 9, 10, 4, 7, 11))

tests/model/test_model.py

@@ -1,4 +1,7 @@
+import numpy
 import pathfinder.model as model
+import scipy.optimize
+import pytest
 
 
 def test_dotdipolemodel_repr():
@@ -7,5 +10,98 @@ def test_dotdipolemodel_repr():
 
 
 def test_dotdipolemodel_m():
-	mod = model.DotDipoleModel([model.Dot(1, (0, 0, 0)), model.Dot(2, (0, 0, 0))], 1)
+	mod = model.DotDipoleModel([model.DotMeasurement(1, (0, 0, 0)), model.DotMeasurement(2, (0, 0, 0))], 1)
 	assert mod.m == 2
+
+
+def print_result(msg, result):
+	print(msg)
+	print(f"\tResult: {result.x}")
+	print(f"\tSuccess: {result.success}. {result.message}")
+	try:
+		print(f"\tFunc evals: {result.nfev}")
+	except AttributeError:
+		pass
+	try:
+		print(f"\tJacb evals: {result.njev}")
+	except AttributeError:
+		pass
+
+
+@pytest.mark.skip(reason="old")
+def test_dot_dipole_model_jac():
+	v1 = -0.05547767706400186526225414
+	v2 = -0.06018573388098888319642888
+	v3 = -0.06364032191901859480476888
+	v4 = -0.06488383879243851188402150
+	v5 = -0.06297148063759813929659130
+	v6 = -0.05735489606460216
+	v7 = -0.07237320672886623
+	v8 = -0.1082531754730548
+
+	c1 = model.DotMeasurement(v1, [0, 0, 1])
+	c2 = model.DotMeasurement(v2, [0, 0, 2])
+	c3 = model.DotMeasurement(v3, [0, 0, 3])
+	c4 = model.DotMeasurement(v4, [0, 0, 4])
+	c5 = model.DotMeasurement(v5, [0, 0, 5])
+	c6 = model.DotMeasurement(v6, [0, 0, 6])
+	c7 = model.DotMeasurement(v7, [1, 1, 7])
+	c8 = model.DotMeasurement(v8, [1, 2, 3])
+
+	mod = model.DotDipoleModel([c1, c2, c3, c4, c5, c6], 1)
+	res = scipy.optimize.least_squares(mod.costs(), numpy.array([1, 2, 3, 4, 5, 6]), jac=mod.jac(), ftol=1e-12)
+	print_result("6 dots, least sq", res)
+	mod2 = model.DotDipoleModel([c1, c2, c3, c4, c5, c6, c7, c8], 1)
+	res2 = scipy.optimize.least_squares(mod2.costs(), numpy.array([0, 0, 0, 0, 0, 0]), jac=mod2.jac(), ftol=1e-12, gtol=3e-16)
+	print_result("7 dots, least squares", res2)
+	print(mod2.costs()(res2.x))
+	print(mod2.costs()(numpy.array([1, 3, 5, 5, 6, 7])))
+
+
+@pytest.mark.skip(reason="bad test")
+def test_dot_2dipoles_model_jac():
+	dots_12andone = [
+		model.DotMeasurement(-0.1978319326584865, [-4, -1, 0]),
+		model.DotMeasurement(-0.1273171638293727, [-4, -5, 0]),
+		model.DotMeasurement(-0.05545025617224288, [-4, -9, 0]),
+		model.DotMeasurement(-0.4960209997369774, [-1, -1, 0]),
+		model.DotMeasurement(-0.1763373289754278, [-1, -5, 0]),
+		model.DotMeasurement(-0.04946346672578462, [-1, -9, 0]),
+		model.DotMeasurement(-0.5633156098386561, [1, -1, 0]),
+		model.DotMeasurement(-0.113765134433174, [1, -5, 0]),
+		model.DotMeasurement(-0.0294893572499722, [1, -9, 0]),
+		model.DotMeasurement(0.7794941616360612, [4, -1, 0]),
+		model.DotMeasurement(0.1110683086477768, [4, -5, 0]),
+		model.DotMeasurement(0.01183272220840589, [4, -9, 0]),
+		model.DotMeasurement(-0.1096485462119833, [1, 1, 1]),
+		model.DotMeasurement(-0.3925851888077783, [1, 1, -1])
+	]
+	# dots = [
+	# 	model.Dot(-0.1978319326584865, [-4, -1, 0]),
+	# 	model.Dot(-0.1273171638293727, [-4, -5, 0]),
+	# 	model.Dot(-0.05545025617224288, [-4, -9, 0]),
+	# 	model.Dot(-0.4960209997369774, [-1, -1, 0]),
+	# 	model.Dot(-0.1763373289754278, [-1, -5, 0]),
+	# 	model.Dot(-0.04946346672578462, [-1, -9, 0]),
+	# 	model.Dot(-0.5633156098386561, [1, -1, 0]),
+	# 	model.Dot(-0.113765134433174, [1, -5, 0]),
+	# 	model.Dot(-0.0294893572499722, [1, -9, 0]),
+	# 	model.Dot(0.7794941616360612, [4, -1, 0]),
+	# 	model.Dot(0.1110683086477768, [4, -5, 0]),
+	# 	model.Dot(0.01183272220840589, [4, -9, 0]),
+	# 	model.Dot(-0.0261092728062841, [-4, -13, 0]),
+	# 	model.Dot(-0.02035559593904894, [-1, -13, 0]),
+	# 	model.Dot(-0.01296894715810522, [1, -13, 0]),
+	# 	model.Dot(0.0003306001626435171, [4, -13, 0]),
+	# 	model.Dot(0.2612817068810759, [7, -1, 0]),
+	# 	model.Dot(0.1203841445911355, [7, -5, 0]),
+	# 	model.Dot(0.03425933872931543, [7, -9, 0]),
+	# 	model.Dot(0.01068547688208644, [7, -13, 0]),
+	# ]
+
+	mod = model.DotDipoleModel(dots_12andone, 2)
+	res = scipy.optimize.least_squares(mod.costs(), numpy.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), jac=mod.jac(), ftol=1e-12)
+	print_result("6 dots, least sq", res)
+	print(mod.costs()(res.x))
+	for val in res.x:
+		print(val)

tests/test_gradient_descent.py

@@ -1,118 +0,0 @@
-import numpy
-import pathfinder.gradient_descent
-
-
-def circ_cost(radius, center=(0, 0)):
-
-	def cf(pt):
-		pt2 = numpy.array(pt) - numpy.array(center)
-		return (radius**2 - pt2.dot(pt2))
-
-	return cf
-
-
-def test_circ_cost():
-	cost = circ_cost(5)
-	actual = cost([3, 4])
-	expected = 0
-	assert actual == expected
-
-	cost = circ_cost(13, [12, 5])
-	actual = cost([0, 0])
-	expected = 0
-	assert actual == expected
-
-
-def test_find_sols():
-	c1 = circ_cost(5)
-	c2 = circ_cost(13, [8, -8])
-
-	def costs(pt):
-		return numpy.array(
-			[c1(pt), c2(pt)]
-		)
-
-	def jac(pt):
-		x, y = pt
-		return numpy.array([[-2 * x, -2 * y], [-2 * (x - 8), -2 * (y + 8)]])
-
-	message, iterations, result = pathfinder.gradient_descent.find_sols(costs, jac, step_size=0.01, max_iterations=5000, initial=(2, 10), desired_cost=1e-6)
-	numpy.testing.assert_almost_equal(
-		result, (3, 4),
-		decimal=7, err_msg='the result was off', verbose=True
-	)
-
-
-def dipole_cost(vn, xn_raw):
-	xn = numpy.array(xn_raw)
-
-	def dc(pt):
-		p = pt[0:3]
-		s = pt[3:6]
-
-		diff = xn - s
-		return (vn * (numpy.linalg.norm(diff)**3)) - p.dot(diff)
-
-	return dc
-
-
-def test_actual_dipole_finding():
-	def c0(pt):
-		p = pt[0:3]
-		return (p.dot(p) - 35)
-
-	v1 = -0.05547767706400186526225414
-	v2 = -0.06018573388098888319642888
-	v3 = -0.06364032191901859480476888
-	v4 = -0.06488383879243851188402150
-	v5 = -0.06297148063759813929659130
-
-	# the 0 here is a red herring for index purposes later
-	vns = [0, v1, v2, v3, v4, v5]
-	# the 0 here is a red herring
-	xns = [numpy.array([0, 0, n]) for n in range(0, 6)]
-
-	c1 = dipole_cost(v1, [0, 0, 1])
-	c2 = dipole_cost(v2, [0, 0, 2])
-	c3 = dipole_cost(v3, [0, 0, 3])
-	c4 = dipole_cost(v4, [0, 0, 4])
-	c5 = dipole_cost(v5, [0, 0, 5])
-
-	def costs(pt):
-		return numpy.array(
-			[c0(pt), c1(pt), c2(pt), c3(pt), c4(pt), c5(pt)]
-		)
-
-	def jac_row(n):
-		def jr(pt):
-			p = pt[0:3]
-			s = pt[3:6]
-			vn = vns[n]
-			xn = xns[n]
-			diff = xn - s
-			return [
-				-diff[0], -diff[1], -diff[2],
-				p[0] - vn * 3 * numpy.linalg.norm(diff) * (diff)[0],
-				p[1] - vn * 3 * numpy.linalg.norm(diff) * (diff)[1],
-				p[2] - vn * 3 * numpy.linalg.norm(diff) * (diff)[2]
-			]
-		return jr
-
-	def jac(pt):
-		return numpy.array([
-			[2 * pt[0], 2 * pt[1], 2 * pt[2], 0, 0, 0],
-			jac_row(1)(pt),
-			jac_row(2)(pt),
-			jac_row(3)(pt),
-			jac_row(4)(pt),
-			jac_row(5)(pt),
-		])
-
-	_, iterations, result = pathfinder.gradient_descent.find_sols(costs, jac, step_size=1, max_iterations=10000, initial=(1, 2, 3, 4, 5, 6), desired_cost=1e-6, step_size_tries=30)
-	print(f"\n\n{iterations} iterations")
-	print(result)
-	print("\n")
-	numpy.testing.assert_allclose(
-		result, (1, 3, 5, 5, 6, 7),
-		rtol=5e-2, err_msg='the result was off', verbose=True
-	)