Commit a2c29670 authored by Reimar H Leike's avatar Reimar H Leike
Browse files

now using realtive error and cleaned code

parent 35a44a2e
...@@ -74,7 +74,6 @@ def check_value_gradient_consistency(E, tol=1e-8, ntries=100): ...@@ -74,7 +74,6 @@ def check_value_gradient_consistency(E, tol=1e-8, ntries=100):
E2 = _get_acceptable_model(E) E2 = _get_acceptable_model(E)
val = E.value val = E.value
dir = E2.position - E.position dir = E2.position - E.position
# Enext = E2
dirnorm = dir.norm() dirnorm = dir.norm()
for i in range(50): for i in range(50):
Emid = E.at(E.position + 0.5*dir) Emid = E.at(E.position + 0.5*dir)
...@@ -82,20 +81,20 @@ def check_value_gradient_consistency(E, tol=1e-8, ntries=100): ...@@ -82,20 +81,20 @@ def check_value_gradient_consistency(E, tol=1e-8, ntries=100):
dirder = Emid.gradient.vdot(dir)/dirnorm dirder = Emid.gradient.vdot(dir)/dirnorm
else: else:
dirder = Emid.gradient(dir)/dirnorm dirder = Emid.gradient(dir)/dirnorm
numgrad = (E2.value-val)/dirnorm
if isinstance(E, Model): if isinstance(E, Model):
# TODO: use relative error here instead xtol = tol*dirder.norm()
if ((E2.value-val)/dirnorm - dirder).norm()/dirder.norm()< tol: if (abs(numgrad-dirder) < xtol).all():
break break
else: else:
xtol = tol*Emid.gradient_norm xtol = tol*Emid.gradient_norm
if abs((E2.value-val)/dirnorm - dirder) < xtol: if abs(numgrad-dirder) < xtol:
break break
dir *= 0.5 dir *= 0.5
dirnorm *= 0.5 dirnorm *= 0.5
E2 = Emid E2 = Emid
else: else:
raise ValueError("gradient and value seem inconsistent") raise ValueError("gradient and value seem inconsistent")
# E = Enext
def check_value_gradient_curvature_consistency(E, tol=1e-8, ntries=100): def check_value_gradient_curvature_consistency(E, tol=1e-8, ntries=100):
...@@ -105,7 +104,6 @@ def check_value_gradient_curvature_consistency(E, tol=1e-8, ntries=100): ...@@ -105,7 +104,6 @@ def check_value_gradient_curvature_consistency(E, tol=1e-8, ntries=100):
E2 = _get_acceptable_energy(E) E2 = _get_acceptable_energy(E)
val = E.value val = E.value
dir = E2.position - E.position dir = E2.position - E.position
# Enext = E2
dirnorm = dir.norm() dirnorm = dir.norm()
for i in range(50): for i in range(50):
Emid = E.at(E.position + 0.5*dir) Emid = E.at(E.position + 0.5*dir)
...@@ -120,4 +118,3 @@ def check_value_gradient_curvature_consistency(E, tol=1e-8, ntries=100): ...@@ -120,4 +118,3 @@ def check_value_gradient_curvature_consistency(E, tol=1e-8, ntries=100):
E2 = Emid E2 = Emid
else: else:
raise ValueError("gradient, value and curvature seem inconsistent") raise ValueError("gradient, value and curvature seem inconsistent")
# E = Enext
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment