import numpy as np 
num_obs = 100
x1 = np.random.uniform(0,10,num_obs)
x2 = np.random.uniform(0,10,num_obs)
const = np.ones(num_obs)
eps = np.random.normal(0,.5,num_obs)

b = 1.5
theta_1 = 2
theta_2 = 5

y = b*const + theta_1*x1 + theta_2*x2 + eps

x_mat = np.array([const,x1,x2]).T

learning_rate = 1e-3
num_iter = 2000
theta_initial = np.array([3,3,3])

def gradient_descent(learning_rate,num_iter,theta_initial):
    ## initialization steps
    theta = theta_initial
    theta_path = np.zeros((num_iter+1,3))  # 两位数组，存放参数轨迹
    theta_path[0,:] = theta_initial

    loss_vec = np.zeros(num_iter)

    for i in range(num_iter):
        y_pred = np.dot(theta.T,x_mat.T)
        loss_vec[i] = np.sum((y-y_pred)**2) # 记录每一轮迭代的loss
        grad_vec = (y-y_pred).dot(x_mat)
        grad_vec /= num_obs
        # 随机梯度下降
        # for j in range(len(y_pred)):
        #     grad_vec = (y[j]-y_pred[j])*(x_mat[j,:])
        #     theta = theta+learning_rate*grad_vec
        theta = theta + learning_rate*grad_vec
        theta_path[i+1,:] = theta
    return theta_path, loss_vec 

import matplotlib.pyplot as plt
theta_path,loss_vect = gradient_descent(learning_rate,num_iter,theta_initial)
plt.plot(loss_vect)
plt.show()

plt.plot(theta_path[:,0])
plt.show()

plt.plot(theta_path[:,1])
plt.show()

plt.plot(theta_path[:,2])
plt.show()

# Sigmoid
def sigmoid(x):
    return 1.0/(1.0+np.exp(-x))

def logic_gate(w1,w2,b):
    return lambda x1,x2:sigmoid(w1*x1+w2*x2+b)

def test(gate):
    for a,b in (0,0),(0,1),(1,0),(1,1):
        print('{} , {} : {}'.format(a,b,np.round(gate(a,b))))

or_gate = logic_gate(20,20,-10)
test(or_gate)
print("#")
nand_gate = logic_gate(-10,-10,14)
test(nand_gate)
print("#")
and_gate = logic_gate(10,10,-19)
test(and_gate)
print("#")

# XOR问题：
def xor_gate(a,b):
    c = or_gate(a,b)
    d = nand_gate(a,b)
    return and_gate(c,d)
test(xor_gate)

W_1 = np.array([[2,-1,1,4],[-1,2,-3,1],[3,-2,-1,5]])
W_2 = np.array([[3,1,-2,1],[-2,4,1,-4],[-1,-3,2,-5],[3,1,1,1]])

W_3 = np.array([[-1,3,-2],[1,-1,-3],[3,-2,2],[1,2,1]])
x_in = np.array([.5,.8,.2])
x_mat_in = np.array([[.5,.8,.2],[.1,.9,.6],[.2,.2,.3],[.6,.1,.9],[.5,.5,.4],[.9,.1,.9],[.1,.8,.7]])

def soft_max_vec(vec):
    return np.exp(vec)/(np.sum(np.exp(vec)))

print("Matrix W_1\n")
print(W_1)
print("-"*30)
print("vector input x_in\n")
print(x_in)
print("-"*30)

z_2 = np.dot(x_in,W_1)
z_2

a_2 = sigmoid(z_2)
a_2

z_3 = np.dot(a_2,W_2)
a_3 = sigmoid(z_3)

z_4 = np.dot(a_3,W_3)
y_out = soft_max_vec(z_4)
y_out

def soft_max_mat(mat):
    return np.exp(mat)/(np.sum(np.exp(mat),axis=1).reshape(-1,1))

print('batch矩阵输入 -- starts with the x_mat_in\n')
print(x_mat_in)

z_2 = np.dot(x_mat_in,W_1)
z_2

a_2 = sigmoid(z_2)
z_3 = np.dot(a_2,W_2)
a_3 = sigmoid(z_3)
z_4 = np.dot(a_3,W_3)
y_out = soft_max_mat(z_4)
y_out

# 生成数据
num_obs = 500
x_mat_1 = np.random.uniform(-1,1,size = (num_obs,2))# -1~1均匀分布
x_mat_bias = np.ones((num_obs,1))
x_mat_full = np.concatenate((x_mat_1,x_mat_bias),axis = 1)# 三维数据

# x,y坐标值绝对值较大的小于0.5，返回True
# y = ((np.maximum(np.abs(x_mat_full[:,0]),np.abs(x_mat_full[:,1])))<.5).astype(int)

# 距离原点距离1的菱形区域
y = ((np.abs(x_mat_full[:,0])+np.abs(x_mat_full[:,1]))<1).astype(int)
print(y)

fig,ax = plt.subplots(figsize = (5,5))
ax.plot(x_mat_full[y==1,0],x_mat_full[y == 1,1],"ro",label = 'class 1',color = "darkslateblue")
ax.plot(x_mat_full[y==0,0],x_mat_full[y == 0,1],"bx",label = 'class 0',color = "chocolate")
ax.grid(True)
ax.legend(loc = "best")
ax.axis('equal')

# 模型：两层网络，sigmoid激活函数
# 向前传播与反向传播
def forward_pass(W_1,W_2):
    global x_mat
    global y
    global num_

    # 向前传播
    z_2 = np.dot(x_mat,W_1)
    a_2 = sigmoid(z_2)
    z_3 = np.dot(a_2,W_2)
    y_pred = sigmoid(z_3).reshape((len(x_mat),))
    # 反向传播
    J_z_3_grad = -y+y_pred
    J_W_2_grad = np.dot(J_z_3_grad,a_2)
    a_2_z_2_grad = sigmoid(z_2)*(1-sigmoid(z_2))
    J_W_1_grad = (np.dot((J_z_3_grad).reshape(-1,1),W_2.reshape(-1,1).T)*a_2_z_2_grad).T.dot(x_mat).T

    # 记录多层的梯度
    gradient = (J_W_1_grad,J_W_2_grad)

    return y_pred,gradient

def loss_fn(y,y_pred):
    return np.log(np.sum((y-y_pred)**2))

def plot_loss_accuracy(loss_vals,accuracies):
    plt.figure(figsize = (13,6))
    plt.figure(1)
   
    ax_1 = plt.subplot(121)
    ax_1.set_title("Accuracies")
    plt.plot(accuracies)

    ax_2 = plt.subplot(122)
    ax_2.set_title("Loss")
    plt.plot(loss_vals)

    plt.suptitle("Accuracies and Loss")

W_1 = np.random.uniform(-1,1,size = (3,4))
W_2 = np.random.uniform(-1,1,size = (4))
num_iter = 20000
learning_rate = .001
x_mat = x_mat_full

loss_vals,accuracies = [],[]
for i in range(num_iter):
    # 向前传播并反向修正
    y_pred,(J_W_1_grad,J_W_2_grad) = forward_pass(W_1,W_2)
    # 调整模型参数
    W_1 = W_1 - learning_rate*J_W_1_grad
    W_2 = W_2 - learning_rate*J_W_2_grad
    # 计算loss和正确率
    curr_loss = loss_fn(y,y_pred)
    loss_vals.append(curr_loss)
    acc = np.sum((y_pred >= .5)==y)/num_obs
    accuracies.append(acc)

    # 每500次迭代打印一次loss和正确率
    if ((i%2000)==0):
        print('iteration {}, log loss is {:.4f}, accuracy is {}'.format(
            i,curr_loss,acc
        ))

plot_loss_accuracy(loss_vals,accuracies)

# 可视化结果
pred1 = (y_pred>=0.5)
pred0 = (y_pred<0.5)
fig,ax = plt.subplots(figsize = (8,8))
# true predictions
ax.plot(x_mat[pred1 & (y==1),0],x_mat[pred1 & (y==1),1],'ro',label = "true positives")
ax.plot(x_mat[pred0 & (y==0),0],x_mat[pred0 & (y==0),1],'bx',label = "true negatives")
# false predictions
ax.plot(x_mat[pred1 & (y==0),0],x_mat[pred1 & (y==0),1],'yx',label = "false positives",markersize=15)
ax.plot(x_mat[pred0 & (y==1),0],x_mat[pred0 & (y==1),1],'yo',label = "false negatives",markersize=15,alpha = .6)
ax.set(title = "Truth vs Prediction")
ax.legend(bbox_to_anchor = (1,0.8),fancybox = True,shadow = True,fontsize = 'x-large')

# Tensor 的声明与初始化
import torch
x = torch.zeros(5,3,dtype = torch.long)
print(x)

x = torch.tensor([5.5,3])
print(x)

# "空Tensor"
x = torch.empty(5,3)
print(x)
y = torch.empty(5,3)
y = torch.randn_like(x,dtype = torch.float)
print(y)
print(id(x))
print(id(y))

x = torch.rand(5,3) # 重新生成
 
print(x)
print(id(x))

# 运算
y = torch.ones(5,3,dtype = torch.int)
print(x+y)
print(x.add_(y))
# 自动转型
# print(y.add_(x)) # result type Float can't be cast to the desired output type Int

# 切片与广播
x[2:] = torch.tensor([1,2,3])
print(x)
y = torch.tensor([0,1,2])
print(x*y)

# 与Numpy ndarray 共享数据定义
a = torch.ones(5,dtype = torch.int)
print(a)
b = a.numpy()
print(b)
b += 1
print(a)
a[3] = 6
print(b)

# 类型自动升级，赋值语句生成新的对象
b = b + 1.0
print(b)
print(a)

# GPU 加速
if torch.cuda.is_available():
    device = torch.device('cuda')
    y = torch.ones_like(x,device = device)
    x = x.to(device)
    z = x+y
    print(z)
    print(z.to('cpu',torch.double))
else:
    print("No cuda avaiaable")

# Autograd, NN模型回归机制,requires_grad = True用于跟踪计算
x = torch.ones(2,2,requires_grad = True)
print(x)
y = x+1
print(y)
z = y*2
print(z)

out = z.mean()
print(out)

# out.backward(retain_graph=True)
# print(x.grad)
# #print(y.grad)