diff --git a/.idea/Numerical.iml b/.idea/Numerical.iml index d06e659..4f21c2d 100644 --- a/.idea/Numerical.iml +++ b/.idea/Numerical.iml @@ -2,7 +2,7 @@ - + \ No newline at end of file diff --git a/chap1/1_3_infinite.py b/chap1/1_3_infinite.py new file mode 100644 index 0000000..048ba4b --- /dev/null +++ b/chap1/1_3_infinite.py @@ -0,0 +1,71 @@ +__author__ = 'chunk' + +import struct +import numpy as np + +from common import * + +def float2bits(f, fmt='bin'): + if fmt == 'hex': + return hex(struct.unpack('!l', struct.pack('!f', f))[0]) + return bin(struct.unpack('!l', struct.pack('!f', f))[0]) + + +def double2bits(d, fmt='bin'): + if fmt == 'hex': + return hex(struct.unpack('!q', struct.pack('!d', d))[0]) + return bin(struct.unpack('!q', struct.pack('!d', d))[0]) + + +def float2bin(num): + # http://stackoverflow.com/questions/16444726/binary-representation-of-float-in-python-bits-not-hex + bits = [bin(ord(c)).replace('0b', '').rjust(8, '0') for c in struct.pack('!f', num)] + print bits + return ''.join(bits) + + +def double2bin(num): + bits = [bin(ord(c)).replace('0b', '').rjust(8, '0') for c in struct.pack('!d', num)] + print bits + return ''.join(bits) + + +def infinite_float(): + s = np.float32(0) + tmp = np.float32(-1) + i = 0 + while s != tmp: + i += 1 + tmp = s + s += np.float32(1.0 / i) + + print i, s, tmp, float2bin(s) + + +def infinite_double(n): + s = np.float64(0) + for i in range(1, n + 1): + s += np.float64(1.0 / i) + return s + + +def testn0(): + for i in range(100): + print i, np.log(i) + if np.log(i) >= 8: + print i + + +if __name__ == '__main__': + # testn0() + # infinite_float() + + # i = 1 + # while infinite_double(i) < 16: + # i += 1 + # print i, infinite_double(i) + timer = Timer() + timer.mark() + infinite_double(10000000) # 1125899906842624 + timer.report() + pass diff --git a/chap1/1_3_infiniteseries.py b/chap1/1_3_infiniteseries.py deleted file mode 100644 index 3b6d040..0000000 --- a/chap1/1_3_infiniteseries.py +++ /dev/null @@ -1,47 +0,0 @@ -__author__ = 'chunk' - -import struct -import numpy as np - -def float2bits(f, fmt='bin'): - if fmt == 'hex': - return hex(struct.unpack('!l', struct.pack('!f', f))[0]) - return bin(struct.unpack('!l', struct.pack('!f', f))[0]) - - -def double2bits(d, fmt='bin'): - if fmt == 'hex': - return hex(struct.unpack('!q', struct.pack('!d', d))[0]) - return bin(struct.unpack('!q', struct.pack('!d', d))[0]) - - -def float2bin(num): - # http://stackoverflow.com/questions/16444726/binary-representation-of-float-in-python-bits-not-hex - bits = [bin(ord(c)).replace('0b', '').rjust(8, '0') for c in struct.pack('!f', num)] - print bits - return ''.join(bits) - - -def double2bin(num): - # http://stackoverflow.com/questions/16444726/binary-representation-of-float-in-python-bits-not-hex - bits = [bin(ord(c)).replace('0b', '').rjust(8, '0') for c in struct.pack('!d', num)] - print bits - return ''.join(bits) - - -def infiniteseries(): - s = np.float32(0) - tmp = np.float32(-1) - i = 0 - while s != tmp: - i += 1 - tmp = s - s += np.float32(1.0 / i) - - print i, s, tmp - - -if __name__ == '__main__': - print float2bin(1.0) - infiniteseries() - pass diff --git a/chap2/__init__.py b/chap2/__init__.py new file mode 100644 index 0000000..a1459cf --- /dev/null +++ b/chap2/__init__.py @@ -0,0 +1 @@ +__author__ = 'chunk' diff --git a/chap2/newton.py b/chap2/newton.py new file mode 100644 index 0000000..3fc75df --- /dev/null +++ b/chap2/newton.py @@ -0,0 +1,67 @@ +__author__ = 'chunk' + +import numpy as np +import math + + +def deriv(fn): + dx = 0.00001 + return lambda x: (fn(x + dx) - fn(x)) / dx + + +def dampnm(fun, x0, lmbda, epsilon, gfun=None): + assert lmbda < 1 and lmbda > 0 + + if gfun == None: + gfun = deriv(fun) + x_new = x0 + x_old = None + # while np.absolute(fun(x_new)) > epsilon: + while x_old == None or np.absolute(fun(x_new) - fun(x_old)) > epsilon: + x_old = x_new + s = fun(x_old) / gfun(x_old) + x_new = x_old - s + while np.absolute(fun(x_new)) >= np.absolute(fun(x_old)): + x_new = x_old - lmbda * s + lmbda /= 2.0 + print lmbda, x_new + return x_new + + +def test(): + fn = lambda x: math.pow(x, 4) - x - 2 + gfn = lambda x: 4 * math.pow(x, 3) - 1 + print dampnm(fn, 1.5, 0.99, 0.00001, gfn) + + fn = lambda x: math.pow(x, 3) - x - 1 + gfn = lambda x: 3 * math.pow(x, 2) - 1 + print dampnm(fn, 0.6, 0.99, 0.00001, gfn) + + fn = lambda x: -math.pow(x, 3) + 5 * x + gfn = lambda x: -3 * math.pow(x, 2) + 5 + print dampnm(fn, 1.2, 0.99, 0.00001, gfn) + + # 0.99 1.375 + # 0.99 1.35377701579 + # 0.99 1.35321036029 + # 0.99 1.3532099642 + # 1.3532099642 + # 0.01546875 1.13521875 + # 0.01546875 1.36976137194 + # 0.01546875 1.32649943078 + # 0.01546875 1.32472090757 + # 0.01546875 1.32471795725 + # 0.01546875 1.32471795724 + # 1.32471795724 + # 0.2475 -1.90976470588 + # 0.2475 -2.34458457506 + # 0.2475 -2.24316474291 + # 0.2475 -2.23610151441 + # 0.2475 -2.23606797825 + # 0.2475 -2.2360679775 + # -2.2360679775 + +if __name__ == '__main__': + test() + + pass diff --git a/chap3/__init__.py b/chap3/__init__.py new file mode 100644 index 0000000..a1459cf --- /dev/null +++ b/chap3/__init__.py @@ -0,0 +1 @@ +__author__ = 'chunk' diff --git a/chap3/cholesky.py b/chap3/cholesky.py new file mode 100644 index 0000000..bcebac4 --- /dev/null +++ b/chap3/cholesky.py @@ -0,0 +1,92 @@ +__author__ = 'chunk' + +import numpy as np +import math + + +def calc_upper(U, b): + n = len(U) + assert len(b) == n + x = [0.0] * n + for i in reversed(range(n)): + if U[i][i] == 0: + raise Exception('u_ii == 0') + x[i] = b[i] + for j in reversed(range(i + 1, n)): + x[i] -= U[i][j] * x[j] + x[i] /= U[i][i] + return x + + +def calc_lower(L, b): + n = len(L) + assert len(b) == n + x = [0.0] * n + for i in range(n): + if L[i][i] == 0: + raise Exception('l_ii == 0') + x[i] = b[i] + for j in range(i): + x[i] -= L[i][j] * x[j] + x[i] /= L[i][i] + return x + + +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, x=None): + if x == None: + x = [1] * n + return np.dot(gen_hilbert(n), x) + + +def cholesky(A): + n = len(A) + for j in range(n): + for k in range(j): + A[j][j] -= (A[j][k] * A[j][k]) + A[j][j] = math.sqrt(A[j][j]) + + for i in range(j + 1, n): + for k in range(j): + A[i][j] -= (A[i][k] * A[j][k]) + A[i][j] /= A[j][j] + + return np.array(A) + + +def calculate(n): + H = gen_hilbert(n) + b = gen_b(n) + L = cholesky(H) + y = calc_lower(L, b) + 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) + + U = [[1, -0.7, 0], [0, 1, -60], [0, 0, 155]] + b = [0.7, -61, 155] + print calc_upper(U, b) + + L = [[1, 0, 0], [-60, 1, 0], [0, -0.7, 1]] + b = [1, -61, 0.7] + print calc_lower(L, b) + + 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) + + +if __name__ == '__main__': + test() + pass diff --git a/chap4/__init__.py b/chap4/__init__.py new file mode 100644 index 0000000..a1459cf --- /dev/null +++ b/chap4/__init__.py @@ -0,0 +1 @@ +__author__ = 'chunk' diff --git a/common.py b/common.py new file mode 100644 index 0000000..9a82dd0 --- /dev/null +++ b/common.py @@ -0,0 +1,76 @@ +""" +Common utils. + +@author: chunk +chunkplus@gmail.com +2014 Dec +""" +__author__ = 'hadoop' + +import os, sys +import time +import StringIO +import ConfigParser + +import numpy as np + +package_dir_imager = os.path.dirname(os.path.abspath(__file__)) + + +class Timer(): + def __init__(self): + self.__newtime = time.time() + self.__oldtime = self.__newtime + + def mark(self): + self.__oldtime = self.__newtime + self.__newtime = time.time() + return self.__newtime - self.__oldtime + + def report(self): + print "%-24s%fs" % ("time elapsed:", self.mark()) + + +def ttimer(): + newtime = time.time() + while True: + oldtime = newtime + newtime = time.time() + yield newtime - oldtime + + +def ctimer(): + newtime = time.clock() + while True: + oldtime = newtime + newtime = time.clock() + yield newtime - oldtime + + +def bytes2bits(arry_bytes): + """ + :param arry_bytes: 1-D np.unit8 array + :return: 1-D 0/1 array + """ + hid_data_bits = [map(int, '{0:08b}'.format(byte)) for byte in arry_bytes] + return np.array(hid_data_bits).ravel() + + +def bits2bytes(arry_bits): + """ + :param arry_bits: 1-D 0/1 array + :return: 1-D np.unit8 array + """ + str_bits = ''.join(map(str, arry_bits)) + arry_bytes = [int(str_bits[i:i + 8], 2) for i in range(0, len(str_bits), 8)] + return np.array(arry_bytes, dtype=np.uint8).ravel() + + +def test_grammer(): + a = 'fsaf' + b = ['dasf', 'dff'] + c = 'dgfsfdg' + # print a + b + print [a] + b # ['fsaf', 'dasf', 'dff'] + print [a] + [b] # ['fsaf', ['dasf', 'dff']] + print [a] + [c] # ['fsaf', 'dgfsfdg'] diff --git a/common.pyc b/common.pyc new file mode 100644 index 0000000..3f6cc80 Binary files /dev/null and b/common.pyc differ -- libgit2 0.21.2