# -*- coding: utf-8 -*-
# author: 'boliang'
# date: 2018/5/15 21:20

import numpy as np

class Solution(object):
    def dfp(self, fun, gfun, x0):
        print('DFP Algorithm')
        rho = 0.55
        sigma = 0.4
        epsilon = 1e-5
        maxk = 10000
        k = 1
        n = np.shape(x0)[0]
        Hk = np.eye(n)
        xk = x0
        while k < maxk:
            gk = np.mat(gfun(xk))
            if np.linalg.norm(gk) <= epsilon:
                break
            dk = -np.mat(Hk)*gk
            m = 0
            mk = 0
            while m < 20:
                tmp_a = fun(xk + (rho**m) * dk)
                tmp_b = fun(xk) + sigma * (rho**m) * (gk.T*dk)[0, 0]
                if tmp_a <= tmp_b:
                    mk = m
                    break
                m += 1

            x = xk + rho**mk * dk
            sk = x - xk
            yk = gfun(x) - gk
            if sk.T*yk > 0:
                Hk = Hk - (Hk*yk*yk.T*Hk)/(yk.T*Hk*yk) + (sk*sk.T)/(sk.T*yk)
            k += 1
            xk = x

        val = fun(xk)
        print('迭代次数: {0}'.format(k))
        print('x1: {0}\nx2: {1}'.format(xk[0, 0], xk[1, 0]))
        print('result: {0}'.format(val))

    def bfgs(self, fun, gfun, x0):
        print('BFGS Algorithm')
        rho = 0.55
        sigma = 0.4
        epsilon = 1e-5
        maxk = 100
        k = 1
        n = np.shape(x0)[0]
        Bk = np.eye(n)
        xk = x0
        while k < maxk:
            gk = np.mat(gfun(xk))
            if np.linalg.norm(gk) <= epsilon:
                break
            dk = np.mat(-np.linalg.solve(Bk, gk))
            m = 0
            mk = 0
            while m < 20:
                tmp_a = fun(xk + (rho ** m) * dk)
                tmp_b = fun(xk) + sigma * (rho ** m) * (gk.T * dk)[0, 0]
                if tmp_a <= tmp_b:
                    mk = m
                    break
                m += 1

            x = xk + rho ** mk * dk
            sk = x - xk
            yk = gfun(x) - gk
            if yk.T*sk > 0:
                Bk = Bk - (Bk*sk*sk.T*Bk) / (sk.T*Bk*sk) + (yk*yk.T) / (yk.T*sk)
            k += 1
            xk = x

        val = fun(xk)
        print('迭代次数: {0}'.format(k))
        print('x1: {0}\nx2: {1}'.format(xk[0, 0], xk[1, 0]))
        print('result: {0}'.format(val))


def fun(x):
    return 100*(x[0, 0]**2-x[1, 0])**2 + (x[0, 0] - 1)**2

def gfun(x):
    a = 400*x[0, 0] * (x[0, 0]**2 - x[1, 0]) + 2 * (x[0, 0] - 1)
    b = -200*(x[0, 0]**2 - x[1, 0])
    result = np.zeros((2, 1))
    result[0, 0] = a
    result[1, 0] = b
    return result

def test(flag):
    sol = Solution()
    if flag == 'dfp':
        sol.dfp(fun, gfun, np.mat([[-1.2], [1]]))
    else:
        sol.bfgs(fun, gfun, np.mat([[-1.2], [1]]))

if __name__ == '__main__':
    test('dfp')
    print()
    test('bfgs')