16 List of examples

List of examples shipped in the distribution of Fusion API for .NET:

Table 16.1 List of distributed examples

File

Description

TrafficNetworkModel.cs

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

alan.cs

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

baker.cs

A small bakery revenue maximization linear problem

breaksolver.cs

Shows how to break a long-running task

callback.cs

An example of data/progress callback

ceo1.cs

A simple conic exponential problem

cqo1.cs

A simple conic quadratic problem

diet.cs

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

duality.cs

Shows how to access the dual solution

elastic.cs

Linear regression with elastic net. Demonstrates model parametrization.

facility_location.cs

Demonstrates a small one-facility location problem (CQO)

gp1.cs

A simple geometric program (GP) in conic form

lo1.cs

A simple linear problem

logistic.cs

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

lownerjohn_ellipsoid.cs

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

lpt.cs

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

mico1.cs

A simple mixed-integer conic problem

milo1.cs

A simple mixed-integer linear problem

mioinitsol.cs

A simple mixed-integer linear problem with an initial guess

modelLib.cs

Library of implementations of basic functions

nearestcorr.cs

Solves the nearest correlation matrix problem (SDO, CQO)

opt_server_sync.cs

Uses MOSEK OptServer to solve an optimization problem synchronously

parameters.cs

Shows how to set optimizer parameters and read information items

portfolio_1_basic.cs

Portfolio optimization - basic Markowitz model

portfolio_2_frontier.cs

Portfolio optimization - efficient frontier

portfolio_3_impact.cs

Portfolio optimization - market impact costs

portfolio_4_transcost.cs

Portfolio optimization - transaction costs

portfolio_5_card.cs

Portfolio optimization - cardinality constraints

pow1.cs

A simple power cone problem

primal_svm.cs

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

qcqp_sdo_relaxation.cs

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

reoptimization.cs

Demonstrate how to modify and re-optimize a linear problem

response.cs

Demonstrates proper response handling

sdo1.cs

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

sdo2.cs

A simple semidefinite problem with two matrix variables

sdo3.cs

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

sospoly.cs

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

sudoku.cs

A SUDOKU solver (MIP)

total_variation.cs

Demonstrates how to solve a total variation problem (CQO)

tsp.cs

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.cs

Listing 16.1 TrafficNetworkModel.cs Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      TrafficNetworkModel.cs
//
// Purpose:   Demonstrates a traffix 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.


using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class TrafficNetworkModel : Model
  {

    public static double[] solve
    ( int      numberOfNodes,
      int      source_idx,
      int      sink_idx,
      int[]    arc_i,
      int[]    arc_j,
      double[] arcSensitivity,
      double[] arcCapacity,
      double[] arcBaseTravelTime,
      double   T)

    {
      using (Model M = new Model("Traffic Network"))
      {
        int n = numberOfNodes;        
        int m = arc_j.Length;

        double[] n_ones = new double[m]; for (int i = 0; i < m; ++i) n_ones[i] = 1.0;
        var NxN = Set.Make(n, n);
        Matrix basetime =
          Matrix.Sparse(n, n, arc_i, arc_j, arcBaseTravelTime);
        Matrix e =
          Matrix.Sparse(n, n, arc_i, arc_j, n_ones);
        Matrix e_e =
          Matrix.Sparse(n, n,
                        new int[]  { sink_idx}, new int[] { source_idx},
                        new double[] { 1.0 });

        int[,] sparsity = new int[m+1,2]; for (int i = 0; i < m; ++i) { sparsity[i,0] = arc_i[i]; sparsity[i,1] = arc_j[i]; }
        sparsity[m,0] = sink_idx;
        sparsity[m,1] = source_idx;

        Console.WriteLine("sparsity dims : {0} x {1}",sparsity.GetLength(0), sparsity.GetLength(1));
        
        double[] cs_inv = new double[m];
        double[] s_inv  = new double[m];
        for (int i = 0; i < m; ++i) cs_inv[i] = 1.0 / (arcSensitivity[i] * arcCapacity[i]);
        for (int i = 0; i < m; ++i) s_inv[i]  = 1.0 / arcSensitivity[i];
        Matrix cs_inv_matrix = Matrix.Sparse(n, n, arc_i, arc_j, cs_inv);
        Matrix s_inv_matrix  = Matrix.Sparse(n, n, arc_i, arc_j, s_inv);

        Variable x = M.Variable("traffic_flow", NxN, Domain.GreaterThan(0.0).Sparse(sparsity));
        Variable d = M.Variable("d",            NxN, Domain.GreaterThan(0.0).Sparse(sparsity));
        Variable z = M.Variable("z",            NxN, Domain.GreaterThan(0.0).Sparse(sparsity));


        // Set the objective:       
        M.Objective("Average travel time", ObjectiveSense.Minimize,
                    Expr.Mul(1.0 / T, Expr.Add(Expr.Dot(basetime, x), Expr.Sum(d))));

        // Set up constraints
        // Constraint (1a)
        M.Constraint("(1a)", Expr.Hstack(d.Pick(sparsity),
                                         z.Pick(sparsity),
                                         x.Pick(sparsity)).Slice(new int[]{0,0},new int[]{m,3} ), Domain.InRotatedQCone());
        
        // 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();
        M.WriteTask("TrafficNetworkModel.task");

        return x.Level();
      }
    }

    public static void Main(string[] args)
    {
      int n             = 4;
      int[]    arc_i    = new int[]    {  0,    0,    2,    1,    2 };
      int[]    arc_j    = new int[]    {  1,    2,    1,    3,    3 };
      double[] arc_base = new double[] {  4.0,  1.0,  2.0,  1.0,  6.0 };
      double[] arc_cap  = new double[] { 10.0, 12.0, 20.0, 15.0, 10.0 };
      double[] arc_sens = new double[] {  0.1,  0.7,  0.9,  0.5,  0.1 };

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

      try {
          var flow = solve(n, source_idx, sink_idx,
                           arc_i, arc_j,
                           arc_sens,
                           arc_cap,
                           arc_base,
                           T);
          System.Console.WriteLine("Optimal flow: {0} {1}",arc_i.Length,flow.Length);
          
          for (int i = 0; i < arc_i.Length; ++i)
            System.Console.WriteLine("\tflow node {0}->node {1} = {2}",
                                     arc_i[i], arc_j[i], flow[arc_i[i] * n + arc_j[i]]);
      }
      catch (mosek.fusion.SolutionError e )
        {
          Console.WriteLine("SolutionError: {0}",e.ToString());
          Console.WriteLine("{0}",e.StackTrace);
          throw;
        }
    }
  }
}

alan.cs

Listing 16.2 alan.cs Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:     alan.cs
//
//  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.
//

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class alan
  {
    /////////////////////////////////////////////////////////////////////
    // Problem data.

    // Security names
    private static string[]
    securities = { "hardware", "software", "show-biz", "t-bills" };
    // Two examples of mean returns on securities
    private static double[]
    mean1      =  {        9.0,        7.0,       11.0,       5.0 };
    private static double[]
    mean2       = {        8.0,        9.0,       12.0,       7.0 };
    // Target mean return
    private static double
    target     = 10.0;

    // Factor of covariance matrix.
    private static Matrix U =
      Matrix.Dense(
        new double[,]
    { {  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 }
    });
    private static int numsec = securities.Length;
    public static void Main(String[] args)
    {
      // Create a parametrized model
      using (Model M = new Model("alan"))
      {
        Variable x = M.Variable("x",        numsec, Domain.GreaterThan(0.0));
        Variable t = M.Variable("variance", Domain.GreaterThan(0.0));
        M.Objective("minvar", ObjectiveSense.Minimize, t.AsExpr());

        // sum securities to 1.0
        M.Constraint("wealth",  Expr.Sum(x), Domain.EqualsTo(1.0));
        // define target expected return
        Parameter mean = M.Parameter(numsec);
        M.Constraint("dmean", Expr.Dot(mean, x), Domain.GreaterThan(target));

        M.Constraint("q",Expr.Vstack(Expr.ConstTerm(1, 0.5),
                                 t.AsExpr(),
                                 Expr.Mul(U, x)),
                     Domain.InRotatedQCone());

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

    public static void solve(Model M, Parameter mean, Variable x, double[] meanVal) {
      Console.Write("Solve with mean = ");
      for (int i = 0; i < numsec; ++i)
        Console.Write("  {0}", meanVal[i]);
      Console.WriteLine();

      mean.SetValue(meanVal);
      M.Solve();        

      double[] solx = x.Level();

      Console.Write("Solution = ");
      for (int i = 0; i < numsec; ++i)
        Console.Write(" {0}", solx[i]);
      Console.WriteLine("");        
    }
  }
}

baker.cs

Listing 16.3 baker.cs Click here to download.
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      baker.cs
//
// 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.

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class baker
  {
    private static string[] ingredientnames = { "Flour", "Sugar", "Butter" };
    private static double[] stock = { 150.0,   22.0,    25.0 };

    private static double[,] recipe_data =
    { { 3.0, 5.0 },
      { 1.0, 0.5 },
      { 1.2, 0.5 }
    };
    private static string[] productnames = { "Cakes", "Breads" };

    private static double[] revenue = { 4.0, 6.0 };

    public static void Main(string[] args)
    {
      Matrix recipe = Matrix.Dense(recipe_data);
      using (Model M = new Model("Recipe"))
      {
        // "production" defines the amount of each product to bake.
        Variable 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(Console.Out);

        // We solve and fetch the solution:
        M.Solve();
        double[] res = production.Level();
        Console.WriteLine("Solution:");
        for (int i = 0; i < res.Length; ++i)
        {
          Console.WriteLine(" Number of {0} : {1}", productnames[i], res[i]);
        }
        Console.WriteLine(" Revenue : ${0}",
                          res[0] * revenue[0] + res[1] * revenue[1]);
      }
    }
  }
}

breaksolver.cs

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


using System;
using System.Threading;
using mosek.fusion;
using System.Diagnostics;

namespace mosek.fusion.example
{
  public class breaksolver
  {
    private static void nshuffle(Random R, int[] a, int n)
    {
      for (int i = 0; i < n; ++i)
      {
        int idx = R.Next(i, a.Length);
        int tmp = a[i]; a[i] = a[idx]; a[idx] = tmp;
      }
    }

    public static void Main(string[] args)
    {
      long timeout = 5;

      int n = 200;    // number of binary variables
      int m = n / 5; // number of constraints
      int p = n / 5; // Each constraint picks p variables and requires that half of them are 1
      Random R = new Random(1234);
      using (Model M = new Model("SolveBinary"))
      {
        M.SetLogHandler(System.Console.Out);

        //Variable x = M.Variable("x", n, Domain.InRange(0,1));
        Variable x = M.Variable("x", n, Domain.Binary());
        M.Objective(ObjectiveSense.Minimize, Expr.Sum(x));
        //M.SetSolverParam("numThreads",1);

        int[] idxs  = new int[n]; for (int i = 0; i < n; ++i) idxs[i] = i;
        int[] cidxs = new int[p];

        for (var i = 0; i < m; ++i)
        {
          nshuffle(R, idxs, p);
          Array.Copy(idxs, cidxs, p);
          M.Constraint(Expr.Sum(x.Pick(cidxs)), Domain.EqualsTo(p / 2));
        }

        var T = new Thread(new ThreadStart(M.Solve));
        T.Start();

        Stopwatch w = new Stopwatch(); w.Start();
        while (true)
        {
          if (w.ElapsedMilliseconds > timeout * 1000)
          {
            Console.WriteLine("Solver terminated due to timeout!");
            M.BreakSolver();
            T.Join();
            break;
          }
          if (! T.IsAlive) {
            Console.WriteLine("Solver terminated before anything happened!");
            T.Join();
            break;
          }
        }

        Console.WriteLine("End.");
      }
    }
  }
}

callback.cs

Listing 16.5 callback.cs Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      callback.cs
//
//   Purpose:   To demonstrate how to use the progress
//              callback.
//
//              Use this script as follows:
//              callback psim
//              callback dsim
//              callback intpnt
//
//              The first argument tells which optimizer to use
//              i.e. psim is primal simplex, dsim is dual simplex
//              and intpnt is interior-point.
//

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  class myCallback : mosek.DataCallback
  {
    double maxtime;
    Model M;

    public myCallback( double maxtime_,
                       Model M_)
    {
      maxtime = maxtime_;
      M = M_;
    }

    public override int callback( callbackcode caller,
                                  double[]     douinf,
                                  int[]        intinf,
                                  long[]       lintinf )
    {
      double opttime = 0.0;
      int itrn;
      double pobj, dobj, stime;

      switch (caller)
      {
        case callbackcode.begin_intpnt:
          Console.WriteLine("Starting interior-point optimizer");
          break;
        case callbackcode.intpnt:
          itrn    = intinf[(int) iinfitem.intpnt_iter      ];
          pobj    = douinf[(int) dinfitem.intpnt_primal_obj];
          dobj    = douinf[(int) dinfitem.intpnt_dual_obj  ];
          stime   = douinf[(int) dinfitem.intpnt_time      ];
          opttime = douinf[(int) dinfitem.optimizer_time   ];

          Console.WriteLine("Iterations: {0,-3}", itrn);
          Console.WriteLine("  Elapsed: Time: {0,6:F2}({1:F2})", opttime, stime);
          Console.WriteLine("  Primal obj.: {0,-18:E6}  Dual obj.: {1,018:E6}e", pobj, dobj);
          break;
        case callbackcode.end_intpnt:
          Console.WriteLine("Interior-point optimizer finished.");
          break;
        case callbackcode.begin_primal_simplex:
          Console.WriteLine("Primal simplex optimizer started.");
          break;
        case callbackcode.update_primal_simplex:
          itrn    = intinf[(int) iinfitem.sim_primal_iter  ];
          pobj    = douinf[(int) dinfitem.sim_obj          ];
          stime   = douinf[(int) dinfitem.sim_time         ];
          opttime = douinf[(int) dinfitem.optimizer_time   ];

          Console.WriteLine("Iterations: {0,-3}}", itrn);
          Console.WriteLine("  Elapsed time: {0,6:F2}({1:F2})", opttime, stime);
          Console.WriteLine("  Obj.: {0,-18:E6}", pobj );
          break;
        case callbackcode.end_primal_simplex:
          Console.WriteLine("Primal simplex optimizer finished.");
          break;
        case callbackcode.begin_dual_simplex:
          Console.WriteLine("Dual simplex optimizer started.");
          break;
        case callbackcode.update_dual_simplex:
          itrn    = intinf[(int) iinfitem.sim_dual_iter    ];
          pobj    = douinf[(int) dinfitem.sim_obj          ];
          stime   = douinf[(int) dinfitem.sim_time         ];
          opttime = douinf[(int) dinfitem.optimizer_time   ];
          Console.WriteLine("Iterations: {0,-3}}", itrn);
          Console.WriteLine("  Elapsed time: {0,6:F2}({1:F2})", opttime, stime);
          Console.WriteLine("  Obj.: {0,-18:E6}", pobj );
          break;
        case callbackcode.end_dual_simplex:
          Console.WriteLine("Dual simplex optimizer finished.");
          break;
        case callbackcode.begin_bi:
          Console.WriteLine("Basis identification started.");
          break;
        case  callbackcode.end_bi:
          Console.WriteLine("Basis identification finished.");
          break;
        default:
          break;
      }
      if (opttime >= maxtime)
      {
        Console.WriteLine("MOSEK is spending too much time. Terminate it.");
        return 1;
      }
      return 0;
    }
  }

  class myProgressCallback : mosek.Progress
  {
    public override int progressCB( callbackcode caller )
    {
      //Handle caller code here
      return 0;
    }
  }

  public class callback
  {
    public static void Main(string[] args)
    {
      string slvr     = "intpnt";

      if (args.Length < 1)
      {
        Console.WriteLine("Usage: callback ( psim | dsim | intpnt )");
      }

      if (args.Length >= 1) slvr = args[0];

      // We create a large linear problem
      int n = 150;
      int m = 700;
      double[] A = new double[m * n];
      double[] b = new double[m];
      double[] c = new double[n];

      Random rnd = new Random();
      for (int i = 0; i < m * n; i++) A[i] = rnd.NextDouble();
      for (int i = 0; i < m; i++)   b[i] = rnd.NextDouble();
      for (int i = 0; i < n; i++)   c[i] = rnd.NextDouble();

      double maxtime = 0.07;

      Model M = new Model("callback");

      Variable 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");

      M.SetDataCallbackHandler( new myCallback(maxtime, M) );

      M.Solve();
    }
  }
}

ceo1.cs

Listing 16.6 ceo1.cs Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      ceo1.cs
//
//   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)
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class ceo1
  {
    public static void Main(string[] args)
    {
      using (Model M = new Model("ceo1"))
      {

        Variable 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 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
        double[] solx = x.Level();
        Console.WriteLine("x1,x2,x3 = {0}, {1}, {2}", solx[0], solx[1], solx[2]);

        // Get conic solution of expc
        double[] expclvl = expc.Level();
        double[] expcsn  = expc.Dual();
        
        Console.Write("expc levels = {0}", expclvl[0]);
        for (int i = 1; i < expclvl.Length; ++i)
          Console.Write(", {0}", expclvl[i]);
        Console.WriteLine();

        Console.Write("expc dual conic var levels = {0}", expcsn[0]);
        for (int i = 1; i < expcsn.Length; ++i)
          Console.Write(", {0}", expcsn[i]);
        Console.WriteLine();
      }
    }
  }
}

cqo1.cs

Listing 16.7 cqo1.cs Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      cqo1.cs
//
//   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              }
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class cqo1
  {
    public static void Main(string[] args)
    {
      using (Model M = new Model("cqo1"))
      {

        Variable x = M.Variable("x", 3, Domain.GreaterThan(0.0));
        Variable 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 z1 = Var.Vstack(y.Index(0),  x.Slice(0, 2));
        Variable z2 = Var.Vstack(y.Slice(1, 3), x.Index(2));

        // Create the constraint
        //      x[0] + x[1] + 2.0 x[2] = 1.0
        double[] aval = new double[] {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 qc1 = M.Constraint("qc1", z1.AsExpr(), Domain.InQCone());
        Constraint qc2 = M.Constraint("qc2", z2.AsExpr(), 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
        double[] solx = x.Level();
        double[] soly = y.Level();
        Console.WriteLine("x1,x2,x3 = {0}, {1}, {2}", solx[0], solx[1], solx[2]);
        Console.WriteLine("y1,y2,y3 = {0}, {1}, {2}", soly[0], soly[1], soly[2]);

        // Get conic solution of qc1
        double[] qc1lvl = qc1.Level();
        double[] qc1sn  = qc1.Dual();
        
        Console.Write("qc1 levels = {0}", qc1lvl[0]);
        for (int i = 1; i < qc1lvl.Length; ++i)
          Console.Write(", {0}", qc1lvl[i]);
        Console.WriteLine();

        Console.Write("qc1 dual conic var levels = {0}", qc1sn[0]);
        for (int i = 1; i < qc1sn.Length; ++i)
          Console.Write(", {0}", qc1sn[i]);
        Console.WriteLine();
      }
    }
  }
}

diet.cs

Listing 16.8 diet.cs Click here to download.
/*
  File : diet.cs

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

  Description :  Implements a basic linear model.
*/
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class diet : Model
  {
    public class DataException : FusionException
    {
      public DataException(string msg) : base(msg) {  }
    }

    public double[] dailyAllowance;
    public Matrix nutrientValue;

    private Variable   dailyPurchase;
    private Constraint dailyNutrients;

    public diet(string    name,
                string[]  foods,
                string[]  nutrients,
                double[]  daily_allowance,
                double[,] nutritive_value)
    : base(name)
    {
      bool finished = false;
      try
      {
        dailyAllowance = daily_allowance;
        nutrientValue = (Matrix.Dense(nutritive_value)).Transpose();

        int M = foods.Length;
        int N = nutrients.Length;

        if (dailyAllowance.Length != N)
          throw new DataException(
            "Length of daily_allowance does not match the number of nutrients");
        if (nutrientValue.NumColumns() != M)
          throw new DataException(
            "Number of rows in nutrient_value does not match the number of foods");
        if (nutrientValue.NumRows() != N)
          throw new DataException(
            "Number of columns in nutrient_value does not match the number of nutrients");

        dailyPurchase = Variable("Daily Purchase",
                                 Set.Make(foods),
                                 Domain.GreaterThan(0.0));
        dailyNutrients =
          Constraint("Nutrient Balance",
                     Expr.Mul(nutrientValue, dailyPurchase),
                     Domain.GreaterThan(dailyAllowance));
        Objective(ObjectiveSense.Minimize, Expr.Sum(dailyPurchase));
        finished = true;
      }
      finally
      {
        if (! finished)
          Dispose();
      }
    }

    public double[] getDailyPurchase()
    {
      return dailyPurchase.Level();
    }

    public double[] getDailyNutrients()
    {
      return dailyNutrients.Level();
    }


    /* Main class with data definitions */
    public static void Main(string[] argv)
    {
      string[] nutrient_unit = {
        "thousands",  "grams",        "grams",
        "milligrams", "thousand ius", "milligrams",
        "milligrams", "milligrams",   "milligrams"
      };
      string[] nutrients = {
        "calorie",    "protein",      "calcium",
        "iron",       "vitamin a",    "vitamin b1",
        "vitamin b2", "niacin",       "vitamin c"
      };
      string[] foods = {
        "wheat",        "cornmeal", "cannedmilk", "margarine", "cheese",
        "peanut butter", "lard",     "liver",      "porkroast", "salmon",
        "greenbeans",   "cabbage",  "onions",     "potatoes",  "spinach",
        "sweet potatos", "peaches",  "prunes",     "limabeans", "navybeans"
      };
      double[,] nutritive_value = {
        //  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,   0.6,     6,   55.8,     0.2,    0,      0,     0},  // margarine
        { 7.4,   448,  16.4,    19,   28.1,     0.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,    0.2,     0,      0.5,    5,     0},  // lard
        { 2.2,   333,   0.2,   139,  169.2,     6.4,   50.8,  316,   525},  // liver
        { 4.4,   249,   0.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
      double[] daily_allowance =
      {   3.0,    70.0, 0.8,    12.0,  5.0,     1.8,    2.7,   18.0,  75.0 };

      using (diet M = new diet("Stinglers Diet Problem",
                               foods,
                               nutrients,
                               daily_allowance,
                               nutritive_value))
      {
        M.Solve();

        try
        {
          Console.WriteLine("Solution (x 365):");
          double[] x = M.getDailyPurchase();
          for (int i = 0; i < x.Length; ++i)
            Console.WriteLine("  {0,15} : {1,10:0.00} $", foods[i], 365 * x[i]);

          Console.WriteLine("Nutrients (daily):");
          double[] y = M.getDailyNutrients();
          for (int i = 0; i < y.Length; ++i)
            Console.WriteLine("  {0,15} : {1,10:0.00} {2,-12} ({3})",
                              nutrients[i],
                              y[i],
                              nutrient_unit[i],
                              daily_allowance[i]);
        }
        catch (SolutionError e)
        {
          Console.WriteLine("Solution error : {0}", e.ToString());
        }
      }
    }
  }
}

duality.cs

Listing 16.9 duality.cs Click here to download.
/*
  File : duality.cs

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

  Description :  Shows how to read the dual value of a constraint
*/
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class duality
  {
    public static void Main(string[] args)
    {
      double[,] A =  new double[,] { { -0.5, 1.0 }  };
      double[]  b =  new double[] { 1.0 };
      double[]  c =  new double[] { 1.0, 1.0 };

      using (Model M = new Model("duality"))
      {
        Variable x = M.Variable("x", 2, Domain.GreaterThan(0.0));

        Constraint con = M.Constraint(Expr.Sub(b, Expr.Mul(Matrix.Dense(A), x)), Domain.EqualsTo(0.0));

        M.Objective("obj", ObjectiveSense.Minimize, Expr.Dot(c, x));

        M.Solve();
        double[] xsol = x.Level();
        double[] ysol = con.Dual();

        Console.WriteLine("x1,x2,y = %{0}, %{1}, %{2}\n", xsol[0], xsol[1], ysol[0]);
      }
    }
  }
}

elastic.cs

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

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class elastic
  {
    // Construct the model with parameters b, lambda1, lambda2
    // and with prescribed matrix A
    public static Model initializeModel(int m, int n, double[,] A) {
      Model M = new Model();
      Variable x = M.Variable("x", n);

      // t >= |Ax-b|_2 where b is a parameter
      Parameter b = M.Parameter("b", m);
      Variable  t = M.Variable();
      M.Constraint(Expr.Vstack(t, Expr.Sub(Expr.Mul(A, x), b)), Domain.InQCone());

      // p_i >= |x_i|, i=1..n
      Variable p = M.Variable(n);
      M.Constraint(Expr.Hstack(p, x), Domain.InQCone());

      // q >= |x|_2
      Variable q = M.Variable();
      M.Constraint(Expr.Vstack(q, x), Domain.InQCone());

      // Objective, parametrized with lambda1, lambda2
      // t + lambda1*sum(p) + lambda2*q
      Parameter lambda1 = M.Parameter("lambda1");
      Parameter lambda2 = M.Parameter("lambda2");
      Expression obj = Expr.Add(new Expression[] {t, Expr.Mul(lambda1, Expr.Sum(p)), Expr.Mul(lambda2, q)});
      M.Objective(ObjectiveSense.Minimize, obj);

      // Return the ready model
      return M;
    }

    public static void smallExample() {
      //Create a small example
      int m = 4;
      int n = 2;
      double[,] A = { {1.0,   2.0},
                      {3.0,   4.0},
                      {-2.0, -1.0},
                      {-4.0, -3.0} };
      double[] sol;
      Model M = initializeModel(m, n, A);

      // For convenience retrieve some elements of the model
      Parameter b = M.GetParameter("b");
      Parameter lambda1 = M.GetParameter("lambda1");
      Parameter lambda2 = M.GetParameter("lambda2");
      Variable x = M.GetVariable("x");

      // First solve
      b.SetValue(new double[]{0.1, 1.2, -1.1, 3.0});
      lambda1.SetValue(0.1);
      lambda2.SetValue(0.01);

      M.Solve();
      sol = x.Level();
      Console.WriteLine("Objective {0}, solution {1}, {2}", M.PrimalObjValue(), sol[0], sol[1]);

      // Increase lambda1
      lambda1.SetValue(0.5);
      
      M.Solve();
      sol = x.Level();
      Console.WriteLine("Objective {0}, solution {1}, {2}", M.PrimalObjValue(), sol[0], sol[1]);

      // Now change the data completely
      b.SetValue(new double[] {1.0, 1.0, 1.0, 1.0});
      lambda1.SetValue(0.0);
      lambda2.SetValue(0.0);
      
      M.Solve();
      sol = x.Level();
      Console.WriteLine("Objective {0}, solution {1}, {2}", M.PrimalObjValue(), sol[0], sol[1]);

      // And increase lamda2
      lambda2.SetValue(1.4145);
      
      M.Solve();
      sol = x.Level();
      Console.WriteLine("Objective {0}, solution {1}, {2}", M.PrimalObjValue(), sol[0], sol[1]);

      M.Dispose();
    }

    public static void Main() {
      smallExample();
    }
  }
}

facility_location.cs

Listing 16.11 facility_location.cs Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      facility_location.cs
//
//  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
////

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class facility_location
  {
    // Customer locations
    private static Matrix  customerloc =
      Matrix.Dense (new double[,]
    { { 12,  2 },
      { 15, 13 },
      { 10,  8 },
      {  0, 10 },
      {  6, 13 },
      {  5,  8 },
      { 10, 12 },
      {  4,  6 },
      {  5,  2 },
      {  1, 10 }
    } );
    private static int N = customerloc.NumRows();
    public static void Main(string[] args)
    {
      using (Model M = new Model("FacilityLocation"))
      {
        // Variable holding the facility location
        Variable f = M.Variable("facility", Set.Make(1, 2), Domain.Unbounded());
        // Variable defining the euclidian distances to each customer
        Variable d = M.Variable("dist", Set.Make(N, 1), Domain.GreaterThan(0.0));
        // Variable defining the x and y differences to each customer;
        Variable t = M.Variable("t", Set.Make(N, 2), Domain.Unbounded());
        M.Constraint("dist measure",
        Var.Hstack(new Variable[] { d, t }),
        Domain.InQCone(N, 3));

        Variable 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();
        M.WriteTask("facility_location.task");
        double[] floc = f.Level();
        Console.WriteLine("Facility location = {0},{1}", floc[0], floc[1]);
      }
    }
  }
}

gp1.cs

Listing 16.12 gp1.cs Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      gp1.cs
//
//   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
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class gp1
  {
    // Models log(sum(exp(Ax+b))) <= 0.
    // Each row of [A b] describes one of the exp-terms
    public static void logsumexp(Model      M, 
                                 double[,]  A, 
                                 Variable   x,
                                 double[]   b)
    {
      int k = A.GetLength(0);
      Variable 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 )
    public static double[] max_volume_box(double Aw, double Af, 
                                          double alpha, double beta, double gamma, double delta)
    {
      using (Model M = new Model("max_vol_box"))
      {
        Variable xyz = M.Variable(3);
        M.Objective("Objective", ObjectiveSense.Maximize, Expr.Sum(xyz));
        
        logsumexp(M, 
                  new double[,] {{1,1,0}, {1,0,1}}, 
                  xyz, 
                  new double[] {System.Math.Log(2.0/Aw), System.Math.Log(2.0/Aw)});
        
        M.Constraint(Expr.Dot(new double[] {0,1,1}, xyz), Domain.LessThan(System.Math.Log(Af)));
        M.Constraint(Expr.Dot(new double[] {1,-1,0}, xyz), Domain.InRange(System.Math.Log(alpha),System.Math.Log(beta)));
        M.Constraint(Expr.Dot(new double[] {0,-1,1}, xyz), Domain.InRange(System.Math.Log(gamma),System.Math.Log(delta)));
        
        M.SetLogHandler(Console.Out);
        M.Solve();
        
        double[] xyzVal = xyz.Level();
        double[] hwdVal = new double[3];
        for(int i=0; i<3; i++) hwdVal[i] = System.Math.Exp(xyzVal[i]);

        return hwdVal;
      }
    }

    public static void Main(String[] args)
    {
      double Aw    = 200.0;
      double Af    = 50.0;
      double alpha = 2.0;
      double beta  = 10.0;
      double gamma = 2.0;
      double delta = 10.0;
      
      double[] hwd = max_volume_box(Aw, Af, alpha, beta, gamma, delta);

      Console.WriteLine("h={0:f4} w={1:f4} d={2:f4}", hwd[0], hwd[1], hwd[2]);
    }
  }
}

lo1.cs

Listing 16.13 lo1.cs Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      lo1.cs
//
//  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
////
using System;
using mosek.fusion;
namespace mosek.fusion.example
{
  public class lo1
  {
    public static void Main(string[] args)
    {
      double[][] A =
      { new double[] { 3.0, 1.0, 2.0, 0.0 },
        new double[] { 2.0, 1.0, 3.0, 1.0 },
        new double[] { 0.0, 2.0, 0.0, 3.0 }
      };
      double[] c = { 3.0, 1.0, 5.0, 1.0  };

      // Create a model with the name 'lo1'
      Model M = new Model("lo1");

      // Create variable 'x' of length 4
      Variable 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(A[0], x), Domain.EqualsTo(30.0));
      M.Constraint("c2", Expr.Dot(A[1], x), Domain.GreaterThan(15.0));
      M.Constraint("c3", Expr.Dot(A[2], 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
      double[] sol = x.Level();
      Console.WriteLine("[x0,x1,x2,x3] = [{0}, {1}, {2}, {3} ]", sol[0], sol[1], sol[2], sol[3]);
    }
  }
}

logistic.cs

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

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class logistic
  {
    // t >= log( 1 + exp(u) ) coordinatewise
    public static void softplus(Model      M,
                                Expression t,
                                Expression u)
    {
      int n = t.GetShape()[0];
      Variable z1 = M.Variable(n);
      Variable 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
    public static Model logisticRegression(double[,]  X, 
                                           bool[]     y,
                                           double     lamb)
    {
      int n = X.GetLength(0);
      int d = X.GetLength(1);       // num samples, dimension
      
      Model M = new Model();   

      Variable theta = M.Variable("theta", d);
      Variable t     = M.Variable(n);
      Variable reg   = M.Variable();

      M.Objective(ObjectiveSense.Minimize, Expr.Add(Expr.Sum(t), Expr.Mul(lamb,reg)));
      M.Constraint(Var.Vstack(reg, theta), Domain.InQCone());
      
      double[] signs = new double[n];
      for(int i = 0; i < n; i++)
        if (y[i]) signs[i] = -1; else signs[i] = 1;

      softplus(M, t, Expr.MulElm(Expr.Mul(X, theta), signs));

      return M;
    }

    public static void Main(String[] args)
    {
      // Test: detect and approximate a circle using degree 2 polynomials
      int n = 30;
      double[,] X = new double[n*n, 6];
      bool[] Y     = new bool[n*n];

      for(int i=0; i<n; i++) 
      for(int j=0; j<n; j++)
      {
        int k = i*n+j;
        double x = -1 + 2.0*i/(n-1);
        double y = -1 + 2.0*j/(n-1);
        X[k,0] = 1.0; X[k,1] = x; X[k,2] = y; X[k,3] = x*y;
        X[k,4] = x*x; X[k,5] = y*y;
        Y[k] = (x*x+y*y>=0.69);
      }

      Model M = logisticRegression(X, Y, 0.1);
      Variable theta = M.GetVariable("theta");

      M.Solve();
      for(int i=0; i<6; i++)
        Console.WriteLine(theta.Level()[i]);
    }
  }
}

lownerjohn_ellipsoid.cs

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

   File:      lownerjohn_ellipsoid.cs

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


   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.


   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

     minimize         t
     subject to       t       >= det(P)^(-1/n)
                 || P*xi + c ||_2 <= 1,  i=1,...,m
                   P is PSD.

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

using System;

namespace mosek.fusion.example
{
  public class lownerjohn_ellipsoid
  {
    /**
        Purpose: Models the convex set

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

        using three-dimensional power cones
    */
    public static void geometric_mean(Model M, Variable x, Variable t)
    {
      int n = (int)x.GetSize();
      if (n==1)
        M.Constraint(Expr.Sub(t, x), Domain.LessThan(0.0));
      else
      {
        Variable t2 = M.Variable();
        M.Constraint(Var.Hstack(t2, x.Index(n-1), t), Domain.InPPowerCone(1-1.0/n));
        geometric_mean(M, x.Slice(0,n-1), t2);
      }
    }

    /**
        Purpose: Models the hypograph of the n-th power of the
        determinant of a positive definite matrix.

        The convex set (a hypograph)

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

        can be modeled as the intersection of a semidefinite cone

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

        and a number of rotated quadratic cones and affine hyperplanes,

          t <= (Z11*Z22*...*Znn)^{1/n}  (see geometric_mean).

        References:
        [1] "Lectures on Modern Optimization", Ben-Tal and Nemirovski, 2000.
     */
    public static Variable det_rootn(Model M, Variable t, int n)
    {
      // Setup variables
      Variable Y = M.Variable(Domain.InPSDCone(2 * n));

      Variable X   = Y.Slice(new int[]{0, 0}, new int[]{n, n});
      Variable Z   = Y.Slice(new int[]{0, n}, new int[]{n, 2 * n});
      Variable DZ  = Y.Slice(new int[]{n, n}, new int[]{2 * n, 2 * n});

      // Z is lower-triangular
      int[,] low_tri = new int[n*(n-1)/2,2];
      int k = 0;
      for(int i = 0; i < n; i++)
        for(int j = i+1; j < n; j++)
          { low_tri[k,0] = i; low_tri[k,1] = j; ++k; }
      M.Constraint(Z.Pick(low_tri), Domain.EqualsTo(0.0));
      // DZ = Diag(Z)
      M.Constraint(Expr.Sub(DZ, Expr.MulElm(Z, Matrix.Eye(n))), Domain.EqualsTo(0.0));

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

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

    public static Tuple<double[], double[]> lownerjohn_inner(double[][] A, double[] b)
    {
      using( Model M = new Model("lownerjohn_inner"))
      {
        int m = A.Length;
        int n = A[0].Length;

        // Setup variables
        Variable t = M.Variable("t", 1, Domain.GreaterThan(0.0));
        Variable C = det_rootn(M, t, n);
        Variable d = M.Variable("d", n, Domain.Unbounded());

        // (bi - ai^T*d, C*ai) \in Q
        for (int i = 0; i < m; ++i)
          M.Constraint("qc" + i, Expr.Vstack(Expr.Sub(b[i], Expr.Dot(A[i], d)), Expr.Mul(C, A[i])),
                       Domain.InQCone().Axis(0) );

        // Objective: Maximize t
        M.Objective(ObjectiveSense.Maximize, t);
        M.Solve();

        return Tuple.Create(C.Level(), d.Level());
      }
    }

    public static Tuple<double[], double[]> lownerjohn_outer(double[,] x)
    {
      using( Model M = new Model("lownerjohn_outer") )
      {
        int m = x.GetLength(0);
        int n = x.GetLength(1);

        // Setup variables
        Variable t = M.Variable("t", 1, Domain.GreaterThan(0.0));
        Variable P = det_rootn(M, t, n);
        Variable c = M.Variable("c", Domain.Unbounded().WithShape(1,n));

        // (1, P(*xi+c)) \in Q
        M.Constraint("qc",
                     Expr.Hstack(Expr.Ones(m), Expr.Sub(Expr.Mul(x,P), Expr.Repeat(c,m,0))),
                     Domain.InQCone().Axis(1));

        // Objective: Maximize t
        M.Objective(ObjectiveSense.Maximize, t);
        M.Solve();

        return Tuple.Create(P.Level(), c.Level());
      }
    }

    public static void Main(string[] argv)
    {
      double[,] p = new double[,] { {0.0, 0.0},
                                    {1.0, 3.0},
                                    {5.5, 4.5},
                                    {7.0, 4.0},
                                    {7.0, 1.0},
                                    {3.0, -2.0} };
      int n = p.GetLength(0);
      double[][] A = new double[n][];
      double[]   b = new double[n];
      for (int i = 0; i < n; ++i)
      {
        A[i] = new double[] { -p[i,1] + p[(i - 1 + n) % n,1], p[i,0] - p[(i - 1 + n) % n,0] };
        b[i] = A[i][0] * p[i,0] + A[i][1] * p[i,1];
      }

      Tuple<double[], double[]> ResInn = lownerjohn_inner(A, b);
      Tuple<double[], double[]> ResOut = lownerjohn_outer(p);
      double[] Ci = ResInn.Item1;
      double[] di = ResInn.Item2;
      double[] Po = ResOut.Item1;
      double[] co = ResOut.Item2;
      Console.WriteLine("Inner:");
      Console.WriteLine("  C = [{0}]", string.Join(", ", Ci));
      Console.WriteLine("  d = [{0}]", string.Join(", ", di));
      Console.WriteLine("Outer:");
      Console.WriteLine("  P = [{0}]", string.Join(", ", Po));
      Console.WriteLine("  c = [{0}]", string.Join(", ", co));
    }
  }
}

lpt.cs

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


using System;
using System.Collections.Generic;

using mosek.fusion;

namespace mosek.fusion.example
{
  public class lpt
  {
    public static void Main(string [] args)
    {
      int n = 30;           //Number of tasks
      int m = 6;              //Number of processors

      double lb = 1.0;         //The range of lengths of short tasks
      double ub = 5.0;

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

      double[] T = new double[n];
      Random gen = new Random(0);
      for (int i = 0; i < n_short; i++)
        T[i] = gen.NextDouble() * (ub - lb) + lb;
      for (int i = n_short; i < n; i++)
        T[i] = 20 * (gen.NextDouble() * (ub - lb) + lb);
      Array.Sort<double>(T);

      Model M = new Model("Multi-processor scheduling");

      Variable x = M.Variable("x", new int[] {m, n}, Domain.Binary());
      Variable t = M.Variable("t", 1, Domain.Unbounded());

      M.Constraint( Expr.Sum(x, 0), Domain.EqualsTo(1.0) );
      M.Constraint( Expr.Sub( Var.Repeat(t, m), Expr.Mul(x, T)), Domain.GreaterThan(0.0) );

      M.Objective( ObjectiveSense.Minimize, t );

      //Computing LPT solution
      double [] init    = new double[n * m];
      double [] schedule = new double[m];

      for (int i = n - 1; i >= 0; i--)
      {
        int next = 0;
        for (int j = 1; j < m; j++)
          if (schedule[j] < schedule[next]) next = j;
        schedule[next] += T[i];
        init[next * n + i] = 1;
      }

      //Comment this line to switch off feeding in the initial LPT solution
      x.SetLevel(init);

      M.SetLogHandler(Console.Out);

      M.SetSolverParam("mioTolRelGap", 0.01);
      M.Solve();

      try
      {
        Console.Write("initial solution: \n");
        for (int i = 0; i < m; i++)
        {
          Console.Write("M {0} [", i);
          for (int y = 0; y < n; y++)
            Console.Write( "{0}, ", (int) init[i * n + y]);
          Console.Write("]\n");
        }
        Console.Write("MOSEK solution:\n");
        for (int i = 0; i < m; i++)
        {
          Console.Write("M {0} [", i);
          for (int y = 0; y < n; y++)
            Console.Write( "{0}, ", (int)x.Index(i, y).Level()[0]);
          Console.Write("]\n");
        }
      }
      catch (SolutionError e) {}
    }
  }
}

mico1.cs

Listing 16.17 mico1.cs Click here to download.
//
//    Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//    File:    mico1.cs
//
//    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
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class mico1
  {
    public static void Main(string[] args)
    {
      using (Model M = new Model("mico1"))
      {
        Variable x = M.Variable(Domain.Integral(Domain.Unbounded()));
        Variable y = M.Variable(Domain.Integral(Domain.Unbounded()));
        Variable 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();

        Console.WriteLine("x, y = {0}, {1}", x.Level()[0], y.Level()[0]);
      }
    }
  }
}

milo1.cs

Listing 16.18 milo1.cs Click here to download.
//
//    Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//    File:    milo1.cs
//
//    Purpose:  Demonstrates how to solve a small mixed
//              integer linear optimization problem.
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class milo1
  {
    public static void Main(string[] args)
    {
      double[][] A =
      { new double[] { 50.0, 31.0 },
        new double[] { 3.0,  -2.0 }
      };
      double[] c = { 1.0, 0.64 };
      using (Model M = new Model("milo1"))
      {
        Variable 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(A[0], x), Domain.LessThan(250.0));
        M.Constraint("c2", Expr.Dot(A[1], 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
        double[] sol = x.Level();
        Console.WriteLine("x1,x2 = {0}, {1}", sol[0], sol[1]);
        double miorelgap = M.GetSolverDoubleInfo("mioObjRelGap");
        double mioabsgap = M.GetSolverDoubleInfo("mioObjAbsGap");
        Console.WriteLine("MIP rel gap = {0} ({0})", miorelgap, mioabsgap);
      }
    }
  }
}

mioinitsol.cs

Listing 16.19 mioinitsol.cs Click here to download.
//
//    Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//    File:    mioinitsol.cs
//
//    Purpose:  Demonstrates how to solve a small mixed
//             integer linear optimization problem
//             providing an initial feasible solution.
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class mioinitsol
  {
    public static void Main(string[] args)
    {
      double[] c = { 7.0, 10.0, 1.0, 5.0 };
      int n = 4;
      using (Model M = new Model("mioinitsol"))
      {

        Variable x = M.Variable("x", n, Domain.GreaterThan(0.0));
        x.Slice(0,3).MakeInteger();

        // Create the constraint
        M.Constraint(Expr.Sum(x), Domain.LessThan(2.5));

        // Set the objective function to (c^T * x)
        M.Objective("obj", ObjectiveSense.Maximize, Expr.Dot(c, x));

        // Assign values to integer variables.
        // We only set a slice of x     
        double[] init_sol = { 1, 1, 0 };
        x.Slice(0,3).SetLevel( init_sol );

        // Solve the problem
        M.Solve();

        // Get the solution values
        double[] sol = x.Level();
        Console.Write("x = [");
        for(int i=0;i<n;i++)
        {
          Console.Write("{0}, ",sol[i]);
        }
        Console.WriteLine("]");

        // Was the initial solution used?
        int constr = M.GetSolverIntInfo("mioConstructSolution");
        double constrVal = M.GetSolverDoubleInfo("mioConstructSolutionObj");
        Console.WriteLine("Initial solution utilization: " + constr);
        Console.WriteLine("Initial solution objective: " +  constrVal);
      }
    }
  }
}

modelLib.cs

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

  File :      modelLib.cs

  Purpose: Library of simple building blocks in Mosek Fusion.
*/
using System;

namespace mosek.fusion.example 
{
  public class modelLib 
  {
    // Duplicate variables
    // x = y
    public static void dup(Model M, Variable x, Variable y) {
      M.Constraint(Expr.Sub(x,y), Domain.EqualsTo(0.0));
    }

    // Absolute value
    // t >= |x|, where t, x have the same shape
    public static void abs(Model M, Variable t, Variable x) {
      M.Constraint(Expr.Add(t,x), Domain.GreaterThan(0.0));
      M.Constraint(Expr.Sub(t,x), Domain.GreaterThan(0.0));
    }

    // 1-norm
    // t >= sum( |x_i| ), x is a vector Variable
    public static void norm1(Model M, Variable t, Variable x) {
      Variable u = M.Variable(x.GetShape(), Domain.Unbounded());
      abs(M, u, x);
      M.Constraint(Expr.Sub(t, Expr.Sum(u)), Domain.GreaterThan(0.0));
    }

    // Square
    // t >= x^2
    public static void sq(Model M, Variable t, Variable x) {
      M.Constraint(Expr.Hstack(0.5, t, x), Domain.InRotatedQCone());
    }

    // 2-norm
    // t >= sqrt(x_1^2 + ... + x_n^2) where x is a vector
    public static void norm2(Model M, Variable t, Variable x) {
      M.Constraint(Expr.Vstack(t, x), Domain.InQCone());
    }

    // Power with exponent > 1
    // t >= |x|^p (where p>1)
    public static void pow(Model M, Variable t, Variable x, double p) {
      M.Constraint(Expr.Hstack(t, 1, x), Domain.InPPowerCone(1.0/p));
    }

    // Inverse of power 
    // t >= 1/|x|^p, x>0 (where p>0)
    public static void pow_inv(Model M, Variable t, Variable x, double p) {
      M.Constraint(Expr.Hstack(t, x, 1), Domain.InPPowerCone(1.0/(1.0+p)));
    }

    // p-norm, p>1
    // t >= \|x\|_p (where p>1), x is a vector Variable
    public static void pnorm(Model M, Variable t, Variable x, double p) {
      int n = (int) x.GetSize();
      Variable r = M.Variable(n);
      M.Constraint(Expr.Sub(t, Expr.Sum(r)), Domain.EqualsTo(0.0));
      M.Constraint(Expr.Hstack(Var.Repeat(t,n), r, x), Domain.InPPowerCone(1.0-1.0/p));
    }

    // Geometric mean
    // |t| <= (x_1...x_n)^(1/n), x_i>=0, x is a vector Variable of length >= 1
    public static void geo_mean(Model M, Variable t, Variable x) {
      int n = (int) x.GetSize();
      if (n==1) {
          abs(M, x, t);
      }
      else {
          Variable t2 = M.Variable();
          M.Constraint(Expr.Hstack(t2, x.Index(n-1), t), Domain.InPPowerCone(1.0-1.0/n));
          geo_mean(M, t2, x.Slice(0,n-1));
      }
    }

    // Logarithm
    // t <= log(x), x>=0
    public static void log(Model M, Variable t, Variable x) {
      M.Constraint(Expr.Hstack(x, 1, t), Domain.InPExpCone());
    }

    // Exponential
    // t >= exp(x)
    public static void exp(Model M, Variable t, Variable x) {
      M.Constraint(Expr.Hstack(t, 1, x), Domain.InPExpCone());
    }

    // Entropy
    // t >= x * log(x), x>=0
    public static void ent(Model M, Variable t, Variable x) {
      M.Constraint(Expr.Hstack(1, x, Expr.Neg(t)), Domain.InPExpCone());
    }

    // Relative entropy
    // t >= x * log(x/y), x,y>=0
    public static void relent(Model M, Variable t, Variable x, Variable y) {
      M.Constraint(Expr.Hstack(y, x, Expr.Neg(t)), Domain.InPExpCone());
    }

    // Log-sum-exp
    // log( sum_i(exp(x_i)) ) <= t, where x is a vector
    public static void logsumexp(Model M, Variable t, Variable x) {
      int n = (int) x.GetSize();
      Variable u = M.Variable(n);
      M.Constraint(Expr.Hstack(u, Expr.ConstTerm(n, 1.0), Expr.Sub(x, Var.Repeat(t, n))), Domain.InPExpCone());
      M.Constraint(Expr.Sum(u), Domain.LessThan(1.0));
    }

    // Semicontinuous variable
    // x = 0 or a <= x <= b
    public static void semicontinuous(Model M, Variable x, double a, double b) {
      Variable u = M.Variable(x.GetShape(), Domain.Binary());
      M.Constraint(Expr.Sub(x, Expr.Mul(a, u)), Domain.GreaterThan(0.0));
      M.Constraint(Expr.Sub(x, Expr.Mul(b, u)), Domain.LessThan(0.0));
    }

    // Indicator variable
    // x!=0 implies t=1. Assumes that |x|<=1 in advance.
    public static void indicator(Model M, Variable t, Variable x) {
      M.Constraint(t, Domain.InRange(0,1));
      t.MakeInteger();
      abs(M, t, x);
    }

    // Logical OR
    // x OR y, where x, y are binary
    public static void logic_or(Model M, Variable x, Variable y) {
      M.Constraint(Expr.Add(x, y), Domain.GreaterThan(1.0));
    }
    // x_1 OR ... OR x_n, where x is a binary vector
    public static void logic_or_vect(Model M, Variable x) {
      M.Constraint(Expr.Sum(x), Domain.GreaterThan(1.0));
    }

    // SOS1 (NAND)
    // at most one of x_1,...,x_n, where x is a binary vector (SOS1 Constraint)
    public static void logic_sos1(Model M, Variable x) {
      M.Constraint(Expr.Sum(x), Domain.LessThan(1.0));
    }
    // NOT(x AND y), where x, y are binary
    public static void logic_nand(Model M, Variable x, Variable y) {
      M.Constraint(Expr.Add(x, y), Domain.LessThan(1.0));
    }

    // Cardinality bound
    // At most k of entries in x are nonzero, assuming in advance |x_i|<=1.
    public static void card(Model M, Variable x, int k) {
      Variable t = M.Variable(x.GetShape(), Domain.Binary());
      abs(M, t, x);
      M.Constraint(Expr.Sum(t), Domain.LessThan(k));
    }

    // This is just a syntactic test without much sense
    public static void testModels() {
      Model M = new Model();
      Variable x = M.Variable();
      Variable y = M.Variable();
      Variable t = M.Variable();
      Variable p = M.Variable(5);
      Variable a = M.Variable(new int[]{10,2});
      Variable b = M.Variable(new int[]{10,2});
      Variable e = M.Variable(Domain.Binary());
      Variable f = M.Variable(Domain.Binary());

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

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

      M.Constraint(Expr.Sub(x, y), Domain.GreaterThan(3.0));
      norm2(M, t.Index(0), Var.Vstack(x,y));
      log  (M, t.Index(1), y);
      pow  (M, t.Index(2), x, 1.5);

      M.Objective(ObjectiveSense.Maximize, Expr.Dot(t, new double[]{-1,1,-1}));

      M.SetLogHandler(Console.Out);
      M.Solve();
    }
      
    // Full test
    public static void Main(String[] argv) {
      testModels();
      testExample();
    }
  }
}

nearestcorr.cs

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

   File:      nearestcorr.cs

   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
*/
using mosek.fusion;
using System;

namespace mosek.fusion.example
{
  public class nearestcorr
  {
    /** Assuming that e is an NxN expression, return the lower triangular part as a vector.
    */
    public static Expression Vec(Expression e)
    {
      int N       = e.GetShape()[0];
      int[] msubi = new int[N * (N + 1) / 2],
      msubj = new int[N * (N + 1) / 2];
      double[] mcof = new double[N * (N + 1) / 2];

      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;
          if (i == j) mcof[k] = 1.0;
          else        mcof[k] = Math.Sqrt(2);
        }

      var S = Matrix.Sparse(N * (N + 1) / 2, N * N, msubi, msubj, mcof);
      return Expr.Mul(S, Expr.Flatten(e));
    }

    public static void Main(string[] argv)
    {
      int N = 5;
      var A = Matrix.Dense( new double[,]
      { { 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}
      });

      var gammas = Range(0.0, 2.0, 0.1);
      var res  = new double[gammas.Length];
      var rank = new int[gammas.Length];

      nearestcorr_frobenius(A);
      nearestcorr_nn(A, gammas, res, rank);

      for (int i = 0; i < gammas.Length; ++i)
        Console.WriteLine("gamma = {0:f1}, res={1:e3}, rank = {2}", gammas[i], res[i], rank[i]);
    }

    public static void nearestcorr_nn(Matrix A, double[] gammas, double[] res, int[] rank)
    {
      int N = A.NumRows();
      using (var M = new Model("NucNorm"))
      {
        // Setup variables
        var t = M.Variable("t", 1, Domain.Unbounded());
        var X = M.Variable("X", Domain.InPSDCone(N));
        var w = M.Variable("w", N, Domain.GreaterThan(0.0));

        // (t, vec (X + diag(w) - A)) in Q
        var D = Expr.MulElm( Matrix.Eye(N), Var.Repeat(w, 1, N) );
        M.Constraint( Expr.Vstack( t, Vec(Expr.Sub(Expr.Add(X, D), A)) ), Domain.InQCone() );

        for (var k = 0; k < gammas.Length; ++k)
        {
          // Objective: Minimize t + gamma*Tr(X)
          var gamm_trX = Expr.Mul( gammas[k], Expr.Sum(X.Diag()) );
          M.Objective(ObjectiveSense.Minimize, Expr.Add(t, gamm_trX));
          M.Solve();

          // Find the eigenvalues of X and approximate rank
          var d = new double[N];
          mosek.LinAlg.syeig(mosek.uplo.lo, N, X.Level(), d);
          var rnk = 0; foreach (var v in d) if (v > 1e-6) ++rnk;

          res[k] = t.Level()[0];
          rank[k] = rnk;
        }
      }
    }

    public static void nearestcorr_frobenius(Matrix A)
    {
      int N = A.NumRows();
      using (var M = new Model("NearestCorrelation"))
      {
        // Setting up the variables
        var X = M.Variable("X", Domain.InPSDCone(N));
        var 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
        Console.WriteLine("X = \n{0}", mattostr(X.Level(), N));
        Console.WriteLine("t = {0}", mattostr(t.Level(), N));
      }
    }

    /* Utility methods */
    private static string mattostr(double[] a, int n)
    {
      var c = new System.Globalization.CultureInfo("en-US");
      var b = new System.Text.StringBuilder();
      if (a.Length == 1 )
      {
        b.Append(String.Format(c, "{0:e3}", a[0]));
      }
      else
      {
        for (int i = 0; i < a.Length; ++i)
        {
          if (  i  % n == 0 )
            b.Append("[");
          b.Append(" ").Append(String.Format(c, "{0:e3}", a[i]));
          if ( ( i + 1 ) % n == 0 )
            b.Append("]\n");
        }
      }
      return b.ToString();
    }

    public static double[] Range (double start, double stop, double step)
    {
      int len;
      if (start < stop && step > 0.0)
        len = 1 + (int)((stop - start - 1) / step);
      else if (stop < start && step < 0)
        len = 1 + (int)((start - stop - 1) / (- step));
      else
        len = 0;
      double[] res = new double[len];
      double v = start;
      for (int i = 0; i < len; ++i, v += step)
        res[i] = v;
      return res;
    }
  }
}

opt_server_sync.cs

Listing 16.22 opt_server_sync.cs Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      opt_server_sync.cs
//
//  Purpose :   Demonstrates how to use MOSEK OptServer
//              to solve optimization problem synchronously
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class opt_server_sync
  {
    public static void Main(string[] args)
    {

      if (args.Length<1) {
        Console.WriteLine("Missing argument, syntax is:");
        Console.WriteLine("   opt_server_sync serveraddr");
        return;
      }

      String addr = args[0];

      // Setup a simple test problem
      Model M = new Model("testOptServer");
      Variable x = M.Variable("x", 3, Domain.GreaterThan(0.0));
      M.Constraint("lc", Expr.Dot(new double[] {1.0, 1.0, 2.0}, x), Domain.EqualsTo(1.0));
      M.Objective("obj", ObjectiveSense.Minimize, Expr.Sum(x));

      // Attach log handler
      M.SetLogHandler(Console.Out);

      // Set OptServer URL
      M.OptserverHost(addr);

      // Solve the problem on the OptServer
      M.Solve();

      // Get the solution
      double[] solx = x.Level();
      Console.WriteLine("x1,x2,x3 = {0}, {1}, {2}", solx[0], solx[1], solx[2]);
    }
  }
}

parameters.cs

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

  File :      parameters.cs

  Purpose :   Demonstrates a very simple example about how to set
              parameters and read information items
              with MOSEK Fusion
*/
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class parameters
  {
    public static void Main(string[] args)
    {
      // Create the Model        
      using (Model M = new Model()) {
        Console.WriteLine("Test MOSEK parameter get/set functions");

        // 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 e) {
          Console.WriteLine("Wrong parameter value");
        }


        // Define and solve an optimization problem here
        // M.Solve()
        // After optimization: 

        Console.WriteLine("Get MOSEK information items");

        double tm = M.GetSolverDoubleInfo("optimizerTime");
        int it = M.GetSolverIntInfo("intpntIter");

        Console.WriteLine("Time: " + tm);
        Console.WriteLine("Iterations: " + it);
      }
    }
  }
}

portfolio_1_basic.cs

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

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

  Description :  Implements a basic portfolio optimization model.
  
  Note:
    This example uses LINQ, which is only available in .NET Framework 3.5 and later.
*/

using System.IO;
using System;
using System.Linq;
using System.Globalization;

namespace mosek.fusion.example
{
  public class portfolio_1_basic
  {

    private static int        n      = 3;
    private static double     w      = 1.0;
    private static double[]   mu     = {0.1073, 0.0737, 0.0627};
    private static double[]   x0     = {0.0, 0.0, 0.0};
    private static double[]   gammas = {0.035, 0.040, 0.050, 0.060, 0.070, 0.080, 0.090};
    private static double[,] GT     = {
      { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
      { 0.0              , 0.102863378954911  , -0.00222873156550421},
      { 0.0              , 0.0                ,  0.0338148677744977 }
    };
    public static double sum(double[] x)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i];
      return r;
    }

    public static double dot(double[] x, double[] y)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++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
    */
    public static double BasicMarkowitz
    ( int n,
      double[] mu,
      double[,]GT,
      double[] x0,
      double   w,
      double   gamma)
    {

      using( Model M = new Model("Basic Markowitz"))
      {

        // Redirect log output from the solver to stdout for debugging.
        // if uncommented.
        //M.SetLogHandler(Console.Out);

        // Defines the variables (holdings). Shortselling is not allowed.
        Variable x = M.Variable("x", n, Domain.GreaterThan(0.0));

        //  Maximize expected return
        M.Objective("obj", ObjectiveSense.Maximize, Expr.Dot(mu, x));

        // The amount invested  must be identical to intial wealth
        M.Constraint("budget", Expr.Sum(x), Domain.EqualsTo(w + sum(x0)));

        // Imposes a bound on the risk
        M.Constraint("risk", Expr.Vstack(gamma, Expr.Mul(GT, x)), Domain.InQCone());

        // Solves the model.
        M.Solve();

        return dot(mu, x.Level());
      }
    }


    /*
      The example. Reads in data and solves the portfolio models.
     */
    public static void Main(string[] argv)
    {

      Console.WriteLine("\n-------------------------------------------------------------------");
      Console.WriteLine("Basic Markowitz portfolio optimization");
      Console.WriteLine("---------------------------------------------------------------------");
      foreach (var gamma in gammas)
      {
        double res = BasicMarkowitz(n, mu, GT, x0, w, gamma);
        Console.WriteLine("Expected return: {0,-12:f4}  St deviation: {1,-12:f4} ", res, gamma);
      }
    }
  }
}

portfolio_2_frontier.cs

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

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

  Description :  Implements a basic portfolio optimization model.
                 Computes points on the efficient frontier.

  Note:
    This example uses LINQ, which is only available in .NET Framework 3.5 and later.
*/

using System.IO;
using System;
using System.Linq;
using System.Globalization;

namespace mosek.fusion.example
{
  public class portfolio_2_frontier
  {

    private static int        n      = 3;
    private static double     w      = 1.0;
    private static double[]   mu     = {0.1073, 0.0737, 0.0627};
    private static double[]   x0     = {0.0, 0.0, 0.0};
    private static double[,]  GT     = {
      { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
      { 0.0              , 0.102863378954911  , -0.00222873156550421},
      { 0.0              , 0.0                ,  0.0338148677744977 }
    };
    public static double sum(double[] x)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i];
      return r;
    }

    public static double dot(double[] x, double[] y)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++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)
     */
    public static void EfficientFrontier
    ( int n,
      double[]    mu,
      double [,]  GT,
      double[]    x0,
      double      w,
      double[]    alphas,
      double[]    frontier_mux,
      double[]    frontier_s)
    {
      using(Model M = new Model("Efficient frontier"))
      {
        //M.SetLogHandler(Console.Out);

        // Defines the variables (holdings). Shortselling is not allowed.
        Variable x = M.Variable("x", n, Domain.GreaterThan(0.0)); // Portfolio variables
        Variable 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 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 (int i = 0; i < alphas.Length; ++i) {
          alpha.SetValue(alphas[i]);
          M.Solve();

          frontier_mux[i] = dot(mu, x.Level());
          frontier_s[i]   = s.Level()[0];
        }
      }
    }


    /*
      The example. Reads in data and solves the portfolio models.
     */
    public static void Main(string[] argv)
    {

      // Some predefined alphas are chosen
      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 };
      int      niter = alphas.Length;
      double[] frontier_mux = new double[niter];
      double[] frontier_s   = new double[niter];

      EfficientFrontier(n, mu, GT, x0, w, alphas, frontier_mux, frontier_s);
      Console.WriteLine("\n-------------------------------------------------------------------------");
      Console.WriteLine("Efficient frontier") ;
      Console.WriteLine("------------------------------------------------------------------------");
      Console.WriteLine("{0,-12}  {1,-12}  {2,-12}", "alpha", "return", "variance") ;
      for (int i = 0; i < frontier_mux.Length; ++i)
        Console.WriteLine("{0,-12:f4}  {1,-12:e4}  {2,-12:e4}", alphas[i], frontier_mux[i], frontier_s[i]);

    }
  }
}

portfolio_3_impact.cs

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

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

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

  Note:
    This example uses LINQ, which is only available in .NET Framework 3.5 and later.
*/

using System.IO;
using System;
using System.Linq;
using System.Globalization;

namespace mosek.fusion.example
{
  public class portfolio_3_impact
  {

    private static int        n      = 3;
    private static double     w      = 1.0;
    private static double[]   mu     = {0.1073, 0.0737, 0.0627};
    private static double[]   x0     = {0.0, 0.0, 0.0};
    private static double[]   gammas = {0.035, 0.040, 0.050, 0.060, 0.070, 0.080, 0.090};
    private static double[,] GT     = {
      { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
      { 0.0              , 0.102863378954911  , -0.00222873156550421},
      { 0.0              , 0.0                ,  0.0338148677744977 }
    };
    public static double sum(double[] x)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i];
      return r;
    }

    public static double dot(double[] x, double[] y)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++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

    */
    public static void MarkowitzWithMarketImpact
    ( int n,
      double[] mu,
      double[,]GT,
      double[] x0,
      double   w,
      double   gamma,
      double[] m,
      double[] xsol,
      double[] tsol)
    {
      using(Model M = new Model("Markowitz portfolio with market impact"))
      {

        //M.SetLogHandler(Console.Out);

        // Defines the variables. No shortselling is allowed.
        Variable x = M.Variable("x", n, Domain.GreaterThan(0.0));

        // Variables computing market impact
        Variable 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();

        if (xsol != null)
          Array.Copy(x.Level(), xsol, n);
        if (tsol != null)
          Array.Copy(t.Level(), tsol, n);
      }
    }


    /*
      The example. Reads in data and solves the portfolio models.
     */
    public static void Main(string[] argv)
    {
      // Somewhat arbirtrary choice of m
      double[] m = new double[n]; for (int i = 0; i < n; ++i) m[i] = 1.0e-2;
      double[] x = new double[n];
      double[] t = new double[n];
      double gamma = gammas[0];

      MarkowitzWithMarketImpact(n, mu, GT, x0, w, gamma, m, x, t);
      Console.WriteLine("\n-----------------------------------------------------------------------");
      Console.WriteLine("Markowitz portfolio optimization with market impact cost");
      Console.WriteLine("------------------------------------------------------------------------");
      Console.WriteLine("Expected return: {0:e4} St deviation: {1:e4} Market impact cost: {2:e4}\n",
                        dot(mu, x),
                        gamma,
                        dot(m, t));

    }
  }
}

portfolio_4_transcost.cs

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

  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.

  Note:
    This example uses LINQ, which is only available in .NET Framework 3.5 and later.
*/

using System.IO;
using System;
using System.Linq;
using System.Globalization;

namespace mosek.fusion.example
{
  public class portfolio_4_transcost
  {

    private static int        n      = 3;
    private static double     w      = 1.0;
    private static double[]   mu     = {0.1073, 0.0737, 0.0627};
    private static double[]   x0     = {0.0, 0.0, 0.0};
    private static double[]   gammas = {0.035, 0.040, 0.050, 0.060, 0.070, 0.080, 0.090};
    private static double[,] GT     = {
      { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
      { 0.0              , 0.102863378954911  , -0.00222873156550421},
      { 0.0              , 0.0                ,  0.0338148677744977 }
    };
    public static double sum(double[] x)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i];
      return r;
    }

    public static double dot(double[] x, double[] y)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++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

    */
    public static double[] MarkowitzWithTransactionsCost
    ( int n,
      double[] mu,
      double[,]GT,
      double[] x0,
      double   w,
      double   gamma,
      double[] f,
      double[] g)
    {

      // Upper bound on the traded amount
      double[] u = new double[n];
      {
        double v = w + sum(x0);
        for (int i = 0; i < n; ++i) u[i] = v;
      }

      using( Model M = new Model("Markowitz portfolio with transaction costs") )
      {
        //M.SetLogHandler(Console.Out);

        // Defines the variables. No shortselling is allowed.
        Variable x = M.Variable("x", n, Domain.GreaterThan(0.0));

        // Addtional "helper" variables
        Variable z = M.Variable("z", n, Domain.Unbounded());
        // Binary varables
        Variable 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));

        //M.constraint("trade", Expr.hstack(z,Expr.sub(x,x0)), Domain.inQcone())"

        // Consraints for turning y off and on. z-diag(u)*y<=0 i.e. z_j <= u_j*y_j
        M.Constraint("y_on_off", Expr.Sub(z, Expr.Mul(Matrix.Diag(u), y)), Domain.LessThan(0.0));

        // Integer optimization problems can be very hard to solve so limiting the
        // maximum amount of time is a valuable safe guard
        M.SetSolverParam("mioMaxTime", 180.0);
        M.Solve();
        return x.Level();
      }
    }


    /*
      The example. Reads in data and solves the portfolio models.
     */
    public static void Main(string[] argv)
    {

      double[] f = new double[n]; for (var i = 0; i < n; ++i) f[i] = 0.01;
      double[] g = new double[n]; for (var i = 0; i < n; ++i) g[i] = 0.001;
      double gamma = gammas[0];

      double[] x = MarkowitzWithTransactionsCost(n, mu, GT, x0, w, gamma, f, g);
      Console.WriteLine("\n-------------------------------------------------------------------------");
      Console.WriteLine("Markowitz portfolio optimization with transaction cost");
      Console.WriteLine("------------------------------------------------------------------------");
      Console.WriteLine("Expected return: {0:e4}", dot(mu, x));

    }
  }
}

portfolio_5_card.cs

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

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

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

  Note:
    This example uses LINQ, which is only available in .NET Framework 3.5 and later.
*/

using System.IO;
using System;
using System.Linq;
using System.Globalization;

namespace mosek.fusion.example
{
  public class portfolio_5_card
  {

    private static int        n      = 3;
    private static double     w      = 1.0;
    private static double[]   mu     = {0.1073, 0.0737, 0.0627};
    private static double[]   x0     = {0.0, 0.0, 0.0};
    private static double     gamma  = 0.035;
    private static double[,]  GT     = {
      { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
      { 0.0              , 0.102863378954911  , -0.00222873156550421},
      { 0.0              , 0.0                ,  0.0338148677744977 }
    };
    private static int[]      kValues = { 1, 2, 3 };

    public static double sum(double[] x)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i];
      return r;
    }

    public static double dot(double[] x, double[] y)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++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

    */
    public static double[,] MarkowitzWithCardinality
    ( int n,
      double[] mu,
      double[,]GT,
      double[] x0,
      double   w,
      double   gamma,
      int[]    kValues)
    {
      // Upper bound on the traded amount
      double[] u = new double[n];
      {
        double v = w + sum(x0);
        for (int i = 0; i < n; ++i) u[i] = v;
      }

      using( Model M = new Model("Markowitz portfolio with cardinality bounds") )
      {
        //M.SetLogHandler(Console.Out);

        // Defines the variables. No shortselling is allowed.
        Variable x = M.Variable("x", n, Domain.GreaterThan(0.0));

        // Addtional "helper" variables
        Variable z = M.Variable("z", n, Domain.Unbounded());
        // Binary varables
        Variable 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
        Parameter 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 parameter k 
        double[,] results = new double[kValues.Length,n];

        for(int i = 0; i < kValues.Length; i++) {
          cardMax.SetValue(kValues[i]);
          M.Solve();
          double[] sol = x.Level();
          for(int j = 0; j < n; j++) results[i,j] = sol[j];
        }   

        return results;
      }
    }


    /*
      The example. Reads in data and solves the portfolio models.
     */
    public static void Main(string[] argv)
    {
      Console.WriteLine("\n-------------------------------------------------------------------------");
      Console.WriteLine("Markowitz portfolio optimization with cardinality constraints");
      Console.WriteLine("------------------------------------------------------------------------");

      double[,] results = MarkowitzWithCardinality(n, mu, GT, x0, w, gamma, kValues);

      for(int k=1; k<=n; k++) 
      {
        Console.Write("Bound:  {0:d}  Solution: ", k);
        for(int i=0; i<n; i++)
        {
          Console.Write("{0:f4} ", results[k-1,i]);
        }
        Console.WriteLine();
      }
    }
  }
}

pow1.cs

Listing 16.29 pow1.cs Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      pow1.cs
//
//   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
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class pow1
  {
    public static void Main(string[] args)
    {
      using (Model M = new Model("pow1"))
      {

        Variable x  = M.Variable("x", 3, Domain.Unbounded());
        Variable x3 = M.Variable();
        Variable x4 = M.Variable();

        // Create the linear constraint
        double[] aval = new double[] {1.0, 1.0, 0.5};
        M.Constraint(Expr.Dot(x, aval), Domain.EqualsTo(2.0));

        // Create the exponential conic constraint
        // 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));      

        // Set the objective function
        double[] cval = new double[] {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
        double[] solx = x.Level();
        Console.WriteLine("x,y,z = {0}, {1}, {2}", solx[0], solx[1], solx[2]);
      }
    }
  }
}

primal_svm.cs

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

using System;
using System.Collections.Generic;

using mosek.fusion;

namespace mosek.fusion.example
{
  public class primal_svm
  {
    public static void Main(string[] args)
    {
      int nc = 10;
      int m = 50;
      int n = 3;
      int seed = 1;
      double variance = 1.0;
      double mean = 1.0;

      Random gen = new Random(seed);

      int nump = gen.Next(m);

      double [] y = new double[m];

      for (int i = 0; i < nump; i++)  y[i] = 1.0;
      for (int i = nump; i < m; i++) y[i] = -1.0;

      double [,] X = new double[m, n];

      for (int i = 0; i < nump; i++)
        for (int j = 0; j < n; j++) {
          double xx = Math.Sqrt( -2.0 * Math.Log(gen.NextDouble()) );
          double yy = Math.Sin( 2.0 * Math.PI * gen.NextDouble() );
          X[i, j] = (xx * yy) * variance + mean;
        }
      Model M = new Model("Primal SVM");
      try {
        Console.WriteLine("Number of data    : {0}\n", m);
        Console.WriteLine("Number of features: {0}\n", n);

        Variable w =  M.Variable( n, Domain.Unbounded() );
        Variable t =  M.Variable( 1, Domain.Unbounded() );
        Variable b =  M.Variable( 1, Domain.Unbounded() );
        Variable xi = M.Variable( m, Domain.GreaterThan(0.0) );

        M.Constraint( Expr.Add(Expr.MulElm( y, Expr.Sub( Expr.Mul(X, w), Var.Repeat(b, m) ) ), xi ),
                      Domain.GreaterThan( 1.0 ) );
        M.Constraint( Expr.Vstack(1.0, t, w) , Domain.InRotatedQCone());

        Console.WriteLine(" c     | b      | w");
        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();

          Console.Write("{0} | {1} |", c, b.Level()[0] );
          for (int j = 0; j < n; j++)
            Console.Write(" {0}", w.Level()[j] );
          Console.WriteLine();
        }
      }
      finally {
        M.Dispose();
      }
    }
  }
}

qcqp_sdo_relaxation.cs

Listing 16.31 qcqp_sdo_relaxation.cs Click here to download.
/*
* Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
*
* File:      qcqp_sdo_relaxation.cs
*
* Purpose:   Demonstrate how to use SDP to solve
*            convex relaxation of a mixed-integer QCQP
**/
using System;

using mosek.fusion;
using mosek;

namespace mosek.fusion.example
{
  public class qcqp_sdo_relaxation
  {
    public static void Main(string[] args)
    {
      Random rnd = new Random();
      // problem dimensions
      int n = 20;
      int m = 2 * n;

      // problem data
      double[] A = new double[m * n];
      double[] b = new double[m];
      double[] c = new double[n];
      double[] P = new double[n * n];
      double[] q = new double[n];

      for (int j = 0; j < n; j++) {
        for (int i = 0; i < m; i++) {
          double u1 = rnd.NextDouble(); //these are uniform(0,1) random doubles
          double u2 = rnd.NextDouble();
          A[i * n + j] = Math.Sqrt(-2.0 * Math.Log(u1)) * Math.Sin(2.0 * Math.PI * u2);
        }
        c[j] = rnd.NextDouble();
      }

      // P = A^T A
      LinAlg.syrk(mosek.uplo.lo, mosek.transpose.yes,
                  n, m, 1.0, A, 0.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
      LinAlg.gemv(mosek.transpose.no, n, n, -1.0, P, c, 0.0, q);
      LinAlg.gemv(mosek.transpose.no, m, n,  1.0, A, c, 0.0, b);

      // Solve the problems
      Model M = miqcqp_sdo_relaxation(n, Matrix.Dense(n, n, P), q);
      Model Mint = int_least_squares(n, Matrix.Dense(n, m, A).Transpose(), b);
      M.Solve();
      Mint.Solve();

      // rounded and optimal solution
      double[] xRound = M.GetVariable("Z").Slice(new int[] {0, n}, new int[] {n, n + 1}).Level();
      for (int i = 0; i < n; i++) xRound[i] = Math.Round(xRound[i]);
      double[] yRound = (double[]) b.Clone();
      double[] xOpt   = Mint.GetVariable("x").Level();
      double[] yOpt   = (double[]) b.Clone();
      LinAlg.gemv(mosek.transpose.no, m, n,  1.0, A, xRound, -1.0, yRound);        // Ax_round-b
      LinAlg.gemv(mosek.transpose.no, m, n,  1.0, A, xOpt, -1.0, yOpt);            // Ax_opt-b

      Console.WriteLine(M.GetSolverDoubleInfo("optimizerTime") + " " + Mint.GetSolverDoubleInfo("optimizerTime"));
      Console.WriteLine(Math.Sqrt(LinAlg.dot(m, yRound, yRound)) + " " +
                        Math.Sqrt(LinAlg.dot(m, yOpt, yOpt)));
    }

    public static Model miqcqp_sdo_relaxation(int n, Matrix P, double[] q) {
      Model M = new Model();

      Variable Z = M.Variable("Z", Domain.InPSDCone(n + 1));

      Variable X = Z.Slice(new int[] {0, 0}, new int[] {n, n});
      Variable x = Z.Slice(new int[] {0, n}, new int[] {n, n + 1});

      M.Constraint( Expr.Sub(X.Diag(), x), Domain.GreaterThan(0.0) );
      M.Constraint( Z.Index(n, n), Domain.EqualsTo(1.0) );

      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|
    public static Model int_least_squares(int n, Matrix A, double[] b) {
      Model M = new Model();

      Variable x = M.Variable("x", n, Domain.Integral(Domain.Unbounded()));
      Variable 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;
    }
  }
}

reoptimization.cs

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

  File :      reoptimization.cs

  Purpose:   Demonstrates how to solve a  linear
             optimization problem using the MOSEK API
             and modify and re-optimize the problem.
*/
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class reoptimization
  {
    public static void printsol(double[] x) {
      Console.Write("x = ");
      foreach(double val in x)
          Console.Write(String.Format("{0:F4} ", val));
      Console.WriteLine();
    }

    public static void Main(string[] args)
    {
      using (Model M = new Model())
      {
        double[]   c = new double[] { 1.5, 2.5, 3.0 };
        double[,]  A = new double[,] { {2, 4, 3},
                                       {3, 2, 3},
                                       {2, 3, 2} };
        double[]   b = new double[] { 100000.0, 50000.0, 60000.0 };
        int   numvar = c.Length;
        int   numcon = b.Length;

        // Create a model and input data
        Variable x = M.Variable(numvar, Domain.GreaterThan(0.0));
        Constraint 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 element of 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
        Variable x3 = M.Variable(Domain.GreaterThan(0.0));
        Variable xNew = Var.Vstack(x, x3);
        // Add to the exising constraint
        con.Update(Expr.Mul(x3, new double[]{4, 0, 1}), x3);
        // Change the objective to include x3
        M.Objective(ObjectiveSense.Maximize, Expr.Dot(new double[]{1.5,2.5,3.0,1.0}, xNew));
        M.Solve();
        printsol(xNew.Level());

        /**************** Add a new constraint *****************************/
        Constraint con2 = M.Constraint(Expr.Dot(xNew, new double[]{1, 2, 1, 1}), Domain.LessThan(30000.0));
        M.Solve();
        printsol(xNew.Level());

        /**************** Change constraint bounds *****************************/
        // Assemble all constraints we previously defined into one
        Constraint cAll = Constraint.Vstack(con, con2);
        // Change bounds by effectively updating fixed terms with the difference
        cAll.Update(new double[]{20000, 10000, 10000, 8000});
        M.Solve();
        printsol(xNew.Level());
      }
    }
  }
}

response.cs

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

  File:      response.cs

  Purpose:   This example demonstrates proper response handling
             for problems solved with the interior-point optimizers.
*/
using System;
using mosek;
using mosek.fusion;
using System.Text;

namespace mosek.fusion.example
{
  public class response
  {
    // Set up a small artificial conic problem to test with
    public static void SetupExample(Model M) {
      Variable x = M.Variable("x", 3, Domain.GreaterThan(0.0));
      Variable y = M.Variable("y", 3, Domain.Unbounded());
      Variable z1 = Var.Vstack(y.Index(0), x.Slice(0, 2));
      Variable z2 = Var.Vstack(y.Slice(1, 3), x.Index(2));
      M.Constraint("lc", Expr.Dot(new double[] {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));
    }

    public static void Main(String[] argv) {
      Model M = new Model();

      // (Optional) set a log stream
      // M.SetLogHandler(Console.Out);
      
      // (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);

        Variable x = M.GetVariable("x");
        long xsize = x.GetSize();
        double[] xVal = x.Level();
        Console.Write("Optimal value of x = ");
        for(int i = 0; i < xsize; ++i)
          Console.Write(xVal[i] + " ");
        Console.WriteLine("\nOptimal primal objective: " + M.PrimalObjValue());
        // .. continue analyzing the solution

      }
      catch (OptimizeError e)
      {
        Console.WriteLine("Optimization failed. Error: " + e.ToString());
      }
      catch (SolutionError)
      {
        // The solution with at least the expected status was not available.
        // We try to diagnoze why.
        Console.WriteLine("Requested solution was not available.");
        ProblemStatus prosta = M.GetProblemStatus();
        switch(prosta)
        {
          case ProblemStatus.DualInfeasible:
            Console.WriteLine("Dual infeasibility certificate found.");
            break;

          case ProblemStatus.PrimalInfeasible:
            Console.WriteLine("Primal infeasibility certificate found.");
            break;

          case ProblemStatus.Unknown:
            // The solutions status is unknown. The termination code
            // indicates why the optimizer terminated prematurely.
            Console.WriteLine("The solution status is unknown.");
            StringBuilder symname = new StringBuilder();
            StringBuilder desc = new StringBuilder();
            Env.getcodedesc((rescode)M.GetSolverIntInfo("optimizeResponse"), symname, desc);
            Console.WriteLine("  Termination code: {0} {1}", symname, desc);
            break;

          default:
            Console.WriteLine("Another unexpected problem status: " + prosta);
            break;
        }
      }
      catch (Exception e)
      {
        Console.WriteLine("Unexpected error: " + e.ToString());
      }

      M.Dispose();
    }
  }
}

sdo1.cs

Listing 16.34 sdo1.cs Click here to download.
//
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      sdo1.cs
//
//  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)
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class sdo1
  {
    public static void Main(string[] args)
    {
      using (Model M  = new Model("sdo1"))
      {
        // Setting up the variables
        Variable X  = M.Variable("X", Domain.InPSDCone(3));
        Variable x  = M.Variable("x", Domain.InQCone(3));

        // Setting up constant coefficient matrices
        Matrix C  = Matrix.Dense ( new double[,] { {2, 1, 0}, {1, 2, 1}, {0, 1, 2}} );
        Matrix A1 = Matrix.Eye(3);
        Matrix 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();

        Console.WriteLine("[{0}]", (new Utils.StringBuffer()).A(X.Level()).ToString());
        Console.WriteLine("[{0}]", (new Utils.StringBuffer()).A(x.Level()).ToString());
      }
    }
  }
}

sdo2.cs

Listing 16.35 sdo2.cs Click here to download.
/*
   Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
 
   File :      sdo2.cs
 
   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.
*/
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class sdo2
  {
    public static void Main(string[] args)
    {

      // Sample data in sparse, symmetric triplet format
      int[]    C1_k = {0, 2};
      int[]    C1_l = {0, 2};
      double[] C1_v = {1, 6};
      int[]    A1_k = {0, 2, 0, 2};
      int[]    A1_l = {0, 0, 2, 2};
      double[] A1_v = {1, 1, 1, 2};
      int[]    C2_k = {0, 1, 0, 1, 2};
      int[]    C2_l = {0, 0, 1, 1, 2};
      double[] C2_v = {1, -3, -3, 2, 1};
      int[]    A2_k = {1, 0, 1, 3};
      int[]    A2_l = {0, 1, 1, 3};
      double[] A2_v = {1, 1, -1, -3};
      double b = 23;
      double k = -3;

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

      using (Model M  = new Model("sdo2"))
      {
        // Two semidefinite variables
        Variable X1 = M.Variable(Domain.InPSDCone(3));
        Variable 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(new int[] {0,1}), Domain.LessThan(k));

        // Solve
        M.SetLogHandler(Console.Out);
        M.Solve();

        // Print solution
        Console.WriteLine("Solution (vectorized):");
        Console.WriteLine("[{0}]", (new Utils.StringBuffer()).A(X1.Level()).ToString());
        Console.WriteLine("[{0}]", (new Utils.StringBuffer()).A(X2.Level()).ToString());
      }
    }
  }
}

sdo3.cs

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

  File :      sdo3.cs

  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.
*/

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class sdo3
  {

    // A helper method computing a semidefinite slice of a 3-dim variable
    public static Variable Slice(Variable X, int d, int j) {
      return
        X.Slice(new int[] {j,0,0}, new int[] {j+1,d,d})
         .Reshape(new int[] {d,d});
    }

    public static void Main(string[] args)
    {

      // Sample input data
      int n = 100;
      int d = 4;
      int k = 3;
      double[] b = {9,10,11};
      double[][,] A = new double[n*k][,];
      var rand = new Random();
      for(int i=0; i<n*k; i++) {
        A[i] = new double[d,d];
        for(int s1=0; s1<d; s1++)
          for(int s2=0; s2<=s1; s2++)
            A[i][s1,s2] = A[i][s2,s1] = rand.NextDouble();
      }

      // Create a model with n semidefinite variables od dimension d x d
      using (Model M  = new Model("sdo3"))
      {
        Variable X = M.Variable(Domain.InPSDCone(d, n));

        // Pick indexes of diagonal entries for the objective
        int[,] alldiag = new int[d*n,3];
        for(int j=0; j<n; j++) for(int s=0; s<d; s++) {
          alldiag[j*d+s,0] = j; 
          alldiag[j*d+s,1] = alldiag[j*d+s,2] = s;
        }
        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++) {
          Expression[] addlist = new Expression[n];
          for(int j=0; j<n; j++) 
            addlist[j] = Expr.Dot(A[i*n+j], Slice(X, d, j));
          M.Constraint(Expr.Add(addlist), Domain.GreaterThan(b[i]));
        }

        // Solve
        M.SetLogHandler(Console.Out);           // Add logging
        M.WriteTask("sdo3.ptf");                // Save problem in readable format
        M.Solve();

        // Get results. Each variable is a slice of X
        Console.WriteLine("Contributing variables:");
        for(int j=0; j<n; j++) {
          double[] Xj = Slice(X, d, j).Level();
          double maxval = 0;
          for(int s=0; s<d*d; s++) maxval = Math.Max(maxval, Xj[s]);
          if (maxval > 1e-6) {
            Console.WriteLine("X" + j + "=");
            for(int s1=0; s1<d; s1++) {
              for(int s2=0; s2<d; s2++)
                Console.Write(Xj[s1*d+s1] + "  ");    
              Console.WriteLine();
            }
          }
        }
      } 
    }
  }
}

sospoly.cs

Listing 16.37 sospoly.cs Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      sospoly.cs
//
//  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.
using mosek.fusion;
using System;

namespace mosek.fusion.example
{
  class sospoly
  {
    // Creates a Hankel matrix of dimension n+1, where
    // H_lk = a if l+k=i, and 0 otherwise.
    public static Matrix 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
    public static void nn_inf(Model M, Variable x)
    {
      int m = (int)x.GetSize() - 1;
      int n = (m / 2); // degree of polynomial is 2*n

      //assert(m == 2*n)

      // Setup variables
      Variable 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)
    public static void nn_semiinf(Model M, Variable x)
    {
      int n = (int)x.GetSize() - 1;
      int n1 = n / 2,
          n2 = (n - 1) / 2;

      // Setup variables
      Variable X1 = M.Variable(Domain.InPSDCone(n1 + 1));
      Variable 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]
    public static void nn_finite(Model M, Variable x, double a, double b)
    {
      //assert(a < b)
      int m = (int)x.GetSize() - 1;
      int n = m / 2;

      if (m == 2 * n)
      {
        Variable X1 = M.Variable(Domain.InPSDCone(n + 1));
        Variable 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 X1 = M.Variable(Domain.InPSDCone(n + 1));
        Variable 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).
    public static Variable diff(Model M, Variable x)
    {
      int n = (int)x.GetSize() - 1;
      Variable u = M.Variable(n, Domain.Unbounded());
      M.Constraint(Expr.Sub(u, Expr.MulElm(Matrix.Dense(n, 1, Range(1.0, n + 1)), x.Slice(1, n + 1))), Domain.EqualsTo(0.0));
      return u;
    }

    public static double[] fitpoly(double[,] data, int n)
    {
      using (var M = new Model("smooth poly"))
      {
        int m = data.GetLength(0);
        double[] Adata = new double[m * (n + 1)];
        for (int j = 0, k = 0; j < m; ++j)
          for (int i = 0; i < n + 1; ++i, ++k)
            Adata[k] = Math.Pow(data[j, 0], i);

        Matrix A = Matrix.Dense(m, n + 1, Adata);
        double[] b = new double[m];
        for (int j = 0; j < m; ++j)
          b[j] = data[j, 1];

        Variable x  = M.Variable("x", n + 1, Domain.Unbounded());
        Variable z  = M.Variable("z", 1,   Domain.Unbounded());
        Variable dx = diff(M, x);

        M.Constraint(Expr.Mul(A, x), Domain.EqualsTo(b));

        // z - f'(t) >= 0, for all t \in [a, b]
        Variable 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[m - 1, 0]);

        // f'(t) + z >= 0, for all t \in [a, b]
        Variable lb = M.Variable(n, Domain.Unbounded());
        M.Constraint(Expr.Sub(lb,
                              Expr.Vstack(Expr.Add(z, dx.Index(0)), dx.Slice(1, n).AsExpr())),
                     Domain.EqualsTo(0.0));
        nn_finite(M, lb, data[0, 0], data[m - 1, 0]);

        M.Objective(ObjectiveSense.Minimize, z);
        M.Solve();
        return x.Level();
      }
    }
    public static void Main(string[] argv)
    {

      double[,] data =
      { { -1.0, 1.0 },
        {  0.0, 0.0 },
        {  1.0, 1.0 }
      };

      double[] x2 = fitpoly(data, 2);
      double[] x4 = fitpoly(data, 4);
      double[] x8 = fitpoly(data, 8);

      Console.WriteLine("fitpoly(data,2) -> {0}", arrtostr(x2));
      Console.WriteLine("fitpoly(data,4) -> {0}", arrtostr(x4));
      Console.WriteLine("fitpoly(data,8) -> {0}", arrtostr(x8));
    }



    // Some utility functions to make things look nicer
    public static int[] Range (int start, int stop, int step)
    {
      int len;
      if (start < stop && step > 0)
        len = 1 + (stop - start - 1) / step;
      else if (stop < start && step < 0)
        len = 1 + (start - stop - 1) / (- step);
      else
        len = 0;
      int[] res = new int[len];
      for (int i = 0, v = start; i < len; ++i, v += step)
        res[i] = v;
      return res;
    }
    public static int[] Range (int start, int stop) { return Range(start, stop, 1); }
    public static int[] Range (int stop) { return Range(0, stop, 1); }
    public static double[] Range (double start, double stop, double step)
    {
      int len;
      if (start < stop && step > 0.0)
        len = 1 + (int)((stop - start - 1) / step);
      else if (stop < start && step < 0)
        len = 1 + (int)((start - stop - 1) / (- step));
      else
        len = 0;
      double[] res = new double[len];
      double v = start;
      for (int i = 0; i < len; ++i, v += step)
        res[i] = v;
      return res;
    }
    public static double[] Range (double start, double stop) { return Range(start, stop, 1.0); }
    public static double[] Range (double stop) { return Range(0.0, stop, 1.0); }

    public static double[] Fill (int num, double v) { double[] r = new double[num]; for (int i = 0; i < num; ++i) r[i] = v; return r;}

    private static string arrtostr(double[] a)
    {
      var c = new System.Globalization.CultureInfo("en-US");
      var b = new System.Text.StringBuilder("[");
      if (a.Length > 0) b.Append(String.Format(c, "{0:e3}", a[0]));
      for (int i = 1; i < a.Length; ++i) b.Append(",").Append(String.Format(c, "{0:e3}", a[i]));
      b.Append("]");
      return b.ToString();
    }
  }
}

sudoku.cs

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

using mosek.fusion;

namespace mosek.fusion.example
{
  public class sudoku
  {
    public static void Main(string[] args)
    {
      int m = 3;
      int n = m * m;

      //fixed cells in human readable (i.e. 1-based) format
      int[,] hr_fixed = {
        {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 nf = hr_fixed.Length / 3;
      int[,] fixed_cells = new int[nf, 3];
      for (int i = 0; i < nf; i++)
        for (int d = 0; d < m; d++)
          fixed_cells[i, d] = hr_fixed[i, d] - 1;

      using (Model M = new Model("SUDOKU")) {
        M.SetLogHandler(Console.Out);
        Variable x = M.Variable(new int[] { 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.0));

        //each number must appears 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++)
            {
              Variable block = x.Slice(new int[] { i * m, j * m, k },
                                       new int[] { (i + 1)*m, (j + 1)*m, k + 1 });
              M.Constraint(Expr.Sum(block), Domain.EqualsTo(1.0));
            }

        M.Constraint(x.Pick(fixed_cells), Domain.EqualsTo(1.0));

        M.Solve();

        //print the solution, if any...
        if (M.GetPrimalSolutionStatus() == SolutionStatus.Optimal)
          print_solution(m, x);
        else
          Console.WriteLine("No solution found!\n");
      }
    }

    public static void print_solution(int m, Variable x)
    {
      int n = m * m;
      Console.WriteLine("\n");
      try
      {
        for (int i = 0; i < n; i++)
        {
          String s = "";
          for (int j = 0; j < n; j++)
          {
            if (j % m == 0) s += " |";
            for (int k = 0; k < n; k++)
            {
              double[] sol = x.Index(new int[] { i, j, k }).Level();
              if (sol[0] > 0.5)
              {
                s += " " + (k + 1);
                break;
              }
            }
          }
          Console.WriteLine(s + " |");
          if ((i + 1) % m == 0)
            Console.WriteLine("\n");
        }
      }
      catch (SolutionError se) { }
    }
  }
}

total_variation.cs

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

namespace mosek.fusion.example
{
  public class total_variation
  {
    public static Model total_var(int n, int m) {
      Model M = new Model("TV");

      Variable u = M.Variable("u", new int[] {n + 1, m + 1}, Domain.InRange(0.0, 1.0));
      Variable t = M.Variable("t", new int[] {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 sigma = M.Parameter("sigma");
      Parameter f = M.Parameter("f", n, m);

      Variable ucore = u.Slice(new int[] {0, 0}, new int[] {n, m});

      Expression deltax = Expr.Sub( u.Slice( new int[] {1, 0}, new int[] {n + 1, m} ), ucore );
      Expression deltay = Expr.Sub( u.Slice( new int[] {0, 1}, new int[] {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;
    }

    public static void Main(string[] args)
    {
      Random randGen = new Random(0);
      
      int n = 100;
      int m = 200;
      double[] sigmas = { 0.0004, 0.0005, 0.0006 };

      // Create a parametrized model with given shape
      Model M = total_var(n, m);
      Parameter sigma = M.GetParameter("sigma");
      Parameter f     = M.GetParameter("f");
      Variable  ucore = M.GetVariable("u").Slice(new int[] {0,0}, new int[] {n,m});

      // Example: Linear signal with Gaussian noise    
      double[,] signal = new double[n,m];
      double[,] noise  = new double[n,m];
      double[,] fVal   = new double[n,m];
      double[,] sol    = new double[n,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] = Math.Sqrt(-2.0 * Math.Log(randGen.NextDouble())) * 
                     Math.Sin(2.0 * Math.PI * randGen.NextDouble()) *
                     0.08;
        fVal[i,j] = Math.Max( Math.Min(1.0, signal[i,j] + noise[i,j]), .0 );
      }
     
      // Set value for f
      f.SetValue(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
        double[] 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
        // ...

        Console.WriteLine("rel_sigma = " + sigmas[iter] + " total_var = " + M.PrimalObjValue());
      }

      M.Dispose();
    }
  }
}

tsp.cs

Listing 16.40 tsp.cs Click here to download.
//
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      tsp.cs
//
//  Purpose: Demonstrates a simple technique to the TSP
//           usign the Fusion API.
//
using System;
using System.Collections.Generic;
using System.Collections;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class tsp {
    public static void Main(string [] args) {
      int []A_i = new int[] {0, 1, 2, 3, 1, 0, 2, 0};
      int []A_j = new int[] {1, 2, 3, 0, 0, 2, 1, 3};
      double []C_v =  new double[] {1.0, 1.0, 1.0, 1.0, 0.1, 0.1, 0.1, 0.1};

      int n = 4;
      Matrix adj = Matrix.Sparse(n, n, A_i, A_j, 1.0);
      Matrix costs = Matrix.Sparse(n, n, A_i, A_j, C_v);

      tsp_fusion(n, adj, costs , true, true);
      tsp_fusion(n, adj, costs , true, false);
    }

    public static void tsp_fusion(int n, Matrix A, Matrix C,
                                  bool remove_loops,
                                  bool remove_2_hop_loops) {
      Model M = new Model();

      Variable x = M.Variable(new int[] {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_2_hop_loops)
        M.Constraint(Expr.Add(x, x.Transpose()), Domain.LessThan(1.0));

      if (remove_loops)
        M.Constraint(x.Diag(), Domain.EqualsTo(0.0));

      for (int it = 0; true; it++)
      {
        M.Solve();
        try
        {
          Console.WriteLine("\nit #{0} - solution cost: {1}", it, M.PrimalObjValue());
        } catch (SolutionError e) {}

        ArrayList cycles = new ArrayList();

        for (int i = 0; i < n; i++) {
          for (int j = 0; j < n; j++) {
            try {
              if ( x.Level()[i * n + j] <= 0.5) continue;
            } catch (SolutionError e) {}

            bool found = false;
            foreach (ArrayList c in cycles) {
              foreach (int[] a in c)
                if ( i == a[0] || j == a[0] ||
                     i == a[1] || j == a[1]  )
                {
                  c.Add( new int[] {i, j} );
                  found = true;
                  break;
                }
              if (found == true) break;
            }

            if (found != true) {
              ArrayList ll = new ArrayList();
              ll.Add( new int[] {i, j} );
              cycles.Add(ll);
            }
          }
        }

        Console.WriteLine("\ncycles:");
        foreach (ArrayList c in cycles) {
          foreach (int[] a in c)
            Console.Write("[{0},{1}] - ", a[0], a[1]);
          Console.Write("\n");
        }

        if (cycles.Count == 1) break;

        foreach (ArrayList c in cycles) {
          int [,] cc = new int[c.Count, 2];
          int i = 0;
          foreach (int[] a in c) {
            cc[i, 0] = a[0];
            cc[i, 1] = a[1];
            i++;
          }
          M.Constraint(Expr.Sum(x.Pick(cc)), Domain.LessThan( 1.0 * c.Count - 1 ));
        }
      }

      try {
        Console.WriteLine("\nsolution:");
        for (int i = 0; i < n; i++) {
          for (int j = 0; j < n; j++)
            Console.Write(" {0} ", ((int) x.Level()[i * n + j]));
          Console.Write("\n");
        }
      } catch (SolutionError e) {}
    }
  }
}