src/LinearProgramming.py

   1 """
   2 Classes to create, solve, and make dinner for linear programs. Handles
   3 integration with lp_solve.
   4 """
   5
   6 import fractions
   7 from numpy import *
   8 import os
   9 import site
  10 import sys
  11
  12 # Add LP_SOLVE_PATH to our path. There is no point to this variable
  13 # other than to make the site.addsitedir() line fit within 80
  14 # characters.
  15 LP_SOLVE_PATH1 = '/../../lib/lp_solve'
  16 LP_SOLVE_PATH2 = '/../lib/lp_solve'
  17 site.addsitedir(os.path.dirname(os.path.abspath(sys.argv[0])) + LP_SOLVE_PATH1)
  18 site.addsitedir(os.path.dirname(os.path.abspath(sys.argv[0])) + LP_SOLVE_PATH2)
  19
  20 from lpsolve55 import *
  21
  22
  23 # Constants representing the two types of linear programs.
  24 # MINIMIZE means that we would like to minimize the objective
  25 # function, and MAXIMIZE means that we would like to maximize it.
  26 MINIMIZE = 0
  27 MAXIMIZE = 1
  28
  29
  30 class LinearProgram(object):
  31     """
  32     Represents an instance of an lp_solve linear program.
  33     The actual lp_solve linear program is only created when it
  34     is needed, and modifications to it are cached beforehand.
  35     """
  36
  37
  38     def get_row_count(self):
  39         """
  40         Return the number of rows in the constraint matrix.
  41         """
  42         return len(self.constraint_matrix)
  43
  44
  45     def get_column_count(self):
  46         """
  47         Return the number of columns in the constraint matrix.
  48         If we don't have any rows yet, claim zero columns as well.
  49         """
  50         if self.get_row_count() == 0:
  51             return 0
  52         else:
  53             return len(self.constraint_matrix[0])
  54
  55
  56
  57     @property
  58     def type(self):
  59         """
  60         A property representing the type of linear program, either
  61         MINIMIZE or MAXIMIZE.
  62         """
  63         return self._type
  64
  65     @type.setter
  66     def type(self, type):
  67         if type == MINIMIZE:
  68             self._type = MINIMIZE
  69             if self._lp != None:
  70                 lpsolve('set_minim', self._lp)
  71         else:
  72             self._type = MAXIMIZE
  73             if self._lp != None:
  74                 lpsolve('set_maxim', self._lp)
  75
  76
  77
  78     @property
  79     def objective_coefficients(self):
  80         return self._objective_coefficients
  81
  82
  83     @objective_coefficients.setter
  84     def objective_coefficients(self, value):
  85         self._objective_coefficients = value
  86
  87         if self._lp != None:
  88             lpsolve('set_obj_fn',
  89                     self._lp,
  90                     self._objective_coefficients)
  91
  92
  93
  94     @property
  95     def constraint_matrix(self):
  96         return self._constraint_matrix
  97
  98     @constraint_matrix.setter
  99     def constraint_matrix(self, value):
 100         self._constraint_matrix = value
 101
 102         if self._lp != None:
 103             lpsolve('set_mat', self._lp, value)
 104
 105
 106
 107     @property
 108     def rhs(self):
 109         return self._rhs
 110
 111     @rhs.setter
 112     def rhs(self, value):
 113         self._rhs = value
 114
 115         if self._lp != None:
 116             lpsolve('set_rh_vec', self._lp, self._rhs)
 117
 118
 119
 120     @property
 121     def inequalities(self):
 122         return self._inequalities
 123
 124     @inequalities.setter
 125     def inequalities(self, value):
 126         self._inequalities = value
 127
 128         if self._lp != None:
 129             for i in range(self.get_row_count()):
 130                 lpsolve('set_constr_type', self._lp, i+1, value[i])
 131
 132
 133     @property
 134     def solution_lower_bounds(self):
 135         return self._solution_lower_bounds
 136
 137     @solution_lower_bounds.setter
 138     def solution_lower_bounds(self, value):
 139         if len(value) != self.get_column_count():
 140             return
 141
 142         self._solution_lower_bounds = value
 143
 144         if self._lp != None:
 145             for i in range(self.get_column_count()):
 146                 lpsolve('set_lowbo', self._lp, i+1, value[i])
 147
 148
 149
 150     @property
 151     def solution_upper_bounds(self):
 152         return self._solution_upper_bounds
 153
 154
 155     @solution_upper_bounds.setter
 156     def solution_upper_bounds(self, value):
 157         if len(value) != self.get_column_count():
 158             return
 159
 160         self._solution_upper_bounds = value
 161
 162         if self._lp != None:
 163             for i in range(self.get_column_count()):
 164                 lpsolve('set_upbo', self._lp, i+1, value[i])
 165
 166
 167
 168     @property
 169     def integer_variables(self):
 170         """
 171         A vector containing the indices of any solution variables
 172         which must be integers.
 173         """
 174         return self._integer_variables
 175
 176     @integer_variables.setter
 177     def integer_variables(self, value):
 178         self._integer_variables = value
 179
 180         if self._lp != None:
 181             for i in range(len(value)):
 182                 lpsolve('set_int', self._lp, value[i], 1)
 183
 184
 185
 186     def make_all_variables_integers(self):
 187         """
 188         Force all solution variables to be integers. This is achieved
 189         by filling the integer_variables vector with all possible
 190         indices.
 191         """
 192         ivs = []
 193
 194         for i in range(self.get_column_count()):
 195             ivs.append(i+1)
 196             if self._lp != None:
 197                 lpsolve('set_int', self._lp, i+1, 1)
 198
 199         self.integer_variables = ivs
 200
 201
 202
 203     @property
 204     def scale_mode(self):
 205         """
 206         The scaling mode used for handling floating point numbers.
 207         See <http://lpsolve.sourceforge.net/5.5/scaling.htm> for more
 208         information.
 209         """
 210         return self._scale_mode
 211
 212
 213     @scale_mode.setter
 214     def scale_mode(self, value):
 215         self._scale_mode = value
 216
 217         if self._lp != None:
 218             lpsolve('set_scaling', self._lp, value)
 219
 220
 221
 222     def print_tableau(self):
 223         """
 224         Tell lp_solve to print its simplex tableau. Only works after
 225         a successful call to solve().
 226         """
 227         lpsolve('set_outputfile', self._lp, '')
 228         lpsolve('print_tableau', self._lp)
 229
 230
 231     def __init__(self):
 232         """
 233         Initialize the object, setting all of the properties
 234         either empty or to sane defaults.
 235
 236         The _lp variable is set to None, initially. All of the
 237         property setters will test for _lp == None, and will refuse
 238         to make calls to lp_solve if that is the case. A new instance
 239         of an lp_solve linear program will be created (and stored in
 240         the _lp variable) on demand.
 241
 242         If the _lp variable is *not* None, the property setters will
 243         make calls to lp_solve, updating the pre-existing linear program.
 244         """
 245
 246         self._lp = None
 247         self._objective_coefficients = []
 248         self._constraint_matrix = []
 249         self._rhs = []
 250         self._inequalities = []
 251         self._integer_variables = []
 252         self._solution_lower_bounds = []
 253         self._solution_upper_bounds = []
 254         self._scale_mode = 0
 255         self._type = MINIMIZE
 256
 257
 258     def set_all_lp_properties(self):
 259         """
 260         Re-set all linear program properties. After a new linear
 261         program is created, it will be 'empty'. We already have
 262         its properties stored in our member variables, however,
 263         we need to make calls to lp_solve to set them on the new
 264         linear program instance.
 265
 266         All of the property setters will check for the existence of
 267         self._lp and make calls to lp_solve as necessary. So, to set
 268         all of our properties, we just have to trigger the property
 269         setters a second time.
 270         """
 271         self.constraint_matrix = self.constraint_matrix
 272         self.rhs = self.rhs
 273         self.objective_coefficients = self.objective_coefficients
 274         self.inequalities = self.inequalities
 275         self.integer_variables = self.integer_variables
 276         self.solution_lower_bounds = self.solution_lower_bounds
 277         self.solution_upper_bounds = self.solution_upper_bounds
 278         self.scale_mode = self.scale_mode
 279         self.type = self.type
 280
 281
 282
 283     def delete(self):
 284         if self._lp != None:
 285             lpsolve('delete_lp', self._lp)
 286
 287
 288
 289     def create_lp_if_necessary(self):
 290         """
 291         If we already have a linear program instance, do nothing.
 292         Otherwise, create one, and set all of the necessary properties.
 293         """
 294         if self._lp != None:
 295             return
 296
 297         self._lp = lpsolve('make_lp',
 298                            self.get_row_count(),
 299                            self.get_column_count())
 300
 301         # This is not critical, but it will encourage lp_solve to
 302         # warn us about potential problems.
 303         lpsolve('set_verbose', self._lp, IMPORTANT)
 304
 305         self.set_all_lp_properties()
 306
 307
 308
 309     def solve(self):
 310         """
 311         Solve the linear program. The lp_solve instance is
 312         created beforehand if necessary.
 313         """
 314         self.create_lp_if_necessary()
 315         result = lpsolve('solve', self._lp)
 316
 317         # Default to empty return values.
 318         obj = []
 319         x = []
 320         duals = []
 321
 322         # See http://lpsolve.sourceforge.net/5.5/solve.htm for a
 323         # description of these constants.
 324         if (result == OPTIMAL    or
 325             result == SUBOPTIMAL or
 326             result == PROCBREAK  or
 327             result == FEASFOUND):
 328
 329             # If the result was "good," i.e. if get_solution will work,
 330             # call it and use its return value as ours.
 331             [obj, x, duals, ret] = lpsolve('get_solution', self._lp)
 332
 333         return [obj, x, duals]
 334
 335
 336     def objective_coefficient_gcd(self):
 337         """
 338         Return the GCD of all objective function coefficients.
 339         """
 340         return reduce(fractions.gcd, self.objective_coefficients)
 341
 342
 343
 344 class SimplexIteration(object):
 345     """
 346     Represents the 'current' iteration of the simplex method at some
 347     point. It needs an A,b,x, and c corresponding to the linear
 348     program in standard form. It can then determine whether or not the
 349     current vertex (x) is optimal, and possible move to a better one.
 350     """
 351
 352     @property
 353     def constraint_matrix(self):
 354         return self._constraint_matrix.tolist()
 355
 356     @constraint_matrix.setter
 357     def constraint_matrix(self, value):
 358         self._constraint_matrix = matrix(value)
 359
 360     @property
 361     def rhs(self):
 362         return self._rhs
 363
 364     @rhs.setter
 365     def rhs(self, value):
 366         self._rhs = value
 367
 368     @property
 369     def objective_coefficients(self):
 370         return self._objective_coefficients
 371
 372     @objective_coefficients.setter
 373     def objective_coefficients(self, value):
 374         self._objective_coefficients = value
 375
 376     @property
 377     def solution_vector(self):
 378         return self._solution_vector.tolist()
 379
 380     @solution_vector.setter
 381     def solution_vector(self, value):
 382         self._solution_vector = array(value)
 383
 384
 385     @property
 386     def basic_variables(self):
 387         # The current set of basic variables. Constructed from the
 388         # "true" source, our list of basic indices.
 389         idxs = self.basic_indices
 390         return map(lambda x: "x" + str(x+1), idxs)
 391
 392
 393     @basic_variables.setter
 394     def basic_variables(self, value):
 395         """
 396         Syntactic sugar to set the basic indices. We take a string
 397         of the form x1,x2,...xn, and subtract one from each of the
 398         subscripts to get the basic indices.
 399         """
 400         basic_indices = []
 401         vars = value.split(',')
 402
 403         for var in vars:
 404             var = var.strip()
 405             var = var.replace('x', '')
 406             basic_indices.append(int(var)-1)
 407
 408         self.basic_indices = basic_indices
 409
 410
 411     @property
 412     def nonbasic_variables(self):
 413         # All elements of the solution vector that have value zero.
 414         idxs = self.nonbasic_indices
 415         return map(lambda x: "x" + str(x+1), idxs)
 416
 417
 418     @property
 419     def basic_indices(self):
 420         return self._basic_indices
 421
 422
 423     @basic_indices.setter
 424     def basic_indices(self, value):
 425         self._basic_indices = value
 426
 427
 428     @property
 429     def nonbasic_indices(self):
 430         all_indices = range(0, len(self.solution_vector))
 431         return list(set(all_indices) - set(self.basic_indices))
 432
 433
 434     @property
 435     def optimal(self):
 436         # True if the current solution is optimal, false if not.
 437         for idx in self.nonbasic_indices:
 438             if self.reduced_cost(idx) < 0.0:
 439                 return False
 440
 441         return True
 442
 443
 444     def constraint_column(self, idx):
 445         # Return the column of the constraint matrix corresponding
 446         # to index idx.
 447         bm = matrix(self.constraint_matrix).transpose().tolist()
 448         return bm[idx]
 449
 450     def negative_constraint_column(self, idx):
 451         # Return the column of the constraint matrix corresponding
 452         # to index idx, multiplied by negative one.
 453         bm = (-matrix(self.constraint_matrix).transpose()).tolist()
 454         return bm[idx]
 455
 456     def basis_matrix(self):
 457         # Return the columns of our constraint matrix corresponding
 458         # to the basic variables.
 459         idxs = self.nonbasic_indices
 460         bm = matrix(self.constraint_matrix)
 461         bm = delete(bm, idxs, axis=1)
 462         return bm.tolist()
 463
 464     def reduced_cost(self, idx):
 465         # Find the reduced cost ofthe variable whose column has index
 466         # idx.
 467         dx = array(self.delta_x(idx))
 468         c = array(self.objective_coefficients)
 469         return dot(dx, c)
 470
 471     def delta_x(self, idx):
 472         # Return the array of deltas that would result from increasing
 473         # the variable with index idx by one.
 474         A = matrix(self.basis_matrix())
 475         b = array(self.negative_constraint_column(idx))
 476         x = linalg.solve(A,b).tolist()
 477
 478         # Now we get to add back the nonbasic columns. The 'idx' column
 479         # should have a value of 1, and the other two non-basics should
 480         # have a value of 0.
 481         xfull = [0.0] * len(self.solution_vector)
 482         bidxs = self.basic_indices
 483
 484         for bidx in bidxs:
 485             xfull[bidx] = x.pop(0)
 486
 487         xfull[idx] = 1.0
 488
 489         return xfull