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 = 0
MAXIMIZE = 1
+# "Epsilon levels"
+# Used to control rounding thresholds.
+EPS_TIGHT = 0
+EPS_MEDIUM = 1
+EPS_LOOSE = 2
+EPS_BAGGY = 3
class LinearProgram(object):
"""
"""
return self._type
-
+ @type.setter
def type(self, type):
if type == MINIMIZE:
self._type = MINIMIZE
lpsolve('set_rh_vec', self._lp, self._rhs)
+ @property
+ def eps_level(self):
+ return self._eps_level
+
+ @eps_level.setter
+ def eps_level(self, value):
+ self._eps_level = value
+
+ if self._lp != None:
+ lpsolve('set_epslevel', self._lp, value)
+
@property
def inequalities(self):
+ 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):
"""
Initialize the object, setting all of the properties
self._solution_upper_bounds = []
self._scale_mode = 0
self._type = MINIMIZE
+ self._eps_level = EPS_MEDIUM
def set_all_lp_properties(self):
self.solution_upper_bounds = self.solution_upper_bounds
self.scale_mode = self.scale_mode
self.type = self.type
+ self.eps_level = self.eps_level
[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 variables(self):
+ vars = []
+ for idx in range(1, len(self.solution_vector)+1):
+ vars.append("x" + str(idx))
+
+ return vars
+
+
+ @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 objective_function_value(self):
+ c = array(self.objective_coefficients)
+ sv = array(self.solution_vector)
+ return dot(sv, c)
+
+
+ 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