16 List of examples

List of examples shipped in the distribution of Fusion API for C++:

Table 16.1 List of distributed examples

File

Description

TrafficNetworkModel.cc

Demonstrates a traffic network problem as a conic quadratic problem (CQO)

alan.cc

A portfolio choice model alan.gms from the GAMS model library

baker.cc

A small bakery revenue maximization linear problem

breaksolver.cc

Shows how to break a long-running task

callback.cc

An example of data/progress callback

ceo1.cc

A simple conic exponential problem

cqo1.cc

A simple conic quadratic problem

diet.cc

Solving Stigler’s Nutrition model diet from the GAMS model library

duality.cc

Shows how to access the dual solution

elastic.cc

Linear regression with elastic net. Demonstrates model parametrization.

facility_location.cc

Demonstrates a small one-facility location problem (CQO)

gp1.cc

A simple geometric program (GP) in conic form

lo1.cc

A simple linear problem

logistic.cc

Implements logistic regression and simple log-sum-exp (CEO)

lownerjohn_ellipsoid.cc

Computes the Lowner-John inner and outer ellipsoidal approximations of a polytope (SDO, Power Cone)

lpt.cc

Demonstrates how to solve the multi-processor scheduling problem and input an integer feasible point (MIP)

mico1.cc

A simple mixed-integer conic problem

milo1.cc

A simple mixed-integer linear problem

mioinitsol.cc

A simple mixed-integer linear problem with an initial guess

modelLib.cc

Library of implementations of basic functions

nearestcorr.cc

Solves the nearest correlation matrix problem (SDO, CQO)

opt_server_sync.cc

Uses MOSEK OptServer to solve an optimization problem synchronously

parameters.cc

Shows how to set optimizer parameters and read information items

portfolio_1_basic.cc

Portfolio optimization - basic Markowitz model

portfolio_2_frontier.cc

Portfolio optimization - efficient frontier

portfolio_3_impact.cc

Portfolio optimization - market impact costs

portfolio_4_transcost.cc

Portfolio optimization - transaction costs

portfolio_5_card.cc

Portfolio optimization - cardinality constraints

pow1.cc

A simple power cone problem

primal_svm.cc

Implements a simple soft-margin Support Vector Machine (CQO)

qcqp_sdo_relaxation.cc

Demonstrate how to use SDP to solve convex relaxation of a mixed-integer QCQO problem

reoptimization.cc

Demonstrate how to modify and re-optimize a linear problem

response.cc

Demonstrates proper response handling

sdo1.cc

A simple semidefinite problem with one matrix variable and a quadratic cone

sdo2.cc

A simple semidefinite problem with two matrix variables

sdo3.cc

A simple semidefinite problem with many matrix variables of the same dimension

sospoly.cc

Models the cone of nonnegative polynomials and nonnegative trigonometric polynomials using Nesterov’s framework

sudoku.cc

A SUDOKU solver (MIP)

total_variation.cc

Demonstrates how to solve a total variation problem (CQO)

tsp.cc

Solves a simple Travelling Salesman Problem and shows how to add constraints to a model and re-optimize (MIP)

Additional examples can be found on the MOSEK website and in other MOSEK publications.

TrafficNetworkModel.cc

Listing 16.1 TrafficNetworkModel.cc Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      TrafficNetworkModel.cc
//
// Purpose:   Demonstrates a traffic network problem as a conic quadratic problem.
//
// Source:    Robert Fourer, "Convexity Checking in Large-Scale Optimization",
//            OR 53 --- Nottingham 6-8 September 2011.
//
// The problem:
//            Given a directed graph representing a traffic network
//            with one source and one sink, we have for each arc an
//            associated capacity, base travel time and a
//            sensitivity. Travel time along a specific arc increases
//            as the flow approaches the capacity.
//
//            Given a fixed inflow we now wish to find the
//            configuration that minimizes the average travel time.

#include <string>
#include <iostream>
#include <iomanip>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

struct TrafficNetworkError : std::runtime_error
{
  TrafficNetworkError(const char * msg) : std::runtime_error(msg) {}
};

std::shared_ptr<ndarray<double, 1>>
                                 TrafficNetworkModel
                                 ( int                                numberOfNodes,
                                   int                                source_idx,
                                   int                                sink_idx,
                                   std::shared_ptr<ndarray<int, 1>>    arc_i,
                                   std::shared_ptr<ndarray<int, 1>>    arc_j,
                                   std::shared_ptr<ndarray<double, 1>> arcSensitivity,
                                   std::shared_ptr<ndarray<double, 1>> arcCapacity,
                                   std::shared_ptr<ndarray<double, 1>> arcBaseTravelTime,
                                   double                             T)
{
  Model::t M = new Model("Traffic Network"); auto _M = finally([&]() { M->dispose(); });
  M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; } );

  int n = numberOfNodes;
  int narcs = arc_j->size(0);

  Matrix::t sens =
    Matrix::sparse(n, n, arc_i, arc_j, arcSensitivity);
  Matrix::t cap =
    Matrix::sparse(n, n, arc_i, arc_j, arcCapacity);
  Matrix::t basetime =
    Matrix::sparse(n, n, arc_i, arc_j, arcBaseTravelTime);
  Matrix::t e =
    Matrix::sparse(n, n, arc_i, arc_j, 1.0);
  Matrix::t e_e =
    Matrix::sparse(n, n,
                   new_array_ptr<int, 1>({ sink_idx }),
                   new_array_ptr<int, 1>({ source_idx }),
                   1.0);
  std::shared_ptr<ndarray<double, 1>>
  cs_inv(new ndarray<double, 1>(narcs, std::function<double(ptrdiff_t)>([&](ptrdiff_t i) { return 1.0 / ((*arcSensitivity)[i] * (*arcCapacity)[i]); }))),
  s_inv (new ndarray<double, 1>(narcs, std::function<double(ptrdiff_t)>([&](ptrdiff_t i) { return 1.0 / (*arcSensitivity)[i]; })));

  Matrix::t cs_inv_matrix = Matrix::sparse(n, n, arc_i, arc_j, cs_inv);
  Matrix::t s_inv_matrix  = Matrix::sparse(n, n, arc_i, arc_j, s_inv);

  Variable::t flow       = M->variable("traffic_flow", Set::make(n,n), Domain::greaterThan(0.0));

  Variable::t x = flow;
  Variable::t t = M->variable("travel_time", Set::make(n,n), Domain::greaterThan(0.0));
  Variable::t d = M->variable("d",           Set::make(n,n), Domain::greaterThan(0.0));
  Variable::t z = M->variable("z",           Set::make(n,n), Domain::greaterThan(0.0));

  // Set the objective:
  M->objective("Average travel time",
               ObjectiveSense::Minimize,
               Expr::mul(1.0 / T, Expr::add(Expr::dot(basetime, x), Expr::dot(e, d))));

  // Set up constraints
  // Constraint (1a)
  M->constraint("(1a)",
                Var::hstack(d->pick(arc_i, arc_j),
                            z->pick(arc_i, arc_j),
                            x->pick(arc_i, arc_j)),
                Domain::inRotatedQCone(narcs, 3));

  // Constraint (1b)
  M->constraint("(1b)",
                Expr::sub(Expr::add(Expr::mulElm(z, e),
                                    Expr::mulElm(x, cs_inv_matrix)),
                          s_inv_matrix),
                Domain::equalsTo(0.0));

  // Constraint (2)
  M->constraint("(2)",
                Expr::sub(Expr::add(Expr::mulDiag(x, e->transpose()),
                                    Expr::mulDiag(x, e_e->transpose())),
                          Expr::add(Expr::mulDiag(x->transpose(), e),
                                    Expr::mulDiag(x->transpose(), e_e))),
                Domain::equalsTo(0.0));
  // Constraint (3)
  M->constraint("(3)",
                x->index(sink_idx, source_idx), Domain::equalsTo(T));
  M->solve();
  return flow->level();
}


int main(int argc, char ** argv)
{
  std::shared_ptr<ndarray<int, 1>>    arc_i    = new_array_ptr<int, 1>   ({  0,    0,    2,    1,    2   });
  std::shared_ptr<ndarray<int, 1>>    arc_j    = new_array_ptr<int, 1>   ({  1,    2,    1,    3,    3   });
  std::shared_ptr<ndarray<double, 1>> arc_base = new_array_ptr<double, 1>({  4.0,  1.0,  2.0,  1.0,  6.0 });
  std::shared_ptr<ndarray<double, 1>> arc_cap  = new_array_ptr<double, 1>({ 10.0, 12.0, 20.0, 15.0, 10.0 });
  std::shared_ptr<ndarray<double, 1>> arc_sens = new_array_ptr<double, 1>({  0.1,  0.7,  0.9,  0.5,  0.1 });

  int      n          = 4;
  double   T          = 20.0;
  int      source_idx = 0;
  int      sink_idx   = 3;

  auto flow = TrafficNetworkModel(
                n, source_idx, sink_idx,
                arc_i, arc_j,
                arc_sens,
                arc_cap,
                arc_base,
                T);
  std::cout << "Optimal flow:" << std::endl;
  for (int i = 0; i < arc_i->size(); ++i)
    std::cout << "\tflow " << (*arc_i)[i] << " -> " << (*arc_j)[i] << " = " << (*flow)[(*arc_i)[i] * n + (*arc_j)[i]] << std::endl;
}

alan.cc

Listing 16.2 alan.cc Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:     alan.cc
//
//  Purpose: This file contains an implementation of the alan.gms (as
//  found in the GAMS online model collection) using Fusion.
//
//  The model is a simple portfolio choice model. The objective is to
//  invest in a number of assets such that we minimize the risk, while
//  requiring a certain expected return.
//
//  We operate with 4 assets (hardware,software, show-biz and treasure
//  bill). The risk is defined by the covariance matrix
//    Q = [[  4.0, 3.0, -1.0, 0.0 ],
//         [  3.0, 6.0,  1.0, 0.0 ],
//         [ -1.0, 1.0, 10.0, 0.0 ],
//         [  0.0, 0.0,  0.0, 0.0 ]]
//
//
//  We use the form Q = U^T * U, where U is a Cholesky factor of Q.
//

#include <iostream>
#include "monty.h"
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

// Security names
std::string securities[] =          { "hardware", "software", "show-biz", "t-bills" };
// Two examples of mean returns on securities
auto mean1 = new_array_ptr<double, 1>({       8.0,        9.0,       12.0,       7.0 });
auto mean2 = new_array_ptr<double, 1>({       9.0,        7.0,       11.0,       5.0 });  
// Target mean return
double target = 10.0;

int numsec = 4;

// Factor of covariance matrix.
auto U = new_array_ptr<double, 2>(
{ { 2.0       ,  1.5       , -0.5       , 0.0 },
  { 0.0       ,  1.93649167,  0.90369611, 0.0 },
  { 0.0       ,  0.0       ,  2.98886824, 0.0 },
  { 0.0       ,  0.0       ,  0.0       , 0.0 }
});

// Solve an instance with given expected return
void solve(Model::t M, 
           Parameter::t mean, 
           Variable::t x, 
           std::shared_ptr<ndarray<double,1>> meanVal) {
  
  std::cout << "Solve with mean = " << std::endl;
  for (int i = 0; i < numsec; ++i)
    std::cout << "  " << securities[i] << " : " << (*meanVal)[i] << std::endl;

  mean->setValue(meanVal);
  M->solve();

  auto solx = x->level();

  std::cout << "Solution = " << std::endl;
  for (int i = 0; i < numsec; ++i)
    std::cout << "  " << securities[i] << " : " << (*solx)[i] << std::endl;
}


int main(int argc, char ** argv)
{
  // Create a parametrized model
  Model::t M = new Model("alan"); auto _M = finally([&]() { M->dispose(); });
  //M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; } );

  Variable::t x = M->variable("x", numsec, Domain::greaterThan(0.0));
  Variable::t t = M->variable("t", 1,      Domain::greaterThan(0.0));
  
  M->objective("minvar", ObjectiveSense::Minimize, t);

  // sum securities to 1.0
  M->constraint("wealth",  Expr::sum(x), Domain::equalsTo(1.0));
  // define target expected return
  Parameter::t mean = M->parameter(numsec);
  M->constraint("dmean", Expr::dot(mean, x), Domain::greaterThan(target));

  M->constraint("t > ||Ux||^2", Expr::vstack(0.5, t, Expr::mul(U, x)), Domain::inRotatedQCone());

  // Solve two instances of the problem
  solve(M, mean, x, mean1);
  solve(M, mean, x, mean2);

  return 0;
}

baker.cc

Listing 16.3 baker.cc Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      baker.cc
//
// Purpose: Demonstrates a small linear problem.
//
// Source: "Linaer Algebra" by Knut Sydsaeter and Bernt Oeksendal.
//
// The problem: A baker has 150 kg flour, 22 kg sugar, 25 kg butter and two
// recipes:
//   1) Cakes, requiring 3.0 kg flour, 1.0 kg sugar and 1.2 kg butter per dozen.
//   2) Breads, requiring 5.0 kg flour, 0.5 kg sugar and 0.5 kg butter per dozen.
// Let the revenue per dozen cakes be $4 and the revenue per dozen breads be $6.
//
// We now wish to compute the combination of cakes and breads that will optimize
// the total revenue.


#include <iostream>
#include <string>
#include "monty.h"
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  auto ingredientnames = new_array_ptr<std::string, 1>({ "Flour", "Sugar", "Butter" });
  auto productnames = new_array_ptr<std::string, 1>({ "Cakes", "Breads" });
  auto stock = new_array_ptr<double, 1>(
  { 150.0,   22.0,    25.0 });
  auto recipe_data = new_array_ptr<double, 2>(
  { { 3.0, 5.0 },
    { 1.0, 0.5 },
    { 1.2, 0.5 }
  });

  auto revenue = new_array_ptr<double, 1>({ 4.0, 6.0 });
  Matrix::t recipe = Matrix::dense(recipe_data);
  Model::t M = new Model("Recipe"); auto _M = finally([&]() { M->dispose(); });
  // "production" defines the amount of each product to bake.
  Variable::t production = M->variable("production",
                                       Set::make(productnames),
                                       Domain::greaterThan(0.0));
  // The objective is to maximize the total revenue.
  M->objective("revenue",
               ObjectiveSense::Maximize,
               Expr::dot(revenue, production));
  // The prodoction is constrained by stock:
  M->constraint(Expr::mul(recipe, production), Domain::lessThan(stock));
  M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; });

  // We solve and fetch the solution:
  M->solve();
  auto res = production->level();
  std::cout << "Solution:" << std::endl;
  for (int i = 0; i < res->size(); ++i)
    std::cout << "  Number of " << (*productnames)[i] << " : " << (*res)[i] << std::endl;
  std::cout << " Revenue : $" << ((*res)[0] * (*revenue)[0] + (*res)[1] * (*revenue)[1]) << std::endl;
}

breaksolver.cc

Listing 16.4 breaksolver.cc Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      breaksolver.cc
//
//  Purpose: Show how to break a long-running task.
//
//  Requires a parameter defining a timeout in seconds.
//

#include <fusion.h>
#include <iostream>
#include <random>
#include <thread>
#include <stdlib.h>
#include <time.h>
#include <functional>
#include <algorithm>
#include <vector>
#include <stdlib.h>


#include "assert.h"

using namespace mosek::fusion;
using namespace monty;


template<class T>
static std::ostream & operator<<(std::ostream & strm, const std::vector<T> & arr)
{
  strm << "[";
  if (arr.size() > 0)
    strm << arr[0];
  for (auto iter = arr.begin() + 1; iter != arr.end(); ++iter)
    strm << " ," << *iter;

  strm << "]";
  return strm;
}

int main(int argc, char ** argv)
{
  int timeout = 5;

  int n = 200;   // number of binary variables
  int m = n / 3; // number of constraints
  int p = n / 5; // Each constraint picks p variables and requires that exactly half of them are 1

  std::cout << "Build problem...\n";


  std::vector<int> idxs(n); for (int i = 0; i < n; ++i) idxs[i] = i;
  std::shared_ptr< ndarray<int> > cidxs(new ndarray<int>(shape(p)));
  //auto rand = std::bind(std::uniform_int_distribution<int>(0,n-1), std::mt19937(0));

  Model::t M = new Model("SolveBinary"); auto _M = finally([&]() { M->dispose(); } );
  M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; } );
  Variable::t x = M->variable("x", n, Domain::binary());

  for (int i = 0; i < m; ++i)
  {
    std::random_shuffle(idxs.begin(), idxs.end(), [ = ](int i) { return std::rand() % i; });
    std::copy(idxs.begin(), idxs.begin() + p, cidxs->begin());
    M->constraint(Expr::sum(x->pick(cidxs)), Domain::equalsTo(p / 2));
  }
  M->objective(ObjectiveSense::Minimize, Expr::sum(x));

  std::cout << "Start thread...\n";
  bool alive = true;
  std::thread T(std::function<void(void)>([&]() { M->solve(); alive = false; }) );

  time_t T0 = time(NULL);
  while (true)
  {
    if (time(NULL) - T0 > timeout)
    {
      std::cout << "Solver terminated due to timeout!\n";
      M->breakSolver();
      T.join();
      break;
    }
    if (! alive)
    {
      std::cout << "Solver terminated before anything happened!\n";
      T.join();
      break;
    }
  }
  return 0;
}

callback.cc

Listing 16.5 callback.cc Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      callback.cc

   Purpose:   To demonstrate how to use the progress
              callback.

              Compile and link the file with MOSEK, then
              use as follows:

              callback psim
              callback dsim
              callback intpnt

              The argument tells which optimizer to use
              i.e. psim is primal simplex, dsim is dual simplex
              and intpnt is interior-point.
 */


#include <memory>
#include <random>
#include <iostream>
#include <string>

#include "mosek.h"
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;


/* Note: This function is declared using MSKAPI,
         so the correct calling convention is
         employed. */

// For demonstration. Actually defined in mosektask.h as datacbhandler_t and cbhandler_t
typedef std::function<bool(MSKcallbackcodee, const double *, const int32_t *, const int64_t *)> callbackHandler_t;
typedef std::function<int(MSKcallbackcodee)> progressHandler_t;

static int MSKAPI usercallback( MSKcallbackcodee caller,
                                const double   * douinf,
                                const int32_t  * intinf,
                                const int64_t  * lintinf,
                                Model::t mod,
                                const double maxtime)
{
  switch ( caller )
  {
    case MSK_CALLBACK_BEGIN_INTPNT:
      std::cerr << "Starting interior-point optimizer\n";
      break;
    case MSK_CALLBACK_INTPNT:
      std::cerr << "Iterations: " << intinf[MSK_IINF_INTPNT_ITER];
      std::cerr << " (" << douinf[MSK_DINF_OPTIMIZER_TIME] << "/";
      std::cerr << douinf[MSK_DINF_INTPNT_TIME] << ")s. \n";
      std::cerr << "Primal obj.: " << douinf[MSK_DINF_INTPNT_PRIMAL_OBJ];
      std::cerr << " Dual obj.: " <<  douinf[MSK_DINF_INTPNT_DUAL_OBJ] << std::endl;
      break;
    case MSK_CALLBACK_END_INTPNT:
      std::cerr << "Interior-point optimizer finished.\n";
      break;
    case MSK_CALLBACK_BEGIN_PRIMAL_SIMPLEX:
      std::cerr << "Primal simplex optimizer started.\n";
      break;
    case MSK_CALLBACK_UPDATE_PRIMAL_SIMPLEX:
      std::cerr << "Iterations: " << intinf[MSK_IINF_SIM_PRIMAL_ITER];
      std::cerr << "  Elapsed time: " << douinf[MSK_DINF_OPTIMIZER_TIME];
      std::cerr << "(" << douinf[MSK_DINF_SIM_TIME] << ")\n";
      std::cerr << "Obj.: " << douinf[MSK_DINF_SIM_OBJ] << std::endl;
      break;
    case MSK_CALLBACK_END_PRIMAL_SIMPLEX:
      std::cerr << "Primal simplex optimizer finished.\n";
      break;
    case MSK_CALLBACK_BEGIN_DUAL_SIMPLEX:
      std::cerr << "Dual simplex optimizer started.\n";
      break;
    case MSK_CALLBACK_UPDATE_DUAL_SIMPLEX:
      std::cerr << "Iterations: " << intinf[MSK_IINF_SIM_DUAL_ITER];
      std::cerr << "  Elapsed time: " << douinf[MSK_DINF_OPTIMIZER_TIME];
      std::cerr << "(" << douinf[MSK_DINF_SIM_TIME] << ")\n";
      std::cerr << "Obj.: " << douinf[MSK_DINF_SIM_OBJ] << std::endl;
      break;
    case MSK_CALLBACK_END_DUAL_SIMPLEX:
      std::cerr << "Dual simplex optimizer finished.\n";
      break;
    case MSK_CALLBACK_BEGIN_BI:
      std::cerr << "Basis identification started.\n";
      break;
    case MSK_CALLBACK_END_BI:
      std::cerr << "Basis identification finished.\n";
      break;
    default:
      break;
  }
  if ( douinf[MSK_DINF_OPTIMIZER_TIME] >= maxtime )
  {
    std::cerr << "MOSEK is spending too much time.Terminate it.\n";
    return 1;
  }
  return 0;
} /* usercallback */

int main(int argc, char ** argv)
{
  std::string slvr("intpnt");

  if (argc <= 1)
  {
    std::cerr << "Usage: ( psim | dsim | intpnt ) \n";
  }

  if (argc >= 2) slvr = argv[1];

  /* Solve a big random linear optimization problem */
  int n = 150, m = 700;
  std::default_random_engine generator;
  std::uniform_real_distribution<double> unif_distr(0., 10.);
  auto c = new_array_ptr<double, 1>(n);
  auto b = new_array_ptr<double, 1>(m);
  auto A = new_array_ptr<double, 1>(m * n);
  std::generate(c->begin(), c->end(), std::bind(unif_distr, generator));
  std::generate(b->begin(), b->end(), std::bind(unif_distr, generator));
  std::generate(A->begin(), A->end(), std::bind(unif_distr, generator));

  Model::t M = new Model();
  auto x = M->variable(n, Domain::unbounded());
  M->constraint(Expr::mul(Matrix::dense(m, n, A), x), Domain::lessThan(b));
  M->objective(ObjectiveSense::Maximize, Expr::dot(c, x));

  if ( slvr == "psim")
    M->setSolverParam("optimizer", "primalSimplex");
  else if ( slvr == "dsim")
    M->setSolverParam("optimizer", "dualSimplex");
  else if ( slvr == "intpnt")
    M->setSolverParam("optimizer", "intpnt");

  double maxtime = 0.07;

  callbackHandler_t cllbck = [&](MSKcallbackcodee caller,
                                 const double * douinf, const int32_t* intinf, const int64_t* lintinf)
  {
    return usercallback(caller, douinf, intinf, lintinf, M, maxtime);
  };

  M->setDataCallbackHandler(cllbck);

  M->setSolverParam("log", 0);
  M->solve();

  return 0;
}

ceo1.cc

Listing 16.6 ceo1.cc Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      ceo1.cc
//
//   Purpose: Demonstrates how to solve the problem
//
//   minimize x1 + x2
//   such that
//            x1 + x2 + x3  = 1.0
//                x1,x2    >= 0.0
//   and      x1 >= x2 * exp(x3/x2)
//
#include <iostream>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  Model::t M = new Model("ceo1"); auto _M = finally([&]() { M->dispose(); });

  Variable::t x  = M->variable("x", 3, Domain::unbounded());

  // Create the constraint
  //      x[0] + x[1] + x[2] = 1.0
  M->constraint("lc", Expr::sum(x), Domain::equalsTo(1.0));

  // Create the exponential conic constraint
  Constraint::t expc = M->constraint("expc", x, Domain::inPExpCone());

  // Set the objective function to (x[0] + x[1])
  M->objective("obj", ObjectiveSense::Minimize, Expr::sum(x->slice(0,2)));

  // Solve the problem
  M->solve();

  // Get the linear solution values
  ndarray<double, 1> xlvl   = *(x->level());
  // Get conic solution of expc1
  ndarray<double, 1> expclvl = *(expc->level());
  ndarray<double, 1> expcdl  = *(expc->dual());

  std::cout << "x1,x2,x3 = " << xlvl << std::endl;
  std::cout << "expc levels = " << expclvl << std::endl;
  std::cout << "expc dual conic var levels = " << expcdl << std::endl;
}

cqo1.cc

Listing 16.7 cqo1.cc Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      cqo1.cc
//
//   Purpose: Demonstrates how to solve the problem
//
//   minimize y1 + y2 + y3
//   such that
//            x1 + x2 + 2.0 x3  = 1.0
//                    x1,x2,x3 >= 0.0
//   and
//            (y1,x1,x2) in C_3,
//            (y2,y3,x3) in K_3
//
//   where C_3 and K_3 are respectively the quadratic and
//   rotated quadratic cone of size 3 defined as
//       C_3 = { z1,z2,z3 :      z1 >= sqrt(z2^2 + z3^2) }
//       K_3 = { z1,z2,z3 : 2 z1 z2 >= z3^2              }
//
#include <iostream>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  Model::t M = new Model("cqo1"); auto _M = finally([&]() { M->dispose(); });

  Variable::t x  = M->variable("x", 3, Domain::greaterThan(0.0));
  Variable::t y  = M->variable("y", 3, Domain::unbounded());

  // Create the aliases
  //      z1 = [ y[0],x[0],x[1] ]
  //  and z2 = [ y[1],y[2],x[2] ]
  Variable::t z1 = Var::vstack(y->index(0),  x->slice(0, 2));
  Variable::t z2 = Var::vstack(y->slice(1, 3), x->index(2));

  // Create the constraint
  //      x[0] + x[1] + 2.0 x[2] = 1.0
  auto aval = new_array_ptr<double, 1>({1.0, 1.0, 2.0});
  M->constraint("lc", Expr::dot(aval, x), Domain::equalsTo(1.0));

  // Create the constraints
  //      z1 belongs to C_3
  //      z2 belongs to K_3
  // where C_3 and K_3 are respectively the quadratic and
  // rotated quadratic cone of size 3, i.e.
  //                 z1[0] >= sqrt(z1[1]^2 + z1[2]^2)
  //  and  2.0 z2[0] z2[1] >= z2[2]^2
  Constraint::t qc1 = M->constraint("qc1", z1, Domain::inQCone());
  Constraint::t qc2 = M->constraint("qc2", z2, Domain::inRotatedQCone());

  // Set the objective function to (y[0] + y[1] + y[2])
  M->objective("obj", ObjectiveSense::Minimize, Expr::sum(y));

  // Solve the problem
  M->solve();

  // Get the linear solution values
  ndarray<double, 1> xlvl   = *(x->level());
  ndarray<double, 1> ylvl   = *(y->level());
  // Get conic solution of qc1
  ndarray<double, 1> qc1lvl = *(qc1->level());
  ndarray<double, 1> qc1dl  = *(qc1->dual());

  std::cout << "x1,x2,x2 = " << xlvl << std::endl;
  std::cout << "y1,y2,y3 = " << ylvl << std::endl;
  std::cout << "qc1 levels = " << qc1lvl << std::endl;
  std::cout << "qc1 dual conic var levels = " << qc1dl << std::endl;
}

diet.cc

Listing 16.8 diet.cc Click here to download.
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      diet.cc
//
// Purpose: Solving Stigler's Nutrition model (DIET,SEQ=7)
//
// Source: GAMS Model library,
//           Dantzig, G B, Chapter 27.1. In Linear Programming and Extensions.
//           Princeton University Press, Princeton, New Jersey, 1963.
//
// Given a set of nutrients, foods with associated nutrient values, allowance of
// nutrients per day, the model will find the cheapest combination of foods
// which will provide the necessary nutrients.
//
//  Arguments for construction:
//  name      - Model name.
//  foods     - List of M names of the foods.
//  nutrients - List of N names of the nutrients.
//  daily_allowance - List of N floats denoting the daily allowance of each
//              nutrient.
//  nutritive_value - Two-dimensional MxN array of floats where each row
//              denotes the nutrient values for a single food per $ spent.

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
#include "fusion.h"
#include <cassert>

using namespace mosek::fusion;
using namespace monty;


struct DataException : std::runtime_error
{
  DataException(const char * msg) : std::runtime_error(msg) {}
};


typedef std::pair< std::shared_ptr<ndarray<double, 1>>, std::shared_ptr<ndarray<double, 1>> > diet_return_t;

diet_return_t diet
( std::shared_ptr<ndarray<std::string, 1>> foods,
  std::shared_ptr<ndarray<std::string, 1>> nutrients,
  std::shared_ptr<ndarray<double, 1>>      dailyAllowance,
  std::shared_ptr<ndarray<double, 2>>      nutritive_value )
{
  Model::t M = new Model("Stinglers Diet Problem"); auto _M = finally([&]() { M->dispose(); });

  Matrix::t nutrientValue = Matrix::t(Matrix::dense(nutritive_value))->transpose();

  int m = foods->size();
  int n = nutrients->size();

  if (dailyAllowance->size() != n)
    throw DataException("Length of dailyAllowance does not match the number of nutrients");
  if (nutrientValue->numColumns() != m)
    throw DataException("Number of rows in nutrient_value does not match the number of foods");
  if (nutrientValue->numRows() != n)
    throw DataException("Number of columns in nutrient_value does not match the number of nutrients");

  
  Variable::t dailyPurchase =
    M->variable(Set::make(foods),
                Domain::greaterThan(0.0));
  
  Constraint::t dailyNutrients =
    M->constraint(Expr::mul(nutrientValue, dailyPurchase),
                  Domain::greaterThan(dailyAllowance));
  M->objective(ObjectiveSense::Minimize, Expr::sum(dailyPurchase));
  M->solve();

  return diet_return_t(dailyPurchase->level(), dailyNutrients->level());
}

/* Main class with data definitions */
int main(int argc, char ** argv)
{
  auto nutrient_unit = new_array_ptr<std::string, 1>({
    "thousands",  "grams",        "grams",
    "milligrams", "thousand ius", "milligrams",
    "milligrams", "milligrams",   "milligrams"
  });
  auto nutrients = new_array_ptr<std::string, 1>({
    "calorie",    "protein",      "calcium",
    "iron",       "vitamin a",    "vitamin b1",
    "vitamin b2", "niacin",       "vitamin c"
  });
  auto foods = new_array_ptr<std::string, 1>({
    "wheat",        "cornmeal", "cannedmilk", "margarine", "cheese",
    "peanut butter", "lard",     "liver",      "porkroast", "salmon",
    "greenbeans",   "cabbage",  "onions",     "potatoes",  "spinach",
    "sweet potatos", "peaches",  "prunes",     "limabeans", "navybeans"
  });

  auto nutritive_value = new_array_ptr<double, 2>(
  {
    //  calorie       calcium      vitamin a        vitamin b2      vitamin c
    //         protein        iron           vitamin b1      niacin
    {44.7,  1411,   2.0,   365,    0,      55.4,   33.3,  441,     0},  // wheat
    {36,     897,   1.7,    99,   30.9,    17.4,    7.9,  106,     0},  // cornmeal
    { 8.4,   422,  15.1,     9,   26,       3,     23.5,   11,    60},  // cannedmilk
    {20.6,    17,    .6,     6,   55.8,      .2,    0,      0,     0},  // margarine
    { 7.4,   448,  16.4,    19,   28.1,      .8,   10.3,    4,     0},  // cheese
    {15.7,   661,   1,      48,     0,      9.6,    8.1,  471,     0},  // peanut butter
    {41.7,     0,   0,       0,     .2,     0,       .5,    5,     0},  // lard
    { 2.2,   333,    .2,   139,  169.2,     6.4,   50.8,  316,   525},  // liver
    { 4.4,   249,    .3,    37,    0,      18.2,    3.6,   79,     0},  // porkroast
    { 5.8,   705,   6.8,    45,    3.5,     1,      4.9,  209,     0},  // salmon
    { 2.4,   138,   3.7,    80,   69,       4.3,    5.8,   37,   862},  // greenbeans
    { 2.6,   125,   4,      36,    7.2,     9,      4.5,   26,  5369},  // cabbage
    { 5.8,   166,   3.8,    59,   16.6,     4.7,    5.9,   21,  1184},  // onions
    {14.3,   336,   1.8,   118,    6.7,    29.4,    7.1,  198,  2522},  // potatoes
    { 1.1,   106,   0.0,   138,  918.4,     5.7,   13.8,   33,  2755},  // spinach
    { 9.6,   138,   2.7,    54,  290.7,     8.4,    5.4,   83,  1912},  // sweet potatos
    { 8.5,    87,   1.7,   173,   86.8,     1.2,    4.3,   55,    57},  // peaches
    {12.8,    99,   2.5,   154,   85.7,     3.9,    4.3,   65,   257},  // prunes
    {17.4,  1055,   3.7,   459,    5.1,    26.9,   38.2,   93,     0},  // limabeans
    {26.9,  1691,  11.4,   792,    0,      38.4,   24.6,  217,     0}
  });// navybeans

  auto daily_allowance = new_array_ptr<double, 1>(
  {   3.,     70.,   .8,    12.,   5.,      1.8,    2.7,   18.,   75. });

  auto res = diet(foods,
                  nutrients,
                  daily_allowance,
                  nutritive_value);

  std::cout << "Solution:" << std::endl;
  for (int i = 0; i < foods->size(); ++i)
    std::cout << (*foods)[i] << ": " << std::setprecision(2) << std::setw(8) << (*res.first)[i] << std::endl;

  std::cout << "Nutrients:" << std::endl;
  for (int i = 0; i < nutrients->size(); ++i)
    std::cout << "\t" << std::setw(15) << (*nutrients)[i] << ": " << std::setprecision(2) << std::setw(8) << (*res.second)[i]
              << "(" << (*daily_allowance)[i] << ")"
              << std::endl;
}

duality.cc

Listing 16.9 duality.cc Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      duality.cc
//
//   Purpose: A small example of how to access dual values.

#include <memory>
#include <iostream>
#include <iomanip>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  auto A = new_array_ptr<double, 2>({ { -0.5, 1.0 }  });
  auto b = new_array_ptr<double, 1>({ 1.0 });
  auto c = new_array_ptr<double, 1>({ 1.0, 1.0 });

  Model::t M = new Model("duality"); auto _M = finally([&]() { M->dispose(); });

  Variable::t x = M->variable("x", 2, Domain::greaterThan(0.0));

  Constraint::t con = M->constraint(Expr::sub(Expr::mul(A, x), b), Domain::equalsTo(0.0));

  M->objective("obj", ObjectiveSense::Minimize, Expr::dot(c, x));

  M->solve();
  auto xsol = x->level();
  auto ssol = x->dual();
  auto ysol = con->dual();

  std::cout << std::setiosflags(std::ios::scientific)
            << "x1 = " << (*xsol)[0] << std::endl
            << "x2 = " << (*xsol)[1] << std::endl
            << "s1 = " << (*ssol)[0] << std::endl
            << "s1 = " << (*ssol)[1] << std::endl
            << "y =  " << (*ysol)[0] << std::endl;
}

elastic.cc

Listing 16.10 elastic.cc Click here to download.
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      elastic.cc
//
// Purpose: Demonstrates model parametrization on the example of an elastic net linear regression:
//
//          min_x  |Ax-b|_2 + lambda1*|x|_1 + lambda2*|x|_2

#include <iostream>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

// Construct the model with parameters b, lambda1, lambda2
// and with prescribed matrix A
Model::t initializeModel(int m, int n, std::shared_ptr<ndarray<double,2>> A) {
  Model::t M = new Model(); 
  auto x = M->variable("x", n);

  // t >= |Ax-b|_2 where b is a parameter
  auto b = M->parameter("b", m);
  auto t = M->variable();
  M->constraint(Expr::vstack(t, Expr::sub(Expr::mul(A, x), b)), Domain::inQCone());

  // p_i >= |x_i|, i=1..n
  auto p = M->variable(n);
  M->constraint(Expr::hstack(p, x), Domain::inQCone());

  // q >= |x|_2
  auto q = M->variable();
  M->constraint(Expr::vstack(q, x), Domain::inQCone());

  // Objective, parametrized with lambda1, lambda2
  // t + lambda1*sum(p) + lambda2*q
  auto lambda1 = M->parameter("lambda1");
  auto lambda2 = M->parameter("lambda2");
  auto obj = Expr::add(new_array_ptr<Expression::t, 1>({t, Expr::mul(lambda1, Expr::sum(p)), Expr::mul(lambda2, q)}));
  M->objective(ObjectiveSense::Minimize, obj);

  // Return the ready model
  return M;
}

int smallExample() {
  //Create a small example
  int m = 4;
  int n = 2;
  auto A = new_array_ptr<double, 2>(
                { {1.0,   2.0},
                  {3.0,   4.0},
                  {-2.0, -1.0},
                  {-4.0, -3.0} });

  auto M = initializeModel(m, n, A);

  // For convenience retrieve some elements of the model
  auto b = M->getParameter("b");
  auto lambda1 = M->getParameter("lambda1");
  auto lambda2 = M->getParameter("lambda2");
  auto x = M->getVariable("x");

  // First solve
  b->setValue(new_array_ptr<double, 1>({0.1, 1.2, -1.1, 3.0}));
  lambda1->setValue(0.1);
  lambda2->setValue(0.01);

  M->solve();
  auto sol = x->level();
  std::cout << "Objective " << M->primalObjValue() << ", solution " << (*sol)[0] << ", " << (*sol)[1] << "\n";

  // Increase lambda1
  lambda1->setValue(0.5);
  
  M->solve();
  sol = x->level();
  std::cout << "Objective " << M->primalObjValue() << ", solution " << (*sol)[0] << ", " << (*sol)[1] << "\n";

  // Now change the data completely
  b->setValue(new_array_ptr<double, 1>({1.0, 1.0, 1.0, 1.0}));
  lambda1->setValue(0.0);
  lambda2->setValue(0.0);
  
  M->solve();
  sol = x->level();
  std::cout << "Objective " << M->primalObjValue() << ", solution " << (*sol)[0] << ", " << (*sol)[1] << "\n";

  // And increase lamda2
  lambda2->setValue(1.4145);
  
  M->solve();
  sol = x->level();
  std::cout << "Objective " << M->primalObjValue() << ", solution " << (*sol)[0] << ", " << (*sol)[1] << "\n";

  M->dispose();
  return 0;
}

int main() {
  return smallExample();
}

facility_location.cc

Listing 16.11 facility_location.cc Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      facility_location.cc
//
//  Purpose: Demonstrates a small one-facility location problem.
//
//  Given 10 customers placed in a grid we wish to place a facility
//  somewhere so that the total sum of distances to customers is
//  minimized.
//
//  The problem is formulated as a conic optimization problem as follows.
//  Let f=(fx,fy) be the (unknown) location of the facility, and let
//  c_i=(cx_i,cy_i) be the (known) customer locations; then we wish to
//  minimize
//      sum_i || f - c_i ||
//  where
//      ||.||
//  denotes the euclidian norm.
//  This is formulated as
//
//  minimize   sum(d_i)
//  such that  d_i ^ 2 > tx_i ^ 2 + ty_i ^ 2, for all i
//             tx_i = cx_i - fx, for all i
//             ty_i = cy_i - fy, for all i
//             d_i > 0, for all i
////

#include <memory>
#include <iostream>
#include <iomanip>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  auto customerloc = new_array_ptr<double, 2>(
  { {12.,  2. },
    {15., 13. },
    {10.,  8. },
    { 0., 10. },
    { 6., 13. },
    { 5.,  8. },
    {10., 12. },
    { 4.,  6. },
    { 5.,  2. },
    { 1., 10. }
  } );
  int N = customerloc->size(0);
  Model::t M = new Model("FacilityLocation"); auto _M = finally([&]() { M->dispose(); });
  // Variable holding the facility location
  Variable::t f = M->variable("facility", Set::make(1, 2), Domain::unbounded());
  // Variable defining the euclidian distances to each customer
  Variable::t d = M->variable("dist", Set::make(N, 1), Domain::greaterThan(0.0));
  // Variable defining the x and y differences to each customer;
  Variable::t t = M->variable("t", Set::make(N, 2), Domain::unbounded());
  M->constraint("dist measure",
                Var::hstack(d, t),
                Domain::inQCone(N, 3));

  Variable::t fxy = Var::repeat(f, N);
  M->constraint("xy diff", Expr::add(t, fxy), Domain::equalsTo(customerloc));


  M->objective("total_dist", ObjectiveSense::Minimize, Expr::sum(d));

  M->solve();

  auto res = f->level();
  std::cout << std::setprecision(2)
            << "Facility location = " << (*res)[0] << "," << (*res)[1] << std::endl;
}

gp1.cc

Listing 16.12 gp1.cc Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      gp1.cc
//
//   Purpose:   Demonstrates how to solve a simple Geometric Program (GP)
//              cast into conic form with exponential cones and log-sum-exp.
//
//              Example from
//                https://gpkit.readthedocs.io/en/latest/examples.html//maximizing-the-volume-of-a-box
//
#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

// Models log(sum(exp(Ax+b))) <= 0.
// Each row of [A b] describes one of the exp-terms
void logsumexp(Model::t                             M, 
               std::shared_ptr<ndarray<double, 2>>  A, 
               Variable::t                          x,
               std::shared_ptr<ndarray<double, 1>>  b)
{
  int k = A->size(0);
  auto u = M->variable(k);
  M->constraint(Expr::sum(u), Domain::equalsTo(1.0));
  M->constraint(Expr::hstack(u,
                             Expr::constTerm(k, 1.0),
                             Expr::add(Expr::mul(A, x), b)), Domain::inPExpCone());
}

// maximize     h*w*d
// subjecto to  2*(h*w + h*d) <= Awall
//              w*d <= Afloor
//              alpha <= h/w <= beta
//              gamma <= d/w <= delta
//
// Variable substitutions:  h = exp(x), w = exp(y), d = exp(z).
//
// maximize     x+y+z
// subject      log( exp(x+y+log(2/Awall)) + exp(x+z+log(2/Awall)) ) <= 0
//                              y+z <= log(Afloor)
//              log( alpha ) <= x-y <= log( beta )
//              log( gamma ) <= z-y <= log( delta )
std::shared_ptr<ndarray<double, 1>> max_volume_box(double Aw, double Af, 
                                                   double alpha, double beta, double gamma, double delta)
{
  Model::t M = new Model("max_vol_box"); auto _M = finally([&]() { M->dispose(); });

  auto xyz = M->variable(3);
  M->objective("Objective", ObjectiveSense::Maximize, Expr::sum(xyz));
    
  logsumexp(M, 
            new_array_ptr<double,2>({{1,1,0}, {1,0,1}}), 
            xyz, 
            new_array_ptr<double,1>({log(2.0/Aw), log(2.0/Aw)}));
    
  M->constraint(Expr::dot(new_array_ptr<double,1>({0,1,1}), xyz), Domain::lessThan(log(Af)));
  M->constraint(Expr::dot(new_array_ptr<double,1>({1,-1,0}), xyz), Domain::inRange(log(alpha),log(beta)));
  M->constraint(Expr::dot(new_array_ptr<double,1>({0,-1,1}), xyz), Domain::inRange(log(gamma),log(delta)));
    
  M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; } );
  M->solve();
  
  return std::make_shared<ndarray<double, 1>>(shape(3), [xyz](ptrdiff_t i) { return exp((*(xyz->level()))[i]); });
}

int main()
{
  double Aw    = 200.0;
  double Af    = 50.0;
  double alpha = 2.0;
  double beta  = 10.0;
  double gamma = 2.0;
  double delta = 10.0;
  
  auto hwd = max_volume_box(Aw, Af, alpha, beta, gamma, delta);

  std::cout << std::setprecision(4);
  std::cout << "h=" << (*hwd)[0] << " w=" << (*hwd)[1] << " d=" << (*hwd)[2] << std::endl;
}

lo1.cc

Listing 16.13 lo1.cc Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      lo1.cc
//
//  Purpose: Demonstrates how to solve the problem
//
//  maximize 3*x0 + 1*x1 + 5*x2 + x3
//  such that
//           3*x0 + 1*x1 + 2*x2        = 30,
//           2*x0 + 1*x1 + 3*x2 + 1*x3 > 15,
//                  2*x1 +      + 3*x3 < 25
//  and
//           x0,x1,x2,x3 > 0,
//           0 < x1 < 10
////

#include <iostream>
#include "fusion.h"
using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  auto A1 = new_array_ptr<double, 1>({ 3.0, 1.0, 2.0, 0.0 });
  auto A2 = new_array_ptr<double, 1>({ 2.0, 1.0, 3.0, 1.0 });
  auto A3 = new_array_ptr<double, 1>({ 0.0, 2.0, 0.0, 3.0 });
  auto c  = new_array_ptr<double, 1>({ 3.0, 1.0, 5.0, 1.0 });

  // Create a model with the name 'lo1'
  Model::t M = new Model("lo1"); auto _M = finally([&]() { M->dispose(); });

  M->setLogHandler([ = ](const std::string & msg) { std::cout << msg << std::flush; } );

  // Create variable 'x' of length 4
  Variable::t x = M->variable("x", 4, Domain::greaterThan(0.0));

  // Create constraints
  M->constraint(x->index(1), Domain::lessThan(10.0));
  M->constraint("c1", Expr::dot(A1, x), Domain::equalsTo(30.0));
  M->constraint("c2", Expr::dot(A2, x), Domain::greaterThan(15.0));
  M->constraint("c3", Expr::dot(A3, x), Domain::lessThan(25.0));

  // Set the objective function to (c^t * x)
  M->objective("obj", ObjectiveSense::Maximize, Expr::dot(c, x));

  // Solve the problem
  M->solve();

  // Get the solution values
  auto sol = x->level();
  std::cout << "[x0,x1,x2,x3] = " << (*sol) << "\n";
}

logistic.cc

Listing 16.14 logistic.cc Click here to download.
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      logistic.cc
//
// Purpose: Implements logistic regression with regulatization.
//
//          Demonstrates using the exponential cone and log-sum-exp in Fusion.

#include <iostream>
#include <random>
#include <sstream>
#include <vector>
#include <utility>
#include "fusion.h"
using namespace mosek::fusion;
using namespace monty;

// t >= log( 1 + exp(u) ) coordinatewise
void softplus(Model::t      M,
              Expression::t t,
              Expression::t u)
{
    int n = (*t->getShape())[0];
    auto z1 = M->variable(n);
    auto z2 = M->variable(n);
    M->constraint(Expr::add(z1, z2), Domain::equalsTo(1));
    M->constraint(Expr::hstack(z1, Expr::constTerm(n, 1.0), Expr::sub(u,t)), Domain::inPExpCone());
    M->constraint(Expr::hstack(z2, Expr::constTerm(n, 1.0), Expr::neg(t)), Domain::inPExpCone());
}

// Model logistic regression (regularized with full 2-norm of theta)
// X - n x d matrix of data points
// y - length n vector classifying training points
// lamb - regularization parameter
std::pair<Model::t, Variable::t>
logisticRegression(std::vector<std::vector<double>> & X, 
                   std::vector<bool>                & y,
                   double                           lamb)
{
    int n = X.size();
    int d = X[0].size();       // num samples, dimension
    
    Model::t M = new Model();   

    auto theta = M->variable(d);
    auto t     = M->variable(n);
    auto reg   = M->variable();

    M->objective(ObjectiveSense::Minimize, Expr::add(Expr::sum(t), Expr::mul(lamb,reg)));
    M->constraint(Var::vstack(reg, theta), Domain::inQCone());

    auto signs = std::make_shared<ndarray<double,1>>(shape(n), [y](ptrdiff_t i) { return y[i] ? -1 : 1; });   

    softplus(M, t, Expr::mulElm(Expr::mul(new_array_ptr<double>(X), theta), signs));

    return std::make_pair(M, theta);
}


int main(int arc, char** argv)
{

  // Test: detect and approximate a circle using degree 2 polynomials
  std::vector<std::vector<double>> X;
  std::vector<bool> Y;

  for(double x = -1; x <= 1; x += 1.0/15)
    for(double y = -1; y <= 1; y += 1.0/15) 
  {
    std::vector<double> p;
    p.push_back(1.0); p.push_back(x); p.push_back(y); p.push_back(x*y);
    p.push_back(x*x); p.push_back(y*y);
    X.push_back(p);
    Y.push_back(x*x+y*y>=0.69);
  }

  auto res = logisticRegression(X, Y, 0.1);
  auto M = res.first;
  auto theta = res.second;

  M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; });
  M->solve();
  std::cout << *(theta->level());
  M->dispose();

  return 0;
}

lownerjohn_ellipsoid.cc

Listing 16.15 lownerjohn_ellipsoid.cc Click here to download.
/*
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:      lownerjohn_ellipsoid.cc

  Purpose:
  Computes the Lowner-John inner and outer ellipsoidal
  approximations of a polytope.

  References:
    [1] "Lectures on Modern Optimization", Ben-Tal and Nemirovski, 2000.
    [2] "MOSEK modeling manual", 2013
*/

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
#include "fusion.h"
#include <cassert>

using namespace mosek::fusion;
using namespace monty;

/**
Models the convex set

  S = { (x, t) \in R^n x R | x >= 0, t <= (x1 * x2 * ... * xn)^(1/n) }

using three-dimensional power cones
*/
void geometric_mean(Model::t M, Variable::t x, Variable::t t)
{
  int n = (int) x->getSize();
  if (n==1)
      M->constraint(Expr::sub(t, x), Domain::lessThan(0.0));
  else
  {
    Variable::t t2 = M->variable();
    M->constraint(Var::hstack(t2, x->index(n-1), t), Domain::inPPowerCone(1-1.0/n));
    geometric_mean(M, x->slice(0,n-1), t2);
  }
}

/**
 Purpose: Models the hypograph of the n-th power of the
 determinant of a positive definite matrix. See [1,2] for more details.

   The convex set (a hypograph)

   C = { (X, t) \in S^n_+ x R |  t <= det(X)^{1/n} },

   can be modeled as the intersection of a semidefinite cone

   [ X, Z; Z^T Diag(Z) ] >= 0

   and a number of rotated quadratic cones and affine hyperplanes,

   t <= (Z11*Z22*...*Znn)^{1/n}  (see geometric_mean).
*/
Variable::t det_rootn(Model::t M, Variable::t t, int n)
{
  // Setup variables
  Variable::t Y = M->variable(Domain::inPSDCone(2 * n));

  Variable::t X   = Y->slice(new_array_ptr<int,1>({0, 0}), new_array_ptr<int,1>({n, n}));
  Variable::t Z   = Y->slice(new_array_ptr<int,1>({0, n}), new_array_ptr<int,1>({n, 2 * n}));
  Variable::t DZ  = Y->slice(new_array_ptr<int,1>({n, n}), new_array_ptr<int,1>({2 * n, 2 * n}));

  // Z is lower-triangular
  std::shared_ptr<ndarray<int,2>> low_tri( new ndarray<int,2>( shape_t<2>(n*(n-1)/2, 2) )); 
  int k = 0;
  for(int i = 0; i < n; i++)
    for(int j = i+1; j < n; j++)
       (*low_tri)(k,0) = i, (*low_tri)(k,1) = j, ++k;
  M->constraint(Z->pick(low_tri), Domain::equalsTo(0.0));
  // DZ = Diag(Z)
  M->constraint(Expr::sub(DZ, Expr::mulElm(Z, Matrix::eye(n))), Domain::equalsTo(0.0));

  // t^n <= (Z11*Z22*...*Znn)
  geometric_mean(M, DZ->diag(), t);

  // Return an n x n PSD variable which satisfies t <= det(X)^(1/n)
  return X;
}

/**
  The inner ellipsoidal approximation to a polytope

     S = { x \in R^n | Ax < b }.

  maximizes the volume of the inscribed ellipsoid,

     { x | x = C*u + d, || u ||_2 <= 1 }.

  The volume is proportional to det(C)^(1/n), so the
  problem can be solved as

    maximize         t
    subject to       t       <= det(C)^(1/n)
                || C*ai ||_2 <= bi - ai^T * d,  i=1,...,m
                C is PSD

  which is equivalent to a mixed conic quadratic and semidefinite
  programming problem.
*/
std::pair<std::shared_ptr<ndarray<double, 1>>, std::shared_ptr<ndarray<double, 1>>>
    lownerjohn_inner
    ( std::shared_ptr<ndarray<double, 2>> A,
      std::shared_ptr<ndarray<double, 1>> b)
{
  Model::t M = new Model("lownerjohn_inner"); auto _M = finally([&]() { M->dispose(); });
  int m = A->size(0);
  int n = A->size(1);

  // Setup variables
  Variable::t t = M->variable("t", 1, Domain::greaterThan(0.0));
  Variable::t C = det_rootn(M, t, n);
  Variable::t d = M->variable("d", n, Domain::unbounded());

  // quadratic cones
  M->constraint(Expr::hstack(Expr::sub(b, Expr::mul(A, d)), Expr::mul(A, C)),
                Domain::inQCone());

  // Objective: Maximize t
  M->objective(ObjectiveSense::Maximize, t);
  M->solve();

  return std::make_pair(C->level(), d->level());
}

/**
  The outer ellipsoidal approximation to a polytope given
  as the convex hull of a set of points

    S = conv{ x1, x2, ... , xm }

  minimizes the volume of the enclosing ellipsoid,

    { x | || P*x-c ||_2 <= 1 }

  The volume is proportional to det(P)^{-1/n}, so the problem can
  be solved as

    maximize         t
    subject to       t       <= det(P)^(1/n)
                || P*xi - c ||_2 <= 1,  i=1,...,m
                P is PSD.
*/
std::pair<std::shared_ptr<ndarray<double, 1>>, std::shared_ptr<ndarray<double, 1>>>
    lownerjohn_outer(std::shared_ptr<ndarray<double, 2>> x)
{
  Model::t M = new Model("lownerjohn_outer");
  int m = x->size(0);
  int n = x->size(1);

  // Setup variables
  Variable::t t = M->variable("t", 1, Domain::greaterThan(0.0));
  Variable::t P = det_rootn(M, t, n);
  Variable::t c = M->variable("c", n, Domain::unbounded());

  // (1, Px-c) \in Q
  M->constraint(Expr::hstack(
                  Expr::ones(m), Expr::sub(Expr::mul(x, P),
                      Var::reshape(Var::repeat(c, m), new_array_ptr<int, 1>({m, n}))) ),
                Domain::inQCone());

  // Objective: Maximize t
  M->objective(ObjectiveSense::Maximize, t);
  M->solve();

  return std::make_pair(P->level(), c->level());
}

std::ostream & operator<<(std::ostream & os, ndarray<double, 1> & a)
{
  os << "[ ";
  if (a.size() > 0)
  {
    os << a[0];
    for (int i = 1; i < a.size(); ++i)
      os << "," << a[i];
  }
  os << " ]";
  return os;
}

/******************************************************************************************/
int main(int argc, char ** argv)
{
  //Vertices of a pentagon in 2D
  int n = 6;
  std::shared_ptr<ndarray<double, 2>> p =
  new_array_ptr<double, 2>({ {0., 0.}, {1., 3.}, {5.5, 4.5}, {7., 4.}, {7., 1.}, {3., -2.} });

  //The h-representation of that polygon
  std::shared_ptr<ndarray<double, 2>> A(
                                     new ndarray<double, 2>(shape_t<2>(n, 2), std::function<double(const shape_t<2> &)>(
                                         [&](const shape_t<2> & ij)
  { int i = ij[0], j = ij[1];
    if (j == 0)
      return -((*p)(i, 1)) + (*p)((i - 1 + n) % n, 1);
    else
      return  ((*p)(i, 0)) - (*p)((i - 1 + n) % n, 0);
  })));
  std::shared_ptr<ndarray<double, 1>> b(
                                     new ndarray<double, 1>(n, std::function<double(ptrdiff_t)>( [&](ptrdiff_t i)
  { return (*A)(i, 0) * (*p)(i, 0) + (*A)(i, 1) * (*p)(i, 1); } )));

  auto Cd = lownerjohn_inner(A, b);
  auto Pc = lownerjohn_outer(p);

  std::cout << "Inner:" << std::endl;
  std::cout << "  C = " << *Cd.first << std::endl;
  std::cout << "  d = " << *Cd.second << std::endl;
  std::cout << "Outer:" << std::endl;
  std::cout << "  P = " << *Pc.first << std::endl;
  std::cout << "  c = " << *Pc.second << std::endl;
}

lpt.cc

Listing 16.16 lpt.cc Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      lpt.cc
//
// Purpose:  Demonstrates how to solve the multi-processor
//           scheduling problem using the Fusion API.


#include <iostream>
#include <random>
#include <sstream>

#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;


int main(int arc, char** argv)
{
  double lb = 1.0;          //Bounds for the length of a short task
  double ub = 5.;

  int     n = 30;           //Number of tasks
  int     m = 6;            //Number of processors

  double  sh = 0.8;         //The proportion of short tasks
  int     n_short = (int)(sh * n);
  int     n_long = n - n_short;

  auto gen = std::bind(std::uniform_real_distribution<double>(lb, ub), std::mt19937(0));

  auto T = std::shared_ptr<ndarray<double, 1> >(new ndarray<double, 1>(n));
  for (int i = 0; i < n_short; i++) (*T)[i] = gen();
  for (int i = n_short; i < n; i++) (*T)[i] = 20 * gen();
  std::sort(T->begin(), T->end(), std::greater<double>());

  Model::t M = new Model("Multi-processor scheduling"); auto _M = finally([&]() { M->dispose(); });

  Variable::t x = M->variable("x", new_array_ptr<int, 1>({m, n}), Domain::binary());
  Variable::t t = M->variable("t", 1, Domain::unbounded());

  M->constraint( Expr::sum(x, 0), Domain::equalsTo(1.) );
  M->constraint( Expr::sub( Var::repeat(t, m), Expr::mul(x, T) ) , Domain::greaterThan(0.) );

  M->objective( ObjectiveSense::Minimize, t );

  //LPT heuristic
  auto schedule = std::shared_ptr<ndarray<double, 1> >(new ndarray<double, 1>(m, 0.));
  auto init = std::shared_ptr<ndarray<double, 1> >(new ndarray<double, 1>(n * m, 0.));

  for (int i = 0; i < n; i++)
  {
    auto pos = std::distance(schedule->begin(), std::min_element(schedule->begin(), schedule->end()));
    (*schedule)[pos] += (*T)[i];
    (*init)[pos * n + i] = 1;
  }

  //Comment this line to switch off feeding in the initial LPT solution
  x->setLevel(init);

  M->setLogHandler([ = ](const std::string & msg) { std::cout << msg << std::flush; } );

  M->setSolverParam("mioTolRelGap", .01);
  M->solve();

  std::cout << "initial solution: \n";
  for (int i = 0; i < m; i++)
  {
    std::cout << "M " << i << " [";
    for (int y = 0; y < n; y++)
      std::cout << int( (*init)[i * n + y] ) << ", ";
    std::cout << "]\n";
  }

  std::cout << "MOSEK solution:\n";
  for (int i = 0; i < m; i++)
  {
    std::cout << "M " << i << " [";
    for (int y = 0; y < n; y++)
      std::cout << int((*(x->index(i, y)->level()))[0]) << ", ";
    std::cout << "]\n";
  }

  return 0;
}

mico1.cc

Listing 16.17 mico1.cc Click here to download.
//
//    Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//    File:    mico1.cc
//
//    Purpose:  Demonstrates how to solve a small mixed
//              integer conic optimization problem.
//
//              minimize    x^2 + y^2
//              subject to  x >= e^y + 3.8
//                          x, y - integer
//
#include <iostream>
#include <iomanip>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  Model::t M = new Model("mico1"); auto _M = finally([&]() { M->dispose(); });

  Variable::t x = M->variable(Domain::integral(Domain::unbounded()));
  Variable::t y = M->variable(Domain::integral(Domain::unbounded()));
  Variable::t t = M->variable();
  
  M->constraint(Expr::vstack(t, x, y), Domain::inQCone());
  M->constraint(Expr::vstack(Expr::sub(x, 3.8), 1, y), Domain::inPExpCone());

  M->objective(ObjectiveSense::Minimize, t);

  M->solve();

  std::cout << std::setprecision(2)
            << "x = " << (*(x->level()))[0] << std::endl
            << "y = " << (*(y->level()))[0] << std::endl ;

  return 0;
}

milo1.cc

Listing 16.18 milo1.cc Click here to download.
//
//    Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//    File:    milo1.cc
//
//    Purpose:  Demonstrates how to solve a small mixed
//              integer linear optimization problem.
//
#include <iostream>
#include <iomanip>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  auto a1 = new_array_ptr<double, 1>({ 50.0, 31.0 });
  auto a2 = new_array_ptr<double, 1>({ 3.0,  -2.0 });
  auto c  = new_array_ptr<double, 1>({  1.0, 0.64 });

  Model::t M = new Model("milo1"); auto _M = finally([&]() { M->dispose(); });
  Variable::t x = M->variable("x", 2, Domain::integral(Domain::greaterThan(0.0)));

  // Create the constraints
  //      50.0 x[0] + 31.0 x[1] <= 250.0
  //       3.0 x[0] -  2.0 x[1] >= -4.0
  M->constraint("c1", Expr::dot(a1, x), Domain::lessThan(250.0));
  M->constraint("c2", Expr::dot(a2, x), Domain::greaterThan(-4.0));

  // Set max solution time
  M->setSolverParam("mioMaxTime", 60.0);
  // Set max relative gap (to its default value)
  M->setSolverParam("mioTolRelGap", 1e-4);
  // Set max absolute gap (to its default value)
  M->setSolverParam("mioTolAbsGap", 0.0);

  // Set the objective function to (c^T * x)
  M->objective("obj", ObjectiveSense::Maximize, Expr::dot(c, x));

  // Solve the problem
  M->solve();

  // Get the solution values
  auto sol = x->level();
  std::cout << std::setiosflags(std::ios::scientific) << std::setprecision(2)
            << "x1 = " << (*sol)[0] << std::endl
            << "x2 = " << (*sol)[1] << std::endl
            << "MIP rel gap = " << M->getSolverDoubleInfo("mioObjRelGap") << " (" << M->getSolverDoubleInfo("mioObjAbsGap") << ")" << std::endl;
}

mioinitsol.cc

Listing 16.19 mioinitsol.cc Click here to download.
//
//    Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//    File:    mioinitsol.cc
//
//    Purpose:  Demonstrates how to solve a small mixed
//              integer linear optimization problem
//              providing an initial feasible solution.
//
#include <memory>
#include <iostream>


#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{

  auto c  = new_array_ptr<double, 1>({7.0, 10.0, 1.0, 5.0});

  Model::t M = new Model("lo1"); auto _M = finally([&]() { M->dispose(); });
  M->setLogHandler([ = ](const std::string & msg) { std::cout << msg << std::flush; } );

  int n = c->size();

  auto x = M->variable("x", n, Domain::greaterThan(0.0));
  x->slice(0,3)->makeInteger();

  M->constraint(Expr::sum(x), Domain::lessThan(2.5));

  M->objective("obj", ObjectiveSense::Maximize, Expr::dot(c, x));

  // Assign values to integer variables.
  // We only set a slice of x     
  auto init_sol = new_array_ptr<double, 1>({ 1.0, 1.0, 0.0 });
  x->slice(0,3)->setLevel( init_sol );

  M->solve();

  auto ss = M->getPrimalSolutionStatus();
  std::cout << "Solution status: " << ss << std::endl;
  auto sol = x->level();
  std::cout << "x = ";

  for (auto s : *sol)
    std::cout << s << ", ";

  // Was the initial solution used?
  int constr = M->getSolverIntInfo("mioConstructSolution");
  double constrVal = M->getSolverDoubleInfo("mioConstructSolutionObj");
  std::cout << "Initial solution utilization: " << constr << std::endl;
  std::cout << "Initial solution objective: " << constrVal << std::endl;

}

modelLib.cc

Listing 16.20 modelLib.cc Click here to download.
/*
  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File :      modelLib.cc

  Purpose: Library of simple building blocks in Mosek Fusion.
*/
#include <iostream>
#include "monty.h"
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

// Duplicate Variable::ts
// x = y
static void dup(Model::t M, Variable::t x, Variable::t y) {
  M->constraint(Expr::sub(x,y), Domain::equalsTo(0.0));
}

// Absolute value
// t >= |x|, where t, x have the same shape
static void abs(Model::t M, Variable::t t, Variable::t x) {
  M->constraint(Expr::add(t,x), Domain::greaterThan(0.0));
  M->constraint(Expr::sub(t,x), Domain::greaterThan(0.0));
}

// 1-norm
// t >= sum( |x_i| ), x is a vector variable
static void norm1(Model::t M, Variable::t t, Variable::t x) {
  Variable::t u = M->variable(x->getShape(), Domain::unbounded());
  abs(M, u, x);
  M->constraint(Expr::sub(t, Expr::sum(u)), Domain::greaterThan(0.0));
}

// Square
// t >= x^2
static void sq(Model::t M, Variable::t t, Variable::t x) {
  M->constraint(Expr::hstack(0.5, t, x), Domain::inRotatedQCone());
}

// 2-norm
// t >= sqrt(x_1^2 + ->->-> + x_n^2) where x is a vector
static void norm2(Model::t M, Variable::t t, Variable::t x) {
  M->constraint(Expr::vstack(t, x), Domain::inQCone());
}

// Power with exponent > 1
// t >= |x|^p (where p>1)
static void pow(Model::t M, Variable::t t, Variable::t x, double p) {
  M->constraint(Expr::hstack(t, 1, x), Domain::inPPowerCone(1.0/p));
}

// Inverse of power 
// t >= 1/|x|^p, x>0 (where p>0)
static void pow_inv(Model::t M, Variable::t t, Variable::t x, double p) {
  M->constraint(Expr::hstack(t, x, 1), Domain::inPPowerCone(1.0/(1.0+p)));
}

// p-norm, p>1
// t >= \|x\|_p (where p>1), x is a vector variable
static void pnorm(Model::t M, Variable::t t, Variable::t x, double p) {
  int n = (int) x->getSize();
  Variable::t r = M->variable(n);
  M->constraint(Expr::sub(t, Expr::sum(r)), Domain::equalsTo(0.0));
  M->constraint(Expr::hstack(Var::repeat(t,n), r, x), Domain::inPPowerCone(1.0-1.0/p));
}

// Geometric mean
// |t| <= (x_1->->->x_n)^(1/n), x_i>=0, x is a vector variable of length >= 1
static void geo_mean(Model::t M, Variable::t t, Variable::t x) {
  int n = (int) x->getSize();
  if (n==1) {
      abs(M, x, t);
  }
  else {
      Variable::t t2 = M->variable();
      M->constraint(Expr::hstack(t2, x->index(n-1), t), Domain::inPPowerCone(1.0-1.0/n));
      geo_mean(M, t2, x->slice(0,n-1));
  }
}

// Logarithm
// t <= log(x), x>=0
static void log(Model::t M, Variable::t t, Variable::t x) {
  M->constraint(Expr::hstack(x, 1, t), Domain::inPExpCone());
}

// Exponential
// t >= exp(x)
static void exp(Model::t M, Variable::t t, Variable::t x) {
  M->constraint(Expr::hstack(t, 1, x), Domain::inPExpCone());
}

// Entropy
// t >= x * log(x), x>=0
static void ent(Model::t M, Variable::t t, Variable::t x) {
  M->constraint(Expr::hstack(1, x, Expr::neg(t)), Domain::inPExpCone());
}

// Relative entropy
// t >= x * log(x/y), x,y>=0
static void relent(Model::t M, Variable::t t, Variable::t x, Variable::t y) {
  M->constraint(Expr::hstack(y, x, Expr::neg(t)), Domain::inPExpCone());
}

// Log-sum-exp
// log( sum_i(exp(x_i)) ) <= t, where x is a vector
static void logsumexp(Model::t M, Variable::t t, Variable::t x) {
  int n = (int) x->getSize();
  Variable::t u = M->variable(n);
  M->constraint(Expr::hstack(u, Expr::constTerm(n, 1.0), Expr::sub(x, Var::repeat(t, n))), Domain::inPExpCone());
  M->constraint(Expr::sum(u), Domain::lessThan(1.0));
}

// Semicontinuous Variable::t
// x = 0 or a <= x <= b
static void semicontinuous(Model::t M, Variable::t x, double a, double b) {
  Variable::t u = M->variable(x->getShape(), Domain::binary());
  M->constraint(Expr::sub(x, Expr::mul(a, u)), Domain::greaterThan(0.0));
  M->constraint(Expr::sub(x, Expr::mul(b, u)), Domain::lessThan(0.0));
}

// Indicator Variable::t
// x!=0 implies t=1. Assumes that |x|<=1 in advance.
static void indicator(Model::t M, Variable::t t, Variable::t x) {
  M->constraint(t, Domain::inRange(0,1));
  t->makeInteger();
  abs(M, t, x);
}

// Logical OR
// x OR y, where x, y are binary
static void logic_or(Model::t M, Variable::t x, Variable::t y) {
  M->constraint(Expr::add(x, y), Domain::greaterThan(1.0));
}
// x_1 OR ->->-> OR x_n, where x is a binary vector
static void logic_or_vect(Model::t M, Variable::t x) {
  M->constraint(Expr::sum(x), Domain::greaterThan(1.0));
}

// SOS1 (NAND)
// at most one of x_1,->->->,x_n, where x is a binary vector (SOS1 constraint)
static void logic_sos1(Model::t M, Variable::t x) {
  M->constraint(Expr::sum(x), Domain::lessThan(1.0));
}
// NOT(x AND y), where x, y are binary
static void logic_nand(Model::t M, Variable::t x, Variable::t y) {
  M->constraint(Expr::add(x, y), Domain::lessThan(1.0));
}

// Cardinality bound
// At most k of entries in x are nonzero, assuming in advance |x_i|<=1.
static void card(Model::t M, Variable::t x, int k) {
  Variable::t t = M->variable(x->getShape(), Domain::binary());
  abs(M, t, x);
  M->constraint(Expr::sum(t), Domain::lessThan(k));
}

// This is just a syntactic test without much sense
static void testModel() {
  Model::t M = new Model();
  Variable::t x = M->variable();
  Variable::t y = M->variable();
  Variable::t t = M->variable();
  Variable::t p = M->variable(5);
  Variable::t a = M->variable(new_array_ptr<int,1>({10,2}));
  Variable::t b = M->variable(new_array_ptr<int,1>({10,2}));
  Variable::t e = M->variable(Domain::binary());
  Variable::t f = M->variable(Domain::binary());

  log(M, t, x);
  exp(M, t, y);
  ent(M, t, x);
  relent(M, t, x, y);
  logsumexp(M, t, p);
  abs(M, a, b);
  norm1(M, t, a);
  sq(M, t, y);
  norm2(M, t, p);
  pow(M, t, x, 1.5);
  pow_inv(M, t, x, 3.3);
  geo_mean(M, t, p);
  semicontinuous(M, y, 1.1, 2.2);
  indicator(M, e, y);
  logic_or(M, e, f);
  logic_nand(M, e, f);
  card(M, b, 5);
}

// A sample problem using functions from the library
//
// max -sqrt(x^2 + y^2) + log(y) - x^1->5
//  st x >= y + 3
//
static void testExample() {
  Model::t M = new Model();
  Variable::t x = M->variable();
  Variable::t y = M->variable();
  Variable::t t = M->variable(3);

  M->constraint(Expr::sub(x, y), Domain::greaterThan(3.0));
  norm2(M, t->index(0), Var::vstack(x,y));
  log  (M, t->index(1), y);
  pow  (M, t->index(2), x, 1.5);

  M->objective(ObjectiveSense::Maximize, Expr::dot(t, new_array_ptr<double,1>({-1,1,-1})));

  M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; } );
  M->solve();
}
  
// Full test
int main() {
  testModel();
  testExample();
  return 0;
}

nearestcorr.cc

Listing 16.21 nearestcorr.cc Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      nearestcorr.cc

   Purpose:
   Solves the nearest correlation matrix problem

     minimize   || A - X ||_F   s.t.  diag(X) = e, X is PSD

   as the equivalent conic program

     minimize     t

     subject to   (t, vec(A-X)) in Q
                  diag(X) = e
                  X >= 0
*/
#include <iostream>
#include <iomanip>
#include <cmath>

#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

//-------------------------------------

void nearestcorr ( std::shared_ptr<ndarray<double, 2>>  A);

void nearestcorr_nn (
  std::shared_ptr<ndarray<double, 2>>  A,
  const std::vector<double>         & gammas,
  std::vector<double>         & res,
  std::vector<double>         & rank);

Expression::t vec(Expression::t e);

void print_mat(std::ostream & os, const std::shared_ptr<ndarray<double, 1>> & a);

int main(int argc, char ** argv)
{
  int N = 5;
  auto A = new_array_ptr<double, 2>(
  { { 0.0,  0.5,  -0.1,  -0.2,   0.5},
    { 0.5,  1.25, -0.05, -0.1,   0.25},
    { -0.1, -0.05,  0.51,  0.02, -0.05},
    { -0.2, -0.1,   0.02,  0.54, -0.1},
    { 0.5,  0.25, -0.05, -0.1,   1.25}
  });

  nearestcorr(A);

  std::vector<double> gammas(11); for (int i = 0; i < gammas.size(); ++i) gammas[i] = 0.1 * i;
  std::vector<double> res(gammas.size());
  std::vector<double> rank(gammas.size());

  nearestcorr_nn(A,
                 gammas,
                 res,
                 rank);

  std::cout << std::setprecision(2) << std::fixed;
  for (int i = 0; i < gammas.size(); ++i)
    std::cout << "gamma = " << gammas[i] << ", rank = " << rank[i] << ", res = " << res[i] << std::endl;
  std::cout.unsetf(std::ios::floatfield);

  return 0;
}


void nearestcorr( std::shared_ptr<ndarray<double, 2>> A)
{
  int N = A->size(0);

  // Create a model
  Model::t M = new Model("NearestCorrelation"); auto _M = finally([&]() { M->dispose(); });

  // Setting up the variables
  Variable::t X = M->variable("X", Domain::inPSDCone(N));
  Variable::t t = M->variable("t", 1, Domain::unbounded());

  // (t, vec (A-X)) \in Q
  M->constraint( Expr::vstack(t, vec(Expr::sub(A, X))), Domain::inQCone() );

  // diag(X) = e
  M->constraint(X->diag(), Domain::equalsTo(1.0));

  // Objective: Minimize t
  M->objective(ObjectiveSense::Minimize, t);

  // Solve the problem
  M->solve();

  // Get the solution values
  std::cout << "X = \n"; print_mat(std::cout, X->level());
  std::cout << "t = " << *(t->level()->begin()) << std::endl;
}

void nearestcorr_nn(
  std::shared_ptr<ndarray<double, 2>>  A,
  const std::vector<double>           & gammas,
  std::vector<double>                 & res,
  std::vector<double>                 & rank)
{
  int N = A->size(0);

  Model::t M = new Model("NucNorm"); auto M_ = monty::finally([&]() { M->dispose(); });

  // Setup variables
  Variable::t t = M->variable("t", 1, Domain::unbounded());
  Variable::t X = M->variable("X", Domain::inPSDCone(N));
  Variable::t w = M->variable("w", N, Domain::greaterThan(0.0));

  // (t, vec (X + diag(w) - A)) in Q
  Expression::t D = Expr::mulElm( Matrix::eye(N), Var::repeat(w, 1, N) );
  M->constraint( Expr::vstack( t, vec(Expr::sub(Expr::add(X, D), A)) ), Domain::inQCone() );

  // Trace(X)
  auto TrX = Expr::sum(X->diag());

  for (int k = 0; k < gammas.size(); ++k)
  {
    // Objective: Minimize t + gamma*Tr(X)
    M->objective(ObjectiveSense::Minimize, Expr::add(t, Expr::mul(gammas[k], TrX )));
    M->solve();

    // Find the eigenvalues of X and approximate its rank
    auto d = new_array_ptr<double, 1>(N);
    mosek::LinAlg::syeig(MSK_UPLO_LO, N, X->level(), d);
    int rnk = 0; for (int i = 0; i < N; ++i) if ((*d)[i] > 1e-6) ++rnk;

    res[k]  = (*(t->level()))[0];
    rank[k] = rnk;
  }
}


/** Assuming that e is an NxN expression, return the lower triangular part as a vector.
*/
Expression::t vec(Expression::t e)
{

  int N = (*e->getShape())[0];
  int dim = N * (N + 1) / 2;

  auto msubi = new_array_ptr<int, 1>(dim);
  auto msubj = new_array_ptr<int, 1>(dim);
  auto mcof  = new_array_ptr<double, 1>(dim);

  for (int i = 0, k = 0; i < N; ++i)
    for (int j = 0; j < i + 1; ++j, ++k)
    {
      (*msubi)[k] = k;
      (*msubj)[k] = i * N + j;
      (*mcof) [k] = (i == j) ? 1.0 : std::sqrt(2.0);
    }

  Matrix::t S = Matrix::sparse(N * (N + 1) / 2, N * N, msubi, msubj, mcof);
  return Expr::mul(S, Expr::reshape(e, N * N));
}

/* Utility methods */
void print_mat(std::ostream & os, const std::shared_ptr<ndarray<double, 1>> & a)
{
  if ( a->size() == 1)
    os << *(a->begin());
  else
  {
    int n = std::sqrt(a->size());
    int i = 0;
    for (auto aa = a->begin(); aa != a->end(); ++aa, ++i)
    {
      if (i % n == 0)
        os << "[";
      os << " " << *aa;
      if ( (i + 1) % n == 0)
        os << "]\n";
    }
  }
}

opt_server_sync.cc

Listing 16.22 opt_server_sync.cc Click here to download.
//
// Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File :      opt_server_sync.cc
//
//  Purpose :   Demonstrates how to use MOSEK OptServer
//              to solve optimization problem synchronously
//
#include <fusion.h>
#include <stdlib.h>
#include <iostream>

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv) {
    if (argc<2) {
        std::cout << "Missing argument, syntax is:" << std::endl;
        std::cout << "  opt_server_sync serveraddr" << std::endl;
        exit(0);
    }

    std::string addr(argv[1]);

    // Setup a simple test problem
    Model::t M = new Model("testOptServer"); auto _M = finally([&]() { M->dispose(); } );
    Variable::t x = M->variable("x", 3, Domain::greaterThan(0.0));
    M->constraint("lc", Expr::dot(new_array_ptr<double, 1>({1.0, 1.0, 2.0}), x), Domain::equalsTo(1.0));
    M->objective("obj", ObjectiveSense::Minimize, Expr::sum(x));

    // Attach log handler
    M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; } );

    // Set OptServer URL
    M->optserverHost(addr);

    // Solve the problem on the OptServer
    M->solve();

    // Get the solution
    std::cout << "x1,x2,x2 = " << *(x->level()) << std::endl;

    return 0;
}

parameters.cc

Listing 16.23 parameters.cc Click here to download.
/*
  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File :      parameters.cc

  Purpose :   Demonstrates a very simple example about how to set
              parameters and read information items
              with MOSEK Fusion
*/
#include <iostream>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  Model::t M = new Model(""); auto _M = finally([&]() { M->dispose(); });
  std::cout << "Test MOSEK parameter get/set functions\n";

  // Set log level (integer parameter)
  M->setSolverParam("log", 1);
  // Select interior-point optimizer... (parameter with symbolic string values)
  M->setSolverParam("optimizer", "intpnt");
  // ... without basis identification (parameter with symbolic string values)
  M->setSolverParam("intpntBasis", "never");
  // Set relative gap tolerance (double parameter)
  M->setSolverParam("intpntCoTolRelGap", 1.0e-7);

  // The same in a different way
  M->setSolverParam("intpntCoTolRelGap", "1.0e-7");

  // Incorrect value
  try {
      M->setSolverParam("intpntCoTolRelGap", -1);
  }
  catch (mosek::fusion::ParameterError) {
      std::cout << "Wrong parameter value\n";
  }

  // Define and solve an optimization problem here
  // M->solve()
  // After optimization: 

  std::cout << "Get MOSEK information items\n";

  double tm = M->getSolverDoubleInfo("optimizerTime");
  int it = M->getSolverIntInfo("intpntIter");

  std::cout << "Time: " << tm << "\nIterations: " << it << "\n";

  return 0;
}

portfolio_1_basic.cc

Listing 16.24 portfolio_1_basic.cc Click here to download.
/*
  File : portfolio_1_basic.cc

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :
    Implements a basic portfolio optimization model.
*/

#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include "monty.h"
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

static double sum(std::shared_ptr<ndarray<double, 1>> x)
{
  double r = 0.0;
  for (auto v : *x) r += v;
  return r;
}

static double dot(std::shared_ptr<ndarray<double, 1>> x,
                  std::shared_ptr<ndarray<double, 1>> y)
{
  double r = 0.0;
  for (int i = 0; i < x->size(); ++i) r += (*x)[i] * (*y)[i];
  return r;
}

static double dot(std::shared_ptr<ndarray<double, 1>> x,
                  std::vector<double> & y)
{
  double r = 0.0;
  for (int i = 0; i < x->size(); ++i) r += (*x)[i] * y[i];
  return r;
}

/*
Purpose:
    Computes the optimal portfolio for a given risk

Input:
    n: Number of assets
    mu: An n dimmensional vector of expected returns
    GT: A matrix with n columns so (GT')*GT  = covariance matrix
    x0: Initial holdings
    w: Initial cash holding
    gamma: Maximum risk (=std. dev) accepted

Output:
    Optimal expected return and the optimal portfolio
*/
double BasicMarkowitz
( int                                n,
  std::shared_ptr<ndarray<double, 1>> mu,
  std::shared_ptr<ndarray<double, 2>> GT,
  std::shared_ptr<ndarray<double, 1>> x0,
  double                             w,
  double                             gamma)
{
  Model::t M = new Model("Basic Markowitz"); auto _M = finally([&]() { M->dispose(); });
  // Redirect log output from the solver to stdout for debugging.
  // M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; } );

  // Defines the variables (holdings). Shortselling is not allowed.
  Variable::t x = M->variable("x", n, Domain::greaterThan(0.0));

  //  Maximize expected return
  M->objective("obj", ObjectiveSense::Maximize, Expr::dot(mu, x));

  // The amount invested  must be identical to intial wealth
  M->constraint("budget", Expr::sum(x), Domain::equalsTo(w + sum(x0)));

  // Imposes a bound on the risk
  M->constraint("risk", Expr::vstack(gamma, Expr::mul(GT, x)), Domain::inQCone());

  // Solves the model.
  M->solve();

  return dot(mu, x->level());
}


/*
  The example reads in data and solves the portfolio models.
 */
int main(int argc, char ** argv)
{

  int        n      = 3;
  double     w      = 1.0;
  auto       mu     = new_array_ptr<double, 1>( {0.1073, 0.0737, 0.0627} );
  auto       x0     = new_array_ptr<double, 1>({0.0, 0.0, 0.0});
  auto       gammas = new_array_ptr<double, 1>({0.035, 0.040, 0.050, 0.060, 0.070, 0.080, 0.090});
  auto       GT     = new_array_ptr<double, 2>({
    { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
    { 0.0              , 0.102863378954911  , -0.00222873156550421},
    { 0.0              , 0.0                ,  0.0338148677744977 }
  });

  std::cout << std::endl << std::endl
            << "================================" << std::endl
            << "Markowitz portfolio optimization" << std::endl
            << "================================" << std::endl;

  std::cout << std::endl
            << "-----------------------------------------------------------------------------------" << std::endl
            << "Basic Markowitz portfolio optimization" << std::endl
            << "-----------------------------------------------------------------------------------" << std::endl;

  std::cout << std::setprecision(4)
            << std::setiosflags(std::ios::scientific);

  for (auto gamma : *gammas)
    std::cout << "Expected return: " << BasicMarkowitz( n, mu, GT, x0, w, gamma) << " St deviation: " << gamma << std::endl;

  return 0;
}

portfolio_2_frontier.cc

Listing 16.25 portfolio_2_frontier.cc Click here to download.
/*
  File : portfolio_2_frontier.cc

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model.
                 Determines points on the efficient frontier.
*/

#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include "monty.h"
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

static double sum(std::shared_ptr<ndarray<double, 1>> x)
{
  double r = 0.0;
  for (auto v : *x) r += v;
  return r;
}

static double dot(std::shared_ptr<ndarray<double, 1>> x,
                  std::shared_ptr<ndarray<double, 1>> y)
{
  double r = 0.0;
  for (int i = 0; i < x->size(); ++i) r += (*x)[i] * (*y)[i];
  return r;
}

static double dot(std::shared_ptr<ndarray<double, 1>> x,
                  std::vector<double> & y)
{
  double r = 0.0;
  for (int i = 0; i < x->size(); ++i) r += (*x)[i] * y[i];
  return r;
}

/*
  Purpose:
      Computes several portfolios on the optimal portfolios by

          for alpha in alphas:
              maximize   expected return - alpha * variance
              subject to the constraints

  Input:
      n: Number of assets
      mu: An n dimmensional vector of expected returns
      GT: A matrix with n columns so (GT')*GT  = covariance matrix
      x0: Initial holdings
      w: Initial cash holding
      alphas: List of the alphas

  Output:
      The efficient frontier as list of tuples (alpha, expected return, variance)
 */
void EfficientFrontier
( int n,
  std::shared_ptr<ndarray<double, 1>> mu,
  std::shared_ptr<ndarray<double, 2>> GT,
  std::shared_ptr<ndarray<double, 1>> x0,
  double w,
  std::vector<double> & alphas,
  std::vector<double> & frontier_mux,
  std::vector<double> & frontier_s)
{

  Model::t M = new Model("Efficient frontier");  auto M_ = finally([&]() { M->dispose(); });

  // Defines the variables (holdings). Shortselling is not allowed.
  Variable::t x = M->variable("x", n, Domain::greaterThan(0.0)); // Portfolio variables
  Variable::t s = M->variable("s", 1, Domain::unbounded());      // Variance variable

  M->constraint("budget", Expr::sum(x), Domain::equalsTo(w + sum(x0)));

  // Computes the risk
  M->constraint("variance", Expr::vstack(s, 0.5, Expr::mul(GT, x)), Domain::inRotatedQCone());

  //  Define objective as a weighted combination of return and variance
  Parameter::t alpha = M->parameter();
  M->objective("obj", ObjectiveSense::Maximize, Expr::sub(Expr::dot(mu, x), Expr::mul(alpha, s)));

  // Solve the same problem for many values of parameter alpha
  for (double a : alphas) {
      alpha->setValue(a);
      M->solve();

      frontier_mux.push_back(dot(mu, x->level()));
      frontier_s.push_back((*s->level())[0]);
  }
}

/*
  The example reads in data and solves the portfolio models.
 */
int main(int argc, char ** argv)
{

  int        n      = 3;
  double     w      = 1.0;
  auto       mu     = new_array_ptr<double, 1>( {0.1073, 0.0737, 0.0627} );
  auto       x0     = new_array_ptr<double, 1>({0.0, 0.0, 0.0});
  auto       GT     = new_array_ptr<double, 2>({
    { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
    { 0.0              , 0.102863378954911  , -0.00222873156550421},
    { 0.0              , 0.0                ,  0.0338148677744977 }
  });

  std::cout << std::endl << std::endl
            << "================================" << std::endl
            << "Markowitz portfolio optimization" << std::endl
            << "================================" << std::endl;

  std::cout << std::setprecision(4)
            << std::setiosflags(std::ios::scientific);

  // Some predefined alphas are chosen
  std::vector<double> alphas{ 0.0, 0.01, 0.1, 0.25, 0.30, 0.35, 0.4, 0.45, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 10.0 };
  std::vector<double> frontier_mux;
  std::vector<double> frontier_s;

  EfficientFrontier(n, mu, GT, x0, w, alphas, frontier_mux, frontier_s);
  std::cout << "\n-----------------------------------------------------------------------------------" << std::endl
            << "Efficient frontier" << std::endl
            << "-----------------------------------------------------------------------------------" << std::endl
            << std::endl;
  std::cout << std::setw(12) << "alpha" << std::setw(12) << "return" << std::setw(12) << "variance" << std::endl;
  for (int i = 0; i < frontier_mux.size(); ++i)
    std::cout << std::setw(12) << alphas[i] << std::setw(12) << frontier_mux[i] << std::setw(12) << frontier_s[i] << std::endl;

  return 0;
}

portfolio_3_impact.cc

Listing 16.26 portfolio_3_impact.cc Click here to download.
/*
  File : portfolio_3_impact.cc

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model
                 with transaction costs of order x^(3/2).
*/

#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include "monty.h"
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

static double sum(std::shared_ptr<ndarray<double, 1>> x)
{
  double r = 0.0;
  for (auto v : *x) r += v;
  return r;
}

static double dot(std::shared_ptr<ndarray<double, 1>> x,
                  std::shared_ptr<ndarray<double, 1>> y)
{
  double r = 0.0;
  for (int i = 0; i < x->size(); ++i) r += (*x)[i] * (*y)[i];
  return r;
}

static double dot(std::shared_ptr<ndarray<double, 1>> x,
                  std::vector<double> & y)
{
  double r = 0.0;
  for (int i = 0; i < x->size(); ++i) r += (*x)[i] * y[i];
  return r;
}


/*
    Description:
        Extends the basic Markowitz model with a market cost term.

    Input:
        n: Number of assets
        mu: An n dimmensional vector of expected returns
        GT: A matrix with n columns so (GT')*GT  = covariance matrix'
        x0: Initial holdings
        w: Initial cash holding
        gamma: Maximum risk (=std. dev) accepted
        m: It is assumed that  market impact cost for the j'th asset is
           m_j|x_j-x0_j|^3/2

    Output:
       Optimal expected return and the optimal portfolio

*/
void MarkowitzWithMarketImpact
( int n,
  std::shared_ptr<ndarray<double, 1>> mu,
  std::shared_ptr<ndarray<double, 2>> GT,
  std::shared_ptr<ndarray<double, 1>> x0,
  double      w,
  double      gamma,
  std::shared_ptr<ndarray<double, 1>> m,
  std::vector<double> & xsol,
  std::vector<double> & tsol)
{
  Model::t M = new Model("Markowitz portfolio with market impact");  auto M_ = finally([&]() { M->dispose(); });

  // Defines the variables. No shortselling is allowed.
  Variable::t x = M->variable("x", n, Domain::greaterThan(0.0));

  // Variables computing the impact cost
  Variable::t t = M->variable("t", n, Domain::unbounded());

  //  Maximize expected return
  M->objective("obj", ObjectiveSense::Maximize, Expr::dot(mu, x));

  // Invested amount + slippage cost = initial wealth
  M->constraint("budget", Expr::add(Expr::sum(x), Expr::dot(m, t)), Domain::equalsTo(w + sum(x0)));

  // Imposes a bound on the risk
  M->constraint("risk", Expr::vstack( gamma, Expr::mul(GT, x)),
                Domain::inQCone());

  // t >= |x-x0|^1.5, using a power cone
  M->constraint("tz", Expr::hstack(t, Expr::constTerm(n, 1.0), Expr::sub(x, x0)), Domain::inPPowerCone(2.0/3.0));

  M->solve();

  xsol.resize(n);
  tsol.resize(n);
  auto xlvl = x->level();
  auto tlvl = t->level();

  std::copy(xlvl->flat_begin(), xlvl->flat_end(), xsol.begin());
  std::copy(tlvl->flat_begin(), tlvl->flat_end(), tsol.begin());
}


/*
  The example reads in data and solves the portfolio models.
 */
int main(int argc, char ** argv)
{

  int        n      = 3;
  double     w      = 1.0;
  auto       mu     = new_array_ptr<double, 1>( {0.1073, 0.0737, 0.0627} );
  auto       x0     = new_array_ptr<double, 1>({0.0, 0.0, 0.0});
  auto       gammas = new_array_ptr<double, 1>({0.035, 0.040, 0.050, 0.060, 0.070, 0.080, 0.090});
  auto       GT     = new_array_ptr<double, 2>({
    { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
    { 0.0              , 0.102863378954911  , -0.00222873156550421},
    { 0.0              , 0.0                ,  0.0338148677744977 }
  });

  std::cout << std::endl << std::endl
            << "================================" << std::endl
            << "Markowitz portfolio optimization" << std::endl
            << "================================" << std::endl;

  std::cout << std::setprecision(4)
            << std::setiosflags(std::ios::scientific);

  // Somewhat arbirtrary choice of m
  std::shared_ptr<ndarray<double, 1>> m(new ndarray<double, 1>(shape_t<1>(n), 1.0e-2));
  double gamma = (*gammas)[0];
  std::vector<double> x;
  std::vector<double> t;

  MarkowitzWithMarketImpact(n, mu, GT, x0, w, gamma, m, x, t);

  std::cout << std::resetiosflags(std::ios::left);
  std::cout << std::endl
            << "-----------------------------------------------------------------------------------" << std::endl
            << "Markowitz portfolio optimization with market impact cost" << std::endl
            << "-----------------------------------------------------------------------------------" << std::endl
            << std::endl
            << "Expected return: " << dot(mu, x) << " St deviation: " << gamma << " Market impact cost: " << dot(m, t) << std::endl;

  return 0;
}

portfolio_4_transcost.cc

Listing 16.27 portfolio_4_transcost.cc Click here to download.
/*
  File : portfolio_4_transcost.cc

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model
                 with fixed setup costs and transaction costs
                 as a mixed-integer problem.
*/

#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include "monty.h"
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

static double sum(std::shared_ptr<ndarray<double, 1>> x)
{
  double r = 0.0;
  for (auto v : *x) r += v;
  return r;
}

static double dot(std::shared_ptr<ndarray<double, 1>> x,
                  std::shared_ptr<ndarray<double, 1>> y)
{
  double r = 0.0;
  for (int i = 0; i < x->size(); ++i) r += (*x)[i] * (*y)[i];
  return r;
}

static double dot(std::shared_ptr<ndarray<double, 1>> x,
                  std::vector<double> & y)
{
  double r = 0.0;
  for (int i = 0; i < x->size(); ++i) r += (*x)[i] * y[i];
  return r;
}


/*
    Description:
        Extends the basic Markowitz model with a market cost term.

    Input:
        n: Number of assets
        mu: An n dimmensional vector of expected returns
        GT: A matrix with n columns so (GT')*GT  = covariance matrix
        x0: Initial holdings
        w: Initial cash holding
        gamma: Maximum risk (=std. dev) accepted
        f: If asset j is traded then a fixed cost f_j must be paid
        g: If asset j is traded then a cost g_j must be paid for each unit traded

    Output:
       Optimal expected return and the optimal portfolio

*/
std::shared_ptr<ndarray<double, 1>> MarkowitzWithTransactionsCost
                                 ( int n,
                                   std::shared_ptr<ndarray<double, 1>> mu,
                                   std::shared_ptr<ndarray<double, 2>> GT,
                                   std::shared_ptr<ndarray<double, 1>> x0,
                                   double                             w,
                                   double                             gamma,
                                   std::shared_ptr<ndarray<double, 1>> f,
                                   std::shared_ptr<ndarray<double, 1>> g)
{
  // Upper bound on the traded amount
  std::shared_ptr<ndarray<double, 1>> u(new ndarray<double, 1>(shape_t<1>(n), w + sum(x0)));

  Model::t M = new Model("Markowitz portfolio with transaction costs");  auto M_ = finally([&]() { M->dispose(); });

  // Defines the variables. No shortselling is allowed.
  Variable::t x = M->variable("x", n, Domain::greaterThan(0.0));

  // Addtional "helper" variables
  Variable::t z = M->variable("z", n, Domain::unbounded());
  // Binary varables
  Variable::t y = M->variable("y", n, Domain::binary());

  //  Maximize expected return
  M->objective("obj", ObjectiveSense::Maximize, Expr::dot(mu, x));

  // Invest amount + transactions costs = initial wealth
  M->constraint("budget", Expr::add(Expr::add(Expr::sum(x), Expr::dot(f, y)), Expr::dot(g, z)),
                Domain::equalsTo(w + sum(x0)));

  // Imposes a bound on the risk
  M->constraint("risk", Expr::vstack( gamma, Expr::mul(GT, x)),
                Domain::inQCone());

  // z >= |x-x0|
  M->constraint("buy", Expr::sub(z, Expr::sub(x, x0)), Domain::greaterThan(0.0));
  M->constraint("sell", Expr::sub(z, Expr::sub(x0, x)), Domain::greaterThan(0.0));

  // Consraints for turning y off and on. z-diag(u)*y<=0 i.e. z_j <= u_j*y_j
  M->constraint("y_on_off", Expr::sub(z, Expr::mul(Matrix::diag(u), y)), Domain::lessThan(0.0));

  // Integer optimization problems can be very hard to solve so limiting the
  // maximum amount of time is a valuable safe guard
  M->setSolverParam("mioMaxTime", 180.0);
  M->solve();

  return x->level();
}


/*
  The example reads in data and solves the portfolio models.
 */
int main(int argc, char ** argv)
{

  int        n      = 3;
  double     w      = 1.0;
  auto       mu     = new_array_ptr<double, 1>( {0.1073, 0.0737, 0.0627} );
  auto       x0     = new_array_ptr<double, 1>({0.0, 0.0, 0.0});
  auto       gammas = new_array_ptr<double, 1>({0.035, 0.040, 0.050, 0.060, 0.070, 0.080, 0.090});
  auto       GT     = new_array_ptr<double, 2>({
    { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
    { 0.0              , 0.102863378954911  , -0.00222873156550421},
    { 0.0              , 0.0                ,  0.0338148677744977 }
  });

  std::cout << std::endl << std::endl
            << "================================" << std::endl
            << "Markowitz portfolio optimization" << std::endl
            << "================================" << std::endl;

  std::cout << std::setprecision(4)
            << std::setiosflags(std::ios::scientific);

  std::shared_ptr<ndarray<double, 1>> f(new ndarray<double, 1>(shape_t<1>(n), 0.01));
  std::shared_ptr<ndarray<double, 1>> g(new ndarray<double, 1>(shape_t<1>(n), 0.001));
  std::cout << std::endl
            << "-----------------------------------------------------------------------------------" << std::endl
            << "Markowitz portfolio optimization with transaction cost" << std::endl
            << "-----------------------------------------------------------------------------------" << std::endl
            << std::endl;

  double gamma = (*gammas)[0];
  auto x = MarkowitzWithTransactionsCost(n, mu, GT, x0, w, gamma, f, g);
  std::cout << "Optimal portfolio:" << std::endl;
  for ( int i = 0; i < x->size(); ++i)
    std::cout << "\tx[" << std::setw(2) << i << "]  " << std::setw(12) << (*x)[i] << std::endl;

  return 0;
}

portfolio_5_card.cc

Listing 16.28 portfolio_5_card.cc Click here to download.
/*
  File : portfolio_5_card.cc

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model
                 with cardinality constraints on number of assets traded.
*/

#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include "monty.h"
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

static double sum(std::shared_ptr<ndarray<double, 1>> x)
{
  double r = 0.0;
  for (auto v : *x) r += v;
  return r;
}

static double dot(std::shared_ptr<ndarray<double, 1>> x,
                  std::shared_ptr<ndarray<double, 1>> y)
{
  double r = 0.0;
  for (int i = 0; i < x->size(); ++i) r += (*x)[i] * (*y)[i];
  return r;
}

static double dot(std::shared_ptr<ndarray<double, 1>> x,
                  std::vector<double> & y)
{
  double r = 0.0;
  for (int i = 0; i < x->size(); ++i) r += (*x)[i] * y[i];
  return r;
}


/*
    Description:
        Extends the basic Markowitz model with cardinality constraints.

    Input:
        n: Number of assets
        mu: An n dimmensional vector of expected returns
        GT: A matrix with n columns so (GT')*GT  = covariance matrix
        x0: Initial holdings
        w: Initial cash holding
        gamma: Maximum risk (=std. dev) accepted
        k: Maximal number of assets in which we allow to change position.

    Output:
       Optimal expected return and the optimal portfolio

*/
std::vector<std::vector<double>> MarkowitzWithCardinality
                                 ( int n,
                                   std::shared_ptr<ndarray<double, 1>> mu,
                                   std::shared_ptr<ndarray<double, 2>> GT,
                                   std::shared_ptr<ndarray<double, 1>> x0,
                                   double                              w,
                                   double                              gamma,
                                   std::vector<int>                    kValues)
{
  // Upper bound on the traded amount
  std::shared_ptr<ndarray<double, 1>> u(new ndarray<double, 1>(shape_t<1>(n), w + sum(x0)));

  Model::t M = new Model("Markowitz portfolio with cardinality constraints");  auto M_ = finally([&]() { M->dispose(); });

  // Defines the variables. No shortselling is allowed.
  Variable::t x = M->variable("x", n, Domain::greaterThan(0.0));

  // Addtional "helper" variables
  Variable::t z = M->variable("z", n, Domain::unbounded());
  // Binary varables
  Variable::t y = M->variable("y", n, Domain::binary());

  //  Maximize expected return
  M->objective("obj", ObjectiveSense::Maximize, Expr::dot(mu, x));

  // The amount invested  must be identical to initial wealth
  M->constraint("budget", Expr::sum(x), Domain::equalsTo(w + sum(x0)));

  // Imposes a bound on the risk
  M->constraint("risk", Expr::vstack( gamma, Expr::mul(GT, x)),
                Domain::inQCone());

  // z >= |x-x0|
  M->constraint("buy", Expr::sub(z, Expr::sub(x, x0)), Domain::greaterThan(0.0));
  M->constraint("sell", Expr::sub(z, Expr::sub(x0, x)), Domain::greaterThan(0.0));

  // Consraints for turning y off and on. z-diag(u)*y<=0 i.e. z_j <= u_j*y_j
  M->constraint("y_on_off", Expr::sub(z, Expr::mul(Matrix::diag(u), y)), Domain::lessThan(0.0));

  // At most k assets change position
  auto cardMax = M->parameter();
  M->constraint("cardinality", Expr::sub(Expr::sum(y), cardMax), Domain::lessThan(0));

  // Integer optimization problems can be very hard to solve so limiting the
  // maximum amount of time is a valuable safe guard
  M->setSolverParam("mioMaxTime", 180.0);

  // Solve multiple instances by varying the cardinality bound
  std::vector<std::vector<double>> results;

  for(auto k : kValues) {
    cardMax->setValue(k);
    M->solve();
    auto sol = x->level();
    results.push_back(new_vector_from_array_ptr(sol));
  }

  return results;
}


/*
  The example reads in data and solves the portfolio models.
 */
int main(int argc, char ** argv)
{

  int        n      = 3;
  double     w      = 1.0;
  auto       mu     = new_array_ptr<double, 1>( {0.1073, 0.0737, 0.0627} );
  auto       x0     = new_array_ptr<double, 1>({0.0, 0.0, 0.0});
  auto       gamma  = 0.035;
  auto       GT     = new_array_ptr<double, 2>({
    { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
    { 0.0              , 0.102863378954911  , -0.00222873156550421},
    { 0.0              , 0.0                ,  0.0338148677744977 }
  });
  std::vector<int> kValues = { 1, 2, 3 };

  std::cout << std::endl << std::endl
            << "================================" << std::endl
            << "Markowitz portfolio optimization" << std::endl
            << "================================" << std::endl;

  std::cout << std::setprecision(4)
            << std::setiosflags(std::ios::scientific);

  std::cout << std::endl
            << "-----------------------------------------------------------------------------------" << std::endl
            << "Markowitz portfolio optimization with cardinality bounds" << std::endl
            << "-----------------------------------------------------------------------------------" << std::endl
            << std::endl;

  auto results = MarkowitzWithCardinality(n, mu, GT, x0, w, gamma, kValues);

  for(int k=1; k<=n; k++)
  {
    std::cout << "Bound " << k << "  Portfolio: ";
    for (int i = 0; i < n; ++i)
      std::cout << results[k-1][i] << " ";
    std::cout << std::endl;
  }

  return 0;
}

pow1.cc

Listing 16.29 pow1.cc Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      pow1.cc
//
//   Purpose: Demonstrates how to solve the problem
//
//     maximize x^0.2*y^0.8 + z^0.4 - x
//           st x + y + 0.5z = 2
//              x,y,z >= 0
//
#include <iostream>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  Model::t M = new Model("pow1"); auto _M = finally([&]() { M->dispose(); });

  Variable::t x  = M->variable("x", 3, Domain::unbounded());
  Variable::t x3 = M->variable();
  Variable::t x4 = M->variable();

  // Create the linear constraint
  auto aval = new_array_ptr<double, 1>({1.0, 1.0, 0.5});
  M->constraint(Expr::dot(x, aval), Domain::equalsTo(2.0));

  // Create the conic constraints
  M->constraint(Var::vstack(x->slice(0,2), x3), Domain::inPPowerCone(0.2));
  M->constraint(Expr::vstack(x->index(2), 1.0, x4), Domain::inPPowerCone(0.4));     

  auto cval = new_array_ptr<double, 1>({1.0, 1.0, -1.0});
  M->objective(ObjectiveSense::Maximize, Expr::dot(cval, Var::vstack(x3, x4, x->index(0))));

  // Solve the problem
  M->solve();

  // Get the linear solution values
  ndarray<double, 1> xlvl   = *(x->level());
  std::cout << "x,y,z = " << xlvl << std::endl;
}

primal_svm.cc

Listing 16.30 primal_svm.cc Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      primal_svm.cc
//
// Purpose: Implements a simple soft-margin SVM
//          using the Fusion API.

#include <random>
#include <memory>
#include <iostream>
#include <iterator>
#include <iomanip>

#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  int m = 50 ;
  int n = 3;
  int nc = 10;

  int nump = m / 2;
  int numm = m - nump;

  auto y = new_array_ptr<double, 1> (m);
  std::fill( y->begin(), y->begin() + nump, 1.);
  std::fill( y->begin() + nump, y->end(), -1.);

  double mean = 1.;
  double var = 1.;

  auto X = std::shared_ptr< ndarray<double, 2> > (  new ndarray<double, 2> ( shape_t<2>(m, n) ) );

  std::mt19937 e2(0);

  for (int i = 0; i < nump; i++)
  {
    auto ram = std::bind(std::normal_distribution<>(mean, var), e2);
    for ( int j = 0; j < n; j++)
      (*X)(i, j) = ram();
  }

  std::cout << "Number of data    : " << m << std::endl;
  std::cout << "Number of features: " << n << std::endl;
  Model::t M = new Model("primal SVM"); auto _M = finally([&]() { M->dispose(); });

  Variable::t w  =  M->variable( n, Domain::unbounded());
  Variable::t t  =  M->variable( 1, Domain::unbounded());
  Variable::t b  =  M->variable( 1, Domain::unbounded());
  Variable::t xi =  M->variable( m, Domain::greaterThan(0.));

  auto ex = Expr::sub( Expr::mul(X, w), Var::repeat(b, m) );
  M->constraint( Expr::add(Expr::mulElm( y, ex ), xi ) , Domain::greaterThan( 1. ) );

  M->constraint( Expr::vstack(1., t, w) , Domain::inRotatedQCone() );

  std::cout << "   c   | b      | w\n";
  for (int i = 0; i < nc; i++)
  {
    double c = 500.0 * i;
    M->objective( ObjectiveSense::Minimize, Expr::add( t, Expr::mul(c, Expr::sum(xi) ) ) );
    M->solve();

    try
    {
      std::cout << std::setw(6) << c << " | " << std::setw(8) << (*(b->level())) [0] << " | ";
      std::cout.width(8);
      auto wlev = w->level();
      std::copy( wlev->begin(), wlev->end() , std::ostream_iterator<double>(std::cout, " ") );
      std::cout << "\n";
    }
    catch (...) {}

  }
  return 0;
}

qcqp_sdo_relaxation.cc

Listing 16.31 qcqp_sdo_relaxation.cc Click here to download.
/*
* Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
*
* File:      qcqp_sdo_relaxation.cc
*
* Purpose:   Demonstrate how to use SDP to solve
*            convex relaxation of a mixed-integer QCQP
**/
#include <random>
#include <iostream>
#include <ctime>
#include "fusion.h"

using namespace mosek::fusion;
using namespace mosek::LinAlg;
using namespace monty;


Model::t miqcqp_sdo_relaxation(int n, Matrix::t P, const std::shared_ptr<ndarray<double, 1>> & q) {
  Model::t M = new Model();

  Variable::t Z = M->variable("Z", Domain::inPSDCone(n + 1));

  Variable::t X = Z->slice(new_array_ptr<int, 1>({0, 0}), new_array_ptr<int, 1>({n, n}));
  Variable::t x = Z->slice(new_array_ptr<int, 1>({0, n}), new_array_ptr<int, 1>({n, n + 1}));

  M->constraint( Expr::sub(X->diag(), x), Domain::greaterThan(0.) );
  M->constraint( Z->index(n, n), Domain::equalsTo(1.) );

  M->objective( ObjectiveSense::Minimize, Expr::add(
                  Expr::sum( Expr::mulElm( P, X ) ),
                  Expr::mul( 2.0, Expr::dot(x, q) )
                ) );

  return M;
}

// A direct integer model for minimizing |Ax-b|
Model::t int_least_squares(int n, Matrix::t A, const std::shared_ptr<ndarray<double, 1>> & b) {
  Model::t M = new Model();

  Variable::t x = M->variable("x", n, Domain::integral(Domain::unbounded()));
  Variable::t t = M->variable("t", 1, Domain::unbounded());

  M->constraint( Expr::vstack(t, Expr::sub(Expr::mul(A, x), b)), Domain::inQCone() );
  M->objective( ObjectiveSense::Minimize, t );

  return M;
}


int main(int argc, char ** argv)
{
  std::default_random_engine generator;
  generator.seed(time(0));
  std::uniform_real_distribution<double> unif_distr(0., 1.);
  std::normal_distribution<double> normal_distr(0., 1.);

  // problem dimensions
  int n = 20;
  int m = 2 * n;

  auto c = new_array_ptr<double, 1>(n);
  auto A = new_array_ptr<double, 1>(n * m);
  auto P = new_array_ptr<double, 1>(n * n);
  auto b = new_array_ptr<double, 1>(m);
  auto q = new_array_ptr<double, 1>(n);


  std::generate(A->begin(), A->end(), std::bind(normal_distr, generator));
  std::generate(c->begin(), c->end(), std::bind(unif_distr, generator));
  std::fill(b->begin(), b->end(), 0.0);
  std::fill(q->begin(), q->end(), 0.0);

  // P = A^T A
  syrk(MSK_UPLO_LO, MSK_TRANSPOSE_YES,
       n, m, 1.0, A, 0., P);
  for (int j = 0; j < n; j++) for (int i = j + 1; i < n; i++) (*P)[i * n + j] = (*P)[j * n + i];

  // q = -P c, b = A c
  gemv(MSK_TRANSPOSE_NO, n, n, -1.0, P, c, 0., q);
  gemv(MSK_TRANSPOSE_NO, m, n, 1.0, A, c, 0., b);

  // Solve the problems
  {
    Model::t M = miqcqp_sdo_relaxation(n, Matrix::dense(n, n, P), q);
    Model::t Mint = int_least_squares(n, Matrix::dense(n, m, A)->transpose(), b);
    M->solve();
    Mint->solve();

    auto xRound = M->getVariable("Z")->
                  slice(new_array_ptr<int, 1>({0, n}), new_array_ptr<int, 1>({n, n + 1}))->level();
    for (int i = 0; i < n; i++) (*xRound)[i] = round((*xRound)[i]);
    auto yRound = new_array_ptr<double, 1>(m);
    auto xOpt   = Mint->getVariable("x")->level();
    auto yOpt   = new_array_ptr<double, 1>(m);
    std::copy(b->begin(), b->end(), yRound->begin());
    std::copy(b->begin(), b->end(), yOpt->begin());
    gemv(MSK_TRANSPOSE_NO, m, n, 1.0, A, xRound, -1.0, yRound);        // Ax_round-b
    gemv(MSK_TRANSPOSE_NO, m, n, 1.0, A, xOpt, -1.0, yOpt);            // Ax_opt-b

    std::cout << M->getSolverDoubleInfo("optimizerTime") << " " << Mint->getSolverDoubleInfo("optimizerTime") << "\n";
    double valRound, valOpt;
    dot(m, yRound, yRound, valRound); dot(m, yOpt, yOpt, valOpt);
    std::cout << sqrt(valRound) << " " << sqrt(valOpt) << "\n";
  }
  return 0;
};

reoptimization.cc

Listing 16.32 reoptimization.cc Click here to download.
/**
  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File :      reoptimization.cc

  Purpose:   Demonstrates how to solve a  linear
             optimization problem using the MOSEK API
             and modify and re-optimize the problem.
*/
#include <iostream>
#include <cmath>

#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

void printsol(const std::shared_ptr<ndarray<double, 1>> & a) {
  std::cout << "x = ";
  for(auto val : *a) std::cout << val << " ";
  std::cout << "\n";
}

int main() {
  auto c = new_array_ptr<double, 1>({ 1.5, 2.5, 3.0 });
  auto A = new_array_ptr<double, 2>({ {2, 4, 3},
                                      {3, 2, 3},
                                      {2, 3, 2} });
  auto b = new_array_ptr<double, 1>({ 100000.0, 50000.0, 60000.0 });
  int numvar = 3;
  int numcon = 3;

  // Create a model and input data
  Model::t M = new Model(); auto M_ = monty::finally([&]() { M->dispose(); });

  auto x = M->variable(numvar, Domain::greaterThan(0.0));
  auto con = M->constraint(Expr::mul(A, x), Domain::lessThan(b));
  M->objective(ObjectiveSense::Maximize, Expr::dot(c, x));
  // Solve the problem
  M->solve();
  printsol(x->level());

  /***************** Change an entry in the A matrix ********************/
  con->index(0)->update(Expr::mul(3.0, x->index(0)), x->index(0));
  M->solve();
  printsol(x->level());

  /*************** Add a new variable ******************************/
  // Create a variable and a compound view of all variables
  auto x3 = M->variable(Domain::greaterThan(0.0));
  auto xNew = Var::vstack(x, x3);
  // Add to the exising constraint
  con->update(Expr::mul(x3, new_array_ptr<double, 1>({ 4, 0, 1 })),x3);
  // Change the objective to include x3
  M->objective(ObjectiveSense::Maximize, Expr::dot(new_array_ptr<double, 1>({1.5,2.5,3.0,1.0}), xNew));
  M->solve();
  printsol(xNew->level());

  /**************** Add a new constraint *****************************/
  auto con2 = M->constraint(Expr::dot(xNew, new_array_ptr<double, 1>({1, 2, 1, 1})), Domain::lessThan(30000.0));
  M->solve();
  printsol(xNew->level());

  /**************** Change constraint bounds *****************************/
  // Assemble all constraints we previously defined into one
  auto cAll = Constraint::vstack(con, con2);
  // Change bounds by effectively updating fixed terms with the difference
  cAll->update(new_array_ptr<double, 1>({20000, 10000, 10000, 8000}));
  M->solve();
  printsol(xNew->level());
}

response.cc

Listing 16.33 response.cc Click here to download.
/*
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:      response.cc

  Purpose:   This example demonstrates proper response handling
             for problems solved with the interior-point optimizers.
*/

#include <iostream>
#include <random>
#include <sstream>
#include <vector>
#include <utility>
#include "fusion.h"
using namespace mosek::fusion;
using namespace monty;

// Set up a small artificial conic problem to test with
void setupExample(Model::t M)
{
  auto x = M->variable("x", 3, Domain::greaterThan(0.0));
  auto y = M->variable("y", 3, Domain::unbounded());
  auto z1 = Var::vstack(y->index(0), x->slice(0, 2));
  auto z2 = Var::vstack(y->slice(1, 3), x->index(2));
  M->constraint("lc", Expr::dot(new_array_ptr<double, 1>({1.0, 1.0, 2.0}), x), Domain::equalsTo(1.0));
  M->constraint("qc1", z1, Domain::inQCone());
  M->constraint("qc2", z2, Domain::inRotatedQCone());
  M->objective("obj", ObjectiveSense::Minimize, Expr::sum(y));
}

int main(int arc, char** argv)
{
  Model::t M = new Model(); auto _M = finally([&]() { M->dispose(); });

  // (Optional) set a log stream
  // M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; });

  // (Optional) uncomment to see what happens when solution status is unknown
  // M->setSolverParam("intpntMaxIterations", 1);

  // In this example we set up a small conic problem
  setupExample(M);

  // Optimize
  try 
  {
    M->solve();

    // We expect solution status OPTIMAL (this is also default)
    M->acceptedSolutionStatus(AccSolutionStatus::Optimal);

    auto x = M->getVariable("x");
    auto xsize = x->getSize();
    auto xVal = x->level();
    std::cout << "Optimal value of x = ";
    for(int i = 0; i < xsize; ++i)
      std::cout << (*xVal)[i] << " ";
    std::cout << "\nOptimal primal objective: " << M->primalObjValue() <<"\n";
    // .. continue analyzing the solution

  }
  catch (const OptimizeError& e)
  {
    std::cout << "Optimization failed. Error: " << e.what() << "\n";
  }
  catch (const SolutionError& e)
  {
    // The solution with at least the expected status was not available.
    // We try to diagnoze why.
    std::cout << "Requested solution was not available.\n";
    auto prosta = M->getProblemStatus();
    switch(prosta)
    {
      case ProblemStatus::DualInfeasible:
        std::cout << "Dual infeasibility certificate found.\n";
        break;

      case ProblemStatus::PrimalInfeasible:
        std::cout << "Primal infeasibility certificate found.\n";
        break;

      case ProblemStatus::Unknown:
        // The solutions status is unknown. The termination code
        // indicates why the optimizer terminated prematurely.
        std::cout << "The solution status is unknown.\n";
        char symname[MSK_MAX_STR_LEN];
        char desc[MSK_MAX_STR_LEN];
        MSK_getcodedesc((MSKrescodee)(M->getSolverIntInfo("optimizeResponse")), symname, desc);
        std::cout << "  Termination code: " << symname << " " << desc << "\n";
        break;

      default:
        std::cout << "Another unexpected problem status: " << prosta << "\n";
    }
  }
  catch (const std::exception& e)
  {
    std::cerr << "Unexpected error: " << e.what() << "\n";
  }

  M->dispose();
  return 0;
 }

sdo1.cc

Listing 16.34 sdo1.cc Click here to download.
//
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      sdo1.cc
//
//  Purpose:
//  Solves the mixed semidefinite and conic quadratic optimization problem
//
//                   [2, 1, 0]
//    minimize    Tr [1, 2, 1] * X + x0
//                   [0, 1, 2]
//
//                   [1, 0, 0]
//    subject to  Tr [0, 1, 0] * X + x0           = 1.0
//                   [0, 0, 1]
//
//                   [1, 1, 1]
//                Tr [1, 1, 1] * X      + x1 + x2 = 0.5
//                   [1, 1, 1]
//
//                   X >> 0,  x0 >= (x1^2 + x2^2) ^ (1/2)
//
#include <iostream>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

int main(int argc, char ** argv)
{
  Model::t M = new Model("sdo1"); auto _M = finally([&]() { M->dispose(); });

  // Setting up the variables
  Variable::t X  = M->variable("X", Domain::inPSDCone(3));
  Variable::t x  = M->variable("x", Domain::inQCone(3));

  // Setting up the constant coefficient matrices
  Matrix::t C  = Matrix::dense ( new_array_ptr<double, 2>({{2., 1., 0.}, {1., 2., 1.}, {0., 1., 2.}}));
  Matrix::t A1 = Matrix::eye(3);
  Matrix::t A2 = Matrix::ones(3, 3);

  // Objective
  M->objective(ObjectiveSense::Minimize, Expr::add(Expr::dot(C, X), x->index(0)));

  // Constraints
  M->constraint("c1", Expr::add(Expr::dot(A1, X), x->index(0)), Domain::equalsTo(1.0));
  M->constraint("c2", Expr::add(Expr::dot(A2, X), Expr::sum(x->slice(1, 3))), Domain::equalsTo(0.5));

  M->solve();

  std::cout << "Solution : " << std::endl;
  std::cout << "  X = " << *(X->level()) << std::endl;
  std::cout << "  x = " << *(x->level()) << std::endl;

  return 0;
}

sdo2.cc

Listing 16.35 sdo2.cc Click here to download.
/*
   Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
 
   File :      sdo2.cc
 
   Purpose :   Solves the semidefinite problem with two symmetric variables:
 
                  min   <C1,X1> + <C2,X2>
                  st.   <A1,X1> + <A2,X2> = b
                              (X2)_{1,2} <= k
                 
                  where X1, X2 are symmetric positive semidefinite,
 
                  C1, C2, A1, A2 are assumed to be constant symmetric matrices,
                  and b, k are constants.
*/
#include <iostream>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

std::shared_ptr<ndarray<int,1>>    nint(const std::vector<int> &X)    { return new_array_ptr<int>(X); }
std::shared_ptr<ndarray<double,1>> ndou(const std::vector<double> &X) { return new_array_ptr<double>(X); }

int main(int argc, char ** argv)
{
    // Sample data in sparse, symmetric triplet format
    std::vector<int>    C1_k = {0, 2};
    std::vector<int>    C1_l = {0, 2};
    std::vector<double> C1_v = {1, 6};
    std::vector<int>    A1_k = {0, 2, 0, 2};
    std::vector<int>    A1_l = {0, 0, 2, 2};
    std::vector<double> A1_v = {1, 1, 1, 2};
    std::vector<int>    C2_k = {0, 1, 0, 1, 2};
    std::vector<int>    C2_l = {0, 0, 1, 1, 2};
    std::vector<double> C2_v = {1, -3, -3, 2, 1};
    std::vector<int>    A2_k = {1, 0, 1, 3};
    std::vector<int>    A2_l = {0, 1, 1, 3};
    std::vector<double> A2_v = {1, 1, -1, -3};
    double b = 23;
    double k = -3;

    // Convert input data into Fusion sparse matrices
    auto C1 = Matrix::sparse(3, 3, nint(C1_k), nint(C1_l), ndou(C1_v));
    auto C2 = Matrix::sparse(4, 4, nint(C2_k), nint(C2_l), ndou(C2_v));
    auto A1 = Matrix::sparse(3, 3, nint(A1_k), nint(A1_l), ndou(A1_v));
    auto A2 = Matrix::sparse(4, 4, nint(A2_k), nint(A2_l), ndou(A2_v));

    // Create model
    Model::t M = new Model("sdo2"); auto _M = finally([&]() { M->dispose(); });

    // Two semidefinite variables
    auto X1 = M->variable(Domain::inPSDCone(3));
    auto X2 = M->variable(Domain::inPSDCone(4));

    // Objective
    M->objective(ObjectiveSense::Minimize, Expr::add(Expr::dot(C1,X1), Expr::dot(C2,X2)));

    // Equality constraint
    M->constraint(Expr::add(Expr::dot(A1,X1), Expr::dot(A2,X2)), Domain::equalsTo(b));

    // Inequality constraint
    M->constraint(X2->index(nint({0,1})), Domain::lessThan(k));

    // Solve
    M->setLogHandler([ = ](const std::string & msg) { std::cout << msg << std::flush; } );
    M->solve();

    // Retrieve solution
    std::cout << "Solution (vectorized) : " << std::endl;
    std::cout << "  X1 = " << *(X1->level()) << std::endl;
    std::cout << "  X2 = " << *(X2->level()) << std::endl;

    return 0;
}

sdo3.cc

Listing 16.36 sdo3.cc Click here to download.
/*
  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File :      sdo3.cc

  Purpose :   Solves the semidefinite problem:

                 min   tr(X_1) + ... + tr(X_n)
                 st.   <A_11,X_1> + ... + <A_1n,X_n> >= b_1
                       ...
                       <A_k1,X_1> + ... + <A_kn,X_n> >= b_k
                
                 where X_i are symmetric positive semidefinite of dimension d,

                 A_ji are constant symmetric matrices and b_i are constant.

              This example is to demonstrate creating and using 
              many matrix variables of the same dimension.
*/

#include <iostream>
#include <random>
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

// A helper function which returns a slice corresponding to j-th variable
Variable::t slice(Variable::t X, int d, int j) {
    return 
        X->slice(new_array_ptr<int,1>({j,0,0}), new_array_ptr<int,1>({j+1,d,d}))
         ->reshape(new_array_ptr<int,1>({d,d}));
}

int main(int argc, char ** argv)
{
    std::random_device rd;
    std::mt19937 e2(rd());
    std::uniform_real_distribution<> dist(0, 1);

    // Sample data 
    int n = 100, d = 4, k = 3;
    std::vector<double> b({9, 10, 11});
    std::vector< std::shared_ptr<ndarray<double,2>> > A;
    for(int i=0; i<k*n; i++) {
        auto Ai = std::make_shared<ndarray<double,2>>(shape(d,d));
        for(int s1=0; s1<d; s1++)
            for(int s2=0; s2<=s1; s2++)
                (*Ai)(s1,s2) = (*Ai)(s2,s1) = dist(e2);
        A.push_back(Ai);
    }

    // Create a model with n semidefinite variables od dimension d x d
    Model::t M = new Model("sdo3"); auto _M = finally([&]() { M->dispose(); });

    Variable::t X = M->variable(Domain::inPSDCone(d, n));

    // Pick indexes (j, s, s), j=0..n-1, s=0..d, of diagonal entries for the objective
    auto alldiag =  std::make_shared<ndarray<int,2>>(
            shape(d*n,3), 
            std::function<int(const shape_t<2> &)>([d](const shape_t<2> & p) { return p[1]==0 ? p[0]/d : p[0]%d; }));

    M->objective(ObjectiveSense::Minimize, Expr::sum( X->pick(alldiag) ));

    // Each constraint is a sum of inner products
    // Each semidefinite variable is a slice of X
    for(int i=0; i<k; i++) {
        std::vector<Expression::t> sumlist;
        for(int j=0; j<n ;j++)
            sumlist.push_back( Expr::dot(A[i*n+j], slice(X, d, j)) );

        M->constraint(Expr::add(new_array_ptr(sumlist)), Domain::greaterThan(b[i]));
    }

    // Solve
    M->setLogHandler([ = ](const std::string & msg) { std::cout << msg << std::flush; } );            // Add logging
    M->writeTask("sdo3.ptf");                // Save problem in readable format
    M->solve();

    // Get results. Each variable is a slice of X
    std::cout << "Contributing variables:" << std::endl;
    for(int j=0; j<n; j++) {
        auto Xj = *(slice(X, d, j)->level());
        double maxval = 0;
        for(int s=0; s<d*d; s++) maxval = std::max(maxval, Xj[s]);
        if (maxval>1e-6) {
            std::cout << "X" << j << "=" << std::endl;
            for(int s1=0; s1<d; s1++) {
                for(int s2=0; s2<d; s2++) std::cout << Xj[s1*d+s2] << "  ";
                std::cout << std::endl;
            }
        }
    }

    return 0;
}

sospoly.cc

Listing 16.37 sospoly.cc Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      sospoly.cc
//
//  Purpose:
//  Models the cone of nonnegative polynomials and nonnegative trigonometric
//  polynomials using Nesterov's framework  [1].
//
//  Given a set of coefficients (x0, x1, ..., xn), the functions model the
//  cone of nonnegative polynomials
//
//  P_m = { x \in R^{n+1} | x0 + x1*t + ... xn*t^n >= 0, forall t \in I }
//
//  where I can be the entire real axis, the semi-infinite interval [0,inf), or
//  a finite interval I = [a, b], respectively.
//
//  References:
//
//  [1] "Squared Functional Systems and Optimization Problems",
//      Y. Nesterov, in High Performance Optimization,
//      Kluwer Academic Publishers, 2000.
#include <string>
#include <iostream>
#include <cmath>
#include "fusion.h"
#include <cassert>

using namespace mosek::fusion;
using namespace monty;

std::shared_ptr<ndarray<int, 1>> range (int start, int stop, int step);
std::shared_ptr<ndarray<int, 1>> range (int start, int stop);
std::shared_ptr<ndarray<int, 1>> range (int stop);
std::shared_ptr<ndarray<double, 1>> range (double start, double stop, double step);
std::shared_ptr<ndarray<double, 1>> range (double start, double stop);
std::shared_ptr<ndarray<double, 1>> range (double stop);
template<typename T> std::ostream & operator<<(std::ostream & os, ndarray<T, 1> & a);

// Creates a Hankel matrix of dimension n+1, where
// H_lk = a if l+k=i, and 0 otherwise.
Matrix::t Hankel(int n, int i, double a)
{
  if (i < n + 1)
    return Matrix::sparse(n + 1, n + 1, range(i, -1, -1),      range(i + 1),      a);
  else
    return Matrix::sparse(n + 1, n + 1, range(n, i - n - 1, -1), range(i - n, n + 1), a);
}

// Models the cone of nonnegative polynomials on the real axis
void nn_inf(Model::t M, Variable::t x)
{
  int m = (int)x->getSize() - 1;
  int n = (m / 2); // degree of polynomial is 2*n

  // Setup variables
  Variable::t X = M->variable(Domain::inPSDCone(n + 1));

  // x_i = Tr H(n, i) * X  i=0,...,m
  for (int i = 0; i < m + 1; ++i)
    M->constraint( Expr::sub(x->index(i), Expr::dot(Hankel(n, i, 1.0), X)), Domain::equalsTo(0.0));
}

// Models the cone of nonnegative polynomials on the semi-infinite interval [0,inf)
void nn_semiinf(Model::t M, Variable::t x)
{
  int n = (int)x->getSize() - 1;
  int n1 = n / 2,
      n2 = (n - 1) / 2;

  // Setup variables
  Variable::t X1 = M->variable(Domain::inPSDCone(n1 + 1));
  Variable::t X2 = M->variable(Domain::inPSDCone(n2 + 1));

  // x_i = Tr H(n1, i) * X1 + Tr H(n2,i-1) * X2, i=0,...,n
  for (int i = 0; i < n + 1; ++i)
    M->constraint( Expr::sub(x->index(i),
                             Expr::add(Expr::dot(Hankel(n1, i, 1.0),  X1),
                                       Expr::dot(Hankel(n2, i - 1, 1.0), X2))), Domain::equalsTo(0.0) );
}

// Models the cone of nonnegative polynomials on the finite interval [a,b]
void nn_finite(Model::t M, Variable::t x, double a, double b)
{
  int m = (int)x->getSize() - 1;
  int n = m / 2;

  if (m == 2 * n)
  {
    Variable::t X1 = M->variable(Domain::inPSDCone(n + 1));
    Variable::t X2 = M->variable(Domain::inPSDCone(n));

    // x_i = Tr H(n,i)*X1 + (a+b)*Tr H(n-1,i-1) * X2 - a*b*Tr H(n-1,i)*X2 - Tr H(n-1,i-2)*X2, i=0,...,m
    for (int i = 0; i < m + 1; ++i)
      M->constraint( Expr::sub(x->index(i),
                               Expr::add(Expr::sub(Expr::dot(Hankel(n, i, 1.0),  X1), Expr::dot(Hankel(n - 1, i, a * b), X2)),
                                         Expr::sub(Expr::dot(Hankel(n - 1, i - 1, a + b), X2), Expr::dot(Hankel(n - 1, i - 2, 1.0),  X2)))),
                     Domain::equalsTo(0.0) );
  }
  else
  {
    Variable::t X1 = M->variable(Domain::inPSDCone(n + 1));
    Variable::t X2 = M->variable(Domain::inPSDCone(n + 1));

    // x_i = Tr H(n,i-1)*X1 - a*Tr H(n,i)*X1 + b*Tr H(n,i)*X2 - Tr H(n,i-1)*X2, i=0,...,m
    for (int i = 0; i < m + 1; ++i)
      M->constraint(Expr::sub(x->index(i),
                              Expr::add(Expr::sub(Expr::dot(Hankel(n, i - 1, 1.0),  X1), Expr::dot(Hankel(n, i, a), X1)),
                                        Expr::sub(Expr::dot(Hankel(n, i, b), X2), Expr::dot(Hankel(n, i - 1, 1.0),  X2)))),
                    Domain::equalsTo(0.0) );
  }
}

// returns variables u representing the derivative of
//  x(0) + x(1)*t + ... + x(n)*t^n,
// with u(0) = x(1), u(1) = 2*x(2), ..., u(n-1) = n*x(n).
Variable::t diff(Model::t M, Variable::t x)
{
  int n = (int)x->getSize() - 1;
  Variable::t u = M->variable(n, Domain::unbounded());
  Matrix::t mx = Matrix::dense(n, 1, range(1.0, n + 1.0));
  M->constraint(Expr::sub(u, Expr::mulElm(mx, x->slice(1, n + 1))), Domain::equalsTo(0.0));
  return u;
}

std::shared_ptr<ndarray<double, 1>> fitpoly(ndarray<double, 2> & data, int n)
{
  Model::t M = new Model("smooth poly"); auto _M = finally([&]() { M->dispose(); });

  int m = data.size(0);
  // create an (m x n+1) array, initialize each (i,j) entry with pow(data[i,0],j)
  std::shared_ptr<ndarray<double, 2>> A(new ndarray<double, 2>(shape_t<2>(m, n + 1), std::function<double(const shape_t<2>&)>([&](const shape_t<2>& p) { return std::pow(data(p[0], 0), p[1]); })));
  std::shared_ptr<ndarray<double, 1>> b(new ndarray<double, 1>(m, std::function<double(ptrdiff_t)>([&](ptrdiff_t i) { return data(i, 1); })));

  Variable::t x  = M->variable("x", n + 1, Domain::unbounded());
  Variable::t z  = M->variable("z", 1,   Domain::unbounded());
  Variable::t dx = diff(M, x);

  M->constraint(Expr::mul(A, x), Domain::equalsTo(b));

  // z - f'(t) >= 0, for all t \in [a, b]
  Variable::t ub = M->variable(n, Domain::unbounded());
  M->constraint(Expr::sub(ub,
                          Expr::vstack(Expr::sub(z, dx->index(0)), Expr::neg(dx->slice(1, n)))),
                Domain::equalsTo(0.0));
  nn_finite(M, ub, data(0, 0), data(data.size(0) - 1, 0));

  // f'(t) + z >= 0, for all t \in [a, b]
  Variable::t lb = M->variable(n, Domain::unbounded());
  M->constraint(Expr::sub(lb, Expr::vstack(Expr::add(z, dx->index(0)), dx->slice(1, n))),
                Domain::equalsTo(0.0));
  nn_finite(M, lb, data(0, 0), data(data.size(0) - 1, 0));

  M->objective(ObjectiveSense::Minimize, z);
  M->solve();
  return x->level();
}


int main(int argc, char ** argv)
{
  ndarray<double, 2> data
  { { -1.0, 1.0 },
    {  0.0, 0.0 },
    {  1.0, 1.0 } };

  auto x2 = fitpoly(data, 2);
  auto x4 = fitpoly(data, 4);
  auto x8 = fitpoly(data, 8);

  std::cout << "fitpoly(data,2) -> " << x2 << std::endl;
  std::cout << "fitpoly(data,4) -> " << x4 << std::endl;
  std::cout << "fitpoly(data,8) -> " << x8 << std::endl;
}


// Some utility functions to make things look nicer

template<typename T>
std::ostream & operator<<(std::ostream & os, ndarray<T, 1> & a)
{
  os << "[";
  for (auto v : a) os << " " << v;
  os << std::endl;
  return os;
}


std::shared_ptr<ndarray<int, 1>> range (int start, int stop, int step)
{
  range_t<int> r(start, stop, step);
  return std::shared_ptr<ndarray<int, 1>>(new ndarray<int, 1>(shape_t<1>(r.size()), iterable(r)));
}
std::shared_ptr<ndarray<int, 1>> range (int start, int stop) { return range(start, stop, 1); }
std::shared_ptr<ndarray<int, 1>> range (int stop) { return range(0, stop, 1); }

std::shared_ptr<ndarray<double, 1>> range (double start, double stop, double step)
{
  ptrdiff_t num = ( start < stop && step > 0 ? floor((stop - start) / step) :
                    start > step && step < 0 ? floor((stop - start) / step) :
                    0 );
  auto res = new_array_ptr<double, 1>(num);
  double v  = start;
  auto   it = res->begin();
  for (ptrdiff_t i = 0; i < num; ++i, v += step, ++it)
    * it = v;
  return res;
}
std::shared_ptr<ndarray<double, 1>> range (double start, double stop) { return range(start, stop, 1.0); }
std::shared_ptr<ndarray<double, 1>> range (double stop) { return range(0.0, stop, 1.0); }

sudoku.cc

Listing 16.38 sudoku.cc Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      sudoku.cc
//
// Purpose:  A MILP-based SUDOKU solver
//
//

#include <iostream>
#include <sstream>
#include <cmath>

#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

void print_solution(int n, Variable::t X)
{
  using namespace std;

  cout << "\n";
  int m( std::sqrt(n) );
  for (int i = 0; i < n; i++)
  {
    stringstream ss;

    for (int j = 0; j < n; j++)
    {
      if (j % m == 0) ss << " |";

      for (int k = 0; k < n; k++)
      {
        auto x = X->index(new_array_ptr<int, 1>({i, j, k}))->level();
        if ( (*x)[0] > 0.5 )
        {
          ss << " " << (k + 1);
          break;
        }
      }
    }
    cout << ss.str() << " |";

    cout << "\n";
    if ((i + 1) % m == 0)
      cout << "\n";
  }
}

int main(int argc, char ** argv)
{

  int m = 3;
  int n = m * m;

  //fixed cells in human readable (i.e. 1-based) format
  auto hr_fixed = new_array_ptr<int, 2>(
  { {1, 5, 4},
    {2, 2, 5}, {2, 3, 8}, {2, 6, 3},
    {3, 2, 1}, {3, 4, 2}, {3, 5, 8}, {3, 7, 9},
    {4, 2, 7}, {4, 3, 3}, {4, 4, 1}, {4, 7, 8}, {4, 8, 4},
    {6, 2, 4}, {6, 3, 1}, {6, 6, 9}, {6, 7, 2}, {6, 8, 7},
    {7, 3, 4}, {7, 5, 6}, {7, 6, 5}, {7, 8, 8},
    {8, 4, 4}, {8, 7, 1}, {8, 8, 6},
    {9, 5, 9}
  }
  );

  int nfixed = hr_fixed->size() / m;

  Model::t M = new Model("SUDOKU"); auto _M = finally([&]() { M->dispose(); });

  M->setLogHandler([ = ](const std::string & msg) { std::cout << msg << std::flush; });

  Variable::t X = M->variable("X", new_array_ptr<int, 1>({n, n, n}), Domain::binary());

  //each value only once per dimension
  for (int d = 0; d < m; d++)
    M->constraint( Expr::sum(X, d), Domain::equalsTo(1.) );

  //each number must appear only once in a block
  for (int k = 0; k < n ; k++)
    for (int i = 0; i < m ; i++)
      for (int j = 0; j < m ; j++)
        M->constraint( Expr::sum( X->slice( new_array_ptr<int, 1>({i * m, j * m, k}),
                                            new_array_ptr<int, 1>({(i + 1)*m, (j + 1)*m, k + 1}) ) ),
                       Domain::equalsTo(1.) );

  auto fixed = std::shared_ptr< ndarray<int, 2> >( new ndarray<int, 2>( shape(nfixed, 3) ) );

  for (int i = 0; i < nfixed; i++)
    for (int d = 0; d < m; d++)
      (*fixed)(i, d) =  (*hr_fixed)(i, d) - 1;

  M->constraint( X->pick( fixed ) , Domain::equalsTo(1.0) ) ;

  M->solve();

  //print the solution, if any...
  if ( M->getPrimalSolutionStatus() == SolutionStatus::Optimal )
    print_solution(n, X);
  else
    std::cout << "No solution found!\n";

  return 0;
}

total_variation.cc

Listing 16.39 total_variation.cc Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      total_variation.cc
//
// Purpose:   Demonstrates how to solve a total
//            variation problem using the Fusion API.

#include <iostream>
#include <vector>
#include <random>

#include "monty.h"
#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;

Model::t total_var(int n, int m) {
  Model::t M = new Model("TV");

  Variable::t u = M->variable("u", new_array_ptr<int, 1>({n + 1, m + 1}), Domain::inRange(0.0, 1.0));
  Variable::t t = M->variable("t", new_array_ptr<int, 1>({n, m}), Domain::unbounded());

  // In this example we define sigma and the input image f as parameters
  // to demonstrate how to solve the same model with many data variants.
  // Of course they could simply be passed as ordinary arrays if that is not needed.
  Parameter::t sigma = M->parameter("sigma");
  Parameter::t f = M->parameter("f", n, m);

  Variable::t ucore = u->slice(new_array_ptr<int, 1>({0, 0}), new_array_ptr<int, 1>({n, m}));

  Expression::t deltax = Expr::sub( u->slice( new_array_ptr<int, 1>({1, 0}), new_array_ptr<int, 1>({n + 1, m}) ), ucore );
  Expression::t deltay = Expr::sub( u->slice( new_array_ptr<int, 1>({0, 1}), new_array_ptr<int, 1>({n, m + 1}) ), ucore );

  M->constraint( Expr::stack(2, t, deltax, deltay), Domain::inQCone()->axis(2) );

  M->constraint( Expr::vstack(sigma, Expr::flatten( Expr::sub(f,  ucore) ) ),
                 Domain::inQCone() );

  M->objective( ObjectiveSense::Minimize, Expr::sum(t) );

  return M;
}

int main(int argc, char ** argv)
{
  std::normal_distribution<double>       ndistr(0., 1.);
  std::mt19937 engine(0);

  int n = 100;
  int m = 200;
  std::vector<double> sigmas({ 0.0004, 0.0005, 0.0006 });

  // Create a parametrized model with given shape
  Model::t M = total_var(n, m);
  Parameter::t sigma = M->getParameter("sigma");
  Parameter::t f     = M->getParameter("f");
  Variable::t  ucore = M->getVariable("u")->slice(new_array_ptr<int, 1>({0, 0}), new_array_ptr<int, 1>({n, m}));

  // Example: Linear signal with Gaussian noise    
  std::vector<std::vector<double>> signal(n, std::vector<double>(m));
  std::vector<std::vector<double>> noise(n, std::vector<double>(m));
  std::vector<std::vector<double>> fVal(n, std::vector<double>(m));
  std::vector<std::vector<double>> sol(n, std::vector<double>(m));

  for(int i=0; i<n; i++) for(int j=0; j<m; j++) {
    signal[i][j] = 1.0*(i+j)/(n+m);
    noise[i][j] = ndistr(engine) * 0.08;
    fVal[i][j] = std::max( std::min(1.0, signal[i][j] + noise[i][j]), .0 );
  }
 
  // Set value for f
  f->setValue(new_array_ptr(fVal));

  for(int iter=0; iter<3; iter++) {
    // Set new value for sigma and solve
    sigma->setValue(sigmas[iter]*n*m);

    M->solve();

    // Retrieve solution from ucore
    auto ucoreLev = *(ucore->level());
    for(int i=0; i<n; i++) for(int j=0; j<m; j++)
      sol[i][j] = ucoreLev[i*n+m];

    // Use the solution
    // ...

    std::cout << "rel_sigma = " << sigmas[iter] << " total_var = " << M->primalObjValue() << std::endl;
  }

  M->dispose();

  return 0;
}

tsp.cc

Listing 16.40 tsp.cc Click here to download.
//
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      tsp.cc
//
//  Purpose: Demonstrates a simple technique to the TSP
//           usign the Fusion API.
//
#include <iostream>
#include <list>
#include <vector>

#include "fusion.h"

using namespace mosek::fusion;
using namespace monty;


void tsp(int n, Matrix::t A, Matrix::t C, bool remove_1_hop_loops, bool remove_2_hop_loops)
{
  Model::t M = new Model();
  auto M_ = finally([&]() { M->dispose(); });

  auto x = M->variable("x", Set::make(n,n), Domain::binary());

  M->constraint(Expr::sum(x, 0), Domain::equalsTo(1.0));
  M->constraint(Expr::sum(x, 1), Domain::equalsTo(1.0));
  M->constraint(x, Domain::lessThan( A ));

  M->objective(ObjectiveSense::Minimize, Expr::dot(C, x));

  if (remove_1_hop_loops)
    M->constraint(x->diag(), Domain::equalsTo(0.));

  if (remove_2_hop_loops)
    M->constraint(Expr::add(x, x->transpose()), Domain::lessThan(1.0));

  int it = 0;
  while (true)
  {
    M->solve();
    it++;

    typedef std::vector< std::tuple<int, int> > cycle_t;
    std::list< cycle_t > cycles;

    auto xlevel = x->level();

    for (int i = 0; i < n; i++)
      for (int j = 0; j < n; j++)
      {

        if ( (*xlevel)[i * n + j] <= 0.5 )
          continue;

        bool found = false;
        for (auto && c : cycles)
        {
          for (auto && cc : c)
          {
            if ( i == std::get<0>(cc) || i == std::get<1>(cc) ||
                 j == std::get<0>(cc) || j == std::get<1>(cc) )
            {
              c.push_back( std::make_tuple(i, j) );
              found = true;
              break;
            }
          }
          if (found) break;
        }

        if (!found)
          cycles.push_back(cycle_t(1, std::make_tuple(i, j)));
      }

    std::cout << "Iteration " << it << "\n";
    for (auto c : cycles) {
      for (auto cc : c)
        std::cout << "(" << std::get<0>(cc) << "," << std::get<1>(cc) << ") ";
      std::cout << "\n";
    }

    if (cycles.size() == 1) break;

    for (auto c : cycles)
    {
      int csize = c.size();

      auto tmp = std::shared_ptr<monty::ndarray<int, 2> >(new  ndarray<int, 2>( shape(csize, 2)) );
      for (auto i = 0; i < csize; ++i)
      {
        (*tmp)(i, 0) = std::get<0>(c[i]);
        (*tmp)(i, 1) = std::get<1>(c[i]);
      }

      M->constraint(Expr::sum(x->pick(tmp)), Domain::lessThan( 1.0 * csize - 1 ));
    }
  }
  try {
    auto xlevel = x->level();
    std::cout << "Solution\n";
    for (int i = 0; i < n; i++) {
      for (int j = 0; j < n; j++)
        std::cout << (int) (*xlevel)(i * n + j);
      std::cout << "\n";
    }
  } catch (...) {}
}

int main()
{
  auto A_i = new_array_ptr<int, 1>({0, 1, 2, 3, 1, 0, 2, 0});
  auto A_j = new_array_ptr<int, 1>({1, 2, 3, 0, 0, 2, 1, 3});

  auto C_v = new_array_ptr<double, 1>({1., 1., 1., 1., 0.1, 0.1, 0.1, 0.1});

  int n = 4;
  tsp(n, Matrix::sparse(n, n, A_i, A_j, 1.), Matrix::sparse(n, n, A_i, A_j, C_v), true, false);
  tsp(n, Matrix::sparse(n, n, A_i, A_j, 1.), Matrix::sparse(n, n, A_i, A_j, C_v), true, true);

  return 0;
}