diff --git a/pathfinder/gradient_descent.py b/pathfinder/gradient_descent.py index 4718a25..3aabf87 100644 --- a/pathfinder/gradient_descent.py +++ b/pathfinder/gradient_descent.py @@ -1,25 +1,29 @@ 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=10): +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) - next = current - step_size * gradient - next_cost = numpy.array(cost_array(next)) + tries = step_size_tries current_step = step_size - while tries > 0 and next_cost.dot(next_cost) > curr_cost.dot(curr_cost): - current_step = current_step / 10 - next = current - current_step / 10 * gradient - next_cost = numpy.array(cost_array(next)) + while total_cost(current - current_step * gradient) > (total_cost(current) - 0.5 * current_step * gradient.dot(gradient)): + current_step = current_step * .8 tries -= 1 - current = next + 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/tests/test_gradient_descent.py b/tests/test_gradient_descent.py index dc7406e..849b8a7 100644 --- a/tests/test_gradient_descent.py +++ b/tests/test_gradient_descent.py @@ -36,9 +36,11 @@ def test_find_sols(): x, y = pt return numpy.array([[-2 * x, -2 * y], [-2 * (x - 8), -2 * (y + 8)]]) - print(costs([3, 4])) - result = pathfinder.gradient_descent.find_sols(costs, jac, step_size=0.01, max_iterations=5000, initial=(2, 10), desired_cost=1e-6) - print(result) + 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): @@ -59,11 +61,11 @@ def test_actual_dipole_finding(): p = pt[0:3] return (p.dot(p) - 35) - v1 = -0.0554777 - v2 = -0.0601857 - v3 = -0.0636403 - v4 = -0.0648838 - v5 = -0.0629715 + 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] @@ -106,5 +108,8 @@ def test_actual_dipole_finding(): jac_row(5)(pt), ]) - result = pathfinder.gradient_descent.find_sols(costs, jac, step_size=0.01, max_iterations=10000, initial=(1, 2, 3, 4, 5, 6), desired_cost=1e-6, step_size_tries=25) - print(result) + _, _, 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) + numpy.testing.assert_allclose( + result, (1, 3, 5, 5, 6, 7), + rtol=5e-2, err_msg='the result was off', verbose=True + )