#coding=utf-8
import numpy as np
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import time
from scipy.optimize import leastsq
from scipy import stats
import scipy.optimize as opt
import matplotlib.pyplot as plt
from scipy.stats import norm, poisson
from scipy.interpolate import BarycentricInterpolator
from scipy.interpolate import CubicSpline
import scipy as sp
import math

def residual(t, x, y):
    return y - (t[0] * x ** 2 + t[1] * x + t[2])

def residual2(t, x, y):
    print(t[0], t[1])
    return y - (t[0]*np.sin(t[1]*x) + t[2])

def f(x):
    y = np.ones_like(x)
    i = x > 0
    y[i] = np.power(x[i], x[i])
    i = x < 0
    y[i] = np.power(-x[i], -x[i])
    return y

a = np.arange(0, 60, 10).reshape((-1, 1)) + np.arange(6)
print(a)

[[ 0  1  2  3  4  5]
 [10 11 12 13 14 15]
 [20 21 22 23 24 25]
 [30 31 32 33 34 35]
 [40 41 42 43 44 45]
 [50 51 52 53 54 55]]

# 通过array函数传递list对象
L = [1, 2, 3, 4, 5, 6]
print('L = ', L)
a = np.array(L)
print('a = ', a)
print(type(a),type(L))

L =  [1, 2, 3, 4, 5, 6]
a =  [1 2 3 4 5 6]
 

# # 若传递的是多层嵌套的list，将创建多维数组
b = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])
print(b)

[[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9 10 11 12]]

# # # # # 数组大小可以通过其shape属性获得
print(a.shape)
print(b.shape)

(6,)
(3, 4)

# # # 也可以强制修改shape
b.shape = 4, 3
print(b)
# # 注：从(3,4)改为(4,3)并不是对数组进行转置，而只是改变每个轴的大小，数组元素在内存中的位置并没有改变

[[ 1  2  3]
 [ 4  5  6]
 [ 7  8  9]
 [10 11 12]]

# # # 当某个轴为-1时，将根据数组元素的个数自动计算此轴的长度
b.shape = 2, -1
print(b)
print(b.shape)

[[ 1  2  3  4  5  6]
 [ 7  8  9 10 11 12]]
(2, 6)

b.shape = 3, 4
print(b)
# # # # 使用reshape方法，可以创建改变了尺寸的新数组，原数组的shape保持不变
c = b.reshape((4, -1))
print('b = \n', b)
print('c = \n', c)

[[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9 10 11 12]]
b = 
 [[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9 10 11 12]]
c = 
 [[ 1  2  3]
 [ 4  5  6]
 [ 7  8  9]
 [10 11 12]]

# # # # 数组b和c共享内存，修改任意一个将影响另外一个
c[1][2] = 600
print('b = \n', b)
print('c = \n', c)

b = 
 [[  1   2   3   4]
 [  5 600   7   8]
 [  9  10  11  12]]
c = 
 [[  1   2   3]
 [  4   5 600]
 [  7   8   9]
 [ 10  11  12]]

# # # # 数组的元素类型可以通过dtype属性获得
print(a.dtype)
print(b.dtype)

int32
int32

 # # # # # 可以通过dtype参数在创建时指定元素类型
d = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]], dtype=np.float)
f = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]], dtype=np.complex)
print(d)
print(f)

[[  1.   2.   3.   4.]
 [  5.   6.   7.   8.]
 [  9.  10.  11.  12.]]
[[  1.+0.j   2.+0.j   3.+0.j   4.+0.j]
 [  5.+0.j   6.+0.j   7.+0.j   8.+0.j]
 [  9.+0.j  10.+0.j  11.+0.j  12.+0.j]]

# # # 如果更改元素类型，可以使用astype安全的转换
f = d.astype(np.int)
print(f)

[[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9 10 11 12]]

# 如果生成一定规则的数据，可以使用NumPy提供的专门函数
# arange函数类似于python的range函数：指定起始值、终止值和步长来创建数组
# 和Python的range类似，arange同样不包括终值；但arange可以生成浮点类型，而range只能是整数类型
a = np.arange(1,20,3)
print(a)
print(type(a),a.dtype)

[ 1  4  7 10 13 16 19]
 int32

# # # # linspace函数通过指定起始值、终止值和元素个数来创建数组，缺省包括终止值
b = np.linspace(1, 10, 10)
print('b = ', b)

b =  [  1.   2.   3.   4.   5.   6.   7.   8.   9.  10.]

# # 可以通过endpoint关键字指定是否包括终值
c = np.linspace(1, 10, 10, endpoint=False)
print('c = ', c)

c =  [ 1.   1.9  2.8  3.7  4.6  5.5  6.4  7.3  8.2  9.1]

# # # 和linspace类似，logspace可以创建等比数列
# # 下面函数创建起始值为10^1，终止值为10^2，有10个数的等比数列
d = np.logspace(1, 4, 4, endpoint=True,base=2)
print('d = ', d)

d =  [  2.   4.   8.  16.]

# # # # # 下面创建起始值为2^0，终止值为2^10(包括)，有10个数的等比数列
np.set_printoptions(suppress=True, linewidth=150)
f = np.logspace(0, 10 ,11, endpoint=True,base=2)
print('f = ', f)

f =  [    1.     2.     4.     8.    16.    32.    64.   128.   256.   512.  1024.]

# # # # 使用 frombuffer, fromstring, fromfile等函数可以从字节序列创建数组
s = 'abcdABCDEF'
g = np.fromstring(s, dtype=np.int8)
print(g)

[ 97  98  99 100  65  66  67  68  69  70]

#3.1常规办法：数组元素的存取方法和Python的标准方法相同
a = np.arange(10)
print(a)

[0 1 2 3 4 5 6 7 8 9]

# 获取某个元素
print(a[3])

3

# # # # 切片[3,6)，左闭右开
print(a[3:6])

[3 4 5]

# # 省略开始下标，表示从0开始
print(a[:5])

[0 1 2 3 4]

# # 从开始到最后
print(a[3:])

[3 4 5 6 7 8 9]

# 步长为2
print(a[1:9:2])

[1 3 5 7]

# # # # # 步长为-1，即翻转
print(a[::-1])

[9 8 7 6 5 4 3 2 1 0]

# 切片数据是原数组的一个视图，与原数组共享内容空间，可以直接修改元素值
a[1:4] = 10, 20, 30
print(a)

[ 0 10 20 30  4  5  6  7  8  9]

# 因此，在实践中，切实注意原始数据是否被破坏，如：
b = a[2:5]
print(b)
b[0] = 100
print(a)

[100  30   4]
[  0  10 100  30   4   5   6   7   8   9]

#根据整数数组存取：当使用整数序列对数组元素进行存取时，
# 将使用整数序列中的每个元素作为下标，整数序列可以是列表(list)或者数组(ndarray)。
# 使用整数序列作为下标获得的数组不和原始数组共享数据空间。
a = np.logspace(0, 9, 10, base=2)
print(a)
i = np.arange(0, 10, 2)
print(i)
# # 利用i取a中的元素
b = a[i]
print(b)
# # b的元素更改，a中元素不受影响
b[2] = 1.6
print(b)
print(a)

[   1.    2.    4.    8.   16.   32.   64.  128.  256.  512.]
[0 2 4 6 8]
[   1.    4.   16.   64.  256.]
[   1.     4.     1.6   64.   256. ]
[   1.    2.    4.    8.   16.   32.   64.  128.  256.  512.]

# 使用布尔数组i作为下标存取数组a中的元素：返回数组a中所有在数组b中对应下标为True的元素
# 生成10个满足[0,1)中均匀分布的随机数
np.set_printoptions(linewidth=200)
a = np.random.rand(10)
print(a)

[ 0.50849212  0.89850527  0.25803128  0.09660402  0.96200096  0.48726367  0.05403987  0.50161325  0.96409502  0.82349272]

# 大于0.5的元素索引
print(a > 0.5)

[ True  True False False  True False False  True  True  True]

# # 大于0.5的元素
print(a[a > 0.5])

[ 0.50849212  0.89850527  0.96200096  0.50161325  0.96409502  0.82349272]

# # # 将原数组中大于0.5的元素截取成0.5
a[a > 0.5] = 0.5
print(a)

[ 0.5         0.5         0.25803128  0.09660402  0.5         0.48726367  0.05403987  0.5         0.5         0.5       ]

# # # # b不受影响
print(b)

[   1.     4.     1.6   64.   256. ]

a = np.arange(0,60,10) # 行向量
print('a = ', a)

a =  [ 0 10 20 30 40 50]

b= a.reshape((-1, 1))# 转换成列向量
print('b = ', b)

b =  [[ 0]
 [10]
 [20]
 [30]
 [40]
 [50]]

c = np.arange(6)
print('c = ', c)
f = b + c
print(f)

c =  [0 1 2 3 4 5]
[[ 0  1  2  3  4  5]
 [10 11 12 13 14 15]
 [20 21 22 23 24 25]
 [30 31 32 33 34 35]
 [40 41 42 43 44 45]
 [50 51 52 53 54 55]]

# 合并上述代码
a = np.arange(0, 60, 10).reshape((-1, 1)) + np.arange(6)
print(a)

[[ 0  1  2  3  4  5]
 [10 11 12 13 14 15]
 [20 21 22 23 24 25]
 [30 31 32 33 34 35]
 [40 41 42 43 44 45]
 [50 51 52 53 54 55]]

# # # 二维数组的切片
print(a[[0, 1, 2], [2, 3, 4]])
print(a[4, [2, 3, 4]])
print(a[4:, [2, 3, 4]])

[ 2 13 24]
[42 43 44]
[[42 43 44]
 [52 53 54]]

i = np.array([True, False, True, False, False, True])
print(a[i])
print(a[i, 3])

[[ 0  1  2  3  4  5]
 [20 21 22 23 24 25]
 [50 51 52 53 54 55]]
[ 3 23 53]

for j in np.logspace(0, 7, 8):
        x = np.linspace(0, 10, j)
        start = time.clock()
        y = np.sin(x)
        t1 = time.clock() - start
        
x = x.tolist()
start = time.clock()
for i, t in enumerate(x):
    x[i] = math.sin(t)
t2 = time.clock() - start
print(j, ": ", t1, t2, t2/t1)

C:\Program Files\Anaconda3\lib\site-packages\ipykernel\__main__.py:2: DeprecationWarning: object of type  cannot be safely interpreted as an integer.
  from ipykernel import kernelapp as app


10000000.0 :  0.17834738633189676 4.184879334368658 23.464764022842356

# 4.2.1直接使用库函数
a = np.array((1, 2, 3, 4, 5, 5, 7, 3, 2, 2, 8, 8))
print('原始数组：', a)

原始数组： [1 2 3 4 5 5 7 3 2 2 8 8]

# # 使用库函数unique
b = np.unique(a)
print('去重后：', b)

去重后： [1 2 3 4 5 7 8]

# # 4.2.2 二维数组的去重，结果会是预期的么？
c = np.array(((1, 2), (3, 4), (5, 6), (1, 3), (3, 4), (7, 6)))
print('二维数组：\n', c)
print('去重后：', np.unique(c))
# # 4.2.3 方案1：转换为虚数
#4.2.3 方案2：利用set
print('去重方案2：\n', np.array(list(set([tuple(t) for t in c]))))

二维数组：
 [[1 2]
 [3 4]
 [5 6]
 [1 3]
 [3 4]
 [7 6]]
去重后： [1 2 3 4 5 6 7]
去重方案2：
 [[1 2]
 [7 6]
 [1 3]
 [3 4]
 [5 6]]

a = np.arange(1, 7).reshape((2, 3))
b = np.arange(11, 17).reshape((2, 3))
c = np.arange(21, 27).reshape((2, 3))
d = np.arange(31, 37).reshape((2, 3))
print('a = \n', a)
print('b = \n', b)
print('c = \n', c)
print('d = \n', d)

a = 
 [[1 2 3]
 [4 5 6]]
b = 
 [[11 12 13]
 [14 15 16]]
c = 
 [[21 22 23]
 [24 25 26]]
d = 
 [[31 32 33]
 [34 35 36]]

s = np.stack((a, b, c, d), axis=0)
print('axis = 0 ', s.shape, '\n', s)

axis = 0  (4, 2, 3) 
 [[[ 1  2  3]
  [ 4  5  6]]

 [[11 12 13]
  [14 15 16]]

 [[21 22 23]
  [24 25 26]]

 [[31 32 33]
  [34 35 36]]]

s = np.stack((a, b, c, d), axis=1)
print('axis = 1 ', s.shape, '\n', s)

axis = 1  (2, 4, 3) 
 [[[ 1  2  3]
  [11 12 13]
  [21 22 23]
  [31 32 33]]

 [[ 4  5  6]
  [14 15 16]
  [24 25 26]
  [34 35 36]]]

s = np.stack((a, b, c, d), axis=2)
print('axis = 2 ', s.shape, '\n', s)

axis = 2  (2, 3, 4) 
 [[[ 1 11 21 31]
  [ 2 12 22 32]
  [ 3 13 23 33]]

 [[ 4 14 24 34]
  [ 5 15 25 35]
  [ 6 16 26 36]]]

a = np.arange(1, 10).reshape(3,3)
print('a = ', a)
b = a + 10
print('b = ', b)
print(np.dot(a, b))
print(a * b)

a =  [[1 2 3]
 [4 5 6]
 [7 8 9]]
b =  [[11 12 13]
 [14 15 16]
 [17 18 19]]
[[ 90  96 102]
 [216 231 246]
 [342 366 390]]
[[ 11  24  39]
 [ 56  75  96]
 [119 144 171]]

a = np.arange(1, 10)
print(a)
b = np.arange(20, 25)
print(b)
print(np.concatenate((a, b)))

[1 2 3 4 5 6 7 8 9]
[20 21 22 23 24]
[ 1  2  3  4  5  6  7  8  9 20 21 22 23 24]

#基本配置
mpl.rcParams['font.sans-serif'] = [u'SimHei']  #FangSong/黑体 FangSong/KaiTi
mpl.rcParams['axes.unicode_minus'] = False

# 5.1 绘制正态分布概率密度函数
mu = 0
sigma = 1
x = np.linspace(mu - 3 * sigma, mu + 3  * sigma, 51)
y = np.exp(-(x - mu) ** 2 / (2 * sigma ** 2)) / (math.sqrt(2 * math.pi) * sigma)
print(x.shape)
print('x = \n', x)
print(y.shape)
print('y = \n', y)
plt.figure(facecolor='w')
plt.plot(x, y, 'ro-', linewidth=2)
plt.xlabel('X', fontsize=15)
plt.ylabel('Y', fontsize=15)
plt.title(u'正态函数分布', fontsize=18)
plt.grid(True)
plt.show()

(51,)
x = 
 [-3.   -2.88 -2.76 -2.64 -2.52 -2.4  -2.28 -2.16 -2.04 -1.92 -1.8  -1.68 -1.56 -1.44 -1.32 -1.2  -1.08 -0.96 -0.84 -0.72 -0.6  -0.48 -0.36 -0.24 -0.12  0.    0.12  0.24  0.36  0.48  0.6   0.72  0.84
  0.96  1.08  1.2   1.32  1.44  1.56  1.68  1.8   1.92  2.04  2.16  2.28  2.4   2.52  2.64  2.76  2.88  3.  ]
(51,)
y = 
 [ 0.00443185  0.00630673  0.00884645  0.01223153  0.0166701   0.02239453  0.02965458  0.03870686  0.04980009  0.06315656  0.07895016  0.09728227  0.1181573   0.14145997  0.16693704  0.19418605
  0.2226535   0.25164434  0.28034381  0.30785126  0.3332246   0.35553253  0.37391061  0.38761662  0.39608021  0.39894228  0.39608021  0.38761662  0.37391061  0.35553253  0.3332246   0.30785126
  0.28034381  0.25164434  0.2226535   0.19418605  0.16693704  0.14145997  0.1181573   0.09728227  0.07895016  0.06315656  0.04980009  0.03870686  0.02965458  0.02239453  0.0166701   0.01223153
  0.00884645  0.00630673  0.00443185]

#5.2 损失函数：Logistic损失(-1,1)/SVM Hinge损失/ 0/1损失
plt.figure(figsize=(10,8), facecolor='w')
x = np.linspace(start=-2, stop=3, num=1001, dtype=np.float)
y_logit = np.log(1 + np.exp(-x)) / math.log(2)
y_boost = np.exp(-x)
y_01 = x < 0
y_hinge = 1.0 - x
y_hinge[y_hinge < 0] = 0
plt.plot(x, y_logit, 'r-', label='Logistic Loss',linewidth=2)
plt.plot(x, y_01, 'g-', label='0/1 Loss', linewidth=2)
plt.plot(x, y_hinge, 'b-', label='Hinge Loss', linewidth=2)
plt.plot(x, y_boost, 'm--', label='Adaboost Loss', linewidth=2)
plt.grid()
plt.legend(loc='upper right')
plt.savefig('1.png')
plt.show()

#5.3 x^x
plt.figure(facecolor='w')
x = np.linspace(-1.3, 1.3, 101)
#y = f(x)
y = x
plt.plot(x, y, 'g-', label='x^x', linewidth=2)
plt.grid()
plt.legend(loc='upper left')
plt.show()

### 5.4 胸型线
x = np.arange(1, 0, -0.01)
y = (-3 * x *np.log(x) + np.exp(-(40 * (x - 1 / np.e)) ** 4) / 25) / 2
plt.figure(figsize=(5, 7),facecolor='w')
plt.plot(y, x, 'r-', linewidth=2)
plt.grid(True)
plt.title(u'胸型线', fontsize=20)
plt.savefig('breast.png')
plt.show()

# 5.5 心形线
t = np.linspace(0, 2*np.pi, 100)
x = 16 * np.sin(t) ** 3
y = 13 * np.cos(t) - 5 * np.cos(2*t) - 2 * np.cos(3*t) - np.cos(4*t)
plt.plot(x, y, 'r-', linewidth=2)
plt.grid(True)
plt.show()

# # 5.6 渐开线
t = np.linspace(0, 50, num=1000)
x = t*np.sin(t) + np.cos(t)
y = np.sin(t) - t*np.cos(t)
plt.plot(x, y, 'r-', linewidth=2)
plt.grid()
plt.show()

# Bar
x = np.arange(0, 10, 0.1)
y = np.sin(x)
plt.bar(x, y, width=0.04, linewidth=0.2)
plt.plot(x, y, 'r--', linewidth=2)
plt.title(u'Sin曲线')
plt.xticks(rotation=0)
plt.xlabel('X')
plt.ylabel('Y')
plt.grid()
plt.show()

# 6.1 均匀分布
x = np.random.rand(10000)
t = np.arange(len(x))
plt.hist(x, 30, color='m', alpha=0.5, label=u'均匀分布')
# plt.plot(t, x, 'g.', label=u'均匀分布')
plt.legend(loc='upper left')
plt.grid()
plt.show()

# # 6.2 验证中心极限定理
t = 1000
a = np.zeros(10000)
for i in range(t):
    a += np.random.uniform(-5, 5, 10000)
a /= t
plt.hist(a, bins=30, color='g', alpha=0.5, normed=True, label=u'均匀分布叠加')
plt.legend(loc='upper left')
plt.grid()
plt.show()

# 6.21 其他分布的中心极限定理
lamda = 7
p = stats.poisson(lamda)
y = p.rvs(size=1000)
mx = 30
r = (0, mx)
bins = r[1] - r[0]
plt.figure(figsize=(15, 8), facecolor='w')
plt.subplot(121)
plt.hist(y, bins=bins, range=r, color='g', alpha=0.8, normed=True)
t = np.arange(0, mx+1)
plt.plot(t, p.pmf(t), 'ro-', lw=2)
plt.grid(True)

N = 1000
M = 10000
plt.subplot(122)
a = np.zeros(M, dtype=np.float)
p = stats.poisson(lamda)
for i in np.arange(N):
    a += p.rvs(size=M)
a /= N
plt.hist(a, bins=20, color='g', alpha=0.8, normed=True)
plt.grid(b=True)
plt.show()

# 6.3 Poisson分布
x = np.random.poisson(lam=5, size=10000)
print(x)
pillar = 15
a = plt.hist(x, bins=pillar, normed=True, range=[0, pillar], color='g', alpha=0.5)
plt.grid()
plt.show()
print(a)
print(a[0].sum())

[2 8 4 ..., 7 3 8]

(array([ 0.0059,  0.0331,  0.0798,  0.1407,  0.1771,  0.1779,  0.1488,  0.1048,  0.066 ,  0.0356,  0.0175,  0.0081,  0.0029,  0.0011,  0.0007]), array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,  11.,  12.,  13.,  14.,  15.]), )
1.0