/home/mark/model/software/ScrumPy/ScrumPy/LP/glpk.py

"""
glpk.py - (c) Mark Poolman 2008 onwards.
A high(ish) level interface between python and the glpk library of
Andrew Makhorin (mao AT gnu DOT org) which can be found at www.gnu.org/software/glpk


This module is still experimental and is

  *************************   NOT TO BE REDISTRIBUTED  ***************************

without the prior consent of the author.

"""




import os, sys, exceptions, random
import glpk_ll # _ll -> low-level interface created with Swig


########################## freindly interface to low level c arrays

from Array import BaseArray
class IntArray(BaseArray):
    _del = glpk_ll.delete_intArray
    _get = glpk_ll.intArray_getitem
    _set = glpk_ll.intArray_setitem
    _new = glpk_ll.new_intArray
    _cast = int

    #_offset=1 by default

    def __init__(self,SizeOrList):
        if type(SizeOrList) ==int:   # we're initing an undetlying malloc'd C array
            SizeOrList *= [0]        # ensure all elements are iniitialised
        BaseArray.__init__(self,SizeOrList)

class DoubleArray(BaseArray):
    _del = glpk_ll.delete_doubleArray
    _get = glpk_ll.doubleArray_getitem
    _set = glpk_ll.doubleArray_setitem
    _new = glpk_ll.new_doubleArray
    _cast = float

    def __init__(self,SizeOrList):
        if type(SizeOrList) ==int:   # we're initing an undetlying malloc'd C array
            SizeOrList *= [0.0]        # ensure all elements are iniitialised
        BaseArray.__init__(self,SizeOrList)
    #_offset=1 by default




########################## mappings between arrays/lists etc

def AddOne(l):
    return map(lambda x: x+1, l)

def SubOne(l):
    return map(lambda x: x-1, l)

def Vec2Lists(vec):
    """ pre: vec indexable type with elements comparable to 0
       post: list of indices and vals of non-zereos in vec """

    d = {}
    for i in range(len(vec)):
        val = vec[i]
        if val != 0:
            d[i] = val
    return d.keys(), d.values()

def Lists2Vec(idxs,vals,size):
    """ pre: len(vals) == len(idxs), max(idx)<size
       post: rv[indxs[i]] = vals[i], i[0..size-1] else 0.0 """


    rv = [0.0] * size
    for i in range(len(idxs)):
        idx = idxs[i]
        if idx >=0:
            rv[idx] = vals[i]
    return rv



#########################


VerySmall = 1e-12
def IsZero(number, tol=VerySmall):
    return abs(number) < tol




########################## mappings between various constants

Sym2Num = {}             # map between symbolic constants and their values
Num2Sym = {}
Str2Num = {}
Num2Str = {             # map constant status values to human strings
    glpk_ll.LPX_MIN : "Min",
    glpk_ll.LPX_MAX : "Max"
}

Status2Str = {           # map symbolic status constants to human strings
    "LPX_OPT"    : "Optimal",
    "LPX_FEAS"   : "Feasible",
    "LPX_INFEAS" : "Infeasible",
    "LPX_NOFEAS" : "No feasible",
    "LPX_UNDEF"  : "Undefined",
    'LPX_UNBND'  : "Unbounded!"
}


for d in dir(glpk_ll):
    if d.startswith("LPX_"):
        num =  getattr(glpk_ll,d)
        if type(num)==int:
            Sym2Num[d] = num
            Num2Sym[num] = d

for k in Status2Str.keys():
    Num2Str[Sym2Num[k]] = Status2Str[k] # generic status

    kp = k.replace("_","_P_")        # primal status if present
    if Sym2Num.has_key(kp):
        Status2Str[kp] = Num2Str[Sym2Num[kp]] = "Primal "+ Status2Str[k]

    kd = k.replace("_","_D_")         # dual status if present
    if Sym2Num.has_key(kd):
        Status2Str[kd] = Num2Str[Sym2Num[kd]] = "Dual " +Status2Str[k]

for k in Num2Str.keys():
    Str2Num[Num2Str[k]] = k


#
# class that executes a glpk_ll function with fixed first arg
# and arbitrary subsequent args
#
class exec_ll:
        def __init__(self,name, arg0):
            self.name = name
            self.arg0 = arg0

        def __call__(self, *args):
            return getattr(glpk_ll, self.name)(self.arg0,*args)




#
# IO utils
#
#
def F2Str(fun):
    """ pre: fun(filename) writes to a file
       post: returns as a string what fun was going to write """

    dest = str(random.randint(0,100000))
    fun(dest)
    rv = open(dest).read()
    os.remove(dest)
    return rv


#
##
### ##################### main user interface to glpk
##
#
class lp:


    _Readers = {      # map humanoid format descriptions to glpk function names
        "FixedMPS"    : glpk_ll._glp_lpx_read_mps,
        "FreeMPS"     : glpk_ll._glp_lpx_read_freemps,
        "CPLEX"       : glpk_ll._glp_lpx_read_cpxlp,
        "GnuMathProg" : glpk_ll._glp_lpx_read_model
    }

    _MethodsMap = {   # map methods to names of functions that have trivial equivalence
                      # methods can't be bound until we have an lpx instance - see __init__()
        "SetName":"_glp_lpx_set_prob_name",
        "GetName":"_glp_lpx_get_prob_name",

        "SetObjDir":"_glp_lpx_set_obj_dir",
        "GetObjDir":"_glp_lpx_get_obj_dir",

        "GetObjVal" :"_glp_lpx_get_obj_val",

        "SetObjName":"_glp_lpx_set_obj_name",
        "GetObjName":"_glp_lpx_get_obj_name",

        "GetNumRows":"_glp_lpx_get_num_rows",
        "GetNumCols":"_glp_lpx_get_num_cols"
    }


    _StatusGettersMap = {  # as above, these query problem status
        "Generic": "_glp_lpx_get_status",
        "Primal" : "_glp_lpx_get_prim_stat",
        "Dual"   : "_glp_lpx_get_dual_stat"
    }
    _StatusGetters = {} # will be bound here on __init__()


    _WritersMap = {  # map nice warm fuzzy pythonic string format descriptions to glpk functions
        "Human" :"_glp_lpx_print_prob",
        "FixedMPS" : "_glp_lpx_write_mps",
        "FreeMPS"  : "_glp_lpx_write_freemps",
        "CPLEX"    : "_glp_lpx_write_cpxlp",
        "Solution" : "_glp_lpx_print_sol"
    }
    _Writers = {} # as previous maps




    def __init__(self, name="Unamed"):
        self.lpx = glpk_ll._glp_lpx_create_prob()
        self.__setup()

        self.SetName(name)


    def __setup(self):
        self.rnames={}
        self.cnames={}
        self.cidxs = {}
        self.ridxs = {}

        self.__map_ll()


        for method in self._MethodsMap:
            setattr(self, method, exec_ll(self._MethodsMap[method],self.lpx))

        for writer in self._WritersMap:
            self._Writers[writer] = exec_ll(self._WritersMap[writer],self.lpx)

        for getter in self._StatusGettersMap:
            self._StatusGetters[getter] =  exec_ll(self._StatusGettersMap[getter],self.lpx)


    def __repr__(self):
        return self.GetName()+" (instance of glpk.lp)"

    def __str__(self):
        return self.Write()


    def __map_ll(self):

        for i in dir(glpk_ll):
            if i.startswith("glp_lpx") and not i.startswith("glp_lpx_read"):
                setattr(self, "_"+i, exec_ll(i,self.lpx))
            elif i.startswith("_glp_lpx"):
                setattr(self, i, exec_ll(i,self.lpx))


    #def __getattr__(self, attr):
    #    return
        # "just-in-time" alternative to __map_ll
        # the approach, and variants thereof, will almost certainly
        # cause more problems than it solves.
        # left here for your entertainment.
        #if attr.startswith("_glp_lpx"):
        #    a = exec_ll(attr[1:], self.lpx)
        #    setattr(self, attr,a)
        #    return a



    def __diccheck(self, dic, k):
        if not k in dic:
            sys.stderr.write(k+" is invalid - options are: "+", ".join(dic.keys()))
            raise exceptions.KeyError
        return dic[k]


    def __rowcheck(self, r):
        if r<0 or r>=self.GetNumRows():
            raise exceptions.IndexError, str(r)
        return r+1

    def __colcheck(self, c):
        if c<0 or c>=self.GetNumCols():
            raise exceptions.IndexError, str(c)
        return c+1


    def __getrow(self,r):
        if type(r)==str:
            if self.ridxs.has_key(r):
                r = self.ridxs[r]
            else:
                raise exceptions.KeyError, str(r)

        return self.__rowcheck(r)


    def __getcol(self,c):
        if type(c)==str:
            if self.cidxs.has_key(c):
                c = self.cidxs[c]
            else:
                raise exceptions.KeyError, str(c)

        return self.__colcheck(c)



    def __getbnds(self, lo, hi):

        if   lo == None and hi == None:
            flag, lo, hi = glpk_ll.LPX_FR, 0.0, 0.0
        elif lo != None and hi == None:
            flag, lo, hi = glpk_ll.LPX_LO, lo, 0.0
        elif lo == None and hi != None:
            flag, lo, hi = glpk_ll.LPX_UP, 0.0, hi
        elif lo != None and hi != None:
            if lo < hi:
                flag =  glpk_ll.LPX_DB
            elif lo==hi:
                flag = glpk_ll.LPX_FX
            else:
                raise exceptions.ValueError, "lo > hi !"

        return flag, lo, hi



    def __RebuildNameDics(self):

        self.rnames = {}
        self.cnames = {}
        self.ridxs = {}
        self.cidxs = {}
        nr,nc = self.GetDims()

        for r in range(nr):
            rname = self._glp_lpx_get_row_name(r+1)
            self.rnames[r] =  rname
            self.ridxs[rname] = r

        for c in range(nc):
            cname = self._glp_lpx_get_col_name(c+1)
            self.cnames[c] =  cname
            self.cidxs[cname] = c


    def __del__(self):
            self._glp_lpx_delete_prob()

    #
    ## # Public methods start here
    #


    #
    ## define the problem
    #

    def SetObjDirec(self, direc="Min"):
        """ pre: direc = "Min" | "Max" """

        if not (direc=="Max" or direc=="Min"):
            raise exceptions.ValueError, direc

        #if direc == "Min":
        #    direc = glpk_ll.LPX_MIN
        #elif direc == "Max":
        #    direc = glpk_ll.LPX_MAX
        #else:
        #    raise exceptions.ValueError, direc
        self._glp_lpx_set_obj_dir(Str2Num[direc])


    def AddRows(self,rows):

        tr = type(rows)
        if tr == int:
            if rows <1:
                raise exceptions.ValueError, "Can't add less then one row !"
            self._glp_lpx_add_rows(rows)

        elif tr == list:
            oldnr = self.GetNumRows()
            newnr = len(rows)
            if True in map(self.HasRow,rows):
                raise exceptions.ValueError, "Can't add duplicate row index !"

            self._glp_lpx_add_rows(newnr)
            for ridx in range(newnr):
                self.SetRowName(ridx+oldnr, rows[ridx])

        else:
            raise execeptions.TypeError, "rows must be int or list"



    def AddCols(self,cols):

        tc = type(cols)
        if tc == int:
            if cols <1:
                raise exceptions.ValueError, "Can't add less then one col !"
            self._glp_lpx_add_cols(cols)

        elif tc == list:
            oldnc = self.GetNumCols()
            newnc = len(cols)
            if True in map(self.HasCol,cols):
                raise exceptions.ValueError, "Can't add duplicate col index !"

            self._glp_lpx_add_cols(newnc)
            for cidx in range(newnc):
                self.SetColName(cidx+oldnc, cols[cidx])

        else:
            raise execeptions.TypeError, "cols must be int or list"


    def GetRowIdxs(self,rows):
        return map(self.__getrow, rows)

    def GetColIdxs(self,cols):
        return map(self.__getcol, cols)


    def DelRows(self,rows):
        n = len(rows)
        rows = self.GetRowIdxs(rows)
        idxs = IntArray(rows)
        self._glp_lpx_del_rows(n,idxs._array)
        self.__RebuildNameDics()


    def DelCols(self,cols):
        n = len(cols)
        cols = self.GetColIdxs(cols)
        idxs = IntArray(cols)
        self._glp_lpx_del_cols(n,idxs._array)
        self.__RebuildNameDics()


    def SetRowName(self, r, name):

        if self.ridxs.has_key(name):
            raise exceptions.KeyError, "duplicate row name, "+name

        gr = self.__getrow(r)
        self._glp_lpx_set_row_name(gr,name)
        self.ridxs[name] = r
        self.rnames[r] = name


    def SetColName(self, c, name):

        if self.cidxs.has_key(name):
            raise exceptions.KeyError, "duplicate col name,"+name

        gc = self.__getcol(c)
        self._glp_lpx_set_col_name(gc,name)
        self.cidxs[name] = c
        self.cnames[c] = name



    def SetRowBounds(self, r, lo=None, hi=None):

        r = self.__getrow(r)
        flag, lo, hi = self.__getbnds(lo,hi)
        self._glp_lpx_set_row_bnds(r, flag, lo, hi)


    def SetColBounds(self, c, lo=None, hi=None):

        c = self.__getcol(c)
        flag, lo, hi = self.__getbnds(lo,hi)
        self._glp_lpx_set_col_bnds(c,flag,lo,hi)



    def SetRowVals(self, row, vals):

        row = self.__getrow(row)
        lenv = len(vals)
        if len(vals) != self.GetNumCols():
            raise exceptions.IndexError, "incorrect row length"
        zs = []
        nzs = []
        for i in range(lenv):
            if vals[i] == 0:
                zs.append(i)
            else:
                nzs.append(i)
        zs.reverse()
        if len(zs) !=0:
            vals = vals[:] # copy, because we are going to change the contents
        for i in zs:
            del vals[i]
        nzs = IntArray(AddOne(nzs))
        vals = DoubleArray(vals)
        self._glp_lpx_set_mat_row(row, len(vals), nzs._array, vals._array)




    def LoadMtxFromILists(self, RowIdxs, ColIdxs, ElVals):

        nels = len(ElVals)
        nr,nc = self.GetDims()

        if not len(RowIdxs) == len(ColIdxs) == nels:
            raise exceptions.IndexError, "Inconsistant indices and/or values"

        try:
            RowIdxs = map(self.__getrow, RowIdxs)
        except:
            raise exceptions.IndexError, "Bad row index"

        try:
            ColIdxs = map(self.__getcol, ColIdxs)
        except:
            raise exceptions.IndexError, "Bad col index"

        tupdic = {}
        for i in range(len(RowIdxs)):
            tup = str((RowIdxs[i],ColIdxs[i]))
            if tupdic.has_key(tup):
                raise exceptions.IndexError, "Duplicate indices not allowed"
                # don't ask me why not, I didn't write gplk.c !
            else:
                tupdic[tup]=1


        zs = []                     # remove any zero element values, which for reasons unknown
        for i in range(nels):       # will cause precipitate termination in glpk
            if ElVals[i] ==0:
                zs.append(i)
        if len(zs) != 0:
            RowIdxs, ColIdxs, ElVals = RowIdxs[:], ColIdxs[:], ElVals[:] # copy, keep orignals untouched
            zs.reverse()
            for i in zs:
                for l in RowIdxs, ColIdxs, ElVals:
                    del l[i]


        ia = IntArray(RowIdxs)
        ja = IntArray(ColIdxs)
        ar = DoubleArray(ElVals)

        self._glp_lpx_load_matrix(
            len(ElVals),
            ia._array,
            ja._array,
            ar._array)


    def LoadMtxFromDic(self,dic, KeysAreRowNames=True):

        if KeysAreRowNames:
            for k in dic.keys():
                self.SetRowVals(k,dic[k])
        else:
            print "D'oh !!"


    def SetObjCoef(self, c, coef):

        if c==-1 or c=="shift":
            self._glp_lpx_set_obj_coef(0, coef)

        c = self.__getcol(c)
        self._glp_lpx_set_obj_coef(c, coef)


    def SetObjCoefsFromLists(self, cs, coefs):

        lenc = len(cs)
        if lenc != len(coefs):
            raise exceptions.IndexError, "cs and coefs not of same length !"

        cs = map(self.__getcol,cs)
        for i in range(lenc):
            self._glp_lpx_set_obj_coef(cs[i], coefs[i])


    def SetObjCoefsFromDic(self, d):

        self.SetObjCoefsFromLists(d.keys(), d.values())



    #
    ##  problem information
    #

    def GetDims(self):
        return self.GetNumRows(), self.GetNumCols()

    def HasRow(self, r):
        """ pre: True
           post: HasRow(r) => r is a valid row index """
        try:
            self.__getrow(r)
            return True
        except:
            return False


    def HasCol(self, c):
        """pre: True
          post: HasCol(c) => c is a valid col index """
        try:
            self.__getcol(c)
            return True
        except:
            return False


    def GetRow(self, r):
        """ self.HasRow(r) => return row r as a list
            else: exception """

        r = self.__getrow(r)
        lenr = self.GetNumCols()
        inds = IntArray(lenr)
        vals = DoubleArray(lenr)
        self._glp_lpx_get_mat_row(r, inds._array, vals._array)
        inds = SubOne(inds.ToList())
        vals = vals.ToList()

        return Lists2Vec(inds, vals, lenr)


    def GetRowsAsLists(self):

        nr = self.GetDims()[0]
        return map(self.GetRow, range(nr))




    def GetCol(self, c):

        c = self.__getcol(c)
        lenc = self.GetNumCols()
        inds = IntArray(lenc)
        vals = DoubleArray(lenc)
        self._glp_lpx_get_mat_col(c, inds._array, vals._array)
        inds = SubOne(inds.ToList())
        vals = vals.ToList()

        return Lists2Vec(inds, vals, lenc)


    def GetRowName(self, r):
        self.__getrow(r) # sanity check, exception if r invalid

        if self.rnames.has_key(r):
            return self.rnames[r]
        return "row_"+str(r)


    def GetRowNames(self):

        return map(self.GetRowName, range(self.GetDims()[0]))

    def GetColName(self, c):
        self.__getcol(c) # sanity check, exception if r invalid

        if self.cnames.has_key(c):
            return self.cnames[c]
        return "col_"+str(c)


    def GetColNames(self):

        return map(self.GetColName, range(self.GetDims()[1]))








    #
    ## Solve the problem
    #



    def Solve(self, Method="Simplex"):

        self._glp_lpx_std_basis()
        self._glp_lpx_simplex()



    #
    ## Status information
    #

    def GetStatus(self, Stype="Generic"):
        """pre: Stype in ["Generic", "Primal", "Dual"]
          post: numerical value of current status """

        return self.__diccheck(self._StatusGetters, Stype)()


    def GetStatusMsg(self,Stype="Generic"):
        """pre: self.GetStatus(Stype)
          post: human string describing current status """
        return Num2Str[self.GetStatus(Stype)]


    def GetStatusSym(self,Stype="Generic"):
        """pre: self.GetStatus(Stype)
          post: symbolic string describing current status """
        return Num2Sym[self.GetStatus(Stype)]


    #
    ## Solution information
    #
    # see also: GetObjVal and GetObjName defined by __init__()

    def GetRowPrimal(self, r):
        r = self.__getrow(r)
        return self._glp_lpx_get_row_prim(r)

    def GetRowDual(self, r):
        r = self.__getrow(r)
        return self._glp_lpx_get_row_dual(r)



    def GetColPrimal(self, c):
        c = self.__getcol(c)
        return self._glp_lpx_get_col_prim(c)

    def GetColDual(self, c):
        c = self.__getcol(c)
        return self._glp_lpx_get_col_dual(c)



    #
    ##  I/O methods
    #

    def Read(self, src, format='CPLEX'):
        """pre: format in ['GnuMathProg', 'FixedMPS', 'CPLEX', 'FreeMPS']
                src is readable filename in format format
          post: self represents problem in src """


        self._glp_lpx_delete_prob()
        self.lpx = self.__diccheck(self._Readers, format)(src)
        self.__setup()
        self.__RebuildNameDics()



    def Write(self, format="Human", dest=None):
        """ pre: format in ['FixedMPS', 'CPLEX', 'Human', 'FreeMPS', 'Solution']
                 Dest writable file name OR None
            post: Dest != None -> problem written in format to Dest
                  Dest == None -> returns proble as a string in format """


        writer = self.__diccheck(self._Writers, format)
        nr,nc = self.GetDims()

        if (nc < 1 or nr <1) and format != "Human":
            raise exceptions.IndexError, "must have at least one row and column"

        if dest==None:
            return F2Str(writer)
        else:
            writer(dest)
            rv = None
        return rv




    def PrintMtx(self):
        """ print current constraint matrix on stdout """

        rv =""
        nr,nc = self.GetDims()
        print " "*10,
        for c in range(nc):
            print self.GetColName(c).ljust(6),
        print
        for r in range(nr):
            print self.GetRowName(r).ljust(10), " ".join(map(lambda x:("%3.3g"%x).ljust(6), self.GetRow(r)))







"""
lp = glpk_ll.glp_lpx_create_prob()
glpk_ll.glp_lpx_set_prob_name(lp, "sample")
glpk_ll.glp_lpx_set_obj_dir(lp, glpk_ll.LPX_MAX)
glpk_ll.glp_lpx_add_rows(lp, 3)
glpk_ll.glp_lpx_set_row_name(lp, 1, "p")
glpk_ll.glp_lpx_set_row_bnds(lp, 1, glpk_ll.LPX_UP, 0.0, 100.0)
glpk_ll.glp_lpx_set_row_name(lp, 2, "q")
glpk_ll.glp_lpx_set_row_bnds(lp, 2, glpk_ll.LPX_UP, 0.0, 600.0)
glpk_ll.glp_lpx_set_row_name(lp, 3, "r")
glpk_ll.glp_lpx_set_row_bnds(lp, 3, glpk_ll.LPX_UP, 0.0, 300.0)
glpk_ll.glp_lpx_add_cols(lp, 3)
glpk_ll.glp_lpx_set_col_name(lp, 1, "x1")
glpk_ll.glp_lpx_set_col_bnds(lp, 1, glpk_ll.LPX_LO, 0.0, 0.0)
glpk_ll.glp_lpx_set_obj_coef(lp, 1, 10.0)
glpk_ll.glp_lpx_set_col_name(lp, 2, "x2")
glpk_ll.glp_lpx_set_col_bnds(lp, 2, glpk_ll.LPX_LO, 0.0, 0.0)
glpk_ll.glp_lpx_set_obj_coef(lp, 2, 6.0)
glpk_ll.glp_lpx_set_col_name(lp, 3, "x3")
glpk_ll.glp_lpx_set_col_bnds(lp, 3, glpk_ll.LPX_LO, 0.0, 0.0)
glpk_ll.glp_lpx_set_obj_coef(lp, 3, 4.0)


ia = IntArray([1,1,1,2,3,2,3,2,3])
ja = IntArray([1,2,3,1,1,2,2,3,3])
ar = DoubleArray([1,1,1,10,2,4,2,5,6])
#ia[1] = 1 ; ja[1] = 1 ; ar[1] =  1.0 # /* a[1,1] =  1 */
#ia[2] = 1 ; ja[2] = 2 ; ar[2] =  1.0 # /* a[1 ;2] =  1 */
#ia[3] = 1 ; ja[3] = 3 ; ar[3] =  1.0 # /* a[1 ;3] =  1 */
#ia[4] = 2 ; ja[4] = 1 ; ar[4] = 10.0 # /* a[2 ;1] = 10 */
#ia[5] = 3 ; ja[5] = 1 ; ar[5] =  2.0 # /* a[3 ;1] =  2 */
#ia[6] = 2 ; ja[6] = 2 ; ar[6] =  4.0 # /* a[2 ;2] =  4 */
#ia[7] = 3 ; ja[7] = 2 ; ar[7] =  2.0 # /* a[3 ;2] =  2 */
#ia[8] = 2 ; ja[8] = 3 ; ar[8] =  5.0 # /* a[2 ;3] =  5 */
#ia[9] = 3 ; ja[9] = 3 ; ar[9] =  6.0 # /* a[3 ;3] =  6 */

glpk_ll.glp_lpx_load_matrix(lp, 9, ia._array, ja._array, ar._array)
glpk_ll.glp_lpx_simplex(lp)
Z = glpk_ll.glp_lpx_get_obj_val(lp)
x1 = glpk_ll.glp_lpx_get_col_prim(lp, 1)
x2 = glpk_ll.glp_lpx_get_col_prim(lp, 2)
x3 = glpk_ll.glp_lpx_get_col_prim(lp, 3)
print "\nZ =", Z, " x1 = ",x1," x2 = ",x2," x3 = ", x3
#glpk_ll.glp_lpx_delete_prob(lp)
 #     return 0
"""