16 List of examples¶
List of examples shipped in the distribution of Fusion API for C++:
File |
Description |
|---|---|
Demonstrates a traffic network problem as a conic quadratic problem (CQO) |
|
A portfolio choice model |
|
A small bakery revenue maximization linear problem |
|
Shows how to break a long-running task |
|
An example of data/progress callback |
|
A simple conic exponential problem |
|
A simple conic quadratic problem |
|
Solving Stigler’s Nutrition model |
|
A simple problem with disjunctive constraints (DJC) |
|
Shows how to access the dual solution |
|
Linear regression with elastic net. Demonstrates model parametrization. |
|
Demonstrates a small one-facility location problem (CQO) |
|
A simple geometric program (GP) in conic form |
|
A simple linear problem |
|
Implements logistic regression and simple log-sum-exp (CEO) |
|
Computes the Lowner-John inner and outer ellipsoidal approximations of a polytope (SDO, Power Cone) |
|
Demonstrates how to solve the multi-processor scheduling problem and input an integer feasible point (MIP) |
|
A simple mixed-integer conic problem |
|
A simple mixed-integer linear problem |
|
A simple mixed-integer linear problem with an initial guess |
|
Solves the nearest correlation matrix problem (SDO, CQO) |
|
Uses MOSEK OptServer to solve an optimization problem synchronously |
|
Demonstrates parallel optimization using a batch method in MOSEK |
|
Shows how to set optimizer parameters and read information items |
|
Shows how to obtain and analyze a primal infeasibility certificate |
|
Portfolio optimization - basic Markowitz model |
|
Portfolio optimization - efficient frontier |
|
Portfolio optimization - market impact costs |
|
Portfolio optimization - transaction costs |
|
Portfolio optimization - cardinality constraints |
|
Portfolio optimization - factor model |
|
A simple power cone problem |
|
Implements a simple soft-margin Support Vector Machine (CQO) |
|
Demonstrate how to use SDP to solve convex relaxation of a mixed-integer QCQO problem |
|
Demonstrate how to modify and re-optimize a linear problem |
|
Demonstrates proper response handling |
|
A simple semidefinite problem with one matrix variable and a quadratic cone |
|
A simple semidefinite problem with two matrix variables |
|
A simple semidefinite problem with many matrix variables of the same dimension |
|
Models the cone of nonnegative polynomials and nonnegative trigonometric polynomials using Nesterov’s framework |
|
A SUDOKU solver (MIP) |
|
Demonstrates how to solve a total variation problem (CQO) |
|
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
//
// 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
//
// 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
//
// 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
//
// 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
/*
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
//
// 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
//
// 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
// 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
////
// 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
//
// 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
// 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
////
// 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
//
// 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
////
// 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
// 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
/*
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
//
// 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 these lines to switch off feeding in the initial LPT solution
x->setLevel(init);
M->setSolverParam("mioConstructSol", "on");
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
//
// 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
//
// 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
//
// 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
/*
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, N, 1) );
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
//
// 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
////
// 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
/*
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
////
// 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
/*
File : portfolio_1_basic.cc
Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
Description :
Implements a basic portfolio optimization model.
*/
#include <iostream>
#include <sstream>
#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();
// Check if the solution is an optimal point
SolutionStatus solsta = M->getPrimalSolutionStatus();
if (solsta != SolutionStatus::Optimal)
{
// See https://docs.mosek.com/latest/cxxfusion/accessing-solution.html about handling solution statuses.
std::ostringstream oss;
oss << "Unexpected solution status: " << solsta << std::endl;
throw SolutionError(oss.str());
}
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
/*
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();
// Check if the solution is an optimal point
SolutionStatus solsta = M->getPrimalSolutionStatus();
if (solsta != SolutionStatus::Optimal)
{
// See https://docs.mosek.com/latest/cxxfusion/accessing-solution.html about handling solution statuses.
std::ostringstream oss;
oss << "Unexpected solution status: " << solsta << std::endl;
throw SolutionError(oss.str());
}
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
/*
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();
// Check if the solution is an optimal point
SolutionStatus solsta = M->getPrimalSolutionStatus();
if (solsta != SolutionStatus::Optimal)
{
// See https://docs.mosek.com/latest/cxxfusion/accessing-solution.html about handling solution statuses.
std::ostringstream oss;
oss << "Unexpected solution status: " << solsta << std::endl;
throw SolutionError(oss.str());
}
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
/*
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();
// Check if the solution is an optimal point
SolutionStatus solsta = M->getPrimalSolutionStatus();
if (solsta != SolutionStatus::Optimal)
{
// See https://docs.mosek.com/latest/cxxfusion/accessing-solution.html about handling solution statuses.
std::ostringstream oss;
oss << "Unexpected solution status: " << solsta << std::endl;
throw SolutionError(oss.str());
}
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
/*
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();
// Check if the solution is an optimal point
SolutionStatus solsta = M->getPrimalSolutionStatus();
if (solsta != SolutionStatus::Optimal)
{
// See https://docs.mosek.com/latest/cxxfusion/accessing-solution.html about handling solution statuses.
std::ostringstream oss;
oss << "Unexpected solution status: " << solsta << std::endl;
throw SolutionError(oss.str());
}
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
/*
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();
// Check if the solution is an optimal point
SolutionStatus solsta = M->getPrimalSolutionStatus();
if (solsta != SolutionStatus::Optimal)
{
// See https://docs.mosek.com/latest/cxxfusion/accessing-solution.html about handling solution statuses.
std::ostringstream oss;
oss << "Unexpected solution status: " << solsta << std::endl;
throw SolutionError(oss.str());
}
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
//
// 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
//
// 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
/*
* 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
/**
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
/*
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
//
// 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
/*
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
/*
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
////
// 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
//
// 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
//
// 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
//
// 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;
}
