Generic Examples¶
These examples use the generic interface, which allow the same set of functions to be called on all supported solvers.
Example: load¶
This example shows how to load a model, solve it and display basic information.
#include "ampls/ampls.h"
#include "test-config.h" // for MODELS_DIR
#include <cassert>
#include <string>
// This example shows how to:
// 1. load a model from an NL file, passing an option at loading time
// 2. solve it
// 3. get basic information about the solution
// 4. write out a solution file ready to be imported in AMPL
// Note that this function could have included the ampls::AMPLModel creation
// but having the base class as a parameter helps intellisense
void doStuff(ampls::AMPLModel& m) {
// Set parameter with common mapping
m.setAMPLParameter(ampls::SolverParams::DBL_MIPGap, 0.1);
// Start the optimization process
m.optimize();
printf("\n%s solution ", m.driver());
// Get the generic status
ampls::Status::SolStatus s = m.getStatus();
assert(s == ampls::Status::OPTIMAL);
switch (s)
{
case ampls::Status::OPTIMAL:
printf("optimal\n");
break;
case ampls::Status::INFEASIBLE:
printf("infeasible\n");
break;
case ampls::Status::UNBOUNDED:
printf("unbounded\n");
break;
default:
printf("Status (%d)\n", s);
}
// Get the objective value
double obj = m.getObj();
printf("%s: Objective=%f\n", m.driver(), obj);
// Get the MIP gap using its generic alias
printf("%s: Relative MIP gap=%f\n", m.driver(), m.getAMPLDoubleAttribute(ampls::SolverAttributes::DBL_RelMIPGap));
// Get the solution vector and count the non-zeros
auto solution = m.getSolutionVector();
int nnz = 0;
for (int i = 0; i < solution.size(); i++)
if (solution[i] != 0) nnz++;
printf("%s: Number of non zeroes=%d\n", m.driver(), nnz);
assert(nnz==51);
// Set an option of the solver driver (in this case, instructing it to
// return the MIP gap when communicating the solution back to AMPL)
m.setOption("return_mipgap", 3);
// Write a solution file ready to be imported in AMPL, that includes the MIP gap
char BUFFER[100];
sprintf(BUFFER, "%s-%s.sol", m.getFileName().c_str(), m.driver());
printf("%s: Writing solution file to: %s\n\n\n", m.driver(), BUFFER);
m.writeSol(BUFFER);
}
template <class T> void example() {
assert(1);
const char* MODELNAME = "tsp.nl";
std::string md(MODELS_DIR);
md += MODELNAME;
const char* options[2] = { "outlev=1", NULL };
T model = ampls::AMPLModel::load<T>(md.c_str(), options);
doStuff(model);
}
int main(int argc, char** argv) {
#ifdef USE_gurobi
example<ampls::GurobiModel>();
#endif
#ifdef USE_xpress
example<ampls::XPRESSModel>();
#endif
#ifdef USE_cplex
example< ampls::CPLEXModel>();
#endif
#ifdef USE_copt
example<ampls::CoptModel>();
#endif
#ifdef USE_cbcmp
example<ampls::CbcModel>();
#endif
return 0;
}
Example: driver error detection¶
This example shows how to catch an error coming from the underlying driver implementation; such errors typically arise when specifying a wrong option when loading the model.
#include <vector>
#include <string>
#include <iostream>
#include <cstring>
#include <exception>
#include "ampls/ampls.h"
#include "test-config.h" // for MODELS_DIR
template <class T> void example()
{
const char* MODELNAME = "tspg96.nl";
std::string md(MODELS_DIR);
md += MODELNAME;
try {
T m=ampls::AMPLModel::load<T>(md.c_str());
m.setOption("nonvalid", 1);
m.optimize();
double obj = m.getObj();
printf("\nObjective (%s)=%f\n", m.driver(), obj);
}
catch (const ampls::AMPLSolverException& e)
{
// Catch error message
std::cout << e.what() << std::endl;
}
}
int main(int argc, char** argv) {
#ifdef USE_xpress
example<ampls::XPRESSModel>();
#endif
#ifdef USE_cplex
example<ampls::CPLEXModel>();
#endif
#ifdef USE_gurobi
example<ampls::GurobiModel>();
#endif
#ifdef USE_cbcmp
example<ampls::CbcModel>();
#endif
#ifdef USE_copt
example<ampls::CoptModel>();
#endif
return 0;
}
Example: information callback¶
This example shows how to monitor the progress of the solution process by registering a callback.
#include "ampls/ampls.h"
#include "test-config.h" // for MODELS_DIR
#include <string>
#include <cassert>
// This example illustrates how to obtain basic information
// during the solution process using generic callbacks.
class MyGenericCallback : public ampls::GenericCallback
{
int nMIPnodes = 0;
int nadd = 0;
int run()
{
// Prints out the name of the solution phase where the solver is called from
// (solver specific)
//printf("\nCalled from %s\n", getWhere());
double obj;
// Get the generic mapping
ampls::Where::CBWhere where = getAMPLWhere();
//printf("Where: %i\n", where);
switch (where)
{
case ampls::Where::MSG:
// Print out all messages coming from the solver
std::cout<<getMessage();
return 0;
case ampls::Where::PRESOLVE:
if((getValue(ampls::Value::PRE_DELROWS).integer+
getValue(ampls::Value::PRE_DELCOLS).integer+
getValue(ampls::Value::PRE_COEFFCHANGED).integer) > 0)
printf("\nRemoved %i rows and %i columns. %i coefficients changed",
getValue(ampls::Value::PRE_DELROWS).integer,
getValue(ampls::Value::PRE_DELCOLS).integer,
getValue(ampls::Value::PRE_COEFFCHANGED).integer);
return 0;
case ampls::Where::MIPNODE:
nMIPnodes++;
return 0;
case ampls::Where::MIPSOL:
try {
obj = getObj();
printf("\nMIP Objective = %f", getObj());
printf("\nRel MIP GAP: %f", getValue(ampls::Value::MIP_RELATIVEGAP).dbl);
return 0;
}
catch (...) {
return 0;
}
}
return 0;
}
};
template<class T> void example()
{
const char* MODELNAME = "queens18.nl";
std::string md(MODELS_DIR);
md += MODELNAME;
T m = ampls::AMPLModel::load<T>(md.c_str());
// Set a (generic) callback
MyGenericCallback cb;
m.setCallback(&cb);
m.setAMPLParameter(ampls::SolverParams::DBL_MIPGap, 0.001);
try {
m.setOption("return_mipgap", 5);
m.setOption("mipstartvalue", 3);
m.setOption("mipstartalg", 2);
m.setOption("mipdisplay", 2);
}
catch (const std::exception& e) {
printf(e.what());
}
// Start the optimization process
m.optimize();
// Get the objective value
double obj = m.getObj();
printf("\nSolution with %s=%f\n", m.driver(), obj);
assert( (obj>= 158-10e-6) && (obj <= 158 + 10e-6) );
ampls::Status::SolStatus s = m.getStatus();
assert(s == ampls::Status::OPTIMAL);
switch (s)
{
case ampls::Status::OPTIMAL:
printf("Optimal.\n");
break;
case ampls::Status::INFEASIBLE:
printf("Infeasible.\n");
break;
case ampls::Status::UNBOUNDED:
printf("Unbounded.\n");
break;
default:
printf("Status: %d\n", s);
}
// Write the AMPL sol file
m.writeSol();
}
int main(int argc, char** argv) {
#ifdef USE_cplex
example<ampls::CPLEXModel >();
#endif
#ifdef USE_gurobi
example<ampls::GurobiModel>();
#endif
#ifdef USE_copt
example<ampls::CoptModel>();
#endif
#ifdef USE_xpress
example<ampls::XPRESSModel>();
#endif
#ifdef USE_cbcmp
example<ampls::CbcModel>();
#endif
return 0;
}
Example: add cuts in callback¶
This example solves a traveling salesman problem by adding lazy cut to exclude subtours:
#include <list>
#include <algorithm>
#include <map>
#include <set>
#include <regex>
#include <iostream>
#include <exception>
#include <tuple>
#include <fstream>
#include "ampls/ampls.h"
#include "test-config.h" // for MODELS_DIR
#include "ampl/ampl.h"
const double INTTOLERANCE = 1e-4;
struct GraphArc {
int from;
int to;
GraphArc() : from(0), to(0) {}
GraphArc(int from, int to)
{
this->from = from;
this->to = to;
}
static GraphArc fromVar(const std::string& var)
{
#if __cplusplus >= 201103L
auto s = split(var, "[A-Za-z0-9]*\\['([A-Za-z0-9]*)','([A-Za-z0-9]*)']");
#else
auto s = split(var, "[A-Za-z0-9]*\\[[']*([A-Za-z0-9]*)[']*,[']*([A-Za-z0-9]*)[']*\\]");
#endif
return GraphArc(std::stoi(s[1]), std::stoi(s[2]));
}
friend std::ostream& operator<<(std::ostream& out, const GraphArc& c)
{
out << '(' << c.from << "->" << c.to << ')';
return out;
}
/**
* Find arc connected to the passed one.
* If it finds one, it changes the arc parameter to the next arc
* in the sequence and returns true. Returns false otherwise.
*/
static bool findArcFromTo(const std::list<GraphArc>& arcs, GraphArc& arc)
{
for (auto other : arcs) // Find next arc in the sequence
// since they are not directed, any "from" can be connected
// to any "to"
{
if ((other.from == arc.to) || (other.to == arc.to) ||
(other.from == arc.from) || (other.to == arc.from))
{
arc = other;
return true;
}
}
return false;
}
static std::list<GraphArc> solutionToArcs(const std::vector<double>& sol,
std::map<int, GraphArc>& xvar, std::size_t startVar)
{
std::list<GraphArc> res;
for (int i = startVar; i < sol.size(); i++)
{
if (sol[i] > INTTOLERANCE)
res.push_back(xvar[i]);
}
return res;
}
static std::vector<std::string> split(const std::string str, const std::string regex_str)
{
const std::regex regex(regex_str);
std::smatch matches;
std::vector<std::string> ret;
if (std::regex_search(str, matches, regex)) {
for (size_t i = 0; i < matches.size(); ++i) {
ret.push_back(matches[i].str());
}
}
else
throw std::runtime_error("Match not found");
return ret;
}
};
inline bool operator<(const GraphArc& c1, const GraphArc& c2)
{
if (c1.from < c2.from)
return true;
if (c1.from > c2.from)
return false;
return (c1.to < c2.to);
}
inline bool operator== (const GraphArc& c1, const GraphArc& c2)
{
return (c1.from == c2.from) && (c1.to == c2.to);
}
class Tour
{
std::vector<GraphArc> arcs;
public:
void add(const GraphArc& arc)
{
arcs.push_back(arc);
}
/**
* Get all different nodes tourched by this tour
*/
std::set<int> getNodes() const {
std::set<int> vertices;
for (auto a : arcs)
{
vertices.insert(a.from);
vertices.insert(a.to);
}
return vertices;
}
std::size_t numNodes() const {
return getNodes().size();
}
friend std::ostream& operator<<(std::ostream& out, const Tour& t)
{
std::set<GraphArc> toVisit(++t.arcs.begin(), t.arcs.end());
GraphArc current = *t.arcs.begin();
out << current.from << "-" << current.to;
int lastTo = current.to;
bool found = false;
while (toVisit.size() > 0)
{
found = false;
for (auto a : toVisit)
{
if ((a.from == lastTo) || (a.to == lastTo))
{
toVisit.erase(a);
current = a;
if (a.from == lastTo)
{
out << "-" << a.to;
lastTo = a.to;
}
else
{
out << "-" << a.from;
lastTo = a.from;
}
found = true;
break;
}
}
if (!found)
{
current = *toVisit.begin();
toVisit.erase(current);
lastTo = current.to;
out << "MALFORMED: (" << current.from << "-" << current.to << ") ";
}
}
return out;
}
static std::vector<Tour> findSubtours(std::list<GraphArc> arcs)
{
std::vector<Tour> subTours;
while (arcs.size() > 0)
{
Tour t;
GraphArc start = arcs.front();
t.add(start);
arcs.remove(start);
while (GraphArc::findArcFromTo(arcs, start))
{
t.add(start);
arcs.remove(start);
}
subTours.push_back(t);
}
return subTours;
}
};
class TSPUtil {
private:
static void declare(ampl::AMPL& a) {
a.eval("set NODES ordered; param hpos{ NODES }; param vpos{NODES};");
a.eval("set PAIRS := {i in NODES, j in NODES : ord(i) < ord(j)};");
a.eval("param distance{ (i,j) in PAIRS }:= sqrt((hpos[j] - hpos[i]) **2 + (vpos[j] - vpos[i]) **2);");
a.eval("var X{ PAIRS } binary;");
a.eval("var Length;");
a.eval("minimize Tour_Length : Length;");
a.eval("subject to Visit_All{i in NODES } : sum{ (i, j) in PAIRS } X[i, j] + sum{ (j, i) in PAIRS } X[j, i] = 2;");
a.eval("c: Length = sum{ (i,j) in PAIRS } distance[i, j] * X[i, j];");
}
static std::string trim(const std::string& str)
{
size_t first = str.find_first_not_of(' ');
if (std::string::npos == first)
{
return str;
}
size_t last = str.find_last_not_of(' ');
return str.substr(first, (last - first + 1));
}
static std::tuple<std::vector<double>,
std::vector<double>, std::vector<double>> readTSP(const char* path) {
std::ifstream infile(path);
std::string line;
int nitems = 0;
bool coords = false;
int node;
double x, y;
std::vector<double> nodes;
std::vector<double> xs, ys;
while (std::getline(infile, line))
{
if (nitems == 0) {
if (line.find("DIMENSION") != std::string::npos)
{
auto ss = trim(line.substr(line.find(":") + 1));
std::istringstream iss(ss);
if (!(iss >> nitems))
break; // errror
}
continue;
}
if (!coords)
{
if (line.find("NODE_COORD_SECTION") != std::string::npos)
{
coords = true;
continue;
}
}
if (coords) {
std::istringstream iss(trim(line));
if (iss >> node >> x >> y) {
nodes.push_back(node);
xs.push_back(x);
ys.push_back(y);
}
}
}
for (auto x : nodes)
{
printf("node: %d\n", x);
}
return std::make_tuple(nodes, xs, ys);
}
public:
static void createModel(ampl::AMPL& a, const char* tspFile) {
declare(a);
auto t = readTSP(tspFile);//
std::vector<double> nodes;
std::vector<double> xs, ys;
std::tie(nodes, xs, ys) = t;
auto df = ampl::DataFrame(1, { "NODES", "hpos", "vpos" });
df.setColumn("NODES", nodes.data(), nodes.size());
df.setColumn("hpos", xs.data(), nodes.size());
df.setColumn("vpos", ys.data(), nodes.size());
a.setData(df, "NODES");
a.eval("display NODES, hpos, vpos;");
}
};
std::set<int> setDiff(std::set<int> s1, std::set<int> s2)
{
std::set<int> result;
std::set_difference(s1.begin(), s1.end(), s2.begin(), s2.end(),
std::inserter(result, result.end()));
return result;
}
class MyGenericCallbackCut : public ampls::GenericCallback
{
std::size_t minXIndex;
std::map<int, GraphArc> xvar;
std::map<GraphArc, int> xinverse;
std::set<int> vertices;
int nrun;
public:
MyGenericCallbackCut() : minXIndex(0), nrun(0) {
}
void setMap(std::map<std::string, int> map)
{
minXIndex = map.size();
for (auto it : map)
{
GraphArc cur = GraphArc::fromVar(it.first);
xvar[it.second] = cur;
xinverse[cur] = it.second;
if (it.second < minXIndex)
minXIndex = it.second;
vertices.insert(cur.from);
vertices.insert(cur.to);
}
}
virtual int run()
{
if (getAMPLWhere() == ampls::Where::MSG)
std::cout << getMessage() << "\n";
// Get the generic mapping
if ((getAMPLWhere() == ampls::Where::MIPSOL) &&
canDo(ampls::CanDo::ADD_LAZY_CONSTRAINT) )
{
std::cout << "Bound=" << getValue(ampls::Value::MIP_OBJBOUND) << "\n";
std::cout << "Obj="<< getValue(ampls::Value::OBJ) << "\n";
nrun++;
// Add the the cut!
auto arcs = GraphArc::solutionToArcs(getSolutionVector(), xvar, minXIndex);
auto sts = Tour::findSubtours(arcs);
std::cout << "Iteration " << nrun << ": Found " << sts.size() << " subtours. \n";
int i = 0;
for(auto st:sts)
std::cout << "Subtour " << i++ << ": (" << st.numNodes() << " nodes)\n";
if (sts.size() > 1)
{
for (auto st : sts)
{
auto stn = st.getNodes();
auto nstn = setDiff(vertices, stn);
std::vector<GraphArc> toAdd;
for (int inside : stn) {
for (int outside : nstn)
{
if (inside < outside)
toAdd.push_back(GraphArc(inside, outside));
else
toAdd.push_back(GraphArc(outside, inside));
}
}
std::vector<int> varsIndexToAdd;
for (GraphArc a : toAdd)
varsIndexToAdd.push_back(xinverse[a]);
std::sort(varsIndexToAdd.begin(), varsIndexToAdd.end());
std::vector<double> coeffs(varsIndexToAdd.size());
std::fill(coeffs.begin(), coeffs.end(), 1);
ampls::Constraint c= addLazyIndices(varsIndexToAdd.size(),
varsIndexToAdd.data(), coeffs.data(),
ampls::CutDirection::GE, 2);
}
std::cout << "Added cuts..";
}
std::cout << "Continue solving.\n";
return 0;
}
return 0;
}
};
double doStuff(ampls::AMPLModel& m)
{
// Set a (generic) callback
MyGenericCallbackCut cb;
cb.setMap(m.getVarMapFiltered("X"));
m.enableLazyConstraints();
m.setCallback(&cb);
// Start the optimization process
m.optimize();
// Get the objective value
double obj = m.getObj();
printf("\nObjective with callback (%s)=%f\n", m.driver(), obj);
// Get the solution vector
std::map<int, GraphArc> xvar;
std::size_t minXIndex = INT_MAX;
for (auto it : m.getVarMapFiltered("X"))
{
GraphArc cur = GraphArc::fromVar(it.first);
xvar[it.second] = cur;
if (it.second < minXIndex)
minXIndex = it.second;
}
auto a = GraphArc::solutionToArcs(m.getSolutionVector(), xvar, minXIndex);
auto sts = Tour::findSubtours(a);
// Print solution
std::cout << "Solution has " << sts.size() << " subtours\n";
int i = 0;
for (auto st : sts)
std::cout << "SUBTOUR " << i++ << " (" << st.numNodes() << " nodes): " << st << "\n";
return obj;
}
template <class T> std::pair<std::string, double> example() {
char buffer[255];
strcpy(buffer, MODELS_DIR);
strcat(buffer, "tsp/gr96.tsp");
ampl::AMPL a;
try {
// Create a TSP model instance using AMPL
TSPUtil::createModel(a, buffer);
}
catch (const std::exception& e) {
std::cout << e.what();
}
// Used to store the results
std::map<std::string, double> res;
double obj;
auto model = ampls::AMPLAPIInterface::exportModel<T>(a);
obj = doStuff(model);
return { model.driver(), obj };
}
int main(int argc, char** argv) {
// Used to store the results
std::map<std::string, double> res;
double obj;
#ifdef USE_gurobi
res.insert(example<ampls::GurobiModel>());
#endif
/*
#ifdef USE_xpress
res.insert(example<ampls::XPRESSModel>());
#endif
*/
#ifdef USE_cplex
res.insert(example< ampls::CPLEXModel>());
#endif
#ifdef USE_copt
res.insert(example<ampls::CoptModel>());
#endif
#ifdef USE_cbcmp
res.insert(example<ampls::CbcModel>());
#endif
// Print out the results
for (auto a : res) {
printf("%s=%f\n", a.first.c_str(), a.second);
}
return 0;
}
Example (requires AMPLAPI): add variables and constraints¶
This example shows how to export a model from AMPLAPI to ampls::, how to modify it at solver interface level by adding a variable and a constraint, optimize it and export the model back to AMPLAPI, including the newly defined entities.
#include <iostream>
#include "ampls/ampls.h"
#include "test-config.h" // for MODELS_DIR
#include "ampl/ampl.h"
void printStatistics(ampl::AMPL& ampl) {
printf("AMPL: I have %d variables and %d constraints\n",
static_cast<int>(ampl.getValue("_nvars").dbl()),
static_cast<int>(ampl.getValue("_ncons").dbl()));
printf("My variables are\n");
ampl.eval("display _VARS;");
}
void printStatistics(ampls::AMPLModel& m, const char* name)
{
printf("\n%s: I have %d variables and %d constraints\n",
name, m.getNumVars(), m.getNumCons());
printf("%s: objective=%f\n", name, m.getObj());
}
template <class T> void example()
{
ampl::AMPL ampl;
ampl.eval("var x >=0 integer; var y>=0 integer;");
ampl.eval("maximize z: x+y + to_come;");
ampl.eval("constraint: x+2*y+ to_come<=4;");
printStatistics(ampl);
T model = ampls::AMPLAPIInterface::exportModel<T>(ampl);
model.optimize();
printStatistics(model, model.driver());
// Create a new constraint using the solver interface
int n = model.getNumVars();
std::vector<int> indices(n);
std::vector<double> coeff(n, 1);
for (int i = 0; i < n; i++)
indices[i] = i;
// Add it to the solver and records it for AMPL
model.record(model.addConstraint(n, indices.data(), coeff.data(), ampls::CutDirection::LE, n));
model.optimize();
printStatistics(model, model.driver());
// Add a variable that does not appear in the constraints matrix
// but with a coefficient of 100 in the objective
model.record(model.addVariable(0, NULL, NULL, 0, 10, 100, ampls::VarType::Integer));
model.optimize();
printStatistics(model, model.driver());
ampls::AMPLAPIInterface::importModel(ampl, model);
printStatistics(ampl);
}
int main(int argc, char** argv) {
#ifdef USE_copt
example<ampls::CoptModel>();
#endif
#ifdef USE_cplexmp
example<ampls::CPLEXModel>();
#endif
#ifdef USE_gurobi
example<ampls::GurobiModel>();
#endif
#ifdef USE_xpress
example<ampls::XPRESSModel>();
#endif
#ifdef USE_cbcmp
example<ampls::CbcModel>();
#endif
return 0;
}