text: // 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:     07/21/05
  //
  // $Id: CglLandP.cpp 1123 2013-04-06 20:47:24Z stefan $
  //
  // This code is licensed under the terms of the Eclipse Public License (EPL).
  //---------------------------------------------------------------------------
  #include "CglLandP.hpp"
  #include "CglLandPSimplex.hpp"
  
  #define INT_INFEAS(value) fabs(value - floor(value+0.5))
  
  #include "CglConfig.h"
  
  #ifdef COIN_HAS_OSICLP
  #include "OsiClpSolverInterface.hpp"
  #endif
  
  #define CLONE_SI //Solver is cloned between two cuts
  
  #include "CoinTime.hpp"
  #include "CglGomory.hpp"
  #include "CoinFactorization.hpp"
  #include <fstream>
  namespace LAP
  {
  //Setup output messages
  LapMessages::LapMessages( )
          :CoinMessages(LAP_MESSAGES_DUMMY_END)
  {
      strcpy(source_,"Lap");
      addMessage(BEGIN_ROUND,CoinOneMessage( 1, 2,"Starting %s round %d variable considered for separation."));
      addMessage(END_ROUND,CoinOneMessage(2, 2,"End ouf %s round %d cut generated in %g seconds."));
      addMessage(DURING_SEP,CoinOneMessage(3,1,"After %g seconds, separated %d cuts."));
      addMessage(CUT_REJECTED, CoinOneMessage(4,1,"Cut rejected for %s."));
      addMessage(CUT_FAILED,CoinOneMessage(5,1,"Generation failed."));
      addMessage(CUT_GAP, CoinOneMessage(7,1,"CUTGAP after %i pass objective is %g"));
      addMessage(LAP_CUT_FAILED_DO_MIG, CoinOneMessage(3006,1,"Failed to generate a cut generate a Gomory cut instead"));
  }
  }
  using namespace LAP;
  CglLandP::Parameters::Parameters():
          CglParam(),
          pivotLimit(20),
          pivotLimitInTree(10),
          maxCutPerRound(5000),
          failedPivotLimit(1),
          degeneratePivotLimit(0),
          extraCutsLimit(5),
          pivotTol(1e-4),
          away(5e-4),
          timeLimit(COIN_DBL_MAX),
          singleCutTimeLimit(COIN_DBL_MAX),
          rhsWeight(1.),
          useTableauRow(true),
          modularize(false),
          strengthen(true),
          countMistakenRc(false),
          sepSpace(Fractional),
          perturb(true),
          normalization(Unweighted),
          rhsWeightType(Fixed),
          lhs_norm(L1),
          generateExtraCuts(none),
          pivotSelection(mostNegativeRc)
  {
      EPS = 1e-08;
  }
  
  CglLandP::Parameters::Parameters(const Parameters &other):
          CglParam(other),
          pivotLimit(other.pivotLimit),
          pivotLimitInTree(other.pivotLimitInTree),
          maxCutPerRound(other.maxCutPerRound),
          failedPivotLimit(other.failedPivotLimit),
          degeneratePivotLimit(other.degeneratePivotLimit),
          extraCutsLimit(other.extraCutsLimit),
          pivotTol(other.pivotTol),
          away(other.away),
          timeLimit(other.timeLimit),
          singleCutTimeLimit(other.singleCutTimeLimit),
          rhsWeight(other.rhsWeight),
          useTableauRow(other.useTableauRow),
          modularize(other.modularize),
          strengthen(other.strengthen),
          countMistakenRc(other.countMistakenRc),
          sepSpace(other.sepSpace),
          perturb(other.perturb),
          normalization(other.normalization),
          rhsWeightType(other.rhsWeightType),
          lhs_norm(other.lhs_norm),
          generateExtraCuts(other.generateExtraCuts),
          pivotSelection(other.pivotSelection)
  {}
  
  CglLandP::Parameters & CglLandP::Parameters::operator=(const Parameters &other)
  {
      if (this != &other)
      {
          CglParam::operator=(other);
          pivotLimit = other.pivotLimit;
          pivotLimitInTree = other.pivotLimitInTree;
          maxCutPerRound = other.maxCutPerRound;
          failedPivotLimit = other.failedPivotLimit;
          degeneratePivotLimit = other.failedPivotLimit;
          extraCutsLimit = other.extraCutsLimit;
          pivotTol = other.pivotTol;
          away = other.away;
          timeLimit = other.timeLimit;
          singleCutTimeLimit = other.singleCutTimeLimit;
          rhsWeight = other.rhsWeight;
          useTableauRow = other.useTableauRow;
          modularize = other.modularize;
          strengthen = other.strengthen;
          countMistakenRc = other.countMistakenRc;
          sepSpace = other.sepSpace;
          perturb = other.perturb;
          normalization = other.normalization;
          rhsWeightType = other.rhsWeightType;
          lhs_norm = other.lhs_norm;
          generateExtraCuts = other.generateExtraCuts;
          pivotSelection = other.pivotSelection;
      }
      return *this;
  }
  
  CglLandP::CachedData::CachedData(int nBasics, int nNonBasics):
          basics_(NULL), nonBasics_(NULL), nBasics_(nBasics),
          nNonBasics_(nNonBasics), basis_(NULL), colsol_(NULL),
          slacks_(NULL), integers_(NULL), solver_(NULL)
  {
      if (nBasics_>0)
      {
          basics_ = new int[nBasics_];
          integers_ = new bool [nNonBasics_ + nBasics_];
      }
      if (nNonBasics_>0)
          nonBasics_ = new int[nNonBasics_];
      if (nBasics_ + nNonBasics_ > 0)
      {
          colsol_ = new double[nBasics_ + nNonBasics_];
          slacks_ = &colsol_[nNonBasics_];
      }
  }
  
  CglLandP::CachedData::CachedData(const CachedData &source):
          basics_(NULL), nonBasics_(NULL), nBasics_(source.nBasics_),
          nNonBasics_(source.nNonBasics_), basis_(NULL),
          colsol_(NULL), slacks_(NULL), integers_(NULL), solver_(NULL)
  {
      if (nBasics_>0)
      {
          basics_ = new int[nBasics_];
          CoinCopyN(source.basics_, nBasics_, basics_);
          integers_ = new bool [nNonBasics_ + nBasics_];
          CoinCopyN(source.integers_, nBasics_ + nNonBasics_, integers_);
      }
      if (nNonBasics_>0)
      {
          nonBasics_ = new int[nNonBasics_];
          CoinCopyN(source.nonBasics_, nBasics_, nonBasics_);
      }
      if (nBasics_ + nNonBasics_ > 0)
      {
          colsol_ = new double[nBasics_ + nNonBasics_];
          slacks_ = &colsol_[nNonBasics_];
          CoinCopyN(source.colsol_, nBasics_ + nNonBasics_, colsol_);
      }
      if (source.basis_!=NULL)
          basis_ = new CoinWarmStartBasis(*source.basis_);
      if (source.solver_!=NULL)
        solver_ = source.solver_->clone();
  }
  
  CglLandP::CachedData& CglLandP::CachedData::operator=(const CachedData &source)
  {
      if (this != &source)
      {
          nBasics_ = source.nBasics_;
          nNonBasics_ = source.nNonBasics_;
          if (basics_ == NULL) delete [] basics_;
          basics_ = NULL;
          if (nonBasics_ == NULL) delete [] nonBasics_;
          nonBasics_ = NULL;
          if (basis_ == NULL) delete [] basis_;
          basis_ = NULL;
          if (colsol_ == NULL) delete [] colsol_;
          colsol_ = NULL;
          if (slacks_ == NULL) delete [] slacks_;
          slacks_ = NULL;
          if (integers_ == NULL) delete [] integers_;
          integers_ = NULL;
          if (nBasics_>0)
          {
              basics_ = new int[nBasics_];
              CoinCopyN(source.basics_, nBasics_, basics_);
              integers_ = new bool [nBasics_ + nNonBasics_];
              CoinCopyN(source.integers_, nBasics_ + nNonBasics_, integers_);
          }
          if (nNonBasics_>0)
          {
              nonBasics_ = new int[nNonBasics_];
              CoinCopyN(source.nonBasics_, nBasics_, nonBasics_);
          }
          if (nBasics_ + nNonBasics_ > 0)
          {
              colsol_ = new double[nBasics_ + nNonBasics_];
              slacks_ = &colsol_[nNonBasics_];
              CoinCopyN(source.colsol_, nBasics_ + nNonBasics_, colsol_);
          }
          if (source.basis_!=NULL)
              basis_ = new CoinWarmStartBasis(*source.basis_);
          delete solver_;
  	if (source.solver_)
  	  solver_ = source.solver_->clone();
      }
      return *this;
  }
  
  void
  CglLandP::CachedData::getData(const OsiSolverInterface &si)
  {
      int nBasics = si.getNumRows();
      int nNonBasics = si.getNumCols();
      if (basis_ != NULL)
          delete basis_;
      basis_ = dynamic_cast<CoinWarmStartBasis *> (si.getWarmStart());
      if (!basis_)
          throw NoBasisError();
  
      if (nBasics_ > 0 || nBasics != nBasics_)
      {
          delete [] basics_;
          basics_ = NULL;
      }
      if (basics_ == NULL)
      {
          basics_ = new int[nBasics];
          nBasics_ = nBasics;
      }
  
      if (nNonBasics_ > 0 || nNonBasics != nNonBasics_)
      {
          delete [] nonBasics_;
          nonBasics_ = NULL;
      }
      if (nonBasics_ == NULL)
      {
          nonBasics_ = new int[nNonBasics];
          nNonBasics_ = nNonBasics;
      }
      int n = nBasics + nNonBasics;
      if ( nBasics_ + nNonBasics_ > 0 || nBasics_ + nNonBasics_ != n)
      {
          delete [] colsol_;
          delete [] integers_;
          integers_ = NULL;
          colsol_ = NULL;
          slacks_ = NULL;
      }
      if (colsol_ == NULL)
      {
          colsol_ = new double[n];
          slacks_ = &colsol_[nNonBasics];
      }
  
      if (integers_ == NULL)
      {
          integers_ = new bool[n];
      }
  
      const double * rowLower = si.getRowLower();
      const double * rowUpper = si.getRowUpper();
      //determine which slacks are integer
      const CoinPackedMatrix * m = si.getMatrixByCol();
      const double * elems = m->getElements();
      const int * inds = m->getIndices();
      const CoinBigIndex * starts = m->getVectorStarts();
      const int * lengths = m->getVectorLengths();
      //    int numElems = m->getNumElements();
      int numCols = m->getNumCols();
      assert(numCols == nNonBasics_);
      //   int numRows = m->getNumRows();
      CoinFillN(integers_ ,n, true);
      for (int i = 0 ;  i < numCols ; i++)
      {
          if (si.isContinuous(i))
              integers_[i] = false;
      }
      bool * integerSlacks = integers_ + numCols;
      for (int i = 0 ; i < nBasics ; i++)
      {
          if (rowLower[i] > -1e50 && INT_INFEAS(rowLower[i]) > 1e-15)
              integerSlacks[i] = false;
          if (rowUpper[i] < 1e50 && INT_INFEAS(rowUpper[i]) > 1e-15)
              integerSlacks[i] = false;
      }
      for (int i = 0 ;  i < numCols ; i++)
      {
          CoinBigIndex end = starts[i] + lengths[i];
          if (integers_[i])
          {
              for (CoinBigIndex k=starts[i] ; k < end; k++)
              {
                  if (integerSlacks[inds[k]] && INT_INFEAS(elems[k])>1e-15 )
                      integerSlacks[inds[k]] = false;
              }
          }
          else
          {
              for (CoinBigIndex k=starts[i] ; k < end; k++)
              {
                  if (integerSlacks[inds[k]])
                      integerSlacks[inds[k]] = false;
              }
          }
      }
  
      CoinCopyN(si.getColSolution(), si.getNumCols(), colsol_);
      CoinCopyN(si.getRowActivity(), si.getNumRows(), slacks_);
      for (int i = 0 ; i < si.getNumRows() ; i++)
      {
          slacks_[i]*=-1;
          if (rowLower[i]>-1e50)
          {
              slacks_[i] += rowLower[i];
          }
          else
          {
              slacks_[i] += rowUpper[i];
          }
      }
      //Now get the fill the arrays;
      nNonBasics = 0;
      nBasics = 0;
  
  
  
      //For having the index variables correctly ordered we need to access to OsiSimplexInterface
      {
          OsiSolverInterface * ncSi = (const_cast<OsiSolverInterface *>(&si));
          ncSi->enableSimplexInterface(0);
          ncSi->getBasics(basics_);
  	// Save enabled solver
  	solver_ = si.clone();
  #ifdef COIN_HAS_OSICLP
  	OsiClpSolverInterface * clpSi = dynamic_cast<OsiClpSolverInterface *>(solver_);
  	const OsiClpSolverInterface * clpSiRhs = dynamic_cast<const OsiClpSolverInterface *>(&si);
  	if (clpSi)
  	  clpSi->getModelPtr()->copyEnabledStuff(clpSiRhs->getModelPtr());;
  #endif
          ncSi->disableSimplexInterface();
      }
  
      int numStructural = basis_->getNumStructural();
      for (int i = 0 ; i < numStructural ; i++)
      {
          if (basis_->getStructStatus(i)== CoinWarmStartBasis::basic)
          {
              nBasics++;
              //Basically do nothing
          }
          else
          {
              nonBasics_[nNonBasics++] = i;
          }
      }
  
      int numArtificial = basis_->getNumArtificial();
      for (int i = 0 ; i < numArtificial ; i++)
      {
          if (basis_->getArtifStatus(i)== CoinWarmStartBasis::basic)
          {
              //Just check number of basics
              nBasics++;
          }
          else
          {
              nonBasics_[nNonBasics++] = i + basis_->getNumStructural();
          }
      }
  }
  void
  CglLandP::CachedData::clean(){
      if (basics_!=NULL)
          delete [] basics_;
      basics_ = NULL;
      if (nonBasics_!=NULL)
          delete [] nonBasics_;
      nonBasics_ = NULL;
      if (colsol_ != NULL)
          delete [] colsol_;
      colsol_ = NULL;
      delete basis_;
      basis_ = NULL;
      if (integers_)
          delete [] integers_;
      integers_ = NULL;
  
     nBasics_ = 0;
     nNonBasics_ = 0;
     delete solver_;
     solver_ = NULL;
  }
  CglLandP::CachedData::~CachedData()
  {
      if (basics_!=NULL)
          delete [] basics_;
      if (nonBasics_!=NULL)
          delete [] nonBasics_;
      if (colsol_ != NULL)
          delete [] colsol_;
      delete basis_;
      if (integers_)
          delete [] integers_;
      delete solver_;
  }
  
  CglLandP::CglLandP(const CglLandP::Parameters &params,
                     const LAP::Validator &validator):
          params_(params), cached_(), validator_(validator), numcols_(-1),
          originalColLower_(NULL), originalColUpper_(NULL),
          canLift_(false),
          extraCuts_()
  {
      handler_ = new CoinMessageHandler();
      handler_->setLogLevel(0);
      messages_ = LapMessages();
  }
  
  
  CglLandP::~CglLandP()
  {
      delete handler_;
      if (originalColLower_ != NULL)
          delete [] originalColLower_;
      if (originalColUpper_ != NULL)
          delete [] originalColUpper_;
  }
  
  CglLandP::CglLandP(const CglLandP & source):
          CglCutGenerator(source),
          params_(source.params_), cached_(source.cached_),
          validator_(source.validator_), numcols_(source.numcols_),
          originalColLower_(NULL), originalColUpper_(NULL),
          canLift_(source.canLift_),
          extraCuts_(source.extraCuts_)
  {
      handler_ = new CoinMessageHandler();
      handler_->setLogLevel(source.handler_->logLevel());
      messages_ = LapMessages();
      if (numcols_ != -1)
      {
          assert(numcols_ > 0);
          assert(originalColLower_!=NULL);
          assert(originalColUpper_!=NULL);
          originalColLower_ = new double[numcols_];
          originalColUpper_ = new double[numcols_];
          CoinCopyN(source.originalColLower_,numcols_,originalColLower_);
          CoinCopyN(source.originalColUpper_,numcols_,originalColUpper_);
      }
  }
  
  /** Assignment operator */
  CglLandP& CglLandP::operator=(const CglLandP &rhs)
  {
      if (this != &rhs)
      {
          params_ = rhs.params_;
          cached_ = rhs.cached_;
          validator_ = rhs.validator_;
          extraCuts_ = rhs.extraCuts_;
      }
      return *this;
  }
  
  
  CglCutGenerator *
  CglLandP::clone() const
  {
      return new CglLandP(*this);
  }
  
  extern double restaurationTime;
  
  struct cutsCos
  {
      int i;
      int j;
      double angle;
      cutsCos(int  i_, int j_ , double angle_):i(i_), j(j_), angle(angle_)
      {
      }
      bool operator<(const cutsCos&other)const
      {
          return angle > other.angle;
      }
  };
  
  
  void
  CglLandP::scanExtraCuts(OsiCuts& cs, const double * colsol) const
  {
      int numAdded = 0;
      for (int i = extraCuts_.sizeRowCuts() - 1; i > -1 ; i--)
      {
          double violation = extraCuts_.rowCut(i).violated(colsol);
          if (violation > 0.)
          {
              cs.insert(extraCuts_.rowCut(i));
              numAdded++;
              //      std::cout<<"A cut computed in a previous iteration is violated by "<<violation<<"."<<std::endl;
              //extraCuts_.eraseRowCut(i);
          }
      }
      //  std::cout<<"Added "<<numAdded<<" previously generated cuts."<<std::endl;
  }
  
  void
  CglLandP::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
                         const CglTreeInfo info )
  {
      if ((info.pass == 0) && !info.inTree)
      {
          numrows_ = si.getNumRows();
      }
  // scanExtraCuts(cs, si.getColSolution());
      Parameters params = params_;
      params.rhsWeight = numrows_ + 2;
  
      handler_->message(CUT_GAP, messages_)<<info.pass<<si.getObjValue() <<CoinMessageEol;
  
      if (info.inTree)   //put lower pivot limit
      {
          params.pivotLimit = std::min(params.pivotLimit, params.pivotLimitInTree);
          params.countMistakenRc = true;
      }
      if (params.timeLimit < 0)
      {
          params.pivotLimit = 0;
      }
  
      assert(si.basisIsAvailable());
  
  
  #ifdef APPEND_ROW
      OsiSolverInterface * t_si = si.clone();
      if (params.modularize)
      {
          int new_idx = si.getNumCols();
          int v_idx[1] = {new_idx};
          double v_val[1] = {-1};
          CoinPackedVector v(1, v_idx, v_val, false);
          t_si->addCol(CoinPackedVector(), 0, 1, 0);
          t_si->setInteger(new_idx);
          t_si->addRow(v,0, 0);
          t_si->resolve();
      }
  #else
      const OsiSolverInterface * t_si = &si;
  #endif
  
      cached_.getData(*t_si);
      CglLandPSimplex landpSi(*t_si, cached_, params, validator_);
      if (params.generateExtraCuts == CglLandP::AllViolatedMigs)
      {
          landpSi.genThisBasisMigs(cached_, params);
      }
      landpSi.setLogLevel(handler_->logLevel());
      int nCut = 0;
  
      std::vector<int> indices;
      getSortedFractionalIndices(indices,cached_, params);
  
  #ifndef NDEBUG
      int numrows = si.getNumRows();
  #endif
  
  #ifdef DO_STAT
      //Get informations on current optimum
      {
          OsiSolverInterface * gapTester = si.clone();
          gapTester->resolve();
  
          roundsStats_.analyseOptimalBasis(gapTester,info.pass, numrows_);
          delete gapTester;
      }
  #endif
  
      params_.timeLimit += CoinCpuTime();
      CoinRelFltEq eq(1e-04);
  
      for (unsigned int i = 0; i < indices.size() && nCut < params.maxCutPerRound &&
              nCut < cached_.nBasics_ ; i++)
      {
  
          //Check for time limit
          int iRow = indices[i];
          assert(iRow < numrows);
          OsiRowCut cut;
          int code=1;
          OsiSolverInterface * ncSi = NULL;
  
          if (params.pivotLimit != 0)
          {
              ncSi = t_si->clone();
              landpSi.setSi(ncSi);
              ncSi->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
              ncSi->messageHandler()->setLogLevel(0);
          }
  
          int generated = 0;
          if (params.pivotLimit == 0)
          {
              generated = landpSi.generateMig(iRow, cut, params);
          }
          else
          {
              generated = landpSi.optimize(iRow, cut, cached_, params);
              if (params.generateExtraCuts == CglLandP::AllViolatedMigs)
              {
                  landpSi.genThisBasisMigs(cached_, params);
              }
              landpSi.resetSolver(cached_.basis_);
          }
          code = 0;
          if (generated)
              code = validator_(cut, cached_.colsol_, si, params, originalColLower_, originalColUpper_);
          if (!generated || code)
          {
              if (params.pivotLimit !=0)
              {
                  handler_->message(LAP_CUT_FAILED_DO_MIG, messages_)<<validator_.failureString(code)<<CoinMessageEol;
                  landpSi.freeSi();
                  OsiSolverInterface * ncSi = t_si->clone();
                  landpSi.setSi(ncSi);
                  params.pivotLimit = 0;
                  if (landpSi.optimize(iRow, cut, cached_, params))
                  {
                      code = validator_(cut, cached_.colsol_, si, params, originalColLower_, originalColUpper_);
                  }
                  params.pivotLimit = params_.pivotLimit;
              }
          }
  
          if (params.pivotLimit != 0)
          {
              landpSi.freeSi();
          }
          if (code)
          {
              handler_->message(CUT_REJECTED, messages_)<<
              validator_.failureString(code)<<CoinMessageEol;
          }
          else
          {
              if (canLift_)
              {
                  cut.setGloballyValid(true);
              }
              cs.insertIfNotDuplicate(cut, eq);
              //cs.insert(cut);
              {
                  //std::cout<<"Violation "<<cut.violated(cached_.colsol_)<<std::endl;
                  nCut++;
              }
          }
      }
  
      Cuts& extra = landpSi.extraCuts();
      for (int i = 0 ; i < cached_.nNonBasics_; i++)
      {
          OsiRowCut * cut = extra.rowCut(i);
          if (cut == NULL) continue;
          int code = validator_(*cut, cached_.colsol_, si, params,
                                originalColLower_, originalColUpper_);
          if (code)
          {
              handler_->message(LAP_CUT_FAILED_DO_MIG, messages_)
              <<validator_.failureString(code)<<CoinMessageEol;
          }
          else
          {
              cs.insertIfNotDuplicate(*cut, eq);
              {
                  nCut++;
              }
          }
          delete cut;
      }
  
      landpSi.outPivInfo(nCut);
      params_.timeLimit -= CoinCpuTime();
  
      cached_.clean();
  #ifdef APPEND_ROW
      assert(t_si != &si);
      delete t_si;
  #endif
  }
  
  
  template < class S, class T, class U >
  class StableCompare
  {
  public:
      inline bool operator()(const CoinTriple<S,T,U>& t1,
                             const CoinTriple<S,T,U>& t2) const
      {
          return (t1.third < t2.third) ||
                 ((t1.third == t2.third) && (t1.second < t2.second));
      }
  
  };
  
  template <class T1,class T2>
  struct StableExternalComp
  {
      const std::vector<T1> &vec_1_;
      const std::vector<T2> &vec_2_;
      StableExternalComp(const std::vector<T1> &vec_1,
                         const std::vector<T2> &vec_2):
              vec_1_(vec_1),
              vec_2_(vec_2)
      {
      }
      CoinRelFltEq eq;
      bool operator()(int i, int j)
      {
          bool result = (vec_1_[i] < vec_1_[j]) ||
                        ( ((vec_1_[i]== vec_1_[j]))
                          && (vec_2_[i] < vec_2_[j]));
          return result;
      }
  
  };
  void
  CglLandP::getSortedFractionalIndices(std::vector<int> &frac_indices,
                                       const CachedData &data,
                                       const CglLandP::Parameters & params) const
  {
      std::vector<int> colIndices;
      std::vector<double> values;
      std::vector<int> indices;
      for (int i = 0 ; i < data.nBasics_ ; i++)
      {
          const int& iCol = data.basics_[i];
          if (iCol >= data.nNonBasics_ ||
                  !data.integers_[iCol] ||
                  INT_INFEAS(data.colsol_[iCol]) <= params.away)
              continue;
          const double value = INT_INFEAS(data.colsol_[iCol]);
  
          frac_indices.push_back(i);
          indices.push_back(static_cast<int>(values.size()));
          values.push_back(- value);
          colIndices.push_back(iCol);
      }
      std::sort(indices.begin(), indices.end(),StableExternalComp<double, int>(values,colIndices));
      colIndices = frac_indices;
      for (unsigned int i = 0; i < indices.size() ; i++)
      {
          frac_indices[i] = colIndices[indices[i]];
      }
  
  }
  
  
  