cplex_interface.cc

// Copyright 2014 IBM Corporation
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

// Initial version of this code was written by Daniel Junglas (IBM)

#include <limits>
#include <memory>

#include "absl/strings/str_format.h"
#include "ortools/base/integral_types.h"
#include "ortools/base/logging.h"
#include "ortools/base/timer.h"
#include "ortools/linear_solver/linear_solver.h"

#if defined(USE_CPLEX)

extern "C" {
#include "ilcplex/cplexx.h"
// This is an undocumented function, setting the objective offset
// is not supported everywhere (for example it may not be exported if a
// model is written to a file), but it works in the cases we need here.
CPXLIBAPI int CPXPUBLIC CPXEsetobjoffset(CPXCENVptr, CPXLPptr, double);
}

// In case we need to return a double but don't have a value for that
// we just return a NaN.
#if !defined(CPX_NAN)
#define CPX_NAN std::numeric_limits<double>::quiet_NaN()
#endif

// The argument to this macro is the invocation of a CPXX function that
// returns a status. If the function returns non-zero the macro aborts
// the program with an appropriate error message.
#define CHECK_STATUS(s)    \
  do {                     \
    int const status_ = s; \
    CHECK_EQ(0, status_);  \
  } while (0)

namespace operations_research {

using std::unique_ptr;

// For a model that is extracted to an instance of this class there is a
// 1:1 corresponence between MPVariable instances and CPLEX columns: the
// index of an extracted variable is the column index in the CPLEX model.
// Similiar for instances of MPConstraint: the index of the constraint in
// the model is the row index in the CPLEX model.
class CplexInterface : public MPSolverInterface {
 public:
  // NOTE: 'mip' specifies the type of the problem (either continuous or
  //       mixed integer. This type is fixed for the lifetime of the
  //       instance. There are no dynamic changes to the model type.
  explicit CplexInterface(MPSolver* const solver, bool mip);
  ~CplexInterface();

  // Sets the optimization direction (min/max).
  virtual void SetOptimizationDirection(bool maximize);

  // ----- Solve -----
  // Solve the problem using the parameter values specified.
  virtual MPSolver::ResultStatus Solve(MPSolverParameters const& param);

  // ----- Model modifications and extraction -----
  // Resets extracted model
  virtual void Reset();

  virtual void SetVariableBounds(int var_index, double lb, double ub);
  virtual void SetVariableInteger(int var_index, bool integer);
  virtual void SetConstraintBounds(int row_index, double lb, double ub);

  virtual void AddRowConstraint(MPConstraint* const ct);
  virtual void AddVariable(MPVariable* const var);
  virtual void SetCoefficient(MPConstraint* const constraint,
                              MPVariable const* const variable,
                              double new_value, double old_value);

  // Clear a constraint from all its terms.
  virtual void ClearConstraint(MPConstraint* const constraint);
  // Change a coefficient in the linear objective
  virtual void SetObjectiveCoefficient(MPVariable const* const variable,
                                       double coefficient);
  // Change the constant term in the linear objective.
  virtual void SetObjectiveOffset(double value);
  // Clear the objective from all its terms.
  virtual void ClearObjective();

  // ------ Query statistics on the solution and the solve ------
  // Number of simplex iterations
  virtual int64_t iterations() const;
  // Number of branch-and-bound nodes. Only available for discrete problems.
  virtual int64_t nodes() const;

  // Returns the basis status of a row.
  virtual MPSolver::BasisStatus row_status(int constraint_index) const;
  // Returns the basis status of a column.
  virtual MPSolver::BasisStatus column_status(int variable_index) const;

  // ----- Misc -----

  // Query problem type.
  // Remember that problem type is a static property that is set
  // in the constructor and never changed.
  virtual bool IsContinuous() const { return IsLP(); }
  virtual bool IsLP() const { return !mMip; }
  virtual bool IsMIP() const { return mMip; }

  virtual void ExtractNewVariables();
  virtual void ExtractNewConstraints();
  virtual void ExtractObjective();

  virtual std::string SolverVersion() const;

  virtual void* underlying_solver() { return reinterpret_cast<void*>(mLp); }

  virtual double ComputeExactConditionNumber() const {
    if (!IsContinuous()) {
      LOG(DFATAL) << "ComputeExactConditionNumber not implemented for"
                  << " CPLEX_MIXED_INTEGER_PROGRAMMING";
      return CPX_NAN;
    }

    if (CheckSolutionIsSynchronized()) {
      double kappa = CPX_NAN;
      CHECK_STATUS(CPXXgetdblquality(mEnv, mLp, &kappa, CPX_EXACT_KAPPA));
      return kappa;
    }
    LOG(DFATAL) << "Cannot get exact condition number without solution";
    return CPX_NAN;
  }

 protected:
  // Set all parameters in the underlying solver.
  virtual void SetParameters(MPSolverParameters const& param);
  // Set each parameter in the underlying solver.
  virtual void SetRelativeMipGap(double value);
  virtual void SetPrimalTolerance(double value);
  virtual void SetDualTolerance(double value);
  virtual void SetPresolveMode(int value);
  virtual void SetScalingMode(int value);
  virtual void SetLpAlgorithm(int value);

  virtual bool ReadParameterFile(std::string const& filename);
  virtual std::string ValidFileExtensionForParameterFile() const;

 private:
  // Mark modeling object "out of sync". This implicitly invalidates
  // solution information as well. It is the counterpart of
  // MPSolverInterface::InvalidateSolutionSynchronization
  void InvalidateModelSynchronization() {
    mCstat = 0;
    mRstat = 0;
    sync_status_ = MUST_RELOAD;
  }

  // Transform CPLEX basis status to MPSolver basis status.
  static MPSolver::BasisStatus xformBasisStatus(int cplex_basis_status);

 private:
  CPXLPptr mLp;
  CPXENVptr mEnv;
  bool const mMip;
  // Incremental extraction.
  // Without incremental extraction we have to re-extract the model every
  // time we perform a solve. Due to the way the Reset() function is
  // implemented, this will lose MIP start or basis information from a
  // previous solve. On the other hand, if there is a significant changes
  // to the model then just re-extracting everything is usually faster than
  // keeping the low-level modeling object in sync with the high-level
  // variables/constraints.
  // Note that incremental extraction is particularly expensive in function
  // ExtractNewVariables() since there we must scan _all_ old constraints
  // and update them with respect to the new variables.
  bool const supportIncrementalExtraction;

  // Use slow and immediate updates or try to do bulk updates.
  // For many updates to the model we have the option to either perform
  // the update immediately with a potentially slow operation or to
  // just mark the low-level modeling object out of sync and re-extract
  // the model later.
  enum SlowUpdates {
    SlowSetCoefficient = 0x0001,
    SlowClearConstraint = 0x0002,
    SlowSetObjectiveCoefficient = 0x0004,
    SlowClearObjective = 0x0008,
    SlowSetConstraintBounds = 0x0010,
    SlowSetVariableInteger = 0x0020,
    SlowSetVariableBounds = 0x0040,
    SlowUpdatesAll = 0xffff
  } const slowUpdates;
  // CPLEX has no method to query the basis status of a single variable.
  // Hence we query the status only once and cache the array. This is
  // much faster in case the basis status of more than one row/column
  // is required.
  unique_ptr<int[]> mutable mCstat;
  unique_ptr<int[]> mutable mRstat;

  // Setup the right-hand side of a constraint from its lower and upper bound.
  static void MakeRhs(double lb, double ub, double& rhs, char& sense,
                      double& range);
};

// Creates a LP/MIP instance.
CplexInterface::CplexInterface(MPSolver* const solver, bool mip)
    : MPSolverInterface(solver),
      mEnv(0),
      mLp(0),
      mMip(mip),
      slowUpdates(static_cast<SlowUpdates>(SlowSetObjectiveCoefficient |
                                           SlowClearObjective)),
      supportIncrementalExtraction(false),
      mCstat(),
      mRstat() {
  int status;

  mEnv = CPXXopenCPLEX(&status);
  CHECK_STATUS(status);
  DCHECK(mEnv != nullptr);  // should not be NULL if status=0

  char const* name = solver_->name_.c_str();
  mLp = CPXXcreateprob(mEnv, &status, name);
  CHECK_STATUS(status);
  DCHECK(mLp != nullptr);  // should not be NULL if status=0

  CHECK_STATUS(CPXXchgobjsen(mEnv, mLp, maximize_ ? CPX_MAX : CPX_MIN));
  if (mMip) CHECK_STATUS(CPXXchgprobtype(mEnv, mLp, CPXPROB_MILP));
}

CplexInterface::~CplexInterface() {
  CHECK_STATUS(CPXXfreeprob(mEnv, &mLp));
  CHECK_STATUS(CPXXcloseCPLEX(&mEnv));
}

std::string CplexInterface::SolverVersion() const {
  // We prefer CPXXversionnumber() over CPXXversion() since the
  // former will never pose any encoding issues.
  int version = 0;
  CHECK_STATUS(CPXXversionnumber(mEnv, &version));

  int const major = version / 1000000;
  version -= major * 1000000;
  int const release = version / 10000;
  version -= release * 10000;
  int const mod = version / 100;
  version -= mod * 100;
  int const fix = version;

  return absl::StrFormat("CPLEX library version %d.%02d.%02d.%02d", major,
                         release, mod, fix);
}

// ------ Model modifications and extraction -----

void CplexInterface::Reset() {
  // Instead of explicitly clearing all modeling objects we
  // just delete the problem object and allocate a new one.
  CHECK_STATUS(CPXXfreeprob(mEnv, &mLp));

  int status;
  const char* const name = solver_->name_.c_str();
  mLp = CPXXcreateprob(mEnv, &status, name);
  CHECK_STATUS(status);
  DCHECK(mLp != nullptr);  // should not be NULL if status=0

  CHECK_STATUS(CPXXchgobjsen(mEnv, mLp, maximize_ ? CPX_MAX : CPX_MIN));
  if (mMip) CHECK_STATUS(CPXXchgprobtype(mEnv, mLp, CPXPROB_MILP));

  ResetExtractionInformation();
  mCstat = 0;
  mRstat = 0;
}

void CplexInterface::SetOptimizationDirection(bool maximize) {
  InvalidateSolutionSynchronization();
  CPXXchgobjsen(mEnv, mLp, maximize ? CPX_MAX : CPX_MIN);
}

void CplexInterface::SetVariableBounds(int var_index, double lb, double ub) {
  InvalidateSolutionSynchronization();

  // Changing the bounds of a variable is fast. However, doing this for
  // many variables may still be slow. So we don't perform the update by
  // default. However, if we support incremental extraction
  // (supportIncrementalExtraction is true) then we MUST perform the
  // update here or we will lose it.

  if (!supportIncrementalExtraction && !(slowUpdates & SlowSetVariableBounds)) {
    InvalidateModelSynchronization();
  } else {
    if (variable_is_extracted(var_index)) {
      // Variable has already been extracted, so we must modify the
      // modeling object.
      DCHECK_LT(var_index, last_variable_index_);
      char const lu[2] = {'L', 'U'};
      double const bd[2] = {lb, ub};
      CPXDIM const idx[2] = {var_index, var_index};
      CHECK_STATUS(CPXXchgbds(mEnv, mLp, 2, idx, lu, bd));
    } else {
      // Variable is not yet extracted. It is sufficient to just mark
      // the modeling object "out of sync"
      InvalidateModelSynchronization();
    }
  }
}

// Modifies integrality of an extracted variable.
void CplexInterface::SetVariableInteger(int var_index, bool integer) {
  InvalidateSolutionSynchronization();

  // NOTE: The type of the model (continuous or mixed integer) is
  //       defined once and for all in the constructor. There are no
  //       dynamic changes to the model type.

  // Changing the type of a variable should be fast. Still, doing all
  // updates in one big chunk right before solve() is usually faster.
  // However, if we support incremental extraction
  // (supportIncrementalExtraction is true) then we MUST change the
  // type of extracted variables here.

  if (!supportIncrementalExtraction &&
      !(slowUpdates && SlowSetVariableInteger)) {
    InvalidateModelSynchronization();
  } else {
    if (mMip) {
      if (variable_is_extracted(var_index)) {
        // Variable is extracted. Change the type immediately.
        // TODO: Should we check the current type and don't do anything
        //       in case the type does not change?
        DCHECK_LE(var_index, CPXXgetnumcols(mEnv, mLp));
        char const type = integer ? CPX_INTEGER : CPX_CONTINUOUS;
        CHECK_STATUS(CPXXchgctype(mEnv, mLp, 1, &var_index, &type));
      } else
        InvalidateModelSynchronization();
    } else {
      LOG(DFATAL)
          << "Attempt to change variable to integer in non-MIP problem!";
    }
  }
}

// Setup the right-hand side of a constraint.
void CplexInterface::MakeRhs(double lb, double ub, double& rhs, char& sense,
                             double& range) {
  if (lb == ub) {
    // Both bounds are equal -> this is an equality constraint
    rhs = lb;
    range = 0.0;
    sense = 'E';
  } else if (lb > -CPX_INFBOUND && ub < CPX_INFBOUND) {
    // Both bounds are finite -> this is a ranged constraint
    // The value of a ranged constraint is allowed to be in
    //   [ rhs[i], rhs[i]+rngval[i] ]
    // see also the reference documentation for CPXXnewrows()
    if (ub < lb) {
      // The bounds for the constraint are contradictory. CPLEX models
      // a range constraint l <= ax <= u as
      //    ax = l + v
      // where v is an auxiliary variable the range of which is controlled
      // by l and u: if l < u then v in [0, u-l]
      //             else          v in [u-l, 0]
      // (the range is specified as the rngval[] argument to CPXXnewrows).
      // Thus CPLEX cannot represent range constraints with contradictory
      // bounds and we must error out here.
      CHECK_STATUS(CPXERR_BAD_ARGUMENT);
    }
    rhs = lb;
    range = ub - lb;
    sense = 'R';
  } else if (ub < CPX_INFBOUND ||
             (std::abs(ub) == CPX_INFBOUND && std::abs(lb) > CPX_INFBOUND)) {
    // Finite upper, infinite lower bound -> this is a <= constraint
    rhs = ub;
    range = 0.0;
    sense = 'L';
  } else if (lb > -CPX_INFBOUND ||
             (std::abs(lb) == CPX_INFBOUND && std::abs(ub) > CPX_INFBOUND)) {
    // Finite lower, infinite upper bound -> this is a >= constraint
    rhs = lb;
    range = 0.0;
    sense = 'G';
  } else {
    // Lower and upper bound are both infinite.
    // This is used for example in .mps files to specify alternate
    // objective functions.
    // Note that the case lb==ub was already handled above, so we just
    // pick the bound with larger magnitude and create a constraint for it.
    // Note that we replace the infinite bound by CPX_INFBOUND since
    // bounds with larger magnitude may cause other CPLEX functions to
    // fail (for example the export to LP files).
    DCHECK_GT(std::abs(lb), CPX_INFBOUND);
    DCHECK_GT(std::abs(ub), CPX_INFBOUND);
    if (std::abs(lb) > std::abs(ub)) {
      rhs = (lb < 0) ? -CPX_INFBOUND : CPX_INFBOUND;
      sense = 'G';
    } else {
      rhs = (ub < 0) ? -CPX_INFBOUND : CPX_INFBOUND;
      sense = 'L';
    }
    range = 0.0;
  }
}

void CplexInterface::SetConstraintBounds(int index, double lb, double ub) {
  InvalidateSolutionSynchronization();

  // Changing rhs, sense, or range of a constraint is not too slow.
  // Still, doing all the updates in one large operation is faster.
  // Note however that if we do not want to re-extract the full model
  // for each solve (supportIncrementalExtraction is true) then we MUST
  // update the constraint here, otherwise we lose this update information.

  if (!supportIncrementalExtraction &&
      !(slowUpdates & SlowSetConstraintBounds)) {
    InvalidateModelSynchronization();
  } else {
    if (constraint_is_extracted(index)) {
      // Constraint is already extracted, so we must update its bounds
      // and its type.
      DCHECK(mLp != NULL);
      char sense;
      double range, rhs;
      MakeRhs(lb, ub, rhs, sense, range);
      CHECK_STATUS(CPXXchgrhs(mEnv, mLp, 1, &index, &lb));
      CHECK_STATUS(CPXXchgsense(mEnv, mLp, 1, &index, &sense));
      CHECK_STATUS(CPXXchgrngval(mEnv, mLp, 1, &index, &range));
    } else {
      // Constraint is not yet extracted. It is sufficient to mark the
      // modeling object as "out of sync"
      InvalidateModelSynchronization();
    }
  }
}

void CplexInterface::AddRowConstraint(MPConstraint* const ct) {
  // This is currently only invoked when a new constraint is created,
  // see MPSolver::MakeRowConstraint().
  // At this point we only have the lower and upper bounds of the
  // constraint. We could immediately call CPXXaddrows() here but it is
  // usually much faster to handle the fully populated constraint in
  // ExtractNewConstraints() right before the solve.
  InvalidateModelSynchronization();
}

void CplexInterface::AddVariable(MPVariable* const ct) {
  // This is currently only invoked when a new variable is created,
  // see MPSolver::MakeVar().
  // At this point the variable does not appear in any constraints or
  // the objective function. We could invoke CPXXaddcols() to immediately
  // create the variable here but it is usually much faster to handle the
  // fully setup variable in ExtractNewVariables() right before the solve.
  InvalidateModelSynchronization();
}

void CplexInterface::SetCoefficient(MPConstraint* const constraint,
                                    MPVariable const* const variable,
                                    double new_value, double) {
  InvalidateSolutionSynchronization();

  // Changing a single coefficient in the matrix is potentially pretty
  // slow since that coefficient has to be found in the sparse matrix
  // representation. So by default we don't perform this update immediately
  // but instead mark the low-level modeling object "out of sync".
  // If we want to support incremental extraction then we MUST perform
  // the modification immediately or we will lose it.

  if (!supportIncrementalExtraction && !(slowUpdates & SlowSetCoefficient)) {
    InvalidateModelSynchronization();
  } else {
    int const row = constraint->index();
    int const col = variable->index();
    if (constraint_is_extracted(row) && variable_is_extracted(col)) {
      // If row and column are both extracted then we can directly
      // update the modeling object
      DCHECK_LE(row, last_constraint_index_);
      DCHECK_LE(col, last_variable_index_);
      CHECK_STATUS(CPXXchgcoef(mEnv, mLp, row, col, new_value));
    } else {
      // If either row or column is not yet extracted then we can
      // defer the update to ExtractModel()
      InvalidateModelSynchronization();
    }
  }
}

void CplexInterface::ClearConstraint(MPConstraint* const constraint) {
  CPXDIM const row = constraint->index();
  if (!constraint_is_extracted(row))
    // There is nothing to do if the constraint was not even extracted.
    return;

  // Clearing a constraint means setting all coefficients in the corresponding
  // row to 0 (we cannot just delete the row since that would renumber all
  // the constraints/rows after it).
  // Modifying coefficients in the matrix is potentially pretty expensive
  // since they must be found in the sparse matrix representation. That is
  // why by default we do not modify the coefficients here but only mark
  // the low-level modeling object "out of sync".

  if (!(slowUpdates & SlowClearConstraint)) {
    InvalidateModelSynchronization();
  } else {
    InvalidateSolutionSynchronization();

    CPXDIM const len = constraint->coefficients_.size();
    unique_ptr<CPXDIM[]> rowind(new CPXDIM[len]);
    unique_ptr<CPXDIM[]> colind(new CPXDIM[len]);
    unique_ptr<double[]> val(new double[len]);
    CPXDIM j = 0;
    const auto& coeffs = constraint->coefficients_;
    for (auto it(coeffs.begin()); it != coeffs.end(); ++it) {
      CPXDIM const col = it->first->index();
      if (variable_is_extracted(col)) {
        rowind[j] = row;
        colind[j] = col;
        val[j] = 0.0;
        ++j;
      }
    }
    if (j)
      CHECK_STATUS(
          CPXXchgcoeflist(mEnv, mLp, j, rowind.get(), colind.get(), val.get()));
  }
}

void CplexInterface::SetObjectiveCoefficient(MPVariable const* const variable,
                                             double coefficient) {
  CPXDIM const col = variable->index();
  if (!variable_is_extracted(col))
    // Nothing to do if variable was not even extracted
    return;

  InvalidateSolutionSynchronization();

  // The objective function is stored as a dense vector, so updating a
  // single coefficient is O(1). So by default we update the low-level
  // modeling object here.
  // If we support incremental extraction then we have no choice but to
  // perform the update immediately.

  if (supportIncrementalExtraction ||
      (slowUpdates & SlowSetObjectiveCoefficient)) {
    CHECK_STATUS(CPXXchgobj(mEnv, mLp, 1, &col, &coefficient));
  } else
    InvalidateModelSynchronization();
}

void CplexInterface::SetObjectiveOffset(double value) {
  // Changing the objective offset is O(1), so we always do it immediately.
  InvalidateSolutionSynchronization();
  CHECK_STATUS(CPXEsetobjoffset(mEnv, mLp, value));
}

void CplexInterface::ClearObjective() {
  InvalidateSolutionSynchronization();

  // Since the objective function is stored as a dense vector updating
  // it is O(n), so we usually perform the update immediately.
  // If we want to support incremental extraction then we have no choice
  // but to perform the update immediately.

  if (supportIncrementalExtraction || (slowUpdates & SlowClearObjective)) {
    CPXDIM const cols = CPXXgetnumcols(mEnv, mLp);
    unique_ptr<CPXDIM[]> ind(new CPXDIM[cols]);
    unique_ptr<double[]> zero(new double[cols]);
    CPXDIM j = 0;
    const auto& coeffs = solver_->objective_->coefficients_;
    for (auto it(coeffs.begin()); it != coeffs.end(); ++it) {
      CPXDIM const idx = it->first->index();
      // We only need to reset variables that have been extracted.
      if (variable_is_extracted(idx)) {
        DCHECK_LT(idx, cols);
        ind[j] = idx;
        zero[j] = 0.0;
        ++j;
      }
    }
    if (j > 0) CHECK_STATUS(CPXXchgobj(mEnv, mLp, j, ind.get(), zero.get()));
    CHECK_STATUS(CPXEsetobjoffset(mEnv, mLp, 0.0));
  } else
    InvalidateModelSynchronization();
}

// ------ Query statistics on the solution and the solve ------

int64_t CplexInterface::iterations() const {
  int iter;
  if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfIterations;
  if (mMip)
    return static_cast<int64_t>(CPXXgetmipitcnt(mEnv, mLp));
  else
    return static_cast<int64_t>(CPXXgetitcnt(mEnv, mLp));
}

int64_t CplexInterface::nodes() const {
  if (mMip) {
    if (!CheckSolutionIsSynchronized()) return kUnknownNumberOfNodes;
    return static_cast<int64_t>(CPXXgetnodecnt(mEnv, mLp));
  } else {
    LOG(DFATAL) << "Number of nodes only available for discrete problems";
    return kUnknownNumberOfNodes;
  }
}

// Transform a CPLEX basis status to an MPSolver basis status.
MPSolver::BasisStatus CplexInterface::xformBasisStatus(int cplex_basis_status) {
  switch (cplex_basis_status) {
    case CPX_AT_LOWER:
      return MPSolver::AT_LOWER_BOUND;
    case CPX_BASIC:
      return MPSolver::BASIC;
    case CPX_AT_UPPER:
      return MPSolver::AT_UPPER_BOUND;
    case CPX_FREE_SUPER:
      return MPSolver::FREE;
    default:
      LOG(DFATAL) << "Unknown CPLEX basis status";
      return MPSolver::FREE;
  }
}

// Returns the basis status of a row.
MPSolver::BasisStatus CplexInterface::row_status(int constraint_index) const {
  if (mMip) {
    LOG(FATAL) << "Basis status only available for continuous problems";
    return MPSolver::FREE;
  }

  if (CheckSolutionIsSynchronized()) {
    if (!mRstat) {
      CPXDIM const rows = CPXXgetnumrows(mEnv, mLp);
      unique_ptr<int[]> data(new int[rows]);
      mRstat.swap(data);
      CHECK_STATUS(CPXXgetbase(mEnv, mLp, 0, mRstat.get()));
    }
  } else
    mRstat = 0;

  if (mRstat)
    return xformBasisStatus(mRstat[constraint_index]);
  else {
    LOG(FATAL) << "Row basis status not available";
    return MPSolver::FREE;
  }
}

// Returns the basis status of a column.
MPSolver::BasisStatus CplexInterface::column_status(int variable_index) const {
  if (mMip) {
    LOG(FATAL) << "Basis status only available for continuous problems";
    return MPSolver::FREE;
  }

  if (CheckSolutionIsSynchronized()) {
    if (!mCstat) {
      CPXDIM const cols = CPXXgetnumcols(mEnv, mLp);
      unique_ptr<int[]> data(new int[cols]);
      mCstat.swap(data);
      CHECK_STATUS(CPXXgetbase(mEnv, mLp, mCstat.get(), 0));
    }
  } else
    mCstat = 0;

  if (mCstat)
    return xformBasisStatus(mCstat[variable_index]);
  else {
    LOG(FATAL) << "Column basis status not available";
    return MPSolver::FREE;
  }
}

// Extract all variables that have not yet been extracted.
void CplexInterface::ExtractNewVariables() {
  // NOTE: The code assumes that a linear expression can never contain
  //       non-zero duplicates.

  InvalidateSolutionSynchronization();

  if (!supportIncrementalExtraction) {
    // Without incremental extraction ExtractModel() is always called
    // to extract the full model.
    CHECK(last_variable_index_ == 0 ||
          last_variable_index_ == solver_->variables_.size());
    CHECK(last_constraint_index_ == 0 ||
          last_constraint_index_ == solver_->constraints_.size());
  }

  int const last_extracted = last_variable_index_;
  int const var_count = solver_->variables_.size();
  CPXDIM newcols = var_count - last_extracted;
  if (newcols > 0) {
    // There are non-extracted variables. Extract them now.

    unique_ptr<double[]> obj(new double[newcols]);
    unique_ptr<double[]> lb(new double[newcols]);
    unique_ptr<double[]> ub(new double[newcols]);
    unique_ptr<char[]> ctype(new char[newcols]);
    unique_ptr<const char*[]> colname(new const char*[newcols]);

    bool have_names = false;
    for (int j = 0, varidx = last_extracted; j < newcols; ++j, ++varidx) {
      MPVariable const* const var = solver_->variables_[varidx];
      lb[j] = var->lb();
      ub[j] = var->ub();
      ctype[j] = var->integer() ? CPX_INTEGER : CPX_CONTINUOUS;
      colname[j] = var->name().empty() ? 0 : var->name().c_str();
      have_names = have_names || var->name().empty();
      obj[j] = solver_->objective_->GetCoefficient(var);
    }

    // Arrays for modifying the problem are setup. Update the index
    // of variables that will get extracted now. Updating indices
    // _before_ the actual extraction makes things much simpler in
    // case we support incremental extraction.
    // In case of error we just reset the indeces.
    std::vector<MPVariable*> const& variables = solver_->variables();
    for (int j = last_extracted; j < var_count; ++j) {
      CHECK(!variable_is_extracted(variables[j]->index()));
      set_variable_as_extracted(variables[j]->index(), true);
    }

    try {
      bool use_newcols = true;

      if (supportIncrementalExtraction) {
        // If we support incremental extraction then we must
        // update existing constraints with the new variables.
        // To do that we use CPXXaddcols() to actually create the
        // variables. This is supposed to be faster than combining
        // CPXXnewcols() and CPXXchgcoeflist().

        // For each column count the size of the intersection with
        // existing constraints.
        unique_ptr<CPXDIM[]> collen(new CPXDIM[newcols]);
        for (CPXDIM j = 0; j < newcols; ++j) collen[j] = 0;
        CPXNNZ nonzeros = 0;
        // TODO: Use a bitarray to flag the constraints that actually
        //       intersect new variables?
        for (int i = 0; i < last_constraint_index_; ++i) {
          MPConstraint const* const ct = solver_->constraints_[i];
          CHECK(constraint_is_extracted(ct->index()));
          const auto& coeffs = ct->coefficients_;
          for (auto it(coeffs.begin()); it != coeffs.end(); ++it) {
            int const idx = it->first->index();
            if (variable_is_extracted(idx) && idx > last_variable_index_) {
              collen[idx - last_variable_index_]++;
              ++nonzeros;
            }
          }
        }

        if (nonzeros > 0) {
          // At least one of the new variables did intersect with an
          // old constraint. We have to create the new columns via
          // CPXXaddcols().
          use_newcols = false;
          unique_ptr<CPXNNZ[]> begin(new CPXNNZ[newcols + 2]);
          unique_ptr<CPXDIM[]> cmatind(new CPXDIM[nonzeros]);
          unique_ptr<double[]> cmatval(new double[nonzeros]);

          // Here is how cmatbeg[] is setup:
          // - it is initialized as
          //     [ 0, 0, collen[0], collen[0]+collen[1], ... ]
          //   so that cmatbeg[j+1] tells us where in cmatind[] and
          //   cmatval[] we need to put the next nonzero for column
          //   j
          // - after nonzeros have been setup the array looks like
          //     [ 0, collen[0], collen[0]+collen[1], ... ]
          //   so that it is the correct input argument for CPXXaddcols
          CPXNNZ* cmatbeg = begin.get();
          cmatbeg[0] = 0;
          cmatbeg[1] = 0;
          ++cmatbeg;
          for (CPXDIM j = 0; j < newcols; ++j)
            cmatbeg[j + 1] = cmatbeg[j] + collen[j];

          for (int i = 0; i < last_constraint_index_; ++i) {
            MPConstraint const* const ct = solver_->constraints_[i];
            CPXDIM const row = ct->index();
            const auto& coeffs = ct->coefficients_;
            for (auto it(coeffs.begin()); it != coeffs.end(); ++it) {
              int const idx = it->first->index();
              if (variable_is_extracted(idx) && idx > last_variable_index_) {
                CPXNNZ const nz = cmatbeg[idx]++;
                cmatind[nz] = idx;
                cmatval[nz] = it->second;
              }
            }
          }
          --cmatbeg;
          CHECK_STATUS(CPXXaddcols(mEnv, mLp, newcols, nonzeros, obj.get(),
                                   cmatbeg, cmatind.get(), cmatval.get(),
                                   lb.get(), ub.get(),
                                   have_names ? colname.get() : 0));
        }
      }
      if (use_newcols) {
        // Either incremental extraction is not supported or none of
        // the new variables did intersect an existing constraint.
        // We can just use CPXXnewcols() to create the new variables.
        CHECK_STATUS(CPXXnewcols(mEnv, mLp, newcols, obj.get(), lb.get(),
                                 ub.get(), mMip ? ctype.get() : 0,
                                 have_names ? colname.get() : 0));
      } else {
        // Incremental extraction: we must update the ctype of the
        // newly created variables (CPXXaddcols() does not allow
        // specifying the ctype)
        if (mMip) {
          // Query the actual number of columns in case we did not
          // manage to extract all columns.
          int const cols = CPXXgetnumcols(mEnv, mLp);
          unique_ptr<CPXDIM[]> ind(new CPXDIM[newcols]);
          for (int j = last_extracted; j < cols; ++j)
            ind[j - last_extracted] = j;
          CHECK_STATUS(CPXXchgctype(mEnv, mLp, cols - last_extracted, ind.get(),
                                    ctype.get()));
        }
      }
    } catch (...) {
      // Undo all changes in case of error.
      CPXDIM const cols = CPXXgetnumcols(mEnv, mLp);
      if (cols > last_extracted)
        (void)CPXXdelcols(mEnv, mLp, last_extracted, cols - 1);
      std::vector<MPVariable*> const& variables = solver_->variables();
      int const size = variables.size();
      for (int j = last_extracted; j < size; ++j)
        set_variable_as_extracted(j, false);
      throw;
    }
  }
}

// Extract constraints that have not yet been extracted.
void CplexInterface::ExtractNewConstraints() {
  // NOTE: The code assumes that a linear expression can never contain
  //       non-zero duplicates.

  if (!supportIncrementalExtraction) {
    // Without incremental extraction ExtractModel() is always called
    // to extract the full model.
    CHECK(last_variable_index_ == 0 ||
          last_variable_index_ == solver_->variables_.size());
    CHECK(last_constraint_index_ == 0 ||
          last_constraint_index_ == solver_->constraints_.size());
  }

  CPXDIM const offset = last_constraint_index_;
  CPXDIM const total = solver_->constraints_.size();

  if (total > offset) {
    // There are constraints that are not yet extracted.

    InvalidateSolutionSynchronization();

    CPXDIM newCons = total - offset;
    CPXDIM const cols = CPXXgetnumcols(mEnv, mLp);
    DCHECK_EQ(last_variable_index_, cols);
    CPXDIM const chunk = 10;  // max number of rows to add in one shot

    // Update indices of new constraints _before_ actually extracting
    // them. In case of error we will just reset the indices.
    for (CPXDIM c = offset; c < total; ++c)
      set_constraint_as_extracted(c, true);

    try {
      unique_ptr<CPXDIM[]> rmatind(new CPXDIM[cols]);
      unique_ptr<double[]> rmatval(new double[cols]);
      unique_ptr<CPXNNZ[]> rmatbeg(new CPXNNZ[chunk]);
      unique_ptr<char[]> sense(new char[chunk]);
      unique_ptr<double[]> rhs(new double[chunk]);
      unique_ptr<char const*[]> name(new char const*[chunk]);
      unique_ptr<double[]> rngval(new double[chunk]);
      unique_ptr<CPXDIM[]> rngind(new CPXDIM[chunk]);
      bool haveRanges = false;

      // Loop over the new constraints, collecting rows for up to
      // CHUNK constraints into the arrays so that adding constraints
      // is faster.
      for (CPXDIM c = 0; c < newCons; /* nothing */) {
        // Collect up to CHUNK constraints into the arrays.
        CPXDIM nextRow = 0;
        CPXNNZ nextNz = 0;
        for (/* nothing */; c < newCons && nextRow < chunk; ++c, ++nextRow) {
          MPConstraint const* const ct = solver_->constraints_[offset + c];

          // Stop if there is not enough room in the arrays
          // to add the current constraint.
          if (nextNz + ct->coefficients_.size() > cols) {
            DCHECK_GT(nextRow, 0);
            break;
          }

          // Setup right-hand side of constraint.
          MakeRhs(ct->lb(), ct->ub(), rhs[nextRow], sense[nextRow],
                  rngval[nextRow]);
          haveRanges = haveRanges || (rngval[nextRow] != 0.0);
          rngind[nextRow] = offset + c;

          // Setup left-hand side of constraint.
          rmatbeg[nextRow] = nextNz;
          const auto& coeffs = ct->coefficients_;
          for (auto it(coeffs.begin()); it != coeffs.end(); ++it) {
            CPXDIM const idx = it->first->index();
            if (variable_is_extracted(idx)) {
              DCHECK_LT(nextNz, cols);
              DCHECK_LT(idx, cols);
              rmatind[nextNz] = idx;
              rmatval[nextNz] = it->second;
              ++nextNz;
            }
          }

          // Finally the name of the constraint.
          name[nextRow] = ct->name().empty() ? 0 : ct->name().c_str();
        }
        if (nextRow > 0) {
          CHECK_STATUS(CPXXaddrows(mEnv, mLp, 0, nextRow, nextNz, rhs.get(),
                                   sense.get(), rmatbeg.get(), rmatind.get(),
                                   rmatval.get(), 0, name.get()));
          if (haveRanges) {
            CHECK_STATUS(
                CPXXchgrngval(mEnv, mLp, nextRow, rngind.get(), rngval.get()));
          }
        }
      }
    } catch (...) {
      // Undo all changes in case of error.
      CPXDIM const rows = CPXXgetnumrows(mEnv, mLp);
      if (rows > offset) (void)CPXXdelrows(mEnv, mLp, offset, rows - 1);
      std::vector<MPConstraint*> const& constraints = solver_->constraints();
      int const size = constraints.size();
      for (int i = offset; i < size; ++i) set_constraint_as_extracted(i, false);
      throw;
    }
  }
}

// Extract the objective function.
void CplexInterface::ExtractObjective() {
  // NOTE: The code assumes that the objective expression does not contain
  //       any non-zero duplicates.

  CPXDIM const cols = CPXXgetnumcols(mEnv, mLp);
  DCHECK_EQ(last_variable_index_, cols);

  unique_ptr<CPXDIM[]> ind(new CPXDIM[cols]);
  unique_ptr<double[]> val(new double[cols]);
  for (CPXDIM j = 0; j < cols; ++j) {
    ind[j] = j;
    val[j] = 0.0;
  }

  const auto& coeffs = solver_->objective_->coefficients_;
  for (auto it = coeffs.begin(); it != coeffs.end(); ++it) {
    CPXDIM const idx = it->first->index();
    if (variable_is_extracted(idx)) {
      DCHECK_LT(idx, cols);
      val[idx] = it->second;
    }
  }

  CHECK_STATUS(CPXXchgobj(mEnv, mLp, cols, ind.get(), val.get()));
  CHECK_STATUS(CPXEsetobjoffset(mEnv, mLp, solver_->Objective().offset()));
}

// ------ Parameters  -----

void CplexInterface::SetParameters(const MPSolverParameters& param) {
  SetCommonParameters(param);
  if (mMip) SetMIPParameters(param);
}

void CplexInterface::SetRelativeMipGap(double value) {
  if (mMip) {
    CHECK_STATUS(CPXXsetdblparam(mEnv, CPX_PARAM_EPGAP, value));
  } else {
    LOG(WARNING) << "The relative MIP gap is only available "
                 << "for discrete problems.";
  }
}

void CplexInterface::SetPrimalTolerance(double value) {
  CHECK_STATUS(CPXXsetdblparam(mEnv, CPX_PARAM_EPRHS, value));
}

void CplexInterface::SetDualTolerance(double value) {
  CHECK_STATUS(CPXXsetdblparam(mEnv, CPX_PARAM_EPOPT, value));
}

void CplexInterface::SetPresolveMode(int value) {
  MPSolverParameters::PresolveValues const presolve =
      static_cast<MPSolverParameters::PresolveValues>(value);

  switch (presolve) {
    case MPSolverParameters::PRESOLVE_OFF:
      CHECK_STATUS(CPXXsetintparam(mEnv, CPX_PARAM_PREIND, CPX_OFF));
      return;
    case MPSolverParameters::PRESOLVE_ON:
      CHECK_STATUS(CPXXsetintparam(mEnv, CPX_PARAM_PREIND, CPX_ON));
      return;
  }
  SetIntegerParamToUnsupportedValue(MPSolverParameters::PRESOLVE, value);
}

// Sets the scaling mode.
void CplexInterface::SetScalingMode(int value) {
  MPSolverParameters::ScalingValues const scaling =
      static_cast<MPSolverParameters::ScalingValues>(value);

  switch (scaling) {
    case MPSolverParameters::SCALING_OFF:
      CHECK_STATUS(CPXXsetintparam(mEnv, CPX_PARAM_SCAIND, -1));
      break;
    case MPSolverParameters::SCALING_ON:
      // TODO: 0 is equilibrium scaling (the default), CPLEX also supports
      //       1 aggressive scaling
      CHECK_STATUS(CPXXsetintparam(mEnv, CPX_PARAM_SCAIND, 0));
      break;
  }
}

// Sets the LP algorithm : primal, dual or barrier. Note that CPLEX offers other
// LP algorithm (e.g. network) and automatic selection
void CplexInterface::SetLpAlgorithm(int value) {
  MPSolverParameters::LpAlgorithmValues const algorithm =
      static_cast<MPSolverParameters::LpAlgorithmValues>(value);

  int alg = CPX_ALG_NONE;

  switch (algorithm) {
    case MPSolverParameters::DUAL:
      alg = CPX_ALG_DUAL;
      break;
    case MPSolverParameters::PRIMAL:
      alg = CPX_ALG_PRIMAL;
      break;
    case MPSolverParameters::BARRIER:
      alg = CPX_ALG_BARRIER;
      break;
  }

  if (alg == CPX_ALG_NONE)
    SetIntegerParamToUnsupportedValue(MPSolverParameters::LP_ALGORITHM, value);
  else {
    CHECK_STATUS(CPXXsetintparam(mEnv, CPX_PARAM_LPMETHOD, alg));
    if (mMip) {
      // For MIP we have to change two more parameters to specify the
      // algorithm that is used to solve LP relaxations.
      CHECK_STATUS(CPXXsetintparam(mEnv, CPX_PARAM_STARTALG, alg));
      CHECK_STATUS(CPXXsetintparam(mEnv, CPX_PARAM_SUBALG, alg));
    }
  }
}

bool CplexInterface::ReadParameterFile(std::string const& filename) {
  // Return true on success and false on error.
  return CPXXreadcopyparam(mEnv, filename.c_str()) == 0;
}

std::string CplexInterface::ValidFileExtensionForParameterFile() const {
  return ".prm";
}

MPSolver::ResultStatus CplexInterface::Solve(MPSolverParameters const& param) {
  int status;

  // Delete chached information
  mCstat = 0;
  mRstat = 0;

  WallTimer timer;
  timer.Start();

  // Set incrementality
  MPSolverParameters::IncrementalityValues const inc =
      static_cast<MPSolverParameters::IncrementalityValues>(
          param.GetIntegerParam(MPSolverParameters::INCREMENTALITY));
  switch (inc) {
    case MPSolverParameters::INCREMENTALITY_OFF:
      Reset(); /* This should not be required but re-extracting everything
                * may be faster, so we do it. */
      CHECK_STATUS(CPXXsetintparam(mEnv, CPX_PARAM_ADVIND, 0));
      break;
    case MPSolverParameters::INCREMENTALITY_ON:
      CHECK_STATUS(CPXXsetintparam(mEnv, CPX_PARAM_ADVIND, 2));
      break;
  }

  // Extract the model to be solved.
  // If we don't support incremental extraction and the low-level modeling
  // is out of sync then we have to re-extract everything. Note that this
  // will lose MIP starts or advanced basis information from a previous
  // solve.
  if (!supportIncrementalExtraction && sync_status_ == MUST_RELOAD) Reset();
  ExtractModel();
  VLOG(1) << absl::StrFormat("Model build in %.3f seconds.", timer.Get());

  // Set log level.
  CHECK_STATUS(
      CPXXsetintparam(mEnv, CPX_PARAM_SCRIND, quiet() ? CPX_OFF : CPX_ON));

  // Set parameters.
  // NOTE: We must invoke SetSolverSpecificParametersAsString() _first_.
  //       Its current implementation invokes ReadParameterFile() which in
  //       turn invokes CPXXreadcopyparam(). The latter will _overwrite_
  //       all current parameter settings in the environment.
  solver_->SetSolverSpecificParametersAsString(
      solver_->solver_specific_parameter_string_);
  SetParameters(param);
  if (solver_->time_limit()) {
    VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms.";
    CHECK_STATUS(
        CPXXsetdblparam(mEnv, CPX_PARAM_TILIM, solver_->time_limit() * 1e-3));
  }

  // Solve.
  // Do not CHECK_STATUS here since some errors (for example CPXERR_NO_MEMORY)
  // still allow us to query useful information.
  timer.Restart();
  if (mMip) {
    status = CPXXmipopt(mEnv, mLp);
  } else {
    status = CPXXlpopt(mEnv, mLp);
  }

  // Disable screen output right after solve
  (void)CPXXsetintparam(mEnv, CPX_PARAM_SCRIND, CPX_OFF);

  if (status) {
    VLOG(1) << absl::StrFormat("Failed to optimize MIP. Error %d", status);
    // NOTE: We do not return immediately since there may be information
    //       to grab (for example an incumbent)
  } else {
    VLOG(1) << absl::StrFormat("Solved in %.3f seconds.", timer.Get());
  }

  int const cpxstat = CPXXgetstat(mEnv, mLp);
  VLOG(1) << absl::StrFormat("CPLEX solution status %d.", cpxstat);

  // Figure out what solution we have.
  int solnmethod, solntype, pfeas, dfeas;
  CHECK_STATUS(CPXXsolninfo(mEnv, mLp, &solnmethod, &solntype, &pfeas, &dfeas));
  bool const feasible = pfeas != 0;

  // Get problem dimensions for solution queries below.
  CPXDIM const rows = CPXXgetnumrows(mEnv, mLp);
  CPXDIM const cols = CPXXgetnumcols(mEnv, mLp);
  DCHECK_EQ(rows, solver_->constraints_.size());
  DCHECK_EQ(cols, solver_->variables_.size());

  // Capture objective function value.
  objective_value_ = CPX_NAN;
  best_objective_bound_ = CPX_NAN;
  if (feasible) {
    CHECK_STATUS(CPXXgetobjval(mEnv, mLp, &objective_value_));
    if (mMip) {
      CHECK_STATUS(CPXXgetbestobjval(mEnv, mLp, &best_objective_bound_));
    }
  }
  VLOG(1) << "objective=" << objective_value_
          << ", bound=" << best_objective_bound_;

  // Capture primal and dual solutions
  if (mMip) {
    // If there is a primal feasible solution then capture it.
    if (feasible) {
      if (cols > 0) {
        unique_ptr<double[]> x(new double[cols]);
        CHECK_STATUS(CPXXgetx(mEnv, mLp, x.get(), 0, cols - 1));
        for (int i = 0; i < solver_->variables_.size(); ++i) {
          MPVariable* const var = solver_->variables_[i];
          var->set_solution_value(x[i]);
          VLOG(3) << var->name() << ": value =" << x[i];
        }
      }
    } else {
      for (int i = 0; i < solver_->variables_.size(); ++i)
        solver_->variables_[i]->set_solution_value(CPX_NAN);
    }

    // MIP does not have duals
    for (int i = 0; i < solver_->variables_.size(); ++i)
      solver_->variables_[i]->set_reduced_cost(CPX_NAN);
    for (int i = 0; i < solver_->constraints_.size(); ++i)
      solver_->constraints_[i]->set_dual_value(CPX_NAN);
  } else {
    // Continuous problem.
    if (cols > 0) {
      unique_ptr<double[]> x(new double[cols]);
      unique_ptr<double[]> dj(new double[cols]);
      if (feasible) CHECK_STATUS(CPXXgetx(mEnv, mLp, x.get(), 0, cols - 1));
      if (dfeas) CHECK_STATUS(CPXXgetdj(mEnv, mLp, dj.get(), 0, cols - 1));
      for (int i = 0; i < solver_->variables_.size(); ++i) {
        MPVariable* const var = solver_->variables_[i];
        var->set_solution_value(x[i]);
        bool value = false, dual = false;

        if (feasible) {
          var->set_solution_value(x[i]);
          value = true;
        } else
          var->set_solution_value(CPX_NAN);
        if (dfeas) {
          var->set_reduced_cost(dj[i]);
          dual = true;
        } else
          var->set_reduced_cost(CPX_NAN);
        VLOG(3) << var->name() << ":"
                << (value ? absl::StrFormat("  value = %f", x[i]) : "")
                << (dual ? absl::StrFormat("  reduced cost = %f", dj[i]) : "");
      }
    }

    if (rows > 0) {
      unique_ptr<double[]> pi(new double[rows]);
      if (dfeas) CHECK_STATUS(CPXXgetpi(mEnv, mLp, pi.get(), 0, rows - 1));
      for (int i = 0; i < solver_->constraints_.size(); ++i) {
        MPConstraint* const ct = solver_->constraints_[i];
        bool dual = false;
        if (dfeas) {
          ct->set_dual_value(pi[i]);
          dual = true;
        } else
          ct->set_dual_value(CPX_NAN);
        VLOG(4) << "row " << ct->index() << ":"
                << (dual ? absl::StrFormat("  dual = %f", pi[i]) : "");
      }
    }
  }

  // Map CPLEX status to more generic solution status in MPSolver
  switch (cpxstat) {
    case CPX_STAT_OPTIMAL:
    case CPXMIP_OPTIMAL:
      result_status_ = MPSolver::OPTIMAL;
      break;
    case CPXMIP_OPTIMAL_TOL:
      // To be consistent with the other solvers.
      result_status_ = MPSolver::OPTIMAL;
      break;
    case CPX_STAT_INFEASIBLE:
    case CPXMIP_INFEASIBLE:
      result_status_ = MPSolver::INFEASIBLE;
      break;
    case CPX_STAT_UNBOUNDED:
    case CPXMIP_UNBOUNDED:
      result_status_ = MPSolver::UNBOUNDED;
      break;
    case CPX_STAT_INForUNBD:
    case CPXMIP_INForUNBD:
      result_status_ = MPSolver::INFEASIBLE;
      break;
    default:
      result_status_ = feasible ? MPSolver::FEASIBLE : MPSolver::ABNORMAL;
      break;
  }

  sync_status_ = SOLUTION_SYNCHRONIZED;
  return result_status_;
}

MPSolverInterface* BuildCplexInterface(bool mip, MPSolver* const solver) {
  return new CplexInterface(solver, mip);
}

}  // namespace operations_research
#endif  // #if defined(USE_CPLEX)