#!/usr/bin/python3
#coding=utf-8
# 使用图解法（穷举法）寻找线性规划问题的最优解
# 其中，线性方程组使用高斯消元法求解

import numpy as np

## ============================================
## ======== 以下为高斯消元法的算法部分 ========
# 参考代码：https://blog.csdn.net/lzyws739307453/article/details/89816311

# 打印矩阵的内容
def printM(a):
    m = len(a)
    n = len(a[0])
    for i in range(0,m):
        for j in range(0,n):
            print("%10f,\t"%a[i][j],end="")
        print()

# 选择列主元并消元
def SelectColE(a): #a:2-dimension matrix; n:int 
    n = len(a) # 参数n应该代表的是矩阵的行数，也就是线性方程组的方程个数
    for i in range(0,n): # 循环处理每一行的信息
        r = i
        for j in range(i,n):
            if(abs(a[j][i])>abs(a[r][i])): r=j
        if(r!=i): # 当r!=i时，交换两行的内容。这里的操作基于以下逻辑：
            # 由于单纯的高斯消元法会将矩阵化为梯形矩阵，靠下方的行中系数为0的项更多
            # 而非0项也会有很多分母，在除法运算中会导致精度降低
            # 因此将系数大的项尽可能的往前移动，可以提高精度
            tmp  = a[r]
            a[r] = a[i]
            a[i] = tmp
        for j in range(i+1,n): # 消元
            temp = a[j][i]/a[i][i]
            for k in range(i,n+1):
                a[j][k] -= a[i][k]*temp
        #print("第{}列消元后：".format(i))
        #printM(a)

# 高斯消元法
def Gauss(a): 
    n = len(a) # 参数n应该代表的是矩阵的行数，也就是线性方程组的方程个数
    SelectColE(a)
    #print("上三角矩阵的结果：")
    #printM(a)
    for i in range(n-1,-1,-1): # 回代求解
        # i 从 n-1递减循环到0，可以取到0
        # 从第n-1行循环到第0行。由于前面已经对矩阵进行了变换，得到了一个倒三角矩阵
        # 所以n-1行可以直接求得第n-1个未知数的值（未知数的下标从0开始）
        # 而前面的行则需要经过一些迭代
        for j in range(i+1,n):  # 对行矩阵进行回代
            # 在此行的最后一个元素处存储回代后的值
            # a[i][n]是行矩阵的最后一个元素，是存储回代后的数值的地方
            # a[i][i]是本行对应的未知数的系数。一共i行，那么就是i个未知数
            # 例如一个3*4的矩阵，有3个变量分别是x0,x1,x2，那么当i=1时，
            # 对应位置就是变量x1的系数
            # 同样的，a[i][j]是本行第j个未知数的系数，
            # 而第j个未知数的值刚刚我们已经求过了，是a[j][n]，
            # 所以a[i][n]减去他俩相乘的值就行
            a[i][n] -= a[i][j]*a[j][n]
        # 回代结束后，对系数进行化简，使得第i个未知数的系数变为1，即可
        # 于是第i个未知数的值就等于a[i][n]/a[i][i]。我们把它存储在a[i][n]的位置上
        a[i][n] /= a[i][i]

## ======== 以上为高斯消元法的算法部分 ========
## ============================================



# 回溯法求解组合情况枚举问题
# 参考自B站
# https://www.bilibili.com/video/BV1Lg411P7k3
def dfs(a,step,n,r,res):
    # a: 一个数组，用于存储组合方案。需要预先对数组大小初始化
    # step: 当前进行到了第几步
    # n: 元素的总体个数
    # r: 选取几个元素。注意，r<=n
    # res: 一个用于存储所有结果的公共变量
    for i in range(a[step-1]+1,n+1):
        a[step] = i
        if(step==r): #到达目的地
            res.append(a.copy())
        else:
            dfs(a,step+1,n,r,res)
def run_dfs(n,r):
    # n: 元素的总体个数
    # r: 选取几个元素。注意，r<=n
    a = np.zeros(r+1,dtype=int)
    res = []
    dfs(a,1,n,r,res)
    return res

# 判断一个解是否为可行解
def feasible_solution(A,b,X):
    # 传入参数：
    #  - A：m*n的矩阵，参考graph_method(A,b,C)的定义
    #  - b：长度为m的向量，参考graph_method(A,b,C)的定义
    #  - X：长度为n的向量，是决策变量向量的一组基解
    # 我们的任务是判断这一组基解是否为可行解
    print('子程序：判断一组基解是否可行')
    #print('A=',A)
    #print('X=',X)
    for x in X:
        if(x<0):return False
    b1 = np.dot(A,X)
    for i in range(len(b)):
        if(b1[i]!=b[i]):return False
    return True

# 图解法求解线性规划问题
def graph_method(A,b,C): 
    # 传入参数包括一个矩阵A和一个向量b，以及一个向量C
    #  - 矩阵A=(P_1,P_2,...,P_n)，这个矩阵大小为m*n，其中m为约束条件方程的个数，
    #    n为决策变量的变量个数（包括松弛变量）。注意m<n
    #  - 向量b是资源向量，其长度为m，与约束条件方程的个数一致
    #  - 向量C是价值向量，其长度为n，与决策变量个数一致
    #
    # 注意，在调用此函数前，必须将线性规划问题化为带有松弛变量的标准形式！
    # 线性规划问题的标准形式的矩阵表达：
    # max z=CX, 
    #        其中：
    #           C为目标函数的系数向量，或者叫做价值向量；
    #           X是未知数向量，或者决策变量向量
    # s.t. sum(P_j*x_j)=b
    #      x_j >= 0, j=1,2,...,n
    #        其中：
    #           P_j是第j个未知数（决策变量）在所有约束条件中的系数向量
    #           x_j是第j个未知数（决策变量）
    #           b是资源向量，可以理解为约束条件不等式中的常数项
    #
    # 我们的求解思路是，从矩阵A中得到所有可能的基矩阵B，后者是一个m*m的矩阵
    # 并利用这个基和资源向量b求出一个基解
    # 之后，使用约束条件判断这个基解是否是可行解，这样我们得到了一个基可行解的集合
    # 我们在这个可行解的集合中找到最大的那个解，返回这个解的决策变量向量和解的数值
    print("\n\n##############\n要求解的约束向量A=\n",A)
    m = len(b) # 约束条件个数
    n = len(C) # 决策变量个数
    r = run_dfs(n,m) # 通过回溯法，寻找所有可能的基矩阵的组合方式
    s = []     # 一个存储基可行解的列表 
    # 下面这个循环从所有基中获得基解，然后判断基解是否是可行解
    for ri in r:
        print("\n\n>>> 尝试计算一个基解")
        tmp = [] # 一个临时变量，用于存储基向量
        for i in ri[1:]: # 注意，上述dfs函数的返回值中，每个组合是以数组的[1:]存储的
            tmp.append(A[:,i-1].copy().T) # 用这个方式获得基向量.i要记得减1
        B  = np.array(tmp).T # 这就获得了一个基
        # np.column_stack((B,b)) 是构建一个增广向量
        Bb = np.column_stack((B,b)).astype(float).tolist()
        print('Bb=',np.array(Bb))
        try: # 如果线性方程组无解：
        #if(True):
            Gauss(Bb)  
            print('Bb=',Bb)
            X = np.zeros(n)
            for j in range(m):
                X[ri[j+1]-1] = Bb[j][-1]
            #X  = np.array(Bb)[:,-1]
            #print('X=',X)
            if(feasible_solution(A,b,X)):
                print("一个可行解为：\n{}".format(X))
                s.append(X)
            else:
                print("这个基解不是可行解：\n{}".format(X))
        except: # 线性方程组无解：
            print("这个基没有解：\n{}".format(B))
        
        # TODO:
        # 1. 基可行解的形式是X=(x_1,x_2,...,x_m,0,0,...0)，这是因为使用一个基（m*m矩阵）能够求出来的解的个数只有m个，然而约束条件共n个（n>m)，所以剩下的要在对应位置填0。代码中这一部分的逻辑有问题，请重新编写【update：已完成】
        # 2. 高斯消元法在之前的测试中确实成功求出来了一些解，不知道为什么现在总是出现NaN的情况，请检查相关的逻辑【update：是numpy和list不兼容的问题。已完成】
        # 3. 在判断一个基解是可行解之后，应该再去进行比较，以获得目标函数z最大的可行解。请完善这一部分的逻辑【这是下一步的目标】






if(__name__=='__main__'):
    # 下面的内容是高斯消元法求解线性方程组的测试
    # 最终的方程组的解为a[:,-1]
    """
    a =[[ 1,-1,-1, 2 ],
        [ 2,-1,-3, 1 ],
        [ 3, 2,-5, 0 ]]
    Gauss(a.copy())
    print("求解结果：")
    printM(a)
    print("数据结构：")
    printM(a)
    for i in range(0,len(a)):
        print("X%d = %9f"%(i,a[i][3]))
    print(np.array(a)[:,-1])
    """

    # 书上的例题1
    A =[[1, 2, 1, 0, 0],
        [4, 0, 0, 1, 0],
        [0, 4, 0, 0, 1]]
    b = [8, 16,12]
    C = [2, 3, 0, 0, 0]
    # 本题的解为 X=[4,2,*,*,*], z_max=14
    graph_method(np.array(A),b,C) # 图解法求解线性规划问题