// Copyright (C) 2003, International Business Machines
// Corporation and others.  All Rights Reserved.

#include "CoinPragma.hpp"
#include "CoinIndexedVector.hpp"
#include "ClpFactorization.hpp"
#include "ClpSimplex.hpp"
#include "ClpLinearObjective.hpp"
#include "CoinHelperFunctions.hpp"

//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor 
//-------------------------------------------------------------------
ClpLinearObjective::ClpLinearObjective () 
: ClpObjective()
{
  type_=1;
  objective_=NULL;
  numberColumns_=0;
}

//-------------------------------------------------------------------
// Useful Constructor 
//-------------------------------------------------------------------
ClpLinearObjective::ClpLinearObjective (const double * objective , 
					int numberColumns) 
  : ClpObjective()
{
  type_=1;
  numberColumns_=numberColumns;
  objective_ = CoinCopyOfArray(objective,numberColumns_,0.0);
}

//-------------------------------------------------------------------
// Copy constructor 
//-------------------------------------------------------------------
ClpLinearObjective::ClpLinearObjective (const ClpLinearObjective & rhs) 
: ClpObjective(rhs)
{  
  numberColumns_=rhs.numberColumns_;
  objective_ = CoinCopyOfArray(rhs.objective_,numberColumns_);
}
/* Subset constructor.  Duplicates are allowed
   and order is as given.
*/
ClpLinearObjective::ClpLinearObjective (const ClpLinearObjective &rhs,
					int numberColumns, 
					const int * whichColumn) 
: ClpObjective(rhs)
{
  objective_=NULL;
  numberColumns_=0;
  if (numberColumns>0) {
    // check valid lists
    int numberBad=0;
    int i;
    for (i=0;i<numberColumns;i++)
      if (whichColumn[i]<0||whichColumn[i]>=rhs.numberColumns_)
	numberBad++;
    if (numberBad)
      throw CoinError("bad column list", "subset constructor", 
		      "ClpLinearObjective");
    numberColumns_ = numberColumns;
    objective_ = new double[numberColumns_];
    for (i=0;i<numberColumns_;i++) 
      objective_[i]=rhs.objective_[whichColumn[i]];
  }
}
  

//-------------------------------------------------------------------
// Destructor 
//-------------------------------------------------------------------
ClpLinearObjective::~ClpLinearObjective ()
{
  delete [] objective_;
}

//----------------------------------------------------------------
// Assignment operator 
//-------------------------------------------------------------------
ClpLinearObjective &
ClpLinearObjective::operator=(const ClpLinearObjective& rhs)
{
  if (this != &rhs) {
    ClpObjective::operator=(rhs);
    numberColumns_=rhs.numberColumns_;
    delete [] objective_;
    objective_ = CoinCopyOfArray(rhs.objective_,numberColumns_);
  }
  return *this;
}

// Returns gradient
double *  
ClpLinearObjective::gradient(const ClpSimplex * model,
			     const double * solution, double & offset,bool refresh,
			     int includeLinear)
{
  // not sure what to do about scaling
  //assert (!model);
  assert (includeLinear==2); //otherwise need to return all zeros
  offset=0.0;
  return objective_;
}
  
/* Returns reduced gradient.Returns an offset (to be added to current one).
 */
double 
ClpLinearObjective::reducedGradient(ClpSimplex * model, double * region,
				 bool useFeasibleCosts)
{
  int numberRows = model->numberRows();
  //work space
  CoinIndexedVector  * workSpace = model->rowArray(0);
  
  CoinIndexedVector arrayVector;
  arrayVector.reserve(numberRows+1);
  
  int iRow;
#ifdef CLP_DEBUG
  workSpace->checkClear();
#endif
  double * array = arrayVector.denseVector();
  int * index = arrayVector.getIndices();
  int number=0;
  const double * cost = model->costRegion();
  assert (!useFeasibleCosts);
  const int * pivotVariable = model->pivotVariable();
  for (iRow=0;iRow<numberRows;iRow++) {
    int iPivot=pivotVariable[iRow];
    double value = cost[iPivot];
    if (value) {
      array[iRow]=value;
      index[number++]=iRow;
    }
  }
  arrayVector.setNumElements(number);

  int numberColumns=model->numberColumns();
  
  // Btran basic costs
  double * work = workSpace->denseVector();
  model->factorization()->updateColumnTranspose(workSpace,&arrayVector);
  ClpFillN(work,numberRows,0.0);
  // now look at dual solution
  double * rowReducedCost = region+numberColumns;
  double * dual = rowReducedCost;
  double * rowCost = model->costRegion(0);
  for (iRow=0;iRow<numberRows;iRow++) {
    dual[iRow]=array[iRow];
  }
  double * dj = region;
  ClpDisjointCopyN(model->costRegion(1),numberColumns,dj);
  model->transposeTimes(-1.0,dual,dj);
  for (iRow=0;iRow<numberRows;iRow++) {
    // slack
    double value = dual[iRow];
    value += rowCost[iRow];
    rowReducedCost[iRow]=value;
  }
  return 0.0;
}
/* Returns step length which gives minimum of objective for
   solution + theta * change vector up to maximum theta.
   
   arrays are numberColumns+numberRows
*/
double 
ClpLinearObjective::stepLength(ClpSimplex * model,
			       const double * solution,
			       const double * change,
			       double maximumTheta,
			       double & currentObj,
			       double & predictedObj,
			       double & thetaObj)
{
  const double * cost = model->costRegion();
  double delta=0.0;
  int numberRows = model->numberRows();
  int numberColumns = model->numberColumns();
  currentObj=0.0;
  thetaObj=0.0;
  for (int iColumn=0;iColumn<numberColumns+numberRows;iColumn++) {
    delta += cost[iColumn]*change[iColumn];
    currentObj += cost[iColumn]*solution[iColumn];
  }
  thetaObj =currentObj + delta*maximumTheta;
  predictedObj =currentObj + delta*maximumTheta;
  if (delta<0.0) {
    return maximumTheta;
  } else {
    printf("odd linear direction %g\n",delta);
    return 0.0;
  }
}
// Return objective value (without any ClpModel offset) (model may be NULL)
double 
ClpLinearObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
{
  const double * cost = objective_;
  if (model&&model->costRegion())
    cost = model->costRegion();
  double currentObj=0.0;
  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
    currentObj += cost[iColumn]*solution[iColumn];
  }
  return currentObj;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpObjective * ClpLinearObjective::clone() const
{
  return new ClpLinearObjective(*this);
}
/* Subset clone.  Duplicates are allowed
   and order is as given.
*/
ClpObjective * 
ClpLinearObjective::subsetClone (int numberColumns, 
			   const int * whichColumns) const
{
  return new ClpLinearObjective(*this, numberColumns, whichColumns);
}
// Resize objective
void 
ClpLinearObjective::resize(int newNumberColumns)
{
  if (numberColumns_!=newNumberColumns) {
    int i;
    double * newArray = new double[newNumberColumns];
    if (objective_)
      CoinMemcpyN(objective_,CoinMin(newNumberColumns,numberColumns_),newArray);
    delete [] objective_;
    objective_ = newArray;
    for (i=numberColumns_;i<newNumberColumns;i++) 
      objective_[i]=0.0;
    numberColumns_ = newNumberColumns;
  } 
  
}
// Delete columns in  objective
void 
ClpLinearObjective::deleteSome(int numberToDelete, const int * which) 
{
  if (objective_) {
    int i ;
    char * deleted = new char[numberColumns_];
    int numberDeleted=0;
    CoinZeroN(deleted,numberColumns_);
    for (i=0;i<numberToDelete;i++) {
      int j = which[i];
      if (j>=0&&j<numberColumns_&&!deleted[j]) {
	numberDeleted++;
	deleted[j]=1;
      }
    }
    int newNumberColumns = numberColumns_-numberDeleted;
    double * newArray = new double[newNumberColumns];
    int put=0;
    for (i=0;i<numberColumns_;i++) {
      if (!deleted[i]) {
	newArray[put++]=objective_[i];
      }
    }
    delete [] objective_;
    objective_ = newArray;
    delete [] deleted;
    numberColumns_ = newNumberColumns;
  }
}
// Scale objective 
void 
ClpLinearObjective::reallyScale(const double * columnScale) 
{
  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
    objective_[iColumn] *= columnScale[iColumn];
  }
}