#!/usr/bin/env python
# coding: utf-8

# In[2]:


get_ipython().run_line_magic('qtconsole', '')
get_ipython().run_line_magic('load_ext', 'autoreload')
get_ipython().run_line_magic('autoreload', '2')


# In[4]:


import numpy as np


# ## The simplex algorithm
# 
# Implementing the algorithm seen during the lectures

# In[1]:


def simplexit(A,b,c, init, debug=True, expectNullOptimum=False):
    """
    Implements the standard matrix form of the simplex algorithm

    Hugues Talbot 2018, from an initial version in R ca. 2009.
    """
    iter = 0
    dim=A.shape[0]
    nbvar=A.shape[1]
    lastx=0
    lastxl=0
    sol=None
    IBVmem=None # to detect cycles

    vars=np.arange(0,nbvar) # indices start at 0
    IBV=init # Initial basis vector
    while (True):
        if debug:
            print("*** Iteration: ", iter)
        NBV = np.setdiff1d(vars,IBV) # Set complement -> non basis variables
        if debug:
            print("    IBV = ", IBV)
            print("    NBV = ", NBV)

        B=A[:,IBV]
        E=A[:,NBV]
        cb = c[IBV]
        ce = c[NBV]

        if debug:
            print("    B = \n", B)
            print("    E = \n", E)

        try:
            iB = np.linalg.inv(B)
        except:
            print("B is not invertible, error")
            break

        if debug:
            print("    iB = \n", iB)

        bbar = iB.dot(b)
        if (np.min(bbar) < -1e-7):
            print("bbar is not positive, not a feasible solution")
            break

        if debug:
            print("    bbar= ", bbar)

        cebart = ce.T - cb.T.dot(iB).dot(E)

        if debug:
            print("  reduced costs", cebart)

        if (np.min(cebart) >= 0):
            print(" optimum found!")
            sol = (IBV, bbar)
            print("IBV = ", IBV)
            print("Bbar= ", bbar)
            if ((np.min(cebart) > 0) or expectNullOptimum):
                break
            if (IBVmem is None):
                IBVmem = IBV
            else:
                if (IBV == IBVmem):
                    print("Cycle detected\n")
                    break

        l = np.argmin(cebart) # may not be unique
        P = iB.dot(A[:,NBV[l]])
        if (debug):
            print("P = ", P)
        if (np.max(P) <= 0):
            print("Solution in unbounded!")
            break
        ratios = bbar/P ## elementwise division
        if (debug):
            print("Ratios = ", ratios)
        ## Ignore negative members f P
        negP = P < 0
        ratios[negP] = np.inf ## this way they will not be selected as the min
        nanP = np.isnan(ratios)
        ratios[nanP] = np.inf ## nan are eliminated as well
        if (debug):
            print("Fixed ratios = ", ratios)
        xl = np.min(ratios)
        j = np.argmin(ratios)
        lastx = xl
        lastxl = NBV[l]
        if (debug):
            print("Variable (", IBV[j],") in the basis replaced by variable (",NBV[l],")\n")
            print("xl = ",xl)

        IBV[j] = NBV[l]
        if (debug):
            print("IBV = ", IBV)

        iter += 1
        ## end while
        
    return sol


# In[2]:


## exercice 2, TD1 (production de A,B,C liée)
def exampleII(debug=True):
    # A <- matrix(data =c(-1, 0, 1, 2, -1, 2, 0, 1, 3, 1, 0, 0, 0, 1, 0, 0, 0, 1), ncol=6, nrow=3)
    A = np.array([[-1,  2, 0, 1, 0, 0], 
                  [ 0, -1, 1, 0, 1, 0], 
                  [ 1,  2, 3, 0, 0, 1]])
    b = np.array([  0,  0, 35])
    C = np.array([-10,-36,-44, 0, 0, 0])
    IBV = np.array([3,4,5]) # index of the IBV
    sol = simplexit(A,b,C,IBV,debug)
    
    # calcul optimum
    z = C[sol[0]].dot(sol[1])
    return(z)


# In[6]:


exampleII(debug=True)


# In[7]:


## Bakery problem
def bakery2(debug=True):
    # A <- matrix(data =c(-1, 0, 1, 2, -1, 2, 0, 1, 3, 1, 0, 0, 0, 1, 0, 0, 0, 1), ncol=6, nrow=3)
    A = np.array([[ 5, -1, 1, 1, 0], 
                  [ 1,  0, 0, 0, 1]])
    b = np.array([30, 5])
    C = np.array([-30, 4, -4, 0, 0])
    IBV = np.array([3,4]) # index of the IBV
    sol = simplexit(A,b,C,IBV,debug)
    
    # calcul optimum
    z = C[sol[0]].dot(sol[1])
    return(z)


# In[8]:


bakery2(debug=False)


# ### Alloy problem

# In[9]:


## Alloy problem
def alloys(debug=True):
    A = np.array([[ 10, 10, 40, 60, 30, 30, 30, 50, 20, 1, 0, 0], 
                  [ 10, 30, 50, 30, 30, 40, 20, 40, 30, 0, 1, 0],
                  [ 80, 60, 10, 10, 40, 30, 50, 10, 50, 0, 0, 1]])
    b = np.array([30, 30, 40])
    C = 100*np.array([4.1, 4.3, 5.8, 6.0, 7.6, 7.5, 7.3, 6.9, 7.3, 100, 100, 100])
    IBV = np.array([9,10,11]) # index of the IBV
    sol = simplexit(A,b,C,IBV,debug)
    
    # calcul optimum
    z = C[sol[0]].dot(sol[1])
    return(z)


# In[10]:


alloys(debug=False)


# ## Collisions
# 
# In the dual, if the solution is unbounded, this means that there is no solution in the primal

# In[11]:


## Alloy problem
def collisions(debug=True):
    A = np.array([[ 1,  1, -2, -1, -1, 2, 1, 0], 
                  [ 1, -2,  1,  1, -1, 1, 0, 1]])
    b = np.array([0, 0])
    C = -1*np.array([5, -10, -4, 6, -10, 10, 0, 0])
    IBV = np.array([6,7]) # index of the IBV
    sol = simplexit(A,b,C,IBV,debug)
    
    # calcul optimum not needed
    # z = C[sol[0]].dot(sol[1])
    return


# In[12]:


collisions(debug=False)


# #### Unbounded solution in the dual, so no solution in the primal, so no collision !

# In[14]:


## Young couple
def youngcouple(debug=False):
    A = np.array([[ 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
                  [ 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0],
                  [ 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0],
                  [ 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0],
                  [ 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0],
                  [ 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
                 ])
    b = np.array([2, 2, 1, 1, 1, 1])
    C = np.array([4.5, 7.8, 3.6, 2.9, 4.9, 7.2, 4.3, 3.1, 100, 100, 100, 100, 100, 100])
    IBV = np.array([8,9,10,11,12,13])
    sol = simplexit(A,b,C,IBV,debug)
    z= C[sol[0]].dot(sol[1])
    return(z)


# In[15]:


print(youngcouple())


# In[ ]:


## tighter problem formulation