From a623ccf342d2083ff0aae09c8797cd1e7de02a61 Mon Sep 17 00:00:00 2001
From: =?utf8?q?Bj=C3=B8rn=20Rustad?= <rustadbjornen@gmail.com>
Date: Tue, 27 Sep 2011 18:05:43 +0200
Subject: [PATCH] Initial commit of working steepest descent method

---
 function.py | 32 ++++++++++++++++++++++++++++++++
 generate.py | 29 +++++++++++++++++++++++++++++
 steepest.py | 46 ++++++++++++++++++++++++++++++++++++++++++++++
 surf.py     | 18 ++++++++++++++++++
 4 files changed, 125 insertions(+)
 create mode 100644 function.py
 create mode 100644 generate.py
 create mode 100644 steepest.py
 create mode 100644 surf.py

diff --git a/function.py b/function.py
new file mode 100644
index 0000000..9d352e0
--- /dev/null
+++ b/function.py
@@ -0,0 +1,32 @@
+import numpy as np
+
+def function_g(x, H, b, c):
+    return -b.T.dot(x) + 0.5 * x.T.dot(H.dot(x)) + (1/12) * x.T.dot(_C(x,c).dot(x))
+
+def _C(x, c):
+    A = np.zeros( (len(c), len(c)) )
+    for i in xrange(0, len(c)):
+        A[i][i] = c[i]*(x[i])**2
+    return A
+
+def gradient_g(x, H, b, c):
+    return -b.T + 0.5 * H.dot(x) + 0.5 * H.T.dot(x) + (1/3) * x.T.dot(_C(x,c))
+
+def hessian_g(x, H, b, c):
+    return 0.5 * H.T + 0.5 * H + _C(x,c)
+
+def main():
+    A = np.array([ (1, 2), (2, 3) ])
+    x = np.array([5, 4])
+    b = np.array([1, 2])
+    c = np.array([7, 8])
+
+    print "Function:"
+    print function_g(x, A, b, c)
+    print "Gradient:"
+    print gradient_g(x, A, b, c)
+    print "Hessian:"
+    print hessian_g(x, A, b, c)
+
+if __name__ == "__main__":
+    main()
diff --git a/generate.py b/generate.py
new file mode 100644
index 0000000..7013e28
--- /dev/null
+++ b/generate.py
@@ -0,0 +1,29 @@
+import numpy as np
+import math
+import random
+
+def spdmatrix(n, limit):
+    A = np.zeros( (n,n) )
+    for y in xrange(0, n):
+        A[y][y] = math.sqrt(random.uniform(0, limit))
+        for x in xrange(y+1, n):
+            A[y][x] = random.uniform(-limit, limit)
+    return A.T.dot(A)
+
+def vector(n, lower, upper):
+    X = np.zeros(n)
+    for x in xrange(0, n):
+        X[x] = random.uniform(lower, upper)
+    return X
+
+def main():
+    A = spdmatrix(5, 10)
+    print "SPD-matrise:"
+    print A
+    print "Egenverdier:"
+    print np.linalg.eig(A)[0]
+    print "Tilfeldig vektor:"
+    print vector(10, -10, 10)
+
+if __name__ == "__main__":
+    main()
diff --git a/steepest.py b/steepest.py
new file mode 100644
index 0000000..04e55ce
--- /dev/null
+++ b/steepest.py
@@ -0,0 +1,46 @@
+from enthought.mayavi import mlab
+from function import function_g, gradient_g
+import generate
+import surf
+import numpy as np
+import matplotlib.pyplot as plt
+from numpy.linalg import norm
+
+def steepest_iter(x, alpha, grad_f):
+    return x - alpha * grad_f(x)
+
+H = generate.spdmatrix(2, 3)
+b = generate.vector(2, -1, 1)
+c = generate.vector(2, -3, 3)
+
+alpha = 0.1
+xiter = np.array([0, 0])
+x0_norm = norm(gradient_g(xiter, H, b, c))
+
+tolerance = 0.01
+
+relative_residual = norm(gradient_g(xiter, H, b, c)) / x0_norm
+relative_residuals = [relative_residual]
+
+while relative_residual > tolerance:
+    relative_residual = norm(gradient_g(xiter, H, b, c)) / x0_norm
+    relative_residuals.append(relative_residual)
+
+    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)
+    print xiter
+
+X, Y, Z = surf.surf(H, b, c, xiter[0]-2, xiter[0]+2, xiter[1]-2, xiter[1]+2, -10, 3)
+
+plt.figure(1)
+ax = plt.subplot(111)
+plt.plot(relative_residuals)
+ax.set_yscale('log')
+plt.show()
+
+#print relative_residuals
+
+s = mlab.mesh(X, Y, Z)
+axes = mlab.axes()
+mlab.show()
diff --git a/surf.py b/surf.py
new file mode 100644
index 0000000..27ccde9
--- /dev/null
+++ b/surf.py
@@ -0,0 +1,18 @@
+from scipy import NaN
+import numpy as np
+import function
+import generate
+
+def surf(H, b, c, xmin, xmax, ymin, ymax, zmin, zmax):
+    X = np.linspace(xmin, xmax)
+    Y = np.linspace(ymin, ymax)
+    X, Y = np.meshgrid(X, Y)
+
+    dim = len(X)
+    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
+    return X, Y, Z
-- 
2.47.3