Prettify

diff --git a/function.py b/function.py
index 5c9cd9fbc3fdca3cebfe2a93affc2bdd7715004f..a7a26a307cec83661e877477b27f4421cb2d6fc7 100644 (file)
--- a/function.py
+++ b/function.py
@@ -19,6 +19,7 @@ def gradient_g(x, H, b, c):
 def hessian_g(x, H, b, c):
     return 0.5 * H.T + 0.5 * H + _C(x,c)
 
+# Small test of the above functions
 def main():
     A = np.array([ (1, 2), (2, 3) ])
     x = np.array([5, 4])

diff --git a/generate.py b/generate.py
index 4af4c0b689d450ee8ad86771437b32130d09db55..7bcef1e109741671430460711fd7adcda33660d1 100644 (file)
--- a/generate.py
+++ b/generate.py
@@ -2,6 +2,7 @@ import numpy as np
 import math
 import random
 
+# Generate nxn symmetric positive definite matrix
 def spdmatrix(n, limit):
     A = np.zeros( (n,n) )
     for y in xrange(0, n):
@@ -10,12 +11,15 @@ def spdmatrix(n, limit):
             A[y][x] = random.uniform(-limit, limit)
     return A.T.dot(A)
 
+# Generate n-dim vector with elements uniformly distributed
+# between lower and upper limit
 def vector(n, lower, upper):
     X = np.zeros(n)
     for x in xrange(0, n):
         X[x] = random.uniform(lower, upper)
     return X
 
+# Test the functions above
 def main():
     A = spdmatrix(5, 10)
     print "SPD-matrise:"

diff --git a/main.py b/main.py
index df0e8820b205cc78c5541c24d01d5b7112224138..8d2dd7a47f959c4bdaf2e1a94dcdd2179ee94a6f 100644 (file)
--- a/main.py
+++ b/main.py
@@ -6,7 +6,7 @@ import surf
 from steepest import steepest
 from function import function_g
 from newton import newton
-from newtsteep import newtsteep
+from steepestnewton import steepestnewton
 
 def main():
     H = np.array([ ( 1.8580682, -1.01197101), (-1.01197101, 0.63456543)] )
@@ -35,7 +35,7 @@ def main():
     print "newton brukte %.4f" % (stop - start)
 
     start = time.clock()
-    nwst_points, nwst_residuals = newtsteep(x0, H, b, c, 0.01, 0.000000001)
+    nwst_points, nwst_residuals = steepestnewton(x0, H, b, c, 0.01, 0.000000001)
     stop = time.clock()
     print "newton + steepest brukte %.4f" % (stop - start)
 

diff --git a/newton.py b/newton.py
index 0be139a65ef87e73257866bfc7c78262f6714f16..c63922946ae76ef6cf461d0c2e7052a6cb4b2775 100644 (file)
--- a/newton.py
+++ b/newton.py
@@ -6,9 +6,11 @@ import numpy as np
 import matplotlib.pyplot as plt
 from numpy.linalg import norm, inv
 
+# Calculate the next iteration of the Newton method
 def newton_iter(x, grad_g, hess_g):
     return x - inv(hess_g).dot(grad_g)
 
+# Run the Newton method from point xk
 def newton(xk, H, b, c, tolerance, maxiter, x0_norm=0):
     if x0_norm == 0:
         x0_norm = norm(gradient_g(xk, H, b, c))
@@ -26,6 +28,7 @@ def newton(xk, H, b, c, tolerance, maxiter, x0_norm=0):
 
     return (points, relative_residuals)
 
+# Test the Newton method
 def main():
     H = np.array([ ( 1.8580682, -1.01197101), (-1.01197101, 0.63456543)] )
     b = np.array([-0.3119187,  0.14508516])

diff --git a/newtsteep.py b/newtsteep.py
deleted file mode 100644 (file)
index 977a02f..0000000
--- a/newtsteep.py
+++ /dev/null
@@ -1,25 +0,0 @@
-import numpy as np
-import generate as g
-import steepest as sd
-import newton
-import surf 
-import matplotlib.pyplot as plt
-from function import gradient_g
-from numpy.linalg import norm
-
-def newtsteep(x0, H, b, c, steep_tol, newt_tol):
-    x = []
-    r = []
-
-    # Finne optimal alpha og kjore en iterasjon med steepest.
-    x0_norm = norm(gradient_g(x0, H, b, c))
-    points, residuals = sd.steepest(x0, H, b, c, steep_tol, 100, 0, x0_norm)
-    x.extend(points)
-    r.extend(residuals)
-
-    # Angir tolerance og kjorer Newtons til vi finner nullpunkt.
-    points, residuals = newton.newton(x[len(x)-1], H, b, c, newt_tol, 100, x0_norm)
-    x.extend(points[1:])
-    r.extend(residuals[1:])
-
-    return (x, r)

diff --git a/steepest.py b/steepest.py
index 49147c48615ef9fa1d02e1b5913a0a21f7448640..976388c886c5db64e78626b3b5d65832743af259 100644 (file)
--- a/steepest.py
+++ b/steepest.py
@@ -7,9 +7,11 @@ import matplotlib.pyplot as plt
 import time
 from numpy.linalg import norm
 
+# Calculate next iteration of steepest method
 def steepest_iter(x, alpha, grad_f):
-    return x - alpha * grad_f(x)
+    return x - alpha * grad_f
 
+# Return optimal step size from point x in direction of -grad_f
 def optimal_step(x, H, b, c, grad_f):
     pk = -grad_f
     Cpk = _C(pk, c)
@@ -17,19 +19,24 @@ def optimal_step(x, H, b, c, grad_f):
 
     coeff = [0] * 4
 
+    # Setup the coefficients, where coeff[0] * a^3 etc.
     coeff[0] = (1.0/3.0) * pk.T.dot(Cpk.dot(pk))
     coeff[1] = x.T.dot(Cpk.dot(pk))
     coeff[2] = pk.T.dot(H.dot(pk)) + pk.T.dot(Cxk.dot(pk))
     coeff[3] = (1.0/3.0) * pk.T.dot(Cxk.dot(x)) + pk.T.dot(H.dot(x)) - b.T.dot(pk)
 
+    # Find roots of polynomial
     roots = np.roots(coeff)
 
+    # Only return the real roots
     for a in roots:
         if abs(a.imag) < 0.0000001:
             return a.real
 
     return 0.1
 
+# Run the steepest method, starting in x. Use optimal step size if step is set
+# to 0
 def steepest(x, H, b, c, tolerance, maxiter, step, x0_norm=0):
     xk = np.array(x)
     if x0_norm == 0:
@@ -51,29 +58,16 @@ def steepest(x, H, b, c, tolerance, maxiter, step, x0_norm=0):
         relative_residuals.append(relative_residual)
         points.append(xk)
     
-        xk = steepest_iter(xk, alpha, lambda x: gradient_g(x, H, b, c))
+        xk = steepest_iter(xk, alpha, gradient_g(xk))
     
-        z = function_g(xk, H, b, c)
-
     return (points, relative_residuals)
     
+# Test the functions above
 def main():
     H = generate.spdmatrix(2, 3)
     b = generate.vector(2, -1, 1)
     c = generate.vector(2, 0, 3)
 
-    #H = np.array([ ( 2.73146076, 0.75347534), ( 0.75347534, 2.56104263) ])
-    #b = np.array([ 0.74290645, 0.64954508])
-    #c = np.array([ 0.27501961, 2.73262633])
-
-    H = np.array([ ( 1.8580682, -1.01197101), (-1.01197101, 0.63456543)] )
-    b = np.array([-0.3119187,  0.14508516])
-    c = np.array([ 0.47736331, 1.48834978])
-
-    #print H
-    #print b
-    #print c
-
     x0 = np.array([1.5, 1.1])
     
     function_values = []

diff --git a/surf.py b/surf.py
index 1c1f2f4a31277841ebf42933d771f6e2daee693c..f0c5b35071f40ded724675df24e7bd98e5f149ef 100644 (file)
--- a/surf.py
+++ b/surf.py
@@ -1,8 +1,8 @@
-from scipy import NaN
 import numpy as np
 import function
 import generate
 
+# Calculate the function over a XY-plane
 def surf(H, b, c, xmin, xmax, ymin, ymax, zmin, zmax):
     X = np.linspace(xmin, xmax)
     Y = np.linspace(ymin, ymax)
@@ -12,7 +12,6 @@ def surf(H, b, c, xmin, xmax, ymin, ymax, zmin, zmax):
     Z = np.zeros( (dim,dim) )
     for i in xrange(dim):
         for j in xrange(dim):
-            Z[i][j] = function.function_g(np.array([X[i][j],Y[i][j]]), H, b, c)
-#            if Z[i][j] > zmax:
-#                Z[i][j] = NaN
+            Z[i][j] = function.function_g(np.array([X[i][j], Y[i][j]]), H, b, c)
+
     return X, Y, Z