Solvers.h

Go to the documentation of this file.

#ifndef LIMBO_SOLVERS_SOLVERS_H
#define LIMBO_SOLVERS_SOLVERS_H

#include <iostream>
#include <limits>
#include <string>
#include <vector>
#include <map>
#include <algorithm>
#include <limbo/preprocessor/Msg.h>
#include <limbo/math/Math.h>
#include <limbo/parsers/lp/bison/LpDriver.h>

namespace limbo 
{
namespace solvers 
{

enum SolverProperty
{
    MIN, 
    MAX, 
    BINARY, 
    INTEGER, 
    CONTINUOUS, 
    OPTIMAL, 
    INFEASIBLE, 
    SUBOPTIMAL, 
    UNBOUNDED 
};

inline std::string toString(SolverProperty sp)
{
    switch (sp)
    {
        case MIN:
            return "MIN";
        case MAX:
            return "MAX";
        case BINARY:
            return "BINARY";
        case INTEGER:
            return "INTEGER";
        case CONTINUOUS:
            return "CONTINUOUS";
        case OPTIMAL:
            return "OPTIMAL";
        case INFEASIBLE:
            return "INFEASIBLE";
        case SUBOPTIMAL:
        default:
            return "SUBOPTIMAL";
    }
}

// forward declaration 
template <typename T>
class Variable;
template <typename T>
class VariableProperty;
template <typename T>
class LinearTerm; 
template <typename T>
class LinearExpression;
template <typename T>
class LinearConstraint;
template <typename T, typename V>
class LinearModel; 
template <typename T, typename I, int StartingIndex>
struct MatrixCSR; 

template <typename T>
class Variable 
{
    public:
        typedef T coefficient_value_type;

        explicit Variable(unsigned int id = std::numeric_limits<unsigned int>::max()) 
            : m_id(id) 
        {
        }
        Variable(Variable const& rhs)
        {
            copy(rhs);
        }
        Variable& operator=(Variable const& rhs)
        {
            if (this != &rhs)
                copy(rhs);
            return *this; 
        }
        ~Variable() {}

        unsigned int id() const {return m_id;}
        void setId(unsigned int id) {m_id = id;}

        bool valid() const {return m_id != std::numeric_limits<unsigned int>::max();}

        bool operator==(Variable const& rhs) const 
        {
            return m_id == rhs.m_id;
        }
        bool operator!=(Variable const& rhs) const 
        {
            return !this->operator==(rhs);
        }
        LinearTerm<coefficient_value_type> operator-() const 
        {
            LinearTerm<coefficient_value_type> term (*this); 
            return term.negate();
        }
        friend LinearTerm<coefficient_value_type> operator*(Variable const& var, coefficient_value_type coef)
        {
            return LinearTerm<coefficient_value_type>(var, coef); 
        }
        friend LinearTerm<coefficient_value_type> operator*(coefficient_value_type coef, Variable const& var)
        {
            return var*coef; 
        }
        friend LinearTerm<coefficient_value_type> operator/(Variable const& var, coefficient_value_type coef)
        {
            return LinearTerm<coefficient_value_type>(var)/coef; 
        }
        friend LinearExpression<coefficient_value_type> operator+(Variable const& var, coefficient_value_type constant)
        {
            return LinearTerm<coefficient_value_type>(var)+constant; 
        }
        friend LinearExpression<coefficient_value_type> operator+(coefficient_value_type constant, Variable const& var)
        {
            return var+constant;
        }
        friend LinearExpression<coefficient_value_type> operator+(Variable const& var1, Variable const& var2)
        {
            return LinearTerm<coefficient_value_type>(var1)+var2; 
        }
        friend LinearExpression<coefficient_value_type> operator-(Variable const& var, coefficient_value_type constant)
        {
            return LinearTerm<coefficient_value_type>(var)-constant; 
        }
        friend LinearExpression<coefficient_value_type> operator-(coefficient_value_type constant, Variable const& var)
        {
            return LinearTerm<coefficient_value_type>(var, -1)+constant; 
        }
        friend LinearExpression<coefficient_value_type> operator-(Variable const& var1, Variable const& var2)
        {
            return LinearTerm<coefficient_value_type>(var1)-var2; 
        }
        friend LinearConstraint<coefficient_value_type> operator<(Variable const& var, coefficient_value_type rhs)
        {
            return (LinearTerm<coefficient_value_type>(var) < rhs);
        }
        friend LinearConstraint<coefficient_value_type> operator<=(Variable const& var, coefficient_value_type rhs)
        {
            return (var < rhs);
        }
        friend LinearConstraint<coefficient_value_type> operator>(Variable const& var, coefficient_value_type rhs)
        {
            return (LinearTerm<coefficient_value_type>(var) > rhs);
        }
        friend LinearConstraint<coefficient_value_type> operator>=(Variable const& var, coefficient_value_type rhs)
        {
            return (var > rhs);
        }
        friend LinearConstraint<coefficient_value_type> operator==(Variable const& var, coefficient_value_type rhs)
        {
            return (LinearTerm<coefficient_value_type>(var) == rhs);
        }
    protected:
        void copy(Variable const& rhs)
        {
            m_id = rhs.m_id; 
        }

        unsigned int m_id; 
};

template <typename V>
class VariableProperty
{
    public:
        typedef V variable_value_type; 

        VariableProperty(variable_value_type lb, variable_value_type ub, SolverProperty nt, std::string const& n)
            : m_lowerBound(lb)
            , m_upperBound(ub)
            , m_numericType(nt)
            , m_name(n)
        {
        }
        VariableProperty(VariableProperty const& rhs)
        {
            copy(rhs);
        }
        VariableProperty& operator=(VariableProperty const& rhs)
        {
            if (this != &rhs)
                copy(rhs);
            return *this;
        }

        variable_value_type lowerBound() const {return m_lowerBound;}
        void setLowerBound(variable_value_type lb) {m_lowerBound = lb;}
        void updateLowerBound(variable_value_type lb) {m_lowerBound = std::max(m_lowerBound, lb);}
        variable_value_type upperBound() const {return m_upperBound;}
        void setUpperBound(variable_value_type ub) {m_upperBound = ub;}
        void updateUpperBound(variable_value_type ub) {m_upperBound = std::min(m_upperBound, ub);}
        SolverProperty numericType() const {return m_numericType;}
        void setNumericType(SolverProperty nt) 
        {
            m_numericType = nt;
            if (m_numericType == BINARY) // update bounds for binary 
            {
                m_lowerBound = 0; 
                m_upperBound = 1; 
            }
        }
        std::string const& name() const {return m_name;}
        void setName(std::string const& n) {m_name = n;}

    protected:
        void copy(VariableProperty const& rhs)
        {
            m_lowerBound = rhs.m_lowerBound;
            m_upperBound = rhs.m_upperBound; 
            m_numericType = rhs.m_numericType; 
            m_name = rhs.m_name;
        }

        variable_value_type m_lowerBound; 
        variable_value_type m_upperBound; 
        SolverProperty m_numericType; 
        std::string m_name; 
};

template <typename T>
class LinearTerm
{
    public:
        typedef T coefficient_value_type;
        typedef Variable<coefficient_value_type> variable_type; 

        LinearTerm(variable_type var = variable_type(), coefficient_value_type coef = 1)
            : m_var(var)
            , m_coef(coef)
        {
        }
        LinearTerm(LinearTerm const& rhs)
        {
            copy(rhs);
        }
        LinearTerm& operator=(LinearTerm const& rhs)
        {
            if (this != &rhs)
                copy(rhs);
            return *this;
        }
        ~LinearTerm()
        {
        }

        variable_type const& variable() const {return m_var;}
        void setVariable(variable_type const& var) {m_var = var;}
        coefficient_value_type coefficient() const {return m_coef;}
        void setCoefficient(coefficient_value_type coef) {m_coef = coef;}

        LinearTerm operator-() const 
        {
            LinearTerm term (*this); 
            return term.negate();
        }
        LinearTerm& negate() 
        {
            m_coef = -m_coef; 
            return *this; 
        }
        friend LinearTerm operator*(LinearTerm const& term, coefficient_value_type coef)
        {
            LinearTerm other (term); 
            other *= coef; 
            return other; 
        }
        friend LinearTerm operator*(coefficient_value_type coef, LinearTerm const& term)
        {
            return term*coef; 
        }
        LinearTerm& operator*=(coefficient_value_type coef)
        {
            m_coef *= coef; 
            return *this;
        }
        friend LinearTerm operator/(LinearTerm const& term, coefficient_value_type coef)
        {
            LinearTerm other (term); 
            other /= coef; 
            return other; 
        }
        LinearTerm& operator/=(coefficient_value_type coef)
        {
            m_coef /= coef; 
            return *this; 
        }
        friend LinearExpression<coefficient_value_type> operator+(LinearTerm const& term, coefficient_value_type constant)
        {
            return LinearExpression<coefficient_value_type>(term, constant);
        }
        friend LinearExpression<coefficient_value_type> operator+(LinearTerm const& term, variable_type const& var)
        {
            return term+LinearTerm(var);
        }
        friend LinearExpression<coefficient_value_type> operator+(coefficient_value_type constant, LinearTerm const& term)
        {
            return LinearExpression<coefficient_value_type>(term, constant);
        }
        friend LinearExpression<coefficient_value_type> operator+(variable_type const& var, LinearTerm const& term)
        {
            return LinearTerm(var)+term;
        }
        friend LinearExpression<coefficient_value_type> operator+(LinearTerm const& term1, LinearTerm const& term2)
        {
            LinearExpression<coefficient_value_type> expr (term1);
            expr += term2; 
            return expr; 
        }
        friend LinearExpression<coefficient_value_type> operator-(LinearTerm const& term, coefficient_value_type constant)
        {
            return LinearExpression<coefficient_value_type>(term, -constant);
        }
        friend LinearExpression<coefficient_value_type> operator-(LinearTerm const& term, variable_type const& var)
        {
            return term-LinearTerm(var);
        }
        friend LinearExpression<coefficient_value_type> operator-(variable_type const& var, LinearTerm const& term)
        {
            return LinearTerm(var)-term;
        }
        friend LinearExpression<coefficient_value_type> operator-(coefficient_value_type constant, LinearTerm const& term)
        {
            return LinearExpression<coefficient_value_type>(-term, constant);
        }
        friend LinearExpression<coefficient_value_type> operator-(LinearTerm const& term1, LinearTerm const& term2)
        {
            LinearExpression<coefficient_value_type> expr (term1);
            expr -= term2; 
            return expr; 
        }
        friend LinearConstraint<coefficient_value_type> operator<(LinearTerm const& term, coefficient_value_type rhs)
        {
            return LinearConstraint<coefficient_value_type>(term, rhs, '<');
        }
        friend LinearConstraint<coefficient_value_type> operator<=(LinearTerm const& term, coefficient_value_type rhs)
        {
            return (term < rhs);
        }
        friend LinearConstraint<coefficient_value_type> operator>(LinearTerm const& term, coefficient_value_type rhs)
        {
            return LinearConstraint<coefficient_value_type>(term, rhs, '>');
        }
        friend LinearConstraint<coefficient_value_type> operator>=(LinearTerm const& term, coefficient_value_type rhs)
        {
            return (term > rhs);
        }
        friend LinearConstraint<coefficient_value_type> operator==(LinearTerm const& term, coefficient_value_type rhs)
        {
            return LinearConstraint<coefficient_value_type>(term, rhs, '=');
        }
        bool operator==(LinearTerm const& rhs) const 
        {
            return m_var == rhs.m_var && m_coef == rhs.m_coef;
        }
    protected:
        void copy(LinearTerm const& rhs)
        {
            m_var = rhs.m_var;
            m_coef = rhs.m_coef; 
        }

        variable_type m_var; 
        coefficient_value_type m_coef; 
};

struct CompareTermByVariable
{
    template <typename TermType>
    bool operator()(TermType const& t1, TermType const& t2) const 
    {
        return t1.variable().id() < t2.variable().id();
    }
};

template <typename T>
class LinearExpression 
{
    public:
        typedef T coefficient_value_type;
        typedef LinearTerm<coefficient_value_type> term_type;
        typedef Variable<coefficient_value_type> variable_type; 

        LinearExpression()
            : m_vTerm()
            , m_constant(0)
        {
        }
        LinearExpression(term_type const& term, coefficient_value_type constant = 0)
        {
            m_vTerm.push_back(term);
            m_constant = constant; 
        }
        LinearExpression(LinearExpression const& rhs)
        {
            copy(rhs);
        }
        LinearExpression& operator=(LinearExpression const& rhs)
        {
            if (this != &rhs)
                copy(rhs);
            return *this;
        }
        ~LinearExpression()
        {
        }

        void swap(LinearExpression& rhs)
        {
            m_vTerm.swap(rhs.m_vTerm); 
            std::swap(m_constant, rhs.m_constant);
        }
        std::vector<term_type> const& terms() const {return m_vTerm;}
        void clear() 
        {
            m_vTerm.clear();
            m_constant = 0; 
        }
        void reserve(unsigned int n) {m_vTerm.reserve(n);}
        coefficient_value_type constant() const {return m_constant;}
        LinearExpression& setConstant(coefficient_value_type constant) 
        {
            m_constant = constant; 
            return *this; 
        }
        LinearExpression& incrementConstant(coefficient_value_type constant) 
        {
            m_constant += constant; 
            return *this; 
        }

        friend LinearExpression operator+(LinearExpression const& expr, coefficient_value_type constant)
        {
            LinearExpression other (expr); 
            other += constant; 
            return other; 
        }
        friend LinearExpression operator+(LinearExpression const& expr, term_type const& term)
        {
            LinearExpression other (expr); 
            other += term; 
            return other; 
        }
        friend LinearExpression operator+(coefficient_value_type constant, LinearExpression const& expr)
        {
            return expr+constant; 
        }
        friend LinearExpression operator+(term_type const& term, LinearExpression const& expr)
        {
            return expr+term; 
        }
        friend LinearExpression operator+(LinearExpression const& expr1, LinearExpression const& expr2)
        {
            LinearExpression other (expr1);
            other += expr2; 
            return other; 
        }
        LinearExpression& operator+=(coefficient_value_type constant)
        {
            this->incrementConstant(constant);
            return *this; 
        }
        LinearExpression& operator+=(term_type const& term)
        {
            m_vTerm.push_back(term); 
            return *this; 
        }
        LinearExpression& operator+=(LinearExpression const& expr)
        {
            m_vTerm.insert(m_vTerm.end(), expr.terms().begin(), expr.terms().end()); 
            this->incrementConstant(expr.constant());
            return *this; 
        }
        friend LinearExpression operator-(LinearExpression const& expr, coefficient_value_type constant)
        {
            LinearExpression other (expr); 
            other -= constant; 
            return other; 
        }
        friend LinearExpression operator-(LinearExpression const& expr, term_type const& term)
        {
            LinearExpression other (expr); 
            other -= term; 
            return other; 
        }
        friend LinearExpression operator-(coefficient_value_type constant, LinearExpression const& expr)
        {
            return (-expr)+constant; 
        }
        friend LinearExpression operator-(term_type const& term, LinearExpression const& expr)
        {
            return (-expr)+term; 
        }
        friend LinearExpression operator-(LinearExpression const& expr1, LinearExpression const& expr2)
        {
            LinearExpression other (expr1); 
            other -= expr2; 
            return other; 
        }
        LinearExpression& operator-=(coefficient_value_type constant)
        {
            this->incrementConstant(-constant);
            return *this;
        }
        LinearExpression& operator-=(term_type term)
        {
            term.setCoefficient(-term.coefficient());
            return this->operator+=(term);
        }
        LinearExpression& operator-=(LinearExpression const& expr)
        {
            for (typename std::vector<term_type>::const_iterator it = expr.terms().begin(), ite = expr.terms().end(); it != ite; ++it)
                this->operator-=(*it);
            this->incrementConstant(-expr.constant());
            return *this;
        }
        friend LinearExpression operator*(LinearExpression const& expr, coefficient_value_type c)
        {
            LinearExpression other (expr);
            other *= c; 
            return other; 
        }
        friend LinearExpression operator*(coefficient_value_type c, LinearExpression const& expr)
        {
            return expr*c; 
        }
        LinearExpression& operator*=(coefficient_value_type c)
        {
            for (typename std::vector<term_type>::iterator it = m_vTerm.begin(); it != m_vTerm.end(); ++it)
                it->setCoefficient(c*it->coefficient());
            m_constant *= c; 
            return *this; 
        }
        friend LinearExpression operator/(LinearExpression const& expr, coefficient_value_type c)
        {
            LinearExpression other (expr);
            other /= c; 
            return other; 
        }
        LinearExpression& operator/=(coefficient_value_type c)
        {
            for (typename std::vector<term_type>::iterator it = m_vTerm.begin(); it != m_vTerm.end(); ++it)
                it->setCoefficient(it->coefficient()/c);
            m_constant /= c; 
            return *this; 
        }
        LinearExpression operator-() const 
        {
            LinearExpression expr (*this); 
            return expr.negate();
        }
        LinearExpression& negate() 
        {
            for (typename std::vector<term_type>::iterator it = m_vTerm.begin(); it != m_vTerm.end(); ++it)
                it->setCoefficient(-it->coefficient());
            m_constant = -m_constant; 
            return *this; 
        }
        friend LinearConstraint<coefficient_value_type> operator<(LinearExpression const& expr, coefficient_value_type rhs)
        {
            return LinearConstraint<coefficient_value_type>(expr, rhs, '<');
        }
        friend LinearConstraint<coefficient_value_type> operator<=(LinearExpression const& expr, coefficient_value_type rhs)
        {
            return (expr < rhs);
        }
        friend LinearConstraint<coefficient_value_type> operator>(LinearExpression const& expr, coefficient_value_type rhs)
        {
            return LinearConstraint<coefficient_value_type>(expr, rhs, '>');
        }
        friend LinearConstraint<coefficient_value_type> operator>=(LinearExpression const& expr, coefficient_value_type rhs)
        {
            return (expr > rhs);
        }
        friend LinearConstraint<coefficient_value_type> operator==(LinearExpression const& expr, coefficient_value_type rhs)
        {
            return LinearConstraint<coefficient_value_type>(expr, rhs, '=');
        }

        void simplify()
        {
            if (m_vTerm.empty())
                return;
            std::sort(m_vTerm.begin(), m_vTerm.end(), CompareTermByVariable()); 
            // merge terms 
            typename std::vector<term_type>::iterator itw = m_vTerm.begin(); // write iterator 
            typename std::vector<term_type>::iterator itr = m_vTerm.begin(); // read iterator 
            typename std::vector<term_type>::iterator ite = m_vTerm.end();
            for (; itr != ite; ++itr)
            {
                if (itr != itw)
                {
                    if (itr->variable() == itw->variable())
                        itw->setCoefficient(itw->coefficient()+itr->coefficient());
                    else
                    {
                        ++itw; 
                        *itw = *itr; 
                    }
                }
            }
            m_vTerm.resize(std::distance(m_vTerm.begin(), itw+1));
            // remove terms with zero coefficients 
            // the order is not maintained 
            for (itr = m_vTerm.begin(); itr != m_vTerm.end(); )
            {
                if (itr->coefficient() == 0)
                {
                    *itr = m_vTerm.back(); 
                    m_vTerm.pop_back(); 
                }
                else 
                    ++itr; 
            }
        }

    protected:
        void copy(LinearExpression const& rhs)
        {
            m_vTerm = rhs.m_vTerm;
            m_constant = rhs.m_constant; 
        }

        std::vector<term_type> m_vTerm; 
        coefficient_value_type m_constant; 
};

template <typename T>
class LinearConstraint
{
    public:
        typedef T coefficient_value_type; 
        typedef Variable<coefficient_value_type> variable_type; 
        typedef LinearTerm<coefficient_value_type> term_type;
        typedef LinearExpression<coefficient_value_type> expression_type;

        LinearConstraint(expression_type expr = expression_type(), coefficient_value_type rhs = 0, char s = '<')
            : m_id (std::numeric_limits<unsigned int>::max())
            , m_expr(expr)
            , m_rhs(rhs)
            , m_sense(s)
        {
            clearConstant();
        }
        LinearConstraint(LinearConstraint const& rhs)
        {
            copy(rhs);
        }
        LinearConstraint& operator=(LinearConstraint const& rhs)
        {
            if (this != &rhs)
                copy(rhs);
            return *this; 
        }
        ~LinearConstraint()
        {
        }

        unsigned int id() const {return m_id;}
        void setId(unsigned int i) {m_id = i;}
        expression_type const& expression() const {return m_expr;}
        void setExpression(expression_type const& expr) {m_expr = expr;}
        void emplaceExpression(expression_type& expr) {m_expr.swap(expr);}
        coefficient_value_type rightHandSide() const {return m_rhs;}
        void setRightHandSide(coefficient_value_type rhs) {m_rhs = rhs;}
        char sense() const {return m_sense;}
        void setSense(char s) {m_sense = s;}
        void swap(LinearConstraint& rhs)
        {
            std::swap(m_id, rhs.m_id);
            m_expr.swap(rhs.m_expr); 
            std::swap(m_rhs, rhs.m_rhs);
            std::swap(m_sense, rhs.m_sense); 
#ifdef DEBUG_SOLVERS
            limboAssertMsg(m_expr.constant() == 0, "expression constant should always be zero in constraint");
#endif
        }
        void simplify()
        {
#ifdef DEBUG_SOLVERS
            limboAssertMsg(m_expr.constant() == 0, "expression constant should always be zero in constraint");
#endif
            m_expr.simplify();
        }
        void normalize(char s)
        {
            simplify();
            if ((m_sense == '<' && s == '>')
                    || (m_sense == '>' && s == '<'))
            {
                m_expr.negate();
                m_rhs = -m_rhs; 
                m_sense = s; 
            }
        }
        void scale(coefficient_value_type factor)
        {
#ifdef DEBUG_SOLVERS
            limboAssertMsg(m_expr.constant() == 0, "expression constant should always be zero in constraint");
#endif
            // in case factor is negative 
            if (factor < 0)
            {
                if (m_sense == '<')
                    m_sense = '>';
                else if (m_sense == '>')
                    m_sense = '<';
            }
            m_expr *= factor; 
            m_rhs *= factor; 
        }
        void clearConstant()
        {
            // make sure the constant is zero in expression 
            m_rhs -= m_expr.constant();
            m_expr.setConstant(0);
        }
        void reserve(unsigned int n) {m_expr.reserve(n);}

        LinearConstraint& operator+=(term_type const& term)
        {
            m_expr += term; 
            return *this; 
        }
        LinearConstraint& operator+=(expression_type const& expr)
        {
            m_expr += expr; 
            clearConstant();
            return *this; 
        }
        LinearConstraint& operator-=(term_type term)
        {
            m_expr -= term; 
            return *this;
        }
        LinearConstraint& operator-=(expression_type const& expr)
        {
            m_expr -= expr; 
            clearConstant();
            return *this;
        }
    protected:
        void copy(LinearConstraint const& rhs)
        {
            m_id = rhs.m_id; 
            m_expr = rhs.m_expr;
            m_rhs = rhs.m_rhs;
            m_sense = rhs.m_sense;
            clearConstant();
        }

        unsigned int m_id; 
        expression_type m_expr; 
        coefficient_value_type m_rhs; 
        char m_sense; 
};

template <typename T, typename V>
class LinearModel : public LpParser::LpDataBase
{
    public:
        typedef T coefficient_value_type; 
        typedef V variable_value_type; 
        typedef LpParser::LpDataBase base_type; 
        typedef Variable<coefficient_value_type> variable_type; 
        typedef LinearTerm<coefficient_value_type> term_type; 
        typedef LinearExpression<coefficient_value_type> expression_type; 
        typedef LinearConstraint<coefficient_value_type> constraint_type; 
        typedef VariableProperty<variable_value_type> property_type; 

        LinearModel() 
            : base_type()
        {
            m_optType = MIN; 
        }
        LinearModel(LinearModel const& rhs) 
            : base_type(rhs)
        {
            copy(rhs);
        }
        LinearModel& operator=(LinearModel const& rhs)
        {
            if (this != &rhs)
            {
                this->base_type::operator=(rhs);
                copy(rhs);
            }
            return *this;
        }
        ~LinearModel()
        {
        }

        std::vector<constraint_type> const& constraints() const {return m_vConstraint;}
        std::vector<constraint_type>& constraints() {return m_vConstraint;}
        bool addConstraint(constraint_type const& constr, std::string name = "") 
        {
            expression_type expr = constr.expression();
            return emplaceConstraint(expr, constr.sense(), constr.rightHandSide(), name);
        }
        bool emplaceConstraint(expression_type& expr, char sense, coefficient_value_type rhs, std::string name = "")
        {
            // simplify expression 
            expr.simplify();
            // verify whether the constraint is actually a bound constraint 
            std::vector<term_type> const& vTerm = expr.terms();
            if (vTerm.empty())
                return false;
            if (vTerm.size() == 1) // one term, bound constraint  
            {
                term_type const& term = vTerm.front();
                switch (sense)
                {
                    case '>':
                        if (term.coefficient() > 0) // lower bound 
                            m_vVariableProperty.at(term.variable().id()).updateLowerBound(rhs/term.coefficient());
                        else if (term.coefficient() < 0) // upper bound 
                            m_vVariableProperty.at(term.variable().id()).updateUpperBound(rhs/term.coefficient());
                        else 
                            return false; 
                        break; 
                    case '<':
                        if (term.coefficient() > 0) // upper bound 
                            m_vVariableProperty.at(term.variable().id()).updateUpperBound(rhs/term.coefficient());
                        else if (term.coefficient() < 0) // lower bound 
                            m_vVariableProperty.at(term.variable().id()).updateLowerBound(rhs/term.coefficient());
                        else 
                            return false; 
                        break; 
                    case '=':
                        if (term.coefficient() != 0)
                        {
                            // lower bound 
                            m_vVariableProperty.at(term.variable().id()).updateLowerBound(rhs/term.coefficient());
                            // upper bound 
                            m_vVariableProperty.at(term.variable().id()).updateUpperBound(rhs/term.coefficient());
                        }
                        else 
                            return false; 
                        break; 
                    default:
                        limboAssertMsg(0, "Unknown sense %c", sense);
                }
            }
            else // actual constraint 
            {
                m_vConstraint.push_back(constraint_type()); 
                m_vConstraintName.push_back(name);
                constraint_type& constr = m_vConstraint.back(); 
                constr.setId(m_vConstraint.size()-1); 
                constr.emplaceExpression(expr); 
                constr.setSense(sense); 
                constr.setRightHandSide(rhs); 
            }
            return true; 
        }
        std::string const& constraintName(constraint_type const& constr) const {return m_vConstraintName[constr.id()];}
        std::vector<std::string> const& constraintNames() const {return m_vConstraintName;}
        std::vector<std::string>& constraintNames() {return m_vConstraintName;}
        expression_type const& objective() const {return m_objective;}
        void setObjective(expression_type const& expr) 
        {
            m_objective = expr;
            m_objective.simplify();
        }
        void emplaceObjective(expression_type& expr)
        {
            m_objective.swap(expr); 
            m_objective.simplify();
        }
        SolverProperty optimizeType() const {return m_optType;}
        void setOptimizeType(SolverProperty optType) {m_optType = optType;}
        std::vector<property_type> const& variableProperties() const {return m_vVariableProperty;}
        std::vector<property_type>& variableProperties() {return m_vVariableProperty;}
        unsigned int numVariables() const {return m_vVariableProperty.size();}
        variable_type variable(unsigned int id) const {return variable_type(id);}
        std::string variableName(variable_type const& var) const
        {
            char buf[256];
            if (m_vVariableProperty.at(var.id()).name().empty())
                limboSPrint(kNONE, buf, "x%u", var.id());
            else 
                limboSPrint(kNONE, buf, "%s", m_vVariableProperty.at(var.id()).name().c_str());
            return buf; 
        }
        void setVariableName(variable_type const& var, std::string const& name) {m_vVariableProperty[var.id()].setName(name);}
        variable_value_type variableLowerBound(variable_type const& var) {return m_vVariableProperty.at(var.id()).lowerBound();}
        void setVariableLowerBound(variable_type const& var, variable_value_type lb) {m_vVariableProperty[var.id()].setLowerBound(lb);}
        void updateVariableLowerBound(variable_type const& var, variable_value_type lb) {m_vVariableProperty[var.id()].updateLowerBound(lb);}
        variable_value_type variableUpperBound(variable_type const& var) {return m_vVariableProperty.at(var.id()).upperBound();}
        void setVariableUpperBound(variable_type const& var, variable_value_type ub) {m_vVariableProperty[var.id()].setUpperBound(ub);}
        void updateVariableUpperBound(variable_type const& var, variable_value_type ub) {m_vVariableProperty[var.id()].updateUpperBound(ub);}
        SolverProperty variableNumericType(variable_type const& var) {return m_vVariableProperty.at(var.id()).numericType();}
        void setVariableNumericType(variable_type const& var, SolverProperty type) {m_vVariableProperty[var.id()].setNumericType(type);}
        variable_type addVariable(variable_value_type lb, variable_value_type ub, SolverProperty nt, std::string name = "") 
        {
            m_vVariableProperty.push_back(property_type(lb, ub, nt, name)); 
            m_vVariableSol.push_back(lb);
            return variable_type(m_vVariableProperty.size()-1);
        }
        std::vector<variable_value_type> const& variableSolutions() const {return m_vVariableSol;}
        std::vector<variable_value_type>& variableSolutions() {return m_vVariableSol;}
        variable_value_type variableSolution(variable_type const& var) const {return m_vVariableSol.at(var.id());}
        void setVariableSolution(variable_type const& var, variable_value_type v) {m_vVariableSol.at(var.id()) = v;}
        void reserveVariables(unsigned int n)
        {
            m_vVariableProperty.reserve(n);
            m_vVariableSol.reserve(n);
        }
        void reserveConstraints(unsigned int n)
        {
            m_vConstraint.reserve(n);
            m_vConstraintName.reserve(n);
        }
        void resizeConstraints(unsigned int n)
        {
            m_vConstraint.resize(n); 
            m_vConstraintName.resize(n);
            for (unsigned int i = 0; i < n; ++i)
                m_vConstraint[i].setId(i);
        }
        void scaleObjective(coefficient_value_type factor) {m_objective *= factor;}
        void scaleConstraint(unsigned int id, coefficient_value_type factor) {m_vConstraint.at(id).scale(factor);}
        coefficient_value_type evaluateExpression(expression_type const& expr, std::vector<variable_value_type> const& vVariableSol) const 
        {
            coefficient_value_type result = 0; 
            for (typename std::vector<term_type>::const_iterator it = expr.terms().begin(), ite = expr.terms().end(); it != ite; ++it)
                result += it->coefficient()*vVariableSol.at(it->variable().id());
            result += expr.constant();
            return result;
        }
        coefficient_value_type evaluateObjective(std::vector<variable_value_type> const& vVariableSol) const 
        {
            return evaluateExpression(m_objective, vVariableSol);
        }
        coefficient_value_type evaluateObjective() const 
        {
            return evaluateExpression(m_objective, m_vVariableSol);
        }
        coefficient_value_type evaluateConstraint(constraint_type const& constr, std::vector<variable_value_type> const& vVariableSol) const 
        {
            coefficient_value_type result = constr.rightHandSide()-evaluateExpression(constr.expression(), vVariableSol);
            if (constr.sense() == '>')
                return -result; 
            else 
                return result; 
        }
        coefficient_value_type evaluateConstraint(constraint_type const& constr) const 
        {
            return evaluateConstraint(constr, m_vVariableSol);
        }
        void evaluateConstraint() const 
        {
            for (unsigned int i = 0, ie = m_vConstraint.size(); i < ie; ++i)
                limboPrint(kNONE, "C[%u] slack = %g\n", i, (double)evaluateConstraint(m_vConstraint.at(i)));
        }

        void read(std::string const& filename) 
        {
            LpParser::read(*this, filename);
        }
        void add_variable(std::string const& vname, double l = limbo::lowest<double>(), double r = std::numeric_limits<double>::max())
        {
            // in case of overflow 
            variable_value_type lb = l; 
            variable_value_type ub = r;
            if (l == -10)
                limboAssert(lb == -10);
            if (l <= (double)limbo::lowest<variable_value_type>())
                lb = limbo::lowest<variable_value_type>();
            if (r >= (double)std::numeric_limits<variable_value_type>::max())
                ub = std::numeric_limits<variable_value_type>::max();
            limboAssertMsg(lb <= ub, "failed to add bound %g <= %s <= %g", l, vname.c_str(), r); 
            if (l == -10)
                limboAssert(lb == -10);

            // no variables with the same name is allowed 
            typename std::map<std::string, variable_type>::const_iterator found = m_mName2Variable.find(vname); 
            if (found == m_mName2Variable.end()) // new variable 
                limboAssertMsg(m_mName2Variable.insert(std::make_pair(vname, addVariable(lb, ub, CONTINUOUS, vname))).second,
                        "failed to insert variable %s to hash table", vname.c_str());
            else // update variable 
            {
                m_vVariableProperty.at(found->second.id()).updateLowerBound(lb);
                m_vVariableProperty.at(found->second.id()).updateUpperBound(ub);
                limboAssertMsg(m_vVariableProperty.at(found->second.id()).lowerBound() <= m_vVariableProperty.at(found->second.id()).upperBound(), 
                        "failed to set bound %g <= %s <= %g", m_vVariableProperty.at(found->second.id()).lowerBound(), vname.c_str(), m_vVariableProperty.at(found->second.id()).upperBound());
            }
        }
        void add_constraint(std::string const& cname, LpParser::TermArray const& terms, char compare, double constant) 
        {
            expression_type expr; 
            for (LpParser::TermArray::const_iterator it = terms.begin(); it != terms.end(); ++it)
            {
                // in case variable is not added yet 
                add_variable(it->var);
                expr += m_mName2Variable.at(it->var)*it->coef; 
            }
            addConstraint(constraint_type(expr, constant, compare), cname);
        }
        void add_objective(bool minimize, LpParser::TermArray const& terms)
        {
            setOptimizeType((minimize)? MIN : MAX); 

            expression_type expr; 
            for (LpParser::TermArray::const_iterator it = terms.begin(); it != terms.end(); ++it)
            {
                // in case variable is not added yet 
                add_variable(it->var);
                expr += m_mName2Variable.at(it->var)*it->coef; 
            }
            setObjective(expr);
        }
        void set_integer(std::string const& vname, bool binary)
        {
            // in case variable is not added yet 
            add_variable(vname);
            variable_type var = m_mName2Variable.at(vname); 
            if (binary)
                m_vVariableProperty.at(var.id()).setNumericType(BINARY); 
            else 
                m_vVariableProperty.at(var.id()).setNumericType(INTEGER); 
        }

        bool print(std::string const& filename) const 
        {
            std::ofstream out (filename.c_str()); 
            if (!out.good())
                return false; 
            print(out); 
            out.close();
            return true; 
        }
        std::ostream& print(std::ostream& os = std::cout) const 
        {
            switch (optimizeType())
            {
                case MIN:
                    os << "Minimize\n";
                    break; 
                case MAX:
                    os << "Maximize\n";
                    break; 
                default:
                    os << "Unknown\n";
                    break; 
            }

            // print objective 
            print(os, objective()); 

            os << "\n\nSubject To\n";

            // print constraints 
            unsigned int i = 0; 
            for (typename std::vector<constraint_type>::const_iterator it = constraints().begin(), ite = constraints().end(); it != ite; ++it, ++i)
            {
                if (constraintName(*it).empty())
                    os << "C" << i << ": ";
                else 
                    os << constraintName(*it) << ": ";
                print(os, *it); 
                os << "\n";
            }

            // print bounds 
            os << "Bounds\n";
            i = 0; 
            for (typename std::vector<property_type>::const_iterator it = m_vVariableProperty.begin(), ite = m_vVariableProperty.end(); it != ite; ++it, ++i)
            {
                if (it->lowerBound() <= limbo::lowest<variable_value_type>() && it->upperBound() >= std::numeric_limits<variable_value_type>::max())
                    os << variableName(variable_type(i)) << " free";
                else if (it->lowerBound() <= limbo::lowest<variable_value_type>())
                    os << variableName(variable_type(i)) << " <= " << it->upperBound();
                else if (it->upperBound() >= std::numeric_limits<variable_value_type>::max())
                    os << variableName(variable_type(i)) << " >= " << it->lowerBound();
                else 
                    os << it->lowerBound() << " <= " << variableName(variable_type(i)) << " <= " << it->upperBound();
                os << "\n";
            }

            // print numeric type 
            os << "Generals\n";
            i = 0; 
            for (typename std::vector<property_type>::const_iterator it = m_vVariableProperty.begin(), ite = m_vVariableProperty.end(); it != ite; ++it, ++i)
            {
                if (it->numericType() == BINARY || it->numericType() == INTEGER)
                    os << variableName(variable_type(i)) << "\n";
            }

            os << "End\n"; 

            return os; 
        }
        std::ostream& print(std::ostream& os, term_type const& term) const 
        {
            os << term.coefficient() << " " << variableName(term.variable());
            return os; 
        }
        std::ostream& print(std::ostream& os, expression_type const& expr) const 
        {
            int i = 0; 
            for (typename std::vector<term_type>::const_iterator it = expr.terms().begin(), ite = expr.terms().end(); it != ite; ++it, ++i)
            {
                if (i)
                    os << " + "; 
                print(os, *it);
                if (i%4 == 3)
                    os << "\n";
            }
            return os; 
        }
        std::ostream& print(std::ostream& os, constraint_type const& constr) const 
        {
            print(os, constr.expression());
            if (constr.sense() == '=')
                os << " " << constr.sense() << " " << constr.rightHandSide();
            else 
                os << " " << constr.sense() << "= " << constr.rightHandSide();
            return os;  
        }
        bool printSolution(std::string const& filename) const 
        {
            std::ofstream out (filename.c_str()); 
            if (!out.good())
                return false; 
            printSolution(out); 
            out.close();
            return true; 
        }
        std::ostream& printSolution(std::ostream& os = std::cout) const 
        {
            coefficient_value_type obj = evaluateObjective(); 
            os << "# Objective " << obj << "\n";
            for (unsigned int i = 0, ie = variableSolutions().size(); i < ie; ++i)
                os << variableName(variable_type(i)) << " " << variableSolution(variable_type(i)) << "\n";
            return os; 
        }
    protected:
        void copy(LinearModel const& rhs)
        {
            m_vConstraint = rhs.m_vConstraint;
            m_objective = rhs.m_objective; 
            m_vVariableProperty = rhs.m_vVariableProperty;
            m_optType = rhs.m_optType;

            m_vVariableSol = rhs.m_vVariableSol;

            m_mName2Variable = rhs.m_mName2Variable; 
        }

        std::vector<constraint_type> m_vConstraint; 
        expression_type m_objective; 
        std::vector<property_type> m_vVariableProperty; 
        SolverProperty m_optType; 
        std::vector<std::string> m_vConstraintName; 

        std::vector<variable_value_type> m_vVariableSol; 

        std::map<std::string, variable_type> m_mName2Variable; 
};

template <typename T, typename I, int StartingIndex = 1>
struct MatrixCSR 
{
    typedef T value_type; 
    typedef I index_type; 

    value_type* vElement; 
    index_type* vColumn; 
    index_type* vRowBeginIndex; 
    index_type numRows; 
    index_type numColumns; 
    index_type numElements; 
    static index_type s_startingIndex; 

    MatrixCSR()
        : vElement(NULL)
        , vColumn(NULL)
        , vRowBeginIndex(NULL)
        , numRows(0)
        , numColumns(0)
        , numElements(0)
    {
    }
    MatrixCSR(MatrixCSR const& rhs)
    {
        copy(rhs);
    }
    MatrixCSR& operator=(MatrixCSR const& rhs)
    {
        if (this != &rhs)
            copy(rhs);
        return *this;
    }
    ~MatrixCSR()
    {
        reset();
    }

    void initialize(index_type nr, index_type nc, index_type nv)
    {
        reset();
        numRows = nr; 
        numColumns = nc; 
        numElements = nv; 
        vElement = new value_type [numElements];
        vColumn = new index_type [numElements];
        vRowBeginIndex = new index_type [numRows+1];
    }

    void reset() 
    {
        if (vElement)
            delete [] vElement; 
        if (vColumn)
            delete [] vColumn; 
        if (vRowBeginIndex)
            delete [] vRowBeginIndex;
        vElement = NULL;
        vColumn = NULL;
        vRowBeginIndex = NULL;
        numRows = 0; 
        numColumns = 0; 
        numElements = 0; 
    }

    void copy(MatrixCSR const& rhs)
    {
        reset(); 
        numRows = rhs.numRows; 
        numColumns = rhs.numColumns; 
        numElements = rhs.numElements;
        if (rhs.vElement)
        {
            vElement = new value_type [numElements];
            vColumn = new index_type [numElements]; 
            vRowBeginIndex = new index_type [numRows+1]; 
            std::copy(rhs.vElement, rhs.vElement+numElements, vElement); 
            std::copy(rhs.vColumn, rhs.vColumn+numElements, vColumn); 
            std::copy(rhs.vRowBeginIndex, rhs.vRowBeginIndex+numRows+1, vRowBeginIndex);
        }
    }

    value_type at(index_type i, index_type j) const
    {
#ifdef DEBUG_SOLVERS 
        limboAssert(i < numRows); 
        limboAssert(i+1 <= numRows);
#endif
        index_type* elementColumnBegin = vColumn+vRowBeginIndex[i]-vRowBeginIndex[0];
        index_type* elementColumnEnd = vColumn+vRowBeginIndex[i+1]-vRowBeginIndex[0];
        // binary search to see if any element match column j 
        index_type* found = std::lower_bound(elementColumnBegin, elementColumnEnd, j+s_startingIndex); 
        if (found != elementColumnEnd && *found == j+s_startingIndex)
        {
#ifdef DEBUG_SOLVERS
            limboAssert(std::distance(vColumn, found)-s_startingIndex < numElements);
#endif
            return vElement[std::distance(vColumn, found)]; 
        }
        else 
            return 0; 
    }

    void set(index_type nr, index_type nc, LinearConstraint<value_type> const* vConstraint) 
    {
        numRows = nr; 
        numColumns = nc; 
        numElements = 0; 
        vRowBeginIndex = new index_type [numRows+1]; 
        vRowBeginIndex[0] = s_startingIndex;

        typedef LinearConstraint<value_type> constraint_type;
        constraint_type const* it = vConstraint; 
        constraint_type const* ite = vConstraint+nr; 
        // initialize vRowBeginIndex 
        index_type i = 1; 
        for (; it != ite; ++it, ++i)
        {
#ifdef DEBUG_SOLVERS
            limboAssert(i <= numRows); 
#endif
            vRowBeginIndex[i] = vRowBeginIndex[i-1]+it->expression().terms().size(); 
        }
        // last element of vRowBeginIndex denotes the total number of elements 
        numElements = vRowBeginIndex[numRows]; 

        // initialize vElement and vColumn 
        vElement = new value_type [numElements]; 
        vColumn = new index_type [numElements];
        i = 0; 
        for (it = vConstraint; it != ite; ++it)
        {
            for (typename std::vector<typename constraint_type::term_type>::const_iterator itt = it->expression().terms().begin(), itte = it->expression().terms().end(); itt != itte; ++itt, ++i)
            {
#ifdef DEBUG_SOLVERS
                limboAssert(i < numElements);
#endif
                vElement[i] = itt->coefficient(); 
                vColumn[i] = itt->variable().id()+s_startingIndex;
            }
        }
    }
};

template <typename T, typename I, int StartingIndex>
typename MatrixCSR<T, I, StartingIndex>::index_type MatrixCSR<T, I, StartingIndex>::s_startingIndex = StartingIndex; 

} // namespace solvers 
} // namespace limbo 

#endif