From d3d7f51dd94ccdd9afbab1f42236cf323614fefe Mon Sep 17 00:00:00 2001
From: =?utf8?q?Bj=C3=B8rn=20Rustad?= <rustadbjornen@gmail.com>
Date: Thu, 29 Sep 2011 10:58:05 +0200
Subject: [PATCH] Refactor newton.py

---
 newton.py | 70 +++++++++++++++++++++++++++++++++++++++++--------------
 1 file changed, 53 insertions(+), 17 deletions(-)

diff --git a/newton.py b/newton.py
index 8ea7538..d626dbb 100644
--- a/newton.py
+++ b/newton.py
@@ -1,3 +1,4 @@
+from enthought.mayavi import mlab
 from function import function_g, gradient_g, hessian_g
 import generate
 import surf
@@ -6,27 +7,62 @@ import matplotlib.pyplot as plt
 from numpy.linalg import norm, inv
 
 def newton_iter(x, grad_g, hess_g):
-    return x - inv(hess_g(x)).dot(grad_g(x))
+    return x - inv(hess_g).dot(grad_g)
 
-H = generate.spdmatrix(2, 3)
-b = generate.vector(2, -1, 1)
-c = generate.vector(2, -3, 3)
+def newton(xk, H, b, c, tolerance, maxiter):
+    x0_norm = norm(gradient_g(xk, H, b, c))
+    
+    relative_residual = norm(gradient_g(xk, H, b, c)) / x0_norm
+    relative_residuals = []
+    points = []
+    
+    while relative_residual > tolerance and maxiter > 0:
+        maxiter -= 1
+        xk = newton_iter(xk, gradient_g(xk, H, b, c), hessian_g(xk, H, b, c))
+        relative_residuals.append(relative_residual)
+        points.append(xk)
+        relative_residual = norm(gradient_g(xk, H, b, c)) / x0_norm
 
-alpha = 0.1
-xiter = np.array([1, 1])
-x0_norm = norm(gradient_g(xiter, H, b, c))
+    return (points, relative_residuals)
 
-tolerance = 0.0001
+def main():
+    H = generate.spdmatrix(2, 3)
+    b = generate.vector(2, -1, 1)
+    c = generate.vector(2, -3, 3)
+    
+    print H
+    print b
+    print c
 
-relative_residual = norm(gradient_g(xiter, H, b, c)) / x0_norm
-relative_residuals = [relative_residual]
+    x0 = np.array([0, 0])
+    x0_norm = norm(gradient_g(x0, H, b, c))
 
-while relative_residual > tolerance:
-    xiter = newton_iter(xiter, lambda x: gradient_g(x, H, b, c), lambda x:
-            hessian_g(x, H, b, c))
-    relative_residual = norm(gradient_g(xiter, H, b, c)) / x0_norm
-    relative_residuals.append(relative_residual)
+    points, residuals = newton(x0, H, b, c, 0.01, 100)
+    function_values = []
 
-    print xiter
+    last = points[len(points)-1]
+    for point in points:
+        print point
+        z = function_g(point, H, b, c)
+        function_values.append(z)
+        if abs(point[0]) < 3 and abs(point[1]) < 3:
+            p = mlab.points3d([point[0]], [point[1]], [z], scale_factor=0.05)
+    
+    X, Y, Z = surf.surf(H, b, c, last[0]-1, last[0]+1, last[1]-1, last[1]+1, -10, 3)
+    
+    plt.figure(1)
+    ax = plt.subplot(211)
+    plt.plot(residuals)
+    ax.set_yscale('log')
+    
+    ax = plt.subplot(212)
+    plt.plot(function_values)
+    
+    plt.show()
+    
+    surface = mlab.mesh(X, Y, Z)
+    axes = mlab.axes()
+    mlab.show()
 
-print "Cirka null: ", gradient_g(xiter, H, b, c)
+if __name__ == "__main__":
+    main()
-- 
2.47.3