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

djc1.java

A simple problem with disjunctive constraints (DJC)

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

nearestcorr.java

Solves the nearest correlation matrix problem (SDO, CQO)

opt_server_sync.java

Uses MOSEK OptServer to solve an optimization problem synchronously

parallel.java

Demonstrates parallel optimization using a batch method in MOSEK

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

portfolio_6_factor.java

Portfolio optimization - factor model

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

djc1.java

Listing 16.9 djc1.java Click here to download.
////
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      djc1.java
//
//   Purpose: Demonstrates how to solve the problem with two disjunctions:
//
//      minimize    2x0 + x1 + 3x2 + x3
//      subject to   x0 + x1 + x2 + x3 >= -10
//                  (x0-2x1<=-1 and x2=x3=0) or (x2-3x3<=-2 and x1=x2=0)
//                  x0=2.5 or x1=2.5 or x2=2.5 or x3=2.5
////
package com.mosek.fusion.examples;
import mosek.fusion.*;

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

    try(Model M = new Model("djc1"))
    {
      // Create variable 'x' of length 4
      Variable x = M.variable("x", 4);

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

      // Second disjunctive constraint
      // Array of terms reading x_i = 2.5 for i = 0,1,2,3
      Term[] terms = new Term[4];
      for(int i = 0; i < 4; i++)
        terms[i] = DJC.term(x.index(i), Domain.equalsTo(2.5));
      // The disjunctive constraint from the array of terms
      M.disjunction(terms);

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

      // Objective
      M.objective(ObjectiveSense.Minimize, Expr.dot(new double[]{2,1,3,1}, x));

      // Useful for debugging
      M.writeTask("djc.ptf");                                 // Save to a readable file
      M.setLogHandler(new java.io.PrintWriter(System.out));   // Enable log output

      // Solve the problem
      M.solve();

      if (M.getPrimalSolutionStatus() == SolutionStatus.Optimal) {
        double[] sol = x.level();
        System.out.printf("[x0,x1,x2,x3] = [%e, %e, %e, %e]\n", sol[0], sol[1], sol[2], sol[3]);
      }
      else {
        System.out.printf("Anoter solution status");
      }
    }
  }
}

duality.java

Listing 16.10 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.11 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.12 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", new int[]{1, 2}, Domain.unbounded());
      // Variable defining the euclidian distances to each customer
      Variable d = M.variable("dist", new int[]{N, 1}, Domain.greaterThan(0.0));
      // Variable defining the x and y differences to each customer;
      Variable t = M.variable("t", new int[]{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.13 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.14 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.15 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.16 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 {
  /**
   Purpose: Models the hypograph of the n-th power of the
   determinant of a positive definite matrix. See [1,2] for more details.

   The convex set (a hypograph)

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

   can be modeled as the intersection of a semidefinite cone

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

   and a geometric mean cone

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

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

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


  /**
  The inner ellipsoidal approximation to a polytope

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

  maximizes the volume of the inscribed ellipsoid,

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

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

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

  which is equivalent to a mixed conic quadratic and semidefinite
  programming problem.
  */
  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.17 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.18 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.19 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.20 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 );

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

      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("Construct solution utilization: " + constr);
      System.out.println("Construct solution objective: " +  constrVal);
    } finally {
      M.dispose();
    }
  }
}

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 addr [certpath]");
        return;
    }

    String serveraddr = args[0];
    String tlscert = (args.length==2) ? args[1] : "";

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

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

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

parallel.java

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

   File:      parallel.java

   Purpose: Demonstrates parallel optimization using solveBatch()
*/
package com.mosek.fusion.examples;
import mosek.fusion.*;

public class parallel
{
  public static Model makeToyParameterizedModel() 
  {
    Model M = new Model();
    Variable x = M.variable("x",3);
    Parameter p = M.parameter("p");
    M.objective(ObjectiveSense.Maximize, Expr.sum(x));
    M.constraint(Expr.vstack(p,x), Domain.inQCone());
    return M;
  }

  /** Example of how to use Model.solveBatch()
   */
  public static void main(String[] argv) throws SolutionError
  {
    // Choose some sample parameters
    int n = 10;                 // Number of models to optimize
    int threadpoolsize = 4;     // Total number of threads available
    int threadspermodel = 1;    // Number of threads per each model

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

    // Set up n copies of the model with different data
    Model[] models = new Model[n];
    for(int i = 0; i < n ; i++)
    {
      models[i] = M.clone();
      models[i].getParameter("p").setValue(i+1);
      // We can set the number of threads individually per model
      models[i].setSolverParam("numThreads", threadspermodel);
    }

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

    // Access the soutions
    for(int i = 0; i < n; i++) 
      if (status[i] == SolverStatus.OK)
        System.out.printf("Model  %d:  Status %s   Solution Status %s   objective %.3f  time %.3f\n", 
          i, 
          status[i],
          models[i].getPrimalSolutionStatus(),
          models[i].primalObjValue(),
          models[i].getSolverDoubleInfo("optimizerTime"));
      else
        System.out.printf("Model  %d:  not solved\n", i);
  }
}

parameters.java

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

    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.26 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;
import java.lang.Math;

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

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

  }
}

portfolio_3_impact.java

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

    // Somewhat arbirtrary choice of m
    double[] m = new double[n]; for (int i = 0; i < n; ++i) m[i] = 0.01;
    double[] x = new double[n];
    double[] t = new double[n];
    double gamma = 0.36;

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

    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 = 0.36;

    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]);
    System.out.format("Expected return: %-12.4e", dot(mu, x));

  }
}

portfolio_5_card.java

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

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

portfolio_6_factor.java

Listing 16.30 portfolio_6_factor.java Click here to download.
/*
  File : portfolio_6_factor.java

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

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

import mosek.fusion.*;
import mosek.LinAlg;
import java.io.FileReader;
import java.io.BufferedReader;
import java.util.Arrays;

public class portfolio_6_factor 
{
  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;
  }

  // Create diagonal matrix from an array.
  public static double[][] diag_matrix(double[] x) 
  {
    int n = x.length;
    double[][] m = new double[n][n];  
    for (double[] row : m) Arrays.fill(row, 0.0);
    for (int i = 0; i < n; ++i) m[i][i] = x[i];
    return m;
  }

  public static double[][] transpose(double[][] m) 
  {
    int ni = m.length; 
    int nj = m[0].length;
    double[][] mt = new double[nj][ni];  
    
    for (int i = 0; i < ni; ++i) 
    {
      for (int j = 0; j < nj; ++j) 
      {
        mt[j][i] = m[i][j];
      }
    }
    return mt;
  }

  public static double[][] array_sqrt(double[][] m) 
  {
    int ni = m.length;
    int nj = m[0].length;
    double[][] sqrtm = new double[ni][nj];  
    
    for (int i = 0; i < ni; ++i) 
    {
      for (int j = 0; j < nj; ++j) 
      {
        sqrtm[i][j] = Math.sqrt(m[i][j]);
      }
    }
    return sqrtm;
  }

  // Create horizontal block matrix
  public static double[][] concat_h(double[][] ml, double[][] mr) 
  {
    int ni = ml.length;
    int njl = ml[0].length;
    int njr = mr[0].length;
    double[][] b = new double[ni][njl + njr];
    
    for (int i = 0; i < ni; ++i) 
    {
      System.arraycopy(ml[i], 0, b[i], 0, njl);
      System.arraycopy(mr[i], 0, b[i], njl, njr);
    }
    return b;
  }

  // Vectorize matrix (column-major order)
  public static double[] mat_to_vec_c(double[][] m) 
  {
    int ni = m.length;
    int nj = m[0].length;
    double[] c = new double[nj * ni];  
    
    for (int j = 0; j < nj; ++j) 
    {
      for (int i = 0; i < ni; ++i) 
      {
        c[j * ni + i] = m[i][j];
      }
    }
    return c;
  }

  // Reshape vector to matrix (column-major order)
  public static double[][] vec_to_mat_c(double[] c, int ni, int nj) 
  {
    double[][] m = new double[ni][nj];
    
    for (int j = 0; j < nj; ++j) 
    {
      for (int i = 0; i < ni; ++i) 
      {
        m[i][j] = c[j * ni + i];
      }
    }
    return m;
  }

  public static double[][] cholesky(double[][] m) 
  {
    int n = m.length;
    double[] vecs = mat_to_vec_c(m);
    LinAlg.potrf(mosek.uplo.lo, n, vecs);
    double[][] s = vec_to_mat_c(vecs, n, n);

    // Zero out upper triangular part (LinAlg.potrf does not use it, original matrix values remain there)
    for (int i = 0; i < n; ++i) 
    {
      for (int j = i+1; j < n; ++j) 
      {
        s[i][j] = 0.0;
      }
    }
    return s;
  }

  // Matrix multiplication
  public static double[][] matrix_mul(double[][] a, double[][] b) {
    int na = a.length;
    int nb = b[0].length;
    int k = b.length;

    double[] vecm = new double[na * nb];
    Arrays.fill(vecm, 0.0);
    LinAlg.gemm(mosek.transpose.no, mosek.transpose.no, na, nb, k, 1.0, mat_to_vec_c(a), mat_to_vec_c(b), 1.0, vecm);
    double[][] m = vec_to_mat_c(vecm, na, nb);     
    
    return m;
  }

  /*
  Purpose:
      Computes the optimal portfolio for a given risk using factor model approximation of the covariance.

  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 FactorMarkowitz
  ( int n,
    double[] mu,
    double[][] GT,
    double[] x0,
    double   w,
    double   gamma)
  throws mosek.fusion.SolutionError {

    Model M = new Model("Factor model 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      = 8;
    double     w      = 1.0;
    double[]   mu     = {0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379};
    double[]   x0     = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    // Factor exposure matrix
    double[][] B = 
    {
      {0.4256, 0.1869},
      {0.2413, 0.3877},
      {0.2235, 0.3697},
      {0.1503, 0.4612},
      {1.5325, -0.2633},
      {1.2741, -0.2613},
      {0.6939, 0.2372},
      {0.5425, 0.2116}
    };

    // Factor covariance matrix
    double[][] S_F = 
    {
      {0.0620, 0.0577},
      {0.0577, 0.0908}
    };

    // Specific risk components
    double[] resid_risk = {0.0720, 0.0508, 0.0377, 0.0394, 0.0663, 0.0224, 0.0417, 0.0459};
    double[][] S_theta = diag_matrix(resid_risk);
    double[][] P = cholesky(S_F);
    double[][] BP = matrix_mul(B, P);  
    double[][] G = concat_h(BP, array_sqrt(S_theta));
    double[][] GT = transpose(G);
    double[]   gammas = {0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48};

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

pow1.java

Listing 16.31 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.32 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.33 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.34 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.35 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.36 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.37 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.38 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.39 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.40 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.41 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.42 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) {}
  }
}