physics/pathfinder
Archived
2
0
Fork 0
Code Issues Pull Requests 2 Projects Releases Wiki Activity

more robust tests and slight tweaks to algorithm

Browse Source
This commit is contained in:
Deepak Mallubhotla 2021-08-04 19:50:54 -05:00
parent 0c8527704d
commit c51c477579
Signed by: deepak
GPG Key ID: 64BF53A3369104E7
2 changed files with 27 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
pathfinder/gradient_descent.py
View File

@ -1,25 +1,29 @@
import numpy 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 desired_cost_squared = desired_cost**2
current = numpy.array(initial) current = numpy.array(initial)
iterations = 0 iterations = 0
curr_cost = numpy.array(cost_array(current)) 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: while iterations < max_iterations:
curr_cost = numpy.array(cost_array(current)) curr_cost = numpy.array(cost_array(current))
if curr_cost.dot(curr_cost) < desired_cost_squared: if curr_cost.dot(curr_cost) < desired_cost_squared:
return ("Finished early", iterations, current) return ("Finished early", iterations, current)
gradient = .5 * numpy.matmul(numpy.transpose(jacobian(current)), curr_cost) 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 tries = step_size_tries
current_step = step_size current_step = step_size
while tries > 0 and next_cost.dot(next_cost) > curr_cost.dot(curr_cost): while total_cost(current - current_step * gradient) > (total_cost(current) - 0.5 * current_step * gradient.dot(gradient)):
current_step = current_step / 10 current_step = current_step * .8
next = current - current_step / 10 * gradient
next_cost = numpy.array(cost_array(next))
tries -= 1 tries -= 1
current = next if tries == 0:
return ("hit minimum step size", iterations, current)
current = current - current_step * gradient
iterations += 1 iterations += 1
return ("Ran out of iterations", iterations, current) return ("Ran out of iterations", iterations, current)

25
tests/test_gradient_descent.py
View File

@ -36,9 +36,11 @@ def test_find_sols():
x, y = pt x, y = pt
return numpy.array([[-2 * x, -2 * y], [-2 * (x - 8), -2 * (y + 8)]]) return numpy.array([[-2 * x, -2 * y], [-2 * (x - 8), -2 * (y + 8)]])
print(costs([3, 4])) message, iterations, result = pathfinder.gradient_descent.find_sols(costs, jac, step_size=0.01, max_iterations=5000, initial=(2, 10), desired_cost=1e-6)
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(
print(result) result, (3, 4),
decimal=7, err_msg='the result was off', verbose=True
)
def dipole_cost(vn, xn_raw): def dipole_cost(vn, xn_raw):
@ -59,11 +61,11 @@ def test_actual_dipole_finding():
p = pt[0:3] p = pt[0:3]
return (p.dot(p) - 35) return (p.dot(p) - 35)
v1 = -0.0554777 v1 = -0.05547767706400186526225414
v2 = -0.0601857 v2 = -0.06018573388098888319642888
v3 = -0.0636403 v3 = -0.06364032191901859480476888
v4 = -0.0648838 v4 = -0.06488383879243851188402150
v5 = -0.0629715 v5 = -0.06297148063759813929659130
# the 0 here is a red herring for index purposes later # the 0 here is a red herring for index purposes later
vns = [0, v1, v2, v3, v4, v5] vns = [0, v1, v2, v3, v4, v5]
@ -106,5 +108,8 @@ def test_actual_dipole_finding():
jac_row(5)(pt), 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) _, _, 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(result) numpy.testing.assert_allclose(
result, (1, 3, 5, 5, 6, 7),
rtol=5e-2, err_msg='the result was off', verbose=True
)