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
  /*
#ifdef USE_scip
  example<ampls::SCIPModel>();
#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
#ifdef USE_scip
  example<ampls::SCIPModel>();
#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.
const double DESIREDGAP = 0.7;

class MyGenericCallback : public ampls::GenericCallback
{
public:
  double lastBound = std::numeric_limits<double>::infinity();

  int  evaluateObjs(bool mipsol) {
    double objBound = getValue(ampls::Value::MIP_OBJBOUND).dbl;
    double gap = getValue(ampls::Value::MIP_RELATIVEGAP).dbl;
    if (!mipsol) // print only if bound has improved
    {
      if (lastBound == objBound)
        return 0;
    }
    lastBound = objBound;
    printf("Objective: %f6.3 - Bound: %6.3f - RelGap: %6.3f%%\n",
      getObj(), objBound, 100*gap);
    if (gap < DESIREDGAP)
    {
      printf("Desired gap reached, terminating");
      return -1;
    }
    return 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:
    //  if (!canDo(ampls::CanDo::GET_MIP_SOLUTION)) return 0; // For some solvers (notably SCIP)
      return evaluateObjs(false);
    case ampls::Where::MIPSOL:
     // if (!canDo(ampls::CanDo::GET_MIP_SOLUTION)) return 0; // For some solvers (notably SCIP)
      double mipgap = getValue(ampls::Value::MIP_RELATIVEGAP).dbl;
      return evaluateObjs(true);
    }
    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);
  try {
    m.setOption("threads", 1);
  }
  catch (...) {}
  m.setAMPLParameter(ampls::SolverParams::DBL_MIPGap, 0.001);
  try {
    m.setOption("outlev", 1);
  //  m.setOption("return_mipgap", 5);
  //  m.setOption("mipstartvalue", 3);
  //  m.setOption("mipstartalg", 2);
  //  m.setOption("mipdisplay", 0);
  }
  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-158*DESIREDGAP) && (obj <= 158 + 158 * DESIREDGAP) );
  ampls::Status::SolStatus s = m.getStatus();
  assert(s == ampls::Status::INTERRUPTED);
  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;
    case ampls::Status::INTERRUPTED:
      printf("Interrupted.\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
  /*
  * COPT does not support stopping the optimization in a callback yet
#ifdef USE_copt
  example<ampls::CoptModel>();
#endif*/

#ifdef USE_scip
  example<ampls::SCIPModel>();
#endif
  
#ifdef USE_gurobi
  example<ampls::GurobiModel>();
#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 <type_traits>
#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);
        }
      }
    }
    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): ";
        for (auto n : st.getNodes())
          std::cout << n << " ";
        std::cout << "\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.. I have now " << getValue(ampls::Value::N_ROWS) << " constraints.\n";

      }
      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);

  #ifdef USE_cplex
  // Implementing parallel MIP callbacks with CPLEX is more involving,
  // see the CPLEX-specific examples:
  // cplex/genericbranch-cplex
  // cplex/genericbranch-ampls
  if (std::is_same<T, ampls::CPLEXModel>::value)
     model.setOption("threads", 1);
  #endif

  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_scip
  example<ampls::SCIPModel>();
#endif
#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);

  #ifdef USE_scip
  if (std::is_same<T, ampls::SCIPModel>::value)
    model.setOption("pre:maxrounds", 0);
  #endif

  model.optimize();
  printStatistics(model, model.driver());
  assert(model.getNumVars() == 2);
  assert(model.getNumCons() == 1);
  assert(model.getObj() == 4.0);


  // 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());
  assert(model.getNumVars() == 2);
  assert(model.getNumCons() == 2);
  assert(model.getObj() == 2.0);


  // 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());
  assert(model.getNumVars() == 3);
  assert(model.getNumCons() == 2);
  assert(model.getObj() == 1002.0);

  ampls::AMPLAPIInterface::importModel(ampl, model);
  printStatistics(ampl);
}


int main(int argc, char** argv) {
#ifdef USE_scip
  example<ampls::SCIPModel>();
#endif
#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;
 
}