Chunk / Numerical

Download as
Browse Code »

Commit e5dbe4161943567f6179545ff3eeb31c50d3e33c

Authored by Chunk
1 parent 8a43d8ec
Exists in master

staged.

Showing 10 changed files with 310 additions and 48 deletions   Show diff stats
.idea/Numerical.iml
  Diff comments   View file @e5dbe41
... ... @@ -2,7 +2,7 @@
2 2 <module type="PYTHON_MODULE" version="4">
3 3 <component name="NewModuleRootManager">
4 4 <content url="file://$MODULE_DIR$" />
5   - <orderEntry type="jdk" jdkName="Python 2.7.8 virtualenv at ~/.virtualenvs/env0" jdkType="Python SDK" />
  5 + <orderEntry type="jdk" jdkName="Python 2.7.6 virtualenv at ~/.virtualenvs/env0" jdkType="Python SDK" />
6 6 <orderEntry type="sourceFolder" forTests="false" />
7 7 </component>
8 8 </module>
9 9 \ No newline at end of file
... ...
chap1/1_3_infinite.py 0 → 100644
  Diff comments   View file @e5dbe41
... ... @@ -0,0 +1,71 @@
  1 +__author__ = 'chunk'
  2 +
  3 +import struct
  4 +import numpy as np
  5 +
  6 +from common import *
  7 +
  8 +def float2bits(f, fmt='bin'):
  9 + if fmt == 'hex':
  10 + return hex(struct.unpack('!l', struct.pack('!f', f))[0])
  11 + return bin(struct.unpack('!l', struct.pack('!f', f))[0])
  12 +
  13 +
  14 +def double2bits(d, fmt='bin'):
  15 + if fmt == 'hex':
  16 + return hex(struct.unpack('!q', struct.pack('!d', d))[0])
  17 + return bin(struct.unpack('!q', struct.pack('!d', d))[0])
  18 +
  19 +
  20 +def float2bin(num):
  21 + # http://stackoverflow.com/questions/16444726/binary-representation-of-float-in-python-bits-not-hex
  22 + bits = [bin(ord(c)).replace('0b', '').rjust(8, '0') for c in struct.pack('!f', num)]
  23 + print bits
  24 + return ''.join(bits)
  25 +
  26 +
  27 +def double2bin(num):
  28 + bits = [bin(ord(c)).replace('0b', '').rjust(8, '0') for c in struct.pack('!d', num)]
  29 + print bits
  30 + return ''.join(bits)
  31 +
  32 +
  33 +def infinite_float():
  34 + s = np.float32(0)
  35 + tmp = np.float32(-1)
  36 + i = 0
  37 + while s != tmp:
  38 + i += 1
  39 + tmp = s
  40 + s += np.float32(1.0 / i)
  41 +
  42 + print i, s, tmp, float2bin(s)
  43 +
  44 +
  45 +def infinite_double(n):
  46 + s = np.float64(0)
  47 + for i in range(1, n + 1):
  48 + s += np.float64(1.0 / i)
  49 + return s
  50 +
  51 +
  52 +def testn0():
  53 + for i in range(100):
  54 + print i, np.log(i)
  55 + if np.log(i) >= 8:
  56 + print i
  57 +
  58 +
  59 +if __name__ == '__main__':
  60 + # testn0()
  61 + # infinite_float()
  62 +
  63 + # i = 1
  64 + # while infinite_double(i) < 16:
  65 + # i += 1
  66 + # print i, infinite_double(i)
  67 + timer = Timer()
  68 + timer.mark()
  69 + infinite_double(10000000) # 1125899906842624
  70 + timer.report()
  71 + pass
... ...
chap1/1_3_infiniteseries.py
View file @8a43d8e
... ... @@ -1,47 +0,0 @@
1   -__author__ = 'chunk'
2   -
3   -import struct
4   -import numpy as np
5   -
6   -def float2bits(f, fmt='bin'):
7   - if fmt == 'hex':
8   - return hex(struct.unpack('!l', struct.pack('!f', f))[0])
9   - return bin(struct.unpack('!l', struct.pack('!f', f))[0])
10   -
11   -
12   -def double2bits(d, fmt='bin'):
13   - if fmt == 'hex':
14   - return hex(struct.unpack('!q', struct.pack('!d', d))[0])
15   - return bin(struct.unpack('!q', struct.pack('!d', d))[0])
16   -
17   -
18   -def float2bin(num):
19   - # http://stackoverflow.com/questions/16444726/binary-representation-of-float-in-python-bits-not-hex
20   - bits = [bin(ord(c)).replace('0b', '').rjust(8, '0') for c in struct.pack('!f', num)]
21   - print bits
22   - return ''.join(bits)
23   -
24   -
25   -def double2bin(num):
26   - # http://stackoverflow.com/questions/16444726/binary-representation-of-float-in-python-bits-not-hex
27   - bits = [bin(ord(c)).replace('0b', '').rjust(8, '0') for c in struct.pack('!d', num)]
28   - print bits
29   - return ''.join(bits)
30   -
31   -
32   -def infiniteseries():
33   - s = np.float32(0)
34   - tmp = np.float32(-1)
35   - i = 0
36   - while s != tmp:
37   - i += 1
38   - tmp = s
39   - s += np.float32(1.0 / i)
40   -
41   - print i, s, tmp
42   -
43   -
44   -if __name__ == '__main__':
45   - print float2bin(1.0)
46   - infiniteseries()
47   - pass
chap2/__init__.py 0 → 100644
  Diff comments   View file @e5dbe41
... ... @@ -0,0 +1 @@
  1 +__author__ = 'chunk'
... ...
chap2/newton.py 0 → 100644
  Diff comments   View file @e5dbe41
... ... @@ -0,0 +1,67 @@
  1 +__author__ = 'chunk'
  2 +
  3 +import numpy as np
  4 +import math
  5 +
  6 +
  7 +def deriv(fn):
  8 + dx = 0.00001
  9 + return lambda x: (fn(x + dx) - fn(x)) / dx
  10 +
  11 +
  12 +def dampnm(fun, x0, lmbda, epsilon, gfun=None):
  13 + assert lmbda < 1 and lmbda > 0
  14 +
  15 + if gfun == None:
  16 + gfun = deriv(fun)
  17 + x_new = x0
  18 + x_old = None
  19 + # while np.absolute(fun(x_new)) > epsilon:
  20 + while x_old == None or np.absolute(fun(x_new) - fun(x_old)) > epsilon:
  21 + x_old = x_new
  22 + s = fun(x_old) / gfun(x_old)
  23 + x_new = x_old - s
  24 + while np.absolute(fun(x_new)) >= np.absolute(fun(x_old)):
  25 + x_new = x_old - lmbda * s
  26 + lmbda /= 2.0
  27 + print lmbda, x_new
  28 + return x_new
  29 +
  30 +
  31 +def test():
  32 + fn = lambda x: math.pow(x, 4) - x - 2
  33 + gfn = lambda x: 4 * math.pow(x, 3) - 1
  34 + print dampnm(fn, 1.5, 0.99, 0.00001, gfn)
  35 +
  36 + fn = lambda x: math.pow(x, 3) - x - 1
  37 + gfn = lambda x: 3 * math.pow(x, 2) - 1
  38 + print dampnm(fn, 0.6, 0.99, 0.00001, gfn)
  39 +
  40 + fn = lambda x: -math.pow(x, 3) + 5 * x
  41 + gfn = lambda x: -3 * math.pow(x, 2) + 5
  42 + print dampnm(fn, 1.2, 0.99, 0.00001, gfn)
  43 +
  44 + # 0.99 1.375
  45 + # 0.99 1.35377701579
  46 + # 0.99 1.35321036029
  47 + # 0.99 1.3532099642
  48 + # 1.3532099642
  49 + # 0.01546875 1.13521875
  50 + # 0.01546875 1.36976137194
  51 + # 0.01546875 1.32649943078
  52 + # 0.01546875 1.32472090757
  53 + # 0.01546875 1.32471795725
  54 + # 0.01546875 1.32471795724
  55 + # 1.32471795724
  56 + # 0.2475 -1.90976470588
  57 + # 0.2475 -2.34458457506
  58 + # 0.2475 -2.24316474291
  59 + # 0.2475 -2.23610151441
  60 + # 0.2475 -2.23606797825
  61 + # 0.2475 -2.2360679775
  62 + # -2.2360679775
  63 +
  64 +if __name__ == '__main__':
  65 + test()
  66 +
  67 + pass
... ...
chap3/__init__.py 0 → 100644
  Diff comments   View file @e5dbe41
... ... @@ -0,0 +1 @@
  1 +__author__ = 'chunk'
... ...
chap3/cholesky.py 0 → 100644
  Diff comments   View file @e5dbe41
... ... @@ -0,0 +1,92 @@
  1 +__author__ = 'chunk'
  2 +
  3 +import numpy as np
  4 +import math
  5 +
  6 +
  7 +def calc_upper(U, b):
  8 + n = len(U)
  9 + assert len(b) == n
  10 + x = [0.0] * n
  11 + for i in reversed(range(n)):
  12 + if U[i][i] == 0:
  13 + raise Exception('u_ii == 0')
  14 + x[i] = b[i]
  15 + for j in reversed(range(i + 1, n)):
  16 + x[i] -= U[i][j] * x[j]
  17 + x[i] /= U[i][i]
  18 + return x
  19 +
  20 +
  21 +def calc_lower(L, b):
  22 + n = len(L)
  23 + assert len(b) == n
  24 + x = [0.0] * n
  25 + for i in range(n):
  26 + if L[i][i] == 0:
  27 + raise Exception('l_ii == 0')
  28 + x[i] = b[i]
  29 + for j in range(i):
  30 + x[i] -= L[i][j] * x[j]
  31 + x[i] /= L[i][i]
  32 + return x
  33 +
  34 +
  35 +def gen_hilbert(n):
  36 + return np.array([1.0 / (i + j + 1) for i in range(n) for j in range(n)]).reshape((n, -1))
  37 +
  38 +
  39 +def gen_b(n, x=None):
  40 + if x == None:
  41 + x = [1] * n
  42 + return np.dot(gen_hilbert(n), x)
  43 +
  44 +
  45 +def cholesky(A):
  46 + n = len(A)
  47 + for j in range(n):
  48 + for k in range(j):
  49 + A[j][j] -= (A[j][k] * A[j][k])
  50 + A[j][j] = math.sqrt(A[j][j])
  51 +
  52 + for i in range(j + 1, n):
  53 + for k in range(j):
  54 + A[i][j] -= (A[i][k] * A[j][k])
  55 + A[i][j] /= A[j][j]
  56 +
  57 + return np.array(A)
  58 +
  59 +
  60 +def calculate(n):
  61 + H = gen_hilbert(n)
  62 + b = gen_b(n)
  63 + L = cholesky(H)
  64 + y = calc_lower(L, b)
  65 + x = calc_upper(L.T, y)
  66 + # validate:
  67 + # print np.dot(H, x) - b
  68 + return x, b - np.dot(H, x), [x[i] - 1.0 for i in range(n)]
  69 +
  70 +
  71 +def test():
  72 + print gen_hilbert(3)
  73 + print gen_b(5)
  74 +
  75 + U = [[1, -0.7, 0], [0, 1, -60], [0, 0, 155]]
  76 + b = [0.7, -61, 155]
  77 + print calc_upper(U, b)
  78 +
  79 + L = [[1, 0, 0], [-60, 1, 0], [0, -0.7, 1]]
  80 + b = [1, -61, 0.7]
  81 + print calc_lower(L, b)
  82 +
  83 + A = [[5, -1, -1], [-1, 3, -1], [-1, -1, 5]]
  84 + print cholesky(A)
  85 +
  86 + x, r, delta = calculate(10)
  87 + print "%s\n%s\n%s\n" % (x, r, delta)
  88 +
  89 +
  90 +if __name__ == '__main__':
  91 + test()
  92 + pass
... ...
chap4/__init__.py 0 → 100644
  Diff comments   View file @e5dbe41
... ... @@ -0,0 +1 @@
  1 +__author__ = 'chunk'
... ...
common.py 0 → 100644
  Diff comments   View file @e5dbe41
... ... @@ -0,0 +1,76 @@
  1 +"""
  2 +Common utils.
  3 +
  4 +@author: chunk
  5 +chunkplus@gmail.com
  6 +2014 Dec
  7 +"""
  8 +__author__ = 'hadoop'
  9 +
  10 +import os, sys
  11 +import time
  12 +import StringIO
  13 +import ConfigParser
  14 +
  15 +import numpy as np
  16 +
  17 +package_dir_imager = os.path.dirname(os.path.abspath(__file__))
  18 +
  19 +
  20 +class Timer():
  21 + def __init__(self):
  22 + self.__newtime = time.time()
  23 + self.__oldtime = self.__newtime
  24 +
  25 + def mark(self):
  26 + self.__oldtime = self.__newtime
  27 + self.__newtime = time.time()
  28 + return self.__newtime - self.__oldtime
  29 +
  30 + def report(self):
  31 + print "%-24s%fs" % ("time elapsed:", self.mark())
  32 +
  33 +
  34 +def ttimer():
  35 + newtime = time.time()
  36 + while True:
  37 + oldtime = newtime
  38 + newtime = time.time()
  39 + yield newtime - oldtime
  40 +
  41 +
  42 +def ctimer():
  43 + newtime = time.clock()
  44 + while True:
  45 + oldtime = newtime
  46 + newtime = time.clock()
  47 + yield newtime - oldtime
  48 +
  49 +
  50 +def bytes2bits(arry_bytes):
  51 + """
  52 + :param arry_bytes: 1-D np.unit8 array
  53 + :return: 1-D 0/1 array
  54 + """
  55 + hid_data_bits = [map(int, '{0:08b}'.format(byte)) for byte in arry_bytes]
  56 + return np.array(hid_data_bits).ravel()
  57 +
  58 +
  59 +def bits2bytes(arry_bits):
  60 + """
  61 + :param arry_bits: 1-D 0/1 array
  62 + :return: 1-D np.unit8 array
  63 + """
  64 + str_bits = ''.join(map(str, arry_bits))
  65 + arry_bytes = [int(str_bits[i:i + 8], 2) for i in range(0, len(str_bits), 8)]
  66 + return np.array(arry_bytes, dtype=np.uint8).ravel()
  67 +
  68 +
  69 +def test_grammer():
  70 + a = 'fsaf'
  71 + b = ['dasf', 'dff']
  72 + c = 'dgfsfdg'
  73 + # print a + b
  74 + print [a] + b # ['fsaf', 'dasf', 'dff']
  75 + print [a] + [b] # ['fsaf', ['dasf', 'dff']]
  76 + print [a] + [c] # ['fsaf', 'dgfsfdg']
... ...
common.pyc 0 → 100644
View file @e5dbe41
No preview for this file type