]>
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 class LinearProgram(object):
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.
38 def get_row_count(self
):
40 Return the number of rows in the constraint matrix.
42 return len(self
.constraint_matrix
)
45 def get_column_count(self
):
47 Return the number of columns in the constraint matrix.
48 If we don't have any rows yet, claim zero columns as well.
50 if self
.get_row_count() == 0:
53 return len(self
.constraint_matrix
[0])
60 A property representing the type of linear program, either
70 lpsolve('set_minim', self
._lp
)
74 lpsolve('set_maxim', self
._lp
)
79 def objective_coefficients(self
):
80 return self
._objective
_coefficients
83 @objective_coefficients.setter
84 def objective_coefficients(self
, value
):
85 self
._objective
_coefficients
= value
90 self
._objective
_coefficients
)
95 def constraint_matrix(self
):
96 return self
._constraint
_matrix
98 @constraint_matrix.setter
99 def constraint_matrix(self
, value
):
100 self
._constraint
_matrix
= value
103 lpsolve('set_mat', self
._lp
, value
)
112 def rhs(self
, value
):
116 lpsolve('set_rh_vec', self
._lp
, self
._rhs
)
121 def inequalities(self
):
122 return self
._inequalities
125 def inequalities(self
, value
):
126 self
._inequalities
= value
129 for i
in range(self
.get_row_count()):
130 lpsolve('set_constr_type', self
._lp
, i
+1, value
[i
])
134 def solution_lower_bounds(self
):
135 return self
._solution
_lower
_bounds
137 @solution_lower_bounds.setter
138 def solution_lower_bounds(self
, value
):
139 if len(value
) != self
.get_column_count():
142 self
._solution
_lower
_bounds
= value
145 for i
in range(self
.get_column_count()):
146 lpsolve('set_lowbo', self
._lp
, i
+1, value
[i
])
151 def solution_upper_bounds(self
):
152 return self
._solution
_upper
_bounds
155 @solution_upper_bounds.setter
156 def solution_upper_bounds(self
, value
):
157 if len(value
) != self
.get_column_count():
160 self
._solution
_upper
_bounds
= value
163 for i
in range(self
.get_column_count()):
164 lpsolve('set_upbo', self
._lp
, i
+1, value
[i
])
169 def integer_variables(self
):
171 A vector containing the indices of any solution variables
172 which must be integers.
174 return self
._integer
_variables
176 @integer_variables.setter
177 def integer_variables(self
, value
):
178 self
._integer
_variables
= value
181 for i
in range(len(value
)):
182 lpsolve('set_int', self
._lp
, value
[i
], 1)
186 def make_all_variables_integers(self
):
188 Force all solution variables to be integers. This is achieved
189 by filling the integer_variables vector with all possible
194 for i
in range(self
.get_column_count()):
197 lpsolve('set_int', self
._lp
, i
+1, 1)
199 self
.integer_variables
= ivs
204 def scale_mode(self
):
206 The scaling mode used for handling floating point numbers.
207 See <http://lpsolve.sourceforge.net/5.5/scaling.htm> for more
210 return self
._scale
_mode
214 def scale_mode(self
, value
):
215 self
._scale
_mode
= value
218 lpsolve('set_scaling', self
._lp
, value
)
222 def print_tableau(self
):
224 Tell lp_solve to print its simplex tableau. Only works after
225 a successful call to solve().
227 lpsolve('set_outputfile', self
._lp
, '')
228 lpsolve('print_tableau', self
._lp
)
233 Initialize the object, setting all of the properties
234 either empty or to sane defaults.
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.
242 If the _lp variable is *not* None, the property setters will
243 make calls to lp_solve, updating the pre-existing linear program.
247 self
._objective
_coefficients
= []
248 self
._constraint
_matrix
= []
250 self
._inequalities
= []
251 self
._integer
_variables
= []
252 self
._solution
_lower
_bounds
= []
253 self
._solution
_upper
_bounds
= []
255 self
._type
= MINIMIZE
258 def set_all_lp_properties(self
):
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.
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.
271 self
.constraint_matrix
= self
.constraint_matrix
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
285 lpsolve('delete_lp', self
._lp
)
289 def create_lp_if_necessary(self
):
291 If we already have a linear program instance, do nothing.
292 Otherwise, create one, and set all of the necessary properties.
297 self
._lp
= lpsolve('make_lp',
298 self
.get_row_count(),
299 self
.get_column_count())
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
)
305 self
.set_all_lp_properties()
311 Solve the linear program. The lp_solve instance is
312 created beforehand if necessary.
314 self
.create_lp_if_necessary()
315 result
= lpsolve('solve', self
._lp
)
317 # Default to empty return values.
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
):
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
)
333 return [obj
, x
, duals
]
336 def objective_coefficient_gcd(self
):
338 Return the GCD of all objective function coefficients.
340 return reduce(fractions
.gcd
, self
.objective_coefficients
)
344 class SimplexIteration(object):
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.
353 def constraint_matrix(self
):
354 return self
._constraint
_matrix
.tolist()
356 @constraint_matrix.setter
357 def constraint_matrix(self
, value
):
358 self
._constraint
_matrix
= matrix(value
)
365 def rhs(self
, value
):
369 def objective_coefficients(self
):
370 return self
._objective
_coefficients
372 @objective_coefficients.setter
373 def objective_coefficients(self
, value
):
374 self
._objective
_coefficients
= value
377 def solution_vector(self
):
378 return self
._solution
_vector
.tolist()
380 @solution_vector.setter
381 def solution_vector(self
, value
):
382 self
._solution
_vector
= array(value
)
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
)
393 @basic_variables.setter
394 def basic_variables(self
, value
):
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.
401 vars = value
.split(',')
405 var
= var
.replace('x', '')
406 basic_indices
.append(int(var
)-1)
408 self
.basic_indices
= basic_indices
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
)
419 def basic_indices(self
):
420 return self
._basic
_indices
423 @basic_indices.setter
424 def basic_indices(self
, value
):
425 self
._basic
_indices
= value
429 def nonbasic_indices(self
):
430 all_indices
= range(0, len(self
.solution_vector
))
431 return list(set(all_indices
) - set(self
.basic_indices
))
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:
444 def constraint_column(self
, idx
):
445 # Return the column of the constraint matrix corresponding
447 bm
= matrix(self
.constraint_matrix
).transpose().tolist()
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()
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)
464 def reduced_cost(self
, idx
):
465 # Find the reduced cost ofthe variable whose column has index
467 dx
= array(self
.delta_x(idx
))
468 c
= array(self
.objective_coefficients
)
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()
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
481 xfull
= [0.0] * len(self
.solution_vector
)
482 bidxs
= self
.basic_indices
485 xfull
[bidx
] = x
.pop(0)