From: Bjørn Rustad Date: Fri, 30 Sep 2011 19:03:35 +0000 (+0200) Subject: Prettify X-Git-Url: http://git.rustad.me/?a=commitdiff_plain;h=7da02d67d11560d7081e346fa3980f1fd080b99a;p=nummat1 Prettify --- diff --git a/function.py b/function.py index 5c9cd9f..a7a26a3 100644 --- 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 4af4c0b..7bcef1e 100644 --- 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 df0e882..8d2dd7a 100644 --- 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 0be139a..c639229 100644 --- 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 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 49147c4..976388c 100644 --- 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 1c1f2f4..f0c5b35 100644 --- 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