"""
Classes to create, solve, and make dinner for linear programs. Handles
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_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 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

# "Epsilon levels"
# Used to control rounding thresholds.
EPS_TIGHT = 0
EPS_MEDIUM = 1
EPS_LOOSE = 2
EPS_BAGGY = 3

class LinearProgram(object):
    """
    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 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):
        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):
        """
        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
        self._eps_level = EPS_MEDIUM


    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
        self.eps_level = self.eps_level



    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):
        """
        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 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