integration with lp_solve.
"""
+import fractions
+from numpy import *
import os
import site
import sys
# Add LP_SOLVE_PATH to our path. There is no point to this variable
# other than to make the site.addsitedir() line fit within 80
# characters.
-LP_SOLVE_PATH = '/../lib/lp_solve'
-site.addsitedir(os.path.dirname(os.path.abspath(sys.argv[0])) + LP_SOLVE_PATH)
+LP_SOLVE_PATH1 = '/../../lib/lp_solve'
+LP_SOLVE_PATH2 = '/../lib/lp_solve'
+site.addsitedir(os.path.dirname(os.path.abspath(sys.argv[0])) + LP_SOLVE_PATH1)
+site.addsitedir(os.path.dirname(os.path.abspath(sys.argv[0])) + LP_SOLVE_PATH2)
-from lp_solve import *
+from lpsolve55 import *
+
+
+# Constants representing the two types of linear programs.
+# MINIMIZE means that we would like to minimize the objective
+# function, and MAXIMIZE means that we would like to maximize it.
+MINIMIZE = 0
+MAXIMIZE = 1
-MINIMUM = 0
-MAXIMUM = 1
class LinearProgram(object):
"""
- Represents an instance of an lp_solve linear program.
+ Represents an instance of an lp_solve linear program.
+ The actual lp_solve linear program is only created when it
+ is needed, and modifications to it are cached beforehand.
"""
-
+
+
+ def get_row_count(self):
+ """
+ Return the number of rows in the constraint matrix.
+ """
+ return len(self.constraint_matrix)
+
+
+ def get_column_count(self):
+ """
+ Return the number of columns in the constraint matrix.
+ If we don't have any rows yet, claim zero columns as well.
+ """
+ if self.get_row_count() == 0:
+ return 0
+ else:
+ return len(self.constraint_matrix[0])
+
+
+
+ @property
+ def type(self):
+ """
+ A property representing the type of linear program, either
+ MINIMIZE or MAXIMIZE.
+ """
+ return self._type
+
+ @type.setter
+ def type(self, type):
+ if type == MINIMIZE:
+ self._type = MINIMIZE
+ if self._lp != None:
+ lpsolve('set_minim', self._lp)
+ else:
+ self._type = MAXIMIZE
+ if self._lp != None:
+ lpsolve('set_maxim', self._lp)
+
+
+
+ @property
+ def objective_coefficients(self):
+ return self._objective_coefficients
+
+
+ @objective_coefficients.setter
+ def objective_coefficients(self, value):
+ self._objective_coefficients = value
+
+ if self._lp != None:
+ lpsolve('set_obj_fn',
+ self._lp,
+ self._objective_coefficients)
+
+
+
+ @property
+ def constraint_matrix(self):
+ return self._constraint_matrix
+
+ @constraint_matrix.setter
+ def constraint_matrix(self, value):
+ self._constraint_matrix = value
+
+ if self._lp != None:
+ lpsolve('set_mat', self._lp, value)
+
+
+
+ @property
+ def rhs(self):
+ return self._rhs
+
+ @rhs.setter
+ def rhs(self, value):
+ self._rhs = value
+
+ if self._lp != None:
+ lpsolve('set_rh_vec', self._lp, self._rhs)
+
+
+
+ @property
+ def inequalities(self):
+ return self._inequalities
+
+ @inequalities.setter
+ def inequalities(self, value):
+ self._inequalities = value
+
+ if self._lp != None:
+ for i in range(self.get_row_count()):
+ lpsolve('set_constr_type', self._lp, i+1, value[i])
+
+
+ @property
+ def solution_lower_bounds(self):
+ return self._solution_lower_bounds
+
+ @solution_lower_bounds.setter
+ def solution_lower_bounds(self, value):
+ if len(value) != self.get_column_count():
+ return
+
+ self._solution_lower_bounds = value
+
+ if self._lp != None:
+ for i in range(self.get_column_count()):
+ lpsolve('set_lowbo', self._lp, i+1, value[i])
+
+
+
+ @property
+ def solution_upper_bounds(self):
+ return self._solution_upper_bounds
+
+
+ @solution_upper_bounds.setter
+ def solution_upper_bounds(self, value):
+ if len(value) != self.get_column_count():
+ return
+
+ self._solution_upper_bounds = value
+
+ if self._lp != None:
+ for i in range(self.get_column_count()):
+ lpsolve('set_upbo', self._lp, i+1, value[i])
+
+
+
+ @property
+ def integer_variables(self):
+ """
+ A vector containing the indices of any solution variables
+ which must be integers.
+ """
+ return self._integer_variables
+
+ @integer_variables.setter
+ def integer_variables(self, value):
+ self._integer_variables = value
+
+ if self._lp != None:
+ for i in range(len(value)):
+ lpsolve('set_int', self._lp, value[i], 1)
+
+
+
+ def make_all_variables_integers(self):
+ """
+ Force all solution variables to be integers. This is achieved
+ by filling the integer_variables vector with all possible
+ indices.
+ """
+ ivs = []
+
+ for i in range(self.get_column_count()):
+ ivs.append(i+1)
+ if self._lp != None:
+ lpsolve('set_int', self._lp, i+1, 1)
+
+ self.integer_variables = ivs
+
+
+
+ @property
+ def scale_mode(self):
+ """
+ The scaling mode used for handling floating point numbers.
+ See <http://lpsolve.sourceforge.net/5.5/scaling.htm> for more
+ information.
+ """
+ return self._scale_mode
+
+
+ @scale_mode.setter
+ def scale_mode(self, value):
+ self._scale_mode = value
+
+ if self._lp != None:
+ lpsolve('set_scaling', self._lp, value)
+
+
+
+ def print_tableau(self):
+ """
+ Tell lp_solve to print its simplex tableau. Only works after
+ a successful call to solve().
+ """
+ lpsolve('set_outputfile', self._lp, '')
+ lpsolve('print_tableau', self._lp)
+
+
def __init__(self):
- self.objective_function_coefficients = []
- self.constraint_matrix = []
- self.rhs = []
- self.inequalities = []
+ """
+ Initialize the object, setting all of the properties
+ either empty or to sane defaults.
+
+ The _lp variable is set to None, initially. All of the
+ property setters will test for _lp == None, and will refuse
+ to make calls to lp_solve if that is the case. A new instance
+ of an lp_solve linear program will be created (and stored in
+ the _lp variable) on demand.
+
+ If the _lp variable is *not* None, the property setters will
+ make calls to lp_solve, updating the pre-existing linear program.
+ """
+
+ self._lp = None
+ self._objective_coefficients = []
+ self._constraint_matrix = []
+ self._rhs = []
+ self._inequalities = []
+ self._integer_variables = []
+ self._solution_lower_bounds = []
+ self._solution_upper_bounds = []
+ self._scale_mode = 0
+ self._type = MINIMIZE
+
+
+ def set_all_lp_properties(self):
+ """
+ Re-set all linear program properties. After a new linear
+ program is created, it will be 'empty'. We already have
+ its properties stored in our member variables, however,
+ we need to make calls to lp_solve to set them on the new
+ linear program instance.
+
+ All of the property setters will check for the existence of
+ self._lp and make calls to lp_solve as necessary. So, to set
+ all of our properties, we just have to trigger the property
+ setters a second time.
+ """
+ self.constraint_matrix = self.constraint_matrix
+ self.rhs = self.rhs
+ self.objective_coefficients = self.objective_coefficients
+ self.inequalities = self.inequalities
+ self.integer_variables = self.integer_variables
+ self.solution_lower_bounds = self.solution_lower_bounds
+ self.solution_upper_bounds = self.solution_upper_bounds
+ self.scale_mode = self.scale_mode
+ self.type = self.type
+
+
+
+ def delete(self):
+ if self._lp != None:
+ lpsolve('delete_lp', self._lp)
+
+
+
+ def create_lp_if_necessary(self):
+ """
+ If we already have a linear program instance, do nothing.
+ Otherwise, create one, and set all of the necessary properties.
+ """
+ if self._lp != None:
+ return
+
+ self._lp = lpsolve('make_lp',
+ self.get_row_count(),
+ self.get_column_count())
+
+ # This is not critical, but it will encourage lp_solve to
+ # warn us about potential problems.
+ lpsolve('set_verbose', self._lp, IMPORTANT)
+
+ self.set_all_lp_properties()
+
def solve(self):
- [v,x,duals] = lp_solve(self.objective_function_coefficients,
- self.constraint_matrix,
- self.rhs,
- self.inequalities)
- return [v,x,duals]
+ """
+ Solve the linear program. The lp_solve instance is
+ created beforehand if necessary.
+ """
+ self.create_lp_if_necessary()
+ result = lpsolve('solve', self._lp)
+
+ # Default to empty return values.
+ obj = []
+ x = []
+ duals = []
+
+ # See http://lpsolve.sourceforge.net/5.5/solve.htm for a
+ # description of these constants.
+ if (result == OPTIMAL or
+ result == SUBOPTIMAL or
+ result == PROCBREAK or
+ result == FEASFOUND):
+
+ # If the result was "good," i.e. if get_solution will work,
+ # call it and use its return value as ours.
+ [obj, x, duals, ret] = lpsolve('get_solution', self._lp)
+
+ return [obj, x, duals]
+
+
+ def objective_coefficient_gcd(self):
+ """
+ Return the GCD of all objective function coefficients.
+ """
+ return reduce(fractions.gcd, self.objective_coefficients)
+
+
+
+class SimplexIteration(object):
+ """
+ Represents the 'current' iteration of the simplex method at some
+ point. It needs an A,b,x, and c corresponding to the linear
+ program in standard form. It can then determine whether or not the
+ current vertex (x) is optimal, and possible move to a better one.
+ """
+
+ @property
+ def constraint_matrix(self):
+ return self._constraint_matrix.tolist()
+
+ @constraint_matrix.setter
+ def constraint_matrix(self, value):
+ self._constraint_matrix = matrix(value)
+
+ @property
+ def rhs(self):
+ return self._rhs
+
+ @rhs.setter
+ def rhs(self, value):
+ self._rhs = value
+
+ @property
+ def objective_coefficients(self):
+ return self._objective_coefficients
+
+ @objective_coefficients.setter
+ def objective_coefficients(self, value):
+ self._objective_coefficients = value
+
+ @property
+ def solution_vector(self):
+ return self._solution_vector.tolist()
+
+ @solution_vector.setter
+ def solution_vector(self, value):
+ self._solution_vector = array(value)
+
+
+ @property
+ def basic_variables(self):
+ # The current set of basic variables. Constructed from the
+ # "true" source, our list of basic indices.
+ idxs = self.basic_indices
+ return map(lambda x: "x" + str(x+1), idxs)
+
+
+ @basic_variables.setter
+ def basic_variables(self, value):
+ """
+ Syntactic sugar to set the basic indices. We take a string
+ of the form x1,x2,...xn, and subtract one from each of the
+ subscripts to get the basic indices.
+ """
+ basic_indices = []
+ vars = value.split(',')
+
+ for var in vars:
+ var = var.strip()
+ var = var.replace('x', '')
+ basic_indices.append(int(var)-1)
+
+ self.basic_indices = basic_indices
+
+
+ @property
+ def nonbasic_variables(self):
+ # All elements of the solution vector that have value zero.
+ idxs = self.nonbasic_indices
+ return map(lambda x: "x" + str(x+1), idxs)
+
+
+ @property
+ def basic_indices(self):
+ return self._basic_indices
+
+
+ @basic_indices.setter
+ def basic_indices(self, value):
+ self._basic_indices = value
+
+
+ @property
+ def nonbasic_indices(self):
+ all_indices = range(0, len(self.solution_vector))
+ return list(set(all_indices) - set(self.basic_indices))
+
+
+ @property
+ def optimal(self):
+ # True if the current solution is optimal, false if not.
+ for idx in self.nonbasic_indices:
+ if self.reduced_cost(idx) < 0.0:
+ return False
+
+ return True
+
+
+ def constraint_column(self, idx):
+ # Return the column of the constraint matrix corresponding
+ # to index idx.
+ bm = matrix(self.constraint_matrix).transpose().tolist()
+ return bm[idx]
+
+ def negative_constraint_column(self, idx):
+ # Return the column of the constraint matrix corresponding
+ # to index idx, multiplied by negative one.
+ bm = (-matrix(self.constraint_matrix).transpose()).tolist()
+ return bm[idx]
+
+ def basis_matrix(self):
+ # Return the columns of our constraint matrix corresponding
+ # to the basic variables.
+ idxs = self.nonbasic_indices
+ bm = matrix(self.constraint_matrix)
+ bm = delete(bm, idxs, axis=1)
+ return bm.tolist()
+
+ def reduced_cost(self, idx):
+ # Find the reduced cost ofthe variable whose column has index
+ # idx.
+ dx = array(self.delta_x(idx))
+ c = array(self.objective_coefficients)
+ return dot(dx, c)
+
+ def delta_x(self, idx):
+ # Return the array of deltas that would result from increasing
+ # the variable with index idx by one.
+ A = matrix(self.basis_matrix())
+ b = array(self.negative_constraint_column(idx))
+ x = linalg.solve(A,b).tolist()
+
+ # Now we get to add back the nonbasic columns. The 'idx' column
+ # should have a value of 1, and the other two non-basics should
+ # have a value of 0.
+ xfull = [0.0] * len(self.solution_vector)
+ bidxs = self.basic_indices
+
+ for bidx in bidxs:
+ xfull[bidx] = x.pop(0)
+
+ xfull[idx] = 1.0
+ return xfull
projects / dead / census-tools.git / blobdiff