#ifndef _PAPILO_INTERFACES_GLOP_INTERFACE_HPP_
#define _PAPILO_INTERFACES_GLOP_INTERFACE_HPP_
#include "papilo/misc/Vec.hpp"
#include <cassert>
#include <stdexcept>
#include "ortools/linear_solver/linear_solver.h"
#include "papilo/core/Problem.hpp"
#include "papilo/core/Solution.hpp"
#include "papilo/interfaces/SolverInterface.hpp"
namespace papilo
{
using operations_research::MPObjective;
using operations_research::MPSolver;
using operations_research::MPConstraint;
using operations_research::MPVariable;
template <typename REAL>
class GlopInterface : public SolverInterface<REAL>
{
private:
int n_cols{};
Vec<MPVariable*> variables;
std::unique_ptr<MPSolver> solver;
std::string solver_id = "PDLP";
int
doSetUp( const Problem<REAL>& problem, const Vec<int>& origRowMap,
const Vec<int>& origColMap )
{
if( !std::is_same<REAL, double>::value )
{
fmt::print( "Please use double precision when solving with "
"OR-TOOLS." );
return -1;
}
auto domains = problem.getVariableDomains();
n_cols = problem.getNCols();
solver = std::unique_ptr<MPSolver>( MPSolver::CreateSolver( solver_id ) );
solver->EnableOutput();
const double infinity = MPSolver::infinity();
auto& row_flags = problem.getRowFlags();
auto lhs = problem.getConstraintMatrix().getLeftHandSides();
auto rhs = problem.getConstraintMatrix().getRightHandSides();
const Vec<String>& varNames = problem.getVariableNames();
auto coefficients = problem.getObjective().coefficients;
variables = Vec<MPVariable*>{};
variables.reserve(problem.getNRows());
MPObjective* const objective = solver->MutableObjective();
for( int i = 0; i < problem.getNCols(); i++ )
{
double lb = domains.flags[i].test( ColFlag::kLbInf )
? -infinity
: domains.lower_bounds[i];
double ub = domains.flags[i].test( ColFlag::kUbInf )
? infinity
: domains.upper_bounds[i];
variables.push_back(solver->MakeNumVar( lb, ub, varNames[origColMap[i]]));
objective->SetCoefficient( variables[i], (double)coefficients[i] );
}
assert( solver->NumVariables() == problem.getNCols() );
for( int i = 0; i < problem.getNRows(); i++ )
{
auto data = problem.getConstraintMatrix().getRowCoefficients( i );
double l = row_flags[i].test( RowFlag::kLhsInf ) ? -infinity : lhs[i];
double r = row_flags[i].test( RowFlag::kRhsInf ) ? infinity : rhs[i];
MPConstraint* const con = solver->MakeRowConstraint( l, r );
for( int j = 0; j < data.getLength(); j++ )
{
double coeff = data.getValues()[j];
int index = data.getIndices()[j];
con->SetCoefficient( variables[index], coeff );
}
}
assert( solver->NumConstraints() == problem.getNRows() );
objective->SetMinimization();
return 0;
}
int
doSetUp( const Problem<REAL>& problem, const Vec<int>& origRowMap,
const Vec<int>& origColMap, const Components& components,
const ComponentInfo& component )
{
if( !std::is_same<REAL, double>::value )
{
fmt::print( "Please use double precision when solving with "
"OR-TOOLS." );
return -1;
}
int ncols = components.getComponentsNumCols( component.componentid );
int nrows = components.getComponentsNumRows( component.componentid );
const int* colset = components.getComponentsCols( component.componentid );
const int* rowset = components.getComponentsRows( component.componentid );
const Vec<String>& varNames = problem.getVariableNames();
const Vec<String>& rowNames = problem.getConstraintNames();
const VariableDomains<REAL>& domains = problem.getVariableDomains();
const Vec<REAL>& obj = problem.getObjective().coefficients;
const Vec<REAL>& rhs = problem.getConstraintMatrix().getRightHandSides();
const Vec<REAL>& lhs = problem.getConstraintMatrix().getLeftHandSides();
const auto consMatrix = problem.getConstraintMatrix();
auto coefficients = problem.getObjective().coefficients;
auto& row_flags = problem.getRowFlags();
n_cols = problem.getNCols();
solver = std::unique_ptr<MPSolver>(MPSolver::CreateSolver( solver_id ));
solver->EnableOutput();
const double infinity = MPSolver::infinity();
variables = Vec<MPVariable*>{};
variables.reserve(problem.getNRows());
MPObjective* const objective = solver->MutableObjective();
for( int j = 0; j < ncols; j++ )
{
int col = colset[j];
double lb = domains.flags[col].test( ColFlag::kLbInf )
? -infinity
: domains.lower_bounds[col];
double ub = domains.flags[col].test( ColFlag::kUbInf )
? infinity
: domains.upper_bounds[col];
variables[col] = solver->MakeNumVar( lb, ub, varNames[origColMap[col]] );
objective->SetCoefficient( variables[col], (double)coefficients[col] );
}
assert( solver->NumVariables() == ncols );
for( int i = 0; i < nrows; i++ )
{
int row = rowset[i];
auto data = problem.getConstraintMatrix().getRowCoefficients( row );
double l = row_flags[row].test( RowFlag::kLhsInf ) ? -infinity : lhs[row];
double r = row_flags[row].test( RowFlag::kRhsInf ) ? infinity : rhs[row];
MPConstraint* const con = solver->MakeRowConstraint( l, r );
for( int j = 0; j < data.getLength(); j++ )
{
double coeff = data.getValues()[j];
int index = components.getColComponentIdx( data.getIndices()[j]);
con->SetCoefficient( variables[index], coeff );
}
}
assert( solver->NumConstraints() == nrows );
objective->SetMinimization();
return 0;
}
public:
GlopInterface() =default;
void
setUp( const Problem<REAL>& prob, const Vec<int>& row_maps,
const Vec<int>& col_maps, const Components& components,
const ComponentInfo& component ) override
{
if( doSetUp( prob, row_maps, col_maps, components, component ) != 0 )
this->status = SolverStatus::kError;
}
void
setNodeLimit( int num ) override
{
}
void
setGapLimit( const REAL& gaplim ) override
{
}
void
setSoftTimeLimit( double tlim ) override
{
}
void
setTimeLimit( double tlim ) override
{
solver->SetTimeLimit(absl::Seconds(tlim));
}
void
setVerbosity( VerbosityLevel verbosity ) override
{
}
void
setUp( const Problem<REAL>& prob, const Vec<int>& row_maps,
const Vec<int>& col_maps ) override
{
if( doSetUp( prob, row_maps, col_maps ) != 0 )
this->status = SolverStatus::kError;
}
void
solve() override
{
MPSolver::ResultStatus result_status = solver->Solve();
assert( this->status != SolverStatus::kError );
switch( result_status )
{
case MPSolver::INFEASIBLE:
this->status = SolverStatus::kInfeasible;
break;
case MPSolver::UNBOUNDED:
this->status = SolverStatus::kUnbounded;
break;
case MPSolver::OPTIMAL:
this->status = SolverStatus::kOptimal;
break;
case MPSolver::FEASIBLE:
case MPSolver::NOT_SOLVED:
this->status = SolverStatus::kInterrupted;
break;
case MPSolver::ABNORMAL:
case MPSolver::MODEL_INVALID:
this->status = SolverStatus::kError;
break;
}
}
REAL
getDualBound() override
{
return 0;
}
bool
is_dual_solution_available() override
{
return false;
}
bool
getSolution( Solution<REAL>& solbuffer ) override
{
Vec<REAL> primal{};
for( int i = 0; i < n_cols; i++ )
primal.push_back( variables[i]->solution_value() );
solbuffer = Solution<REAL>( primal );
return true;
}
bool
getSolution( const Components& components, int component,
Solution<REAL>& solbuffer ) override
{
const int* colset = components.getComponentsCols( component );
for( std::size_t i = 0; i < n_cols; ++i )
solbuffer.primal[colset[i]] = variables[i]->solution_value();
return true;
}
SolverType
getType() override
{
return SolverType::MIP;
}
String
getName() override
{
return "or-tools";
}
void
printDetails() override
{
}
void
addParameters( ParameterSet& paramSet ) override
{
paramSet.addParameter( "ortools.solver_id", "solver id", solver_id );
}
~GlopInterface() = default;
};
template <typename REAL>
class GlopFactory : public SolverFactory<REAL>
{
GlopFactory() {
solver = new GlopInterface<REAL>();
}
GlopInterface<REAL>* solver;
public:
virtual std::unique_ptr<SolverInterface<REAL>>
newSolver( VerbosityLevel verbosity ) const
{
auto glop = std::unique_ptr<SolverInterface<REAL>>( solver );
return std::move( glop );
}
virtual void
add_parameters( ParameterSet& parameter ) const
{
solver->addParameters( parameter );
}
static std::unique_ptr<SolverFactory<REAL>>
create()
{
return std::unique_ptr<SolverFactory<REAL>>( new GlopFactory<REAL>() );
}
};
}
#endif