From c279f1b470178a5953ab5beddd4dffcd5daf9ee7 Mon Sep 17 00:00:00 2001 From: Deepak Date: Sun, 29 Aug 2021 10:54:36 -0500 Subject: [PATCH] Removes gradient descent and begins refactor --- pathfinder/gradient_descent.py | 29 -------- pathfinder/model/__init__.py | 5 +- pathfinder/model/dot.py | 24 +++++-- pathfinder/model/model.py | 10 ++- pathfinder/model/staticdipole.py | 31 ++++++++ tests/model/test_dot.py | 6 +- tests/model/test_model.py | 98 ++++++++++++++++++++++++- tests/test_gradient_descent.py | 118 ------------------------------- 8 files changed, 161 insertions(+), 160 deletions(-) delete mode 100644 pathfinder/gradient_descent.py create mode 100644 pathfinder/model/staticdipole.py delete mode 100644 tests/test_gradient_descent.py diff --git a/pathfinder/gradient_descent.py b/pathfinder/gradient_descent.py deleted file mode 100644 index 3aabf87..0000000 --- a/pathfinder/gradient_descent.py +++ /dev/null @@ -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) diff --git a/pathfinder/model/__init__.py b/pathfinder/model/__init__.py index 662e43b..65ff72d 100644 --- a/pathfinder/model/__init__.py +++ b/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'] diff --git a/pathfinder/model/dot.py b/pathfinder/model/dot.py index 9e6ed6d..a59b23d 100644 --- a/pathfinder/model/dot.py +++ b/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 - return p.dot(diff) / (numpy.linalg.norm(diff)**3) + dipole = StaticDipole(p, s) + + return dipole.v_at_position(self.r) def cost(self, pts: numpy.ndarray) -> float: # 6 because dipole in 3d has 6 degrees of freedom. diff --git a/pathfinder/model/model.py b/pathfinder/model/model.py index 6584fe5..ac246b6 100644 --- a/pathfinder/model/model.py +++ b/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 diff --git a/pathfinder/model/staticdipole.py b/pathfinder/model/staticdipole.py new file mode 100644 index 0000000..0e9c09a --- /dev/null +++ b/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) diff --git a/tests/model/test_dot.py b/tests/model/test_dot.py index 858d958..837fc49 100644 --- a/tests/model/test_dot.py +++ b/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)) diff --git a/tests/model/test_model.py b/tests/model/test_model.py index 6a6e1c7..8eaf597 100644 --- a/tests/model/test_model.py +++ b/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) diff --git a/tests/test_gradient_descent.py b/tests/test_gradient_descent.py deleted file mode 100644 index 47ee030..0000000 --- a/tests/test_gradient_descent.py +++ /dev/null @@ -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 - )