Refactor steepest.py

diff --git a/steepest.py b/steepest.py
index 857cb15e48596d37c2d2bd46f76c94c4c7ee686c..fd8452da10363f77585d032f06d21256134b2998 100644 (file)
--- a/steepest.py
+++ b/steepest.py
@@ -11,7 +11,7 @@ def steepest_iter(x, alpha, grad_f):
     return x - alpha * grad_f(x)
 
 def optimal_step(x, H, b, c, grad_f):
-    pk = -grad_f(x)
+    pk = -grad_f
     Cpk = _C(pk, c)
     Cxk = _C(x, c)
 
@@ -30,54 +30,66 @@ def optimal_step(x, H, b, c, grad_f):
 
     return 0.1
 
-#H = generate.spdmatrix(2, 3)
-#b = generate.vector(2, 0.01, 1)
-#c = generate.vector(2, 0.01, 3)
-
-H = np.array([ (2.06810668,-1.51172618),  (-1.51172618, 3.95629853)])
-b = np.array([-0.09229406, 0.38657143])
-c = np.array([ 0.12659127, 2.42818636])
-
-xiter = np.array([1.5, 1.5])
-x0_norm = norm(gradient_g(xiter, H, b, c))
-
-tolerance = 0.001
-
-relative_residual = norm(gradient_g(xiter, H, b, c)) / x0_norm
-relative_residuals = []
-function_values = []
-
-maxiter = 0
-
-while relative_residual > tolerance and maxiter < 100:
-    maxiter += 1
-
-    alpha = optimal_step(xiter, H, b, c, lambda x: gradient_g(x, H, b, c))
-
-    relative_residual = norm(gradient_g(xiter, H, b, c)) / x0_norm
-    relative_residuals.append(relative_residual)
-    function_values.append(function_g(xiter, H, b, c))
-
-    xiter = steepest_iter(xiter, alpha, lambda x: gradient_g(x, H, b, c))
-
-    z = function_g(xiter, H, b, c)
-    p = mlab.points3d([xiter[0]], [xiter[1]], [z], scale_factor=0.05,
-            color=(maxiter % 2, maxiter % 2, maxiter % 2))
-
-    print xiter
-
-X, Y, Z = surf.surf(H, b, c, xiter[0]-1, xiter[0]+1, xiter[1]-1, xiter[1]+1, -10, 3)
-
-plt.figure(1)
-ax = plt.subplot(211)
-plt.plot(relative_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()
+def steepest(x, H, b, c, tolerance, maxiter, step):
+    xk = np.array(x)
+    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
+ 
+        if step == 0:
+            alpha = optimal_step(xk, H, b, c, gradient_g(xk, H, b, c))
+        else:
+            alpha = step
+    
+        relative_residual = norm(gradient_g(xk, H, b, c)) / x0_norm
+        relative_residuals.append(relative_residual)
+        points.append(xk)
+    
+        xk = steepest_iter(xk, alpha, lambda x: gradient_g(x, H, b, c))
+    
+        z = function_g(xk, H, b, c)
+
+    return (points, relative_residuals)
+    
+def main():
+    H = generate.spdmatrix(2, 3)
+    b = generate.vector(2, 0.01, 1)
+    c = generate.vector(2, 0.01, 3)
+    
+    x0 = np.array([1, 1])
+    
+    function_values = []
+    
+    points, residuals = steepest(x0, H, b, c, 0.01, 100, 0)
+    function_values = []
+
+    for point in points:
+        print point
+        z = function_g(point, H, b, c)
+        function_values.append(z)
+        p = mlab.points3d([point[0]], [point[1]], [z], scale_factor=0.05)
+    
+    last = points[len(points)-1]
+    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()
+
+if __name__ == "__main__":
+    main()