]>
gitweb.michael.orlitzky.com - dead/census-tools.git/blob - src/LinearProgramming.py
2 Classes to create, solve, and make dinner for linear programs. Handles
3 integration with lp_solve.
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
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
)
20 from lpsolve55
import *
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.
30 # Used to control rounding thresholds.
36 class LinearProgram(object):
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.
44 def get_row_count(self
):
46 Return the number of rows in the constraint matrix.
48 return len(self
.constraint_matrix
)
51 def get_column_count(self
):
53 Return the number of columns in the constraint matrix.
54 If we don't have any rows yet, claim zero columns as well.
56 if self
.get_row_count() == 0:
59 return len(self
.constraint_matrix
[0])
66 A property representing the type of linear program, either
76 lpsolve('set_minim', self
._lp
)
80 lpsolve('set_maxim', self
._lp
)
85 def objective_coefficients(self
):
86 return self
._objective
_coefficients
89 @objective_coefficients.setter
90 def objective_coefficients(self
, value
):
91 self
._objective
_coefficients
= value
96 self
._objective
_coefficients
)
101 def constraint_matrix(self
):
102 return self
._constraint
_matrix
104 @constraint_matrix.setter
105 def constraint_matrix(self
, value
):
106 self
._constraint
_matrix
= value
109 lpsolve('set_mat', self
._lp
, value
)
118 def rhs(self
, value
):
122 lpsolve('set_rh_vec', self
._lp
, self
._rhs
)
127 return self
._eps
_level
130 def eps_level(self
, value
):
131 self
._eps
_level
= value
134 lpsolve('set_epslevel', self
._lp
, value
)
138 def inequalities(self
):
139 return self
._inequalities
142 def inequalities(self
, value
):
143 self
._inequalities
= value
146 for i
in range(self
.get_row_count()):
147 lpsolve('set_constr_type', self
._lp
, i
+1, value
[i
])
151 def solution_lower_bounds(self
):
152 return self
._solution
_lower
_bounds
154 @solution_lower_bounds.setter
155 def solution_lower_bounds(self
, value
):
156 if len(value
) != self
.get_column_count():
159 self
._solution
_lower
_bounds
= value
162 for i
in range(self
.get_column_count()):
163 lpsolve('set_lowbo', self
._lp
, i
+1, value
[i
])
168 def solution_upper_bounds(self
):
169 return self
._solution
_upper
_bounds
172 @solution_upper_bounds.setter
173 def solution_upper_bounds(self
, value
):
174 if len(value
) != self
.get_column_count():
177 self
._solution
_upper
_bounds
= value
180 for i
in range(self
.get_column_count()):
181 lpsolve('set_upbo', self
._lp
, i
+1, value
[i
])
186 def integer_variables(self
):
188 A vector containing the indices of any solution variables
189 which must be integers.
191 return self
._integer
_variables
193 @integer_variables.setter
194 def integer_variables(self
, value
):
195 self
._integer
_variables
= value
198 for i
in range(len(value
)):
199 lpsolve('set_int', self
._lp
, value
[i
], 1)
203 def make_all_variables_integers(self
):
205 Force all solution variables to be integers. This is achieved
206 by filling the integer_variables vector with all possible
211 for i
in range(self
.get_column_count()):
214 lpsolve('set_int', self
._lp
, i
+1, 1)
216 self
.integer_variables
= ivs
221 def scale_mode(self
):
223 The scaling mode used for handling floating point numbers.
224 See <http://lpsolve.sourceforge.net/5.5/scaling.htm> for more
227 return self
._scale
_mode
231 def scale_mode(self
, value
):
232 self
._scale
_mode
= value
235 lpsolve('set_scaling', self
._lp
, value
)
239 def print_tableau(self
):
241 Tell lp_solve to print its simplex tableau. Only works after
242 a successful call to solve().
244 lpsolve('set_outputfile', self
._lp
, '')
245 lpsolve('print_tableau', self
._lp
)
250 Initialize the object, setting all of the properties
251 either empty or to sane defaults.
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.
259 If the _lp variable is *not* None, the property setters will
260 make calls to lp_solve, updating the pre-existing linear program.
264 self
._objective
_coefficients
= []
265 self
._constraint
_matrix
= []
267 self
._inequalities
= []
268 self
._integer
_variables
= []
269 self
._solution
_lower
_bounds
= []
270 self
._solution
_upper
_bounds
= []
272 self
._type
= MINIMIZE
273 self
._eps
_level
= EPS_MEDIUM
276 def set_all_lp_properties(self
):
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.
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.
289 self
.constraint_matrix
= self
.constraint_matrix
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
304 lpsolve('delete_lp', self
._lp
)
308 def create_lp_if_necessary(self
):
310 If we already have a linear program instance, do nothing.
311 Otherwise, create one, and set all of the necessary properties.
316 self
._lp
= lpsolve('make_lp',
317 self
.get_row_count(),
318 self
.get_column_count())
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
)
324 self
.set_all_lp_properties()
330 Solve the linear program. The lp_solve instance is
331 created beforehand if necessary.
333 self
.create_lp_if_necessary()
334 result
= lpsolve('solve', self
._lp
)
336 # Default to empty return values.
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
):
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
)
352 return [obj
, x
, duals
]
355 def objective_coefficient_gcd(self
):
357 Return the GCD of all objective function coefficients.
359 return reduce(fractions
.gcd
, self
.objective_coefficients
)
363 class SimplexIteration(object):
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.
372 def constraint_matrix(self
):
373 return self
._constraint
_matrix
.tolist()
375 @constraint_matrix.setter
376 def constraint_matrix(self
, value
):
377 self
._constraint
_matrix
= matrix(value
)
384 def rhs(self
, value
):
388 def objective_coefficients(self
):
389 return self
._objective
_coefficients
391 @objective_coefficients.setter
392 def objective_coefficients(self
, value
):
393 self
._objective
_coefficients
= value
396 def solution_vector(self
):
397 return self
._solution
_vector
.tolist()
399 @solution_vector.setter
400 def solution_vector(self
, value
):
401 self
._solution
_vector
= array(value
)
407 for idx
in range(1, len(self
.solution_vector
)+1):
408 vars.append("x" + str(idx
))
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
)
421 @basic_variables.setter
422 def basic_variables(self
, value
):
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.
429 vars = value
.split(',')
433 var
= var
.replace('x', '')
434 basic_indices
.append(int(var
)-1)
436 self
.basic_indices
= basic_indices
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
)
447 def basic_indices(self
):
448 return self
._basic
_indices
451 @basic_indices.setter
452 def basic_indices(self
, value
):
453 self
._basic
_indices
= value
457 def nonbasic_indices(self
):
458 all_indices
= range(0, len(self
.solution_vector
))
459 return list(set(all_indices
) - set(self
.basic_indices
))
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:
472 def constraint_column(self
, idx
):
473 # Return the column of the constraint matrix corresponding
475 bm
= matrix(self
.constraint_matrix
).transpose().tolist()
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()
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)
493 def objective_function_value(self
):
494 c
= array(self
.objective_coefficients
)
495 sv
= array(self
.solution_vector
)
499 def reduced_cost(self
, idx
):
500 # Find the reduced cost ofthe variable whose column has index
502 dx
= array(self
.delta_x(idx
))
503 c
= array(self
.objective_coefficients
)
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()
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
516 xfull
= [0.0] * len(self
.solution_vector
)
517 bidxs
= self
.basic_indices
520 xfull
[bidx
] = x
.pop(0)