// $Id: CglFlowCover.cpp 664 2008-07-27 02:11:59Z ladanyi $
//-----------------------------------------------------------------------------
// name:     Cgl Lifted Simple Generalized Flow Cover Cut Generator
// author:   Yan Xu                email: yan.xu@sas.com
//           Jeff Linderoth        email: jtl3@lehigh.edu
//           Martin Savelsberg     email: martin.savelsbergh@isye.gatech.edu
// date:     05/01/2003
// comments: please scan this file for '???' and read the comments
//-----------------------------------------------------------------------------
// Copyright (C) 2003, Yan Xu, Jeff Linderoth, Martin Savelsberg and others. 
// All Rights Reserved.
// This code is published under the Common Public License.

#if defined(_MSC_VER)
// Turn off MS VS compiler warning about long names
#  pragma warning(disable:4786)
#endif

#include <cstdlib>
#include <cmath>

#include "CoinHelperFunctions.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinPackedVector.hpp"

#include "CglFlowCover.hpp"


#define CGLFLOW_DEBUG 0

//-------------------------------------------------------------------
// Overloaded operator<< for printing VUB and VLB.
//-------------------------------------------------------------------  
std::ostream& operator<<( std::ostream& os, const CglFlowVUB &v ) 
{ 
    os << " VAR = " << v.getVar() << "\t VAL = " << v.getVal() << std::endl; 
    return os; 
}

// Initialize static memeber
int CglFlowCover::numFlowCuts_ = 0;

//-------------------------------------------------------------------
// Determine row types. Find the VUBS and VLBS. 
//-------------------------------------------------------------------  
void 
CglFlowCover::flowPreprocess(const OsiSolverInterface& si) const
{
    CoinPackedMatrix matrixByRow(*si.getMatrixByRow());

    int numRows = si.getNumRows();
    int numCols = si.getNumCols();
  
    const char* sense        = si.getRowSense();
    const double* RHS        = si.getRightHandSide();

    const double* coefByRow  = matrixByRow.getElements();
    const int* colInds       = matrixByRow.getIndices();
    const int* rowStarts     = matrixByRow.getVectorStarts();
    const int* rowLengths    = matrixByRow.getVectorLengths();
    int iRow      = -1; 
    int iCol      = -1;

    numCols_ = numCols;     // Record col and row numbers for copy constructor
    numRows_ = numRows;

    if (rowTypes_ != 0) {
	delete [] rowTypes_; rowTypes_ = 0;
    }
    rowTypes_ = new CglFlowRowType [numRows];// Destructor will free memory
    // Get integer types
    const char * columnType = si.getColType (true);
    
    // Summarize the row type infomation.
    int numUNDEFINED   = 0;
    int numVARUB       = 0;
    int numVARLB       = 0;
    int numVAREQ       = 0;
    int numMIXUB       = 0;
    int numMIXEQ       = 0;
    int numNOBINUB     = 0;
    int numNOBINEQ     = 0;
    int numSUMVARUB    = 0;
    int numSUMVAREQ    = 0;
    int numUNINTERSTED = 0;

    for (iRow = 0; iRow < numRows; ++iRow) {

	int rowLen   = rowLengths[iRow];
	char sen     = sense[iRow];
	double rhs   = RHS[iRow];
	int* ind     = new int [rowLen];
	double* coef = new double [rowLen];


	CoinDisjointCopyN(colInds + rowStarts[iRow], rowLen, ind);
	CoinDisjointCopyN(coefByRow + rowStarts[iRow], rowLen, coef);
 
	CglFlowRowType rowType = determineOneRowType(si, rowLen, ind, coef, 
						     sen, rhs);

	rowTypes_[iRow] = rowType;

	switch(rowType) {
	case  CGLFLOW_ROW_UNDEFINED:
	    ++numUNDEFINED; 
	    break;
	case  CGLFLOW_ROW_VARUB:
	    ++numVARUB; 
	    break;
	case  CGLFLOW_ROW_VARLB:
	    ++numVARLB; 
	    break;
	case  CGLFLOW_ROW_VAREQ:
	    ++numVAREQ; 
	    break;
	case  CGLFLOW_ROW_MIXUB:
	    ++numMIXUB; 
	    break;
	case  CGLFLOW_ROW_MIXEQ:
	    ++numMIXEQ; 
	    break;
	case  CGLFLOW_ROW_NOBINUB:
	    ++numNOBINUB; 
	    break;
	case  CGLFLOW_ROW_NOBINEQ:
	    ++numNOBINEQ; 
	    break;
	case  CGLFLOW_ROW_SUMVARUB:
	    ++numSUMVARUB; 
	    break;
	case  CGLFLOW_ROW_SUMVAREQ:
	    ++numSUMVAREQ; 
	    break;
	case  CGLFLOW_ROW_UNINTERSTED:
	    ++numUNINTERSTED;
	    break;
	default:
	    throw CoinError("Unknown row type", "flowPreprocess",
			    "CglFlowCover");
	}
    
	delete [] ind;  ind  = 0;
	delete [] coef; coef = 0;
    }

#if CGLFLOW_DEBUG
    std::cout << "The num of rows = "  << numRows        << std::endl;
    std::cout << "Summary of Row Type" << std::endl;
    std::cout << "numUNDEFINED     = " << numUNDEFINED   << std::endl;
    std::cout << "numVARUB         = " << numVARUB       << std::endl;
    std::cout << "numVARLB         = " << numVARLB       << std::endl;
    std::cout << "numVAREQ         = " << numVAREQ       << std::endl;
    std::cout << "numMIXUB         = " << numMIXUB       << std::endl;
    std::cout << "numMIXEQ         = " << numMIXEQ       << std::endl;
    std::cout << "numNOBINUB       = " << numNOBINUB     << std::endl;
    std::cout << "numNOBINEQ       = " << numNOBINEQ     << std::endl;
    std::cout << "numSUMVARUB      = " << numSUMVARUB    << std::endl;
    std::cout << "numSUMVAREQ      = " << numSUMVAREQ    << std::endl;
    std::cout << "numUNINTERSTED   = " << numUNINTERSTED << std::endl;
#endif

  //---------------------------------------------------------------------------
  // Setup  vubs_ and vlbs_
    if (vubs_ != 0) { delete [] vubs_; vubs_ = 0; }
    vubs_ = new CglFlowVUB [numCols];      // Destructor will free memory
    if (vlbs_ != 0) { delete [] vlbs_; vlbs_ = 0; }
    vlbs_ = new CglFlowVLB [numCols];      // Destructor will free memory

    for (iCol = 0; iCol < numCols; ++iCol) {   // Initilized in constructor
	vubs_[iCol].setVar(UNDEFINED_);     // but, need redo since may call
	vlbs_[iCol].setVar(UNDEFINED_);     // preprocess(...) more than once
    }
  
    for (iRow = 0; iRow < numRows; ++iRow) {
	
	CglFlowRowType rowType2 = rowTypes_[iRow];
    
	if ( (rowType2 == CGLFLOW_ROW_VARUB) || 
	     (rowType2 == CGLFLOW_ROW_VARLB) || 
	     (rowType2 == CGLFLOW_ROW_VAREQ) )  { 
      
	    int startPos = rowStarts[iRow];
	    int index0   = colInds[startPos];
	    int index1   = colInds[startPos + 1];
	    double coef0 = coefByRow[startPos];
	    double coef1 = coefByRow[startPos + 1];
	    
	    int    xInd,  yInd;   // x is binary
	    double xCoef, yCoef;

	    if ( columnType[index0]==1 ) {
		xInd  = index0;   yInd  = index1;
		xCoef = coef0;    yCoef = coef1;
	    }
	    else {
		xInd  = index1;   yInd  = index0;
		xCoef = coef1;    yCoef = coef0;
	    }

	    switch (rowType2) {
	    case CGLFLOW_ROW_VARUB:       // Inequality: y <= ? * x
		vubs_[yInd].setVar(xInd);
		vubs_[yInd].setVal(-xCoef / yCoef);
		break;
	    case CGLFLOW_ROW_VARLB:       // Inequality: y >= ? * x
		vlbs_[yInd].setVar(xInd);
		vlbs_[yInd].setVal(-xCoef / yCoef);
		break;
	    case CGLFLOW_ROW_VAREQ:       // Inequality: y >= AND <= ? * x
		vubs_[yInd].setVar(xInd);
		vubs_[yInd].setVal(-xCoef / yCoef);
		vlbs_[yInd].setVar(xInd);
		vlbs_[yInd].setVal(-xCoef / yCoef);
		break;
	    default:
		throw CoinError("Unknown row type: impossible", 
				"flowPreprocess", "CglFlowCover");
	    }
	}
    }

#if CGLFLOW_DEBUG
    printVubs(std::cout);
#endif
}


//-----------------------------------------------------------------------------
// Generate LSGFC cuts
//------------------------------------------------------------------- 
void CglFlowCover::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
				const CglTreeInfo info) const 
{
    if (getMaxNumCuts() <= 0) return;
    
    if (getNumFlowCuts() >= getMaxNumCuts()) return;

#if 0
    bool preInit = false;
    bool preReso = false;
    si.getHintParam(OsiDoPresolveInInitial, preInit);
    si.getHintParam(OsiDoPresolveInResolve, preReso);

    if (preInit == false &&  preReso == false) { // Do once
	if (doneInitPre_ == false) {   
	    flowPreprocess(si);
	    doneInitPre_ = true;
	}
    }
    else
#endif
    
    flowPreprocess(si);

    CoinPackedMatrix matrixByRow(*si.getMatrixByRow());
    const char* sense = si.getRowSense();
    const double* rhs = si.getRightHandSide();

    const double* elementByRow = matrixByRow.getElements();
    const int* colInd = matrixByRow.getIndices();
    const CoinBigIndex* rowStart = matrixByRow.getVectorStarts();
    const int* rowLength = matrixByRow.getVectorLengths();
    
    int* ind        = 0;
    double* coef    = 0;
    int iRow, iCol;

    CglFlowRowType rType;

    for (iRow = 0; iRow < numRows_; ++iRow) {
	rType = getRowType(iRow);
	if( ( rType != CGLFLOW_ROW_MIXUB ) &&
	    ( rType != CGLFLOW_ROW_MIXEQ ) &&
	    ( rType != CGLFLOW_ROW_NOBINUB ) &&
	    ( rType != CGLFLOW_ROW_NOBINEQ ) &&
	    ( rType != CGLFLOW_ROW_SUMVARUB ) &&
	    ( rType != CGLFLOW_ROW_SUMVAREQ ) )
	    continue;  

	const int sta = rowStart[iRow];     // Start position of iRow
	const int rowLen = rowLength[iRow]; // iRow length / non-zero elements

	if (ind != 0) { delete [] ind; ind = 0; }
	ind = new int [rowLen];
	if (coef != 0) { delete [] coef; coef = 0; }
	coef = new double [rowLen];

	int lastPos = sta + rowLen;
	for (iCol = sta; iCol < lastPos; ++iCol) {
	    ind[iCol - sta]  = colInd[iCol];
	    coef[iCol - sta] = elementByRow[iCol];
	}

	OsiRowCut flowCut1, flowCut2, flowCut3;
	double violation = 0.0;
	bool hasCut = false;

	if (sense[iRow] == 'E') {
	    hasCut = generateOneFlowCut(si, rowLen, ind, coef, 'L', 
					rhs[iRow], flowCut1, violation);
	    if (hasCut)  {                         // If find a cut
		cs.insert(flowCut1);
		incNumFlowCuts();
		if (getNumFlowCuts() >= getMaxNumCuts())
		    break;
	    }
	    hasCut = false;
	    hasCut = generateOneFlowCut(si, rowLen, ind, coef, 'G', 
					rhs[iRow], flowCut2, violation);
	    if (hasCut)  {
		cs.insert(flowCut2);
		incNumFlowCuts();
		if (getNumFlowCuts() >= getMaxNumCuts())
		    break;
	    }
	}
	if (sense[iRow] == 'L' || sense[iRow] == 'G') {
	    hasCut = generateOneFlowCut(si, rowLen, ind, coef, sense[iRow], 
					rhs[iRow], flowCut3, violation);
	    if (hasCut)  {
		cs.insert(flowCut3);
		incNumFlowCuts();
		if (getNumFlowCuts() >= getMaxNumCuts())
		    break;
	    }
	}
    }


#if CGLFLOW_DEBUG
    std::cout << "\nnumFlowCuts = "<< getNumFlowCuts()  << std::endl;
    std::cout << "CGLFLOW_COL_BINNEG = "<< CGLFLOW_COL_BINNEG  << std::endl;
#endif

    if (ind != 0)  { delete [] ind; ind = 0; }
    if (coef != 0) { delete [] coef; coef = 0; }
}

//-------------------------------------------------------------------
// Default Constructor 
//-------------------------------------------------------------------
CglFlowCover::CglFlowCover()
    :
    CglCutGenerator(),
    maxNumCuts_(2000),
    EPSILON_(1.0e-6),
    UNDEFINED_(-1),
    INFTY_(1.0e30),
    TOLERANCE_(0.05),
    firstProcess_(true),
    numRows_(0),
    numCols_(0),
    doneInitPre_(false),
    vubs_(0),
    vlbs_(0),
    rowTypes_(0)
{ 
    // DO NOTHING
}


//-------------------------------------------------------------------
// Copy constructor 
//-------------------------------------------------------------------
CglFlowCover::CglFlowCover (const CglFlowCover & source)
    :
    CglCutGenerator(source), 
    maxNumCuts_(source.maxNumCuts_),
    EPSILON_(source.EPSILON_),
    UNDEFINED_(source.UNDEFINED_),
    INFTY_(source.INFTY_),
    TOLERANCE_(source.TOLERANCE_),
    firstProcess_(true),
    numRows_(source.numRows_),
    numCols_(source.numCols_),
    doneInitPre_(source.doneInitPre_)
{ 
    setNumFlowCuts(source.numFlowCuts_);
    if (numCols_ > 0) {
	vubs_ = new CglFlowVUB [numCols_];
	vlbs_ = new CglFlowVLB [numCols_];
	CoinDisjointCopyN(source.vubs_, numCols_, vubs_);
	CoinDisjointCopyN(source.vlbs_, numCols_, vlbs_);
    }
    else {
	vubs_ = 0;
	vlbs_ = 0;
    }
    if (numRows_ > 0) {
	rowTypes_ = new CglFlowRowType [numRows_];
	CoinDisjointCopyN(source.rowTypes_, numRows_, rowTypes_);
    }
    else {
	rowTypes_ = 0;
    }
}


//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
CglCutGenerator *
CglFlowCover::clone() const
{
  return new CglFlowCover(*this);
}

//------------------------------------------------------------------
// Assignment operator 
//-------------------------------------------------------------------
CglFlowCover &
CglFlowCover::operator=(const CglFlowCover& rhs)
{
    if (this != &rhs) {
	CglCutGenerator::operator=(rhs);
	maxNumCuts_ = rhs.maxNumCuts_;
	EPSILON_ = rhs.EPSILON_;
	UNDEFINED_ = rhs.UNDEFINED_;
	INFTY_ = rhs.INFTY_;
	TOLERANCE_ = rhs.TOLERANCE_;
	numRows_ = rhs.numRows_;
	numCols_ = rhs.numCols_;
	//    numFlowCuts_ = rhs.numFlowCuts_;
	setNumFlowCuts(rhs.numFlowCuts_);
	doneInitPre_ = rhs.doneInitPre_;
	if (numCols_ > 0) {
	    vubs_ = new CglFlowVUB [numCols_];
	    vlbs_ = new CglFlowVLB [numCols_];
	    CoinDisjointCopyN(rhs.vubs_, numCols_, vubs_);
	    CoinDisjointCopyN(rhs.vlbs_, numCols_, vlbs_);
	}
	if (numRows_ > 0) {
	    rowTypes_ = new CglFlowRowType [numRows_];
	    CoinDisjointCopyN(rhs.rowTypes_, numRows_, rowTypes_);
	}
    }
    return *this;
}


//-------------------------------------------------------------------
// Destructor 
//-------------------------------------------------------------------  
CglFlowCover::~CglFlowCover ()
{
    if (vubs_ != 0) { delete [] vubs_; vubs_ = 0; }
    if (vlbs_ != 0) { delete [] vlbs_; vlbs_ = 0; }
    if (rowTypes_ != 0) { delete [] rowTypes_; rowTypes_ = 0; } 
}


//-------------------------------------------------------------------
//  Given the model data, a row of the model, and a LP solution, 
//  this function tries to generate a violated lifted simple generalized 
//  flow cover.
//-------------------------------------------------------------------  
bool 
CglFlowCover::generateOneFlowCut( const OsiSolverInterface & si, 
				  const int rowLen,
				  int* ind,
				  double* coef,
				  char sense,
				  double rhs,
				  OsiRowCut& flowCut,
				  double& violation ) const
{
    bool generated       = false;
    const double* xlp    = si.getColSolution();
    const int numCols    = si.getNumCols();
    
    double* up           = new double [rowLen];
    double* x            = new double [rowLen];
    double* y            = new double [rowLen];
    CglFlowColType* sign = new CglFlowColType [rowLen];
    
    int i, j;  
    double value, LB, UB;
    
    CglFlowVLB VLB;
    CglFlowVUB VUB;

    // Get integer types
    const char * columnType = si.getColType ();
    for (i = 0; i < rowLen; ++i) {
	if ( xlp[ind[i]] - floor(xlp[ind[i]]) > EPSILON_ && ceil(xlp[ind[i]]) - xlp[ind[i]] > EPSILON_ )
	    break;
    }

    if (i == rowLen)  {
	delete [] sign;
	delete [] up; 
	delete [] x;   
	delete [] y; 
	return generated;
    }

    //-------------------------------------------------------------------------

    if (sense == 'G') flipRow(rowLen, coef, rhs); // flips everything,
						  // but the sense
					  

#if CGLFLOW_DEBUG
    std::cout << "***************************" << std::endl;
    std::cout << "Generate Flow cover -- initial constraint, converted to L sense..." << std::endl;
    std::cout << "Rhs = " << rhs << std::endl;
    std::cout << "coef [var_index]" << " -- " <<  "xlp[var_index]" << '\t' << "vub_coef[vub_index] vub_lp_value OR var_index_col_ub" << std::endl;
   
    for(int iD = 0; iD < rowLen; ++iD) {
	VUB = getVubs(ind[iD]);

	std::cout << coef[iD] << "["  << ind[iD] << "] -- " 
		  << xlp[ind[iD]] << '\t';
	if (VUB.getVar() != UNDEFINED_) {  
	    std::cout << VUB.getVal() << "[" << VUB.getVar() << "]" 
		      << xlp[VUB.getVar()] << std::endl; 
	}
	else
	    std::cout << si.getColUpper()[ind[iD]] << std::endl;
	
    }
#endif

    //-------------------------------------------------------------------------
    // Generate conservation inequality and capacity equalities from 
    // the given row.
  
    for (i = 0; i < rowLen; ++i) {
	
	VLB = getVlbs(ind[i]);
	LB = ( VLB.getVar() != UNDEFINED_ ) ? 
	    VLB.getVal() : si.getColLower()[ind[i]];

	VUB = getVubs(ind[i]);
	UB = ( VUB.getVar() != UNDEFINED_ ) ? 
	    VUB.getVal() : si.getColUpper()[ind[i]];

	if (LB < -EPSILON_) {   // Only consider rows whose variables are all
	    delete [] sign;       // non-negative (LB>= 0). 
	    delete [] up; 
	    delete [] x;   
	    delete [] y;  
	    return generated;     
	}

	if ( columnType[ind[i]]==1 ) {   // Binary variable
	    value = coef[i];
	    if (value > EPSILON_)
		sign[i] = CGLFLOW_COL_BINPOS;
	    else {
		sign[i] = CGLFLOW_COL_BINNEG;
		value = -value;
	    }    
	    up[i] = value;
	    x[i] =  xlp[ind[i]];
	    y[i] = value * x[i];
	}
	else {   
	    value = coef[i];
	    if (value > EPSILON_)
		sign[i] = CGLFLOW_COL_CONTPOS;
	    else {
		sign[i] = CGLFLOW_COL_CONTNEG;
		value = -value;
	    }
	    up[i] = value* UB;
	    x[i] = (VUB.getVar() != UNDEFINED_) ? xlp[VUB.getVar()] : 1.0;
	    y[i] = value * xlp[ind[i]];
	}
    }

    //-------------------------------------------------------------------------
    // Find a initial cover (C+, C-) in (N+, N-)
    double  knapRHS   = rhs;
    double  tempSum   = 0.0;
    double  tempMin   = INFTY_;
    int*    candidate = new int [rowLen];
    int*    label     = new int [rowLen];
    double* ratio     = new double [rowLen];
    int t = -1;
    for (i = 0; i < rowLen; ++i) {
	candidate[i] = label[i] = CGLFLOW_COL_OUTCUT;
	ratio[i] = INFTY_;
	
	switch(sign[i]) {
	case CGLFLOW_COL_CONTPOS:
	case CGLFLOW_COL_BINPOS:
	    if( y[i] > EPSILON_ ) {
		ratio[i] = (1.0 - x[i]) / up[i];
		if( y[i] > up[i] * x[i] - EPSILON_ ) {       // Violated
		    candidate[i] = CGLFLOW_COL_PRIME;
		    tempSum += up[i];
		}
		else {
		    candidate[i] = CGLFLOW_COL_SECONDARY;
		}
	    }
	    break;
	case CGLFLOW_COL_CONTNEG:
	case CGLFLOW_COL_BINNEG:
	    if( up[i] > ( (1.0 - EPSILON_) * INFTY_ ) ) {  // UB is infty
		label[i] = CGLFLOW_COL_INCUT;
	    }
	    else {
		knapRHS += up[i];
		if( y[i] < up[i] ) {
		    candidate[i] = CGLFLOW_COL_PRIME;
		    ratio[i] = x[i] / up[i];
		    tempSum += up[i];
		}
	    }
	    break;
	}    
    }
    
    double diff, tempD, lambda;
    int xID = -1;
    if (knapRHS >1.0e10) {
#if CGLFLOW_DEBUG
      std::cout << "knapsack RHS too large. RETURN." << std::endl; 
#endif
      delete [] sign;                              
      delete [] up; 
      delete [] x;   
      delete [] y;  
      delete [] candidate;
      delete [] label;
      delete [] ratio;
      return generated;
    }

    while (tempSum < knapRHS + EPSILON_) { // Not a cover yet
	diff = INFTY_;
	for (i = 0; i < rowLen; ++i) {
	    if (candidate[i] == CGLFLOW_COL_SECONDARY) {
		tempD = up[i] * x[i] - y[i];
		if (tempD < diff - EPSILON_) {
		    diff = tempD;
		    xID = i;
		}
	    }
	}
    
	if( diff > (1.0 - EPSILON_) * INFTY_  ) {   // NO cover exits.
	    delete [] sign;                              
	    delete [] up; 
	    delete [] x;   
	    delete [] y;  
	    delete [] candidate;
	    delete [] label;
	    delete [] ratio;
	    return generated;
	}
	else {
	    tempSum += up[xID];
	    candidate[xID] = CGLFLOW_COL_PRIME;
	}
    }

    // Solve the knapsack problem to get an initial cover
    tempSum = 0.0;
    for (i = 0; i < rowLen; ++i) {
	if (candidate[i] == CGLFLOW_COL_PRIME && ratio[i] < EPSILON_) {
	    //Zero ratio
	    label[i] = CGLFLOW_COL_INCUT;
	    tempSum += up[i];
	}
    }
  
    while (tempSum < knapRHS + EPSILON_) {
	tempMin = INFTY_;
	xID=-1;
	for (i = 0; i < rowLen; i++) {   // Search the col with  minimum ratio
	    if (candidate[i] == CGLFLOW_COL_PRIME && label[i] == 0 && 
		ratio[i] < tempMin) {
		tempMin = ratio[i];  
		xID = i; 
	    }
	}
	if (xID>=0) {
	  label[xID] = CGLFLOW_COL_INCUT;
	  tempSum += up[xID];
	} else {
#if CGLFLOW_DEBUG
	  std::cout << "knapsack RHS too large B. RETURN." << std::endl; 
#endif
	  delete [] sign;                              
	  delete [] up; 
	  delete [] x;   
	  delete [] y;  
	  delete [] candidate;
	  delete [] label;
	  delete [] ratio;
	  return generated;
	}
    }
  
    // Reduce to a minimal cover
    for (i = 0; i < rowLen; ++i) {
	if (label[i] == CGLFLOW_COL_INCUT && ratio[i] > EPSILON_) {
	    if (tempSum - up[i] > knapRHS + EPSILON_) {
		label[i] = CGLFLOW_COL_OUTCUT;
		tempSum -= up[i];
	    }
	}
    }
    for (i = 0; i < rowLen; ++i) {
	if (label[i] == CGLFLOW_COL_INCUT && ratio[i] < EPSILON_) {
	    if (tempSum - up[i] > knapRHS + EPSILON_) {
		label[i] = CGLFLOW_COL_OUTCUT;
		tempSum -= up[i];
	    }
	}
    }
    
    // Due to the way to handle N-
    for(i = 0; i < rowLen; ++i) {
	if( sign[i] < 0 ) 
	    label[i] = 1 - label[i];
    }

   // No cover, no cut. 
    bool emptyCover = true; 
    for (i = 0; i < rowLen; ++i) {
      if (label[i] == CGLFLOW_COL_INCUT) {
	emptyCover = false; 
	break;
      }
     }
    if (emptyCover) {	
#if CGLFLOW_DEBUG
      std::cout << "No cover. RETURN." << std::endl; 
#endif
      delete [] sign;                              
      delete [] up; 
      delete [] x;   
      delete [] y;  
      delete [] candidate;
      delete [] label;
      delete [] ratio;
      return generated;  
    }

    lambda = tempSum - knapRHS;



#if CGLFLOW_DEBUG
    double sum_mj_Cplus = 0.0;
    double sum_mj_Cminus= 0.0;
    double checkLambda;
    // print out the knapsack variables
    std::cout << "Knapsack Cover: C+" << std::endl;
    for (i = 0; i < rowLen; ++i) { 
	if ( label[i] == CGLFLOW_COL_INCUT && sign[i] > 0 ) {
	  std::cout << ind[i] << '\t' << up[i] << std::endl;
	  sum_mj_Cplus += up[i];
	}
    } 
    std::cout << "Knapsack Cover: C-" << std::endl;
    for (i = 0; i < rowLen; ++i) { 
	if ( label[i] == CGLFLOW_COL_INCUT && sign[i] < 0 ) {
	  std::cout << ind[i] << '\t' << up[i] << std::endl;
	  sum_mj_Cminus += up[i];
	}
    }

    // rlh: verified "lambda" is lambda in the paper.
    // lambda = (sum coefficients in C+) - (sum of VUB
    // coefficients in C-) - rhs-orig-constraint
    std::cout << "lambda = " << lambda << std::endl;
#endif

    //-------------------------------------------------------------------------
    // Generate a violated SGFC

    int numCMinus = 0;
    int numPlusPlus = 0;
    double* rho     = new double [rowLen];
    double* xCoef   = new double [rowLen]; 
    double* yCoef   = new double [rowLen];
    double cutRHS   = rhs;
    double temp     = 0.0;
    double sum      = 0.0;
    double minPlsM  = INFTY_;
    double minNegM  = INFTY_;

    for(i = 0; i < rowLen; ++i) {
	rho[i]   = 0.0;
	xCoef[i] = 0.0;
	yCoef[i] = 0.0;
    }
    
    // Project out variables in C-
    // d^' = d + sum_{i in C^-} m_i. Now cutRHS = d^'
    for (i = 0; i < rowLen; ++i) { 
	if ( label[i] == CGLFLOW_COL_INCUT && sign[i] < 0 ) {
	    cutRHS += up[i];
	    ++numCMinus;
	}
    }

    // (1) Compute the coefficients of the simple generalized flow cover
    // (2) Compute minPlsM, minNegM and sum
    //
    // sum = sum_{i in C+\C++} m_i + sum_{i in L--} m_i = m. Page 15.
    // minPlsM = min_{i in C++} m_i
    // minNegM = min_{i in L-} m_i

    temp = cutRHS;
    
    for (i = 0; i < rowLen; ++i) {
	if (label[i] == CGLFLOW_COL_INCUT  && sign[i] > 0) { // C+
	    yCoef[i] = 1.0;
	    if ( up[i] > lambda + EPSILON_ ) { // C++
		++numPlusPlus;
		xCoef[i] = lambda - up[i];
		cutRHS += xCoef[i];
		if( up[i] < minPlsM ) {
		    minPlsM = up[i];
		}
	    }
	    else {  // C+\C++
	        xCoef[i] = 0.0;  // rlh: is this necesarry? (xCoef initialized to zero)
		sum += up[i];
	    } 
	}
	
	if (label[i] != CGLFLOW_COL_INCUT && sign[i] < 0) { // N-\C-
	    temp += up[i];
	    if ( up[i] > lambda) {      // L-
#if CGLFLOW_DEBUG
	        std::cout << "Variable " << ind[i] << " is in L-" << std::endl;
#endif
		yCoef[i] = 0.0;
		xCoef[i] = -lambda;
		label[i] = CGLFLOW_COL_INLMIN;
		if ( up[i] < minNegM ) { 
		    minNegM = up[i];
		}
	    }
	    else  {        // L--
#if CGLFLOW_DEBUG
	      std::cout << "Variable " << ind[i] << " is in L-- " << std::endl;
#endif
		yCoef[i] = -1.0;
		xCoef[i] = 0.0; // rlh: is this necesarry? (xCoef initialized to zero)
		label[i] = CGLFLOW_COL_INLMINMIN;
		sum += up[i];
	    }
	}
    }
   
    // Sort the upper bounds (m_i) of variables in C++ and L-.

    int     ix;
    int     index  = 0;
    double  temp1  = 0.0;
    double* mt     = new double [rowLen];
    double* M      = new double [rowLen + 1];
    
    while (true) {
	ix = UNDEFINED_;
	value = lambda;
	for (i = 0; i < rowLen; ++i) {
	    if ( (label[i] == CGLFLOW_COL_INCUT && sign[i] > 0) || 
		 label[i] == CGLFLOW_COL_INLMIN ) {     //  C+ || L- 
		if ( up[i] > value ) {       // C++ || L-(up[i] > lambda)
		    ix = i;
		    value = up[i];
		}
	    }
	}
      
	if( ix == UNDEFINED_ )  break;
      
	mt[index++] = up[ix];  // Record m_i in C++ and L-(not all) descently
	
	if( label[ix] == CGLFLOW_COL_INLMIN )  
	    label[ix] = CGLFLOW_COL_INLMINDONE;
	else
	    label[ix] = CGLFLOW_COL_INCUTDONE;
    }
    //printf("mins %g %g\n",minNegM,minPlsM);
    if( index == 0 || numPlusPlus == 0) {
	// No column in C++ and L-(not all). RETURN.
#if CGLFLOW_DEBUG
	std::cout << "index = 0. RETURN." << std::endl; 
#endif
	delete [] sign;
	delete [] up; 
	delete [] x;   
	delete [] y;  
	delete [] candidate;
	delete [] label;
	delete [] ratio;
	delete [] rho;
	delete [] xCoef;
	delete [] yCoef;
	delete [] mt; 
	delete [] M; 
	return generated;
    }

    for ( i = 0; i < rowLen; i++ ) {
	switch( label[i] ) {
	case  CGLFLOW_COL_INCUTDONE:
	    label[i] = CGLFLOW_COL_INCUT;
	    break;
	case  CGLFLOW_COL_INLMIN:
	case  CGLFLOW_COL_INLMINDONE:
	case  CGLFLOW_COL_INLMINMIN:
	    label[i] = 0;
	    break;
	}
    }
    
    temp1 = temp;
    
    /* Get t */
    t = 0;
    for ( i = 0; i < index; ++i ) {
	if ( mt[i] < minPlsM ) {
	    t = i;
	    break;
	} 
    }

    if (i == index) {
	t = index;
    }
    
    /* Compute M_i */
    M[0] = 0.0;
    for ( i = 1; i <= index; ++i ) {
	M[i] = M[(i-1)] + mt[(i-1)];
#if CGLFLOW_DEBUG
	std::cout << "t = " << t << std::endl; 
	printf("mt[%i]=%f, M[%i]=%f\n",i-1, mt[i-1], i, M[i]);
#endif
    }
    // Exit if very big M
    if (M[index]>1.0e30) { // rlh: should test for huge col UB earler 
			   // no sense doing all this work in that case.
#if CGLFLOW_DEBUG
      std::cout << "M[index]>1.0e30. RETURN." << std::endl; 
#endif
      delete [] sign;
      delete [] up; 
      delete [] x;   
      delete [] y;  
      delete [] candidate;
      delete [] label;
      delete [] ratio;
      delete [] rho;
      delete [] xCoef;
      delete [] yCoef;
      delete [] mt; 
      delete [] M; 
      return generated;
    }

    /* Get ml */
    double ml = CoinMin(sum, lambda);
#if CGLFLOW_DEBUG
    // sum = sum_{i in C+\C++} m_i + sum_{i in L--} m_i = m. Page 15.
    std::cout << "ml = min(m, lambda) = min(" << sum << ", " << lambda << ") =" << ml << std::endl; 
#endif
    /* rho_i = max[0, m_i - (minPlsM - lamda) - ml */
    if (t < index ) { /* rho exits only for t <= index-1 */
	value = (minPlsM - lambda) + ml;
	for (i = t; i < index; ++i) {
	    rho[i] =  CoinMax(0.0, mt[i] - value);
#if CGLFLOW_DEBUG
	    printf("rho[%i]=%f\n", i, rho[i]);
#endif
	}
    }

#if 1 // LIFTING
    double estY, estX;
    double movement = 0.0;
    double dPrimePrime = temp + cutRHS; 
    int lifted = false;
    for( i = 0; i < rowLen; ++i ) {
	if ( (label[i] != CGLFLOW_COL_INCUT) && (sign[i] > 0) ) {/* N+\C+*/
	    lifted = liftPlus(estY, estX,
                              index, up[i],
                              lambda,
                              y[i], x[i], 
                              dPrimePrime, M);
	    
	    xCoef[i] = -estX;
	    yCoef[i] = estY;
#if CGLFLOW_DEBUG
	    if (lifted == true) {
		printf("Success: Lifted col %i (up_i=%f,yCoef[i]=%f,xCoef[i]=%f) in N+\\C+\n", 
		       ind[i], up[i], yCoef[i], xCoef[i]);
	    }
	    else {
		printf("Failed to Lift col %i (m_i=%f) in N+\\C+\n", 
		       ind[i], up[i]);
	    }       
#endif
	}
	if (label[i] == CGLFLOW_COL_INCUT && sign[i] < 0) { 
	    /* C- */
	    liftMinus(movement, t,
		      index, up[i], 
		      dPrimePrime, 
		      lambda, ml,
		      M, rho);
                
	    if(movement > EPSILON_) {
#if CGLFLOW_DEBUG
		printf("Success: Lifted col %i in C-, movement=%f\n", 
		       ind[i], movement);
#endif
		lifted = true;
		xCoef[i] = -movement;
		cutRHS -= movement;
	    }
	    else {
#if CGLFLOW_DEBUG
		printf("Failed to Lift col %i in C-, g=%f\n",
		       ind[i], movement);
#endif
	    }
	}
    }
#endif
    //-------------------------------------------------------------------

    
    // Calculate the violation
    violation = -cutRHS;
    for ( i = 0; i < rowLen; ++i ) {
#if CGLFLOW_DEBUG
	std::cout << "i = " << i << " ind = " << ind[i] << " sign = " 
		  << sign[i] 
		  << " coef = " << coef[i] << " x = " << x[i] << " xCoef = " 
		  << xCoef[i] << " y = " << y[i] << " yCoef = " << yCoef[i] 
		  << " up = " << up[i] << " label = " << label[i] << std::endl;
#endif

	violation += y[i] * yCoef[i] + x[i] * xCoef[i];
    }

#if CGLFLOW_DEBUG
    std::cout << "violation = " << violation << std::endl;
#endif
    
    int     cutLen     = 0;
    char    cutSense   = 'L';
    int*    cutInd     = 0;
    double* cutCoef    = 0;

    // If violated, transform the inequality back to original system
    if ( violation > TOLERANCE_ ) {
	cutLen = 0;
	cutSense = 'L';
	cutInd  = new int [numCols];
	cutCoef = new double [numCols];
      
	for ( i = 0; i < rowLen; ++i )  {
	    VUB = getVubs(ind[i]);
      
	    if ( ( sign[i] == CGLFLOW_COL_CONTPOS ) || 
		 ( sign[i] == CGLFLOW_COL_CONTNEG ) ) {

		if ( fabs( yCoef[i] ) > EPSILON_ ) {
		    
		    if ( sign[i] == CGLFLOW_COL_CONTPOS ) 
			cutCoef[cutLen] = coef[i] * yCoef[i];
		    else 
			cutCoef[cutLen] = -coef[i] * yCoef[i];
		  
		    cutInd[cutLen++] = ind[i];
		}

		if ( fabs( xCoef[i] ) > EPSILON_ ) {
		    if ( VUB.getVar() != UNDEFINED_ ) {
			cutCoef[cutLen] = xCoef[i];
			cutInd[cutLen++] = VUB.getVar();
		    }
		    else
			cutRHS -= xCoef[i];
		}
	    }
            
	    if ( ( sign[i] == CGLFLOW_COL_BINPOS ) || 
		 ( sign[i] == CGLFLOW_COL_BINNEG ) ) {
		if (fabs(yCoef[i]) > EPSILON_ || fabs(xCoef[i]) > EPSILON_) {
		    if (sign[i] == CGLFLOW_COL_BINPOS) 
			cutCoef[cutLen] = coef[i] * yCoef[i] + xCoef[i];
		    else 
			cutCoef[cutLen] = -coef[i] * yCoef[i] + xCoef[i];
		    cutInd[cutLen++] = ind[i];
		}
	    }
	}
    
	for ( i = 0; i < cutLen; ++i ) {
	    for ( j = 0; j < i; j++ ) {
		if ( cutInd[j] == cutInd[i] ) { /* Duplicate*/
		    cutCoef[j] += cutCoef[i];
		    cutInd[i] = -1;
		}
	    }
	}

	for ( j = 0, i = 0; i < cutLen; ++i ) {
	    if ( ( cutInd[i] == -1 ) || ( fabs( cutCoef[i]) < EPSILON_ ) ){
		/* Small coeff*/
	    }
	    else {
		cutCoef[j] = cutCoef[i];
		cutInd[j] = cutInd[i];
		j++;
	    }
	}
    
	cutLen = j;
	// Skip if no elements ? - bug somewhere
	assert (cutLen);
        
	// Recheck the violation.
	violation = 0.0;
	for (i = 0; i < cutLen; ++i) 
	    violation += cutCoef[i] * xlp[cutInd[i]];
    
	violation -= cutRHS;

	if ( violation > TOLERANCE_ ) {
	    flowCut.setRow(cutLen, cutInd, cutCoef);
	    flowCut.setLb(-1.0 * si.getInfinity());
	    flowCut.setUb(cutRHS);
	    flowCut.setEffectiveness(violation);
	    generated = true;

#if CGLFLOW_DEBUG
	    std::cout << "generateOneFlowCover(): Found a cut" << std::endl;
#endif
	}
	else {
#if CGLFLOW_DEBUG
	    std::cout << "generateOneFlowCover(): Lost a cut" << std::endl;
#endif
	}
    }

    //-------------------------------------------------------------------------
    delete [] sign;
    delete [] up; 
    delete [] x;   
    delete [] y;  
    delete [] candidate;
    delete [] label;
    delete [] ratio;
    delete [] rho;
    delete [] xCoef;
    delete [] yCoef;
    delete [] mt; 
    delete [] M; 
    delete [] cutInd;
    delete [] cutCoef;
    
    return generated;
}


//-------------------------------------------------------------------
// Flip a row from ">=" to "<=", and vice versa. 
//-------------------------------------------------------------------
void 
CglFlowCover::flipRow(int rowLen, double* coef, double& rhs) const
{
    for(int i = 0; i < rowLen; ++i) coef[i] = -coef[i]; 
    rhs = -rhs;  
}

//-------------------------------------------------------------------
// Flip a row from ">=" to "<=", and vice versa. Have 'sense'.
//-------------------------------------------------------------------
void 
CglFlowCover::flipRow(int rowLen, double* coef, char& sen,
		      double& rhs) const
{
    for(int i = 0; i < rowLen; ++i) coef[i] = -coef[i]; 
    sen = (sen == 'G') ?  'L' : 'G';
    rhs = -rhs;  
}

//-------------------------------------------------------------------
// Determine the type of a given row 
//-------------------------------------------------------------------
CglFlowRowType
CglFlowCover::determineOneRowType(const OsiSolverInterface& si,
				  int rowLen, int* ind, 
				  double* coef, char sense, 
				  double rhs) const
{
    if (rowLen == 0) 
	return CGLFLOW_ROW_UNDEFINED;
    
    CglFlowRowType rowType = CGLFLOW_ROW_UNDEFINED;
    // Get integer types
    const char * columnType = si.getColType ();
    
    int  numPosBin = 0;      // num of positive binary variables
    int  numNegBin = 0;      // num of negative binary variables
    int  numBin    = 0;      // num of binary variables
    int  numPosCol = 0;      // num of positive variables
    int  numNegCol = 0;      // num of negative variables
    int  i;
    bool flipped = false;

    // Range row will only consider as 'L'
    if (sense == 'G') {        // Transform to " <= "
	flipRow(rowLen, coef, sense, rhs);                
	flipped = true;
    }
    
    // Summarize the variable types of the given row.
    for ( i = 0; i < rowLen; ++i ) {
	if ( coef[i] < -EPSILON_ ) {
	    ++numNegCol;
	    if( columnType[ind[i]]==1 )
		++numNegBin;
	}
	else {
	    ++numPosCol;
	    if( columnType[ind[i]]==1 )
		++numPosBin;    
	}
    }
    numBin = numNegBin + numPosBin;

#if CGLFLOW_DEBUG
    std::cout << "numNegBin = " << numNegBin << std::endl;
    std::cout << "numPosBin = " << numPosBin << std::endl;
    std::cout << "numBin = " << numBin << std::endl;
    std::cout << "rowLen = " << rowLen << std::endl;
#endif
    
    
    //------------------------------------------------------------------------
    // Classify row type based on the types of variables.
    
    // All variables are binary. NOT interested in this type of row right now
    if (numBin == rowLen) 
	rowType = CGLFLOW_ROW_UNINTERSTED;

    // All variables are NOT binary
    if (rowType == CGLFLOW_ROW_UNDEFINED && numBin == 0) {
	if (sense == 'L')
	    rowType = CGLFLOW_ROW_NOBINUB;
	else 
	    rowType = CGLFLOW_ROW_NOBINEQ;
    }

    // There are binary and other types of variables   
    if (rowType == CGLFLOW_ROW_UNDEFINED) {  
	if ((rhs < -EPSILON_) || (rhs > EPSILON_) || (numBin != 1)) {
	    if (sense == 'L')
		rowType = CGLFLOW_ROW_MIXUB;
	    else 
		rowType = CGLFLOW_ROW_MIXEQ;
	}
	else {                               // EXACTLY one binary
	    if (rowLen == 2) {               // One binary and one other type
		if (sense == 'L') {
		    if (numNegCol == 1 && numNegBin == 1)
			rowType = CGLFLOW_ROW_VARUB;
		    if (numPosCol == 1 && numPosBin == 1)
			rowType = CGLFLOW_ROW_VARLB;
		}
		else
		    rowType = CGLFLOW_ROW_VAREQ;
	    }
	    else {               // One binary and 2 or more other types
		if (numNegCol==1 && numNegBin==1) {// Binary has neg coef and 
		    if (sense == 'L')  // other are positive
			rowType = CGLFLOW_ROW_SUMVARUB;
		    else
			rowType = CGLFLOW_ROW_SUMVAREQ;
		}
	    }
	}
    }
  
    // Still undefined
    if (rowType == CGLFLOW_ROW_UNDEFINED) {
	if (sense == 'L') 
	    rowType = CGLFLOW_ROW_MIXUB;
	else
	    rowType = CGLFLOW_ROW_MIXEQ;
    }
    if (flipped == true) {
	flipRow(rowLen, coef, sense, rhs);                
    }

    return rowType;
}

/*===========================================================================*/

void
CglFlowCover::liftMinus(double &movement, /* Output */ 
			int t,
			int r,
			double z,
			double dPrimePrime, 
			double lambda,
			double ml,
			double *M,
			double *rho) const
{
    int i;
    movement = 0.0;
    
    if (z > dPrimePrime) {
        movement = z - M[r] + r * lambda;
    }
    else {
        for (i = 0; i < t; ++i) {
            if ( (z >= M[i]) && (z <= M[(i+1)] - lambda) ) {
		movement = i * lambda;
		return;
            }
        }
        
        for (i = 1; i < t; ++i) {
            if ( (z >= M[i] - lambda) && (z <= M[i]) ) {
                movement = z - M[i] + i * lambda;
                return;
            }
        }
        
        for (i = t; i < r; ++i) {
            if ( (z >= M[i] - lambda) && (z <= M[i] - lambda + ml + rho[i]) ) {
                movement = z - M[i] + i * lambda;
                return;
            }
        }
        
        for (i = t; i < r; ++i) {
            if ( (z >= M[i]-lambda+ml+rho[i]) && (z <= M[(i+1)]-lambda) ) {
                movement = i * lambda;
                return;
            }
        }
    
        if ((z >= M[r] - lambda) && z <= dPrimePrime) {
            movement = z - M[r] + r * lambda;
        }
    }
    
}

/*===========================================================================*/

int
CglFlowCover::liftPlus(double &alpha, 
			     double &beta,
			     int r,
			     double m_j, 
			     double lambda,
			     double y_j,
			     double x_j,
			     double dPrimePrime,
			     double *M) const
{
    int i;
    int status = false;  /* Defalut: fail to lift */
    double value;
    alpha = 0.0;
    beta = 0.0;
    
    if (m_j > M[r] - lambda + EPSILON_) {
        if (m_j < dPrimePrime - EPSILON_) {
            if ((m_j > (M[r] - lambda)) && (m_j <= M[r])){ /* FIXME: Test */
                value = y_j - x_j * (M[r] - r * lambda);
                
                /* FIXME: Is this "if" useful */
                if (value > 0.0) {
                    status = true;
                    alpha = 1.0;
                    beta = M[r] - r * lambda;
#if CGLFLOW_DEBUG
                    printf("liftPlus:1: value=%f, alpah=%f, beta=%f\n",
                           value, alpha,beta);
#endif
                }
                else {
#if CGLFLOW_DEBUG
                    printf("liftPlus:1: value=%f, become worst\n",value);
#endif
                }
            }
        }
	else {    
#if CGLFLOW_DEBUG
	    printf("liftPlus:1: too big number\n");
#endif
	}
    }
    else {
        for (i = 1; i <= r; ++i) {
            if ((m_j > (M[i] - lambda)) && (m_j <= M[i])){ /* FIXME: Test */
                
                value = y_j - x_j * (M[i] - i * lambda);

                /* FIXME: Is this "if" useful */
                if (value > 0.0) {
                    status = true;
                    alpha = 1.0;
                    beta = M[i] - i * lambda;
#if CGLFLOW_DEBUG
                    printf("liftPlus:2: value=%f, alpah=%f, beta=%f\n",
                           value, alpha, beta);
#endif
                }
                else {
#if CGLFLOW_DEBUG
                    printf("liftPlus:2: value=%f, become worst\n",value);
#endif
                }
                return status;
            }
        }
    }
    
    return status;
}
// Create C++ lines to get to current state
std::string
CglFlowCover::generateCpp( FILE * fp) 
{
  CglFlowCover other;
  fprintf(fp,"0#include \"CglFlowCover.hpp\"\n");
  fprintf(fp,"3  CglFlowCover flowCover;\n");
  if (maxNumCuts_!=other.maxNumCuts_)
    fprintf(fp,"3  flowCover.setMaxNumCuts(%d);\n",maxNumCuts_);
  else
    fprintf(fp,"4  flowCover.setMaxNumCuts(%d);\n",maxNumCuts_);
  if (getAggressiveness()!=other.getAggressiveness())
    fprintf(fp,"3  flowCover.setAggressiveness(%d);\n",getAggressiveness());
  else
    fprintf(fp,"4  flowCover.setAggressiveness(%d);\n",getAggressiveness());
  return "flowCover";
}

