MultiKnapsackLagRelax.h

Go to the documentation of this file.

#ifndef LIMBO_SOLVERS_MULTIKNAPSACKLAGRELAX_H
#define LIMBO_SOLVERS_MULTIKNAPSACKLAGRELAX_H

#include <cmath>
#include <limbo/solvers/Solvers.h>
#include <limbo/solvers/Numerical.h>

namespace limbo 
{
namespace solvers 
{

// forward declaration 
template <typename T>
class LagMultiplierUpdater;
template <typename T>
class SubGradientDescent; 
template <typename T, typename V>
class ProblemScaler; 
template <typename T, typename V>
class L2NormScaler; 
template <typename T, typename V>
class FeasibleSearcher;
template <typename T, typename V>
class SearchByAdjustCoefficient;

template <typename T, typename V>
class MultiKnapsackLagRelax
{
    public:
        typedef LinearModel<T, V> model_type; 
        typedef typename model_type::coefficient_value_type coefficient_value_type; 
        typedef typename model_type::variable_value_type variable_value_type; 
        typedef typename model_type::variable_type variable_type;
        typedef typename model_type::constraint_type constraint_type; 
        typedef typename model_type::expression_type expression_type; 
        typedef typename model_type::term_type term_type; 
        typedef typename model_type::property_type property_type;
        typedef MatrixCSR<coefficient_value_type, int, 1> matrix_type; 
        typedef LagMultiplierUpdater<coefficient_value_type> updater_type; 
        typedef ProblemScaler<coefficient_value_type, variable_value_type> scaler_type; 
        typedef FeasibleSearcher<coefficient_value_type, variable_value_type> searcher_type; 
        
        MultiKnapsackLagRelax(model_type* model)
        {
            // T must be a floating point number 
            limboStaticAssert(!std::numeric_limits<T>::is_integer);
            // V must be a floating point number 
            //limboStaticAssert(!std::numeric_limits<V>::is_integer);

            m_model = model; 

            m_vObjCoef = NULL; 
            m_constrMatrix.reset();
            m_vConstrRhs = NULL; 

            m_vGroupedVariable = NULL;
            m_vVariableGroupBeginIndex = NULL; 
            m_numGroups = 0; 

            m_vLagMultiplier = NULL;
            m_vNewLagMultiplier = NULL;
            m_vSlackness = NULL;

            m_objConstant = 0; 
            m_lagObj = 0; 
            m_lagObjFlag = false; 

            m_iter = 0; 
            m_maxIters = 1000; 
            m_bestObj = std::numeric_limits<coefficient_value_type>::max();
            m_useInitialSol = false; 
        }
        MultiKnapsackLagRelax(MultiKnapsackLagRelax const& rhs)
        {
            copy(rhs);
        }
        MultiKnapsackLagRelax& operator=(MultiKnapsackLagRelax const& rhs)
        {
            if (this != &rhs)
                copy(rhs);
            return *this;
        }
        ~MultiKnapsackLagRelax()
        {
            // recycle model 
            destroy(); 
        }
        
        SolverProperty operator()(updater_type* updater = NULL, scaler_type* scaler = NULL, searcher_type* searcher = NULL)
        {
            return solve(updater, scaler, searcher);
        }

        unsigned int maxIterations() const;
        void setMaxIterations(unsigned int maxIter);
        bool useInitialSolutions() const;
        void setUseInitialSolutions(bool f); 
        bool lagObjFlag() const; 
        void setLagObjFlag(bool f); 

        unsigned int numNegativeSlackConstraints(bool evaluateFlag); 

        bool printVariableGroup(std::string const& filename) const; 
        std::ostream& printVariableGroup(std::ostream& os = std::cout) const; 
        bool printObjCoef(std::string const& filename) const; 
        std::ostream& printObjCoef(std::ostream& os = std::cout) const; 
        bool printLagMultiplier(std::string const& filename) const; 
        std::ostream& printLagMultiplier(std::ostream& os = std::cout) const; 
        std::ostream& printConstraint(std::ostream& os, std::string const& name) const; 
        template <typename TT>
        void printArray(unsigned int n, TT const* array, bool nonzeroFlag) const; 

        struct CapacityConstraintPred
        {
            std::vector<constraint_type> const& vConstraint; 

            CapacityConstraintPred(std::vector<constraint_type> const& v) : vConstraint(v) {}

            bool operator()(constraint_type const& constr) const 
            {
                return constr.sense() != '='; 
            }
            bool operator()(unsigned int id) const 
            {
                return this->operator()(vConstraint[id]); 
            }
        };

        struct CompareVariableByCoefficient
        {
            coefficient_value_type const* vObjCoef; 

            CompareVariableByCoefficient(coefficient_value_type const* v) 
                : vObjCoef(v)
            {
            }

            bool operator()(variable_type const& v1, variable_type const& v2) const 
            {
                return vObjCoef[v1.id()] < vObjCoef[v2.id()]; 
            }
        };

    protected:
        void copy(MultiKnapsackLagRelax const& rhs);
        void destroy();
        SolverProperty solve(updater_type* updater, scaler_type* scaler, searcher_type* searcher);
        SolverProperty solveSubproblems(updater_type* updater, unsigned int beginIter, unsigned int endIter);
        void scale(scaler_type* scaler); 
        void unscale(); 
        void prepare();
        void updateLagMultipliers(updater_type* updater); 
        void solveLag(); 
        void computeSlackness(); 
        coefficient_value_type evaluateLagObjective() const;
        SolverProperty converge(); 
        SolverProperty postProcess(updater_type* updater, searcher_type* searcher, SolverProperty status); 
        template <typename TT, typename VV>
        void bMinusAx(matrix_type const& A, VV const* x, TT const* b, TT* y) const; 

        model_type* m_model; 

        coefficient_value_type* m_vObjCoef; 
        matrix_type m_constrMatrix; 
        coefficient_value_type* m_vConstrRhs; 

        variable_type* m_vGroupedVariable; 
        unsigned int* m_vVariableGroupBeginIndex; 
        unsigned int m_numGroups; 
        std::vector<unsigned int> m_vConstraintPartition; 
        coefficient_value_type* m_vLagMultiplier; 
        coefficient_value_type* m_vNewLagMultiplier; 
        coefficient_value_type* m_vSlackness; 
        std::vector<coefficient_value_type> m_vScalingFactor; 

        coefficient_value_type m_objConstant; 
        coefficient_value_type m_lagObj; 
        bool m_lagObjFlag; 

        unsigned int m_iter; 
        unsigned int m_maxIters; 
        bool m_useInitialSol; 

        std::vector<variable_value_type> m_vBestVariableSol; 
        coefficient_value_type m_bestObj; 

    private:
        friend class FeasibleSearcher<coefficient_value_type, variable_value_type>; 
};

template <typename T, typename V>
unsigned int MultiKnapsackLagRelax<T, V>::maxIterations() const 
{
    return m_maxIters;
}
template <typename T, typename V>
void MultiKnapsackLagRelax<T, V>::setMaxIterations(unsigned int maxIters) 
{
    m_maxIters = maxIters; 
}
template <typename T, typename V>
bool MultiKnapsackLagRelax<T, V>::useInitialSolutions() const 
{
    return m_useInitialSol;
}
template <typename T, typename V>
void MultiKnapsackLagRelax<T, V>::setUseInitialSolutions(bool f) 
{
    m_useInitialSol = f; 
}
template <typename T, typename V>
bool MultiKnapsackLagRelax<T, V>::lagObjFlag() const
{
    return m_lagObjFlag;
}
template <typename T, typename V>
void MultiKnapsackLagRelax<T, V>::setLagObjFlag(bool f)
{
    m_lagObjFlag = f; 
}
template <typename T, typename V>
unsigned int MultiKnapsackLagRelax<T, V>::numNegativeSlackConstraints(bool evaluateFlag) 
{
    unsigned int result = 0; 
    if (evaluateFlag)
        computeSlackness();
    for (typename matrix_type::index_type i = 0; i < m_constrMatrix.numRows; ++i)
        result += (m_vSlackness[i] < 0)? 1 : 0; 
    return result; 
}
template <typename T, typename V>
void MultiKnapsackLagRelax<T, V>::copy(MultiKnapsackLagRelax<T, V> const& rhs)
{
    m_model = rhs.m_model; 

    // recycle model 
    destroy(); 
    // copy problem model 
    if (rhs.m_vObjCoef)
    {
        m_vObjCoef = new coefficient_value_type [m_model->numVariables()]; 
        std::copy(rhs.m_vObjCoef, rhs.m_vObjCoef+m_model->numVariables(), m_vObjCoef);
    }
    m_constrMatrix = rhs.m_constrMatrix; 
    if (rhs.m_vConstrRhs)
    {
        m_vConstrRhs = new coefficient_value_type [m_constrMatrix.numRows];
        std::copy(rhs.m_vConstrRhs, rhs.m_vConstrRhs+m_constrMatrix.numRows, m_vConstrRhs);
    }

    // copy variable groups 
    if (rhs.m_vGroupedVariable)
    {
        m_vGroupedVariable = new variable_type [rhs.m_vVariableGroupBeginIndex[rhs.m_numGroups]]; 
        m_numGroups = rhs.m_numGroups;
        m_vVariableGroupBeginIndex = new unsigned int [m_numGroups+1];
        std::copy(rhs.m_vGroupedVariable, rhs.m_vGroupedVariable+rhs.m_vVariableGroupBeginIndex[rhs.m_numGroups], m_vGroupedVariable); 
        std::copy(rhs.m_vVariableGroupBeginIndex, rhs.m_vVariableGroupBeginIndex+rhs.m_numGroups+1, m_vVariableGroupBeginIndex);
    }

    m_vConstraintPartition = rhs.m_vConstraintPartition;
    if (rhs.m_vLagMultiplier)
    {
        m_vLagMultiplier = new coefficient_value_type [m_constrMatrix.numRows]; 
        std::copy(rhs.m_vLagMultiplier, rhs.m_vLagMultiplier+m_constrMatrix.numRows, m_vLagMultiplier); 
        m_vNewLagMultiplier = new coefficient_value_type [m_constrMatrix.numRows];
        std::copy(rhs.m_vNewLagMultiplier, rhs.m_vNewLagMultiplier+m_constrMatrix.numRows, m_vNewLagMultiplier);
        m_vSlackness = new coefficient_value_type [m_constrMatrix.numRows]; 
        std::copy(rhs.m_vSlackness, rhs.m_vSlackness+m_constrMatrix.numRows, m_vSlackness); 
    }
    m_objConstant = rhs.m_objConstant; 
    m_lagObj = rhs.m_lagObj; 
    m_lagObjFlag = rhs.m_lagObjFlag;
    m_iter = rhs.m_iter; 
    m_maxIters = rhs.m_maxIters;
    m_useInitialSol = rhs.m_useInitialSol; 

    m_vBestVariableSol = rhs.m_vBestVariableSol;
    m_bestObj = rhs.m_bestObj;
}
template <typename T, typename V>
void MultiKnapsackLagRelax<T, V>::destroy() 
{
    // recycle variable groups 
    if (m_vGroupedVariable)
    {
        delete [] m_vGroupedVariable; 
        delete [] m_vVariableGroupBeginIndex;
        m_vGroupedVariable = NULL;
        m_vVariableGroupBeginIndex = NULL;
    }
    if (m_vObjCoef)
    {
        delete [] m_vObjCoef; 
        m_vObjCoef = NULL;
    }
    m_constrMatrix.reset();
    if (m_vConstrRhs)
    {
        delete [] m_vConstrRhs;
        m_vConstrRhs = NULL;
    }
    // recycle lagrangian data 
    if (m_vLagMultiplier)
    {
        delete [] m_vLagMultiplier;
        delete [] m_vNewLagMultiplier;
        delete [] m_vSlackness;
        m_vLagMultiplier = NULL;
        m_vNewLagMultiplier = NULL;
        m_vSlackness = NULL;
    }
}
template <typename T, typename V>
SolverProperty MultiKnapsackLagRelax<T, V>::solve(typename MultiKnapsackLagRelax<T, V>::updater_type* updater, typename MultiKnapsackLagRelax<T, V>::scaler_type* scaler, typename MultiKnapsackLagRelax<T, V>::searcher_type* searcher)
{
    bool defaultUpdater = false; 
    bool defaultScaler = false; 
    bool defaultSearcher = false; 
    // use default updater if NULL 
    if (updater == NULL)
    {
        updater = new SubGradientDescent<coefficient_value_type>(); 
        defaultUpdater = true; 
    }
    // use default scaler if NULL 
    if (scaler == NULL)
    {
        scaler = new L2NormScaler<coefficient_value_type, variable_value_type>(); 
        defaultScaler = true; 
    }
    // use default searcher if NULL 
    if (searcher == NULL)
    {
        searcher = new SearchByAdjustCoefficient<coefficient_value_type, variable_value_type>(this); 
        defaultSearcher = true; 
    }

    // recycle old model 
    destroy(); 

    // scale problem 
    scale(scaler); 
    // prepare data structure 
    prepare(); 

    // solve lagrangian subproblem 
    SolverProperty status = solveSubproblems(updater, 0, m_maxIters); 

    // try to find a feasible solution by post refinement 
    status = postProcess(updater, searcher, status); 

    // recover problem from scaling 
    unscale();

    // recycle default updater 
    if (defaultUpdater)
        delete updater; 
    // recycle default scaler 
    if (defaultScaler)
        delete scaler;
    // recycle default searcher 
    if (defaultSearcher)
        delete searcher; 

    return status;
}
template <typename T, typename V>
SolverProperty MultiKnapsackLagRelax<T, V>::solveSubproblems(typename MultiKnapsackLagRelax<T, V>::updater_type* updater, unsigned int beginIter, unsigned int endIter)
{
    // solve lagrangian subproblem 
    SolverProperty status = INFEASIBLE; 
    for (m_iter = beginIter; m_iter < endIter; ++m_iter)
    {
        if (!useInitialSolutions() || m_iter != 0)
            solveLag(); 
        computeSlackness();
#ifdef DEBUG_MULTIKNAPSACKLAGRELAX
        limboPrint(kDEBUG, "iteration %u with %u negative slacks, %g lagrangian objective, %g objective\n", m_iter, numNegativeSlackConstraints(false), evaluateLagObjective(), m_model->evaluateObjective());
        char buf[64];
        limboSPrint(kNONE, buf, "lag%u", m_iter);
        std::ofstream out (buf);
        printObjCoef(out);
        printLagMultiplier(out);
        out.close();
#endif
        if ((status = converge()) == OPTIMAL || m_iter+1 == endIter)
            break; 

        updateLagMultipliers(updater);
    }

    return status; 
}
template <typename T, typename V>
void MultiKnapsackLagRelax<T, V>::scale(typename MultiKnapsackLagRelax<T, V>::scaler_type* scaler)
{
    // compute scaling factors and perform scale 
    m_vScalingFactor.resize(m_model->constraints().size()+1);
    // constraints 
    for (unsigned int i = 0, ie = m_model->constraints().size(); i < ie; ++i)
    {
        m_vScalingFactor[i] = scaler->operator()(m_model->constraints()[i]); 
        m_model->scaleConstraint(i, 1.0/m_vScalingFactor[i]);
    }
    // objective 
    m_vScalingFactor.back() = scaler->operator()(m_model->objective()); 
    m_model->scaleObjective(1.0/m_vScalingFactor.back());
}
template <typename T, typename V>
void MultiKnapsackLagRelax<T, V>::unscale()
{
    // constraints 
    for (unsigned int i = 0, ie = m_model->constraints().size(); i < ie; ++i)
        m_model->scaleConstraint(i, m_vScalingFactor[i]);
    // objective 
    m_model->scaleObjective(m_vScalingFactor.back());
}
template <typename T, typename V>
void MultiKnapsackLagRelax<T, V>::prepare() 
{
    // initialize weights of variables in objective 
    m_vObjCoef = new coefficient_value_type [m_model->numVariables()]; 
    std::fill(m_vObjCoef, m_vObjCoef+m_model->numVariables(), 0);
    for (typename std::vector<term_type>::const_iterator it = m_model->objective().terms().begin(), ite = m_model->objective().terms().end(); it != ite; ++it)
        m_vObjCoef[it->variable().id()] += it->coefficient();

    // for fixed variables, set correct solutions 
    unsigned int i = 0; 
    for (unsigned int ie = m_model->numVariables(); i < ie; ++i)
    {
        variable_type var (i); 
        variable_value_type lowerBound = m_model->variableLowerBound(var); 
        variable_value_type upperBound = m_model->variableUpperBound(var); 
        if (lowerBound == upperBound)
            m_model->setVariableSolution(var, lowerBound); 
    }

    // partition all capacity constraints to the beginning of the array 
    m_vConstraintPartition.resize(m_model->constraints().size()); 
    i = 0; 
    for (unsigned int ie = m_model->constraints().size(); i < ie; ++i)
        m_vConstraintPartition[i] = i; 
    std::vector<unsigned int>::iterator bound = std::partition(m_vConstraintPartition.begin(), m_vConstraintPartition.end(), CapacityConstraintPred(m_model->constraints()));
    unsigned int numCapacityConstraints = std::distance(m_vConstraintPartition.begin(), bound);
    m_vLagMultiplier = new coefficient_value_type [numCapacityConstraints]; 
    std::fill(m_vLagMultiplier, m_vLagMultiplier+numCapacityConstraints, 0);
    m_vNewLagMultiplier = new coefficient_value_type [numCapacityConstraints];
    std::fill(m_vNewLagMultiplier, m_vNewLagMultiplier+numCapacityConstraints, 0);
    m_vSlackness = new coefficient_value_type [numCapacityConstraints]; 
    std::fill(m_vSlackness, m_vSlackness+numCapacityConstraints, 0);

    // change the sense of '>' to '<'
    for (std::vector<unsigned int>::iterator it = m_vConstraintPartition.begin(); it != bound; ++it)
        m_model->constraints().at(*it).normalize('<');

    // initialize constraint matrix 
    m_constrMatrix.numRows = numCapacityConstraints; 
    m_constrMatrix.numColumns = m_model->numVariables(); 
    m_constrMatrix.numElements = 0; 
    m_constrMatrix.vRowBeginIndex = new typename matrix_type::index_type [m_constrMatrix.numRows+1]; 
    m_constrMatrix.vRowBeginIndex[0] = matrix_type::s_startingIndex;

    // initialize vRowBeginIndex 
    i = 1; 
    for (std::vector<unsigned int>::iterator it = m_vConstraintPartition.begin(); it != bound; ++it, ++i)
    {
#ifdef DEBUG_MULTIKNAPSACKLAGRELAX
        limboAssert(i <= m_constrMatrix.numRows); 
#endif
        constraint_type const& constr = m_model->constraints()[*it];
        m_constrMatrix.vRowBeginIndex[i] = m_constrMatrix.vRowBeginIndex[i-1]+constr.expression().terms().size(); 
    }
    // last element of vRowBeginIndex denotes the total number of elements 
    m_constrMatrix.numElements = m_constrMatrix.vRowBeginIndex[m_constrMatrix.numRows]-matrix_type::s_startingIndex; 

    // initialize vElement and vColumn 
    m_constrMatrix.vElement = new coefficient_value_type [m_constrMatrix.numElements]; 
    m_constrMatrix.vColumn = new typename matrix_type::index_type [m_constrMatrix.numElements];
    i = 0; 
    for (std::vector<unsigned int>::iterator it = m_vConstraintPartition.begin(); it != bound; ++it)
    {
        constraint_type const& constr = m_model->constraints()[*it];
        for (typename std::vector<typename constraint_type::term_type>::const_iterator itt = constr.expression().terms().begin(), itte = constr.expression().terms().end(); itt != itte; ++itt, ++i)
        {
#ifdef DEBUG_MULTIKNAPSACKLAGRELAX
            limboAssert(i < m_constrMatrix.numElements); 
#endif
            m_constrMatrix.vElement[i] = itt->coefficient(); 
            m_constrMatrix.vColumn[i] = itt->variable().id()+matrix_type::s_startingIndex;
        }
    }
    // initialize right hand side of constraints 
    m_vConstrRhs = new coefficient_value_type [m_constrMatrix.numRows];
    i = 0; 
    for (std::vector<unsigned int>::iterator it = m_vConstraintPartition.begin(); it != bound; ++it, ++i)
    {
        constraint_type const& constr = m_model->constraints()[*it];
        m_vConstrRhs[i] = constr.rightHandSide();
    }

    // group variables according items 
    // the variables for one item will be grouped 
    // I assume the rest constraints are all single item constraints 
    m_numGroups = std::distance(bound, m_vConstraintPartition.end());
    unsigned int numGroupElements = 0; // it may be smaller than number of variables due to the existence of fixed variables 
    for (std::vector<unsigned int>::iterator it = bound, ite = m_vConstraintPartition.end(); it != ite; ++it, ++i)
    {
        constraint_type const& constr = m_model->constraints()[*it];
        numGroupElements += constr.expression().terms().size();
    }
    m_vGroupedVariable = new variable_type [numGroupElements];
    m_vVariableGroupBeginIndex = new unsigned int [m_numGroups+1]; // append a dummy begin index as the ending index for last group 
    // append a dummy begin index as the ending index for last group 
    m_vVariableGroupBeginIndex[m_numGroups] = numGroupElements; 
    i = 0; // use as group index 
    unsigned int j = 0; // use as group variable index 
    for (std::vector<unsigned int>::iterator it = bound, ite = m_vConstraintPartition.end(); it != ite; ++it, ++i)
    {
        constraint_type const& constr = m_model->constraints()[*it];
#ifdef DEBUG_MULTIKNAPSACKLAGRELAX
        limboAssert(constr.sense() == '='); 
#endif
        expression_type const& expr = constr.expression();
        m_vVariableGroupBeginIndex[i] = j; 
        for (typename std::vector<term_type>::const_iterator itt = expr.terms().begin(), itte = expr.terms().end(); itt != itte; ++itt, ++j)
        {
            m_vGroupedVariable[j] = itt->variable(); 
        }
    }
}
template <typename T, typename V>
void MultiKnapsackLagRelax<T, V>::updateLagMultipliers(typename MultiKnapsackLagRelax<T, V>::updater_type* updater)
{
#if 1
    // update lagrangian multiplier 
    // \lambda^{k+1} = \lambda^{k} + t_k (Ax-b) 
    // vector scaling, the slackness we computed is b-Ax  
    updater->operator()(m_iter, m_constrMatrix.numRows, m_vSlackness, m_vLagMultiplier, m_vNewLagMultiplier); 

    // update objective coefficients
    // c^{k+1, T} = c^{0, T} + \lambda^{k+1, T} (Ax-b) = c^{0, T} + \lambda^{k+1, T} A - \lambda^{k+1, T} b
    // c^{k+1} = c^{0} + A^T \lambda^{k+1} - b^T \lambda^{k+1} = c^{k} + A^T \Delta \lambda^{k+1} - b^T \lambda^{k+1}
    axpy(m_constrMatrix.numRows, (coefficient_value_type)-1, m_vLagMultiplier, m_vNewLagMultiplier); // m_vNewLagMultiplier becomes delta multipliers 
    ATxPlusy((coefficient_value_type)1, m_constrMatrix, m_vNewLagMultiplier, m_vObjCoef); 
    axpy(m_constrMatrix.numRows, (coefficient_value_type)1, m_vNewLagMultiplier, m_vLagMultiplier); // update delta multipliers of m_vNewLagMultiplier to m_vLagMultiplier
#else 
    updater->operator()(m_iter, m_constrMatrix.numRows, m_vSlackness, m_vLagMultiplier, m_vLagMultiplier); 

    std::fill(m_vObjCoef, m_vObjCoef+m_model->numVariables(), 0);
    for (typename std::vector<term_type>::const_iterator it = m_model->objective().terms().begin(), ite = m_model->objective().terms().end(); it != ite; ++it)
        m_vObjCoef[it->variable().id()] += it->coefficient();
    ATxPlusy((coefficient_value_type)1, m_constrMatrix, m_vLagMultiplier, m_vObjCoef); 
#endif

    if (m_lagObjFlag)
    {
        m_objConstant = -dot(m_constrMatrix.numRows, m_vConstrRhs, m_vLagMultiplier);
    }
}
template <typename T, typename V>
void MultiKnapsackLagRelax<T, V>::solveLag()
{
    // for each item 
    CompareVariableByCoefficient helper (m_vObjCoef);
    variable_value_type* vVariableSol = &m_model->variableSolutions()[0];
    unsigned int i = 0; 
    variable_type const* variableGroupBegin = m_vGroupedVariable; 
    variable_type const* variableGroupEnd = m_vGroupedVariable; 
    variable_type const* variable; 

    // reset variables in this group 
    std::fill(vVariableSol, vVariableSol+m_model->numVariables(), 0);
    for (i = 0; i < m_numGroups; ++i)
    {
        variableGroupBegin = variableGroupEnd;
        variableGroupEnd = m_vGroupedVariable+m_vVariableGroupBeginIndex[i+1];
        // find the bin with minimum cost for each item 
        variable = std::min_element(variableGroupBegin, variableGroupEnd, helper);
        vVariableSol[variable->id()] = 1; 
    }
    // evaluate current objective 
    if (m_lagObjFlag)
        m_lagObj = evaluateLagObjective();

}
template <typename T, typename V>
void MultiKnapsackLagRelax<T, V>::computeSlackness() 
{
    // s = b-Ax
    bMinusAx(m_constrMatrix, &m_model->variableSolutions()[0], m_vConstrRhs, m_vSlackness);
}
template <typename T, typename V>
typename MultiKnapsackLagRelax<T, V>::coefficient_value_type MultiKnapsackLagRelax<T, V>::evaluateLagObjective() const 
{
    // evaluate current objective 
    coefficient_value_type objValue = m_objConstant; 
    for (unsigned int i = 0, ie = m_model->numVariables(); i < ie; ++i)
        objValue += m_vObjCoef[i]*m_model->variableSolution(variable_type(i));
    return objValue; 
}
template <typename T, typename V>
SolverProperty MultiKnapsackLagRelax<T, V>::converge()
{
    bool feasibleFlag = true; 
    bool convergeFlag = true; 
    for (typename matrix_type::index_type i = 0; i < m_constrMatrix.numRows; ++i)
    {
        coefficient_value_type multiplier = m_vLagMultiplier[i];
        coefficient_value_type slackness = m_vSlackness[i];
        // must satisfy KKT condition 
        // lambda >= 0
        // g(x) <= 0
        // lambda * g(x) = 0
        //
        // g(x) = Ax-b = -slackness 
        if (slackness < 0)
        {
            convergeFlag = feasibleFlag = false; 
            break; 
        }
        else if (multiplier*slackness != 0)
        {
            convergeFlag = false; 
        }
    }
    SolverProperty status = INFEASIBLE; 
    // store feasible solutions with better objective 
    if (feasibleFlag)
    {
        coefficient_value_type obj = m_model->evaluateObjective(); 
        if (obj < m_bestObj)
        {
            m_vBestVariableSol = m_model->variableSolutions();
            m_bestObj = obj; 
        }
        status = SUBOPTIMAL;
    }
    if (convergeFlag)
        status = OPTIMAL;

    return status;
}
template <typename T, typename V>
SolverProperty MultiKnapsackLagRelax<T, V>::postProcess(typename MultiKnapsackLagRelax<T, V>::updater_type* updater, typename MultiKnapsackLagRelax<T, V>::searcher_type* searcher, SolverProperty status)
{
    if (status == OPTIMAL) // already OPTIMAL
        return status; 
    else if (m_bestObj != std::numeric_limits<coefficient_value_type>::max()) // there is a best feasible solution in store 
    {
        m_model->variableSolutions() = m_vBestVariableSol; 
        return SUBOPTIMAL;
    }
    else // no best solutions in store 
    {
        // try to search feasible solution with heuristic  
        return searcher->operator()(updater); 
    }
}
template <typename T, typename V>
template <typename TT, typename VV>
void MultiKnapsackLagRelax<T, V>::bMinusAx(typename MultiKnapsackLagRelax<T, V>::matrix_type const& A, VV const* x, TT const* b, TT* y) const
{
    // s = b-Ax
    vcopy(A.numRows, b, y);
    AxPlusy((coefficient_value_type)-1, A, x, y);
}
template <typename T, typename V>
bool MultiKnapsackLagRelax<T, V>::printVariableGroup(std::string const& filename) const 
{
    std::ofstream out (filename.c_str()); 
    if (!out.good())
    {
        limbo::limboPrint(kWARN, "failed to open %s for write", filename.c_str()); 
        return false; 
    }
    printVariableGroup(out);
    out.close();
    return true;
}
template <typename T, typename V>
std::ostream& MultiKnapsackLagRelax<T, V>::printVariableGroup(std::ostream& os) const 
{
    os << __func__ << " iteration " << m_iter << "\n";
    for (unsigned int i = 0; i < m_numGroups; ++i)
    {
        os << "[" << i << "]";
        for (unsigned int j = m_vVariableGroupBeginIndex[i]; j < m_vVariableGroupBeginIndex[i+1]; ++j)
        {
            variable_type* variable = m_vGroupedVariable[j]; 
            if (m_model->variableSolution(*variable) == 1)
                os << " *" << m_model->variableName(*variable) << "*";
            else 
                os << " " << m_model->variableName(*variable);
        }
        os << "\n";
    }
    return os; 
}
template <typename T, typename V>
bool MultiKnapsackLagRelax<T, V>::printObjCoef(std::string const& filename) const 
{
    std::ofstream out (filename.c_str()); 
    if (!out.good())
    {
        limbo::limboPrint(kWARN, "failed to open %s for write", filename.c_str()); 
        return false; 
    }
    printObjCoef(out);
    out.close();
    return true;
}
template <typename T, typename V>
std::ostream& MultiKnapsackLagRelax<T, V>::printObjCoef(std::ostream& os) const 
{
    os << __func__ << " iteration " << m_iter << "\n";
    os << "lagrangian objective = " << m_lagObj << "\n";
    os << "objective = " << m_model->evaluateObjective() << "\n";
    for (unsigned int i = 0, ie = m_model->numVariables(); i < ie; ++i)
        os << m_model->variableName(variable_type(i)) << " = " << m_vObjCoef[i] << "\n";
    return os; 
}
template <typename T, typename V>
bool MultiKnapsackLagRelax<T, V>::printLagMultiplier(std::string const& filename) const 
{
    std::ofstream out (filename.c_str()); 
    if (!out.good())
    {
        limbo::limboPrint(kWARN, "failed to open %s for write", filename.c_str()); 
        return false; 
    }
    printLagMultiplier(out);
    out.close();
    return true;
}
template <typename T, typename V>
std::ostream& MultiKnapsackLagRelax<T, V>::printLagMultiplier(std::ostream& os) const 
{
    os << __func__ << " iteration " << m_iter << "\n";
    for (int i = 0; i < m_constrMatrix.numRows; ++i)
        os << "[" << m_model->constraintName(m_model->constraints().at(m_vConstraintPartition[i])) << "] = " << m_vLagMultiplier[i] 
            << " slack = " << m_vSlackness[i]
            << "\n";
    return os; 
}
template <typename T, typename V>
std::ostream& MultiKnapsackLagRelax<T, V>::printConstraint(std::ostream& os, std::string const& name) const 
{
    for (int i = 0; i < this->m_constrMatrix.numRows; ++i)
    {
        constraint_type const& constr = this->m_model->constraints()[i];
        if (m_model->constraintName(constr) == name) 
        {
            os << "constraint " << name << ": rhs = " << constr.rightHandSide() << "\n";
            for (typename std::vector<term_type>::const_iterator it = constr.expression().terms().begin(), ite = constr.expression().terms().end(); it != ite; ++it)
                os << m_model->variableName(it->variable()) << ": " << m_model->variableSolution(it->variable()) << "*" << it->coefficient() << "\n"; 
        }
    }
    return os; 
}
template <typename T, typename V>
template <typename TT>
void MultiKnapsackLagRelax<T, V>::printArray(unsigned int n, TT const* array, bool nonzeroFlag) const
{
    limboPrint(kNONE, "array of length %u = {", n); 
    for (unsigned int i = 0; i < n; ++i)
    {
        if (nonzeroFlag)
        {
            if (array[i] != 0)
            {
                if (i != 0)
                    limboPrint(kNONE, ", ");
                limboPrint(kNONE, "%u:%g", i, (double)array[i]);
            }
        }
        else 
        {
            if (i != 0)
                limboPrint(kNONE, ", ");
            limboPrint(kNONE, "%g", (double)array[i]); 
        }
    }
    limboPrint(kNONE, "}\n"); 
}

template <typename T>
class LagMultiplierUpdater
{
    public:
        typedef T value_type; 

        LagMultiplierUpdater() {}
        virtual ~LagMultiplierUpdater() {}

        virtual value_type operator()(unsigned int iter, value_type multiplier, value_type slackness) = 0;
        virtual void operator()(unsigned int iter, unsigned int n, value_type const* vSlackness, value_type const* vLagMultiplier, value_type* vNewLagMultiplier) = 0; 
};

template <typename T>
class SubGradientDescent : public LagMultiplierUpdater<T>
{
    public:
        typedef LagMultiplierUpdater<T> base_type;
        typedef typename base_type::value_type value_type;

        SubGradientDescent(value_type alpha = 1, value_type beta = 1)
            : SubGradientDescent::base_type()
            , m_alpha(alpha)
            , m_beta(beta)
            , m_iter(std::numeric_limits<unsigned int>::max())
            , m_scalingFactor(1)
        {
        }
        SubGradientDescent(SubGradientDescent const& rhs)
        {
            copy(rhs);
        }
        SubGradientDescent& operator=(SubGradientDescent const& rhs)
        {
            if (this != &rhs)
            {
                this->base_type::operator=(rhs);
                copy(rhs);
            }
            return *this; 
        }
        ~SubGradientDescent()
        {
        }

        value_type operator()(unsigned int iter, value_type multiplier, value_type slackness)
        {
            // compute scaling factor 
            computeScalingFactor(iter);
            value_type result = std::max((value_type)0, multiplier-m_scalingFactor*slackness);
            return result; 
        }
        void operator()(unsigned int iter, unsigned int n, value_type const* vSlackness, value_type const* vLagMultiplier, value_type* vNewLagMultiplier)
        {
            // compute scaling factor 
            computeScalingFactor(iter);
            value_type const* slackness = vSlackness; 
            value_type const* multiplier = vLagMultiplier; 
            value_type* newMultiplier = vNewLagMultiplier;
            for (unsigned int i = 0; i < n; ++i)
            {
                *newMultiplier = std::max((value_type)0, *multiplier-m_scalingFactor*(*slackness));
                ++multiplier;
                ++slackness;
                ++newMultiplier;
            }
        }
    protected:
        void copy(SubGradientDescent const& rhs)
        {
            m_alpha = rhs.m_alpha;
            m_beta = rhs.m_beta;
            m_iter = rhs.m_iter; 
            m_scalingFactor = rhs.m_scalingFactor; 
        }
        void computeScalingFactor(unsigned int iter)
        {
            // avoid frequent computation of scaling factor  
            if (m_iter != iter)
            {
                m_iter = iter; 
                m_scalingFactor = pow((value_type)m_iter+1, -m_alpha)*m_beta;
            }
        }

        value_type m_alpha; 
        value_type m_beta; 
        unsigned int m_iter; 
        value_type m_scalingFactor; 
};

template <typename T, typename V>
class ProblemScaler
{
    public:
        typedef LinearModel<T, V> model_type; 
        typedef typename model_type::coefficient_value_type value_type; 
        typedef typename model_type::expression_type expression_type; 
        typedef typename model_type::constraint_type constraint_type; 
        typedef typename model_type::term_type term_type; 

        ProblemScaler() {}
        virtual ~ProblemScaler() {}

        virtual value_type operator()(expression_type const& /*expr*/) const
        {
            return 1; 
        }
        virtual value_type operator()(constraint_type const& constr) const
        {
            return this->operator()(constr.expression());
        }
};

template <typename T, typename V>
class MinCoefficientScaler : public ProblemScaler<T, V>
{
    public:
        typedef ProblemScaler<T, V> base_type;
        typedef typename base_type::model_type model_type; 
        typedef typename base_type::value_type value_type;
        typedef typename base_type::expression_type expression_type; 
        typedef typename base_type::constraint_type constraint_type; 
        typedef typename base_type::term_type term_type; 

        MinCoefficientScaler(value_type factor = 1) : base_type(), m_scalingFactor(factor) {}
        ~MinCoefficientScaler() {}

        value_type operator()(expression_type const& expr) const
        {
            value_type result = std::numeric_limits<value_type>::max(); 
            for (typename std::vector<term_type>::const_iterator it = expr.terms().begin(), ite = expr.terms().end(); it != ite; ++it)
                result = std::min(result, it->coefficient()); 
            return result*m_scalingFactor;
        }
        value_type operator()(constraint_type const& constr) const
        {
            return this->operator()(constr.expression());
        }

    protected:
        value_type m_scalingFactor; 
};

template <typename T, typename V>
class L2NormScaler : public ProblemScaler<T, V>
{
    public:
        typedef ProblemScaler<T, V> base_type;
        typedef typename base_type::model_type model_type; 
        typedef typename base_type::value_type value_type;
        typedef typename base_type::expression_type expression_type; 
        typedef typename base_type::constraint_type constraint_type; 
        typedef typename base_type::term_type term_type; 

        L2NormScaler() : base_type() {}
        ~L2NormScaler() {}

        value_type operator()(expression_type const& expr) const
        {
            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()*it->coefficient(); 
            return sqrt(result/expr.terms().size()); 
        }
        value_type operator()(constraint_type const& constr) const
        {
            return this->operator()(constr.expression());
        }
};

template <typename T, typename V>
class FeasibleSearcher 
{
    public:
        typedef LinearModel<T, V> model_type; 
        typedef MultiKnapsackLagRelax<T, V> solver_type; 
        typedef typename solver_type::updater_type updater_type;
        typedef typename model_type::coefficient_value_type coefficient_value_type; 
        typedef typename model_type::variable_value_type variable_value_type; 
        typedef typename model_type::variable_type variable_type; 
        typedef typename model_type::expression_type expression_type; 
        typedef typename model_type::constraint_type constraint_type; 
        typedef typename model_type::term_type term_type; 
        typedef typename solver_type::matrix_type matrix_type;

        FeasibleSearcher(solver_type* solver) 
            : m_solver(solver)
            , m_model(solver->m_model)
            , m_vObjCoef(solver->m_vObjCoef)
            , m_constrMatrix(solver->m_constrMatrix)
            , m_vConstrRhs(solver->m_vConstrRhs)
            , m_vGroupedVariable(solver->m_vGroupedVariable)
            , m_vVariableGroupBeginIndex(solver->m_vVariableGroupBeginIndex)
            , m_numGroups(solver->m_numGroups)
            , m_vConstraintPartition(solver->m_vConstraintPartition)
            , m_vLagMultiplier(solver->m_vLagMultiplier)
            , m_vSlackness(solver->m_vSlackness)
            , m_vScalingFactor(solver->m_vScalingFactor)
            , m_objConstant(solver->m_objConstant)
            , m_lagObj(solver->m_lagObj)
            , m_iter(solver->m_iter)
            , m_maxIters(solver->m_maxIters)
            , m_useInitialSol(solver->m_useInitialSol)
            , m_vBestVariableSol(solver->m_vBestVariableSol)
            , m_bestObj(solver->m_bestObj)
        {
        }
        virtual ~FeasibleSearcher() {}

        virtual SolverProperty operator()(updater_type* /*updater*/) {return INFEASIBLE;} 

    protected:
        void computeSlackness()
        {
            m_solver->computeSlackness(); 
        }
        SolverProperty solveSubproblems(updater_type* updater, unsigned int beginIter, unsigned int endIter)
        {
            return m_solver->solveSubproblems(updater, beginIter, endIter);
        }

        solver_type* m_solver; 
        model_type* const& m_model; 

        coefficient_value_type*& m_vObjCoef; 
        matrix_type const& m_constrMatrix; 
        coefficient_value_type* const& m_vConstrRhs; 

        variable_type* const& m_vGroupedVariable; 
        unsigned int* const& m_vVariableGroupBeginIndex; 
        unsigned int const& m_numGroups; 
        std::vector<unsigned int> const& m_vConstraintPartition; 
        coefficient_value_type*& m_vLagMultiplier; 
        coefficient_value_type*& m_vSlackness; 
        std::vector<coefficient_value_type> const& m_vScalingFactor; 
        coefficient_value_type& m_objConstant; 
        coefficient_value_type& m_lagObj; 
        unsigned int& m_iter; 
        unsigned int& m_maxIters; 
        bool& m_useInitialSol; 

        std::vector<variable_value_type>& m_vBestVariableSol; 
        coefficient_value_type& m_bestObj; 
};

template <typename T, typename V>
class SearchByAdjustCoefficient : public FeasibleSearcher<T, V>
{
    public:
        typedef FeasibleSearcher<T, V> base_type; 
        typedef typename base_type::model_type model_type; 
        typedef typename base_type::solver_type solver_type; 
        typedef typename solver_type::updater_type updater_type;
        typedef typename model_type::coefficient_value_type coefficient_value_type; 
        typedef typename model_type::variable_value_type variable_value_type; 
        typedef typename model_type::variable_type variable_type; 
        typedef typename model_type::expression_type expression_type; 
        typedef typename model_type::constraint_type constraint_type; 
        typedef typename model_type::term_type term_type; 
        typedef typename solver_type::matrix_type matrix_type;

        struct VariableMoveCost 
        {
            variable_type variable; 
            coefficient_value_type moveCost; 
            coefficient_value_type capacity; 

            VariableMoveCost(variable_type var, coefficient_value_type mc, coefficient_value_type cap)
                : variable(var)
                , moveCost(mc)
                , capacity(cap)
            {
            }
        };

        struct CompareVariableMoveCost
        {
            bool operator()(VariableMoveCost const& c1, VariableMoveCost const& c2) const 
            {
                return c1.moveCost/(1+c1.capacity) < c2.moveCost/(1+c2.capacity); 
            }
        };

        SearchByAdjustCoefficient(solver_type* solver, coefficient_value_type convergeRatio = 0.1) : base_type(solver), m_convergeRatio(convergeRatio) {}
        ~SearchByAdjustCoefficient() {}

        virtual SolverProperty operator()(updater_type* updater)
        {
            SolverProperty status = INFEASIBLE; 
            // map variable to group 
            this->mapVariable2Group(); 
            // record number of bins with negative slack 
            unsigned int numNegativeSlacksPrev = std::numeric_limits<unsigned int>::max();
            unsigned int numNegativeSlacks = std::numeric_limits<unsigned int>::max();
            // mark whether a variable has been processed 
            std::vector<bool> vVariableProcess (this->m_model->numVariables(), false);
            // record move cost of items 
            std::vector<VariableMoveCost> vVariableMoveCost;
            // large cost 
            coefficient_value_type largeCost = 0; 
            for (coefficient_value_type const* it = this->m_vObjCoef, * ite = this->m_vObjCoef+this->m_model->numVariables(); it != ite; ++it)
                largeCost += (*it > 0)? *it : 0; 
            do 
            {
                vVariableProcess.assign(this->m_model->numVariables(), false);
                // update number of bins with negative slack for previous iteration 
                numNegativeSlacksPrev = numNegativeSlacks; 
                numNegativeSlacks = 0; 
                // go through all bins with negative slack 
                for (typename matrix_type::index_type i = 0; i < this->m_constrMatrix.numRows; ++i)
                {
                    coefficient_value_type slackness = this->m_vSlackness[i];
                    if (slackness < 0)
                    {
                        constraint_type const& constr = this->m_model->constraints()[i];
                        // compute move cost for items in this bin 
                        computeMoveCost(constr, vVariableProcess, vVariableMoveCost);  
                        // try to forbid some items to be assigned to this bin 
                        // sort items by its cost of moving to other bins 
                        std::sort(vVariableMoveCost.begin(), vVariableMoveCost.end(), CompareVariableMoveCost());

                        // the total capacity of items moving out 
                        for (typename std::vector<VariableMoveCost>::const_iterator it = vVariableMoveCost.begin(), ite = vVariableMoveCost.end(); it != ite; ++it)
                        {
                            if (slackness < 0)
                            {
                                // increase the cost of assigning the item to this bin 
                                //this->m_vObjCoef[it->variable.id()] += (it->moveCost+1)*100; 
                                this->m_vObjCoef[it->variable.id()] = largeCost; 
                                // update total capacity of items moving out 
                                slackness += it->capacity; 
                                // mark as processed 
                                vVariableProcess[it->variable.id()] = true; 
                            }
                            else break; 
                        }
                        // compute number of bins with negative slack 
                        ++numNegativeSlacks; 
                    }
                }
                if (numNegativeSlacks == 0)
                    break; 
                // run more iterations with the updated coefficients 
                status = this->solveSubproblems(updater, this->m_iter+1, this->m_iter+this->m_maxIters/2);
            } while (numNegativeSlacks && numNegativeSlacks*(1+m_convergeRatio) < numNegativeSlacksPrev);

            return status; 
        }
    protected:
        void mapVariable2Group()
        {
            // be careful that some variables may not belong to any group 
            // I do not handle this since it should not be a problem if the variables are accessed from constraints 
            this->m_vVariable2Group.assign(this->m_model->numVariables(), std::numeric_limits<unsigned int>::max());
            for (unsigned int i = 0; i < this->m_numGroups; ++i)
            {
                variable_type const* it = this->m_vGroupedVariable+this->m_vVariableGroupBeginIndex[i];
                variable_type const* ite = this->m_vGroupedVariable+this->m_vVariableGroupBeginIndex[i+1];
                for (; it != ite; ++it)
                {
                    this->m_vVariable2Group[it->id()] = i; 
                }
            }
        }
        void computeMoveCost(constraint_type const& constr, std::vector<bool> const& vVariableProcess, std::vector<VariableMoveCost>& vVariableMoveCost) const 
        {
            // extract items that are currently assigned to this bin 
            vVariableMoveCost.clear();
            unsigned int j = 0; 
            for (typename std::vector<term_type>::const_iterator it = constr.expression().terms().begin(), ite = constr.expression().terms().end(); it != ite; ++it, ++j)
            {
                // selected but not processed yet 
                if (this->m_model->variableSolution(it->variable()) && !vVariableProcess[it->variable().id()]) 
                {
                    // compute move cost of moving item to other bins 
                    coefficient_value_type moveCost = std::numeric_limits<coefficient_value_type>::max(); 
                    unsigned int group = this->m_vVariable2Group[it->variable().id()];
                    variable_type const* itv = this->m_vGroupedVariable+this->m_vVariableGroupBeginIndex[group];
                    variable_type const* itve = this->m_vGroupedVariable+this->m_vVariableGroupBeginIndex[group+1];
                    for (; itv != itve; ++itv)
                    {
                        if (it->variable() != *itv) // different variables 
                        {
                            moveCost = std::min(moveCost, this->m_vObjCoef[itv->id()]-this->m_vObjCoef[it->variable().id()]);
                        }
                    }
                    // use the ratio of move cost to capacity 
                    vVariableMoveCost.push_back(VariableMoveCost(it->variable(), moveCost, it->coefficient()));
                }
            }
        }

        std::vector<unsigned int> m_vVariable2Group; 
        coefficient_value_type m_convergeRatio; 
};

template <typename T, typename V>
class SearchByBinSmoothing : public SearchByAdjustCoefficient<T, V>
{
    public:
        typedef SearchByAdjustCoefficient<T, V> base_type; 
        typedef typename base_type::model_type model_type; 
        typedef typename base_type::solver_type solver_type; 
        typedef typename solver_type::updater_type updater_type;
        typedef typename model_type::coefficient_value_type coefficient_value_type; 
        typedef typename model_type::variable_type variable_type; 
        typedef typename model_type::expression_type expression_type; 
        typedef typename model_type::constraint_type constraint_type; 
        typedef typename model_type::term_type term_type; 
        typedef typename solver_type::matrix_type matrix_type;

        struct VariableMoveCost 
        {
            variable_type variable; 
            variable_type targetVariable; 
            coefficient_value_type moveCost; 

            VariableMoveCost(variable_type var, variable_type targetVar, coefficient_value_type mc)
                : variable(var)
                , targetVariable(targetVar)
                , moveCost(mc)
            {
            }
        };

        struct CompareVariableMoveCost
        {
            bool operator()(VariableMoveCost const& c1, VariableMoveCost const& c2) const 
            {
                return c1.moveCost < c2.moveCost; 
            }
        };


        struct CompareConstraintSlack
        {
            coefficient_value_type const* vSlackness; 

            CompareConstraintSlack(coefficient_value_type const* v) : vSlackness(v) {}
            bool operator()(unsigned int i, unsigned int j) const 
            {
                return vSlackness[i] < vSlackness[j];
            }
        };

        SearchByBinSmoothing(solver_type* solver) : base_type(solver) {}
        ~SearchByBinSmoothing() {}

        virtual SolverProperty operator()(updater_type* /*updater*/) 
        {
            SolverProperty status = INFEASIBLE; 
            // map variable to group 
            this->mapVariable2Group();
            // map variable to constraint 
            this->mapVariable2Constraint(); 
            // record move cost of items 
            std::vector<VariableMoveCost> vVariableMoveCost;
            // collect bins with negative slack 
            std::vector<unsigned int> vNegativeSlackConstr; 

            // go through all bins with negative slack 
            for (typename matrix_type::index_type i = 0; i < this->m_constrMatrix.numRows; ++i)
            {
                coefficient_value_type slackness = this->m_vSlackness[i];
                if (slackness < 0)
                    vNegativeSlackConstr.push_back(i);
            }

            // sort bins from small slackness to large 
            std::sort(vNegativeSlackConstr.begin(), vNegativeSlackConstr.end(), CompareConstraintSlack(this->m_vSlackness)); 
                    
            // go through all bins with negative slack 
            for (std::vector<unsigned int>::const_iterator it = vNegativeSlackConstr.begin(), ite = vNegativeSlackConstr.end(); it != ite; ++it)
            {
                constraint_type const& constr = this->m_model->constraints()[*it]; 

                while (this->m_vSlackness[*it] < 0)
                {
                    // compute move cost for items in this bin 
                    vVariableMoveCost.clear(); 
                    computeMoveCost(constr, vVariableMoveCost); 
                    // sort items by its cost of moving to other bins 
                    typename std::vector<VariableMoveCost>::const_iterator itm = std::min_element(vVariableMoveCost.begin(), vVariableMoveCost.end(), CompareVariableMoveCost());

                    // move out at most one item once 
                    // so that we can always get the accurate capacity of all bins 
                    if (itm->moveCost != std::numeric_limits<coefficient_value_type>::max())
                    {
                        // update total capacity of all bins when the item moves out 
                        variable_type const& var = itm->variable;
                        this->m_model->setVariableSolution(var, 0); 
                        for (std::vector<std::pair<unsigned int, unsigned int> >::const_iterator itc = this->m_mVariable2Constr[var.id()].begin(), itce = this->m_mVariable2Constr[var.id()].end(); itc != itce; ++itc)
                        {
                            constraint_type const& curConstr = this->m_model->constraints()[itc->first]; 
                            term_type const& curTerm = curConstr.expression().terms()[itc->second]; 
                            coefficient_value_type& curSlackness = this->m_vSlackness[itc->first]; 

                            curSlackness += curTerm.coefficient(); 
                        }
                        // update capacity of all bins when the item moves into other bins 
                        variable_type const& targetVar = itm->targetVariable; 
                        this->m_model->setVariableSolution(targetVar, 1); 
                        for (std::vector<std::pair<unsigned int, unsigned int> >::const_iterator itc = this->m_mVariable2Constr[targetVar.id()].begin(), itce = this->m_mVariable2Constr[targetVar.id()].end(); itc != itce; ++itc)
                        {
                            constraint_type const& curConstr = this->m_model->constraints()[itc->first]; 
                            term_type const& curTerm = curConstr.expression().terms()[itc->second]; 
                            coefficient_value_type& curSlackness = this->m_vSlackness[itc->first]; 

                            curSlackness -= curTerm.coefficient(); 
                        }
                    }
                    else break; 
                }
            }
            // go through all bins with negative slack and compute status  
            status = SUBOPTIMAL;
            for (std::vector<unsigned int>::const_iterator it = vNegativeSlackConstr.begin(), ite = vNegativeSlackConstr.end(); it != ite; ++it)
            {
                if (this->m_vSlackness[*it] < 0)
                {
                    status = INFEASIBLE;
                    break; 
                }
            }
            return status; 
        }
    protected:
        void mapVariable2Constraint()
        {
            this->m_mVariable2Constr.resize(this->m_model->numVariables());
            for (typename matrix_type::index_type i = 0; i < this->m_constrMatrix.numRows; ++i)
            {
                constraint_type const& constr = this->m_model->constraints()[i];
                unsigned int j = 0; 
                for (typename std::vector<term_type>::const_iterator it = constr.expression().terms().begin(), ite = constr.expression().terms().end(); it != ite; ++it, ++j)
                    this->m_mVariable2Constr[it->variable().id()].push_back(std::make_pair(i, j)); 
            }

            m_vObjCoefOrig.resize(this->m_model->numVariables(), 0); 
            for (typename std::vector<term_type>::const_iterator it = this->m_model->objective().terms().begin(), ite = this->m_model->objective().terms().end(); it != ite; ++it)
                m_vObjCoefOrig[it->variable().id()] = it->coefficient();
        }
        void computeMoveCost(constraint_type const& constr, std::vector<VariableMoveCost>& vVariableMoveCost) const 
        {
            // extract items that are currently assigned to this bin 
            vVariableMoveCost.clear();
            unsigned int j = 0; 
            for (typename std::vector<term_type>::const_iterator it = constr.expression().terms().begin(), ite = constr.expression().terms().end(); it != ite; ++it, ++j)
            {
                variable_type const& var = it->variable(); 
                // selected items 
                if (this->m_model->variableSolution(var)) 
                {
                    // compute move cost of moving item to other bins 
                    coefficient_value_type moveCost = std::numeric_limits<coefficient_value_type>::max(); 
                    variable_type targetVar; 
                    unsigned int group = this->m_vVariable2Group[it->variable().id()];
                    variable_type const* itv = this->m_vGroupedVariable+this->m_vVariableGroupBeginIndex[group];
                    variable_type const* itve = this->m_vGroupedVariable+this->m_vVariableGroupBeginIndex[group+1];
                    for (; itv != itve; ++itv)
                    {
                        if (var != *itv) // different variables 
                        {
                            // find the minimum slackness of target bins 
                            bool enoughCapacityFlag = true; 
                            for (std::vector<std::pair<unsigned int, unsigned int> >::const_iterator itc = this->m_mVariable2Constr[itv->id()].begin(), itce = this->m_mVariable2Constr[itv->id()].end(); itc != itce; ++itc)
                            {
                                coefficient_value_type slackness = this->m_vSlackness[itc->first]; 
                                coefficient_value_type coeff = this->m_model->constraints()[itc->first].expression().terms()[itc->second].coefficient(); 
                                if (slackness < coeff)
                                {
                                    enoughCapacityFlag = false; 
                                    break; 
                                }
                            }
                            if (!enoughCapacityFlag)
                                continue; 
                            coefficient_value_type targetCost = this->m_vObjCoefOrig[itv->id()];
                            if (moveCost > targetCost)
                            {
                                moveCost = targetCost; 
                                targetVar = *itv; 
                            }
                        }
                    }
                    // use the ratio of move cost to capacity 
                    vVariableMoveCost.push_back(VariableMoveCost(var, targetVar, moveCost));
                }
            }
        }

        std::vector<std::vector<std::pair<unsigned int, unsigned int> > > m_mVariable2Constr; 
        std::vector<coefficient_value_type> m_vObjCoefOrig; 
};

template <typename T, typename V>
class SearchByCombinedStrategy : public FeasibleSearcher<T, V>
{
    public:
        typedef FeasibleSearcher<T, V> base_type; 
        typedef typename base_type::model_type model_type; 
        typedef typename base_type::solver_type solver_type; 
        typedef typename solver_type::updater_type updater_type;
        typedef typename model_type::coefficient_value_type coefficient_value_type; 
        typedef typename model_type::variable_value_type variable_value_type; 
        typedef typename model_type::variable_type variable_type; 
        typedef typename model_type::expression_type expression_type; 
        typedef typename model_type::constraint_type constraint_type; 
        typedef typename model_type::term_type term_type; 

        SearchByCombinedStrategy(solver_type* solver, coefficient_value_type convergeRatio = 0.1) 
            : base_type(solver)
            , m_searcherCoeff(solver, convergeRatio)
            , m_searcherSmoothing(solver)
        {
        }
        ~SearchByCombinedStrategy() {}

        virtual SolverProperty operator()(updater_type* updater) 
        {
            SolverProperty status = m_searcherCoeff(updater); 
            if (status == INFEASIBLE)
                status = m_searcherSmoothing(updater);
            return status; 
        }

    protected:
        SearchByAdjustCoefficient<coefficient_value_type, variable_value_type> m_searcherCoeff; 
        SearchByBinSmoothing<coefficient_value_type, variable_value_type> m_searcherSmoothing; 
};

} // namespace solvers 
} // namespace limbo 

#endif