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

djc1.cc

A simple problem with disjunctive constraints (DJC)

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

nearestcorr.cc

Solves the nearest correlation matrix problem (SDO, CQO)

opt_server_sync.cc

Uses MOSEK OptServer to solve an optimization problem synchronously

parallel.cc

Demonstrates parallel optimization using a batch method in MOSEK

parameters.cc

Shows how to set optimizer parameters and read information items

pinfeas.cc

Shows how to obtain and analyze a primal infeasibility certificate

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

portfolio_6_factor.cc

Portfolio optimization - factor model

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

djc1.cc

Listing 16.9 djc1.cc Click here to download.
////
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      djc1.cc
//
//   Purpose: Demonstrates how to solve the problem with two disjunctions:
//
//      minimize    2x0 + x1 + 3x2 + x3
//      subject to   x0 + x1 + x2 + x3 >= -10
//                  (x0-2x1<=-1 and x2=x3=0) or (x2-3x3<=-2 and x1=x2=0)
//                  x0=2.5 or x1=2.5 or x2=2.5 or x3=2.5
////
#include <iostream>
#include "fusion.h"
using namespace mosek::fusion;
using namespace monty;

std::shared_ptr<ndarray<double,1>> dblarray(std::initializer_list<double> x) { 
  return new_array_ptr<double,1>(x); 
}

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

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

  // First disjunctive constraint
  M->disjunction( DJC::AND( DJC::term(Expr::dot(dblarray({1,-2,0,0}), x), Domain::lessThan(-1)), // x0 - 2x1 <= -1  
                            DJC::term(x->slice(2, 4), Domain::equalsTo(0)) ),                    // x2 = x3 = 0
                  DJC::AND( DJC::term(Expr::dot(dblarray({0,0,1,-3}), x), Domain::lessThan(-2)), // x2 - 3x3 <= -2
                            DJC::term(x->slice(0, 2), Domain::equalsTo(0)) ) );                  // x0 = x1 = 0

  // Second disjunctive constraint
  // Disjunctive constraint from an array of terms reading x_i = 2.5 for i = 0,1,2,3
  M->disjunction(std::make_shared<ndarray<Term::t,1>>(shape(4), [x](int i) { return DJC::term(x->index(i), Domain::equalsTo(2.5)); }));

  // The linear constraint
  M->constraint(Expr::sum(x), Domain::greaterThan(-10));

  // Objective
  M->objective(ObjectiveSense::Minimize, Expr::dot(dblarray({2,1,3,1}), x));

  // Useful for debugging
  M->writeTask("djc1.ptf");
  M->setLogHandler([ = ](const std::string & msg) { std::cout << msg << std::flush; } );

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

  // Get the solution values
  if (M->getPrimalSolutionStatus() == SolutionStatus::Optimal) {
    auto sol = x->level();
    std::cout << "[x0,x1,x2,x3] = " << (*sol) << std::endl;
  }
  else {
    std::cout << "Another solution status" << std::endl;
  }
}

duality.cc

Listing 16.10 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.11 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.12 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.13 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.14 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.15 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->writeTask("dump.ptf");
  M->solve();
  std::cout << *(theta->level());
  M->dispose();

  return 0;
}

lownerjohn_ellipsoid.cc

Listing 16.16 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;

/**
 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 geometric mean cone

   t <= (Z11*Z22*...*Znn)^{1/n}
*/
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));

  // (Z11*Z22*...*Znn) >= t^n
  M->constraint(Expr::vstack(DZ->diag(), t), Domain::inPGeoMeanCone());

  // 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 polygon 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.17 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.18 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.19 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.20 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 );

  // Request constructing the solution from integer variable values
  M->setSolverParam("mioConstructSol", "on");

  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 << "Construct solution utilization: " << constr << std::endl;
  std::cout << "Construct solution objective: " << constrVal << std::endl;

}

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 addr [certpath]" << std::endl;
        exit(0);
    }

    std::string serveraddr(argv[1]);
    std::string tlscert(argc==3 ? argv[2] : "");

    // 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(serveraddr);

    // Path to certificate, if any
    M->setSolverParam("remoteTlsCertPath", tlscert);

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

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

    return 0;
}

parallel.cc

Listing 16.23 parallel.cc Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      parallel.cc
//
//  Purpose: Demonstrates parallel optimization using solveBatch()
////

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

Model::t makeToyParameterizedModel() 
{
  auto M = new Model();
  auto x = M->variable("x",3);
  auto p = M->parameter("p");
  M->objective(ObjectiveSense::Maximize, Expr::sum(x));
  M->constraint(Expr::vstack(p,x), Domain::inQCone());
  return M;
}

/** Example of how to use Model.solveBatch()
*/
int main(int argc, char ** argv)
{
  // Choose some sample parameters
  int n = 10;                 // Number of models to optimize
  int threadpoolsize = 4;     // Total number of threads available
  int threadspermodel = 1;    // Number of threads per each model

  // Create a toy model for this example
  auto M = makeToyParameterizedModel();

  // Set up n copies of the model with different data
  auto models = std::make_shared<ndarray<Model::t,1>>(shape(n));

  for(int i = 0; i < n ; i++)
  {
    (*models)[i] = M->clone();
    (*models)[i]->getParameter("p")->setValue(i+1);
    // We can set the number of threads individually per model
    (*models)[i]->setSolverParam("numThreads", threadspermodel);
  }

  // Solve all models in parallel
  auto status = Model::solveBatch(false,         // No race
                                  -1.0,          // No time limit
                                  threadpoolsize,
                                  models);       // Array of Models to solve

  // Access the solutions
  for(int i = 0; i < n; i++) 
    if ((*status)[i] == SolverStatus::OK)
      std::cout << "Model "            <<  i << ":  "
                << "  Status "          <<  (*status)[i]
                << "  Solution Status " <<  (*models)[i]->getPrimalSolutionStatus()
                << "  Objective "       <<  (*models)[i]->primalObjValue()
                << "  Time "            <<  (*models)[i]->getSolverDoubleInfo("optimizerTime") << std::endl;
    else
      std::cout << "Model "           <<  i << ": not solved" << std::endl;
}

parameters.cc

Listing 16.24 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;
}

pinfeas.cc

Listing 16.25 pinfeas.cc Click here to download.
////
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      pinfeas.cc
//
//
//   Purpose: Demonstrates how to fetch a primal infeasibility certificate
//            for a linear problem
////
#include <iostream>
#include "fusion.h"
using namespace mosek::fusion;
using namespace monty;

//Analyzes and prints infeasibility certificate for a single object,
//which can be a variable or constraint
static void analyzeCertificate(std::string name,                              // name of the analyzed object
                               long size,                                     // size of the object
                               std::shared_ptr<ndarray<double, 1>> duals,     // actual dual values
                               double eps)                                    // tolerance determining when a dual value is considered important
{
  for(int i = 0; i < size; i++) {
    if (abs((*duals)[i]) > eps) 
      std::cout << name << "[" << i << "],   dual = " << (*duals)[i] << std::endl;
  }
}


int main(int argc, char ** argv)
{
  // Construct the sample model from the example in the manual
  auto sMat = Matrix::sparse(3, 7, new_array_ptr<int,1>({0,0,1,1,2,2,2}),
                                   new_array_ptr<int,1>({0,1,2,3,4,5,6}),
                                   new_array_ptr<double,1>({1,1,1,1,1,1,1}));
  auto sBound = new_array_ptr<double,1>({200, 1000, 1000});
  auto dMat = Matrix::sparse(4, 7, new_array_ptr<int,1>({0,0,1,2,2,3,3}),
                                   new_array_ptr<int,1>({0,4,1,2,5,3,6}),
                                   new_array_ptr<double,1>({1,1,1,1,1,1,1}));
  auto dBound = new_array_ptr<double,1>({1100, 200, 500, 500});
  auto c = new_array_ptr<double,1>({1, 2, 5, 2, 1, 2, 1});

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

  Variable::t x = M->variable("x", 7, Domain::greaterThan(0));
  Constraint::t s = M->constraint("s", Expr::mul(sMat, x), Domain::lessThan(sBound));
  Constraint::t d = M->constraint("d", Expr::mul(dMat, x), Domain::equalsTo(dBound));
  M->objective(ObjectiveSense::Minimize, Expr::dot(c,x));

  // Useful for debugging
  M->writeTask("pinfeas.ptf");
  M->setLogHandler([ = ](const std::string & msg) { std::cout << msg << std::flush; } );

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

  // Check problem status
  if (M->getProblemStatus() == ProblemStatus::PrimalInfeasible) {
    // Set the tolerance at which we consider a dual value as essential
    double eps = 1e-7;

    // We want to retrieve infeasibility certificates
    M->acceptedSolutionStatus(AccSolutionStatus::Certificate);

    // Go through variable bounds
    std::cout << "Variable bounds important for infeasibility: " << std::endl;
    analyzeCertificate("x", x->getSize(), x->dual(), eps);

    // Go through constraint bounds
    std::cout << "Constraint bounds important for infeasibility: " << std::endl;
    analyzeCertificate("s", s->getSize(), s->dual(), eps);
    analyzeCertificate("d", d->getSize(), d->dual(), eps);
  }
  else {
    std::cout << "The problem is not primal infeasible, no certificate to show" << std::endl;
  }
}

portfolio_1_basic.cc

Listing 16.26 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      = 8;
  auto       w      = 59.0;
  auto       mu     = new_array_ptr<double, 1>({0.07197349, 0.15518171, 0.17535435, 0.0898094 , 0.42895777, 0.39291844, 0.32170722, 0.18378628});
  auto       x0     = new_array_ptr<double, 1>({8.0, 5.0, 3.0, 5.0, 2.0, 9.0, 3.0, 6.0});
  auto       gammas = new_array_ptr<double, 1>({36});
  auto       GT     = new_array_ptr<double, 2>({
        {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
        {0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
        {0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
        {0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
        {0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 },
        {0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187},
        {0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327},
        {0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 }
  });

  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.27 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      = 8;
  double     w      = 1.0;
  auto       mu     = new_array_ptr<double, 1>({0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379});
  auto       x0     = new_array_ptr<double, 1>({0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
  auto       GT     = new_array_ptr<double, 2>({
    {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
    {0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
    {0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
    {0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
    {0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 },
    {0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187},
    {0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327},
    {0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 }
  });

  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) << "std. dev." << 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) << sqrt(frontier_s[i]) << std::endl;

  return 0;
}

portfolio_3_impact.cc

Listing 16.28 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      = 8;
  double     w      = 1.0;
  auto       mu     = new_array_ptr<double, 1>( {0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379} );
  auto       x0     = new_array_ptr<double, 1>({0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
  auto       GT     = new_array_ptr<double, 2>({
    {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
    {0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
    {0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
    {0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
    {0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 },
    {0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187},
    {0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327},
    {0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 }
  });

  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), 0.01));
  double gamma = 0.36;
  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.29 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      = 8;
  double     w      = 1.0;
  auto       mu     = new_array_ptr<double, 1>( {0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379} );
  auto       x0     = new_array_ptr<double, 1>({0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
  auto       GT     = new_array_ptr<double, 2>({
    {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
    {0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
    {0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
    {0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
    {0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 },
    {0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187},
    {0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327},
    {0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 }
  });

  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 = 0.36;
  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;
  std::cout << "Expected return:" << dot(mu, x) << std::endl;
  return 0;
}

portfolio_5_card.cc

Listing 16.30 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      = 8;
  double     w      = 1.0;
  auto       mu     = new_array_ptr<double, 1>( {0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379} );
  auto       x0     = new_array_ptr<double, 1>({0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
  auto       GT     = new_array_ptr<double, 2>({
    {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
    {0.     , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
    {0.     , 0.     , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
    {0.     , 0.     , 0.     , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
    {0.     , 0.     , 0.     , 0.     , 0.36096, 0.12574, 0.10157, 0.0571 },
    {0.     , 0.     , 0.     , 0.     , 0.     , 0.21552, 0.05663, 0.06187},
    {0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.22514, 0.03327},
    {0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.     , 0.2202 }
  });
  auto       gamma  = 0.25;
  
  std::vector<int> kValues = { 1, 2, 3, 4, 5, 6, 7, 8 };

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

portfolio_6_factor.cc

Listing 16.31 portfolio_6_factor.cc Click here to download.
/*
  File : portfolio_6_factor.cc

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

  Description :
    Implements a portfolio optimization model using factor model.
*/

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

using namespace mosek::fusion;
using namespace monty;

template<class T, int N>
using farray = std::shared_ptr<monty::ndarray<T, N>>;

// These are for debug purposes
void print_farray1(const farray<double, 1> &a)
{
  for (int j = 0; j < a->shape[0]; ++j)
  {
    std::cout << (*a)(j) << ",\t";
  }    
  std::cout << std::endl;
}

void print_farray2(const farray<double, 2> &a)
{
  for (int i = 0; i < a->shape[0]; ++i)
  {
    for (int j = 0; j < a->shape[1]; ++j)
    {
      std::cout << (*a)(i, j) << ",\t";
    }    
    std::cout << std::endl;
  }
}

static double sum(const farray<double, 1> &x)
{
  double r = 0.0;
  for (auto v : *x) r += v;
  return r;
}

static double dot(const farray<double, 1> &x, const farray<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(const farray<double, 1> &x, const std::vector<double> &y)
{
  double r = 0.0;
  for (int i = 0; i < x->size(); ++i) r += (*x)[i] * y[i];
  return r;
}

static farray<double, 2> transpose(const farray<double, 2> &m) 
{
  int ni = m->shape[0];
  int nj = m->shape[1];
  auto mt = std::make_shared<ndarray<double, 2>>(shape(nj, ni), 0.0);

  for (int i = 0; i < ni; ++i)
  {
    for (int j = 0; j < nj; ++j)
    {
      (*mt)(j, i) = (*m)(i, j);
    }
  }
  return mt;
}

static farray<double, 1> vector_sqrt(const farray<double, 1> &m) 
{
  int ni = m->shape[0];
  farray<double, 1> sqrtm = std::make_shared<ndarray<double, 1>>(shape(ni));

  for (int i = 0; i < ni; ++i)
  {
    (*sqrtm)(i) = sqrt((*m)(i));
  }
  return sqrtm;
}

// Vectorize matrix (column-major order)
static farray<double, 1> mat_to_vec_c(const farray<double, 2> &m) 
{
  int ni = m->shape[0];
  int nj = m->shape[1];
  auto c = std::make_shared<ndarray<double, 1>>(shape(nj * ni));

  for (int j = 0; j < nj; ++j)
  {
    for (int i = 0; i < ni; ++i)
    {
      (*c)(j * ni + i) = (*m)(i, j);
    }
  }
  return c;
}

// Reshape vector to matrix (column-major order)
static farray<double, 2> vec_to_mat_c(const farray<double, 1> &c, int ni, int nj) 
{
  auto m = std::make_shared<ndarray<double, 2>>(shape(ni, nj));

  for (int j = 0; j < nj; ++j)
  {
    for (int i = 0; i < ni; ++i)
    {
      (*m)(i, j) = (*c)(j * ni + i);
    }
  }
  return m;
}

static farray<double, 2> cholesky(const farray<double, 2> &m) 
{
  int n = m->shape[0];
  auto vecs = mat_to_vec_c(m);
  mosek::LinAlg::potrf(MSK_UPLO_LO, n, vecs);
  auto s = vec_to_mat_c(vecs, n, n);
  
  // Zero out upper triangular part (LinAlg::potrf does not use it, original matrix values remain there)
  for (int i = 0; i < n; ++i)
  {
    for (int j = i+1; j < n; ++j)
    {
      (*s)(i, j) = 0.0;
    }
  }
  return s;
}

static farray<double, 2> matrix_mul(const farray<double, 2> &a, const farray<double, 2> &b) 
{
  int na = a->shape[0];
  int nb = b->shape[1];
  int k = b->shape[0];

  auto vecm = std::make_shared<ndarray<double, 1>>(shape(na * nb), 0.0);
  mosek::LinAlg::gemm(MSK_TRANSPOSE_NO, MSK_TRANSPOSE_NO, na, nb, k, 1.0, mat_to_vec_c(a), mat_to_vec_c(b), 1.0, vecm);
  auto m = vec_to_mat_c(vecm, na, nb);

  return m;
}

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

Input:
    n: Number of assets
    mu: An n dimmensional vector of expected returns
    G_factor_T: The factor (dense) part of the factorized risk
    theta: specific risk vector
    x0: Initial holdings
    w: Initial cash holding
    gamma: Maximum risk (=std. dev) accepted

Output:
    Optimal expected return and the optimal portfolio
*/
double FactorMarkowitz( 
  int                 n,
  farray<double, 1>   mu,
  farray<double, 2>   G_factor_T,
  farray<double, 1>   theta,
  farray<double, 1>   x0,
  double              w,
  double              gamma
)
{
  Model::t M = new Model("Factor 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(new_array_ptr<Expression::t, 1>({Expr::constTerm(gamma),
                                                                      Expr::mul(G_factor_T, x),
                                                                      Expr::mulElm(vector_sqrt(theta), 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      = 8;
  double     w      = 1.0;
  auto       mu     = new_array_ptr<double, 1>({0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379});
  auto       x0     = new_array_ptr<double, 1>({0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
  // Factor exposure matrix
  auto B = new_array_ptr<double, 2>({
    {0.4256, 0.1869},
    {0.2413, 0.3877},
    {0.2235, 0.3697},
    {0.1503, 0.4612},
    {1.5325, -0.2633},
    {1.2741, -0.2613},
    {0.6939, 0.2372},
    {0.5425, 0.2116}
  });

  // Factor covariance matrix
  auto S_F = new_array_ptr<double, 2>({
      {0.0620, 0.0577},
      {0.0577, 0.0908}
  });

  // Specific risk components
  auto theta = new_array_ptr<double, 1>({0.0720, 0.0508, 0.0377, 0.0394, 0.0663, 0.0224, 0.0417, 0.0459});

  auto P = cholesky(S_F);
  auto G_factor = matrix_mul(B, P);
  auto G_factor_T = transpose(G_factor);

  auto gammas = new_array_ptr<double, 1>({0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48});

  std::cout << std::endl
            << "-----------------------------------------------------------------------------------" << std::endl
            << "Markowitz portfolio optimization based on a factor model." << std::endl
            << "-----------------------------------------------------------------------------------" << std::endl;

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

  for (auto gamma : *gammas)
    std::cout << "Expected return: " << FactorMarkowitz(n, mu, G_factor_T, theta, x0, w, gamma) << " St deviation: " << gamma << std::endl;
  return 0;
}

pow1.cc

Listing 16.32 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.33 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.34 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.35 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.36 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.37 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.38 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.39 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.40 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.41 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.42 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.43 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;
}