CglLandPSimplex.hpp

Go to the documentation of this file.

// Copyright (C) 2005-2009 Pierre Bonami and others.  All Rights Reserved.
// Author:   Pierre Bonami
//           Tepper School of Business
//           Carnegie Mellon University, Pittsburgh, PA 15213
// Date:     21/07/05
//
// $Id: CglLandPSimplex.hpp 1122 2013-04-06 20:39:53Z stefan $
//
// This code is licensed under the terms of the Eclipse Public License (EPL).
//---------------------------------------------------------------------------
#ifndef CglLandPSimplex_H
#define CglLandPSimplex_H

#include <iostream>
#include <vector>

#include "CglConfig.h"
#include "CglLandP.hpp"

#include "OsiSolverInterface.hpp"
#include "CoinMessage.hpp"
#include "CoinMessageHandler.hpp"
#include "CoinWarmStartBasis.hpp"
#include "CoinPackedMatrix.hpp"

#ifdef COIN_HAS_OSICLP
#include "OsiClpSolverInterface.hpp"
#endif

#include "CglLandPTabRow.hpp"
#include "CglLandPUtils.hpp"
#include "CglLandPMessages.hpp"

//#define APPEND_ROW
#define OLD_COMPUTATION

namespace LAP
{
class DebugData;

class CglLandPSimplex
{
public:
    CglLandPSimplex(const OsiSolverInterface &si,
                    const CglLandP::CachedData &cached,
                    const CglLandP::Parameters &params,
                    Validator &validator);
    ~CglLandPSimplex();
    void cacheUpdate(const CglLandP::CachedData &cached, bool reducedSpace = 0);
    bool resetSolver(const CoinWarmStartBasis * basis);
    bool optimize(int var, OsiRowCut & cut, const CglLandP::CachedData &cached, const CglLandP::Parameters & params);
    bool generateMig(int row, OsiRowCut &cut, const CglLandP::Parameters & params);

    int generateExtraCuts(const CglLandP::CachedData &cached, const CglLandP::Parameters & params);
    int generateExtraCut(int i, const CglLandP::CachedData & cached,
                         const CglLandP::Parameters& params);

    void genThisBasisMigs(const CglLandP::CachedData &cached,
                          const CglLandP::Parameters & params) ;

    int insertAllExtr(OsiCuts & cs, CoinRelFltEq eq);

    void setLogLevel(int level)
    {
        handler_->setLogLevel(level);
    }


    void setSi(OsiSolverInterface *si)
    {
        si_ = si;
#ifdef COIN_HAS_OSICLP
        OsiClpSolverInterface * clpSi = dynamic_cast<OsiClpSolverInterface *>(si_);
        if (clpSi)
        {
            clp_ = clpSi;
        }
#endif
    }
    void freeSi()
    {
        assert(si_ != NULL);
        delete si_;
        si_ = NULL;
#ifdef COIN_HAS_OSICLP
        clp_ = NULL;
#endif
    }

    Cuts& extraCuts()
    {
        return cuts_;
    }
    void loadBasis(const OsiSolverInterface &si,
                   std::vector<int> &M1, std::vector<int> &M2,
                   int k);


    int getNumCols() const
    {
        return ncols_;
    }

    int getNumRows() const
    {
        return nrows_;
    }

    const CoinWarmStartBasis  * getBasis() const
    {
        return basis_;
    }
    const int * getNonBasics() const
    {
        return nonBasics_;
    }

    const int * getBasics() const
    {
        return basics_;
    }

    void outPivInfo(int ncuts)
    {
        handler_->message(RoundStats, messages_)<<ncuts<<numPivots_
        <<numSourceRowEntered_
        <<numIncreased_<<CoinMessageEol;
    }
#ifdef APPEND_ROW

    void append_row(int row_num, bool modularize) ;
    void update_row(TabRow &row);

    void check_mod_row(TabRow &row);
#endif
protected:
    bool changeBasis(int incoming, int leaving, int direction, 
#ifndef OLD_COMPUTATION
                     bool recompute_source_row,
#endif
                     bool modularize);
    int fastFindCutImprovingPivotRow( int &direction, int &gammaSign, double tolerance, bool flagPositiveRows);
    int rescanReducedCosts( int &direction, int &gammaSign, double tolerance);
    int fastFindBestPivotColumn(int direction, int gammaSign,
                                double pivotTol, double rhsTol,
                                bool reducedSpace,
                                bool allowNonImproving,
                                double &bestSigma, bool modularize);

    int findBestPivot(int &leaving, int & direction,
                      const CglLandP::Parameters & params);


    double computeCglpObjective(const TabRow & row, bool modularize = false) const;
    inline double strengthenedIntersectionCutCoef(int i, double alpha_i, double beta) const;
    inline double newRowCoefficient(int j, double gamma) const;
    void createIntersectionCut(TabRow & row, OsiRowCut &cut) const;
    double normalizationFactor(const TabRow & row) const;
    void scaleCut(OsiRowCut & cut, double factor) const;
    //  void createIntersectionCut(double * row);
    void createMIG( TabRow & row, OsiRowCut &cut) const;
    void pullTableauRow(TabRow & row) const;
    void adjustTableauRow(int var, TabRow & row, int direction);
    void resetOriginalTableauRow(int var, TabRow & row, int direction);
    inline double getLoBound(int index) const
    {
        return lo_bounds_[original_index_[index]];
    }
    inline double getUpBound(int index) const
    {
        return up_bounds_[original_index_[index]];
    }
    inline double getColsolToCut(int index) const
    {
        return colsolToCut_[original_index_[index]];
    }
    bool isGtConst(int index) const
    {
        return (index >= ncols_ && lo_bounds_[original_index_[index]] < -1e-10 && up_bounds_[original_index_[index]] <= 1e-09);
    }
    inline void setColsolToCut(int index, double value)
    {
        colsolToCut_[original_index_[index]] = value;
    }
    inline CoinWarmStartBasis::Status getStatus(int index) const
    {
        if (index < ncols_) return basis_->getStructStatus(index);
        return basis_->getArtifStatus(index - ncols_);
    }
    inline bool isInteger(int index) const
    {
        return integers_[original_index_[index]];
    }
    void computeWeights(CglLandP::LHSnorm norm, CglLandP::Normalization type,
                        CglLandP::RhsWeightType rhs);
    double normedCoef(double a, int ii) const
    {
        if (norm_weights_.empty())
        {
            return a;
        }
        else
        {
            return a*norm_weights_[ii];
        }
    }
    void printTableau(std::ostream & os);

  void printEverything();
    void printTableauLateX(std::ostream & os);
    void printRowLateX(std::ostream & os, int i);
    void printCutLateX(std::ostream & os, int i);

    void printCglpBasis(std::ostream& os = std::cout);
    void get_M1_M2_M3(const TabRow & row,
                      std::vector<int> &M1,
                      std::vector<int> &M2,
                      std::vector<int> &M3);
    void eliminate_slacks(double * vec) const;
private:
    CglLandPSimplex();
    CglLandPSimplex(const CglLandPSimplex&);
    CglLandPSimplex& operator=(const CglLandPSimplex&);
#ifdef COIN_HAS_OSICLP

    OsiClpSolverInterface * clp_;
#endif

    void updateM1_M2_M3(TabRow & row, double tolerance, bool alwaysComputeCheap);
    void removeRows(int nDelete, const int * rowsIdx);


    void compute_p_q_r_s(double gamma, int gammaSign, double &p, double & q, double & r , double &s);



    TabRow row_k_;
    TabRow original_row_k_;
    TabRow row_i_;
#ifndef NDBEUG
    TabRow new_row_;
#endif

    CoinPackedVector gammas_;
    std::vector<double> rWk1_;
    std::vector<double> rWk2_;
    std::vector<double> rWk3_;
    std::vector<double> rWk4_;
    std::vector<int> rIntWork_;
    bool * rowFlags_;
    std::vector<bool> col_in_subspace;
    bool *colCandidateToLeave_;
    int * basics_;
    int * nonBasics_;
    std::vector<int> M1_;
    std::vector<int> M2_;
    std::vector<int> M3_;
    double sigma_;
    CoinWarmStartBasis * basis_;
    double * colsolToCut_;
    double * colsol_;
    int ncols_orig_;
    int nrows_orig_;
    int ncols_;
    int nrows_;
    // for fast access to lower bounds (both cols and rows)
    std::vector<double> lo_bounds_;
    // for fast access to upper bounds (both cols and rows)
    std::vector<double> up_bounds_;
    bool inDegenerateSequence_;
    double chosenReducedCostVal_;
    const bool * integers_;
    std::vector<int> original_index_;
    Cuts cuts_;

    OsiSolverInterface * si_;
    bool own_;
    Validator & validator_;
    std::vector<double> norm_weights_;
    double rhs_weight_;

    int nNegativeRcRows_;
    bool checkBasis();

    int numPivots_;
    int numSourceRowEntered_;
    int numIncreased_;

    CoinMessageHandler * handler_;
    CoinMessages messages_;
#ifndef NDEBUG
    double bestSigma_;
#endif

protected:
    double computeCglpRedCost(int direction, int gammaSign, double tau);
    double computeRedCostConstantsInRow();
    double computeCglpObjective(double gamma, bool strengthen, TabRow &row);
    double computeCglpObjective(double gamma, bool strengthen);
    int findCutImprovingPivotRow( int &direction, int &gammaSign, double tolerance);
    int findBestPivotColumn(int direction,
                            double pivotTol, bool reducedSpace, bool allowDegeneratePivot,
                            bool modularize);
#if 1
    int plotCGLPobj(int direction, double gammaTolerance,
                    double pivotTol, bool reducedSpace, bool allowDegenerate, bool modularize);
#endif

};


double CglLandPSimplex::strengthenedIntersectionCutCoef(int i, double alpha_i, double beta) const
{
    //  double ratio = beta/(1-beta);
    if ( (!integers_[i]))
        return intersectionCutCoef(alpha_i, beta);
    else
    {
        double f_i = alpha_i - floor(alpha_i);
        if (f_i < beta)
            return f_i*(1- beta);
        else
            return (1 - f_i)*beta;//(1-beta);
    }
}

double
CglLandPSimplex::newRowCoefficient(int j, double gamma) const
{
    return row_k_[j] + gamma * row_i_[j];
}

}
#endif