/* -*- mode: C++; indent-tabs-mode: nil; -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library.

 *

 * Copyright (C) 2003-2013

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#ifndef LEMON_LP_BASE_H

#define LEMON_LP_BASE_H

#include<iostream>

#include<vector>

#include<map>

#include<limits>

#include<lemon/math.h>

#include<lemon/error.h>

#include<lemon/assert.h>

#include<lemon/core.h>

#include<lemon/bits/solver_bits.h>

#include<lemon/bits/stl_iterators.h>

///\file

///\brief The interface of the LP solver interface.

///\ingroup lp_group

namespace lemon {

  ///Common base class for LP and MIP solvers

  ///Usually this class is not used directly, please use one of the concrete

  ///implementations of the solver interface.

  ///\ingroup lp_group

  class LpBase {

  protected:

    _solver_bits::VarIndex _rows;

    _solver_bits::VarIndex _cols;

  public:

    ///Possible outcomes of an LP solving procedure

    enum SolveExitStatus {

      /// = 0. It means that the problem has been successfully solved: either

      ///an optimal solution has been found or infeasibility/unboundedness

      ///has been proved.

      SOLVED = 0,

      /// = 1. Any other case (including the case when some user specified

      ///limit has been exceeded).

      UNSOLVED = 1

    };

    ///Direction of the optimization

    enum Sense {

      /// Minimization

      MIN,

      /// Maximization

      MAX

    };

    ///Enum for \c messageLevel() parameter

    enum MessageLevel {

      /// No output (default value).

      MESSAGE_NOTHING,

      /// Error messages only.

      MESSAGE_ERROR,

      /// Warnings.

      MESSAGE_WARNING,

      /// Normal output.

      MESSAGE_NORMAL,

      /// Verbose output.

      MESSAGE_VERBOSE

    };

    ///The floating point type used by the solver

    typedef double Value;

    ///The infinity constant

    static const Value INF;

    ///The not a number constant

    static const Value NaN;

    friend class Col;

    friend class ColIt;

    friend class Row;

    friend class RowIt;

    ///Refer to a column of the LP.

    ///This type is used to refer to a column of the LP.

    ///

    ///Its value remains valid and correct even after the addition or erase of

    ///other columns.

    ///

    ///\note This class is similar to other Item types in LEMON, like

    ///Node and Arc types in digraph.

    class Col {

      friend class LpBase;

    protected:

      int _id;

      explicit Col(int id) : _id(id) {}

    public:

      typedef Value ExprValue;

      typedef True LpCol;

      /// Default constructor

      /// \warning The default constructor sets the Col to an

      /// undefined value.

      Col() {}

      /// Invalid constructor \& conversion.

      /// This constructor initializes the Col to be invalid.

      /// \sa Invalid for more details.

      Col(const Invalid&) : _id(-1) {}

      /// Equality operator

      /// Two \ref Col "Col"s are equal if and only if they point to

      /// the same LP column or both are invalid.

      bool operator==(Col c) const  {return _id == c._id;}

      /// Inequality operator

      /// \sa operator==(Col c)

      ///

      bool operator!=(Col c) const  {return _id != c._id;}

      /// Artificial ordering operator.

      /// To allow the use of this object in std::map or similar

      /// associative container we require this.

      ///

      /// \note This operator only have to define some strict ordering of

      /// the items; this order has nothing to do with the iteration

      /// ordering of the items.

      bool operator<(Col c) const  {return _id < c._id;}

    };

    ///Iterator for iterate over the columns of an LP problem

    /// Its usage is quite simple, for example, you can count the number

    /// of columns in an LP \c lp:

    ///\code

    /// int count=0;

    /// for (LpBase::ColIt c(lp); c!=INVALID; ++c) ++count;

    ///\endcode

    class ColIt : public Col {

      const LpBase *_solver;

    public:

      /// Default constructor

      /// \warning The default constructor sets the iterator

      /// to an undefined value.

      ColIt() {}

      /// Sets the iterator to the first Col

      /// Sets the iterator to the first Col.

      ///

      ColIt(const LpBase &solver) : _solver(&solver)

      {

        _solver->_cols.firstItem(_id);

      }

      /// Invalid constructor \& conversion

      /// Initialize the iterator to be invalid.

      /// \sa Invalid for more details.

      ColIt(const Invalid&) : Col(INVALID) {}

      /// Next column

      /// Assign the iterator to the next column.

      ///

      ColIt &operator++()

      {

        _solver->_cols.nextItem(_id);

        return *this;

      }

    };

    /// \brief Gets the collection of the columns of the LP problem.

    ///

    /// This function can be used for iterating on

    /// the columns of the LP problem. It returns a wrapped ColIt, which looks

    /// like an STL container (by having begin() and end())

    /// which you can use in range-based for loops, STL algorithms, etc.

    /// For example you can write:

    ///\code

    /// for(auto c: lp.cols())

    ///   doSomething(c);

    ///\endcode

    LemonRangeWrapper1<ColIt, LpBase> cols() {

      return LemonRangeWrapper1<ColIt, LpBase>(*this);

    }

    /// \brief Returns the ID of the column.

    static int id(const Col& col) { return col._id; }

    /// \brief Returns the column with the given ID.

    ///

    /// \pre The argument should be a valid column ID in the LP problem.

    static Col colFromId(int id) { return Col(id); }

    ///Refer to a row of the LP.

    ///This type is used to refer to a row of the LP.

    ///

    ///Its value remains valid and correct even after the addition or erase of

    ///other rows.

    ///

    ///\note This class is similar to other Item types in LEMON, like

    ///Node and Arc types in digraph.

    class Row {

      friend class LpBase;

    protected:

      int _id;

      explicit Row(int id) : _id(id) {}

    public:

      typedef Value ExprValue;

      typedef True LpRow;

      /// Default constructor

      /// \warning The default constructor sets the Row to an

      /// undefined value.

      Row() {}

      /// Invalid constructor \& conversion.

      /// This constructor initializes the Row to be invalid.

      /// \sa Invalid for more details.

      Row(const Invalid&) : _id(-1) {}

      /// Equality operator

      /// Two \ref Row "Row"s are equal if and only if they point to

      /// the same LP row or both are invalid.

      bool operator==(Row r) const  {return _id == r._id;}

      /// Inequality operator

      /// \sa operator==(Row r)

      ///

      bool operator!=(Row r) const  {return _id != r._id;}

      /// Artificial ordering operator.

      /// To allow the use of this object in std::map or similar

      /// associative container we require this.

      ///

      /// \note This operator only have to define some strict ordering of

      /// the items; this order has nothing to do with the iteration

      /// ordering of the items.

      bool operator<(Row r) const  {return _id < r._id;}

    };

    ///Iterator for iterate over the rows of an LP problem

    /// Its usage is quite simple, for example, you can count the number

    /// of rows in an LP \c lp:

    ///\code

    /// int count=0;

    /// for (LpBase::RowIt c(lp); c!=INVALID; ++c) ++count;

    ///\endcode

    class RowIt : public Row {

      const LpBase *_solver;

    public:

      /// Default constructor

      /// \warning The default constructor sets the iterator

      /// to an undefined value.

      RowIt() {}

      /// Sets the iterator to the first Row

      /// Sets the iterator to the first Row.

      ///

      RowIt(const LpBase &solver) : _solver(&solver)

      {

        _solver->_rows.firstItem(_id);

      }

      /// Invalid constructor \& conversion

      /// Initialize the iterator to be invalid.

      /// \sa Invalid for more details.

      RowIt(const Invalid&) : Row(INVALID) {}

      /// Next row

      /// Assign the iterator to the next row.

      ///

      RowIt &operator++()

      {

        _solver->_rows.nextItem(_id);

        return *this;

      }

    };

    /// \brief Gets the collection of the rows of the LP problem.

    ///

    /// This function can be used for iterating on

    /// the rows of the LP problem. It returns a wrapped RowIt, which looks

    /// like an STL container (by having begin() and end())

    /// which you can use in range-based for loops, STL algorithms, etc.

    /// For example you can write:

    ///\code

    /// for(auto c: lp.rows())

    ///   doSomething(c);

    ///\endcode

    LemonRangeWrapper1<RowIt, LpBase> rows() {

      return LemonRangeWrapper1<RowIt, LpBase>(*this);

    }

    /// \brief Returns the ID of the row.

    static int id(const Row& row) { return row._id; }

    /// \brief Returns the row with the given ID.

    ///

    /// \pre The argument should be a valid row ID in the LP problem.

    static Row rowFromId(int id) { return Row(id); }

  public:

    ///Linear expression of variables and a constant component

    ///This data structure stores a linear expression of the variables

    ///(\ref Col "Col"s) and also has a constant component.

    ///

    ///There are several ways to access and modify the contents of this

    ///container.

    ///\code

    ///e[v]=5;

    ///e[v]+=12;

    ///e.erase(v);

    ///\endcode

    ///or you can also iterate through its elements.

    ///\code

    ///double s=0;

    ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)

    ///  s+=*i * primal(i);

    ///\endcode

    ///(This code computes the primal value of the expression).

    ///- Numbers (<tt>double</tt>'s)

    ///and variables (\ref Col "Col"s) directly convert to an

    ///\ref Expr and the usual linear operations are defined, so

    ///\code

    ///v+w

    ///2*v-3.12*(v-w/2)+2

    ///v*2.1+(3*v+(v*12+w+6)*3)/2

    ///\endcode

    ///are valid expressions.

    ///The usual assignment operations are also defined.

    ///\code

    ///e=v+w;

    ///e+=2*v-3.12*(v-w/2)+2;

    ///e*=3.4;

    ///e/=5;

    ///\endcode

    ///- The constant member can be set and read by dereference

    ///  operator (unary *)

    ///

    ///\code

    ///*e=12;

    ///double c=*e;

    ///\endcode

    ///

    ///\sa Constr

    class Expr {

      friend class LpBase;

    public:

      /// The key type of the expression

      typedef LpBase::Col Key;

      /// The value type of the expression

      typedef LpBase::Value Value;

    protected:

      Value const_comp;

      std::map<int, Value> comps;

    public:

      typedef True SolverExpr;

      /// Default constructor

      /// Construct an empty expression, the coefficients and

      /// the constant component are initialized to zero.

      Expr() : const_comp(0) {}

      /// Construct an expression from a column

      /// Construct an expression, which has a term with \c c variable

      /// and 1.0 coefficient.

      Expr(const Col &c) : const_comp(0) {

        typedef std::map<int, Value>::value_type pair_type;

        comps.insert(pair_type(id(c), 1));

      }

      /// Construct an expression from a constant

      /// Construct an expression, which's constant component is \c v.

      ///

      Expr(const Value &v) : const_comp(v) {}

      /// Returns the coefficient of the column

      Value operator[](const Col& c) const {

        std::map<int, Value>::const_iterator it=comps.find(id(c));

        if (it != comps.end()) {

          return it->second;

        } else {

          return 0;

        }

      }

      /// Returns the coefficient of the column

      Value& operator[](const Col& c) {

        return comps[id(c)];

      }

      /// Sets the coefficient of the column

      void set(const Col &c, const Value &v) {

        if (v != 0.0) {

          typedef std::map<int, Value>::value_type pair_type;

          comps.insert(pair_type(id(c), v));

        } else {

          comps.erase(id(c));

        }

      }

      /// Returns the constant component of the expression

      Value& operator*() { return const_comp; }

      /// Returns the constant component of the expression

      const Value& operator*() const { return const_comp; }

      /// \brief Removes the coefficients which's absolute value does

      /// not exceed \c epsilon. It also sets to zero the constant

      /// component, if it does not exceed epsilon in absolute value.

      void simplify(Value epsilon = 0.0) {

        std::map<int, Value>::iterator it=comps.begin();

        while (it != comps.end()) {

          std::map<int, Value>::iterator jt=it;

          ++jt;

          if (std::fabs((*it).second) <= epsilon) comps.erase(it);

          it=jt;

        }

        if (std::fabs(const_comp) <= epsilon) const_comp = 0;

      }

      void simplify(Value epsilon = 0.0) const {

        const_cast<Expr*>(this)->simplify(epsilon);

      }

      ///Sets all coefficients and the constant component to 0.

      void clear() {

        comps.clear();

        const_comp=0;

      }

      ///Compound assignment

      Expr &operator+=(const Expr &e) {

        for (std::map<int, Value>::const_iterator it=e.comps.begin();

             it!=e.comps.end(); ++it)

          comps[it->first]+=it->second;

        const_comp+=e.const_comp;

        return *this;

      }

      ///Compound assignment

      Expr &operator-=(const Expr &e) {

        for (std::map<int, Value>::const_iterator it=e.comps.begin();

             it!=e.comps.end(); ++it)

          comps[it->first]-=it->second;

        const_comp-=e.const_comp;

        return *this;

      }

      ///Multiply with a constant

      Expr &operator*=(const Value &v) {

        for (std::map<int, Value>::iterator it=comps.begin();

             it!=comps.end(); ++it)

          it->second*=v;

        const_comp*=v;

        return *this;

      }

      ///Division with a constant

      Expr &operator/=(const Value &c) {

        for (std::map<int, Value>::iterator it=comps.begin();

             it!=comps.end(); ++it)

          it->second/=c;

        const_comp/=c;

        return *this;

      }

      ///Iterator over the expression

      ///The iterator iterates over the terms of the expression.

      ///

      ///\code

      ///double s=0;

      ///for(LpBase::Expr::CoeffIt i(e);i!=INVALID;++i)

      ///  s+= *i * primal(i);

      ///\endcode

      class CoeffIt {

      private:

        std::map<int, Value>::iterator _it, _end;

      public:

        /// Sets the iterator to the first term

        /// Sets the iterator to the first term of the expression.

        ///

        CoeffIt(Expr& e)

          : _it(e.comps.begin()), _end(e.comps.end()){}

        /// Convert the iterator to the column of the term

        operator Col() const {

          return colFromId(_it->first);

        }

        /// Returns the coefficient of the term

        Value& operator*() { return _it->second; }

        /// Returns the coefficient of the term

        const Value& operator*() const { return _it->second; }

        /// Next term

        /// Assign the iterator to the next term.

        ///

        CoeffIt& operator++() { ++_it; return *this; }

        /// Equality operator

        bool operator==(Invalid) const { return _it == _end; }

        /// Inequality operator

        bool operator!=(Invalid) const { return _it != _end; }

      };

      /// Const iterator over the expression

      ///The iterator iterates over the terms of the expression.

      ///

      ///\code

      ///double s=0;

      ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)

      ///  s+=*i * primal(i);

      ///\endcode

      class ConstCoeffIt {

      private:

        std::map<int, Value>::const_iterator _it, _end;

      public:

        /// Sets the iterator to the first term

        /// Sets the iterator to the first term of the expression.

        ///

        ConstCoeffIt(const Expr& e)

          : _it(e.comps.begin()), _end(e.comps.end()){}

        /// Convert the iterator to the column of the term

        operator Col() const {

          return colFromId(_it->first);

        }

        /// Returns the coefficient of the term

        const Value& operator*() const { return _it->second; }

        /// Next term

        /// Assign the iterator to the next term.

        ///

        ConstCoeffIt& operator++() { ++_it; return *this; }

        /// Equality operator

        bool operator==(Invalid) const { return _it == _end; }

        /// Inequality operator

        bool operator!=(Invalid) const { return _it != _end; }

      };

    };

    ///Linear constraint

    ///This data stucture represents a linear constraint in the LP.

    ///Basically it is a linear expression with a lower or an upper bound

    ///(or both). These parts of the constraint can be obtained by the member

    ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),

    ///respectively.

    ///There are two ways to construct a constraint.

    ///- You can set the linear expression and the bounds directly

    ///  by the functions above.

    ///- The operators <tt>\<=</tt>, <tt>==</tt> and  <tt>\>=</tt>

    ///  are defined between expressions, or even between constraints whenever

    ///  it makes sense. Therefore if \c e and \c f are linear expressions and

    ///  \c s and \c t are numbers, then the followings are valid expressions

    ///  and thus they can be used directly e.g. in \ref addRow() whenever

    ///  it makes sense.

    ///\code

    ///  e<=s

    ///  e<=f

    ///  e==f

    ///  s<=e<=t

    ///  e>=t

    ///\endcode

    ///\warning The validity of a constraint is checked only at run

    ///time, so e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will

    ///compile, but will fail an assertion.

    class Constr

    {

    public:

      typedef LpBase::Expr Expr;

      typedef Expr::Key Key;

      typedef Expr::Value Value;

    protected:

      Expr _expr;

      Value _lb,_ub;

    public:

      ///\e

      Constr() : _expr(), _lb(NaN), _ub(NaN) {}

      ///\e

      Constr(Value lb, const Expr &e, Value ub) :

        _expr(e), _lb(lb), _ub(ub) {}

      Constr(const Expr &e) :

        _expr(e), _lb(NaN), _ub(NaN) {}

      ///\e

      void clear()

      {

        _expr.clear();

        _lb=_ub=NaN;

      }

      ///Reference to the linear expression

      Expr &expr() { return _expr; }

      ///Cont reference to the linear expression

      const Expr &expr() const { return _expr; }

      ///Reference to the lower bound.

      ///\return

      ///- \ref INF "INF": the constraint is lower unbounded.

      ///- \ref NaN "NaN": lower bound has not been set.

      ///- finite number: the lower bound

      Value &lowerBound() { return _lb; }

      ///The const version of \ref lowerBound()

      const Value &lowerBound() const { return _lb; }

      ///Reference to the upper bound.

      ///\return

      ///- \ref INF "INF": the constraint is upper unbounded.

      ///- \ref NaN "NaN": upper bound has not been set.

      ///- finite number: the upper bound

      Value &upperBound() { return _ub; }

      ///The const version of \ref upperBound()

      const Value &upperBound() const { return _ub; }

      ///Is the constraint lower bounded?

      bool lowerBounded() const {

        return _lb != -INF && !isNaN(_lb);

      }

      ///Is the constraint upper bounded?

      bool upperBounded() const {

        return _ub != INF && !isNaN(_ub);

      }

    };

    ///Linear expression of rows

    ///This data structure represents a column of the matrix,

    ///thas is it strores a linear expression of the dual variables

    ///(\ref LpBase::Row "Row"s).

    ///

    ///There are several ways to access and modify the contents of this

    ///container.

    ///\code

    ///e[v]=5;

    ///e[v]+=12;

    ///e.erase(v);

    ///\endcode

    ///or you can also iterate through its elements.

    ///\code

    ///double s=0;

    ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)

    ///  s+=*i;

    ///\endcode

    ///(This code computes the sum of all coefficients).

    ///- Numbers (<tt>double</tt>'s)

    ///and variables (\ref LpBase::Row "Row"s) directly convert to an

    ///\ref DualExpr and the usual linear operations are defined, so

    ///\code

    ///v+w

    ///2*v-3.12*(v-w/2)

    ///v*2.1+(3*v+(v*12+w)*3)/2

    ///\endcode

    ///are valid \ref DualExpr dual expressions.

    ///The usual assignment operations are also defined.

    ///\code

    ///e=v+w;

    ///e+=2*v-3.12*(v-w/2);

    ///e*=3.4;

    ///e/=5;

    ///\endcode

    ///

    ///\sa Expr

    class DualExpr {

      friend class LpBase;

    public:

      /// The key type of the expression

      typedef LpBase::Row Key;

      /// The value type of the expression

      typedef LpBase::Value Value;

    protected:

      std::map<int, Value> comps;

    public:

      typedef True SolverExpr;

      /// Default constructor

      /// Construct an empty expression, the coefficients are

      /// initialized to zero.

      DualExpr() {}

      /// Construct an expression from a row

      /// Construct an expression, which has a term with \c r dual

      /// variable and 1.0 coefficient.

      DualExpr(const Row &r) {

        typedef std::map<int, Value>::value_type pair_type;

        comps.insert(pair_type(id(r), 1));

      }

      /// Returns the coefficient of the row

      Value operator[](const Row& r) const {

        std::map<int, Value>::const_iterator it = comps.find(id(r));

        if (it != comps.end()) {

          return it->second;

        } else {

          return 0;

        }

      }

      /// Returns the coefficient of the row

      Value& operator[](const Row& r) {

        return comps[id(r)];

      }

      /// Sets the coefficient of the row

      void set(const Row &r, const Value &v) {

        if (v != 0.0) {

          typedef std::map<int, Value>::value_type pair_type;

          comps.insert(pair_type(id(r), v));

        } else {

          comps.erase(id(r));

        }

      }

      /// \brief Removes the coefficients which's absolute value does

      /// not exceed \c epsilon.

      void simplify(Value epsilon = 0.0) {

        std::map<int, Value>::iterator it=comps.begin();

        while (it != comps.end()) {

          std::map<int, Value>::iterator jt=it;

          ++jt;

          if (std::fabs((*it).second) <= epsilon) comps.erase(it);

          it=jt;

        }

      }

      void simplify(Value epsilon = 0.0) const {

        const_cast<DualExpr*>(this)->simplify(epsilon);

      }

      ///Sets all coefficients to 0.

      void clear() {

        comps.clear();

      }

      ///Compound assignment

      DualExpr &operator+=(const DualExpr &e) {

        for (std::map<int, Value>::const_iterator it=e.comps.begin();

             it!=e.comps.end(); ++it)

          comps[it->first]+=it->second;

        return *this;

      }

      ///Compound assignment

      DualExpr &operator-=(const DualExpr &e) {

        for (std::map<int, Value>::const_iterator it=e.comps.begin();

             it!=e.comps.end(); ++it)

          comps[it->first]-=it->second;

        return *this;

      }

      ///Multiply with a constant

      DualExpr &operator*=(const Value &v) {

        for (std::map<int, Value>::iterator it=comps.begin();

             it!=comps.end(); ++it)

          it->second*=v;

        return *this;

      }

      ///Division with a constant

      DualExpr &operator/=(const Value &v) {

        for (std::map<int, Value>::iterator it=comps.begin();

             it!=comps.end(); ++it)

          it->second/=v;

        return *this;

      }

      ///Iterator over the expression

      ///The iterator iterates over the terms of the expression.

      ///

      ///\code

      ///double s=0;

      ///for(LpBase::DualExpr::CoeffIt i(e);i!=INVALID;++i)

      ///  s+= *i * dual(i);

      ///\endcode

      class CoeffIt {

      private:

        std::map<int, Value>::iterator _it, _end;

      public:

        /// Sets the iterator to the first term

        /// Sets the iterator to the first term of the expression.

        ///

        CoeffIt(DualExpr& e)

          : _it(e.comps.begin()), _end(e.comps.end()){}

        /// Convert the iterator to the row of the term

        operator Row() const {

          return rowFromId(_it->first);

        }

        /// Returns the coefficient of the term

        Value& operator*() { return _it->second; }

        /// Returns the coefficient of the term

        const Value& operator*() const { return _it->second; }

        /// Next term

        /// Assign the iterator to the next term.

        ///

        CoeffIt& operator++() { ++_it; return *this; }

        /// Equality operator

        bool operator==(Invalid) const { return _it == _end; }

        /// Inequality operator

        bool operator!=(Invalid) const { return _it != _end; }

      };

      ///Iterator over the expression

      ///The iterator iterates over the terms of the expression.

      ///

      ///\code

      ///double s=0;

      ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)

      ///  s+= *i * dual(i);

      ///\endcode

      class ConstCoeffIt {

      private:

        std::map<int, Value>::const_iterator _it, _end;

      public:

        /// Sets the iterator to the first term

        /// Sets the iterator to the first term of the expression.

        ///

        ConstCoeffIt(const DualExpr& e)

          : _it(e.comps.begin()), _end(e.comps.end()){}

        /// Convert the iterator to the row of the term

        operator Row() const {

          return rowFromId(_it->first);

        }

        /// Returns the coefficient of the term

        const Value& operator*() const { return _it->second; }

        /// Next term

        /// Assign the iterator to the next term.

        ///

        ConstCoeffIt& operator++() { ++_it; return *this; }

        /// Equality operator

        bool operator==(Invalid) const { return _it == _end; }

        /// Inequality operator

        bool operator!=(Invalid) const { return _it != _end; }

      };

    };

  protected:

    class InsertIterator {

    private:

      std::map<int, Value>& _host;

      const _solver_bits::VarIndex& _index;

    public:

      typedef std::output_iterator_tag iterator_category;

      typedef void difference_type;

      typedef void value_type;

      typedef void reference;

      typedef void pointer;

      InsertIterator(std::map<int, Value>& host,

                   const _solver_bits::VarIndex& index)

        : _host(host), _index(index) {}

      InsertIterator& operator=(const std::pair<int, Value>& value) {

        typedef std::map<int, Value>::value_type pair_type;

        _host.insert(pair_type(_index[value.first], value.second));

        return *this;

      }

      InsertIterator& operator*() { return *this; }

      InsertIterator& operator++() { return *this; }

      InsertIterator operator++(int) { return *this; }

    };

    class ExprIterator {

    private:

      std::map<int, Value>::const_iterator _host_it;

      const _solver_bits::VarIndex& _index;

    public:

      typedef std::bidirectional_iterator_tag iterator_category;

      typedef std::ptrdiff_t difference_type;

      typedef const std::pair<int, Value> value_type;

      typedef value_type reference;

      class pointer {

      public:

        pointer(value_type& _value) : value(_value) {}

        value_type* operator->() { return &value; }

      private:

        value_type value;

      };

      ExprIterator(const std::map<int, Value>::const_iterator& host_it,

                   const _solver_bits::VarIndex& index)

        : _host_it(host_it), _index(index) {}

      reference operator*() {

        return std::make_pair(_index(_host_it->first), _host_it->second);

      }

      pointer operator->() {

        return pointer(operator*());

      }

      ExprIterator& operator++() { ++_host_it; return *this; }

      ExprIterator operator++(int) {

        ExprIterator tmp(*this); ++_host_it; return tmp;

      }

      ExprIterator& operator--() { --_host_it; return *this; }

      ExprIterator operator--(int) {

        ExprIterator tmp(*this); --_host_it; return tmp;

      }

      bool operator==(const ExprIterator& it) const {

        return _host_it == it._host_it;

      }

      bool operator!=(const ExprIterator& it) const {

        return _host_it != it._host_it;

      }

    };

  protected:

    //Abstract virtual functions

    virtual int _addColId(int col) { return _cols.addIndex(col); }

    virtual int _addRowId(int row) { return _rows.addIndex(row); }

    virtual void _eraseColId(int col) { _cols.eraseIndex(col); }

    virtual void _eraseRowId(int row) { _rows.eraseIndex(row); }

    virtual int _addCol() = 0;

    virtual int _addRow() = 0;

    virtual int _addRow(Value l, ExprIterator b, ExprIterator e, Value u) {

      int row = _addRow();

      _setRowCoeffs(row, b, e);

      _setRowLowerBound(row, l);

      _setRowUpperBound(row, u);

      return row;

    }

    virtual void _eraseCol(int col) = 0;

    virtual void _eraseRow(int row) = 0;

    virtual void _getColName(int col, std::string& name) const = 0;

    virtual void _setColName(int col, const std::string& name) = 0;

    virtual int _colByName(const std::string& name) const = 0;