diff --git a/OR/BigM.py b/OR/BigM.py index 92c3fe1..e770168 100644 --- a/OR/BigM.py +++ b/OR/BigM.py @@ -1,73 +1,47 @@ from pandas import DataFrame import numpy as np -from .Format import * + from .TableauBase import TableauBase +from .LP import LP +from .ObjectiveFunction import ObjectiveFunction +from .NonNeg import NonNeg -def getRep(num,M_val): - # solve for integer M values, the rest... soz - round5 = lambda x: round(x,5) - if abs(num%M_val) < abs(num%-M_val): - #kM+x - mNo = int((num - num%M_val)/M_val) - defNo = num%M_val - else: - #KM-x - mNo = int((num - num%-M_val)/M_val) - defNo = num%-M_val - - mNo = round5(mNo) - defNo = round5(defNo) +from typing import TypeVar +BigM = TypeVar("BigM") - if mNo == 1: - mString = " M " - elif mNo == -1: - mString = "-M " - elif mNo == 0: - mString = "" - else: - mString = f"{mNo}M " - - if defNo == 0: - defString = "" - elif defNo > 0: - defString = f"+{defNo}" - else: - defString = f"{defNo}" - if mString + defString == "": - return 0 - else: - return mString + defString class BigM(TableauBase): @classmethod - def fromLP(cls,inLP,M_Val = 1e6): - constraintSigns = inLP.posRHS().body.loc[:,"signs"] - baseT = super().fromLP(inLP) - body = baseT.dataFrame - basicVar = list(baseT.basicVar) - aVals = [] - for i in range(1,len(constraintSigns)): - if constraintSigns[i] in (">","="): - newCol = [0] * body.shape[0] - newCol[0] = M_Val - newCol[i] = 1 - body.insert(body.shape[1]-1,f"a{subscript(i)}",newCol) - basicVar[i] = f"a{subscript(i)}" - - return BigM(body,tuple(basicVar),aVals,M_Val=M_Val) - - def fromTableauBase(self,base): - return BigM(base.dataFrame,base.basicVar,self.aVals) + def fromBase(cls,inLP:LP) -> BigM: + return cls(inLP.dataFrame,inLP.basicVar) + + @classmethod + def fromLP(cls,inLP:LP,M_Val: float = 1e6) -> BigM: + inLP = inLP.standardForm() + for i,c in enumerate(inLP.constraints): + if not isinstance(c,NonNeg) and c.sign in ["=",">"]: + inLP.objFunc[f"A{i+1}"] = M_Val + c[f"A{i+1}"] = 1 + baseLP = super().fromLP(inLP) + return cls.fromBase(baseLP) - def __init__(self,dataFrame,basicVar,aVals,M_Val=1e6): - super().__init__(dataFrame,basicVar,dispFunc = lambda x: getRep(x,M_Val)) - self.aVals = aVals # list of a factors - # self.dispFunc = + def __init__(self,dataFrame:DataFrame,basicVar:DataFrame,M_Val: float=1e6): + self.dataFrame = dataFrame + self.basicVar = basicVar + self.M_Val = M_Val - def canonical(self,echo=False): + def canonical(self,echo=False) -> BigM: t = super().canonical(echo=echo) - return self.fromTableauBase(t) + return self.fromBase(t) + + def isInfeasible(self): + if not self.isOptimal(): + return False + elif "A" in "".join(self.getBFS().keys()): + return False + else: + return True def solve(self,echo=False): - t = super().solve(echo=echo) - return self.fromTableauBase(t) \ No newline at end of file + t = super().canonical(echo=echo).solve(echo=echo) + return self.fromBase(t) \ No newline at end of file diff --git a/OR/CommonLP.py b/OR/CommonLP.py index b459da2..3d7b3bb 100644 --- a/OR/CommonLP.py +++ b/OR/CommonLP.py @@ -1,46 +1,11 @@ +from .ObjectiveFunction import ObjectiveFunction +from .Constraint import Constraint +from .NonNeg import NonNeg from .LP import LP -class CommonLP(): - def __init__(self): - pass - - @classmethod - def bevco(cls): - return LP.new( - ["min", 2, 3], - [[0.5, 0.25, "<", 4], - [1, 3, ">", 20], - [1, 1, "=", 10], ], - [">", ">"], - factorNames=False) - - @classmethod - def bevcoInf(cls): - return LP.new( - ["min", 2, 3], - [[0.5, 0.25, "<", 4], - [1, 3, ">", 36], - [1, 1, "=", 10], ], - [">", ">"], - factorNames=False) - - @classmethod - def dakota(cls): - return LP.new( - ["max", 60, 30, 20], - [[8, 6, 1, "<", 48], - [4, 2, 1.5, "<", 20], - [2, 1.5, 0.5, "<", 8], - [0, 1, 0, "<", 5]], - [">", ">", ">"], - factorNames=False) - - @classmethod - def dakotaDual(cls): - return LP.new( - ["min", 48, 20, 8, 5], - [[8, 4, 2, 0, ">", 60], - [6, 2, 1.5, 1, ">", 30], - [1, 1.5, 0.5, 0, ">", 20], ], - [">", ">", ">", ">"], - factorNames=["Y₁", "Y₂", "Y₃", "Y₄"]) +bevCo = LP( + ObjectiveFunction.new("min", 2, 3), + Constraint.new(0.5, 0.25, "<", 4), + Constraint.new(1, 3, ">", 20), + Constraint.new(1, 1, "=", 10), + *NonNeg.fromArray(">",">",)) \ No newline at end of file diff --git a/OR/Constraint.py b/OR/Constraint.py new file mode 100644 index 0000000..75ab6cc --- /dev/null +++ b/OR/Constraint.py @@ -0,0 +1,119 @@ +from pandas import DataFrame + +from typing import List,Tuple,TypeVar +Constraint = TypeVar("Constraint") + + +class Constraint(): + invertSign = {">": "<", "=": "=", "<": ">"} + + @classmethod + def new(cls, *args, label:str="X", labels=None): + ''' + construct a new Constraint from array of values + ''' + #label tests + if labels == None: + # labels overrides label + if label.upper() in ["A", "S"]: + raise ValueError(f"A and S are reserved variables") + labels = [f"{label.upper()}{i+1}" for i in range(len(args)-2)] + elif len(labels) != len(args)-2: + raise ValueError(f"provided {len(labels)} labels for {len(args)-2} variables") + else: + label = labels[0][0] + #sign tests + if args[-2] not in ["<", ">", "="]: + raise ValueError(f"{args[-2]} is not a valid constraint sign") + + LHS = {k: v for k, v in zip(labels, args[:-2])} + sign = args[-2] + RHS = args[-1] + return cls(label, LHS, sign, RHS) + + def __init__(self, label, LHS, sign, RHS): + self.label = label.upper() # keep track of what variable we are using + self.LHS = LHS + self.sign = sign + self.RHS = RHS + # self.iterCount for iterable + + def sign(self) -> str: + return self.sign + + def RHS(self) -> float: + ''' + all numbera are subtypes of float i think>< + ''' + return self.RHS + + def __getitem__(self, key:str) -> float: + if self.LHS.get(key): + return self.LHS.get(key) + else: + return 0 + + def __setitem__(self, key:str, value: float): + self.LHS[key] = value + + def clear(self) -> Constraint: + return Constraint(self.label,dict(),self.sign,0) + + def add_S(self, sIndex: int) -> Constraint: + newLHS = self.LHS.copy() + if self.sign == ">": + newLHS[f"S{sIndex}"] = -1 + elif self.sign == "<": + newLHS[f"S{sIndex}"] = 1 + return Constraint(self.label,newLHS,'=',self.RHS) + + def invert(self) -> Constraint: + return Constraint( + self.label, + {k:-v for k,v in self.LHS.items()}, + self.invertSign[self.sign], + -self.RHS) + + def invertVar(self, variableName:str) -> Constraint: + newLHS = self.LHS.copy() + if newLHS.get(variableName): + newLHS[variableName] = -newLHS[variableName] + return Constraint( + self.label, + newLHS, + self.sign, + self.RHS) + + def addVar(self,newVars: List[str]) -> Constraint: + newLHS = self.LHS.copy() + for i in newVars: + newLHS[i] = 0 + return Constraint( + self.label, + newLHS, + self.sign, + self.RHS) + + def toDF(self,index=0)->DataFrame: + return DataFrame(self.LHS,index=[index]).join(DataFrame({"sign":self.sign,"RHS":self.RHS},index=[index])) + + def get_variables(self) -> List[str]: + return list(self.LHS.keys()) + + def __len__(self) -> int: + return len(self.LHS) + + def __repr__(self): + outString = "" + factors = sorted([k for k in self.LHS.keys() if self.label in k]) + S_factors = sorted([k for k in self.LHS.keys() if "S" in k]) + A_factors = sorted([k for k in self.LHS.keys() if "A" in k]) + + for k in factors+S_factors+A_factors: + v = self.LHS[k] + if v > 0: + outString += f" +{v}{k}" + elif v < 0: + outString += f" {v}{k}" + outString += f" {self.sign} {self.RHS}" + return outString[1:] \ No newline at end of file diff --git a/OR/Format.py b/OR/Format.py index 60d562a..08d912f 100644 --- a/OR/Format.py +++ b/OR/Format.py @@ -1,20 +1,33 @@ -try: - from ipython import display -except: - pass -def subscript(i): - # ima let you toggle between subscript and full height numbers - # return i - out = "" - for j in str(i): - out +=chr(ord(u'\u2080')+ int(j)) - return out -# LP Object -# duality finding method -def displayHelper(item): - #display works on jupyter but might not always be included? +def display2(item): try: display(item) except: print(item) - print() \ No newline at end of file + +class InfeasibleException(Exception): + pass + +class OptimalException(Exception): + pass + + +# try: +# from ipython import display +# except: +# pass +# def subscript(i): +# # ima let you toggle between subscript and full height numbers +# return i +# out = "" +# for j in str(i): +# out +=chr(ord(u'\u2080')+ int(j)) +# return out +# # LP Object +# # duality finding method +# def displayHelper(item): +# #display works on jupyter but might not always be included? +# try: +# display(item) +# except: +# print(item) +# print() \ No newline at end of file diff --git a/OR/Hungarian.py b/OR/Hungarian.py index 99729c9..1d5d5b1 100644 --- a/OR/Hungarian.py +++ b/OR/Hungarian.py @@ -1,76 +1,81 @@ +from .Format import display2 +from .LP import LP +from .Constraint import Constraint +from .ObjectiveFunction import ObjectiveFunction +from .NonNeg import NonNeg - -from .Format import displayHelper import numpy as np -from pandas import DataFrame +from pandas import DataFrame,Series from itertools import combinations +from typing import TypeVar,List,Tuple,Dict +Hungarian = TypeVar("Hungarian") + class Hungarian(): @classmethod - def new(cls,table): - df = DataFrame(table) - df.columns = [i+1 for i in range(df.shape[0])] - df.index = [i+1 for i in range(df.shape[0])] - return Hungarian(df) - - def __init__(self,costDF,reducedDF=False,assignedDF=False): - self.costDF = DataFrame(costDF).copy() - - if type(reducedDF) != bool: - self.reducedDF = DataFrame(reducedDF).copy() - else: - self.reducedDF = (DataFrame(costDF).copy() - .apply(lambda x:x-min(x),axis=1) - .apply(lambda x:x-min(x),axis=0)) + def new(cls,inArray:List[List[float]]) -> Hungarian: + ''' + convert list of lists to DataFrame for Hungarian + ''' + #data validity + if not all([len(i)-1 == len(inArray[-1]) for i in inArray[:-1]]): + raise ValueError("invalid input array") + elif not all([len(i) == len(inArray) for i in inArray[:-1]]): + raise ValueError("not square") - if type(assignedDF) != bool: - self.assignedDF = DataFrame(assignedDF).copy() - else: - self.assignedDF = False + aggCol = DataFrame(inArray[-1]).T + aggRow = DataFrame([i[-1] for i in inArray[:-1]]) + if not aggCol.applymap(lambda x: x in (0,1)).all().all(): + raise ValueError("invalid data in column subtotal, must be 1 or 0") + elif not aggRow.applymap(lambda x: x in (0,1)).all().all(): + raise ValueError("invalid data in row subtotals, must be 1 or 0") + costDF = DataFrame([i[:-1] for i in inArray[:-1]]) - def solve(self,echo=False): - if type(self.isOptimal()) != bool: - return self.isOptimal() - coveredDF = self.coverLines(echo=echo) - return self.reduceDF(coveredDF,echo=echo).solve(echo=echo) - + for i in (costDF,aggCol,aggRow): + i.columns = [j+1 for j in i.columns] + i.index = [j+1 for j in i.index] - def isOptimal(self): - zeros = () - # locate zeros - for colIndex in self.reducedDF.columns: - for rowIndex in self.reducedDF.index: - if self.reducedDF.loc[rowIndex,colIndex] == 0: - zeros += ((rowIndex,colIndex),) + return cls(costDF,aggCol,aggRow,dict()) - # check that there is a feasible combination - for zero in combinations(zeros,len(self.reducedDF.columns)): - x = set(i[0] for i in zero) - y = set(i[1] for i in zero) - if (len(x) == len(self.reducedDF.columns) - and len(y) == len(self.reducedDF.columns)): - assignedDF = DataFrame(np.zeros_like(self.costDF)) - assignedDF.columns = self.costDF.columns - assignedDF.index = self.costDF.columns - for coord in zero: - assignedDF.loc[coord[0],coord[1]] = 1 - return Hungarian(self.costDF,self.reducedDF,assignedDF) - return False - - def coverLines(self,echo=False): + def __init__(self,costDF: DataFrame,aggCol: DataFrame,aggRow: DataFrame,BFS: Dict): + self.costDF = costDF.copy() + self.aggCol = aggCol.copy() # aggregate of Columns + self.aggRow = aggRow.copy() # aggregate of Rows + self.BFS = BFS.copy() #dict key is coord tuple, value is assigned val + + def subtractMinimums(self,echo: bool=False) -> Hungarian: ''' - return a DataFrame where crossed out columns are represented by np.inf + for each row, subtract minimum + for each col, subtract minimum ''' - workingDF = self.reducedDF.copy() + nCostDF = (self.costDF.copy() + .apply(lambda x:x-min(x),axis=1) + .apply(lambda x:x-min(x),axis=0)) + return Hungarian(nCostDF,self.aggCol,self.aggRow,self.BFS) - #axis{index (0), columns (1)} - rows = (workingDF == 0).astype(int).sum(axis=1) - cols = (workingDF == 0).astype(int).sum(axis=0) - max0 = max(rows.append(cols)) + def isOptimal(self) -> bool: + ''' + if there is a feasible assignment with 0 entries of the zero entries, + this assignment is optimal + ''' + if self.BFS == dict(): + return False + else: + return True + def coverLines(self, echo:str=False) -> Tuple[List[int],List[int]]: + ''' + return a Tuple of (list of crossed rows, list of crossed cols) + ''' + workingDF = self.costDF.copy() + max0 = 1 while max0 >0: + rows = (workingDF == 0).sum(axis=1) + cols = (workingDF == 0).sum(axis=0) + max0 = max(rows.append(cols)) + if max0 in rows.array: for index, value in rows.items(): if value == max0: @@ -81,50 +86,104 @@ def coverLines(self,echo=False): if value == max0: workingDF.loc[:,index] = [np.inf] * workingDF.shape[0] break - rows = (workingDF == 0).astype(int).sum(axis=1) - cols = (workingDF == 0).astype(int).sum(axis=0) - max0 = max(rows.append(cols)) - # if echo: - # echoDisp = workingDF.copy() - # for rowIndex in echoDisp.index: - # for colIndex in echoDisp.columns: - # if echoDisp.loc[rowIndex,colIndex] - # idea to do - # - for row only - # | for col only - # + for row and col - - if echo:{displayHelper(workingDF.applymap(lambda x:"-" if np.isposinf(x) else x))} - return workingDF - - def reduceDF(self,coveredDF,echo=False): - rows = (coveredDF == np.inf).astype(int).sum(axis=1) - cols = (coveredDF == np.inf).astype(int).sum(axis=0) - maskedRows = [index for index, value in rows.items() if value == coveredDF.shape[1]] - maskedCols = [index for index, value in cols.items() if value == coveredDF.shape[1]] - minUncovered = coveredDF.min().min() + # if echo:{display2(workingDF.applymap(lambda x:"-" if np.isposinf(x) else x))} + rows = (workingDF == np.inf).all(axis=1) + rows = rows[rows==True].index.tolist() + cols = (workingDF == np.inf).all(axis=0) + cols = cols[cols==True].index.tolist() + return (rows,cols) - workingDF = self.reducedDF.copy() + def pivot(self, crossedLines:Tuple[List[int],List[int]],echo: bool=False): + if self.isOptimal(): + # return copy if already optimal + return self.copy() + + rows,cols = crossedLines + maskedDF = self.costDF.copy() + for row in rows: + maskedDF.loc[row,:] = np.inf + for col in cols: + maskedDF.loc[:,col] = np.inf + minVal = maskedDF.min().min() + if echo:{display2(maskedDF.applymap(lambda x: "-" if x==np.inf else x))} + workingDF = self.costDF.copy() for colIndex in workingDF.columns: for rowIndex in workingDF.index: - if colIndex in maskedCols and rowIndex in maskedRows: - workingDF.loc[rowIndex,colIndex] += minUncovered - elif colIndex not in maskedCols and rowIndex not in maskedRows: - workingDF.loc[rowIndex,colIndex] -= minUncovered - if echo:{displayHelper(workingDF)} - return Hungarian(self.costDF,workingDF) + if colIndex in cols and rowIndex in rows: + workingDF.loc[rowIndex,colIndex] += minVal + elif colIndex not in cols and rowIndex not in rows: + workingDF.loc[rowIndex,colIndex] -= minVal + if echo:{display2(workingDF)} + + base = Hungarian(workingDF,self.aggCol,self.aggRow,self.BFS) + zeros = base.costDF[base.costDF == base.costDF.min().min()].stack().index.tolist() + for zero in combinations(zeros,len(self.costDF)): + x = set(i[0] for i in zero) + y = set(i[1] for i in zero) + if (len(x) == len(self.costDF) + and len(y) == len(self.costDF)): + base.BFS = {(row,col):1 for row,col in zero} + return base + + def solve(self,echo:bool=False): + if self.isOptimal(): + return self + else: + return (self.subtractMinimums(echo=echo) + .pivot(self.coverLines(echo=echo),echo=echo) + .solve(echo=echo)) + + # getters + def to_LP(self,label:str="X"): + def varRepr(row,col,label): + return f"{label}_{row},{col}" + + # objective Function + allVariables = [varRepr(row,col,label) for row in self.costDF.index for col in self.costDF.columns] + objFunc = ObjectiveFunction.new("min", + *[self.costDF.loc[row,col] for row in self.costDF.index for col in self.costDF.columns], + labels=allVariables + ) + + rowConstraints = [Constraint.new( + *[self.costDF.loc[row,col] for row in self.costDF.index], + "<",self.aggCol.loc[1,col], + labels = [varRepr(row,col,label) for row in self.costDF.index] + ) for col in self.costDF.columns] + + colConstraints = [Constraint.new( + *[self.costDF.loc[row,col] for col in self.costDF.columns], + "<",self.aggRow.loc[row,1], + labels = [varRepr(row,col,label) for col in self.costDF.columns] + ) for row in self.costDF.index] + + nonNegConstraints = NonNeg.fromArray(*([">"]*len(allVariables)),labels=allVariables) + + return LP( + objFunc, + *rowConstraints, + *colConstraints, + *nonNegConstraints + ) + def objectiveValue(self): ''' return calculated objective value of this transport tableau using current assignments ''' - if type(self.assignedDF) != bool: - return self.costDF.mul(self.assignedDF).sum().sum() - else: - raise ValueError("assignments have not been calculated") + return sum([value * self.costDF.loc[row,col] for (row,col),value in self.BFS.items()]) + + def bfsDF(self) -> DataFrame: + baseDF = self.costDF.copy().applymap(lambda x:"") + for (row,col),value in self.BFS.items(): + baseDF.loc[row,col] = value + return baseDF + + def copy(self): + return Hungarian(self.costDF,self.aggCol,self.aggRow,self.BFS) def __repr__(self): - if type(self.assignedDF) == bool: - return f"Hungarian\n{str(self.costDF)}\n" - else: - return f"Assignments\n{str(self.assignedDF)}\n" \ No newline at end of file + baseDF = self.costDF.join(Series(self.aggRow.loc[:,1],name="row")).append(Series(self.aggCol.loc[1,:],name="col")).fillna("") + b = [f"X_{row},{col} = {c}, " for (row,col),c in self.BFS.items()] + b = "".join(b)[:-1] + return f"{baseDF}\nBFS :\n{self.bfsDF()}" diff --git a/OR/LP.py b/OR/LP.py index 3303799..e944352 100644 --- a/OR/LP.py +++ b/OR/LP.py @@ -1,155 +1,96 @@ -from pandas import DataFrame -from pandas import Series -import numpy as np -from .Format import * +from .ObjectiveFunction import ObjectiveFunction +from .Constraint import Constraint +from .NonNeg import NonNeg -class Constraint(): - # unused - def __init__(self): - self.variables = variables # 1d dataFrame - self.sign = sign # String - self.RHS = RHS # number - - @classmethod - def fromArray(cls,variables,sign,RHS,names=False): - if not names: - names = [f"X{subscript(i+1)}"for i in range(variables)] - variables = DataFrame(variables,columns=names) - return Constraint(variables,sign,RHS) +# LP is +# a objective Function +# an objective funtion type +# an array of constraints class LP(): - @classmethod # raw constructor - def new(cls,objectiveFn,constraints,signs,factorNames = False): - #throw in some error checking to ensure valid LP - objectiveType = objectiveFn[0] - if not factorNames: - factorNames = [f"X{subscript(i+1)}" for i in range(len(objectiveFn)-1)] - factorNames = factorNames + ["signs","RHS"] - body = DataFrame( - [objectiveFn[1:]+["",""]]+ - constraints, - columns=factorNames) - return LP(objectiveType,body,signs) + dualConvert = { + "min":{ + "newObj":"max", + "c":{">":">","=":"URS","<":"<"}, + "v":{">":"<","URS":"=","<":">"} + },"max":{ + "newObj":"min", + "c":{"<":">","=":"URS",">":"<"}, + "v":{">":">","URS":"=","<":"<"} + } + } + def __init__(self, objectiveFunction, *constraints): + if not isinstance(objectiveFunction, ObjectiveFunction): + raise ValueError("invalid ObjectiveFunction") + for i, c in enumerate(constraints): + if not isinstance(c, Constraint): + raise ValueError(f"invalid constraint {i} of {c}") + self.objFunc = objectiveFunction + self.constraints = constraints - def __init__(self,objectiveType,body,signs): - self.objectiveType = objectiveType # String - self.body = body # DataFrame - self.signs = signs # List + def standardForm(self): + # add slack surplus variables + # invert constraint row if RHS < 0 + # invert variable row if x1 <= 0 is a constraint + newObjective = self.objFunc + newConstraints = list(self.constraints) - def getDual(self,factorPrefix = "Y",factorNames=[]): - #factorNames is some List /nparray? + # non-negative RHS + for i, c in enumerate(self.constraints): + if c.RHS < 0: + newConstraints[i] = c.invert() - if self.objectiveType == "min": - constraint_sign = {">":">","=":"urs","<":"<"} - variable_sign = {">":"<","urs":"=","<":">"} - else: - constraint_sign = {"<":">","=":"urs",">":"<"} - variable_sign = {">":">","urs":"=","<":"<"} - objectiveType = "min"if self.objectiveType == "max" else "max" - - body = self.body.copy() - signs = self.signs.copy() + # non-neg constraints >= 0 + for i, c in enumerate(self.constraints): + if isinstance(c, NonNeg) and c.sign == "<": + # invert col + newObjective = newObjective.invertVar(c.name) + for j, d in enumerate(newConstraints): + newConstraints[j] = newConstraints[j].invertVar(c.name) + # invert row + newConstraints[i] = newConstraints[i].invert() + return LP(newObjective, *newConstraints) + + def add_S(self): + newConstraints = [c.add_S(i+1) for i,c in enumerate(self.constraints)] + return LP(self.objFunc,*newConstraints) + + def dual(self,label="Y",labels = None): + converter = self.dualConvert[self.objFunc.getObjType()] + newObjective = ObjectiveFunction.new( + converter['newObj'], + *[i.RHS for i in self.constraints if not isinstance(i,NonNeg)], + label='Y',labels=labels) + # newObjective + newVars = set() + for c in self.constraints: + newVars.update(c.LHS.keys()) + newVars = sorted(newVars) + newConstraints = [] + for var in newVars: + newConstraints.append(Constraint.new( + *[c[var] for c in self.constraints if not isinstance(c,NonNeg)], + converter['v'][self.nonNeg(var)],self.objFunc[var], + label=label,labels=labels + )) + newConstraints+= NonNeg.fromArray( + *[converter['c'][i.sign] for i in self.constraints if not isinstance(i,NonNeg)], + label=label,labels=labels + ) + return LP(newObjective,*newConstraints) - obj = body.iloc[:,-1:].T.drop(0,axis=1) - mid = body.iloc[:,:-2].drop(0).T - obj = obj.append(mid).reset_index().drop("index",axis=1) - - if len(factorNames)== len(obj.columns): - obj.columns = factorNames + def nonNeg(self,variable): + temp = [c for c in self.constraints if isinstance(c,NonNeg) and c.name == variable] + if len(temp) == 0: + return "URS" + elif len(temp) > 1: + raise ValueError(f"multiple non-negative constraints for {variable} found") else: - obj.columns = [f"{factorPrefix}{subscript(i)}" for i in obj.columns] - - obj.insert(obj.shape[1],"signs",np.array([""] + [variable_sign[i] for i in signs])) - obj.insert(obj.shape[1],"RHS",np.append([""],body.iloc[0,:-2])) - - newSigns = list(body.loc[:,"signs"])[1:] - newSigns = [constraint_sign[i] for i in newSigns] - return LP(objectiveType,obj,newSigns) - - def getStandardForm(self): - #set positive signs - posSigns = self.posSigns() - #modify urs factors - noURS = posSigns.clearURS() - #add surplus/slack variables - SFactors = noURS.addSFactors() - - #set positive RHS - return SFactors.posRHS() - def posSigns(self): - '''Set the X<0 factors to X>0''' - objectiveType = self.objectiveType - body = self.body.copy() - signs = self.signs.copy() - for i in range(len(signs)): - if signs[i] == "<": - signs[i] = ">" - body.iloc[:,i] = body.iloc[:,i].apply(lambda x: x*-1) - return LP(objectiveType,body,signs) - def clearURS(self): - objectiveType = self.objectiveType - body = self.body.copy() - signs = self.signs.copy() - i = 0 - while i < len(signs): - if signs[i] == "urs": - signs[i] = ">" - signs.insert(i,">") - body.insert(i+1,body.columns[i]+"'",body.iloc[:,i]) - body.insert(i+2,body.columns[i]+"\"",body.iloc[:,i]*-1) - body = body.drop(body.columns[i],axis=1) - i+=1 - return LP(objectiveType,body,signs) - def addSFactors(self): - objectiveType = self.objectiveType - body = self.body.copy() - signs = self.signs.copy() - for i in range(1,body.shape[0]): - if body.iloc[i].loc["signs"] == "<": - newCol = [0] * body.shape[0] - newCol[i] = 1 - body.insert(body.shape[1]-2,f"s{subscript(i)}",np.array(newCol)) - signs.insert(body.shape[1]-2,">") - elif body.iloc[i].loc["signs"] == ">": - newCol = [0] * body.shape[0] - newCol[i] = -1 - body.insert(body.shape[1]-2,f"s{subscript(i)}",np.array(newCol)) - signs.insert(body.shape[1]-2,">") - body.iloc[i] = body.iloc[i].apply(lambda x: "=" if x in [">","<"] else x) - return LP(objectiveType,body,signs) - def posRHS(self): - objectiveType = self.objectiveType - body = self.body.copy() - signs = self.signs.copy() - for i in range(1,body.shape[0]): - row = body.iloc[i] - if row.loc["RHS"]<=0: - body.iloc[i] = body.iloc[i].apply(lambda x: x*-1 if type(x)!=str else x) - return LP(objectiveType,body,signs) - - def __eq__(self,other): - if all(self.body.eq (other.body)): - return True - else: - return False - - def display(self): - outDF = self.body.copy() - outDF.insert(0,"",[1]+[""]*(len(self.body)-1)) - outDF.iloc[0,0] = self.objectiveType - for i in outDF.columns[1:-2]: - outDF.loc[:,i] = (outDF.loc[:,i].apply(lambda x: "" if x==0 else x)) - outDF.loc[:,i] = (outDF.loc[:,i].apply(str)) - outDF.loc[:,i] = (outDF.loc[:,i].apply(lambda x: x if (len(x)==0 or x[0]=="-") else "+"+x)) - outDF.loc[:,i] = (outDF.loc[:,i].apply(lambda x: x if len(x)==0 else x+i)) - - signsDF = DataFrame([[""] + self.signs+["",""]],columns = outDF.columns) - for i in signsDF.columns[1:-2]: - signsDF.loc[:,i] = (signsDF.loc[:,i].apply(lambda x: i + x + " 0" if x in [">","<"] else x)) - signsDF.loc[:,i] = (signsDF.loc[:,i].apply(lambda x: i + " " + x if x =="urs" else x)) - outDF = outDF.append(signsDF) - outDF = outDF.reset_index().drop(["index"],axis=1) - return outDF + return temp[0].sign + def __repr__(self): - return str(self.display()) + outString = str(self.objFunc) + for i in self.constraints: + outString += "\n" + str(i) + return outString diff --git a/OR/NonNeg.py b/OR/NonNeg.py new file mode 100644 index 0000000..e49d138 --- /dev/null +++ b/OR/NonNeg.py @@ -0,0 +1,31 @@ +from .Constraint import Constraint + +class NonNeg(Constraint): + @classmethod + def fromArray(cls, *args, label="X", labels=None): + outList = [] + if labels == None: + labels = [f"{label}{i+1}" for i in range(len(args))] + for i, v in enumerate(args): + if v.upper() != "URS": + outList += [cls(labels[i], v)] + return outList + + def __init__(self, variable, sign, RHS=0): + self.name = variable + self.label = variable[0] + self.LHS = {variable: 1} + self.sign = sign + self.RHS = RHS + + def invertVar(self, variableName): + if variableName == self.name: + return NonNeg(self.name,self.sign, RHS=self.RHS) + else: + return self + + def add_S(self,index): + return self + + def invert(self): + return NonNeg(self.name, self.invertSign[self.sign], RHS=-self.RHS) diff --git a/OR/ObjectiveFunction.py b/OR/ObjectiveFunction.py new file mode 100644 index 0000000..0ca5dce --- /dev/null +++ b/OR/ObjectiveFunction.py @@ -0,0 +1,75 @@ +from pandas import DataFrame + +from .Constraint import Constraint + +from typing import List,Tuple,TypeVar +ObjectiveFunction = TypeVar("ObjectiveFunction") + + +class ObjectiveFunction(Constraint): + @classmethod + def new(cls, objType:str, *args,label:str="X", labels: List[str]=None) -> ObjectiveFunction: + if labels == None: + # labels overrides label + if label.upper() in ["A", "S"]: + raise ValueError(f"A and S are reserved variables") + labels = [f"{label.upper()}{i+1}" for i in range(len(args))] + elif len(labels) != len(args): + raise ValueError( + f"provided {len(labels)} labels for {len(args)} variables") + else: + label = labels[0][0] + + if objType not in ["min", "max"]: + raise ValueError(f"{objType} is not a valid objective") + + LHS = {k: v for k, v in zip(labels, args)} + sign = objType + RHS = None + return cls(label, LHS, objType, RHS) + + def __init__(self, label, LHS, objType, RHS): + super().__init__(label, LHS, objType, RHS) + + def invertVar(self, variableName:str) -> ObjectiveFunction: + newLHS = self.LHS.copy() + if newLHS.get(variableName): + newLHS[variableName] = -newLHS[variableName] + return ObjectiveFunction(self.label, newLHS, self.sign, self.RHS) + + def getObjType(self) -> str: + return self.sign + + def toDF(self,target:str="Z") -> DataFrame: + if self.getObjType() == "max": + return ( + DataFrame({f"{target}":1},index=[0]) + .join(-1*DataFrame(self.LHS,index=[0]) + .join(DataFrame({"sign":"=","RHS":0},index=[0]))) + ) + + else: # objType == "min" + return ( + DataFrame({f"-{target}":1},index=[0]) + .join(DataFrame(self.LHS,index=[0]) + .join(DataFrame({"sign":"=","RHS":0},index=[0]))) + ) + + return target.join(-1*DataFrame(self.LHS,index=[0]).join(DataFrame({"sign":"=","RHS":0},index=[0]))) + + def clear(self) -> ObjectiveFunction: + return ObjectiveFunction(self.label,dict(),self.sign,None) + + def __repr__(self) -> str: + outString = f"{self.sign}" + factors = sorted([k for k in self.LHS.keys() if self.label in k]) + S_factors = sorted([k for k in self.LHS.keys() if "S" in k]) + A_factors = sorted([k for k in self.LHS.keys() if "A" in k]) + + for k in factors+S_factors+A_factors: + v = self.LHS[k] + if v > 0: + outString += f" +{v}{k}" + elif v < 0: + outString += f" {v}{k}" + return outString diff --git a/OR/TableauBase.py b/OR/TableauBase.py index e33c0c9..59d63eb 100644 --- a/OR/TableauBase.py +++ b/OR/TableauBase.py @@ -1,58 +1,63 @@ from pandas import DataFrame import pandas as pd import numpy as np -from .Format import * -''' -keep it nice and simple -just solve a tableau -until row 0 has no more negative numbers -or we face some other issue -''' +from .Format import OptimalException, display2 +from .LP import LP +from .NonNeg import NonNeg -def round5(x): - # default dispFunc - return round(x,5) +from typing import List,Tuple,Dict,TypeVar +TableauBase = TypeVar("TableauBase") class TableauBase(): - def __init__(self,dataFrame,basicVar,dispFunc): + @classmethod + def fromLP(cls,inLP:LP) -> TableauBase: + inLP = inLP.standardForm() + constraints = [(i+1,c.add_S(i+1)) for i,c in enumerate(inLP.constraints) if not isinstance(c,NonNeg)] + target = inLP.objFunc.toDF().columns[0] + label = inLP.objFunc.label + + uniqueVars = set() + uniqueVars.update(inLP.objFunc.get_variables()) + for i,c in constraints: + uniqueVars.update(c.get_variables()) + + baseDF = pd.concat( [inLP.objFunc.toDF()] + [c.toDF(i) for i,c in constraints] ,sort=False) + baseDF = baseDF.fillna(0).drop(columns="sign") + #add Basic var Column + + basicVar = [""] + for i in range(len(constraints)): + if f"A{i+1}" in uniqueVars: + basicVar += [f"A{i+1}"] + elif f"S{i+1}" in uniqueVars: + basicVar += [f"S{i+1}"] + else: + basicVar += [f"{label}{i+1}"] + basicVar = pd.DataFrame( basicVar,columns=["Basic"]) + + cols = ([target]+# if inLP.objFunc.sign == "max" else f"-{target}"]+ + sorted([i for i in uniqueVars if label in i])+ + sorted([i for i in uniqueVars if "S" in i])+ + sorted([i for i in uniqueVars if "A" in i])+ + ["RHS"]) + + return TableauBase(baseDF[cols],basicVar) + + def __init__(self,dataFrame: DataFrame,basicVar: DataFrame): self.dataFrame = dataFrame - self.basicVar = tuple(basicVar) - self.dispFunc = dispFunc - - def setBFS(self,bfs): - # takes in BFS in sequence down the tableau - out = self - for i in range(len(bfs)): - out = out.pivot(bfs[i],i+1) - return out + self.basicVar = basicVar - def canonical(self,echo=False): - if echo: - print("pivoting on all basic Vars") - for var in self.basicVar[1:]: - self = self.pivot(var,self.pickPivotRow(var,echo=echo)) - if echo:{print(f"making {var} canonical\n")} - if echo: - self.display() - print("In canonical form\n") - return self - - #pivot action - def autoPivot(self,echo=False): - ''' - stack a search and a pivot together, single cycle - ''' - #echo is to toggle printing of steps - # set other non-basic factors to 0 - # solve for row maximum value of interested factor where basic variables are 0 to inf - t = self.dataFrame - factor,row = self.pickPivot(echo) - if row == None: - print("unbounded lp") - return self - return self.pivot(factor,row) - def pivot(self, factor, row): + def canonical(self,echo:bool = False) -> TableauBase: + return self.setBFS(list(self.basicVar.loc[:,"Basic"])[1:],echo=echo) + + def asDF(self) -> DataFrame: + return self.dataFrame.join(self.basicVar) + + def __repr__(self) -> str: + return str(self.asDF()) + + def pivot(self, factor: str, row: int,echo: bool=False) -> TableauBase: ''' core pivot action using matrix math ''' @@ -65,36 +70,51 @@ def pivot(self, factor, row): result = self.dataFrame - np.dot(DataFrame(pivot_col),DataFrame(enter_row).T) result.iloc[row] = enter_row - newVar = self.basicVar[:row]+(factor,)+self.basicVar[row+1:] - return TableauBase(result,newVar,self.dispFunc) - def pickPivot(self,echo=False): + newVar = self.basicVar.copy() + newVar.loc[row] = factor + + out = TableauBase(result,newVar) + if echo:{display2(out)} + return out + + def getBFS(self) -> Dict[str,float]: + return {k:v for k,v in zip(self.basicVar.loc[1:,"Basic"],self.dataFrame.loc[1:,"RHS"])} + + def getObjectiveVal(self) -> Dict[str,float]: + if self.dataFrame.columns[0][0] == "-": + return {self.dataFrame.columns[0][1:] : -1* self.dataFrame.loc[0,"RHS"]} + else: + return {self.dataFrame.columns[0] : self.dataFrame.loc[0,"RHS"]} + + def setBFS(self,bfs:List[str],echo=False) -> TableauBase: ''' - combining pickRow and pickCol + pivot tableau to use new BFS ''' - t = self.dataFrame - factor = self.pickPivotCol(echo) - # if factor == None: # unbounded - # return None,None - row = self.pickPivotRow(factor,echo) - if row == None: # unbounded - return None, None - if echo: - print(f"pivoting on {factor} and row {row}\n") - return factor,row - def pickPivotCol(self,echo=False): + k = self + if echo:{display2(k)} + for i,variable in enumerate(bfs): + k = k.pivot(variable,i+1,echo=echo) + return k + + def isOptimal(self) -> bool: + if all(self.dataFrame.loc[0].drop("RHS") >= 0): + return True + else: + return False + + def pickPivotCol(self,echo=False) -> str: ''' get the most negative column ''' - t = self.dataFrame - t2 = t.iloc[:,:-1].iloc[0] # objective function - minT = min(t2) # value of most negative factor - factor = None - for i in t.columns: - if t[i][0] == minT: # get entering variable with the lowest index - factor = i # aka closest to the left of tableau - break # Bland's rule - return factor - def pickPivotRow(self,factor,echo=False): + minVal = self.dataFrame.loc[0].iloc[1:-1].idxmin() + if echo:{display2(DataFrame(self.dataFrame.loc[0]).T)} + if self.dataFrame.loc[0,minVal] < 0 : + return minVal + else: + raise OptimalException("row 0 all positive, Tableau is optimal") + # return None + + def pickPivotRow(self,factor: str,echo=False) -> int: ''' ratio test 1) solution divide entering col must be positive @@ -102,139 +122,46 @@ def pickPivotRow(self,factor,echo=False): 3) take smallest of all remaining positive non inf 4) --> such that pivoting on this will cause next solution set to remain positive no. ''' - t = self.dataFrame - # ratio = t.iloc[1:,-1]/t.loc[1:,factor] - top = t.iloc[1:,-1] - bot = t.loc[1:,factor] - ratio = top.copy() - - for i in range(len(ratio)): - if bot.iloc[i] == 0: - ratio.iloc[i] = np.inf - else: - ratio.iloc[i] = top.iloc[i]/bot.iloc[i] - - for i in range(len(ratio)): - if ratio.iloc[i] < 0: # replace if opposite sign --> values can go to inf - ratio.iloc[i] = np.inf # replace -inf from x/-0.0 - elif ratio.iloc[i] == -np.inf: - ratio.iloc[i] = np.inf - row = pd.to_numeric(ratio).idxmin() #returns first index of smallest value, follows bland's rule + col = self.dataFrame.loc[:,factor] + RHS = self.dataFrame.loc[:,"RHS"] - enterCol = np.array([""]*len(self.basicVar)) - if ratio[row] == np.inf: - row = None - else: - enterCol[row] = "*" - if echo: - t = t.applymap(self.dispFunc) - displayHelper( - t.assign(basicVar=self.basicVar).join(ratio.to_frame(name="ratioTest")).assign(Entering=enterCol) - ) - return row - - @classmethod #constructor from LP - def fromLP(cls,inLP): - LP2 = inLP.getStandardForm() - body = LP2.body - body = body.drop("signs",axis=1) - if LP2.objectiveType == "min": - target = "-Z" + l = RHS/col + l = l[1:][l > 0][l != np.nan] + if echo:{display2(self.asDF().join(DataFrame(l,columns=["Ratio Test"])))} + if len(l) >= 1: + return l.idxmin() else: - target = "Z" - body.iloc[0] = body.iloc[0]*-1 - body.insert(0,target,[1]+[0]*(body.shape[0]-1)) - body.iloc[0,-1] = 0 - - basicVars = ["Z"] + [f"s{subscript(i)}" for i in range(1,body.shape[0])] - return TableauBase(body,basicVars,round5) - - #displays/getters - def display(self): - return displayHelper(self.displayfmt()) - def displayfmt(self): #default round to 5 decimal places for floating point sanity - out = self.dataFrame - out = out.applymap(self.dispFunc) - df = self.dataFrame.assign(Basic = self.basicVar) - return out.assign(Basic = self.basicVar).loc[:, (df != 0).any(axis=0)] - def __repr__(self): - return str(self.displayfmt()) - def getBFS(self): - allVar = self.dataFrame.columns - bfs = "{" - for var in allVar[:-1]: - value = 0 - if var in self.basicVar: - value = self.dataFrame.iloc[self.basicVar.index(var),-1] - value = self.dispFunc(value) - - if var[0] =="-": - var = var[1:] - value*=-1 - - - - bfs += f" {var} = {value}," - bfs = bfs[:-1] - bfs += " }" - return bfs - #comparisons - def equals(self,other): - diff = self.dataFrame-other.dataFrame - return diff.apply(lambda x: abs(x)<1e-15).all().all() - - def solve(self,echo=False,past=[]): - if self in past: - #looped - if echo: - self.display() - print("Solved") - return self - elif all( [i>=0 for i in self.dataFrame.iloc[0,:-1]] ): - # no more negative row 0 factors - if echo: - self.display() - print("Solved") - return self - else: - past +=[self] - t1 = self.autoPivot(echo=echo) - return t1.solve(echo=echo,past=past) + raise OptimalException("no valid pivot row found, Tableau is optimal") - # graveyard - @classmethod - def genericTableau(cls,objective,constraints,factorNames=None): - factorNo = len(objective) - 1 - if not all([len(i)-2 == factorNo for i in constraints]): - raise ValueError("inconsistent number of factors") - elif objective[0] not in ["max","min"]: - raise ValueError("invalid objective function min or max") - elif not all([i[-2] in ["<","=",">"] for i in constraints]): - raise ValueError("constraints must have =, < or >") - elif factorNames: #validate col if it exists - if not( - (type(factorNames)==list) and - (len(factorNames)==factorNo) and - (all([type(i)==str for i in factorNames])) - ): - raise ValueError("invalid factorNames") - - if factorNames == None: - factorNames = [f"X{subscript(i+1)}" for i in range(len(factorNo))] - - #basic factors - factors = np.array([np.array(objective[1:])]+[i[:-2] for i in constraints]) - if objective[0] == "max": - factors[0] *=-1 - #non-basic variables - nonBasic = np.array([[0] + [0 for i in constraints]]).T - solution = np.array([[0] + [i[-1] for i in constraints]]).T - if objective[0] == "min": - target = "-z" - else: - target = "z" - - header = [target] + factorNames + ["Solution"] - tabData = np.hstack((nonBasic,factors,solution)) - basicVars = (target,)+tuple(factorNames) - return header,tabData,basicVars + def pickPivot(self,echo=False) -> Tuple[str,str]: + ''' + return (row,column) + combining pickRow and pickCol + ''' + try: + col = self.pickPivotCol(echo=echo) + row = self.pickPivotRow(col,echo=echo) + return col,row + except Exception as e: + raise e + + def autoPivot(self,echo=False) -> TableauBase: + ''' + stack a search and a pivot together, single cycle + ''' + return self.pivot(*self.pickPivot(echo=echo),echo=echo) + + def solve(self,echo=False,past=[]) -> TableauBase: + ''' + loop autosolve until result produced + or error encountered + ''' + k = self + try: + while True: + k = k.autoPivot(echo=echo) + if echo:{display2(k.asDF())} + except Exception as e: + if echo:{print(e)} + finally: + return k diff --git a/OR/TableauNormal.py b/OR/TableauNormal.py deleted file mode 100644 index f5f55ab..0000000 --- a/OR/TableauNormal.py +++ /dev/null @@ -1,10 +0,0 @@ -from .Tableau import Tableau -import numpy as np - -''' -Not really sure what to do with this one -''' - -class TableauNormal(Tableau): - def __init__(self): - pass \ No newline at end of file diff --git a/OR/TransportTable.py b/OR/TransportTable.py index fa8445c..f5893e0 100644 --- a/OR/TransportTable.py +++ b/OR/TransportTable.py @@ -1,427 +1,380 @@ -from .Format import * -from .LP import LP import numpy as np from pandas import DataFrame,Series +import pandas as pd + +from .LP import LP +from .ObjectiveFunction import ObjectiveFunction +from .Constraint import Constraint +from .NonNeg import NonNeg +from .Format import Optimal, OptimalException + +from typing import List,Tuple,Dict,TypeVar +TransportTable = TypeVar("TransportTable") class TransportTable(): + @classmethod - #factory to produce BFS-less tptTable - def new(cls,costs,RHS,BHS,RHSName = "RHS",BHSName = "BHS"): - costDF = DataFrame(costs).copy() - costDF.columns = [i+1 for i in range(len(costDF.columns))] - costDF.index = [i+1 for i in range(len(costDF.index))] - - RHS = Series(RHS,name=RHSName).copy() - RHS.index = [i+1 for i in range(len(RHS.index))] - BHS = Series(BHS,name=BHSName).copy() - BHS.index = [i+1 for i in range(len(BHS.index))] - BFS = dict() - return TransportTable(costDF,BFS,RHS,BHS) + def new(cls,inArray:List[List[float]]) -> TransportTable: + ''' + convert list of lists to DataFrame for TransportTable + ''' + #data validity + if not all([len(i)-1 == len(inArray[-1]) for i in inArray[:-1]]): + raise ValueError("invalid input array") + + aggCol = DataFrame(inArray[-1]).T + aggRow = DataFrame([i[-1] for i in inArray[:-1]]) + costDF = DataFrame([i[:-1] for i in inArray[:-1]]) - def __init__(self,costDF,BFS,RHS,BHS): + for i in (costDF,aggCol,aggRow): + i.columns = [j+1 for j in i.columns] + i.index = [j+1 for j in i.index] + + return cls(costDF,aggCol,aggRow,dict()) + + def __init__(self,costDF: DataFrame,aggCol: DataFrame,aggRow: DataFrame,BFS: Dict): self.costDF = costDF.copy() - self.BFS = BFS #dict key is coord tuple, value is assigned val - self.BHS = BHS.copy() - self.RHS = RHS.copy() - - # table stuff - def assign(self,rowIndex,colIndex,amount): # TODO - if (amount <= self.BHS.loc[colIndex] # amount + col sum - and amount <= self.RHS.loc[rowIndex]): # amount + row sum - - newRHS = self.RHS.copy() - newRHS.loc[rowIndex] -= amount - - newBHS = self.BHS.copy() - newBHS.loc[colIndex] -= amount - - BFS = self.BFS.copy() - BFS[(rowIndex,colIndex)] = amount - - return TransportTable(self.costDF.copy(),BFS, newRHS,newBHS) - else: - raise ValueError("too large of an assignment") - def assignMax(self,rowIndex,colIndex): - amount = min(self.RHS.loc[rowIndex],self.BHS.loc[colIndex]) - return self.assign(rowIndex,colIndex,amount) - - # LP stuff - def to_LP(self): - variablesNo = self.costDF.shape[0]*self.costDF.shape[1] - #objectiveFn - objFn = ["min"] + list(self.costDF.to_numpy().flatten()) - signs = [">"] * variablesNo - factorNames = [] - constraints = [] - - #row constraints - for rowIndex in self.costDF.index: - template = [0]*variablesNo + ["=",self.RHS.loc[rowIndex]] - for colIndex in self.costDF.columns: - factorNames +=[f"X{subscript(rowIndex)},{subscript(colIndex)}"] - template[(rowIndex-1)*self.costDF.shape[1] + colIndex-1] = 1 - constraints.append(template) - - #col constraints - for colIndex in self.costDF.columns: - template = [0]*variablesNo + ["=",self.BHS.loc[colIndex]] - for rowIndex in self.costDF.index: - template[(rowIndex-1)*self.costDF.shape[1] + colIndex-1] = 1 - constraints.append(template) - return LP.new(objFn,constraints,signs,factorNames=factorNames) - - def to_DualLP(self): - baseLP = self.to_LP() - factorNames = [f"u{i}" for i in self.costDF.index] + [f"v{i}" for i in self.costDF.columns] - return baseLP.getDual(factorNames=factorNames) - - #BFS - - def northWest(self,echo=False,origBHS=False,origRHS=False,curDir="col",checking=(1,1)): - #checking is tuple of coordinates - if echo == True: - self.display() - if self.RHS.sum() + self.BHS.sum() == 0: - #all assignments done - return TransportTable(self.costDF.copy(),self.BFS.copy(),origRHS,origBHS) - if type(origBHS) == bool: - # only copy on first loop - # otherwise pass through orig RHS and BHS - origBHS = self.BHS.copy() - origRHS = self.RHS.copy() + self.aggCol = aggCol.copy() # aggregate of Columns + self.aggRow = aggRow.copy() # aggregate of Rows + self.BFS = BFS.copy() #dict key is coord tuple, value is assigned val + + # LP stuff + def to_LP(self,label: str = "X") -> LP: + ''' + return the LP representing this tableau + with factors as X_row,col + ''' + def varRepr(row,col,label): + return f"{label}_{row},{col}" + + # objective Function + allVariables = [varRepr(row,col,label) for row in self.costDF.index for col in self.costDF.columns] + objFunc = ObjectiveFunction.new("min", + *[self.costDF.loc[row,col] for row in self.costDF.index for col in self.costDF.columns], + labels=allVariables + ) + + rowConstraints = [Constraint.new( + *[self.costDF.loc[row,col] for row in self.costDF.index], + "<",self.aggCol.loc[1,col], + labels = [varRepr(row,col,label) for row in self.costDF.index] + ) for col in self.costDF.columns] - if curDir =="col": - nextState = self.assignMax(checking[0],checking[1]) - rowSum = origRHS.loc[checking[0]] - for k,v in nextState.BFS.items(): - if k[0] == checking[0]: - rowSum -= v - if rowSum == 0: - # row Cleared, assign next row - return nextState.northWest(echo=echo,origBHS=origBHS,origRHS=origRHS,curDir="row",checking=(checking[0]+1,checking[1])) - else: - # col uncleared assign next col - return nextState.northWest(echo=echo,origBHS=origBHS,origRHS=origRHS,curDir="col",checking=(checking[0],checking[1]+1)) - - else: # curDir = "row" - # check (checking[0]+1,checking[1]) - nextState = self.assignMax(checking[0],checking[1]) - colSum = origBHS.loc[checking[1]] - for k,v in nextState.BFS.items(): - if k[1] == checking[1]: - colSum -= v - if colSum == 0: - # col Cleared, assign next col - return nextState.northWest(echo=echo,origBHS=origBHS,origRHS=origRHS,curDir="col",checking=(checking[0],checking[1]+1)) - else: - # col uncleared assign next row - return nextState.northWest(echo=echo,origBHS=origBHS,origRHS=origRHS,curDir="row",checking=(checking[0]+1,checking[1])) + colConstraints = [Constraint.new( + *[self.costDF.loc[row,col] for col in self.costDF.columns], + "<",self.aggRow.loc[row,1], + labels = [varRepr(row,col,label) for col in self.costDF.columns] + ) for row in self.costDF.index] + + nonNegConstraints = NonNeg.fromArray(*([">"]*len(allVariables)),labels=allVariables) + + return LP( + objFunc, + *rowConstraints, + *colConstraints, + *nonNegConstraints + ) - def minimumCost(self,echo=False,prevCost=0,origBHS=False,origRHS=False,crossedOut=[[],[]]): - if echo == True: - self.display() - if self.RHS.sum() + self.BHS.sum() == 0: - #all assignment done - if len(self.BFS) + 1 == len(self.RHS) + len(self.BHS): - # if enough BFS found - return TransportTable(self.costDF.copy(),self.BFS.copy(),origRHS,origBHS) - if type(origBHS) == bool: - # only copy on first loop - # otherwise pass through orig RHS and BHS - origBHS = self.BHS.copy() - origRHS = self.RHS.copy() - - costs = np.sort(np.unique(self.costDF.to_numpy().flatten())) - filter_arr = costs > prevCost-1 - costs = costs[filter_arr] - # start searching from equal to previous cost - for cost in costs: - for rowIndex in self.costDF.index: - if rowIndex in crossedOut[0]: - continue - for colIndex in self.costDF.columns: - if colIndex in crossedOut[1]: - continue - if(self.costDF.loc[rowIndex,colIndex] == cost): - nextState = self.assignMax(rowIndex,colIndex) - - rowSum = origRHS.loc[rowIndex] - colSum = origBHS.loc[colIndex] - for k,v in nextState.BFS.items(): - if k[0] == rowIndex: - rowSum -= v - if k[1] == colIndex: - colSum -= v - #only cross out one row/col - if (rowSum == 0 - and len(self.RHS) - len(crossedOut[0]) > 1): - crossedOut[0].append(rowIndex) - elif (colSum == 0 # TODO testing if last case assigned is value 0 works - and len(self.BHS)-len(crossedOut[1]) > 1): # len-len is supposed to not eliminate the last in each dir - crossedOut[1].append(colIndex) - return (nextState.minimumCost(echo=echo,prevCost=cost,origBHS=origBHS,origRHS=origRHS,crossedOut=crossedOut)) + #BFS stuff + def northWest(self,echo: bool=False) -> TransportTable: + ''' + assign a BFS using northWest Method + ''' + def NW(row:int,col:int)-> Dict: + ''' + private recursive function to do allocations + ''' + if echo: + df = (assignDF.join(Series(rowAgg, name="remain")) + .append(Series(colAgg, name="remain")).fillna("")) + try: + display(df) + except: + print(df) + + # BFS complete + if len(BFS) == len(colAgg) + len(rowAgg) - 1: + return None + + # allocate to current checking + minVal = min(colAgg.loc[col],rowAgg.loc[row]) + assignDF.loc[row,col] = minVal + BFS[(row,col)] = minVal + colAgg.loc[col] -= minVal + rowAgg.loc[row] -= minVal + + # go to next + if colAgg.loc[col] == 0 and col < len(assignDF.columns): + colAgg.loc[col] = "X" + return NW(row,col+1 ) + else: # rowAgg.loc[row] == 0 and row < len(assignDF.index) + rowAgg.loc[row] = "X" + return NW(row+1,col) + + assignDF = self.costDF.copy().applymap(lambda x:"") + colAgg = self.aggCol.copy().loc[1,:] + rowAgg = self.aggRow.copy().loc[:,1] + BFS = dict() + NW(1,1) - raise RuntimeWarning("minimum Cost fell through, likely casued by an inbalanced Transport Table") + return TransportTable(self.costDF,self.aggCol,self.aggRow,BFS) + + def minimumCost(self,echo: bool=False) -> TransportTable: + ''' + assign a BFS using minimum cost Method + ''' + def MC(): + ''' + private recursive function to do allocations + ''' + if echo: + df = (assignDF + .join(Series(rowAgg,name="remain")) + .append(Series(colAgg,name="remain")) + .fillna("") + ) + try: + display(df) + except: + print(df) + + if len(BFS) == len(colAgg) + len(rowAgg) - 1: + #BFS complete + return None + + # find smallest cost + row,col = tCostDF[tCostDF == tCostDF.min().min()].stack().index.tolist()[0] + + + # allocate + minVal = min(rowAgg.loc[row],colAgg.loc[col]) + assignDF.loc[row,col] = minVal + BFS[(row,col)] = minVal + colAgg.loc[col] -= minVal + rowAgg.loc[row] -= minVal + + # cross out row/col by changing it to np.inf + if colAgg.loc[col] == 0: + tCostDF.loc[:,col] = np.inf + colAgg.loc[col] = "X" + return MC() + elif rowAgg.loc[row] == 0: + tCostDF.loc[row,:] = np.inf + rowAgg.loc[row] = "X" + return MC() + + tCostDF = self.costDF.copy() + assignDF = self.costDF.copy().applymap(lambda x:"") + colAgg = self.aggCol.copy().loc[1,:] + rowAgg = self.aggRow.copy().loc[:,1] + BFS = dict() + MC() + return TransportTable(self.costDF,self.aggCol,self.aggRow,BFS) - def vogel(self,echo=False,origBHS=False,origRHS=False,crossedOut=[[],[]]): - if echo: - self.display() - if self.RHS.sum() + self.BHS.sum() ==0: - #all assignment done - if len(self.BFS) + 1 == len(self.RHS) + len(self.BHS): - return TransportTable(self.costDF.copy(),self.BFS.copy(),origRHS,origBHS) - if type(origBHS) == bool: - # only copy on first loop - # otherwise pass through orig RHS and BHS - origBHS = self.BHS.copy() - origRHS = self.RHS.copy() - - RHSOppCost = Series(dtype=float) - for i in self.RHS.index: - rowMin = self.costDF.loc[i,:].min() - rowMax = self.costDF.loc[i,:].max() - rowArr = np.sort(np.unique(self.costDF.loc[i,:])) - if (i in crossedOut[0] - or rowMin == rowMax): - RHSOppCost = RHSOppCost.append(Series([0],index=[i])) - else: - RHSOppCost = RHSOppCost.append(Series([rowArr[1] - rowArr[0]],index=[i])) - - BHSOppCost = Series(dtype=float) - for i in self.BHS.index: - colMin = self.costDF.loc[:,i].min() - colMax = self.costDF.loc[:,i].max() - colArr = np.sort(np.unique(self.costDF.loc[:,i])) - if (i in crossedOut[1] - or colMin == colMax): - BHSOppCost = BHSOppCost.append(Series([0],index=[i])) - else: - BHSOppCost = BHSOppCost.append(Series([colArr[1] - colArr[0]],index=[i])) - - maxCost = RHSOppCost.append(BHSOppCost).max() - if maxCost in RHSOppCost.array: - for rowIndex,value in RHSOppCost.items(): - if value == maxCost and rowIndex not in crossedOut[0]: - toSearch = self.costDF.loc[rowIndex].copy() - for i in self.BHS.index: - if i in crossedOut[1]: - toSearch.loc[i] = np.inf - for colIndex in self.BHS.index: - if toSearch.loc[colIndex] == min(toSearch): - nextState = self.assignMax(rowIndex,colIndex) - rowSum = origRHS.loc[rowIndex] - colSum = origBHS.loc[colIndex] - for k,v in nextState.BFS.items(): - if k[0] == rowIndex: - rowSum -= v - if k[1] == colIndex: - colSum -= v - #only cross out one row/col - if (rowSum == 0 - and len(self.RHS) - len(crossedOut[0]) > 1): - crossedOut[0].append(rowIndex) - elif (colSum == 0 # TODO testing if last case assigned is value 0 works - and len(self.BHS)-len(crossedOut[1]) > 1): # len-len is supposed to not eliminate the last in each dir - crossedOut[1].append(colIndex) - return nextState.vogel(echo=echo,origRHS=origRHS,origBHS=origBHS,crossedOut=crossedOut) - else: #maxCost in BHSOppCost - for colIndex,value in BHSOppCost.items(): - if value == maxCost and colIndex not in crossedOut[1]: - toSearch = self.costDF.loc[:,colIndex].copy() - for i in self.RHS.index: - if i in crossedOut[0]: - toSearch.loc[i] = np.inf - for rowIndex in self.RHS.index: - if toSearch.loc[rowIndex] == min(toSearch): - nextState = self.assignMax(rowIndex,colIndex) - rowSum = origRHS.loc[rowIndex] - colSum = origBHS.loc[colIndex] - for k,v in nextState.BFS.items(): - if k[0] == rowIndex: - rowSum -= v - if k[1] == colIndex: - colSum -= v - #only cross out one row/col - if (rowSum == 0 - and len(self.RHS) - len(crossedOut[0]) > 1): - crossedOut[0].append(rowIndex) - elif (colSum == 0 # TODO testing if last case assigned is value 0 works - and len(self.BHS)-len(crossedOut[1]) > 1): # len-len is supposed to not eliminate the last in each dir - crossedOut[1].append(colIndex) - return nextState.vogel(echo=echo,origRHS=origRHS,origBHS=origBHS,crossedOut=crossedOut) + def vogel(self,echo: bool=False) -> TransportTable: + ''' + assign a BFS using vogel method + ''' + def oppCost(series:Series) -> float: + s = sorted(series) + return s[1] - s[0] + def VG(): + ''' + private recursive function to do allocations + ''' + if echo: + df = (assignDF + .join(Series(rowAgg,name="remain")) + .append(Series(colAgg,name="remain")) + .fillna("") + ) + try: + display(df) + except: + print(df) + + if len(BFS) == len(colAgg) + len(rowAgg) - 1: + #BFS complete + return None + + # select grid to maximise + rowOpp = [oppCost(tCostDF.loc[i,:]) if not rowAgg[i]=="X" else 0 for i in rowAgg.index] + colOpp = [oppCost(tCostDF.loc[:,i]) if not colAgg[i]=="X" else 0 for i in colAgg.index] + maxOpp = max(rowOpp + colOpp) + if maxOpp in colOpp: + col = colOpp.index(maxOpp) + 1 + row = tCostDF.loc[:,col].idxmin() + else: # maxOpp in rowOpp + row = rowOpp.index(maxOpp) + 1 + col = tCostDF.loc[row,:].idxmin() + + # allocate + minVal = min(rowAgg.loc[row],colAgg.loc[col]) + assignDF.loc[row,col] = minVal + BFS[(row,col)] = minVal + colAgg.loc[col] -= minVal + rowAgg.loc[row] -= minVal - raise RuntimeWarning("Vogel fell through, likely casued by an inbalanced Transport Table") - - # colIndex = BHSOppCost.index(maxCost)+1 - # toSearch = self.costDF.loc[:,colIndex].copy() - # for i in self.RHS.index: - # if self.RHS.loc[i] == 0: - # toSearch.loc[i] = np.inf - # for i in self.RHS.index: - # if toSearch.loc[i] == min(toSearch): - # return (self.assignMax(i,colIndex) - # .vogel(echo=echo,origRHS=origRHS,origBHS=origBHS)) - raise RuntimeWarning("Vogel fell through, likely casued by an inbalanced Transport Table") + # cross out the row/col + if colAgg.loc[col] == 0: + tCostDF.loc[:,col] = np.inf + colAgg.loc[col] = "X" + return VG() + elif rowAgg.loc[row] == 0: + tCostDF.loc[row,:] = np.inf + rowAgg.loc[row] = "X" + return VG() + + + tCostDF = self.costDF.copy() + assignDF = self.costDF.copy().applymap(lambda x:"") + rowAgg = self.aggRow.copy().loc[:,1] + colAgg = self.aggCol.copy().loc[1,:] + BFS = dict() + VG() + return TransportTable(self.costDF,self.aggCol,self.aggRow,BFS) # solving - def solve(self,echo=False): - if echo: - self.display() - eV = self.getEnteringVar(echo=echo) - if eV == False: # is better to not run the calculation twice - return self - loop = self.findLoop((eV,),echo=echo) - return self.loopPivot(loop,echo=echo).solve(echo=echo) - - def isOptimal(self): - # reuse duality code here - return bool(self.getEnteringVar) - - def getEnteringVar(self,echo=False): - # SOLVE FOR U and V + def getEnteringVar(self,echo=False) -> Tuple[int,int]: + ''' + find the entering var + ''' + variables = [f"U{i}" for i in self.costDF.index]+[f"V{i}" for i in self.costDF.columns] + LHS = DataFrame(columns=variables) + body = [DataFrame([[1]],columns=[LHS.columns[0]])] + body += [DataFrame([[1,1]],columns=[f"U{row}",f"V{col}"]) for (row,col),c in self.BFS.items()] - matrixRHS = [] - LHSIndex = Series(self.RHS.index).copy().apply(lambda x:"U" + str(x)) - LHSIndex = LHSIndex.append(Series(self.BHS.index).copy().apply(lambda x:"V" + str(x))) - LHS_Template = Series([0]*len(LHSIndex),index=LHSIndex,dtype=int) - matrixLHS = DataFrame(columns=LHSIndex,dtype=int) - - for k,v in self.BFS.items(): - LHSrow = LHS_Template.copy() - LHSrow.loc["U" + str(k[0])] = 1 - LHSrow.loc["V" + str(k[1])] = 1 - matrixLHS = matrixLHS.append(LHSrow,ignore_index=True) - matrixRHS += [self.costDF.loc[k[0],k[1]]] - matrixRHS = Series(matrixRHS) - matrixLHS = matrixLHS.drop(columns=["U1"]) - - soln = np.linalg.solve(matrixLHS,matrixRHS) - soln = Series(soln,index=matrixLHS.columns) - soln = soln.append(Series([0],index=["U1"])) - soln = soln[LHSIndex] - - - basicVarVals = dict() - w_dict = dict() - for colIndex in self.BHS.index: - for rowIndex in self.RHS.index: - uVal = soln.loc[f'U{rowIndex}'] - vVal = soln.loc[f'V{colIndex}'] - cVal = self.costDF.loc[rowIndex,colIndex] - - uVal = int(uVal) if int(uVal) == uVal else uVal - vVal = int(vVal) if int(vVal) == vVal else vVal - cVal = int(cVal) if int(cVal) == cVal else cVal - - - basicVarVals[f"U{subscript(rowIndex)}"] = uVal - basicVarVals[f"V{subscript(colIndex)}"] = vVal - - if (rowIndex,colIndex) not in list(self.BFS.keys()): - if echo: - terms = f"U{subscript(rowIndex)} + V{subscript(colIndex)} - C{subscript(rowIndex)},{subscript(colIndex)}" - middle = f" = {uVal} + {vVal} - {cVal}" - calcValue = f" = {uVal + vVal - cVal}" - print(terms + middle + calcValue) - w_dict[rowIndex,colIndex] = soln.loc[f"U{rowIndex}"] + soln.loc[f"V{colIndex}"] - self.costDF.loc[rowIndex,colIndex] - if echo:{print(basicVarVals)} - maxW = max(w_dict.values()) + LHS = LHS.append(pd.concat(body,sort=False),sort=False).fillna(0) + LHS.index = variables + RHS = DataFrame([0]+[self.costDF.loc[row,col] for (row,col),v in self.BFS.items()],columns=["RHS"]) + RHS.index = variables + soln = np.linalg.solve(LHS,RHS) + soln = Series([i[0] for i in soln],index=variables) - if maxW == 0: - if echo:{print(f"all W values are <= 0, optimal solution found\nobjective value is {self.objectiveValue()}\n")} - return False + if echo: + print("From complementary slackness property:") + for i in LHS.index: + row = LHS.loc[i] + out = "".join([f" + {i}" for i in variables if row[i] == 1]) + f" = {RHS.loc[i,'RHS']}" + print(out[3:]) + print("\nA solution is:") + for i in soln.index: + print(f"{i} = {soln[i]}") + print("\nComputing Ui + Vj – Cij for the nonbasic variables:") + + enteringVar = (None,0) + for row in self.costDF.index: + for col in self.costDF.columns: + if (row,col) in self.BFS.keys(): + continue + calc = soln[f"U{row}"] + soln[f"V{col}"] - self.costDF.loc[row,col] + if echo:{print(f"U{row} + V{col} - C{row}{col} : X{row}{col} = {calc}")} + if calc >= enteringVar[1]: + enteringVar = ((row,col),calc) + + if enteringVar[0] == None: + # optimal solution, no complementary slackness + raise OptimalException("Optimal Solution, no Entering Variable") + if echo:{print(f"\nEntering Variable is X_{row},{col}")} + return enteringVar[0] + + def findLoop(self,start: Tuple[int,int],echo=False) -> List[Tuple[int,int]]: + ''' + recursively find _a_ valid closed loop + initial path is [(startRow,startCol)] + ''' + def FL(path: Tuple[Tuple[int,int]],axis="col") -> Tuple[Tuple[int,int]]: + if ( len(path) > 3 and path[-1] == path[0]): + if echo:{print(f"Pivot on loop of {path[:-1]}")} + return path[:-1] # drop the repeated start + + prevRow,prevCol = path[-1] + if axis == "col": + nextSteps = tuple((i,prevCol) for i in self.costDF.index) + nextSteps = tuple(i for i in nextSteps if i not in path[1:]) + nextSteps = tuple(i for i in nextSteps if i in list(self.BFS.keys()) or i == path[0]) + out = tuple(FL(path+(i,),axis="row") for i in nextSteps) + else: #axis = "row" + nextSteps = tuple((prevRow,i) for i in self.costDF.columns) + nextSteps = tuple(i for i in nextSteps if i not in path[1:]) + nextSteps = tuple(i for i in nextSteps if i in list(self.BFS.keys()) or i == path[0]) + out = tuple(FL(path+(i,),axis="col") for i in nextSteps) + out = tuple(i for i in out if i!=None) + if len(out) == 1: + return out[0] + elif len(out) > 1: + return out - for colIndex in [i[1] for i in self.BFS.keys()]: - for rowIndex in [i[0] for i in self.BFS.keys()]: - if w_dict.get((rowIndex,colIndex)) == maxW: - if echo:{print(f"Entering Variable is {(rowIndex,colIndex)}\n")} - return (rowIndex,colIndex) - - def findLoop(self,path,axis="col",echo=False): - # path is a tuple of - # tuples of (row,col) - #base case when loop closes - if ( len(path) > 3 - and path[-1] == path[0]): - if echo:{print(f"Pivot on loop of {path[:-1]}")} - return path[:-1] # drop the repeated start - prevRow,prevCol = path[-1] - - if axis == "col": - nextSteps = tuple((i,prevCol) for i in self.RHS.index) - nextSteps = tuple(i for i in nextSteps if i not in path[1:]) - nextSteps = tuple(i for i in nextSteps if i in list(self.BFS.keys()) or i == path[0]) - out = tuple(self.findLoop(path+(i,),axis="row",echo=echo) for i in nextSteps) - else: #axis = "row" - nextSteps = tuple((prevRow,i) for i in self.BHS.index) - nextSteps = tuple(i for i in nextSteps if i not in path[1:]) - nextSteps = tuple(i for i in nextSteps if i in list(self.BFS.keys()) or i == path[0]) - out = tuple(self.findLoop(path+(i,),axis="col",echo=echo) for i in nextSteps) - - out = tuple(i for i in out if i!=None) - - if len(out) == 1: - return out[0] - elif len(out) > 1: - return out - + return FL((start,)) + def loopPivot(self,loopSeq,echo=False): ''' - given that we have somehow found the even and odd squares, - now we run the loopPivot on the Transport Table + do the Pivot using loopSeq ''' - if echo:{self.display()} - evenIndexTuple = () - oddIndexTuple = () + evenCells = () # contain tuples of indexes of even cells + oddCells = () # contain tuples of indexes of odd cells i = 0 while i < len(loopSeq): - evenIndexTuple += (loopSeq[i],) - oddIndexTuple += (loopSeq[i+1],) + evenCells += (loopSeq[i],) + oddCells += (loopSeq[i+1],) i+=2 newBFS = self.BFS.copy() newBFS[loopSeq[0]] = 0 - # evenIndexTuple is an array of even squares - # oddIndexTuple is an array of odd squares - # members are tuples of (row,col) similar to x(row,col) notation - oddVals = [newBFS[i] for i in oddIndexTuple] + oddVals = [newBFS[i] for i in oddCells] theta = min(oddVals) - leavingVariable = oddIndexTuple[oddVals.index(theta)] #TODO potential degeneracy when mutiple leaving variables can exist? + leavingVariable = oddCells[oddVals.index(theta)] + #TODO potential degeneracy when mutiple leaving variables can exist? - if echo:{print(f"Leaving Variable is: X{leavingVariable}")} + if echo:{print(f"Leaving Variable is: X_{leavingVariable[0]},{leavingVariable[1]}")} - for i in oddIndexTuple: + for i in oddCells: newBFS[i] -= theta - for i in evenIndexTuple: + for i in evenCells: newBFS[i] += theta del newBFS[leavingVariable] - return TransportTable(self.costDF.copy(),newBFS,self.RHS.copy(),self.BHS.copy()) + return TransportTable(self.costDF,self.aggCol,self.aggRow,newBFS) + + def solve(self,echo: bool=False): + ''' + solve the transport tableau + ''' + newTab = self + while True: + try: + enteringVar = newTab.getEnteringVar(echo=echo) + except OptimalException: + return newTab + loop = newTab.findLoop(enteringVar,echo=echo) + newTab = newTab.loopPivot(loop,echo=echo) + if echo:{print("-"*20)} + #TODO check for degenerate looping + + def isOptimal(self) -> bool: + ''' + check if this is optimal Tableau + ''' + try: + self.getEnteringVar() + except OptimalException: + return True + return False # getters def objectiveValue(self): ''' return calculated objective value of this transport tableau using current assignments ''' - dispDF = DataFrame(np.zeros_like(self.costDF),index=self.costDF.index,columns=self.costDF.columns) - for k,v in self.BFS.items(): - dispDF.loc[k[0],k[1]] = v - dispDF.insert(dispDF.shape[1],self.RHS.name,self.RHS) - return self.costDF.mul(dispDF).sum().sum() - - #display section - def displayFormat(self): - dispDF = DataFrame(np.zeros_like(self.costDF),index=self.costDF.index,columns=self.costDF.columns) - for k,v in self.BFS.items(): - dispDF.loc[k[0],k[1]] = v - dispDF.insert(dispDF.shape[1],self.RHS.name,self.RHS) - dispColVars = self.BHS.copy() - dispColVars.name = self.BHS.name - dispDF = dispDF.append(dispColVars) - dispDF = dispDF.applymap(lambda x:"" if np.isnan(x) else x) - return dispDF - - def display(self): - displayHelper(self.displayFormat()) - + return sum([value * self.costDF.loc[row,col] for (row,col),value in self.BFS.items()]) + + def bfsDF(self) -> DataFrame: + baseDF = self.costDF.copy().applymap(lambda x:"") + for (row,col),value in self.BFS: + baseDF.loc[row,col] = value + return baseDF + def __repr__(self): - return str(self.displayFormat()) + "\n" \ No newline at end of file + baseDF = self.costDF.join(Series(self.aggRow.loc[:,1],name="row")).append(Series(self.aggCol.loc[1,:],name="col")).fillna("") + b = [f"X_{row},{col} = {c}, " for (row,col),c in self.BFS.items()] + b = "".join(b)[:-1] + return f"{baseDF}\nBFS : {b}" diff --git a/OR/TwoFactorGraph.py b/OR/TwoFactorGraph.py deleted file mode 100644 index 7e8e880..0000000 --- a/OR/TwoFactorGraph.py +++ /dev/null @@ -1,146 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -from matplotlib.patches import Polygon - -SMALL_SIZE = 8 -MEDIUM_SIZE = 10 -BIGGER_SIZE = 12 - -plt.rc('font', size=BIGGER_SIZE) # controls default text sizes -plt.rc('axes', titlesize=BIGGER_SIZE) # fontsize of the axes title -plt.rc('axes', labelsize=BIGGER_SIZE) # fontsize of the x and y labels -plt.rc('xtick', labelsize=BIGGER_SIZE) # fontsize of the tick labels -plt.rc('ytick', labelsize=BIGGER_SIZE) # fontsize of the tick labels -plt.rc('legend', fontsize=MEDIUM_SIZE) # legend fontsize -plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title - -def get_MinMaxY(myLimits): - store = dict() - for i in myLimits: - if i[0] in store.keys(): - store[i[0]] +=[i] - else: - store[i[0]] =[i] - minVals = [] - maxVals = [] - for k,v in store.items(): - maxY = max([i[1] for i in v]) - minY = min([i[1] for i in v]) - minVals.append((k,minY)) - maxVals.append((k,maxY)) - minVals.sort(key=lambda x:x[0]) - maxVals.sort(key=lambda x:x[0]) - return minVals,maxVals - -def get_MinMaxX(myLimits): - store = dict() - for i in myLimits: - if i[1] in store.keys(): - store[i[1]] +=[i] - else: - store[i[1]] =[i] - minVals = [] - maxVals = [] - for k,v in store.items(): - maxX = max([i[0] for i in v]) - minX = min([i[0] for i in v]) - minVals.append((minX,k)) - maxVals.append((maxX,k)) - - minVals.sort(key=lambda x:x[1]) - maxVals.sort(key=lambda x:x[1]) - return minVals,maxVals - -def graphical_LP(constraints,xMax=10,yMax=10,integer=False,shading=True,objectiveFunction=False): - nonegative = True - left = 0 - bottom = 0 - right=xMax - top=yMax - if integer: - yDomain = list(range(bottom,yMax+1)) #(bottom,yMax) - xDomain = list(range(left,xMax+1)) #(left,xMax) - else: - yDomain = np.linspace(bottom,yMax,1000) #(bottom,yMax) - xDomain = np.linspace(left,xMax,1000) #(left,xMax) - - fig,axes = plt.subplots() - feasible = set((x,y) for x in xDomain for y in yDomain) - - colors = plt.rcParams['axes.prop_cycle']() - for constraint,color in zip(constraints,colors): - myColor = color['color'] - if constraint[1] in (">",">=") : - myLimits = set((x,y) for x in xDomain for y in yDomain if constraint[0](x,y) >= constraint[2] ) - feasible = feasible.intersection(myLimits) - elif constraint[1] in ("<","<="): - myLimits = set((x,y) for x in xDomain for y in yDomain if constraint[0](x,y) <= constraint[2] ) - feasible = feasible.intersection(myLimits) - elif constraint[1] == "=": - myLimits = set((x,y) for x in xDomain for y in yDomain if abs(constraint[0](x,y) - constraint[2])<0.1) - feasible = feasible.intersection(myLimits) - - #feasible area for this constraint - minY,maxY = get_MinMaxY(myLimits) - minX,maxX = get_MinMaxX(myLimits) - axes.plot([i[0] for i in minY],[i[1] for i in minY],color=myColor)# ,label = constraint) - axes.plot([i[0] for i in maxY],[i[1] for i in maxY],color=myColor) - axes.plot([i[0] for i in minX],[i[1] for i in minX],color=myColor) - axes.plot([i[0] for i in maxX],[i[1] for i in maxX],color=myColor) - - if shading: - allVals = minY+maxX+list(reversed(maxY))+list(reversed(minX)) - axes.add_patch(Polygon(allVals,alpha=0.25/len(constraints),color=myColor)) - if integer: - axes.scatter( - [i[0] for i in feasible], - [i[1] for i in feasible], - color="green" - ) - - if shading: - minY,maxY = get_MinMaxY(feasible) - minX,maxX = get_MinMaxX(feasible) - allVals = minY+maxX+list(reversed(maxY))+list(reversed(minX)) - if len(allVals)>0: - axes.add_patch(Polygon(allVals,alpha=0.25,color="green")) - - - - # if objectiveFunction: - # MaxVal = 0 - # MaxCoords = [] - # for x,y in feasible: - # test = objectiveFunction(x,y) - # if test>MaxVal: - # MaxVal = test - # MaxCoords = [(x,y)] - # elif abs(test - MaxVal)<0.1: - # MaxCoords += [(x,y)] - - # xPlot = [i[0] for i in MaxCoords] - # yPlot = [i[1] for i in MaxCoords] - # plt.plot(xPlot,yPlot) - - axes.set_ylim(bottom=bottom,top=top) - axes.set_xlim(left=left,right=right) - axes.set_xlabel(u"X\u2081") - axes.set_ylabel(u"X\u2082") - axes.grid() - return fig,axes - - - -# Sample -# hw3 - -# constraints = [ -# (lambda x1,x2 :x1+2*x2, "<",350), -# (lambda x1,x2: 2*x1+x2, "<",400) -# ] -# fig,axes = graphical_LP(constraints,xMax=220,yMax=200) -# axes.plot(-1,-1,color="#1f77b4",label = u"X\u2081+2X\u2082 ≤ 350") -# axes.plot(-1,-1,color="#ff7f0e",label = u"2X\u2081+X\u2082 ≤ 400") -# axes.legend() -# fig.show() -# input("exit") \ No newline at end of file diff --git a/OR/TwoPhase.py b/OR/TwoPhase.py index 0ec7c05..913c705 100644 --- a/OR/TwoPhase.py +++ b/OR/TwoPhase.py @@ -1,65 +1,70 @@ -from .TableauBase import TableauBase -from .Format import * import numpy as np +from pandas import DataFrame + +from .TableauBase import TableauBase +from .LP import LP +from .NonNeg import NonNeg +from .ObjectiveFunction import ObjectiveFunction +from .Format import InfeasibleException + +from typing import Dict, TypeVar +TwoPhase = TypeVar("TwoPhase") class TwoPhase(TableauBase): + @classmethod - def fromLP(cls,inLP): - constraintSigns = inLP.posRHS().body.loc[:,"signs"] - baseT = super().fromLP(inLP) - body = baseT.dataFrame - origZ = body.loc[0] - newW = np.zeros_like(origZ) - newW[0] = 1 - body.loc[0] = newW - basicVar = list(baseT.basicVar) - basicVar[0] = "-W" - body.columns = ["-W"] + list(body.columns)[1:] - aVals = [] - - for i in range(1,len(constraintSigns)): - if constraintSigns[i] in (">","="): - newCol = [0] * body.shape[0] - newCol[0] = 1 - newCol[i] = 1 - body.insert(body.shape[1]-1,f"a{subscript(i)}",newCol) - basicVar[i] = f"a{subscript(i)}" - aVals.append(f"a{subscript(i)}") - - #LP to phase 1 - return TwoPhase(body,tuple(basicVar),origZ,tuple(aVals)) + def fromLP(cls,inLP:LP) -> TwoPhase: + inLP = inLP.standardForm() + # origZ = super().fromLP(inLP).dataFrame.loc[0] + # print(origZ) + origZ = inLP.objFunc + inLP.objFunc = ObjectiveFunction.new("min") + for i,c in enumerate(inLP.constraints): + if not isinstance(c,NonNeg) and c.sign in ["=",">"]: + inLP.objFunc[f"A{i+1}"] = 1 + c[f"A{i+1}"] = 1 + baseLP = super().fromLP(inLP) + return cls.fromBase(baseLP,origZ) - def fromTableauBase(self,base): - return TwoPhase(base.dataFrame,base.basicVar,self.origZ,self.aVals) + @classmethod + def fromBase(cls,base:TableauBase,origZ) -> TwoPhase: + return cls(base.dataFrame,base.basicVar,origZ) - def __init__(self,dataFrame,basicVar,origZ,aVals):#,origObjRow):#phase1z,phase2z): - super().__init__(dataFrame,basicVar,lambda x:round(x,5)) - self.origZ = origZ - self.aVals = aVals + def canonical(self,echo=False) -> TwoPhase: + t = super().canonical(echo=echo) + return self.fromBase(t,self.origZ) - def canonical(self,echo=False): + def __init__(self,dataFrame:DataFrame,basicVar:DataFrame,origZ:ObjectiveFunction): + self.dataFrame = dataFrame + self.basicVar = basicVar + self.origZ = origZ + + def canonical(self,echo=False)-> TwoPhase: t = super().canonical(echo=echo) - return self.fromTableauBase(t) + return self.fromBase(t,self.origZ) + + def phase1(self,echo=False) -> TwoPhase: + ''' + use phase 1 to assign a BFS to the new TwoPhase Object + ''' + self = self.canonical(echo=echo) + newBody = super().solve(echo=echo) + if newBody.dataFrame.loc[0,"RHS"] < 0 : + raise InfeasibleException(f"infeasible LP\n{DataFrame(newBody.dataFrame.loc[0]).T}") + return self.fromBase(newBody,self.origZ) + + def phase2(self,echo=False) -> TwoPhase: + ''' + solve the tableau normally + ''' + p2Body = self.dataFrame.copy() + p2Body = p2Body.drop([i for i in p2Body.columns if i[0] == "A"],axis=1) + p2Body = p2Body.drop([0],axis=0) + newZ = self.origZ.toDF().drop("sign",axis=1) + p2 = p2Body.append(newZ).sort_index().fillna(0) + + soln = TableauBase(p2,self.basicVar).canonical(echo=echo).solve(echo=echo) + return self.fromBase(soln,self.origZ) def solve(self,echo=False): - # create Phase1 - # done by constructor - #solve Phase1 - t = super().canonical(echo=echo) - t = t.solve(echo=echo) - # t = t.rawSolve(echo=echo) - if t.dataFrame.iloc[0,-1]>1e-15: - #infeasible - print("Infeasible") - return t - #create Phase2 - t2 = t.dataFrame.copy() - for a in self.aVals: - t2 = t2.drop(a,axis=1) - t2.columns = self.origZ.index - t2.iloc[0] = self.origZ - t2 = TableauBase(t2,(self.origZ.index[0],)+t.basicVar[1:],self.dispFunc) - #solve Phase2 - t2 = t2.canonical(echo=echo) - t2 = t2.solve(echo=echo) - return self.fromTableauBase(t2) \ No newline at end of file + return self.phase1(echo=echo).phase2(echo=echo) \ No newline at end of file diff --git a/OR/__init__.py b/OR/__init__.py index d4a6b0f..d1010e2 100644 --- a/OR/__init__.py +++ b/OR/__init__.py @@ -1,12 +1,29 @@ +#basic LP +from .Constraint import * +from .ObjectiveFunction import * +from .NonNeg import * from .LP import * + +# Simplex from .TableauBase import * from .BigM import * -from .Format import * -from .CommonLP import * -from .TwoPhase import * +from .TwoPhase import * # incomplete + +# Transport from .TransportTable import * -from .Hungarian import * -# from .TwoFactorGraph import * +from .Hungarian import * +# functional but can be improved +# planned rewrite +# inherit from transport base? + +# Network and graphs +# from .Network import * +# networkX stuff + +#Visualisation +# from .TwoFactorGraph import * -# from . import as \ No newline at end of file +# Helper Functions +# from .Format import * +# from .CommonLP import * \ No newline at end of file diff --git a/setup.py b/setup.py index e725526..3d5eb9c 100644 --- a/setup.py +++ b/setup.py @@ -22,5 +22,5 @@ "Operating System :: OS Independent", ], python_requires='>=3.6', - install_requires=['pandas',"matplotlib"] + install_requires=['pandas',"matplotlib","networkx"] ) diff --git a/tests/test_NLP.py b/tests/test_NLP.py new file mode 100644 index 0000000..d12d8c9 --- /dev/null +++ b/tests/test_NLP.py @@ -0,0 +1,13 @@ +# -*- coding: utf-8 -*- + +from context import OR +import unittest + +class BasicTestSuite(unittest.TestCase): + """Basic test cases.""" + def test_some(self): + pass + def test_other(self): + pass +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_Network.py b/tests/test_Network.py new file mode 100644 index 0000000..196df43 --- /dev/null +++ b/tests/test_Network.py @@ -0,0 +1,14 @@ +# -*- coding: utf-8 -*- + +from context import OR +# from numpy import irr,arange +import unittest + +class BasicTestSuite(unittest.TestCase): + """Basic test cases.""" + def test_some(self): + OR.Network() + def test_other(self): + pass +if __name__ == '__main__': + unittest.main()