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