//===========================================================================//
// This file is part of the DIP Solver Framework.                            //
//                                                                           //
// DIP is distributed under the Eclipse Public License as part of the        //
// COIN-OR repository (http://www.coin-or.org).                              //
//                                                                           //
// Authors: Matthew Galati, SAS Institute Inc. (matthew.galati@sas.com)      //
//          Ted Ralphs, Lehigh University (ted@lehigh.edu)                   //
//          Jiadong Wang, Lehigh University (jiw408@lehigh.edu)              //
//                                                                           //
// Copyright (C) 2002-2015, Lehigh University, Matthew Galati, Ted Ralphs    //
// All Rights Reserved.                                                      //
//===========================================================================//

//===========================================================================//
#include "DecompApp.h"
#include "DecompVar.h"
#include "DecompAlgoRC.h"

using namespace std;

//===========================================================================//
DecompPhase DecompAlgoRC::phaseInit()
{
   //THINK: should base have some container for master and sub
   //since in RC, don't need OSI??
   //so base should never call m_masterSI, but m_masterContainer
   //or something
   if (m_param.LogDumpModel > 1)
      printCurrentProblem(m_masterSI,
                          "masterProb",
                          m_nodeStats.nodeIndex,
                          m_nodeStats.cutCallsTotal,
                          m_nodeStats.priceCallsTotal);

   //---
   //--- update primal/dual vectors
   //---
   m_status = STAT_FEASIBLE;
   //TODO: what about the INF case?? artificial columns? DC-ABCC version
   // ---
   // --- update the phase
   // ---
   return PHASE_PRICE2;
}

//===========================================================================//
void DecompAlgoRC::phaseDone()
{
   //take the current set of variables and solve DW master to get primal
   //TODO: right now, creating from scratch each time -- really need
   //      to append and warm start - esp if doing alot of branching
   //delete m_masterSI;
   //m_masterSI = NULL;
   //m_masterSI = new OsiLpSolverInterface();
   //CoinAssertHint(m_masterSI, "Error: Out of Memory");
   //m_masterSI->messageHandler()->setLogLevel(m_param.LogLpLevel);
   DecompConstraintSet*           modelCore   = m_modelCore.getModel();
#if 0
   //---
   //--- Initialize the solver interface for the master problem.
   //--- PC: min c(s lam)
   //---     A''s   lam   >= b'',
   //---     sum{s} lam_s  = 1  ,
   //---            lam_s >= 0  , s in F'[0]
   //--- modelCore contains [A'', b''], from the original model in
   //--- terms of x. In this function we create the DW-LP in terms of
   //--- lambda, [[A''s, 1], [b'', 1]] and load that into the OSI
   //--- interface m_masterSI.
   //---
   //--- NOTE: if 0 is feasible to subproblem, we can relax convexity to <= 1
   //---
   UtilPrintFuncBegin(m_osLog, m_classTag,
                      "phaseDone()", m_param.LogDebugLevel, 2);
   CoinAssert(m_vars.size() > 0);
   int nColsCore = modelCore->getNumCols();
   int nRowsCore = modelCore->getNumRows();
   modelCore->nBaseRowsOrig = modelCore->nBaseRows;
   modelCore->nBaseRows     = nRowsCore;
   //THINK? should initVars be in pool and do addVarsFromPool here?
   CoinPackedMatrix* M = new CoinPackedMatrix(true, 0, 0);
   M->setDimensions(nRowsCore + 1, 0);
   const int n_cols  = static_cast<int>(m_vars.size());
   double* colLB    = new double[n_cols];
   double* colUB    = new double[n_cols];
   double* obj      = new double[n_cols];
   double* denseCol = new double[nRowsCore + 1];
   CoinAssertHint(colLB && colUB && obj && denseCol, "Error: Out of Memory");
   int col_index     = 0;
   DecompVarList::iterator li;

   for (li = m_vars.begin(); li != m_vars.end(); li++) {
      UTIL_DEBUG(m_param.LogDebugLevel, 5,
                 (*li)->print(m_osLog, m_app);
                );
      //---
      //--- get dense column = A''s, append convexity constraint on end
      //---
      modelCore->M->times((*li)->m_s, denseCol);
      denseCol[nRowsCore] = 1.0;
      //---
      //--- create a sparse column from the dense column
      //---
      // THINK: do i need a DecompCol?
      // THINK: does this allocate memory for coinpackedvec twice?
      CoinPackedVector* sparseCol
      = UtilPackedVectorFromDense(nRowsCore + 1,
                                  denseCol, m_param.TolZero);
      UTIL_DEBUG(m_param.LogDebugLevel, 5,
                 (*m_osLog) << "\nSparse Col: \n";
                 UtilPrintPackedVector(*sparseCol, m_osLog);
                );
      //TODO: do in const blocks
      // ---
      // --- append the sparse column to the matrix
      // ---
      M->appendCol(*sparseCol);
      colLB[col_index] = 0.0; //THINK: (*li)->getLowerBound();
      colUB[col_index] = DecompInf; //THINK: (*li)->getUpperBound(); //FIX!!
      obj[col_index]   = (*li)->getOriginalCost(); //c.s
      col_index++;
      UTIL_DELPTR(sparseCol); //THINK
   }

   //---
   //--- THINK: do we want to adjust modelCore directly here?
   //--- adjust row bounds for convexity constraint
   //---
   //TODO: in memory
   double* zeroSol    = new double[nColsCore];
   CoinAssertHint(zeroSol, "Error: Out of Memory");
   UtilFillN(zeroSol, nColsCore, 0.0);
   //TODO - REVISIT - that's not the right check
   //  needs to be feasible to subproblem?
   bool     isZeroFeas = isIPFeasible(zeroSol);
   UTIL_DEBUG(m_param.LogDebugLevel, 5,

              if (isZeroFeas)
              (*m_osLog) << "Zero Sol is Feasible - relax convexity con.\n";
             );

   vector<double> masterRowLB(modelCore->rowLB);
   vector<double> masterRowUB(modelCore->rowUB);

   if (isZeroFeas) {
      masterRowLB.push_back(-DecompInf);
      masterRowUB.push_back(1.0);
   } else {
      masterRowLB.push_back(1.0);
      masterRowUB.push_back(1.0);
   }

   //---
   //--- load the problem into master's solver interface
   //---
   assert(M->getNumRows() == static_cast<int>(masterRowLB.size()));
   assert(M->getNumRows() == static_cast<int>(masterRowUB.size()));
   m_masterSI->loadProblem(*M,  colLB, colUB, obj,
                           &masterRowLB[0],
                           &masterRowUB[0]);
   // ---
   // --- free local memory
   // ---
   UTIL_DELPTR(M);
   UTIL_DELARR(denseCol);
   UTIL_DELARR(colLB);
   UTIL_DELARR(colUB);
   UTIL_DELARR(obj);
   UTIL_DELARR(zeroSol);
#endif
   m_status = DecompAlgo::solutionUpdate(PHASE_UNKNOWN, 99999, 99999);

   //---
   //--- check if IP feasible (are we done?)
   //--- TODO: for nonexplicity, also check user app isfeasible
   //---
   //TODO: should this whole section be phaseDone?
   if (m_status != STAT_INFEASIBLE) {
      DecompAlgo::recomposeSolution(m_masterSI->getColSolution(), m_xhat);
      UTIL_DEBUG(m_param.LogDebugLevel, 4,
                 m_app->printOriginalSolution(modelCore->getNumCols(),
                                              modelCore->getColNames(),
                                              m_xhat);
                );

      if (isIPFeasible(m_xhat)) {
         if (m_app->APPisUserFeasible(m_xhat,
                                      modelCore->getNumCols(),
                                      m_param.TolZero)) {
            DecompSolution* decompSol
            = new DecompSolution(modelCore->getNumCols(),
                                 m_xhat, m_masterSI->getObjValue());
            m_xhatIPFeas.push_back(decompSol);
         }
      }

      vector<DecompSolution*>::iterator vi;
      DecompSolution* viBest = NULL;
      double bestBoundUB = m_nodeStats.objBest.second;

      for (vi = m_xhatIPFeas.begin(); vi != m_xhatIPFeas.end(); vi++) {
         const DecompSolution* xhatIPFeas = *vi;

         if (isIPFeasible(xhatIPFeas->getValues())) {
            if (xhatIPFeas->getQuality() <= bestBoundUB) {
               bestBoundUB = xhatIPFeas->getQuality();
               viBest = *vi;
            }
         }
      }

      if (viBest) {
         //save the best
         setObjBoundIP(bestBoundUB);
         m_xhatIPBest = viBest;
      }
   }

   UtilPrintFuncEnd(m_osLog, m_classTag,
                    "phaseDone()", m_param.LogDebugLevel, 2);
}

//===========================================================================//
void DecompAlgoRC::recomposeSolution(const double* solution,
                                     double*        rsolution)
{
   printf("RC recomposeSolution does nothing\n");
}

//===========================================================================//
void DecompAlgoRC::createMasterProblem(DecompVarList& initVars)
{
   //---
   //--- there is no master LP in RC, just initialize the dual vector
   //---
   UtilPrintFuncBegin(m_osLog, m_classTag,
                      "createMasterProblem()", m_param.LogDebugLevel, 2);
   DecompAlgo::createMasterProblem(initVars);
   CoinAssert(initVars.size() > 0);
   //---
   //--- In order to implement simple branching, we are going to
   //--- treat all column bounds as explicit constraints. Then branching
   //--- for DW can be done in the same way it is done for regular CPM.
   //---
   //TODO: looks like volume doesn't like R rows... change and split up?
   //printf("Volume Algorithm can't work if there is a non ELG row\n");
   //   coreMatrixAppendColBounds();
   //THINK: is this the right place for this
   //TODO: give user option to feed in a good starting dual vector?
   //you might want to do that if switch between DW and RC... etc... THINK
   DecompConstraintSet*           modelCore   = m_modelCore.getModel();
   fill_n(back_inserter(m_u), modelCore->getNumRows(), 0.0);
   //TODO
   m_rc = new double[modelCore->getNumCols()]; //better in constructor?
   CoinAssertHint(m_rc, "Error: Out of Memory");
   //m_vars here will contain all the shat's used - we can only
   //use one of them in subgradient, so just save the last one
   //DecompVarList::iterator it = initVars.begin();
   //assert(*it);
   //think
   //double redCost = (*it)->getOriginalCost() - (*it)->m_s.dotProduct(&m_u[0]);
   //(*it)->setReducedCost(redCost);
   //for rc we want to calc the reduced cost
   //m_vars.push_back(*it);
   //UTIL_DEBUG(m_param.LogDebugLevel, 3,
   //	      (*it)->print(m_osLog);
   //	      );
   UtilPrintFuncEnd(m_osLog, m_classTag,
                    "createMasterProblem()", m_param.LogDebugLevel, 2);
}

// ------------------------------------------------------------------------- //
bool DecompAlgoRC::isDone()
{
   //iter count is checked by phaseUpdate
   //need to check step limit in here
   //need to check if ub-lb gap is small? isn't that always checked?
   //printf("\nm_UB: %12.10f, m_LB: %12.10f", m_UB, m_LB);
   if ((m_step < 1.0e-3) ||             //step length too small
         m_zeroSub         ||              //0 subgradient
         UtilIsZero(m_UB - m_LB, 1.0e-3)) { //gap is small
      return true;
   }

   return false;
}

//do we need a var pool at all here?

// ------------------------------------------------------------------------- //
int DecompAlgoRC::addCutsFromPool()
{
   int nNewRows = DecompAlgo::addCutsFromPool();
   m_u.reserve(m_u.size() + nNewRows);
   UtilFillN(m_u, nNewRows, 0.0);

   //is this the right place to do this?
   //we were pricing out, then cutting, but when we price out,
   //step=0, so we need to start over with step size
   //best place for this would be in phaseUpdatE?
   if (nNewRows > 0) {
      m_step = 2.0;
   }

   return nNewRows;
}

// ------------------------------------------------------------------------- //
int DecompAlgoRC::generateVars(DecompVarList&     newVars,
                               double&            mostNegReducedCost)
{
   UtilPrintFuncBegin(m_osLog, m_classTag,
                      "generateVars()", m_param.LogDebugLevel, 2);
   DecompConstraintSet*           modelCore   = m_modelCore.getModel();
   //really only returning one var here...
   //for RC, doesn't have to be negative??
   mostNegReducedCost = DecompInf;//bad name here
   //TODO: whenever a cut is added, if doing RC, you need to add an
   //element to u.... do we need to override gen cuts just for that?
   //gen cuts (x) is wrong here anyway... need gen cuts (s)
   assert(static_cast<int>(m_u.size()) == modelCore->getNumRows());
   //THINK:
   // if we overload getDual method, then this can be same as PC genvars?
   // in PC, we are letting OSI house the duals, in RC we do it ourself
   // seems silly to use OSI at all for RC?
   //
   // almost want Decomp to house dual even for PC...
   //TODO: m_app->m_model seems dumb
   //reduced cost = c - uA
   const double* origObjective = getOrigObjective();
   modelCore->M->transposeTimes(&m_u[0], m_rc);

   for (int c = 0; c < modelCore->getNumCols(); c++) {
      printf("RC[%d] -> c: %g - uA: %g = m_rc: %g\n",
             c, origObjective[c], m_rc[c],
             origObjective[c] - m_rc[c]);
      m_rc[c] = origObjective[c] - m_rc[c];
   }

   //double alpha = 0.0;
   DecompVarList potentialVars;
   //TODO: stat return, restrict how many? pass that in to user?
   //only take those with negative reduced cost?
   //check for dups here
   //TODO: blocks!
   //WRONG...
#if 0
   solveRelaxed(
      0,
      m_rc, origObjective, alpha,
      modelCore->getNumCols(), false, true, true, false,
      m_subprobSI[0],
      NULL,//WILL FAIL
      potentialVars);//NO CHECK RC??
#endif
   //another way to do this is to just collect all m_vars,
   //not worrying about duplicates -- and when we get to DW
   //strip out the dups before constructing the master formulation
   DecompVarList::iterator it;
   double varRedCost;

   for (it = potentialVars.begin(); it != potentialVars.end(); it++) {
      varRedCost = (*it)->getReducedCost();
      newVars.push_back(*it);

      if (varRedCost < mostNegReducedCost) {
         mostNegReducedCost = varRedCost;
         //TODO: if this winds up being a dup in addVarsToPool
         //this memory will be erased -- which will cause a problem
         //because we need this var
         //TODO: FUGLY - make a copy so if dup doesn't cause prob
         //i don't like this! ugh
         //m_shatVar = *it; //THINK
         m_shatVar = *(*it);
         m_shatVar.fillDenseArr(modelCore->getNumCols(), m_xhat);
      }
   }

   potentialVars.clear(); //THINK? what does clear do exactly ?
#if 0

   for (it = newVars.begin(); it != newVars.end(); it++)
      UTIL_DEBUG(m_app->m_param.LogDebugLevel, 3,
                 (*it)->print(m_osLog);
                );

#endif
   UtilPrintFuncEnd(m_osLog, m_classTag,
                    "generateVars()", m_param.LogDebugLevel, 2);
   return static_cast<int>(newVars.size());
}

// ------------------------------------------------------------------------- //
DecompStatus DecompAlgoRC::solutionUpdate(const DecompPhase phase,
					  const int         maxInnerIter,
					  const int         maxOuterIter)
{
   //---
   //--- C, PC: This step solves (or takes a few steps to solve) master LP
   //---        which updates both the primal (x,lambda) and dual(u) vectors.
   //--- RC   : This does one step of subgradient, which updates the dual (u)
   //---        vector, given some shat (the last variable in m_var).
   //generalize this for the many variants of subgradient?
   //TODO: might make life easier to flip all inequalities to <= or >= !
   UtilPrintFuncBegin(m_osLog, m_classTag,
                      "solutionUpdate()", m_param.LogDebugLevel, 2);
   m_UB = 1.05 * m_nodeStats.objBest.second; //TODO heuristics 1.05??
   //TODO: tols
   //how to allow hooks to other stabilization methods, bundle, etc... user
   //can simply derive from RC or base and recode this method...
   UTIL_DEBUG(m_app->m_param.LogDebugLevel, 3,
              (*m_osLog) << "\nVARS m_vars:\n";
              printVars(m_osLog);
             );
   DecompVar* shatVar = &m_shatVar; //m_vars.back();
   //DecompVar * shatVar = m_vars.back();
   CoinAssert(shatVar);
   UTIL_DEBUG(m_app->m_param.LogDebugLevel, 3,
              (*m_osLog) << "\nshat: ";
              shatVar->print(m_osLog);
             );
   //make this part of class, else realloc every iter...
   //use vector, let it grow?
   int      r;
   //char     sense;
   //double   range;
   DecompConstraintSet*           modelCore   = m_modelCore.getModel();
   int      n_coreRows = modelCore->getNumRows();
   //const double * rhs        = getRightHandSide();
   //const char   * sense      = getRowSense();
   const double* rhs        = &modelCore->rowRhs[0];
   const char*    sense      = &modelCore->rowSense[0];
   double* violation  = new double[n_coreRows];
   double* activity   = new double[n_coreRows];
   modelCore->M->times(shatVar->m_s, activity); //As
   assert(static_cast<int>(m_u.size()) == n_coreRows);
   // =, b  - Ax or Ax - b
   //>=, b  - Ax > 0 is a violation
   //<=, b  - Ax < 0 is a violation
   m_zeroSub = true;

   for (r = 0; r < n_coreRows; r++) {
      violation[r] = rhs[r] - activity[r];
      UTIL_DEBUG(m_app->m_param.LogDebugLevel, 3,
                 (*m_osLog) << setprecision(8);
                 (*m_osLog) << "r: " << r << " vio: " << violation[r]
                 << " rhs: " << rhs[r] << " act: " << activity[r]
                 << " u: " << m_u[r] << " sense: " << sense[r];
                );

      //beasley suggestion... this is just another basic variant of SG?
      // =, b - Ax != 0 is a violation
      //>=, b - Ax > 0  is a violation
      //<=, Ax - b > 0  is a violation
      switch (sense[r]) {
      case 'E':
         //violation[i] = rhs[i] - activity[i];
         break;
      case 'G':

         //violation[i] = rhs[i] - activity[i];
         //TODO: use tol
         if (violation[r] < 0.0 && m_u[r] >= -1.0e-4 && m_u[r] <= 1.0e-4) {
            violation[r] = 0.0;
         }

         break;
      case 'L':

         //violation[i] = rhs[i] - activity[i];
         if (violation[r] > 0.0 && m_u[r] >= -1.0e-4 && m_u[r] <= 1.0e-4) {
            violation[r] = 0.0;
         }

         break;
      }

      UTIL_DEBUG(m_app->m_param.LogDebugLevel, 3,
                 (*m_osLog) << " -> vio: " << violation[r] << "\n";
                );

      //?? shouldn't it be if all are feasible?
      if (fabs(violation[r]) > 0.0001) {
         m_zeroSub = false;
      }
   }

   //when to half the step size?
   //same stuff as m_tlb?? setTrueLowerBound?
   double bound    = shatVar->getReducedCost();//c - uA (is this set?)
   //double constant = calcConstant(m_u.size(), &m_u[0]);
   double constant = 0.0;

   for (r = 0; r < n_coreRows; r++) {
      constant += m_u[r] * rhs[r];
   }

   //needs rhs from OSI - this is why it was better to fake it
   //and have the information in OSI... but then carrying around
   //modelCore and OSI for no good reason... but... this is messier
   //LR Bound = (c - uA)shat + ub, assumes u >= 0
   //first iter, LR Bound = cshat - is that a valid LB?
   //yes, actually is a LB for any u >= 0
   //this assumes u is optimal? or just dual feasible
   //is u = 0 dual feasible?
   //but only a valid bound if shat has the lowest reduced cost
   //for a given u... so, for first iter of smallip, should have given
   //(2,1) with LB = 2... RC don't do genInitVars?
   //TODO: think initial dual vector - solve an LP to get started?
   if (bound + constant > m_LB + m_app->m_param.TolZero) {
      m_LB  = bound + constant;
      m_cntSameLB  = 0;
      //count_sameLB = 0;
      //make param to do or not
      //reducedCostFixing(reducedCost);
   } else {
      m_cntSameLB++;
   }

   //if(count_sameLB >= m_app->m_param.RC_sameLBLimit){
   if (m_cntSameLB >= 10) {
      m_step /= 2.0;
      cout << "LB has not changed in " << m_cntSameLB
           << " iterations - halve step: " << m_step << endl;
      m_cntSameLB = 0;
   }

   printf("m_UB: %12.10f, m_LB: %12.10f\n", m_UB, m_LB);
   assert((m_UB - m_LB) > -0.0001);
   double theta = 0.0;
   double denom = 0.0;

   for (r = 0; r < n_coreRows; r++) {
      denom += violation[r] * violation[r];
   }

   if (denom > 0.0) {
      theta = m_step * ((1.00 * m_UB) - m_LB) / denom;
   }

   //TODO: debug util to print a list of values in a nice format to log?
   UTIL_DEBUG(m_app->m_param.LogDebugLevel, 3,
              (*m_osLog) << "m_UB: " << m_UB << " m_LB: " << m_LB
              << " denom: " << denom << " m_step: " << m_step
              << " theta: " << theta << "\n";
             );

   //STOP 10/6/07
   //How do we deal with range constraints? What does volume do, for example?
   for (r = 0; r < n_coreRows; r++) {
      switch (sense[r]) {
      case 'E':
         m_u[r] += theta * violation[r];
         break;
      case 'G':
         //u > 0, g_i > 0 for violations
         m_u[r] = max(0.0, m_u[r] + (theta * violation[r]));
         break;
      case 'L':
         //u < 0, g_i < 0 for violatoins
         m_u[r] = max(0.0, m_u[r] - (theta * violation[r]));
         break;
      default:
         assert(0);
      }

      UTIL_DEBUG(m_app->m_param.LogDebugLevel, 3,
                 (*m_osLog) << "r: " << r << " m_u: " << m_u[r] << "\n";
                );
   }

   m_iter++;
   //if no var pool buffer is used, make sure you clean up manually
   /*it++;
   UtilDeleteListPtr(potentialVars, it, potentialVars.end());
   potentialVars.clear();*/
   //TODO - back in algo, reduced cost fixing
   //TODO: temp memory don't alloc/free each time!
   //UTIL_DELARR(rhs);
   UTIL_DELARR(violation);
   UTIL_DELARR(activity);
   UtilPrintFuncEnd(m_osLog, m_classTag,
                    "solutionUpdate()", m_param.LogDebugLevel, 2);
   return STAT_FEASIBLE;
}

// ------------------------------------------------------------------------- //
bool DecompAlgoRC::updateObjBound(const double mostNegRC)
{
   //---
   //--- C    : LB = masterLP obj
   //--- PC   : LB = zDW_RMP + RC* <= zDW <= zDW_RMP
   //---    where RC* is the most negative reduced cost
   //---    assuming the relaxation subproblem was solved exactly
   //---
   UtilPrintFuncBegin(m_osLog, m_classTag,
                      "updateObjBound()", m_param.LogDebugLevel, 2);
   DecompConstraintSet*           modelCore   = m_modelCore.getModel();
   //mostNegRC not used?
   //DecompVar * shatVar = m_shatVar;//m_vars.back();
   //DecompVar * shatVar = m_vars.back();
   //CoinAssert(shatVar);
   int r;
   const int n_coreRows = modelCore->getNumRows();
   //double bound = shatVar->getReducedCost();//c - uA (is this set?)
   double constant = 0.0;
   const double* rhs = &modelCore->rowRhs[0];

   for (r = 0; r < n_coreRows; r++) {
      constant += m_u[r] * rhs[r];
   }

   //double thisBoundLB = shatVar->getReducedCost() + constant;
   double thisBoundLB = mostNegRC + constant;
   setObjBound(thisBoundLB, constant);
   UTIL_DEBUG(m_param.LogDebugLevel, 5,
              (*m_osLog)
              << "ThisLB = " << UtilDblToStr(thisBoundLB) << "\t"
              << "BestLB = " << UtilDblToStr(m_nodeStats.objBest.first)
              << "\n";
             );
   UtilPrintFuncEnd(m_osLog, m_classTag,
                    "updateObjBound()", m_param.LogDebugLevel, 2);
   return false;
}
