16 List of examples

List of examples shipped in the distribution of Fusion API for Java:

Table 16.1 List of distributed examples

File

Description

TrafficNetworkModel.java

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

alan.java

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

baker.java

A small bakery revenue maximization linear problem

breaksolver.java

Shows how to break a long-running task

callback.java

An example of data/progress callback

ceo1.java

A simple conic exponential problem

cqo1.java

A simple conic quadratic problem

diet.java

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

duality.java

Shows how to access the dual solution

elastic.java

Linear regression with elastic net. Demonstrates model parametrization.

facility_location.java

Demonstrates a small one-facility location problem (CQO)

gp1.java

A simple geometric program (GP) in conic form

lo1.java

A simple linear problem

logistic.java

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

lownerjohn_ellipsoid.java

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

lpt.java

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

mico1.java

A simple mixed-integer conic problem

milo1.java

A simple mixed-integer linear problem

mioinitsol.java

A simple mixed-integer linear problem with an initial guess

modelLib.java

Library of implementations of basic functions

nearestcorr.java

Solves the nearest correlation matrix problem (SDO, CQO)

opt_server_sync.java

Uses MOSEK OptServer to solve an optimization problem synchronously

parameters.java

Shows how to set optimizer parameters and read information items

portfolio_1_basic.java

Portfolio optimization - basic Markowitz model

portfolio_2_frontier.java

Portfolio optimization - efficient frontier

portfolio_3_impact.java

Portfolio optimization - market impact costs

portfolio_4_transcost.java

Portfolio optimization - transaction costs

portfolio_5_card.java

Portfolio optimization - cardinality constraints

pow1.java

A simple power cone problem

primal_svm.java

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

qcqp_sdo_relaxation.java

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

reoptimization.java

Demonstrate how to modify and re-optimize a linear problem

response.java

Demonstrates proper response handling

sdo1.java

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

sdo2.java

A simple semidefinite problem with two matrix variables

sdo3.java

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

sospoly.java

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

sudoku.java

A SUDOKU solver (MIP)

total_variation.java

Demonstrates how to solve a total variation problem (CQO)

tsp.java

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

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

package com.mosek.fusion.examples;
import mosek.fusion.*;

public class TrafficNetworkModel {

  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) throws SolutionError {
    
    try(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;
      int[]    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;


      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.dot(e, 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();
      
      return x.level();
    }
  }

  public static void main(String[] args)
  throws SolutionError {
    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;


    double[]flow = solve(n, source_idx, sink_idx,
                         arc_i, arc_j,
                         arc_sens,
                         arc_cap,
                         arc_base,
                         T);

    System.out.println("Optimal flow:");
    for (int i = 0; i < arc_i.length; ++i)
      System.out.println("\tflow node" + arc_i[i] + "->node" + arc_j[i] + " = " + flow[arc_i[i] * n + arc_j[i]]);
  }
}

alan.java

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

package com.mosek.fusion.examples;
import mosek.fusion.*;

public class alan {
  /////////////////////////////////////////////////////////////////////
  // Problem data.
  // Security names
  private static final String[]
    securities = { "hardware", "software", "show-biz", "t-bills" };
  // Two sample mean returns on securities
  private static final double[]
    mean1      = {        9.0,        7.0,       11.0,       5.0 };
  private static final double[]
    mean2      = {        8.0,        9.0,       12.0,       7.0 };

  // Target mean return
  private static final double
    target     = 10.0;

  private static final int numsec = securities.length;

  // Factor of covariance matrix.
  private static final double[][] U_data = {
    { 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 final Matrix U = Matrix.dense(U_data);

  public static void main(String[] args)
  throws SolutionError {
    // Create a parametrized model
    Model M = new Model("alan");
    try {

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

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

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

      // Solve with two values of mean
      solve(M, mean, x, mean1);
      solve(M, mean, x, mean2);
    } finally {
      M.dispose();
    }
  }

  // Solve the model given some expected return vector
  public static void solve(Model M, Parameter mean, Variable x, double[] meanVal)
  throws SolutionError {
    mean.setValue(meanVal);

    System.out.print("Solve with mean = ");
    for (int i = 0; i < numsec; ++i)
      System.out.printf(" %e", meanVal[i]);
    System.out.print("\n");

    M.solve();

    double[] solx = x.level();
    System.out.print("Solution = ");
    for (int i = 0; i < numsec; ++i)
      System.out.printf(" %e", solx[i]);
    System.out.print("\n");
  }

}

baker.java

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

package com.mosek.fusion.examples;

import mosek.fusion.*;
import java.io.PrintWriter;

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)
  throws SolutionError {
    Matrix recipe = Matrix.dense(recipe_data);
    Model M = new Model("Recipe");
    try {
      // "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(new PrintWriter(System.out));

      // We solve and fetch the solution:
      M.solve();
      double[] res = production.level();
      System.out.println("Solution:");
      for (int i = 0; i < res.length; ++i) {
        System.out.println(" Number of " + productnames[i] + " : " + res[i]);
      }
      System.out.println(" Revenue : $" +
                         (res[0] * revenue[0] + res[1] * revenue[1]));
    } finally {
      M.dispose();
    }
  }
}

breaksolver.java

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


package com.mosek.fusion.examples;
import mosek.fusion.*;
import java.util.Random;

public class breaksolver {
  private static void nshuffle(Random R, int[] a, int n) {
    for (int i = 0; i < n; ++i) {
      int idx = (int)(R.nextDouble() * (a.length - i)) + i;
      int tmp = a[i]; a[i] = a[idx]; a[idx] = tmp;
    }
  }

  public static void main(String[] args)
  throws java.lang.InterruptedException, SolutionError {
    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);

    System.out.println("Build problem...");
    try (Model M = new Model("SolveBinary")) {
      M.setLogHandler(new java.io.PrintWriter(System.out));

      Variable x = M.variable("x", n, Domain.binary());
      M.objective(ObjectiveSense.Minimize, Expr.sum(x));

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

      for (int i = 0; i < m; ++i) {
        nshuffle(R, idxs, p);
        M.constraint(Expr.sum(x.pick(java.util.Arrays.copyOf(idxs, p))),
                     Domain.equalsTo(p / 2));
      }

      System.out.println("Start thread...");
      Thread T = new Thread() { public void run() { M.solve(); } };
      T.start();

      long T0 = System.currentTimeMillis();
      while (true) {
        if (System.currentTimeMillis() - T0 > timeout * 1000) {
          System.out.println("Solver terminated due to timeout!");
          M.breakSolver();
          T.join();
          break;
        }
        if (! T.isAlive()) {
          System.out.println("Solver terminated before anything happened!");
          T.join();
          break;
        }
      }
    } finally {
      System.out.println("End.");
    }
  }
}

callback.java

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

   File :      callback.java

   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.
*/
package com.mosek.fusion.examples;

import mosek.*;
import mosek.fusion.*;
import java.util.Formatter;
import java.util.*;
import java.lang.*;

public class callback {
  private static DataCallback makeUserCallback(final double maxtime,
                                               final Model mod) {
    return new DataCallback() {
      public int callback(callbackcode caller,
                          double[]     douinf,
                          int[]        intinf,
                          long[]       lintinf) 
      {
        double opttime = 0.0;
        int itrn;
        double pobj, dobj, stime;

        Formatter f = new Formatter(System.out);
        switch (caller) {
          case begin_intpnt:
            f.format("Starting interior-point optimizer\n");
            break;
          case intpnt:
            itrn    = intinf[iinfitem.intpnt_iter.value      ];
            pobj    = douinf[dinfitem.intpnt_primal_obj.value];
            dobj    = douinf[dinfitem.intpnt_dual_obj.value  ];
            stime   = douinf[dinfitem.intpnt_time.value      ];
            opttime = douinf[dinfitem.optimizer_time.value   ];
            f.format("Iterations: %-3d\n", itrn);
            f.format("  Elapsed time: %6.2f(%.2f)\n", opttime, stime);
            f.format("  Primal obj.: %-18.6e  Dual obj.: %-18.6e\n", pobj, dobj);
            break;
          case end_intpnt:
            f.format("Interior-point optimizer finished.\n");
            break;
          case begin_primal_simplex:
            f.format("Primal simplex optimizer started.\n");
            break;
          case update_primal_simplex:
            itrn    = intinf[iinfitem.sim_primal_iter.value  ];
            pobj    = douinf[dinfitem.sim_obj.value          ];
            stime   = douinf[dinfitem.sim_time.value         ];
            opttime = douinf[dinfitem.optimizer_time.value   ];

            f.format("Iterations: %-3d\n", itrn);
            f.format("  Elapsed time: %6.2f(%.2f)\n", opttime, stime);
            f.format("  Obj.: %-18.6e\n", pobj );
            break;
          case end_primal_simplex:
            f.format("Primal simplex optimizer finished.\n");
            break;
          case begin_dual_simplex:
            f.format("Dual simplex optimizer started.\n");
            break;
          case update_dual_simplex:
            itrn    = intinf[iinfitem.sim_dual_iter.value    ];
            pobj    = douinf[dinfitem.sim_obj.value          ];
            stime   = douinf[dinfitem.sim_time.value         ];
            opttime = douinf[dinfitem.optimizer_time.value   ];
            f.format("Iterations: %-3d\n", itrn);
            f.format("  Elapsed time: %6.2f(%.2f)\n", opttime, stime);
            f.format("  Obj.: %-18.6e\n", pobj);
            break;
          case end_dual_simplex:
            f.format("Dual simplex optimizer finished.\n");
            break;
          case begin_bi:
            f.format("Basis identification started.\n");
            break;
          case end_bi:
            f.format("Basis identification finished.\n");
            break;
          default:
        }
        System.out.flush();
        if (opttime >= maxtime) {
          f.format("MOSEK is spending too much time. Terminate it.\n");
          System.out.flush();
          return 1;
        }
        return 0;
      }
    };
  }

  private static Progress makeProgressCallback() {
    return new Progress() {
      public int progress(callbackcode caller) 
      {
        // Process the caller code here
        return 0;
      }
    };
  }
  

  public static void main(String[] args) {
    String slvr     = "intpnt";
    if (args.length < 1) {
      System.out.println("Usage: ( psim | dsim | intpnt ) ");
    }

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

    // Create a random 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();

    try(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");

      double maxtime = 0.07;

      M.setDataCallbackHandler( makeUserCallback(maxtime, M) );

      M.solve();
    } catch (mosek.Exception e) {
      System.out.println ("An error/warning was encountered");
      System.out.println (e.toString());
      throw e;
    }
  }
}

ceo1.java

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

package com.mosek.fusion.examples;
import mosek.fusion.*;

public class ceo1 {
  public static void main(String[] args)
  throws SolutionError {
    Model M = new Model("ceo1");
    try {
      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 conic exponential 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();
      System.out.printf("x1,x2,x3 = %e, %e, %e\n", solx[0], solx[1], solx[2]);

      // Get conic solution of expc
      double[] expclvl = expc.level();
      double[] expcsn  = expc.dual();
      
      System.out.printf("expc levels = %e", expclvl[0]);
      for (int i = 1; i < expclvl.length; ++i)
        System.out.printf(", %e", expclvl[i]);
      System.out.print("\n");

      System.out.printf("expc dual conic var levels = %e", expcsn[0]);
      for (int i = 1; i < expcsn.length; ++i)
        System.out.printf(", %e", expcsn[i]);
      System.out.print("\n");

    } finally {
      M.dispose();
    }
  }
}

cqo1.java

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

package com.mosek.fusion.examples;
import mosek.fusion.*;

public class cqo1 {
  public static void main(String[] args)
  throws SolutionError {
    Model M = new Model("cqo1");
    try {
      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, Domain.inQCone());
      Constraint qc2 = M.constraint("qc2", z2, Domain.inRotatedQCone());

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

      // Solve the problem
      M.writeTask("cqo1.task");
      M.solve();


      // Get the linear solution values
      double[] solx = x.level();
      double[] soly = y.level();
      System.out.printf("x1,x2,x3 = %e, %e, %e\n", solx[0], solx[1], solx[2]);
      System.out.printf("y1,y2,y3 = %e, %e, %e\n", soly[0], soly[1], soly[2]);

      // Get conic solution of qc1
      double[] qc1lvl = qc1.level();
      double[] qc1sn  = qc1.dual();
      
      System.out.printf("qc1 levels = %e", qc1lvl[0]);
      for (int i = 1; i < qc1lvl.length; ++i)
        System.out.printf(", %e", qc1lvl[i]);
      System.out.print("\n");

      System.out.printf("qc1 dual conic var levels = %e", qc1sn[0]);
      for (int i = 1; i < qc1sn.length; ++i)
        System.out.printf(", %e", qc1sn[i]);
      System.out.print("\n");

    } finally {
      M.dispose();
    }
  }
}

diet.java

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

public class diet extends Model {
  public class DataException extends RuntimeException {
    public DataException(String msg) { super(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) {
    super(name);
    boolean 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[] dailyPurchase()
  throws SolutionError {
    return dailyPurchase.level();
  }

  public double[] dailyNutrients()
  throws SolutionError {
    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,    .6,     6,   55.8,      .2,    0,      0,     0},  // margarine
      { 7.4,   448,  16.4,    19,   28.1,      .8,   10.3,    4,     0},  // cheese
      {15.7,   661,   1,      48,     0,      9.6,    8.1,  471,     0},  // peanut butter
      {41.7,     0,   0,       0,     .2,     0,       .5,    5,     0},  // lard
      { 2.2,   333,    .2,   139,  169.2,     6.4,   50.8,  316,   525},  // liver
      { 4.4,   249,    .3,    37,    0,      18.2,    3.6,   79,     0},  // porkroast
      { 5.8,   705,   6.8,    45,    3.5,     1,      4.9,  209,     0},  // salmon
      { 2.4,   138,   3.7,    80,   69,       4.3,    5.8,   37,   862},  // greenbeans
      { 2.6,   125,   4,      36,    7.2,     9,      4.5,   26,  5369},  // cabbage
      { 5.8,   166,   3.8,    59,   16.6,     4.7,    5.9,   21,  1184},  // onions
      {14.3,   336,   1.8,   118,    6.7,    29.4,    7.1,  198,  2522},  // potatoes
      { 1.1,   106,   0.0,   138,  918.4,     5.7,   13.8,   33,  2755},  // spinach
      { 9.6,   138,   2.7,    54,  290.7,     8.4,    5.4,   83,  1912},  // sweet potatos
      { 8.5,    87,   1.7,   173,   86.8,     1.2,    4.3,   55,    57},  // peaches
      {12.8,    99,   2.5,   154,   85.7,     3.9,    4.3,   65,   257},  // prunes
      {17.4,  1055,   3.7,   459,    5.1,    26.9,   38.2,   93,     0},  // limabeans
      {26.9,  1691,  11.4,   792,    0,      38.4,   24.6,  217,     0}   // navybeans
    };

    double[] daily_allowance =
    {   3.,     70.,   .8,    12.,   5.,      1.8,    2.7,   18.,   75. };
    diet M = new diet("Stinglers Diet Problem",
                      foods,
                      nutrients,
                      daily_allowance,
                      nutritive_value);
    try {
      M.solve();
      try {
        System.out.println("Solution:");
        double[] x = M.dailyPurchase();
        for (int i = 0; i < x.length; ++i)
          System.out.println(foods[i] + " : $" + x[i]);

        System.out.println("Nutrients:");
        double[] y = M.dailyNutrients();
        for (int i = 0; i < y.length; ++i)
          System.out.println(nutrients[i] + " : " + y[i] + " " + nutrient_unit[i] +
                             " (" + daily_allowance[i] + ")");
      } catch (SolutionError e) {
        System.out.println("Solution error : " + e.getMessage());
      }
    } finally {
      M.dispose();
    }
  }
}

duality.java

Listing 16.9 duality.java Click here to download.
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      duality.java
//
// Purpose: Show how to read the dual value of a constraint.
package com.mosek.fusion.examples;
import mosek.fusion.*;

public class duality {
  public static void main(String[] args)
  throws SolutionError {
    double[][] A =  { { -0.5, 1.0 }  };
    double[]   b =  { 1.0 };
    double[]   c =  { 1.0, 1.0 };

    Model M = new Model("duality");
    try {
      Variable x = M.variable("x", 2, Domain.greaterThan(0.0));

      Constraint con = M.constraint(Expr.sub(Expr.mul(Matrix.dense(A), x), b), Domain.equalsTo(0.0));

      M.objective("obj", ObjectiveSense.Minimize, Expr.dot(c, x));

      M.solve();
      double[] xsol = x.level();
      double[] ssol = x.dual();
      double[] ysol = con.dual();

      System.out.printf("x1,x2,s1,s2,y = %e, %e, %e, %e, %e\n", xsol[0], xsol[1], ssol[0], ssol[1], ysol[0]);
    } finally {
      M.dispose();
    }
  }
}

elastic.java

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

package com.mosek.fusion.examples;
import mosek.fusion.*;

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() throws SolutionError {
    //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");

    try {
      // 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();
      System.out.printf("Objective %.5f, solution %.3f, %.3f\n", M.primalObjValue(), sol[0], sol[1]);

      // Increase lambda1
      lambda1.setValue(0.5);
      
      M.solve();
      sol = x.level();
      System.out.printf("Objective %.5f, solution %.3f, %.3f\n", 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();
      System.out.printf("Objective %.5f, solution %.3f, %.3f\n", M.primalObjValue(), sol[0], sol[1]);

      // And increase lamda2
      lambda2.setValue(1.4145);
      
      M.solve();
      sol = x.level();
      System.out.printf("Objective %.5f, solution %.3f, %.3f\n", M.primalObjValue(), sol[0], sol[1]);
    } 
    catch (SolutionError e) {
      System.out.println("Solution not available");
      throw e;
    }
    finally {
      M.dispose();
    }
  }

  public static void main(String[] argv) throws SolutionError {
    smallExample();
  }
}

facility_location.java

Listing 16.11 facility_location.java Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      facility_location.java
//
//  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
////
package com.mosek.fusion.examples;
import mosek.fusion.*;

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

  static int N = customerloc.numRows();
  public static void main(String[] args)
  throws SolutionError {
    Model M = new Model("FacilityLocation");
    // Variable holding the facility location
    try {
      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();

      double[] res = f.level();
      System.out.println("Facility location = " + res[0] + "," + res[1]);
    } finally {
      M.dispose();
    }
  }
}

gp1.java

Listing 16.12 gp1.java Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      gp1.java
//
//   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
//
package com.mosek.fusion.examples;
import mosek.fusion.*;
import java.lang.Math;

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.length;
    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)
  throws SolutionError
  {
    Model M = new Model("max_vol_box");
    try {
      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[] {Math.log(2.0/Aw), Math.log(2.0/Aw)});
      
      M.constraint(Expr.dot(new double[] {0,1,1}, xyz), Domain.lessThan(Math.log(Af)));
      M.constraint(Expr.dot(new double[] {1,-1,0}, xyz), Domain.inRange(Math.log(alpha),Math.log(beta)));
      M.constraint(Expr.dot(new double[] {0,-1,1}, xyz), Domain.inRange(Math.log(gamma),Math.log(delta)));
      
      M.setLogHandler(new java.io.PrintWriter(System.out));
      M.solve();
      
      double[] xyzVal = xyz.level();
      double[] hwdVal = new double[3];
      for(int i=0; i<3; i++) hwdVal[i] = Math.exp(xyzVal[i]);

      return hwdVal;
    } finally {
      M.dispose();
    }
  }

  public static void main(String[] args)
  throws SolutionError
  {
    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);

    System.out.format("h=%.4f w=%.4f d=%.4f\n", hwd[0], hwd[1], hwd[2]);
  }
}

lo1.java

Listing 16.13 lo1.java Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      lo1.java
//
//  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
////
package com.mosek.fusion.examples;
import mosek.fusion.*;

public class lo1 {
  public static void main(String[] args)
  throws SolutionError {
    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'
    try(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();
      System.out.printf("[x0,x1,x2,x3] = [%e, %e, %e, %e]\n", sol[0], sol[1], sol[2], sol[3]);
    }
  }
}

logistic.java

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

package com.mosek.fusion.examples;
import mosek.fusion.*;

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, 
                                         boolean[]  y,
                                         double     lamb)
  {
    int n = X.length;
    int d = X[0].length;       // 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)
  throws SolutionError {
    // Test: detect and approximate a circle using degree 2 polynomials
    int n = 30;
    double[][] X = new double[n*n][6];
    boolean[] Y = new boolean[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.setLogHandler(new java.io.PrintWriter(System.out));
    M.solve();
    for(int i=0;i<6;i++)
      System.out.println(theta.level()[i]);
  }
}

lownerjohn_ellipsoid.java

Listing 16.15 lownerjohn_ellipsoid.java Click here to download.
package com.mosek.fusion.examples;
/**
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:      lownerjohn_ellipsoid.java

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

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

import java.util.Arrays;
import mosek.fusion.*;

public class lownerjohn_ellipsoid {
  /**
    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. See [1,2] for more details.

   The convex set (a hypograph)

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

   can be modeled as the intersection of a semidefinite cone

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

   and a number of rotated quadratic cones and affine hyperplanes,

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


  /**
  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.
  */
  public static Object[] lownerjohn_inner(double[][] A, double[] b)
  throws SolutionError {
    Model M = new Model("lownerjohn_inner");
    try {
      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());

      // (b-Ad, AC) generate cones
      M.constraint("qc", Expr.hstack(Expr.sub(b, Expr.mul(A, d)), Expr.mul(A, C)),
                   Domain.inQCone());

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

      return new Object[] {C.level(), d.level()};
    } finally {
      M.dispose();
    }
  }

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

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

  minimizes the volume of the enclosing ellipsoid,

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

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

    maximize         t
    subject to       t       <= det(P)^(1/n)
                || P*xi - c ||_2 <= 1,  i=1,...,m
                P is PSD.
  */
  public static Object[] lownerjohn_outer(double[][] x)
  throws SolutionError {
    Model M = new Model("lownerjohn_outer");
    try {
      int m = x.length;
      int n = x[0].length;

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

      //(1, Px-c) in cone
      M.constraint("qc",
                   Expr.hstack(Expr.ones(m),
                               Expr.sub(Expr.mul(x, P),
                                        Var.reshape(Var.repeat(c, m), new int[] {m, n})
                                       )
                              ),
                   Domain.inQCone());

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

      return new Object[] {P.level(), c.level()};
    } finally {
      M.dispose();
    }
  }

  /**************************************************************************************/
  public static void main(String[] argv)
  throws SolutionError {
    //The list of vertices of a polygon
    final double[][] p = { {0., 0.}, {1., 3.}, {5.5, 4.5}, {7., 4.}, {7., 1.}, {3., -2.}  };
    double[][] A = new double[p.length][];
    double[]   b = new double[p.length];
    int n = p.length;

    //The h-representation of that polygon
    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];
    }

    Object[] ResInn = lownerjohn_inner(A, b);
    Object[] ResOut = lownerjohn_outer(p);
    double[] Ci = (double[]) ResInn[0];
    double[] di = (double[]) ResInn[1];
    double[] Po = (double[]) ResOut[0];
    double[] co = (double[]) ResOut[1];
    System.out.println("Inner:");
    System.out.println("  C = " + Arrays.toString(Ci));
    System.out.println("  d = " + Arrays.toString(di));
    System.out.println("Outer:");
    System.out.println("  P = " + Arrays.toString(Po));
    System.out.println("  c = " + Arrays.toString(co));
  }
}

lpt.java

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

import java.util.Random;

import java.util.*;
import java.io.*;
import java.lang.*;
import java.lang.Math.*;
import java.lang.Number.*;
import mosek.*;
import mosek.fusion.*;

public class lpt {
  public static void main(String [] args) {
    int n = 30;             //Number of tasks
    int m = 6;              //Number of processors

    double lb = 1.;         //The range of lengths of short tasks
    double ub = 5.;

    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);
    Arrays.sort(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.) );
    M.constraint( Expr.sub( Var.repeat(t, m), Expr.mul(x, T)), Domain.greaterThan(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(new PrintWriter(System.out));

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

    try {
      System.out.printf("initial solution: \n");
      for (int i = 0; i < m; i++) {
        System.out.printf("M %d [", i);
        for (int y = 0; y < n; y++)
          System.out.printf( "%d, ", (int) init[i * n + y]);
        System.out.printf("]\n");
      }
      System.out.print("MOSEK solution:\n");
      for (int i = 0; i < m; i++) {
        System.out.printf("M %d [", i);
        for (int y = 0; y < n; y++)
          System.out.printf( "%d, ", (int)x.index(i, y).level()[0]);
        System.out.printf("]\n");
      }
    } catch (SolutionError e) {}
  }
}

mico1.java

Listing 16.17 mico1.java Click here to download.
//
//    Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//    File:    mico1.java
//
//    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
//
package com.mosek.fusion.examples;
import mosek.fusion.*;

public class mico1 {
  public static void main(String[] args)
  throws SolutionError {

    Model M = new Model("mico1");

    try {
      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();

      System.out.printf("x, y = %f, %f\n", x.level()[0], y.level()[0]);

    } finally {
      M.dispose();
    }
  }
}

milo1.java

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

public class milo1 {
  public static void main(String[] args)
  throws SolutionError {
    double[][] A = {
      { 50.0, 31.0 },
      { 3.0,  -2.0 }
    };
    double[] c = { 1.0, 0.64 };

    Model M = new Model("milo1");
    try {
      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();
      System.out.printf("x1,x2 = %e, %e\n", sol[0], sol[1]);
      System.out.printf("MIP rel gap = %.2f (%f)\n",
                        M.getSolverDoubleInfo("mioObjRelGap"),
                        M.getSolverDoubleInfo("mioObjAbsGap"));
    } finally {
      M.dispose();
    }
  }
}

mioinitsol.java

Listing 16.19 mioinitsol.java Click here to download.
//
//    Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//    File:    mioinitsol.java
//
//    Purpose:  Demonstrates how to solve a small mixed
//              integer linear optimization problem.
//
package com.mosek.fusion.examples;
import mosek.fusion.*;

public class mioinitsol {
  public static void main(String[] args)
  throws SolutionError {
    int n = 4;
    double[] c = { 7.0, 10.0, 1.0, 5.0 };

    Model M = new Model("mioinitsol");
    try {
      Variable x = M.variable("x", n, Domain.greaterThan(0.0));
      x.slice(0,3).makeInteger();

      M.constraint(Expr.sum(x), Domain.lessThan(2.5));

      M.objective("obj", ObjectiveSense.Maximize, Expr.dot(c, x));

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

      M.solve();

      // Get the solution values
      double[] sol = x.level();
      System.out.printf("x = [");
      for (int i = 0; i < n; i++) {
        System.out.printf("%e, ", sol[i]);
      }
      System.out.printf("]\n");
      
      // Was the initial solution used?
      int constr = M.getSolverIntInfo("mioConstructSolution");
      double constrVal = M.getSolverDoubleInfo("mioConstructSolutionObj");
      System.out.println("Initial solution utilization: " + constr);
      System.out.println("Initial solution objective: " +  constrVal);

    } finally {
      M.dispose();
    }
  }
}

modelLib.java

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

  File :      modelLib.java

  Purpose: Library of simple building blocks in Mosek Fusion.
*/
package com.mosek.fusion.examples;
import mosek.*;
import mosek.fusion.*;

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(new java.io.PrintWriter(System.out));
    M.solve();
  }
    
  // Full test
  public static void main(String[] argv) {
    testModels();
    testExample();
  }
}

nearestcorr.java

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

   File:      nearestcorr.java

   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
*/
package com.mosek.fusion.examples;
import mosek.*;
import mosek.fusion.*;

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];
    int[] 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);
      }

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

  /* Nearest correlation with nuclear norm penalty */
  private static void nearestcorr_nn(Matrix A, double[] gammas, double[] res, int[] rank)
  throws SolutionError {
    int N = A.numRows();
    Model M = new Model("NucNorm");
    try {
      // Setup variables
      Variable t = M.variable("t", 1, Domain.unbounded());
      Variable X = M.variable("X", Domain.inPSDCone(N));
      Variable w = M.variable("w", N, Domain.greaterThan(0.0));

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

      // Trace(X)
      Expression TX = Expr.sum(X.diag());

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

        // Get the eigenvalues of X and approximate its rank
        double[] d = new double[N];
        LinAlg.syeig(mosek.uplo.lo, N, X.level(), d);
        int rnk = 0; for (int i = 0; i < d.length; ++i) if (d[i] > 1e-6) ++rnk;

        res[k] = t.level()[0];
        rank[k] = rnk;
      }
    } finally {
      M.dispose();
    }
  }

  private static void nearestcorr(Matrix A)
  throws SolutionError {
    int N = A.numRows();

    Model M = new Model("NearestCorrelation");
    try {
      // Setting up the variables
      Variable X = M.variable("X", Domain.inPSDCone(N));
      Variable 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
      System.out.println("X = \n" + mattostr(X.level(), N));
      System.out.println("t = \n" + mattostr(t.level(), N));

    } finally {
      M.dispose();
    }
  }

  public static void main(String[] argv)
  throws SolutionError {
    int N = 5;
    Matrix A = Matrix.dense(N, N,
                            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
                                         });
    double[] gammas = range(0.0, 2.0, 0.1);
    double[] res    = new double[gammas.length];
    int[]    rank   = new int[gammas.length];

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

    for (int i = 0; i < gammas.length; ++i)
      System.out.printf("gamma = %.1f, res=%.3e, rank = %d\n", gammas[i], res[i], rank[i]);
  }

  /* Utility functions */
  private static String mattostr(double[] a, int n) {
    StringBuilder b = new StringBuilder();
    java.util.Formatter f = new java.util.Formatter(b, java.util.Locale.US);

    if ( a.length == 1) {
      f.format("%f\n", a[0]);
      return b.toString();
    }

    for (int i = 0; i < a.length; ++i) {
      if (  i  % n == 0 )
        b.append("[");
      f.format(" %.3e ", a[i]);
      if ( ( i + 1 ) % n == 0 )
        b.append("]\n");
    }
    return b.toString();
  }


  private 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.java

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

package com.mosek.fusion.examples;
import mosek.fusion.*;

public class opt_server_sync {

  public static void main(String[] args) throws SolutionError {
    if (args.length<1) {
        System.out.println("Missing argument, syntax is:");
        System.out.println("   java com.mosek.fusion.examples.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(new java.io.PrintWriter(System.out));

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

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

    // Get the solution
    double[] solx = x.level();
    System.out.printf("x1,x2,x3 = %e, %e, %e\n", solx[0], solx[1], solx[2]);
  }
}

parameters.java

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

  File :      parameters.java

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

public class parameters {
    public static void main(String[] argv) {
        // Create the Model        
        Model M = new Model();
        System.out.println("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) {
            System.out.println("Wrong parameter value");
        }


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

        System.out.println("Get MOSEK information items");

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

        System.out.println("Time: " + tm);
        System.out.println("Iterations: " + it);
    }
}

portfolio_1_basic.java

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

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

  Purpose :   Implements a basic portfolio optimization model.
*/
package com.mosek.fusion.examples;

import mosek.fusion.*;
import java.io.FileReader;
import java.io.BufferedReader;
import java.util.ArrayList;

public class portfolio_1_basic {
  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)
  throws mosek.fusion.SolutionError {

    Model M = new Model("Basic Markowitz");
    try {
      // Redirect log output from the solver to stdout for debugging.
      // if uncommented.
      //M.setLogHandler(new java.io.PrintWriter(System.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());
    } finally {
      M.dispose();
    }
  }

  /*
    The example. Reads in data and solves the portfolio models.
   */
  public static void main(String[] argv)
  throws java.io.IOException,
         java.io.FileNotFoundException,
         mosek.fusion.SolutionError {

    int        n      = 3;
    double     w      = 1.0;
    double[]   mu     = {0.1073, 0.0737, 0.0627};
    double[]   x0     = {0.0, 0.0, 0.0};
    double[]   gammas = {0.035, 0.040, 0.050, 0.060, 0.070, 0.080, 0.090};
    double[][] GT     = {
      { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
      { 0.0              , 0.102863378954911  , -0.00222873156550421},
      { 0.0              , 0.0                ,  0.0338148677744977 }
    };


    System.out.println("\n-----------------------------------------------------------------------------------");
    System.out.println("Basic Markowitz portfolio optimization");
    System.out.println("-----------------------------------------------------------------------------------\n");
    for ( int i = 0; i < gammas.length; ++i) {
      double expret = BasicMarkowitz( n, mu, GT, x0, w, gammas[i]);
      System.out.format("Expected return: %.4e Std. deviation: %.4e\n",
                        expret,
                        gammas[i]);
    }
  }
}

portfolio_2_frontier.java

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

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

  Purpose :   Implements a basic portfolio optimization model.
              Computes points on the efficient frontier.
*/
package com.mosek.fusion.examples;

import mosek.fusion.*;
import java.io.FileReader;
import java.io.BufferedReader;
import java.util.ArrayList;

public class portfolio_2_frontier {
  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)
  throws mosek.fusion.SolutionError {

    Model M = new Model("Efficient frontier");
    try {
      //M.setLogHandler(new java.io.PrintWriter(System.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 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];
      }
    } finally {
      M.dispose();
    }
  }

  /*
    The example. Reads in data and solves the portfolio models.
   */
  public static void main(String[] argv)
  throws java.io.IOException,
         java.io.FileNotFoundException,
         mosek.fusion.SolutionError {

    int        n      = 3;
    double     w      = 1.0;
    double[]   mu     = {0.1073, 0.0737, 0.0627};
    double[]   x0     = {0.0, 0.0, 0.0};
    double[][] GT     = {
      { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
      { 0.0              , 0.102863378954911  , -0.00222873156550421},
      { 0.0              , 0.0                ,  0.0338148677744977 }
    };


    // 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);
    System.out.println("\n-----------------------------------------------------------------------------------");
    System.out.println("Efficient frontier") ;
    System.out.println("-----------------------------------------------------------------------------------\n");
    System.out.format("%-12s  %-12s  %-12s\n", "alpha", "return", "variance") ;
    for (int i = 0; i < frontier_mux.length; ++i)
      System.out.format("\t%-12.4f  %-12.4e  %-12.4e\n", alphas[i], frontier_mux[i], frontier_s[i]);

  }
}

portfolio_3_impact.java

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

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

  Purpose :   Implements a basic portfolio optimization model
              with transaction costs of the form x^(3/2).
*/
package com.mosek.fusion.examples;

import mosek.fusion.*;
import java.io.FileReader;
import java.io.BufferedReader;
import java.util.ArrayList;

public class portfolio_3_impact {
  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)
  throws mosek.fusion.SolutionError {
    Model M = new Model("Markowitz portfolio with market impact");
    try {
      //M.setLogHandler(new java.io.PrintWriter(System.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)
        System.arraycopy(x.level(), 0, xsol, 0, n);
      if (tsol != null)
        System.arraycopy(t.level(), 0, tsol, 0, n);
    } finally {
      M.dispose();
    }
  }


  /*
    The example. Reads in data and solves the portfolio models.
   */
  public static void main(String[] argv)
  throws java.io.IOException,
         java.io.FileNotFoundException,
         mosek.fusion.SolutionError {

    int        n      = 3;
    double     w      = 1.0;
    double[]   mu     = {0.1073, 0.0737, 0.0627};
    double[]   x0     = {0.0, 0.0, 0.0};
    double[]   gammas = {0.035, 0.040, 0.050, 0.060, 0.070, 0.080, 0.090};
    double[][] GT     = {
      { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
      { 0.0              , 0.102863378954911  , -0.00222873156550421},
      { 0.0              , 0.0                ,  0.0338148677744977 }
    };


    // 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);
    System.out.println("\n-----------------------------------------------------------------------------------");
    System.out.println("Markowitz portfolio optimization with market impact cost");
    System.out.println("-----------------------------------------------------------------------------------\n");
    System.out.format("Expected return: %.4e Std. deviation: %.4e Market impact cost: %.4e\n",
                      dot(mu, x),
                      gamma,
                      dot(m, t));
  }
}

portfolio_4_transcost.java

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

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

  Purpose :   Implements a basic portfolio optimization model
              with fixed setup costs and transaction costs
              as a mixed-integer problem.
*/
package com.mosek.fusion.examples;

import mosek.fusion.*;
import java.io.FileReader;
import java.io.BufferedReader;
import java.util.ArrayList;

public class portfolio_4_transcost {
  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)
  throws mosek.fusion.SolutionError {

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

    Model M = new Model("Markowitz portfolio with transaction costs");
    try {
      //M.setLogHandler(new java.io.PrintWriter(System.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();
    } finally {
      M.dispose();
    }
  }


  /*
    The example. Reads in data and solves the portfolio models.
   */
  public static void main(String[] argv)
  throws java.io.IOException,
         java.io.FileNotFoundException,
         mosek.fusion.SolutionError {

    int        n      = 3;
    double     w      = 1.0;
    double[]   mu     = {0.1073, 0.0737, 0.0627};
    double[]   x0     = {0.0, 0.0, 0.0};
    double[]   gammas = {0.035, 0.040, 0.050, 0.060, 0.070, 0.080, 0.090};
    double[][] GT     = {
      { 0.166673333200005, 0.0232190712557243 ,  0.0012599496030238 },
      { 0.0              , 0.102863378954911  , -0.00222873156550421},
      { 0.0              , 0.0                ,  0.0338148677744977 }
    };

    double[] f = new double[n]; java.util.Arrays.fill(f, 0.01);
    double[] g = new double[n]; java.util.Arrays.fill(g, 0.001);
    double gamma = gammas[0];

    System.out.println("\n-----------------------------------------------------------------------------------");
    System.out.println("Markowitz portfolio optimization with transaction cost");
    System.out.println("-----------------------------------------------------------------------------------\n");

    double[] x = new double[n];
    x = MarkowitzWithTransactionsCost(n, mu, GT, x0, w, gamma, f, g);
    System.out.println("Optimal portfolio: \n");
    for ( int i = 0; i < x.length; ++i)
      System.out.format("\tx[%-2d]  %-12.4e\n", i, x[i]);

  }
}

portfolio_5_card.java

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

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

  Description :  Implements a basic portfolio optimization model
                 with cardinality constraints on number of assets traded.
*/
package com.mosek.fusion.examples;

import mosek.fusion.*;
import java.io.FileReader;
import java.io.BufferedReader;
import java.util.ArrayList;

public class portfolio_5_card {
  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: Maximum number of assets on 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)
  throws mosek.fusion.SolutionError {

    // 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;
    }
 
    Model M = new Model("Markowitz portfolio with cardinality bounds");
    try {
      //M.setLogHandler(new java.io.PrintWriter(System.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);
      M.solve();

      // 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;
    } finally {
      M.dispose();
    }
  }

  /*
    The example. Reads in data and solves the portfolio models.
   */
  public static void main(String[] argv)
  throws java.io.IOException,
         java.io.FileNotFoundException,
         mosek.fusion.SolutionError {

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

    double[][] results;
    System.out.println("\n-----------------------------------------------------------------------------------");
    System.out.println("Markowitz portfolio optimization with cardinality constraints");
    System.out.println("-----------------------------------------------------------------------------------\n");

    results = MarkowitzWithCardinality(n, mu, GT, x0, w, gamma, kValues);
    for(int k = 1; k <= n; k++) {
      System.out.format("Bound %d   Solution: ", k);
      for (int i = 0; i < n; ++i)
        System.out.format("%-12.4e ", results[k-1][i]);
      System.out.println();
    }
  }
}

pow1.java

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

package com.mosek.fusion.examples;
import mosek.fusion.*;

public class pow1 {
  public static void main(String[] args)
  throws SolutionError {
    Model M = new Model("pow1");
    try {
      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 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();
      System.out.printf("x, y, z = %e, %e, %e\n", solx[0], solx[1], solx[2]);
    } finally {
      M.dispose();
    }
  }
}

primal_svm.java

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

package com.mosek.fusion.examples;

import java.util.Random;
import java.util.*;
import java.io.*;
import java.lang.*;
import java.lang.Math.*;
import java.lang.Number.*;

import mosek.fusion.*;

public class primal_svm {

  private static final int nc = 10;
  private static final int m = 50;
  private static final int n = 3;
  private static final int seed = 0;
  private static final double var = 1.;
  private static final double mean = 1.;

  public static void main(String[] args)
  throws SolutionError {

    Random gen = new Random(seed);

    int nump = gen.nextInt(m);
    int numm = m - nump;

    double [] y = new double[m];

    Arrays.fill(y, 0, nump, 1.);
    Arrays.fill(y, nump, m, -1.);

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

    for (int i = 0; i < nump; i++)
      for (int j = 0; j < n; j++)
        X[i][j] = gen.nextGaussian() + 1.;

    for (int i = nump; i < m; i++)
      for (int j = 0; j < n; j++)
        X[i][j] = gen.nextGaussian() - 1.;

    Model M = new Model("Primal SVM");
    try {
      System.out.format("Number of data    : %d\n", m);
      System.out.format("Number of features: %d\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.));

      M.constraint(
        Expr.add( xi ,
                  Expr.mulElm( y,
                               Expr.sub( Expr.mul(X, w), Var.repeat(b, m) )
                             )
                ),
        Domain.greaterThan( 1. )
      );

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

      System.out.println("   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();

        try {
          System.out.format("%4f | %8f", c, b.level()[0] );
          for (int j = 0; j < n; j++)
            System.out.format(" | %8f", w.level()[j] );

          System.out.print("\n");
        } catch (FusionException e) {}
      }
    }
    finally {
      M.dispose();
    }
  }
}

qcqp_sdo_relaxation.java

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

import java.util.*;
import java.lang.*;
import java.lang.Math.*;
import java.io.*;

import mosek.*;
import mosek.fusion.*;

public class qcqp_sdo_relaxation {
  public static void main(String [] args)
  throws SolutionError {
    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++)
        A[i * n + j] = rnd.nextGaussian();
      c[j] = rnd.nextDouble();
    }

    // P = A^T A
    LinAlg.syrk(mosek.uplo.lo, mosek.transpose.yes,
                n, m, 1.0, A, 0., P);
    for (int j = 0; j < n; j++) for (int i = j + 1; i < n; i++) P[i * n + j] = P[j * n + i];

    // q = -P c, b = A c
    LinAlg.gemv(mosek.transpose.no, n, n, -1.0, P, c, 0., q);
    LinAlg.gemv(mosek.transpose.no, m, n,  1.0, A, c, 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] = java.lang.Math.round(xRound[i]);
    double[] yRound = b.clone();
    double[] xOpt   = Mint.getVariable("x").level();
    double[] yOpt   = 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

    System.out.println(M.getSolverDoubleInfo("optimizerTime") + " " + Mint.getSolverDoubleInfo("optimizerTime"));
    System.out.println(java.lang.Math.sqrt(LinAlg.dot(m, yRound, yRound)) + " " +
                       java.lang.Math.sqrt(LinAlg.dot(m, yOpt, yOpt)));
  }

  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.) );
    M.constraint( Z.index(n, n), Domain.equalsTo(1.) );

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

    return M;
  }

  // A direct integer model for minimizing |Ax-b|
  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.java

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

  File :      reoptimization.java

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

public class reoptimization {
  public static void printsol(double[] x) {
    System.out.print("x = ");
    for(double val : x)
        System.out.printf("%.2f  ", val);
    System.out.println();   
  }

  public static void main(String[] args)
  throws SolutionError {
    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
    Model M = new Model();
    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.java

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

   File :      response.java

   Purpose :   This example demonstrates proper response handling
               for problems solved with the interior-point optimizers.
*/
package com.mosek.fusion.examples;
import mosek.*;
import mosek.fusion.*;

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
    /*
      task.set_Stream(
        streamtype.log,
        new mosek.Stream()
      { public void stream(String msg) { System.out.print(msg); }});
    */

    // (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();
      System.out.print("Optimal value of x = ");
      for(int i = 0; i < xsize; ++i)
        System.out.print(xVal[i] + " ");
      System.out.println("\nOptimal primal objective: " + M.primalObjValue());
      // .. continue analyzing the solution

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

        case PrimalInfeasible:
          System.out.println("Primal infeasibility certificate found.");
          break;

        case Unknown:
          // The solutions status is unknown. The termination code
          // indicates why the optimizer terminated prematurely.
          System.out.println("The solution status is unknown.");
          StringBuffer symname = new StringBuffer();
          StringBuffer desc = new StringBuffer();
          Env.getcodedesc(rescode.fromValue(M.getSolverIntInfo("optimizeResponse")), symname, desc);
          System.out.printf("   Termination code: %s %s\n", symname, desc);
          break;

        default:
          System.out.println("Another unexpected problem status: " + prosta);
      }
    }
    catch (java.lang.Exception e)
    {
      System.out.println("Unexpected error: " + e.toString());
    }

    M.dispose();
  }
}

sdo1.java

Listing 16.34 sdo1.java Click here to download.
//
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      sdo1.java
//
//  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)
//
package com.mosek.fusion.examples;
import mosek.fusion.*;

public class sdo1 {
  public static void main(String[] args) throws SolutionError {
    Model M  = new Model("sdo1");
    try {
      // 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();

      System.out.println(java.util.Arrays.toString( X.level() ));
      System.out.println(java.util.Arrays.toString( x.level() ));
    } finally {
      M.dispose();
    }
  }
}

sdo2.java

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

public class sdo2 {
  public static void main(String[] args) throws SolutionError {

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

    Model M  = new Model("sdo2");
    try {
      // 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(new java.io.PrintWriter(System.out));
      M.solve();

      // Print solution
      System.out.println("Solution (vectorized):");
      System.out.println(java.util.Arrays.toString( X1.level() ));
      System.out.println(java.util.Arrays.toString( X2.level() ));
    } finally {
      M.dispose();
    }
  }
}

sdo3.java

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

  File :      sdo3.java

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

package com.mosek.fusion.examples;
import mosek.fusion.*;

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) throws SolutionError {

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

    // Create a model with n semidefinite variables od dimension d x d
    Model M  = new Model("sdo3");
    try {
      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(new java.io.PrintWriter(System.out)); // Add logging
      M.writeTask("sdo3.ptf");                              // Save problem in readable format
      M.solve();

      // Get results. Each variable is a slice of X
      System.out.println("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) {
          System.out.println("X" + j + "=");
          for(int s1=0; s1<d; s1++) {
            for(int s2=0; s2<d; s2++)
              System.out.print(Xj[s1*d+s1] + "  ");    
            System.out.println();
          }
        }
      }
    } 
    finally {
      M.dispose();
    }
  }
}

sospoly.java

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

import mosek.fusion.*;

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), fill(i + 1, a));
    else
      return Matrix.sparse(n + 1, n + 1, range(n, i - n - 1, -1), range(i - n, n + 1), fill(2 * n + 1 - i, 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)
  throws SolutionError {
    Model M = new Model("smooth poly");

    try {
      int m = data.length;
      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[data.length - 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))),
                   Domain.equalsTo(0.0));
      nn_finite(M, lb, data[0][0], data[data.length - 1][0]);

      M.objective(ObjectiveSense.Minimize, z);
      M.solve();
      return x.level();
    } finally {
      M.dispose();
    }
  }
  public static void main(String[] argv)
  throws SolutionError {

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

    System.out.println("fitpoly(data,2) -> " + arrtostr(x2));
    System.out.println("fitpoly(data,4) -> " + arrtostr(x4));
    System.out.println("fitpoly(data,8) -> " + 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) {
    StringBuilder b = new StringBuilder("[");
    java.util.Formatter f = new java.util.Formatter(b, java.util.Locale.US);
    if (a.length > 0) f.format(", %.3e", a[0]);
    for (int i = 1; i < a.length; ++i) f.format(", %.3e", a[i]);
    b.append("]");
    return b.toString();
  }
}

sudoku.java

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

package com.mosek.fusion.examples;

import java.lang.*;
import java.util.*;
import java.io.*;
import mosek.fusion.*;

public class sudoku {
  //fixed cells in human readable (i.e. 1-based) format
  private static final 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}
  };

  public static void main(String[] args)
  throws SolutionError {

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

    int nf = hr_fixed.length;

    int [][] fixed = new int[nf][3];

    for (int i = 0; i < nf; i++)
      for (int d = 0; d < m; d++)
        fixed[i][d] =  hr_fixed[i][d] - 1;

    try(Model M = new Model("SUDOKU")) {
      M.setLogHandler(new PrintWriter(System.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.) );

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

      M.constraint( x.pick(fixed) , Domain.equalsTo(1.0) ) ;

      M.solve();

      //print the solution, if any...
      if ( M.getPrimalSolutionStatus() == SolutionStatus.Optimal )
        print_solution(m, x);
      else
        System.out.println("No solution found!\n");

    }
  }


  static void print_solution(int m, Variable x) {
    int n = m * m;
    System.out.println("\n");
    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++)
          try {
            double sol[] = x.index( new int[] {i, j, k}).level();
            if (sol[0] > 0.5) {
              s += " " + (k + 1);
              break;
            }
          } catch (SolutionError se) {}
      }
      System.out.println(s + " |");

      if ((i + 1) % m == 0)
        System.out.println("\n");
    }
  }
}

total_variation.java

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

package com.mosek.fusion.examples;
import java.util.Random;
import mosek.fusion.*;

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., 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)
  throws SolutionError {
    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] = randGen.nextGaussian() * 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
      // ...

      System.out.printf("rel_sigma = %f  total_var = %.3f\n", sigmas[iter], M.primalObjValue());
    }

    M.dispose();
  }
}

tsp.java

Listing 16.40 tsp.java Click here to download.
//
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      tsp.java
//
//  Purpose: Demonstrates a simple technique to the TSP
//           usign the Fusion API.
//
package com.mosek.fusion.examples;

import java.util.*;
import java.lang.*;
import java.io.*;
import java.lang.Math.*;
import mosek.*;
import mosek.fusion.*;


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., 1., 1., 1., 0.1, 0.1, 0.1, 0.1};

    int n = 4;

    Matrix costs = Matrix.sparse(n, n, A_i, A_j, C_v);
    Matrix adj   = Matrix.sparse(n, n, A_i, A_j, 1.);

    tsp_fusion(n, adj, costs, true, false);     //we do not remove loops of lenght 2
    tsp_fusion(n, adj, costs, true, true);      //we include all constraints
  }

  public static void tsp_fusion(int n,
                                Matrix A,
                                Matrix C,
                                boolean remove_loops,
                                boolean 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_loops)
      M.constraint(x.diag(), Domain.equalsTo(0.));

    if (remove_2_hop_loops)
      M.constraint(Expr.add(x, x.transpose()), Domain.lessThan(1.0));

    int it = 0;

    while (true) {
      it = it + 1;
      M.solve();

      try {
        System.out.printf("\nit #%d - solution cost: %f\n", it, M.primalObjValue());
      } catch (SolutionError e) {}

      ArrayList<ArrayList< int[] > > cycles = new ArrayList<ArrayList< int[] > >();

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

          boolean found = false;
          for (ArrayList<int[]> c : cycles) {
            for (int[] a : 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< int[] > ll = new ArrayList<int[]>();
            ll.add( new int[] {i, j} );
            cycles.add(ll);
          }
        }
      }

      System.out.printf("\ncycles:\n");
      for (ArrayList<int[]> c : cycles) {
        for (int[] a : c)
          System.out.printf("[%d,%d] - ", a[0], a[1]);
        System.out.printf("\n");
      }

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

      for (ArrayList<int []> c : cycles) {
        int [][] cc = new int[c.size()][2];
        c.toArray(cc);
        M.constraint(Expr.sum(x.pick(cc)), Domain.lessThan( 1.0 * c.size() - 1 ));
      }
    }

    try {
      System.out.printf("\nsolution:\n");
      for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++)
          System.out.printf(" %d ", ((int) x.level()[i * n + j]) );
        System.out.printf("\n");
      }
    } catch (SolutionError e) {}
  }
}