from __future__ import unicode_literals
from __future__ import division
%matplotlib inline
from IPython.display import *
import numpy as np
import scipy.optimize as op
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import sympy as sp

sp.init_printing(use_latex=True, use_unicode=True)
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']
plt.rcParams['font.size'] = 14
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.figsize'] = 8, 6

def ex1(data_file_path):
    data = np.loadtxt(data_file_path, delimiter=',')
    x, y = data[:,0:-1], data[:, -1:]
    bias = np.ones((x.shape[0], 1))
    x_ = np.hstack((bias, x))
    theta = np.zeros(x_.shape[1])

    # 算法
    def h(x, theta):
        '''hypothesis'''
        return np.dot(x, theta).reshape(-1, 1)


    def J(x, y, theta):
        '''cost'''
        return np.power(h(x, theta) - y, 2).mean()/2


    def gradient_descent(x, y, theta, alpha, iterations):
        thetas = [theta]
        costs = []
        for i in range(iterations):
            delta = alpha * np.multiply(h(x, theta)-y, x).mean(0)
            theta = theta - delta
            thetas.append(theta)
            cost = J(x, y, theta)
            costs.append(cost)
        return thetas, costs


    def normal_equation(x, y):
        x = np.matrix(x)
        y = np.matrix(y)
        return np.linalg.pinv(x.T * x) * x.T * y

    # 数据结果图
    thetas, costs = gradient_descent(x_, y, theta, 0.01, 1500)
    theta_gradient_descent = thetas[-1]
    theta_normal_equation = normal_equation(x_, y)
    theta_minimize = op.minimize(lambda theta: J(x_, y, theta), [0, 0])['x']

    plt.plot(costs)
    plt.title('代价迭代收敛')
    plt.xlabel('迭代次数')
    plt.ylabel('代价')
    plt.show()


    display(Math(r'$\theta_{梯度下降} = \begin{bmatrix}%.2f\\%.2f\end{bmatrix}\
                 ,\ \theta_{正规函数} = \begin{bmatrix}%.2f\\%.2f\end{bmatrix}\
                 ,\ \theta_{最小化BFGS} = \begin{bmatrix}%.2f\\%.2f\end{bmatrix}$' % 
                 (theta_gradient_descent[0], theta_gradient_descent[1], 
                  theta_normal_equation[0], theta_normal_equation[1],
                  theta_minimize[0], theta_minimize[1])))

    y_35 = h([1, 3.5], theta_gradient_descent)[0, 0]
    y_7 = h([1, 7], theta_gradient_descent)[0, 0]

    plt.plot(x, y, 'x', label='样本')
    px = range(2, 30, 2)
    py = [h(np.array([1, i]), theta_gradient_descent)[0, 0] for i in px]
    plt.plot(px, py, label='线性回归模型')
    plt.plot([3.5, 7], [y_35, y_7], 'o', label='预测')
    plt.title('线性回归')
    plt.xlabel('食品车利润(万)')
    plt.ylabel('城市人口(万)')
    plt.legend(loc=4)
    plt.show()
    display(Math(r'h(%s) = %.2f' % (3.5, y_35)))
    display(Math(r'h(%s) = %.2f' % (7, y_7)))

    # 代价函数图
    theta0s = np.linspace(-10, 10, 100)
    theta1s = np.linspace(-1, 4, 100)
    theta0s, theta1s = np.meshgrid(theta0s, theta1s)
    m, n = theta0s.shape
    hs = np.array([J(x_, y, [theta0s[i, j], theta1s[i, j]]) 
                   for i in range(m) for j in range(n)]).reshape(m, n)

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(theta0s, theta1s, hs,
                    rstride=3, cstride=3, cmap=mpl.cm.jet, linewidth=0.3, antialiased=True)
    plt.xlabel(r'$\theta_0$')
    plt.ylabel(r'$\theta_1$')
    plt.title(r'代价函数 $J(\theta)$')
    plt.show()

    plt.contour(theta0s, theta1s, hs, np.logspace(-2, 3, 30))
    plt.xlabel(r'$\theta_0$')
    plt.ylabel(r'$\theta_1$')
    plt.title(r'代价函数等高线')
    plt.plot([i[0] for i in thetas], [i[1] for i in thetas], 'r-')
    plt.plot(theta_gradient_descent[0], theta_gradient_descent[1], 'rx')
    plt.plot(theta_normal_equation[0, 0], theta_normal_equation[1, 0], 'k+')
    plt.plot(theta_minimize[0], theta_minimize[1], 'yx')
    plt.show()


ex1('ex1data1.txt')

z = sp.symbols('z')
sigmoid = 1/(1+sp.exp(-z))
sp.plotting.plot(sigmoid, (z, -8, 8))
display(sigmoid)

def ex2(data_file_path):
    data = np.loadtxt(data_file_path, delimiter=',')
    x, y = data[:,0:-1], data[:, -1:]
    bias = np.ones((x.shape[0], 1))
    x_ = np.hstack((bias, x))
    theta = np.zeros(x_.shape[1])


    # 算法
    def sigmoid(z):
        return 1 / (1+np.exp(-z))

    def h(x, theta):
        return sigmoid(np.dot(x, theta).reshape(-1, 1))

    def J(x, y, theta):
        h_val = h(x, theta)
        return (-y * np.log(h_val) - (1-y)*np.log(1-h_val)).mean();

    def gradient(x, y, theta):
        return np.multiply(h(x, theta)-y, x).mean(0)

    print('初始代价 %.3f' % J(x_, y, theta))

    cost = lambda t: J(x_, y, t)
    gradient_ = lambda t: gradient(x_, y, t)

    # TODO: make it work with other minimize algorithm
    #op.minimize(lambda t: J(x_, y, t), 
    #            np.zeros(x_.shape[1]), method = 'TNC', jac=lambda t: gradient(x_, y, t))
    #print gradient_(init)
    #print 'fmin_cg()\n', op.fmin_cg(cost, init, fprime=gradient_)
    #op.check_grad(cost, gradient_, init)

    theta_fmin = op.fmin(cost, theta, disp=False)
    display(Math(r'\theta_{fmin} = \begin{bmatrix}%.3f\\%.3f\\%.3f\end{bmatrix}' % tuple(theta_fmin)))

    pos = np.where(y == 0)[0]
    neg = np.where(y != 0)[0]
    plt.plot(x[pos, 0], x[pos, 1], 'yo', label='拒绝')
    plt.plot(x[neg, 0], x[neg, 1], 'k+', label='录取')
    plt.xlabel('甲科考分')
    plt.ylabel('乙科考分')
    plot_x = np.array([x[:,0].min()-2, x[:,0].max()+2])
    plot_y =  (theta_fmin[0] + theta_fmin[1]*plot_x) / -theta_fmin[2]
    plt.plot(plot_x, plot_y, label='决策边界') # TODO how to plot decision bondary generally?
    plt.legend(loc=3)
    plt.show()

    print('两门成绩 45 分和 85 分的考生录取率为 {:.1%}'.format(h([1, 45, 85], theta_fmin)[0][0]))
    print('总拟合准确率 {:.1%}'.format(((h(x_, theta_fmin) >= 0.5) == y).mean()))


ex2('ex2data1.txt')