From 7d0297f3f8739a34b4b89b3ecc5deb3862e92d1a Mon Sep 17 00:00:00 2001
From: Chunk <chunkplus@gmail.com>
Date: Sun, 21 Jun 2015 17:03:54 +0800
Subject: [PATCH] 6 finished.

---
 chap3/__init__.pyc | Bin 0 -> 180 bytes
 chap3/cholesky.py  | 49 +++++++++++++++++++++++++++++++++++++++++++++----
 chap3/cholesky.pyc | Bin 0 -> 4621 bytes
 chap4/iterative.py | 87 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 chap5/__init__.py  |  1 +
 chap5/power.py     | 62 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 chap6/__init__.py  |  1 +
 chap6/fitting.py   | 89 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 8 files changed, 285 insertions(+), 4 deletions(-)
 create mode 100644 chap3/__init__.pyc
 create mode 100644 chap3/cholesky.pyc
 create mode 100644 chap4/iterative.py
 create mode 100644 chap5/__init__.py
 create mode 100644 chap5/power.py
 create mode 100644 chap6/__init__.py
 create mode 100644 chap6/fitting.py

diff --git a/chap3/__init__.pyc b/chap3/__init__.pyc
new file mode 100644
index 0000000..ab57761
Binary files /dev/null and b/chap3/__init__.pyc differ
diff --git a/chap3/cholesky.py b/chap3/cholesky.py
index bcebac4..356b27c 100644
--- a/chap3/cholesky.py
+++ b/chap3/cholesky.py
@@ -42,6 +42,14 @@ def gen_b(n, x=None):
     return np.dot(gen_hilbert(n), x)
 
 
+def norm(l):
+    s = 0
+    for item in l:
+        if np.absolute(item) > s:
+            s = np.absolute(item)
+    return s
+
+
 def cholesky(A):
     n = len(A)
     for j in range(n):
@@ -57,7 +65,7 @@ def cholesky(A):
     return np.array(A)
 
 
-def calculate(n):
+def calculate0(n):
     H = gen_hilbert(n)
     b = gen_b(n)
     L = cholesky(H)
@@ -68,6 +76,19 @@ def calculate(n):
     return x, b - np.dot(H, x), [x[i] - 1.0 for i in range(n)]
 
 
+def calculate1(n):
+    H = gen_hilbert(n)
+    b = gen_b(n)
+    b_t = np.copy(b)
+    b_t[-1] += math.pow(10, -7)
+    L = cholesky(H)
+    y = calc_lower(L, b_t)
+    x = calc_upper(L.T, y)
+    # validate:
+    # print np.dot(H, x) - b
+    return x, b - np.dot(H, x), [x[i] - 1.0 for i in range(n)]
+
+
 def test():
     print gen_hilbert(3)
     print gen_b(5)
@@ -83,10 +104,30 @@ def test():
     A = [[5, -1, -1], [-1, 3, -1], [-1, -1, 5]]
     print cholesky(A)
 
-    x, r, delta = calculate(10)
-    print "%s\n%s\n%s\n" % (x, r, delta)
+
+def test1():
+    print "[origin(n=10)]"
+    x, r, delta = calculate0(10)
+    print "x:%s\nr:%s\ndelta:%s\n" % (x, r, delta)
+    print "norm(r):%s\nnorm(delta):%s\n" % (norm(r), norm(delta))
+
+    print "[disturbed(n=10)]"
+    x, r, delta = calculate1(10)
+    print "x:%s\nr:%s\ndelta:%s\n" % (x, r, delta)
+    print "norm(r):%s\nnorm(delta):%s\n" % (norm(r), norm(delta))
+
+    print "[n=8]"
+    x, r, delta = calculate0(8)
+    print "x:%s\nr:%s\ndelta:%s\n" % (x, r, delta)
+    print "norm(r):%s\nnorm(delta):%s\n" % (norm(r), norm(delta))
+
+    print "[n=12]"
+    x, r, delta = calculate0(12)
+    print "x:%s\nr:%s\ndelta:%s\n" % (x, r, delta)
+    print "norm(r):%s\nnorm(delta):%s\n" % (norm(r), norm(delta))
 
 
 if __name__ == '__main__':
-    test()
+    # test()
+    test1()
     pass
diff --git a/chap3/cholesky.pyc b/chap3/cholesky.pyc
new file mode 100644
index 0000000..c3e3de2
Binary files /dev/null and b/chap3/cholesky.pyc differ
diff --git a/chap4/iterative.py b/chap4/iterative.py
new file mode 100644
index 0000000..4c16853
--- /dev/null
+++ b/chap4/iterative.py
@@ -0,0 +1,87 @@
+__author__ = 'chunk'
+
+import numpy as np
+import math
+
+
+def gen_hilbert(n):
+    return np.array([1.0 / (i + j + 1) for i in range(n) for j in range(n)]).reshape((n, -1))
+
+
+def gen_b(n):
+    return np.array([1.0 / i for i in range(1, n + 1)])
+
+
+def norm(l):
+    s = 0
+    for item in l:
+        if np.absolute(item) > s:
+            s = np.absolute(item)
+    return s
+
+
+def jacobi(A, b, epsilon=0.0001, x=None):
+    n = len(A)
+    assert n == len(b)
+
+    if x is None:
+        x = np.zeros(n, dtype=np.float32)
+    x_old = None
+
+    D = np.diag(A)
+    R = A - np.diagflat(D)
+
+    # for i in range(N):
+    iter = 0
+    while x_old is None or norm(x - x_old) >= epsilon:
+        iter += 1
+        x_old = np.copy(x)
+        x = (b - np.dot(R, x_old)) / D
+        print iter, x
+    return x
+
+
+def SOR(A, b, omega, epsilon=0.0001, x=None):
+    n = len(A)
+    assert n == len(b)
+
+    if x is None:
+        x = np.zeros(n, dtype=np.float32)
+    x_old = None
+    iter = 0
+    while x_old is None or norm(x - x_old) >= epsilon:
+        iter += 1
+        x_old = np.copy(x)
+        for i in range(n):
+            s = 0.0
+            for j in range(n):
+                if j != i:
+                    s += A[i][j] * x[j]
+            x[i] = (1 - omega) * x[i] + omega * (b[i] - s) / A[i][i]
+
+            # print iter, x
+    return iter, x
+
+
+def test():
+    # A = [[10, 3, 1], [2, -10, 3], [1, 3, 10]]
+    # b = [14, -5, 14]
+    # print SOR(A, b, 1.25)
+    # print jacobi(A, b)
+
+    n = 10
+    A = gen_hilbert(n)
+    b = gen_b(n)
+    # print jacobi(A, b)
+    print SOR(A, b, 1)
+
+    for w in np.linspace(0.5, 1.5, num=11):
+        iter, x = SOR(A, b, w)
+        # print "w:%f\titer:%d\tx:%s\t" % (w, iter, x)
+        # print "norm(delta):%s\n" % (norm(np.array([1.0] + [0.0] * (n - 1)) - x))
+        print "w:%f\titer:%d\tnorm(delta):%s\n" % (
+        w, iter, norm(np.array([1.0] + [0.0] * (n - 1)) - x))
+
+
+if __name__ == '__main__':
+    test()
diff --git a/chap5/__init__.py b/chap5/__init__.py
new file mode 100644
index 0000000..a1459cf
--- /dev/null
+++ b/chap5/__init__.py
@@ -0,0 +1 @@
+__author__ = 'chunk'
diff --git a/chap5/power.py b/chap5/power.py
new file mode 100644
index 0000000..f7bc8bf
--- /dev/null
+++ b/chap5/power.py
@@ -0,0 +1,62 @@
+__author__ = 'chunk'
+
+import numpy as np
+import math
+
+
+def norm(l):
+    s = 0
+    for item in l:
+        if np.absolute(item) > s:
+            s = np.absolute(item)
+    return s
+
+
+def norm_vec(l):
+    s = 0
+    loc = 0
+    i = -1
+    for item in l:
+        i += 1
+        if np.absolute(item) > s:
+            s = np.absolute(item)
+            loc = i
+    return l[loc], np.array(l, dtype=np.float32) / l[loc]
+
+
+def power_method(A, epsilon=0.00001, v=None):
+    n = len(A)
+    if v == None:
+        v = np.random.rand(n)
+    u = np.copy(v)
+    lmbda = None
+    lmbda_old = None
+
+    iter = 0
+    while lmbda_old is None or np.absolute(lmbda - lmbda_old) >= epsilon:
+        iter += 1
+        lmbda_old = lmbda
+        v = np.dot(A, u)
+        lmbda, u = norm_vec(v)
+        print iter, lmbda
+
+    return lmbda, u
+
+
+def test():
+    ll = [0, 2, 0.5, 1, -3, 0.2]
+    ll2 = [i * 2 for i in ll]
+    print norm_vec(ll)
+    print norm_vec(ll2)
+
+    A = [[3, 1], [1, 3]]
+    print power_method(A)
+
+    A = [[5, -4, 1], [-4, 6, -4], [1, -4, 7]]
+    B = [[25, -41, 10, -6], [-41, 68, -17, 10], [10, -17, 5, -3], [-6, 10, -3, 2]]
+    print power_method(A)
+    print power_method(B)
+
+
+if __name__ == '__main__':
+    test()
diff --git a/chap6/__init__.py b/chap6/__init__.py
new file mode 100644
index 0000000..a1459cf
--- /dev/null
+++ b/chap6/__init__.py
@@ -0,0 +1 @@
+__author__ = 'chunk'
diff --git a/chap6/fitting.py b/chap6/fitting.py
new file mode 100644
index 0000000..e467452
--- /dev/null
+++ b/chap6/fitting.py
@@ -0,0 +1,89 @@
+__author__ = 'chunk'
+
+import numpy as np
+import math
+import matplotlib.pyplot as plt
+import seaborn as sns
+
+from chap3.cholesky import calc_upper, calc_lower, cholesky
+
+plt.ticklabel_format(style='sci', axis='both')
+
+
+def fiiting(t, f, phi):
+    m, n = len(t), len(phi)
+    assert m == len(f)
+
+    A = np.array([phi[j](t[i]) for i in range(m) for j in range(n)]).reshape(m, -1)
+    G = np.dot(A.T, A)
+    b = np.dot(A.T, f)
+    L = cholesky(G)
+    y = calc_lower(L, b)
+    x = calc_upper(L.T, y)
+    return x
+
+
+def test0():
+    t0 = [1, 2, 3, 4, 5]
+    f0 = [4, 4.5, 6, 8, 8.5]
+    phi = [lambda x: 1, lambda x: x]
+    coef = fiiting(t0, f0, phi)
+    print coef
+
+    fit = lambda x: coef[0] + coef[1] * x * x
+    t1 = np.linspace(1, 8, num=1000).tolist()
+    f1 = [fit(i) for i in t1]
+
+    plt.scatter(t0, f0)
+    plt.plot(t1, f1)
+    plt.show()
+
+
+def test():
+    # print  np.linspace(1, 8, num=15).tolist()
+
+
+
+    t = [1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0]
+    f = [33.40, 79.50, 122.65, 159.05, 189.15, 214.15, 238.65, 252.2, 267.55, 280.50, 296.65,
+         301.65, 310.40,
+         318.15, 325.15]
+    phi = [lambda x: 1, lambda x: x, lambda x: x * x]
+    coef = fiiting(t, f, phi)
+    print coef
+    fit = lambda x: coef[0] + coef[1] * x + coef[2] * x * x
+    t1 = np.linspace(1, 8, num=1000).tolist()
+    f1 = [fit(i) for i in t1]
+
+
+    ff = [np.log(i) for i in f]
+    print ff
+    phi2 = [lambda x: 1, lambda x: x]
+    coef2 = fiiting(t, ff, phi2)
+    print coef2
+    fit2 = lambda x: np.exp(coef2[0]) * np.exp(coef2[1] *  x)
+    t2 = np.linspace(1, 8, num=1000).tolist()
+    f2 = [fit2(i) for i in t2]
+
+    tt = [1.0/i for i in t]
+    ff = [np.log(i) for i in f]
+    print ff
+    phi2 = [lambda x: 1, lambda x: x]
+    coef2 = fiiting(tt, ff, phi2)
+    print coef2
+    fit2 = lambda x: np.exp(coef2[0]) * np.exp(coef2[1] * 1.0 / x)
+    t3 = np.linspace(1, 8, num=1000).tolist()
+    f3 = [fit2(i) for i in t3]
+
+    plt.scatter(t, f)
+    plt.plot(t1, f1, label='polynomial')
+    plt.plot(t2, f2, label='exponential', linestyle='--')
+    plt.plot(t3, f3, label='exponential 1/t', linestyle='-')
+    plt.xlabel("t")
+    plt.ylabel("y")
+    plt.legend(loc=2)
+    plt.show()
+
+
+if __name__ == '__main__':
+    test()
--
libgit2 0.21.2