17 List of examples¶
List of examples shipped in the distribution of Optimizer API for Rust:
File |
Description |
|---|---|
A simple problem with one affine conic constraint (ACC) |
|
A simple problem with two affine conic constraints (ACC) |
|
Demonstrates the MOSEK interface to BLAS/LAPACK linear algebra routines |
|
An example of data/progress callback |
|
A simple conic exponential problem |
|
Implementation of a concurrent optimizer for linear and mixed-integer problems |
|
A simple conic quadratic problem |
|
A simple problem with disjunctive constraints (DJC) |
|
A simple example of how to repair an infeasible problem |
|
A simple geometric program (GP) in conic form |
|
A Hello World example |
|
A simple linear problem |
|
A simple linear problem |
|
Implements logistic regression and simple log-sum-exp (CEO) |
|
A simple mixed-integer conic problem |
|
A simple mixed-integer linear problem |
|
A simple mixed-integer linear problem with an initial guess |
|
Uses MOSEK OptServer to solve an optimization problem asynchronously |
|
Uses MOSEK OptServer to solve an optimization problem synchronously |
|
Demonstrates parallel optimization using a batch method in MOSEK |
|
Shows how to set optimizer parameters and read information items |
|
Shows how to obtain and analyze a primal infeasibility certificate |
|
Portfolio optimization - basic Markowitz model |
|
Portfolio optimization - efficient frontier |
|
Portfolio optimization - market impact costs |
|
Portfolio optimization - transaction costs |
|
Portfolio optimization - cardinality constraints |
|
Portfolio optimization - factor model |
|
A simple power cone problem |
|
A simple quadratically constrained quadratic problem |
|
A simple quadratic problem |
|
Demonstrate how to modify and re-optimize a linear problem |
|
Demonstrates proper response handling |
|
A simple semidefinite problem with one matrix variable and a quadratic cone |
|
A simple semidefinite problem with two matrix variables |
|
A simple semidefinite problem with an LMI using the SVEC domain. |
|
Sensitivity analysis performed on a small linear problem |
|
A simple I/O example: read problem from a file, solve and write solutions |
|
Demonstrates how to examine the quality of a solution |
|
Demonstrates solving a linear system with the basis matrix |
|
Demonstrates solving a general linear system |
|
Shows how to find a Cholesky factorization of a sparse matrix |
Additional examples can be found on the MOSEK website and in other MOSEK publications.
acc1.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : acc1.rs
//!
//! Purpose : Tutorial example for affine conic constraints.
//! Models the problem:
//!
//! maximize c^T x
//! subject to sum(x) = 1
//! gamma >= |Gx+h|_2
//!
extern crate mosek;
extern crate itertools;
use mosek::{Task,Objsense,Streamtype,Solsta,Soltype,Boundkey};
// Define problem data
#[allow(non_upper_case_globals)]
const n : i32 = 3;
#[allow(non_upper_case_globals)]
const k : i64 = 2;
#[allow(non_snake_case)]
fn main() -> Result<(),String> {
// Create a task
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
// Attach a printer to the task
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
// Create n free variables
task.append_vars(n)?;
let x : Vec<i32> = (0..n).collect();
task.put_var_bound_slice_const(0, n, Boundkey::FR, 0.0, 0.0)?;
// Set up the objective
let c = &[2.0, 3.0, -1.0];
task.put_obj_sense(Objsense::MAXIMIZE)?;
task.put_c_list(x.as_slice(), c)?;
// One linear constraint - sum(x) = 1
task.append_cons(1)?;
task.put_a_row(0,x.as_slice(), vec![1.0; n as usize].as_slice())?;
task.put_con_bound(0, Boundkey::FX, 1.0, 1.0)?;
// Append empty AFE rows for affine expression storage
task.append_afes(k + 1)?;
// G matrix in sparse form
let Gsubi : &[i64] = &[0, 0, 1, 1];
let Gsubj : &[i32] = &[0, 1, 0, 2];
let Gval = &[1.5, 0.1, 0.3, 2.1];
// Other data
let h = &[0.0, 0.1];
let gamma = 0.03;
// Construct F matrix in sparse form
let Fsubi : Vec<i64> = Gsubi.iter().map(|i| *i+1).collect(); // G will be placed from row number 1 in F
let Fsubj = Gsubj;
let Fval = Gval;
// Fill in F storage
task.put_afe_f_entry_list(Fsubi.as_slice(), Fsubj, Fval)?;
// Fill in g storage
task.put_afe_g(0, gamma)?;
task.put_afe_g_slice(1, k+1, h)?;
// Define a conic quadratic domain
let quadDom = task.append_quadratic_cone_domain(k + 1)?;
// Create the ACC
task.append_acc(quadDom, // Domain index
(0..k+1).collect::<Vec<i64>>().as_slice(), // Indices of AFE rows [0,...,k]
vec![0.0; (k+1) as usize].as_slice())?; // Ignored
// Solve and retrieve solution
let _ = task.optimize()?;
task.write_data("acc1.ptf")?;
let mut xx = vec![0.0; n as usize];
task.get_xx(Soltype::ITR,xx.as_mut_slice())?;
assert!(task.get_sol_sta(Soltype::ITR)? == Solsta::OPTIMAL);
println!("Solution: {:?}",xx);
// Demonstrate retrieving activity of ACC
let mut activity = vec![0.0; (k+1) as usize];
let mut doty = vec![0.0; (k+1) as usize];
task.evaluate_acc(Soltype::ITR,0,activity.as_mut_slice())?;
println!("Activity of ACC:: {:?}",activity);
// Demonstrate retrieving the dual of ACC
task.get_acc_dot_y(Soltype::ITR,0,doty.as_mut_slice())?;
println!("Dual of ACC:: {:?}",doty);
Ok(())
}
fn maxgap(a : &[f64], b : &[f64]) -> f64 {
a.iter().zip(b.iter()).map(|(&a,&b)| (a - b).abs()).max_by(|a,b| if a < b { std::cmp::Ordering::Less } else if b < a { std::cmp::Ordering::Greater } else {std::cmp::Ordering::Equal} ).unwrap()
}
fn dot(a : &[f64], b : &[f64]) -> f64 {
a.iter().zip(b.iter()).map(|(&a,&b)| (a * b)).sum()
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
acc2.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : acc2.rs
//!
//! Purpose : Tutorial example for affine conic constraints.
//! Models the problem:
//!
//! maximize c^T x
//! subject to sum(x) = 1
//! gamma >= |Gx+h|_2
//!
//! This version inputs the linear constraint as an affine conic constraint.
//!
extern crate mosek;
extern crate itertools;
use mosek::{Task,Objsense,Streamtype,Solsta,Soltype,Boundkey};
#[allow(non_upper_case_globals)]
fn main() -> Result<(),String> {
// Define problem data
const n : i32 = 3;
const k : i64 = 2;
// Create a task
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
// Attach a printer to the task
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
// Create n free variables
task.append_vars(n)?;
let x : Vec<i32> = (0..n).collect();
task.put_var_bound_slice_const(0, n, Boundkey::FR, 0.0,0.0)?;
// Set up the objective
let c = &[2.0, 3.0, -1.0];
task.put_obj_sense(Objsense::MAXIMIZE)?;
task.put_c_list(x.as_slice(), c)?;
// Set AFE rows representing the linear constraint
task.append_afes(1)?;
task.put_afe_f_row(0, x.as_slice(), vec![1.0; n as usize].as_slice())?;
task.put_afe_g(0, -1.0)?;
// Set AFE rows representing the quadratic constraint
task.append_afes(k + 1)?;
task.put_afe_f_row(2, // afeidx, row number
&[0i32, 1], // varidx, column numbers
&[1.5, 0.1])?; // values
task.put_afe_f_row(3, // afeidx, row number
&[0i32, 2], // varidx, column numbers
&[0.3, 2.1])?; // values
let h = &[0.0, 0.1];
let gamma = 0.03;
task.put_afe_g(1, gamma)?;
task.put_afe_g_slice(2, k+2, h)?;
// Define domains
let zero_dom = task.append_rzero_domain(1)?;
let quad_dom = task.append_quadratic_cone_domain(k + 1)?;
// Append affine conic constraints
task.append_acc(zero_dom, // Domain index
&[0i64], // Indices of AFE rows
&[0.0])?; // Ignored
task.append_acc(quad_dom, // Domain index
&[1i64,2,3], // Indices of AFE rows
&[0.0,0.0,0.0])?; // Ignored
// Solve and retrieve solution
let _ = task.optimize()?;
task.write_data("acc2.ptf")?;
let mut xx = vec![0.0; n as usize];
task.get_xx(Soltype::ITR,xx.as_mut_slice())?;
assert! (task.get_sol_sta(Soltype::ITR)? == Solsta::OPTIMAL);
println!("Solution: {:?}",xx);
// Demonstrate retrieving activity of ACC
let mut activity = vec![0.0; 3];
let mut doty = vec![0.0; 3];
task.evaluate_acc(Soltype::ITR,1,activity.as_mut_slice())?;
println!("Activity of quadratic ACC:: {:?}",activity);
// Demonstrate retrieving the dual of ACC
task.get_acc_dot_y(Soltype::ITR,1,doty.as_mut_slice())?;
println!("Dual of quadratic ACC:: {:?}",doty);
Ok(())
}
fn maxgap(a : &[f64], b : &[f64]) -> f64 {
a.iter().zip(b.iter()).map(|(&a,&b)| (a - b).abs()).max_by(|a,b| if a < b { std::cmp::Ordering::Less } else if b < a { std::cmp::Ordering::Greater } else {std::cmp::Ordering::Equal} ).unwrap()
}
fn dot(a : &[f64], b : &[f64]) -> f64 {
a.iter().zip(b.iter()).map(|(&a,&b)| (a * b)).sum()
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
blas_lapack.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : blas_lapack.rs
//!
//! Purpose: To demonstrate how to call BLAS/LAPACK routines for
//! which MOSEK provides simplified interfaces.
//!
extern crate mosek;
fn print_matrix(x : &[f64], r : i32, c : i32) {
let mut k = 0usize;
print!("[");
println!(" {:?}",&x[k..k+c as usize]); k += c as usize;
for _i in 1..r {
println!(" {:?}",&x[k..k+c as usize]);
}
println!("]");
}
#[allow(non_snake_case)]
fn main() -> Result<(),String> {
let n = 3i32;
let m = 2i32;
let k = 3i32;
let alpha = 2.0;
let beta = 0.5;
let x = &[1.0, 1.0, 1.0];
let mut y = vec![1.0, 2.0, 3.0];
let mut z = vec![1.0, 1.0];
let A = &[1.0, 1.0, 2.0, 2.0, 3.0, 3.0];
let B = &[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
let mut C = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
let mut D = vec![1.0, 1.0, 1.0, 1.0];
let mut Q = vec![1.0, 0.0, 0.0, 2.0];
let mut v = vec![0.0, 0.0, 0.0];
let mut xy : f64 = 0.0;
/* BLAS routines*/
println!("n={} m={} k={}", m, n, k);
println!("alpha={}", alpha);
println!("beta={}", beta);
mosek::dot(n,x,y.as_slice(),& mut xy)?;
println!("dot results = {}\n", xy);
print_matrix(x, 1, n);
print_matrix(y.as_slice(), 1, n);
mosek::axpy(n, alpha, x, y.as_mut_slice())?;
println!("axpy results is:\n");
print_matrix(y.as_slice(), 1, n);
mosek::gemv(mosek::Transpose::NO, m, n, alpha, A, x, beta, z.as_mut_slice())?;
println!("gemv results is:");
print_matrix(z.as_slice(), 1, m);
mosek::gemm(mosek::Transpose::NO, mosek::Transpose::NO, m, n, k, alpha, A, B, beta, C.as_mut_slice())?;
println!("gemm results is");
print_matrix(C.as_slice(), m, n);
mosek::syrk(mosek::Uplo::LO, mosek::Transpose::NO, m, k, 1., A, beta, D.as_mut_slice())?;
println!("syrk results is");
print_matrix(D.as_slice(), m, m);
/* LAPACK routines*/
mosek::potrf(mosek::Uplo::LO, m, Q.as_mut_slice())?;
println!("potrf results is");
print_matrix(Q.as_slice(), m, m);
mosek::syeig(mosek::Uplo::LO, m, Q.as_slice(), & mut v[0..m as usize])?;
println!("syeig results is");
print_matrix(v.as_slice(), 1, m);
mosek::syevd(mosek::Uplo::LO, m, Q.as_mut_slice(), &mut v[0..m as usize])?;
println!("syevd results is");
print_matrix(v.as_slice(), 1, m);
print_matrix(Q.as_slice(), m, m);
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
callback.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : callback.rs
//!
//! Purpose : To demonstrate how to use the progress
//! callback.
//!
//! Use this script as follows:
//!
//! callback psim 25fv47.mps
//! callback dsim 25fv47.mps
//! callback intpnt 25fv47.mps
//!
//! The first argument tells which optimizer to use
//! i.e. psim is primal simplex, dsim is dual simplex
//! and intpnt is interior-point.
extern crate mosek;
use std::env;
use mosek::{Task,Streamtype,Iparam,Optimizertype,Callbackcode,Dinfitem,Iinfitem};
const MAXTIME : f64 = 0.05;
fn callback(caller : i32, dinf : &[f64], iinf : &[i32], _linf : &[i64]) -> bool {
let mut opttime = 0.0;
match caller {
Callbackcode::BEGIN_INTPNT =>
println!("Starting interior-point optimizer"),
Callbackcode::INTPNT => {
let itrn = iinf[Iinfitem::INTPNT_ITER as usize];
let pobj = dinf[Dinfitem::INTPNT_PRIMAL_OBJ as usize];
let dobj = dinf[Dinfitem::INTPNT_DUAL_OBJ as usize];
let stime = dinf[Dinfitem::INTPNT_TIME as usize];
opttime = dinf[Dinfitem::OPTIMIZER_TIME as usize];
println!("Iterations: {:-3}",itrn);
println!(" Elapsed time: {:6.2}({:.2})",opttime, stime);
println!(" Primal obj.: {:-18.6e} Dual obj.: {:-18.6e}",pobj, dobj);
},
Callbackcode::END_INTPNT =>
println!("Interior-point optimizer finished."),
Callbackcode::BEGIN_PRIMAL_SIMPLEX =>
println!("Primal simplex optimizer started."),
Callbackcode::UPDATE_PRIMAL_SIMPLEX => {
let itrn = iinf[Iinfitem::SIM_PRIMAL_ITER as usize];
let pobj = dinf[Dinfitem::SIM_OBJ as usize];
let stime = dinf[Dinfitem::SIM_TIME as usize];
opttime = dinf[Dinfitem::OPTIMIZER_TIME as usize];
println!("Iterations: {:-3}",itrn);
println!(" Elapsed time: {:6.2}({:.2})",opttime, stime);
println!(" Obj.: {:-18.6e}",pobj);
},
Callbackcode::END_PRIMAL_SIMPLEX =>
println!("Primal simplex optimizer finished."),
Callbackcode::BEGIN_DUAL_SIMPLEX =>
println!("Dual simplex optimizer started."),
Callbackcode::UPDATE_DUAL_SIMPLEX => {
let itrn = iinf[Iinfitem::SIM_DUAL_ITER as usize];
let pobj = dinf[Dinfitem::SIM_OBJ as usize];
let stime = dinf[Dinfitem::SIM_TIME as usize];
opttime = dinf[Dinfitem::OPTIMIZER_TIME as usize];
println!("Iterations: {:-3}", itrn);
println!(" Elapsed time: {:6.2}({:.2})",opttime,stime);
println!(" Obj.: {:-18.6e}",pobj);
},
Callbackcode::END_DUAL_SIMPLEX =>
println!("Dual simplex optimizer finished."),
Callbackcode::NEW_INT_MIO => {
println!("New integer solution has been located.");
// let mut xx = vec![0.0; ]
// xx = task.get_xx(Soltype::ITG);
// println!(xx);
println!("Obj.: {}",dinf[Dinfitem::MIO_OBJ_INT as usize]);
}
_ => {
}
}
if opttime >= MAXTIME {
// mosek is spending too much time. Terminate it.
println!("Terminating.");
false
}
else {
true
}
}
enum FileOrText<'a> {
File(& 'a str),
Text(&'a str)
}
fn main() -> Result<(),String> {
let args : Vec<String> = env::args().collect();
if args.len() < 3 {
println!("Syntax: callback (psim|dsim|intpnt) FILENAME");
}
callbackmain(args[1].as_str(),FileOrText::File(args[2].as_str()))
}
fn callbackmain(which : &str, data : FileOrText) -> Result<(),String> {
/* Create the optimization task. */
let mut task = Task::new().unwrap();
match data {
FileOrText::Text(data) => { task.read_ptf_string(data)? },
FileOrText::File(fname) => { task.read_data(fname)? }
}
task.write_data("callback.ptf")?;
match which {
"psim" => task.put_int_param(Iparam::OPTIMIZER,Optimizertype::PRIMAL_SIMPLEX)?,
"dsim" => task.put_int_param(Iparam::OPTIMIZER,Optimizertype::DUAL_SIMPLEX)?,
"intpnt" => task.put_int_param(Iparam::OPTIMIZER,Optimizertype::INTPNT)?,
s => return Err(format!("Invalid argument '{}'",s))
}
/* Directs the log task stream to the 'printstr' function. */
task.with_stream_callback(
Streamtype::LOG,
& mut |msg| print!("{}",msg),
|task|
task.with_info_callback(
& mut callback,
|task|
task.optimize()
)
)?;
Result::Ok(())
}
#[cfg(test)]
mod tests {
const DFLT_FILE : &str = "Task
Objective
Maximize + 3 @x0 + @x1 + 5 @x2 + @x3
Constraints
@c0 [30] + 3 @x0 + @x1 + 2 @x2
@c1 [15;+inf] + 2 @x0 + @x1 + 3 @x2 + @x3
@c2 [-inf;25] + 2 @x1 + 3 @x3
Variables
@x0 [0;+inf]
@x1 [0;10]
@x2 [0;+inf]
@x3 [0;+inf]
";
#[test]
fn test() {
super::callbackmain("psim", super::FileOrText::Text(DFLT_FILE) ).unwrap();
super::callbackmain("dsim", super::FileOrText::Text(DFLT_FILE) ).unwrap();
super::callbackmain("intpnt", super::FileOrText::Text(DFLT_FILE) ).unwrap();
}
}
ceo1.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : ceo1.rs
//!
//! Description:
//! Demonstrates how to solve a small conic exponential
//! optimization problem using the MOSEK API.
//!
//! Min x1 + x2
//! Such that
//! x1 + x2 + x3 = 1.0
//! |1 | |x1|
//! | 1 | x |x2| in K_exp
//! | 1| |x3|
//! x1,x2,x3 are free
//!
extern crate mosek;
use mosek::{Task,Boundkey,Objsense,Streamtype,Solsta,Soltype};
const INF : f64 = 0.0;
fn main() -> Result<(),String> {
let numcon = 1;
let numvar = 3;
let bkc = mosek::Boundkey::FX;
let blc = 1.0;
let buc = 1.0;
let bkx = vec![ Boundkey::FR,
Boundkey::FR,
Boundkey::FR ];
let blx = vec![ -INF, -INF, -INF ];
let bux = vec![ INF, INF, INF ];
let c = vec![ 1.0, 1.0, 0.0 ];
let a = vec![ 1.0, 1.0, 1.0 ];
let asub = vec![0, 1, 2];
//let csub = new int[numvar];
//double[] xx = new double[numvar];
/* Create the optimization task. */
Task::new().expect("Failed to create task")
.with_stream_callback(
Streamtype::LOG,
&mut|msg| print!("{}",msg),
|task| task.with_callback(
&mut |caller| { println!("caller = {}",caller); false },
|task| {
/* Append 'numcon' empty constraints.
The constraints will initially have no bounds. */
task.append_cons(numcon)?;
/* Append 'numvar' variables.
The variables will initially be fixed at zero (x=0). */
task.append_vars(numvar)?;
/* Define the linear part of the problem */
task.put_c_slice(0, numvar, c.as_slice())?;
task.put_a_row(0, asub.as_slice(), a.as_slice())?;
task.put_con_bound(0, bkc, blc, buc)?;
task.put_var_bound_slice(0, numvar, bkx.as_slice(), blx.as_slice(), bux.as_slice())?;
/* Add a conic constraint */
task.append_afes(3)?;
let afeidxs = vec![0, 1, 2 ];
let b = vec![0.0,0.0,0.0];
let domidx = task.append_primal_exp_cone_domain()?;
task.put_afe_f_row_list(afeidxs.as_slice(),
vec![1,1,1].as_slice(),
vec![0,1,2].as_slice(),
vec![0,1,2].as_slice(),
vec![1.0,1.0,1.0].as_slice())?;
task.append_acc(domidx,afeidxs.as_slice(),b.as_slice())?;
task.put_obj_sense(Objsense::MINIMIZE)?;
println!("optimize");
/* Solve the problem */
task.optimize()?;
// Print a summary containing information
// about the solution for debugging purposes
task.solution_summary(Streamtype::MSG)?;
/* Get status information about the solution */
let solsta = task.get_sol_sta(Soltype::ITR)?;
assert!(solsta == Solsta::OPTIMAL);
let mut xx = vec![0.0; numvar as usize];
task.get_xx(Soltype::ITR, & mut xx[..])?;
println!("Optimal primal solution");
for j in 0..numvar as usize {
println!("x[{}]: {:.4}",j,xx[j]);
}
Ok(())
}))
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
concurrent1.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : concurrent1.rs
//!
//! Purpose: Demonstrates a simple implementation of a concurrent optimizer.
//!
//! The concurrent optimizer starts a few parallel optimizations
//! of the same problem using different algorithms, and reports
//! a solution when the first optimizer is ready.
//!
//! This example also demonstrates how to define a simple callback handler
//! that stops the optimizer when requested.
extern crate mosek;
use mosek::{Task,Objsense,Streamtype,Solsta,Soltype};
use std::sync::{Arc,Mutex};
use std::cmp::Ordering;
use std::thread;
use std::env;
fn optimize(mut t : mosek::Task, stop : Arc<Mutex<bool>>) -> Option<(i32,mosek::Task)> {
let cbstop = Arc::clone(&stop);
if let Some(trm) = t.with_callback(
&mut |_| cbstop.lock().and_then(|p| Ok(*p)).unwrap_or(false),
|task|
if let Ok(trm) = task.optimize() {
let mut st = stop.lock().unwrap();
*st = true;
Some(trm)
} else {
None
}) {
Some((trm,t))
}
else {
None
}
}
fn optimize_concurrent(task : &mut mosek::Task,
optimizers : &[i32]) -> Vec<(usize,i32,mosek::Task)> {
let stop = Arc::new(Mutex::new(false));
optimizers.iter().enumerate()
.filter_map(|(i,&ot)|
if let Some(mut t) = task.clone() {
if let Err(_) = t.put_int_param(mosek::Iparam::OPTIMIZER, ot as i32) { None }
else {
let stopopt = Arc::clone(&stop);
Some((i,thread::spawn(move || optimize(t,stopopt))))
}
}
else { None })
.filter_map(|(i,th)|
match th.join().unwrap() {
None => None,
Some((r,t)) => Some((i,r,t)) } )
.collect()
}
fn optimize_concurrent_mio(task : & mut mosek::Task,
seeds : &[i32]) -> Vec<(usize,i32,mosek::Task)> {
let stop = Arc::new(Mutex::new(false));
seeds.iter().enumerate()
.filter_map(|(i,&seed)| {
if let Some(mut t) = task.clone() {
if let Err(_) = t.put_int_param(mosek::Iparam::MIO_SEED, seed) { None }
else {
let stopopt = Arc::clone(&stop);
Some((i,thread::spawn(move || optimize(t,stopopt))))
}
}
else { None }})
.filter_map(|(i,th)|
match th.join().unwrap() {
None => None,
Some((r,t)) => Some((i,r,t)) } )
.collect()
}
enum FileOrText {
File(String),
Text(String)
}
fn main() -> Result<(),String> {
let mut args = env::args();
if args.len() < 2 {
println!("Syntax: concurrent1 FILENAME [ TIMELIMIT ]");
return Err("Invalid argument list".to_string());
}
let _ = args.next();
let filename = args.next().unwrap();
let timelimit = args.next();
concurrent1(FileOrText::File(filename),
timelimit)
}
fn concurrent1(data : FileOrText, timelimit : Option<String>) -> Result<(),String> {
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
};
match data {
FileOrText::File(fname) => task.read_data(fname.as_str())?,
FileOrText::Text(text) => task.read_ptf_string(text.as_str())?
}
if let Some(timelimit) = timelimit {
task.put_dou_param(mosek::Dparam::OPTIMIZER_MAX_TIME, timelimit.parse().unwrap())?;
}
let numintvar = task.get_num_int_var()?;
let r = if numintvar == 0 {
let optimizers = &[mosek::Optimizertype::CONIC,
mosek::Optimizertype::DUAL_SIMPLEX,
mosek::Optimizertype::PRIMAL_SIMPLEX];
optimize_concurrent(& mut task, optimizers)
}
else {
let seeds = &[ 42, 13, 71749373 ];
optimize_concurrent_mio(& mut task, seeds)
};
let sense = task.get_obj_sense()?;
// Pick the feasible result. For non-integer problems all
// solutions should be the same if more than one is returned, but
// for integer problems tasks may have hit the time limit and
// returned non-optimal solutions.
let n = r.len();
if n == 0 {
println!("All optimizers failed.");
}
else if numintvar > 0 {
let (ii,_,tasks) = split3vec(r);
let pobjs : Vec<(usize,f64)> =
ii.iter().zip(tasks.iter()).enumerate()
.filter_map(|(k,(_i,t))|
match (*t).solution_def(Soltype::ITG) {
Ok(true) => match t.get_sol_sta(Soltype::ITG).unwrap() {
Solsta::PRIM_FEAS|Solsta::INTEGER_OPTIMAL => Some((k,t.get_primal_obj(Soltype::ITG).unwrap())),
_ => None
},
_ => None
})
.collect();
let &(besti,bestobj) = pobjs.iter()
.max_by(|(_,o1),(_,o2)|
match sense {
Objsense::MAXIMIZE => if o1 < o2 {Ordering::Less} else if o2 < o1 {Ordering::Greater} else {Ordering::Equal},
_ => if o1 > o2 {Ordering::Less} else if o2 > o1 {Ordering::Greater} else {Ordering::Equal}
}).unwrap();
drop_except(tasks,besti).unwrap()
.with_stream_callback(
Streamtype::LOG,
&mut |msg| print!("{}",msg),
|t| {
t.optimizer_summary(mosek::Streamtype::LOG)?;
t.solution_summary(mosek::Streamtype::LOG)?;
Ok::<(),String>(())
})?;
println!("{} optimizers succeeded:",pobjs.len());
for &(k,v) in pobjs.iter() {
println!("Optimizer with seed #{} produced result : {:.5e}",ii[k],v);
}
println!("\tBest solution is #{}: {:.5e}",ii[besti],bestobj);
}
else {
let (ii,_,tasks) = split3vec(r);
let pobjs : Vec<(usize,f64)> =
tasks.iter().enumerate()
.filter_map(|(k,t)|
match t.get_sol_sta(Soltype::BAS) {
Ok(Solsta::PRIM_FEAS)|Ok(Solsta::OPTIMAL) => Some((k,t.get_primal_obj(Soltype::BAS).unwrap())),
_ => None
}.or_else(|| match t.get_sol_sta(Soltype::ITR) {
Ok(Solsta::PRIM_FEAS)|Ok(Solsta::OPTIMAL) => Some((k,t.get_primal_obj(Soltype::ITR).unwrap())),
_ => None
}))
.collect();
let &(besti,bestobj) = pobjs.iter()
.max_by(|(_,o1),(_,o2)|
match sense {
Objsense::MAXIMIZE => if o1 < o2 {Ordering::Less} else if o2 < o1 {Ordering::Greater} else {Ordering::Equal},
_ => if o1 > o2 {Ordering::Less} else if o2 > o1 {Ordering::Greater} else {Ordering::Equal}
}).unwrap();
drop_except(tasks,besti).unwrap()
.with_stream_callback(
Streamtype::LOG,
&mut|msg| print!("{}",msg),
|t| {
t.optimizer_summary(mosek::Streamtype::LOG)?;
t.solution_summary(mosek::Streamtype::LOG)?;
Ok::<(),String>(())
})?;
println!("{} optimizers succeeded:",pobjs.len());
for &(k,v) in pobjs.iter() {
println!("Optimizer with seed #{} produced result : {:.5e}",ii[k],v);
}
println!("\tBest solution is #{}: {:.5e}",ii[besti],bestobj);
}
Result::Ok(())
}
fn drop_except<A>(mut a : Vec<A>, idx : usize ) -> Option<A> {
if idx >= a.len() { None }
else {
for _i in 0..a.len()-idx-1 { let _ = a.pop(); }
a.pop()
}
}
fn split3vec<A,B,C>(mut v : Vec<(A,B,C)>) -> (Vec<A>,Vec<B>,Vec<C>) {
let mut ra = Vec::with_capacity(v.len());
let mut rb = Vec::with_capacity(v.len());
let mut rc = Vec::with_capacity(v.len());
while let Some((va,vb,vc)) = v.pop() {
ra.push(va);
rb.push(vb);
rc.push(vc);
}
ra.reverse();
rb.reverse();
rc.reverse();
(ra,rb,rc)
}
#[cfg(test)]
mod tests {
const DFLT_FILE1 : &str = "Task
Objective
Maximize + 2 @x0 + 3 @x1 - @x2
Constraints
@c0 [1] + @x0 + @x1 + @x2
@C0 [QUAD(3)]
@ac1: + 0.03
@ac2: + 1.5 @x0 + 0.1 @x1
@ac3: + 0.3 @x0 + 2.1 @x2 + 0.1
Variables
@x0
@x1
@x2
";
const DFLT_FILE2 : &str = "Task
Objective
Maximize + @x0 + 0.64 @x1
Constraints
@c0 [-inf;250] + 50 @x0 + 31 @x1
@c1 [-4;+inf] + 3 @x0 - 2 @x1
Variables
@x0[0;+inf]
@x1[0;+inf]
Integers
@x0 @x1
";
#[test]
fn test() {
super::concurrent1(super::FileOrText::Text(DFLT_FILE1.to_string()),None).unwrap();
super::concurrent1(super::FileOrText::Text(DFLT_FILE2.to_string()),Some("100.0".to_string())).unwrap();
}
}
cqo1.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : cqo1.rs
//!
//! Purpose: To demonstrate how to solve a small conic quadratic
//! optimization problem using the MOSEK API.
//!
extern crate mosek;
extern crate itertools;
const INF : f64 = 0.0;
use mosek::{Task,Streamtype,Solsta,Soltype};
use itertools::izip;
fn main() -> Result<(),String>
{
let numvar : i32 = 6;
let numcon : i32 = 1;
let bkc = &[ mosek::Boundkey::FX ];
let blc = &[ 1.0 ];
let buc = &[ 1.0 ];
let bkx = &[ mosek::Boundkey::LO,
mosek::Boundkey::LO,
mosek::Boundkey::LO,
mosek::Boundkey::FR,
mosek::Boundkey::FR,
mosek::Boundkey::FR];
let blx = &[ 0.0, 0.0, 0.0, -INF, -INF, -INF ];
let bux = &[ INF, INF, INF, INF, INF, INF ];
let c = &[ 0.0, 0.0, 0.0, 1.0, 1.0, 1.0 ];
let asub = &[0, 1, 2];
let aval = &[1.0, 1.0, 1.0];
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* Append 'numcon' empty constraints.
* The constraints will initially have no bounds. */
task.append_cons(numcon)?;
/* Append 'numvar' variables.
* The variables will initially be fixed at zero (x=0). */
task.append_vars(numvar)?;
for (j,&cj,&bkj,&blj,&buj) in izip!(0..numvar,c,bkx,blx,bux) {
/* Set the linear term c_j in the objective.*/
task.put_c_j(j,cj)?;
/* Set the bounds on variable j.
* blx[j] <= x_j <= bux[j] */
task.put_var_bound( j, /* Index of variable.*/
bkj, /* Bound key.*/
blj, /* Numerical value of lower bound.*/
buj)?; /* Numerical value of upper bound.*/
}
/* Input columns of A */
task.put_a_row(0, asub, aval)?;
/* Set the bounds on constraints.
* for i=1, ...,numcon : blc[i] <= constraint i <= buc[i] */
for (i,&bki,&bli,&bui) in izip!(0..numcon,bkc,blc,buc) {
task.put_con_bound( i, /* Index of constraint.*/
bki, /* Bound key.*/
bli, /* Numerical value of lower bound.*/
bui)?; /* Numerical value of upper bound.*/
}
/* Append the first cone. */
// Create a matrix F such that F * x = [x(3),x(0),x(1),x(4),x(5),x(2)]
{
task.append_afes(6)?;
task.put_afe_f_entry_list(&[0,1,2,3,4,5], // Rows
&[3, 0, 1, 4, 5, 2], // Columns
&[1.0,1.0,1.0,1.0,1.0,1.0])?;
// Quadratic cone (x(3),x(0),x(1)) \in QUAD_3
let qconedom = task.append_quadratic_cone_domain(3)?;
task.append_acc(qconedom, // Domain
&[0, 1, 2], // Rows from F
&[0.0,0.0,0.0])?; // Unused
// Rotated quadratic cone (x(4),x(5),x(2)) \in RQUAD_3
let rqconedom = task.append_r_quadratic_cone_domain(3)?;
task.append_acc(rqconedom, // Domain
&[3, 4, 5], // Rows from F
&[0.0,0.0,0.0])?; // Unused
}
/* Run optimizer */
let trm = task.optimize()?;
task.write_data("cqo1.ptf")?;
/* Print a summary containing information
* about the solution for debugging purposes*/
task.solution_summary (Streamtype::MSG)?;
let solsta = task.get_sol_sta(Soltype::ITR)?;
match solsta
{
Solsta::OPTIMAL =>
{
let mut xx = vec![0.0; numvar as usize];
task.get_xx(Soltype::ITR, /* Request the basic solution. */
xx.as_mut_slice())?;
println!("Optimal primal solution");
for (j,xj) in izip!(0..numvar,xx) {
println!("x[{}]: {}",j,xj);
}
}
Solsta::DUAL_INFEAS_CER |
Solsta::PRIM_INFEAS_CER =>
{
println!("Primal or dual infeasibility certificate found.");
}
Solsta::UNKNOWN => {
/* If the solutions status is unknown, print the termination code
* indicating why the optimizer terminated prematurely. */
println!("The solution status is unknown.");
println!("The optimizer terminitated with code: {}",trm);
}
_ =>
{
println!("Other solution status.");
}
}
Ok(())
} /* main */
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
djc1.rs
//
// Copyright: MOSEK ApS
//
// File: djc1.rs
//
// 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
//
extern crate mosek;
use mosek::{Task,Boundkey,Objsense,Streamtype,Solsta,Soltype};
// Since the value of infinity is ignored, we define it solely
// for symbolic purposes
const INF : f64 = 0.0;
fn main() -> Result<(),String> {
// Create a task object
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
// Append free variables
let numvar : i32 = 4;
task.append_vars(numvar)?;
let x : Vec<i32> = (0..numvar).collect();
task.put_var_bound_slice_const(0, numvar, Boundkey::FR, -INF, INF)?;
// The linear part: the linear constraint
task.append_cons(1)?;
task.put_a_row(0, x.as_slice(), vec![1.0; numvar as usize].as_slice())?;
task.put_con_bound(0, Boundkey::LO, -10.0, -10.0)?;
// The linear part: objective
task.put_obj_sense(Objsense::MINIMIZE)?;
task.put_c_list(x.as_slice(), &[2.0, 1.0, 3.0, 1.0])?;
// Fill in the affine expression storage F, g
let numafe : i64 = 10;
task.append_afes(numafe)?;
let fafeidx : &[i64] = &[0, 0, 1, 1, 2, 3, 4, 5, 6, 7, 8, 9];
let fvaridx : &[i32] = &[0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3];
let fval = &[1.0, -2.0, 1.0, -3.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
let g = &[1.0, 2.0, 0.0, 0.0, 0.0, 0.0, -2.5, -2.5, -2.5, -2.5];
task.put_afe_f_entry_list(fafeidx, fvaridx, fval)?;
task.put_afe_g_slice(0, numafe, g)?;
// Create domains
let zero1 = task.append_rzero_domain(1)?;
let zero2 = task.append_rzero_domain(2)?;
let rminus1 = task.append_rminus_domain(1)?;
// Append disjunctive constraints
let numdjc : i64 = 2;
task.append_djcs(numdjc)?;
// First disjunctive constraint
task.put_djc(0, // DJC index
&[rminus1, zero2, rminus1, zero2], // Domains (domidxlist)
&[0, 4, 5, 1, 2, 3], // AFE indices (afeidxlist)
&[0.0,0.0,0.0,0.0,0.0,0.0], // Unused
&[2, 2])?; // Term sizes (termsizelist)
// Second disjunctive constraint
task.put_djc(1, // DJC index
&[zero1, zero1, zero1, zero1], // Domains (domidxlist)
&[6, 7, 8, 9], // AFE indices (afeidxlist)
&[0.0,0.0,0.0,0.0], // Unused
&[1, 1, 1, 1])?; // Term sizes (termidxlist)
// Useful for debugging
task.write_data("djc1.ptf")?; // Write file in human-readable format
// Solve the problem
let _ = task.optimize()?;
// Print a summary containing information
// about the solution for debugging purposes
task.solution_summary(Streamtype::MSG)?;
// Get status information about the solution
let sta = task.get_sol_sta(Soltype::ITG)?;
let mut xx = vec![0.0; numvar as usize];
task.get_xx(Soltype::ITG,xx.as_mut_slice())?;
println!("Optimal solution: ");
for (i,&xi) in xx.iter().enumerate() {
println!("x[{}]={}",i,xi);
}
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
feasrepairex1.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : feasrepairex1.rs
//!
//! Purpose : To demonstrate how to use the MSK_relaxprimal function to
//! locate the cause of an infeasibility.
//!
//! Syntax : On command line
//!
//! feasrepairex1 feasrepair.lp
//!
//! feasrepair.lp is located in mosek/<version>/tools/examples/data
//!*/
extern crate mosek;
use std::env;
use mosek::{Task,Streamtype};
enum FileOrText {
File(String),
Text(String)
}
fn main() -> Result<(),String> {
let mut args = env::args();
if args.len() < 2 {
println!("Syntax: feasrepairex1 FILENAME");
return Err("Invalid argument list".to_string());
}
let _ = args.next();
feasrepairex1(FileOrText::File(args.next().unwrap()))
}
fn feasrepairex1(filename : FileOrText) -> Result<(),String> {
let mut task = Task::new().unwrap().with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
match filename {
FileOrText::File(fname) => task.read_data(fname.as_str())?,
FileOrText::Text(data) => task.read_lp_string(data.as_str())?
}
task.put_int_param(mosek::Iparam::LOG_FEAS_REPAIR, 3)?;
let wc = vec![1.0; task.get_num_con()? as usize];
let wx = vec![1.0; task.get_num_var()? as usize];
task.primal_repair(wc.as_slice(),wc.as_slice(),wx.as_slice(),wx.as_slice())?;
let sum_viol = task.get_dou_inf(mosek::Dinfitem::PRIMAL_REPAIR_PENALTY_OBJ)?;
println!("Minimized sum of violations = {}", sum_viol);
let _ = task.optimize()?;
task.solution_summary(mosek::Streamtype::MSG)?;
Ok(())
}
#[cfg(test)]
mod tests {
/// Small infeasible linear problem
const FEASREPAIR_LP : &str = "
minimize
obj: - 10 x1 - 9 x2
st
c1: + 7e-01 x1 + x2 <= 630
c2: + 5e-01 x1 + 8.333333333e-01 x2 <= 600
c3: + x1 + 6.6666667e-01 x2 <= 708
c4: + 1e-01 x1 + 2.5e-01 x2 <= 135
bounds
x2 >= 650
end
";
#[test]
fn test() {
super::feasrepairex1(super::FileOrText::Text(FEASREPAIR_LP.to_string())).unwrap();
}
}
gp1.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : gp1.rs
//!
//! 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
//!
extern crate mosek;
use mosek::{Task,Boundkey,Objsense,Streamtype,Soltype};
// Since the value of infinity is ignored, we define it solely
// for symbolic purposes
const INF : f64 = 0.0;
#[allow(non_snake_case)]
fn max_volume_box(Aw : f64,
Af : f64,
alpha : f64,
beta : f64,
gamma : f64,
delta : f64) -> Result<Vec<f64>,String>
{
let numvar = 3i32; // Variables in original problem
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
// Directs the log task stream to the user specified
// method task_msg_obj.stream
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
// Add variables and constraints
task.append_vars(numvar)?;
let x = 0i32;
let y = 1i32;
let z = 2i32;
// Objective is the sum of three first variables
task.put_obj_sense(Objsense::MAXIMIZE)?;
task.put_c_slice(0, numvar, &[1.0,1.0,1.0])?;
task.put_var_bound_slice_const(0, numvar, Boundkey::FR, -INF, INF)?;
task.append_cons(3)?;
// s0+s1 < 1 <=> log(s0+s1) < 0
task.put_aij_list(&[0,0,1,1,2,2],
&[y, z, x, y, z, y],
&[1.0, 1.0, 1.0, -1.0, 1.0, -1.0])?;
task.put_con_bound(0,Boundkey::UP,-INF,Af.ln())?;
task.put_con_bound(1,Boundkey::RA,alpha.ln(),beta.ln())?;
task.put_con_bound(2,Boundkey::RA,gamma.ln(),delta.ln())?;
{
let afei = task.get_num_afe()?;
let u1 = task.get_num_var()?;
let u2 = u1+1;
let afeidx = &[0, 1, 2, 2, 3, 3, 5, 5];
let varidx = &[u1, u2, x, y, x, z, u1, u2];
let fval = &[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
let gfull = &[0.0, 0.0, (2.0/Aw).ln(), (2.0/Aw).ln(), 1.0, -1.0];
task.append_vars(2)?;
task.append_afes(6)?;
task.put_var_bound_slice_const(u1, u1+2, Boundkey::FR, -INF, INF)?;
// Affine expressions appearing in affine conic constraints
// in this order:
// u1, u2, x+y+log(2/Awall), x+z+log(2/Awall), 1.0, u1+u2-1.0
task.put_afe_f_entry_list(afeidx, varidx, fval)?;
task.put_afe_g_slice(afei, afei+6, gfull)?;
{
let dom = task.append_primal_exp_cone_domain()?;
// (u1, 1, x+y+log(2/Awall)) \in EXP
task.append_acc(dom, &[0, 4, 2], &[0.0,0.0,0.0])?;
// (u2, 1, x+z+log(2/Awall)) \in EXP
task.append_acc(dom, &[1, 4, 3], &[0.0,0.0,0.0])?;
}
{
let dom = task.append_rzero_domain(1)?;
// The constraint u1+u2-1 \in \ZERO is added also as an ACC
task.append_acc(dom, &[5], &[0.0])?;
}
}
let _trm = task.optimize()?;
task.write_data("gp1.ptf")?;
let mut xyz = vec![0.0; 3];
task.get_xx_slice(Soltype::ITR, 0i32, numvar, xyz.as_mut_slice())?;
// task.write_data("gp1.ptf")?;
Ok(xyz.iter().map(|v| v.exp()).collect())
}
#[allow(non_snake_case)]
fn main() -> Result<(),String> {
// 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 )
//
// Finally, the model we will implement:
//
// maximize x+y+z
// subject to s0 > exp(x+y+log(2/Awall); (s0,1,x+y+log(2/Awall)) in PEXP
// s1 > exp(x+z+log(2/Awall); (s1,1,x+z+log(2/Awall)) in PEXP
// s0+s1 < 1
//
// y+z < log Afloor
//
// x-y in [log alpha; log beta]
// z-y in [log gamma; log delta]
//
// (x,y,z) in pexp : x0 > x1 * exp(x2/x1)
let Aw = 200.0;
let Af = 50.0;
let alpha = 2.0;
let beta = 10.0;
let gamma = 2.0;
let delta = 10.0;
let hwd = max_volume_box(Aw, Af, alpha, beta, gamma, delta)?;
println!("h={:.4} w={:.4} d={:.4}\n", hwd[0], hwd[1], hwd[2]);
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
helloworld.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : helloworld.rs
//!
//! The most basic example of how to get started with MOSEK.
//!
extern crate mosek;
use mosek::{Task,Boundkey,Objsense,Soltype};
fn main() -> Result<(),String> {
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
};
task.append_vars(1)?; // 1 variable x
task.put_c_j(0, 1.0)?; // c_0 = 1.0
task.put_var_bound(0, Boundkey::RA, 2.0, 3.0)?; // 2.0 <= x <= 3.0
task.put_obj_sense(Objsense::MINIMIZE)?; // minimize
task.optimize()?; // Optimize
let mut x = vec![0.0; 1];
task.get_xx(Soltype::ITR, x.as_mut_slice())?; // Get solution
println!("Solution x = {}", x[0]); // Print solution
return Result::Ok(());
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
lo1.rs
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : lo1.rs
//!
//! Purpose: To demonstrate how to solve a small linear
//! optimization problem using the MOSEK C API,
//! and handle the solver result and the problem
//! solution.
//!
//! ```
//! Maximize 3 x1 + x2 + 5 x3 + x4
//! Such that
//! 3 x1 + x2 + 2 x3 = 39
//! 2 x1 + x2 + 3 x3 + x4 >= 15
//! + 2 x2 + 3 x4 <= 25
//!
//! x1,x3,x4 >= 0
//! 0 <= x2 <= 10
//! ```
extern crate mosek;
use mosek::{Task,Boundkey,Objsense,Streamtype,Solsta,Soltype};
const INF : f64 = 0.0;
fn main() -> Result<(),String> {
let numvar = 4;
let numcon = 3;
let c = vec![3.0, 1.0, 5.0, 1.0];
/* Below is the sparse representation of the A
* matrix stored by column. */
let aptrb = vec![ 0, 2, 5, 7 ];
let aptre = vec![ 2, 5, 7, 9 ];
let asub = vec![ 0, 1,
0, 1, 2,
0, 1,
1, 2 ];
let aval = vec![ 3.0, 2.0,
1.0, 1.0, 2.0,
2.0, 3.0,
1.0, 3.0 ];
/* Bounds on constraints. */
let bkc = vec![ Boundkey::FX, Boundkey::LO, Boundkey::UP ];
let blc = vec![ 30.0, 15.0, -INF ];
let buc = vec![ 30.0, INF, 25.0 ];
/* Bounds on variables. */
let bkx = vec![ Boundkey::LO, Boundkey::RA, Boundkey::LO, Boundkey::LO ];
let blx = vec![ 0.0, 0.0, 0.0, 0.0 ];
let bux = vec![ INF, 10.0, INF, INF ];
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
/* Directs the log task stream to the 'printstr' function. */
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* Append 'numcon' empty constraints.
* The constraints will initially have no bounds. */
task.append_cons(numcon as i32)?;
/* Append 'numvar' variables.
* The variables will initially be fixed at zero (x=0). */
task.append_vars(numvar as i32)?;
for j in 0..numvar
{
/* Set the linear term c_j in the objective.*/
task.put_c_j(j as i32,c[j])?;
/* Set the bounds on variable j.
* blx[j] <= x_j <= bux[j] */
task.put_var_bound(j as i32, /* Index of variable.*/
bkx[j], /* Bound key.*/
blx[j], /* Numerical value of lower bound.*/
bux[j])?; /* Numerical value of upper bound.*/
/* Input column j of A */
task.put_a_col(j as i32, /* Variable (column) index.*/
& asub[aptrb[j]..aptre[j]], /* Pointer to row indexes of column j.*/
& aval[aptrb[j]..aptre[j]])?; /* Pointer to Values of column j.*/
}
/* Set the bounds on constraints.
* for i=1, ...,numcon : blc[i] <= constraint i <= buc[i] */
for i in 0..numcon {
task.put_con_bound(i as i32, /* Index of constraint.*/
bkc[i], /* Bound key.*/
blc[i], /* Numerical value of lower bound.*/
buc[i])?; /* Numerical value of upper bound.*/
}
/* Maximize objective function. */
task.put_obj_sense(Objsense::MAXIMIZE)?;
/* Run optimizer */
let _trmcode = task.optimize()?;
/* Print a summary containing information
* about the solution for debugging purposes. */
task.solution_summary(Streamtype::LOG)?;
let solsta = task.get_sol_sta(Soltype::BAS)?;
match solsta
{
Solsta::OPTIMAL =>
{
let mut xx = vec![0.0,0.0,0.0,0.0];
task.get_xx(Soltype::BAS, /* Request the basic solution. */
& mut xx[..])?;
println!("Optimal primal solution");
for j in 0..numvar as usize
{
println!("x[{}]: {}",j,xx[j]);
}
}
Solsta::DUAL_INFEAS_CER |
Solsta::PRIM_INFEAS_CER =>
{
println!("Primal or dual infeasibility certificate found.");
}
Solsta::UNKNOWN =>
{
/* If the solutions status is unknown, print the termination code
* indicating why the optimizer terminated prematurely. */
println!("The solution status is unknown.");
println!("The optimizer terminitated with code: {}",solsta);
}
_ =>
{
println!("Other solution status.");
}
}
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
lo2.rs
// Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File : lo2.rs
//
// Purpose : Demonstrates how to solve a small linear
// optimization problem using the MOSEK Java API.
extern crate mosek;
extern crate itertools;
use mosek::{Task,Boundkey,Objsense,Streamtype,Solsta,Soltype};
use itertools::izip;
const INF : f64 = 0.0;
fn main() -> Result<(),String> {
let numcon : i32 = 3;
let numvar : i32 = 4;
let c = [3.0, 1.0, 5.0, 1.0];
let asub = [0, 1, 2,
0, 1, 2, 3,
1, 3];
let aval = [3.0, 1.0, 2.0,
2.0, 1.0, 3.0, 1.0,
2.0, 3.0];
let aptr = [0,3,7,9];
let bkc = [Boundkey::FX,
Boundkey::LO,
Boundkey::UP];
let blc = [30.0,
15.0,
-INF];
let buc = [30.0,
INF,
25.0];
let bkx = [Boundkey::LO,
Boundkey::RA,
Boundkey::LO,
Boundkey::LO];
let blx = [0.0,
0.0,
0.0,
0.0];
let bux = [INF,
10.0,
INF,
INF];
// Create a task object linked with the environment env.
let mut task = Task::new().unwrap().with_callbacks();
// Directs the log task stream to the user specified
// method task_msg_obj.streamCB
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* Give MOSEK an estimate of the size of the input data.
This is done to increase the speed of inputting data.
However, it is optional. */
/* Append 'numcon' empty constraints.
The constraints will initially have no bounds. */
task.append_cons(numcon)?;
/* Append 'numvar' variables.
The variables will initially be fixed at zero (x=0). */
task.append_vars(numvar)?;
/* Optionally add a constant term to the objective. */
task.put_cfix(0.0)?;
for (j,(&cj,&bkj,&blj,&buj)) in izip!(c.iter(),bkx.iter(),blx.iter(),bux.iter()).enumerate() {
/* Set the linear term c_j in the objective.*/
task.put_c_j(j as i32, cj)?;
/* Set the bounds on variable j.
blx[j] <= x_j <= bux[j] */
task.put_var_bound(j as i32, bkj, blj, buj)?;
};
/* Set the bounds on constraints.
for i=1, ...,numcon : blc[i] <= constraint i <= buc[i] */
for (i,(&bki,&bli,&bui,&ptrb,&ptre)) in izip!(bkc.iter(),blc.iter(),buc.iter(),aptr[..aptr.len()-1].iter(),aptr[1..].iter()).enumerate() {
task.put_con_bound(i as i32, bki, bli, bui)?;
/* Input row i of A */
task.put_a_row(i as i32, /* Row index.*/
&asub[ptrb..ptre], /* Column indexes of non-zeros in row i.*/
&aval[ptrb..ptre])?; /* Non-zero Values of row i. */
}
task.put_obj_sense(Objsense::MAXIMIZE)?;
task.optimize()?;
// Print a summary containing information
// about the solution for debugging purposes
task.solution_summary(Streamtype::MSG)?;
/* Get status information about the solution */
let solsta = task.get_sol_sta(Soltype::BAS)?;
let mut xx = vec![0.0; numvar as usize];
task.get_xx(Soltype::BAS, // Basic solution.
xx.as_mut_slice());
match solsta {
Solsta::OPTIMAL =>
println!("Optimal primal solution = {:?}",xx),
Solsta::DUAL_INFEAS_CER|Solsta::PRIM_INFEAS_CER =>
println!("Primal or dual infeasibility."),
Solsta::UNKNOWN =>
println!("Unknown solution status."),
_ =>
println!("Other solution status")
}
Ok(())
}
logistic.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : logistic.rs
//!
//! Purpose: Implements logistic regression with regulatization.
//!
//! Demonstrates using the exponential cone and log-sum-exp in Optimizer API.
//!
//! Plots an example for 2D datasets.
extern crate mosek;
extern crate itertools;
use mosek::{Task,TaskCB,Boundkey,Objsense,Soltype,Streamtype};
use itertools::{izip,iproduct};
const INF : f64 = 0.0;
#[allow(non_snake_case)]
fn softplus(task : & mut TaskCB, d : i32, n : i32, theta : i32, t : i32, X : &[f64], Y : &[bool]) -> Result<(),String> {
let nvar = task.get_num_var()?;
let ncon = task.get_num_con()?;
let nafe = task.get_num_afe()?;
task.append_vars(2*n)?; // z1, z2
task.append_cons(n)?; // z1 + z2 = 1
task.append_afes(4*n as i64)?; //theta * X[i] - t[i], -t[i], z1[i], z2[i]
let z1 = nvar;
let z2 = nvar+n;
let zcon = ncon;
let thetaafe = nafe;
let tafe = nafe+n as i64;
let z1afe = tafe+n as i64;
let z2afe = z1afe+n as i64;
// Linear constraints
{
let mut subi = vec![0i32; 2*n as usize];
let mut subj = vec![0i32; 2*n as usize];
let mut aval = vec![0.0; 2*n as usize];
for i in 0..n {
task.put_var_name(z1+i,format!("z1[{}]",i).as_str())?;
task.put_var_name(z2+i,format!("z2[{}]",i).as_str())?;
}
for ((i,&zx),si, sj, av) in izip!(iproduct!(0..n,&[z1,z2]),
subi.iter_mut(),
subj.iter_mut(),
aval.iter_mut()) {
*si = zcon+i;
*sj = zx+i as i32;
*av = 1.0;
}
task.put_aij_list(subi.as_slice(), subj.as_slice(), aval.as_slice())?;
task.put_con_bound_slice_const(zcon, zcon+n, Boundkey::FX, 1.0, 1.0)?;
task.put_var_bound_slice_const(nvar, nvar+2*n, Boundkey::FR, -INF, INF)?;
}
// Affine conic expressions
let mut afeidx = vec![0i64; (d*n+4*n) as usize];
let mut varidx = vec![0i32; (d*n+4*n) as usize];
let mut fval = vec![0.0; (d*n+4*n) as usize];
// Thetas
let mut k : usize = 0;
for ((i,j),afei,vari) in izip!(iproduct!(0..n,0..d),
& mut afeidx[k..k+(n*d) as usize],
& mut varidx[k..k+(n*d) as usize]) {
*afei = thetaafe + i as i64;
*vari = theta + j;
}
for ((&yi,_j),&xij,fv) in izip!(iproduct!(Y,0..d), X, & mut fval[k..k+(n*d) as usize]) {
*fv = (if yi {-1.0} else {1.0}) * xij;
}
k += (n*d) as usize;
for fv in fval[k..k+(2*n) as usize].iter_mut() { *fv = -1.0; }
for (i,afei,vari) in izip!(0..n,
&mut afeidx[k..k+(2*n) as usize].iter_mut().step_by(2),
&mut varidx[k..k+(2*n) as usize].iter_mut().step_by(2)) {
*afei = thetaafe+i as i64; *vari = t+i;
}
for (i,afei,vari) in izip!(0..n,
&mut afeidx[k+1..k+(2*n) as usize].iter_mut().step_by(2),
&mut varidx[k+1..k+(2*n) as usize].iter_mut().step_by(2)) {
*afei = tafe+i as i64; *vari = t+i;
}
k += (n*2) as usize;
for fv in fval[k..k+(2*n) as usize].iter_mut() { *fv = 1.0; }
for (i,afei,vari) in izip!(0..n,
&mut afeidx[k..k+(2*n) as usize].iter_mut().step_by(2),
&mut varidx[k..k+(2*n) as usize].iter_mut().step_by(2)) {
*afei = z1afe+i as i64; *vari = z1+i;
}
for (i,afei,vari) in izip!(0..n,
&mut afeidx[k+1..k+(2*n) as usize].iter_mut().step_by(2),
&mut varidx[k+1..k+(2*n) as usize].iter_mut().step_by(2)) {
*afei = z2afe+i as i64; *vari = z2+i;
}
// Add the expressions
task.put_afe_f_entry_list(afeidx.as_slice(), varidx.as_slice(), fval.as_slice())?;
{
// Add a single row with the constant expression "1.0"
let oneafe = task.get_num_afe()?;
task.append_afes(1)?;
task.put_afe_g(oneafe, 1.0)?;
// Add an exponential cone domain
let dom = task.append_primal_exp_cone_domain()?;
// Conic constraints
let zeros = &[0.0,0.0,0.0];
let mut acci = task.get_num_acc()?;
for i in 0..n as i64 {
task.append_acc(dom, &[z1afe+i, oneafe, thetaafe+i], zeros)?;
task.append_acc(dom, &[z2afe+i, oneafe, tafe+i], zeros)?;
task.put_acc_name(acci,format!("z1:theta[{}]",i).as_str())?;
task.put_acc_name(acci+1,format!("z2:t[{}]",i).as_str())?;
acci += 2;
}
}
Ok(())
}
// 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
#[allow(non_snake_case)]
fn logistic_regression(X : &[f64],
Y : &[bool],
lamb : f64) -> Result<Vec<f64>,String> {
let n = Y.len() as i32;
let d = (X.len()/Y.len()) as i32; // num samples, dimension
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
// Variables [r; theta; t]
let nvar : i32 = 1+d+n;
task.append_vars(nvar)?;
task.put_var_bound_slice_const(0, nvar, Boundkey::FR, -INF, INF)?;
let (r,theta,t) = (0i32,1i32,1+d);
task.put_var_name(r,"r")?;
for j in 0..d { task.put_var_name(theta+j,format!("theta[{}]",j).as_str())?; }
for j in 0..n { task.put_var_name(t+j,format!("t[{}]",j).as_str())?; }
// Objective lambda*r + sum(t)
task.put_c_j(r, lamb)?;
for i in 0..n { task.put_c_j(t+i, 1.0)?; }
task.put_obj_sense(Objsense::MINIMIZE)?;
// Softplus function constraints
softplus(& mut task, d, n, theta, t, X, Y)?;
// Regularization
// Append a sequence of linear expressions (r, theta) to F
let numafe = task.get_num_afe()?;
task.append_afes(1+d as i64)?;
task.put_afe_f_entry(numafe, r, 1.0)?;
for i in 0..d {
task.put_afe_f_entry(numafe + i as i64 + 1, theta + i, 1.0)?;
}
// Add the constraint
{
let dom = task.append_quadratic_cone_domain((1+d) as i64)?;
task.append_acc_seq(dom,
numafe,
vec![0.0; 1+d as usize].as_slice())?;
}
// Solution
task.write_data("logistic.ptf")?;
task.optimize()?;
let mut xx = vec![0.0; d as usize];
task.get_xx_slice(Soltype::ITR, theta, theta+d as i32,xx.as_mut_slice())?;
Ok(xx)
}
#[allow(non_snake_case)]
fn main() -> Result<(),String> {
// Test: detect and approximate a circle using degree 2 polynomials
let n : i32 = 30;
let d : i32 = 2;
let (X,Y) = {
let mut Y = vec![false;(n*n) as usize];
let mut x1 = vec![0.0;(n*n) as usize];
let mut x2 = vec![0.0;(n*n) as usize];
let step = 2.0/((n-1) as f64);
for (x1i,x2i,yi,(z1,z2)) in izip!(x1.iter_mut(),x2.iter_mut(),Y.iter_mut(),
iproduct!((0..30).map(|i| -1.0 + (i as f64) * step),
(0..30).map(|i| -1.0 + (i as f64) * step))) {
*x1i = z1; *x2i = z2; *yi = z1*z1+z2*z2 >= 0.69;
}
let X : Vec<f64> = iproduct!(izip!(x1,x2),iproduct!(0..d+1,0..d+1)).filter_map(|((p0,p1),(a,b))| if a+b <= d {Some(f64::powi(p0,a)*f64::powi(p1,b))} else {None}).collect();
(X,Y)
};
let theta = logistic_regression( X.as_slice(), Y.as_slice(), 0.1)?;
println!("theta = {:?}",theta);
return Result::Ok(());
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
mico1.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : mico1.rs
//!
//! 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
//!
extern crate mosek;
use mosek::{Task,Boundkey,Objsense,Streamtype,Soltype,Variabletype};
fn main() -> Result<(),String> {
/* Create the optimization task. */
let mut task = match Task::new() {
Some(t) => t,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
let infinity = 0.0; // for symbolic use, value is irrelevant
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
task.append_vars(6)?;
task.append_cons(3)?;
task.put_var_bound_slice_const(0, 6, Boundkey::FR, -infinity, infinity)?;
// Integrality constraints
task.put_var_type_list(vec![1i32,2i32].as_slice(),
vec![Variabletype::TYPE_INT, Variabletype::TYPE_INT].as_slice())?;
// Set up the three auxiliary linear constraints
task.put_aij_list(vec![0i32,0i32,1i32,2i32,2i32].as_slice(),
vec![1i32,3i32,4i32,2i32,5i32].as_slice(),
vec![-1.0,1.0,1.0,1.0,-1.0].as_slice())?;
task.put_con_bound_slice(0, 3,
vec![Boundkey::FX, Boundkey::FX, Boundkey::FX].as_slice(),
vec![-3.8, 1.0, 0.0].as_slice(),
vec![-3.8, 1.0, 0.0].as_slice())?;
// Objective
task.put_obj_sense(Objsense::MINIMIZE)?;
task.put_c_j(0, 1.0)?;
// Conic part of the problem
task.append_afes(6)?;
for i in 0..6 {
task.put_afe_f_entry(i as i64, i as i32, 1.0)?;
}
{
let domidx = task.append_quadratic_cone_domain(3)?;
task.append_acc(domidx,
vec![0i64,1i64,2i64].as_slice(),
vec![0.0,0.0,0.0].as_slice())?;
}
{
let domidx = task.append_primal_exp_cone_domain()?;
task.append_acc(domidx,vec![3i64,4i64,5i64].as_slice(),vec![0.0,0.0,0.0].as_slice())?;
}
// Optimize the task
let _trm = task.optimize()?;
task.solution_summary(Streamtype::MSG)?;
let mut xx = vec![0.0; 2];
task.get_xx_slice(Soltype::ITG, 1, 3, xx.as_mut_slice())?;
println!("x = {} y = {}",xx[0],xx[1]);
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
milo1.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : milo1.rs
//!
//! Purpose : Demonstrates how to solve a small mixed
//! integer linear optimization problem using the MOSEK Java API.
//!
extern crate mosek;
use mosek::{Task,Boundkey,Objsense,Streamtype,Solsta,Prosta,Soltype,Variabletype,Dparam};
fn main() -> Result<(),String> {
let numcon : i32 = 2;
let numvar : i32 = 2;
let infinity = 0.0; // only for symbolic purposes, value never used
let bkc = vec![Boundkey::UP, Boundkey::LO];
let blc = vec![ -infinity, -4.0 ];
let buc = vec![ 250.0, infinity ];
let bkx = vec![ Boundkey::LO, Boundkey::LO ];
let blx = vec![ 0.0, 0.0 ];
let bux = vec![ infinity, infinity ];
let c = vec![1.0, 0.64];
let asub = vec![0, 1,
0, 1];
let aval = vec![50.0, 3.0, 31.0, -2.0];
let ptrb : Vec<usize> = vec![ 0, 2 ];
let ptre : Vec<usize> = vec![ 2, 4 ];
/* Create the optimization task. */
Task::new().expect("Failed to create task")
.with_stream_callback(
Streamtype::LOG,
&mut |msg| print!("{}",msg),
|task| task.with_itg_sol_callback(
&mut |xx| { println!("Found a new solution = {:?}",xx); false },
|task| {
/* Append 'numcon' empty constraints.
The constraints will initially have no bounds. */
task.append_cons(numcon)?;
/* Append 'numvar' variables.
The variables will initially be fixed at zero (x=0). */
task.append_vars(numvar)?;
for ((((j,cj),bk),bl),bu) in (0..numvar).zip(c.iter()).zip(bkx.iter()).zip(blx.iter()).zip(bux.iter()) {
/* Set the linear term c_j in the objective.*/
task.put_c_j(j, *cj)?;
/* Set the bounds on variable j.
blx[j] <= x_j <= bux[j] */
task.put_var_bound(j, *bk, *bl, *bu)?;
/* Input column j of A */
task.put_a_col(j, /* Variable (column) index.*/
&asub[ptrb[j as usize]..ptre[j as usize]], /* Row index of non-zeros in column j.*/
&aval[ptrb[j as usize]..ptre[j as usize]])?; /* Non-zero Values of column j. */
}
// Set the bounds on constraints.
// for i=1, ...,numcon : blc[i] <= constraint i <= buc[i]
for (((i,bk),bl),bu) in (0..numcon).zip(bkc.iter()).zip(blc.iter()).zip(buc.iter()) {
task.put_con_bound(i, *bk, *bl, *bu)?;
}
/* Specify integer variables. */
for j in 0..numvar {
task.put_var_type(j, Variabletype::TYPE_INT)?;
}
/* Set max solution time */
task.put_dou_param(Dparam::MIO_MAX_TIME, 60.0)?;
/* A maximization problem */
task.put_obj_sense(Objsense::MAXIMIZE)?;
/* Solve the problem */
let _trm = task.optimize()?;
// Print a summary containing information
// about the solution for debugging purposes
task.solution_summary(Streamtype::MSG)?;
let mut xx = vec![0.0; numvar as usize];
task.get_xx(Soltype::ITG, xx.as_mut_slice())?;
/* Get status information about the solution */
match task.get_sol_sta(Soltype::ITG)? {
Solsta::INTEGER_OPTIMAL => {
println!("Optimal solution");
for (j,xj) in (0..numvar).zip(xx.iter()) {
println!("x[{}]: {}", j,xj);
}
}
Solsta::PRIM_FEAS => {
println!("Feasible solution");
for (j,xj) in (0..numvar).zip(xx.iter()) {
println!("x[{}]: {}", j,xj);
}
}
Solsta::UNKNOWN => {
match task.get_pro_sta(Soltype::ITG)? {
Prosta::PRIM_INFEAS_OR_UNBOUNDED => {
println!("Problem status Infeasible or unbounded");
}
Prosta::PRIM_INFEAS => {
println!("Problem status Infeasible.");
}
Prosta::UNKNOWN => {
println!("Problem status unknown.");
}
_ => {
println!("Other problem status.");
}
}
}
_ => {
println!("Other solution status");
}
}
Ok(())
}))
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
mioinitsol.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : mioinitsol.rs
//!
//! Purpose : Demonstrates how to solve a MIP with a start guess.
//!
extern crate mosek;
use mosek::{Task,Boundkey,Variabletype,Objsense,Streamtype,Soltype};
const INF : f64 = 0.0;
fn main() -> Result<(),String> {
let numvar : i32 = 4;
let numcon : i32 = 1;
let c = &[ 7.0, 10.0, 1.0, 5.0 ];
let bkc = &[Boundkey::UP];
let blc = &[ -INF ];
let buc = &[2.5];
let bkx = &[Boundkey::LO,
Boundkey::LO,
Boundkey::LO,
Boundkey::LO ];
let blx = &[0.0,
0.0,
0.0,
0.0 ];
let bux = &[INF,
INF,
INF,
INF];
let ptrb = &[0i64, 1, 2, 3];
let ptre = &[1i64, 2, 3, 4];
let aval = &[1.0, 1.0, 1.0, 1.0];
let asub = &[0i32, 0, 0, 0 ];
let intsub = &[0i32, 1, 2];
let inttype = &[Variabletype::TYPE_INT,
Variabletype::TYPE_INT,
Variabletype::TYPE_INT ];
let mut task = Task::new().unwrap().with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
task.input_data(numcon, numvar,
c, 0.0,
ptrb, ptre,
asub, aval,
bkc, blc, buc,
bkx, blx, bux)?;
task.put_var_type_list(intsub, inttype)?;
/* A maximization problem */
task.put_obj_sense(Objsense::MAXIMIZE)?;
// Assign values to integer variables
// We only set that slice of xx
task.put_xx_slice(Soltype::ITG, 0, 3, &[1.0,1.0,0.0])?;
// Request constructing the solution from integer variable values
task.put_int_param(mosek::Iparam::MIO_CONSTRUCT_SOL, mosek::Onoffkey::ON)?;
// solve
let _ = task.optimize()?;
task.solution_summary(mosek::Streamtype::LOG)?;
// Read and print solution
let mut xx = vec![0.0; numvar as usize];
task.get_xx(mosek::Soltype::ITG, xx.as_mut_slice())?;
println!("Optimal solution:");
for (i,xi) in xx.iter().enumerate() {
println!("x[{}] = {}",i,*xi);
}
// Was the initial solution used?
let constr = task.get_int_inf(mosek::Iinfitem::MIO_CONSTRUCT_SOLUTION)?;
let constr_val = task.get_dou_inf(mosek::Dinfitem::MIO_CONSTRUCT_SOLUTION_OBJ)?;
println!("Construct solution utilization: {}", constr);
println!("Construct solution objective: {}", constr_val);
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
opt_server_async.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : opt_server_async.rs
//!
//! Purpose : Demonstrates how to use MOSEK OptServer
//! to solve optimization problem asynchronously
extern crate mosek;
use mosek::{Task,Streamtype,Sparam};
use std::env;
use std::time::Duration;
use std::thread::sleep;
#[derive(Debug)]
enum FileOrText {
File(String),
Text(String)
}
fn main() {
let mut args = env::args();
if args.len() < 3 {
println!("Missing argument, syntax is:");
println!(" opt_server_async inputfile http[s]://HOSTNAME:PORT numpolls [certfile]");
panic!("Missing arguments");
}
let _ = args.next();
opt_server_async(FileOrText::File(args.next().unwrap()),
args.next().unwrap(),
args.next().unwrap().parse().unwrap(),
args.next()).unwrap();
}
fn opt_server_async(inputfile : FileOrText, addr : String, numpolls : usize, cert : Option<String>) -> Result<(),String> {
// Path to certificate, if any
let token = {
Task::new().unwrap()
.with_stream_callback(
Streamtype::LOG,
&mut |msg| print!("{}",msg),
|task| {
match inputfile {
FileOrText::File(ref filename) => task.read_data(filename.as_str()).unwrap(),
FileOrText::Text(ref data) => task.read_ptf_string(data.as_str()).unwrap()
}
if let Some(ref cert) = cert {
task.put_str_param(Sparam::REMOTE_TLS_CERT_PATH,cert.as_str())?;
}
task.async_optimize(addr.as_str(),"")
}).expect("Failed to submit async optimization")
};
println!("Task token = '{}'", token);
println!("Setting log stream...");
Task::new().unwrap().with_stream_callback(
Streamtype::LOG,
& mut |msg| print!("{}",msg),
|task| task.with_callback(
&mut|caller| { println!("caller = {}",caller); false },
|task| {
println!("Reading input file '{:?}'...",inputfile);
match inputfile {
FileOrText::File(ref filename) => task.read_data(filename.as_str()).unwrap(),
FileOrText::Text(ref data) => task.read_ptf_string(data.as_str()).unwrap()
}
if let Some(ref cert) = cert {
task.put_str_param(Sparam::REMOTE_TLS_CERT_PATH,cert.as_str())?;
}
println!("Starting polling loop...");
for i in 0..numpolls {
sleep(Duration::new(1,0));
println!("\tpoll {}...", i);
let mut trm : i32 = 0;
let mut resp : i32 = 0;
let respavailable = task.async_poll(addr.as_str(),
"",
token.as_str(),
& mut resp,
& mut trm)?;
if respavailable {
println!("solution available!");
task.async_get_result(addr.as_str(),
"",
token.as_str(),
& mut resp,
& mut trm)?;
task.solution_summary (Streamtype::LOG)?;
return Ok(());
}
}
println!("max num polls reached, stopping host.");
task.async_stop (addr.as_str(), "", token.as_str())?;
Err("Max num polls".to_string())
}))
}
#[cfg(test)]
mod tests {
const DFLT_FILE : &str = "Task
Objective
Maximize + 2 @x0 + 3 @x1 - @x2
Constraints
@c0 [1] + @x0 + @x1 + @x2
@C0 [QUAD(3)]
@ac1: + 0.03
@ac2: + 1.5 @x0 + 0.1 @x1
@ac3: + 0.3 @x0 + 2.1 @x2 + 0.1
Variables
@x0
@x1
@x2
";
#[test]
fn test() {
super::opt_server_async(super::FileOrText::Text(DFLT_FILE.to_string()),
"http://solve.mosek.com:30080".to_string(),
100,
None).unwrap();
}
}
opt_server_sync.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : opt_server_sync.rs
//!
//! Purpose : Demonstrates how to use MOSEK OptServer
//! to solve optimization problem synchronously
//!
extern crate mosek;
use mosek::{Task,Streamtype,Sparam};
use std::env;
enum FileOrText {
File(String),
Text(String)
}
fn main() -> Result<(),String> {
let mut args = env::args();
if args.len() < 3 {
println!("Missing argument, syntax is:");
println!(" opt_server_sync inputfile http[s]://HOSTNAME:PORT [certfile]");
return Err("Missing arguments".to_string())
}
let _ = args.next();
opt_server_sync(FileOrText::File(args.next().unwrap()),
args.next().unwrap(),
args.next())
}
fn opt_server_sync(inputfile : FileOrText, addr : String, cert : Option<String>) -> Result<(),String> {
let mut task = Task::new().unwrap().with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
// Load some data into the task
match inputfile {
FileOrText::File(filename) => task.read_data(filename.as_str())?,
FileOrText::Text(data) => task.read_ptf_string(data.as_str())?
}
// Set OptServer URL
task.put_optserver_host(addr.as_str())?;
// Path to certificate, if any
if let Some(cert) = cert {
task.put_str_param(Sparam::REMOTE_TLS_CERT_PATH, cert.as_str())?;
}
// Optimize remotely, no access token
let _trm = task.optimize()?;
task.solution_summary(Streamtype::LOG)?;
Ok(())
}
#[cfg(test)]
mod tests {
const DFLT_FILE : &str = "Task
Objective
Maximize + 2 @x0 + 3 @x1 - @x2
Constraints
@c0 [1] + @x0 + @x1 + @x2
@C0 [QUAD(3)]
@ac1: + 0.03
@ac2: + 1.5 @x0 + 0.1 @x1
@ac3: + 0.3 @x0 + 2.1 @x2 + 0.1
Variables
@x0
@x1
@x2
";
#[test]
fn test() {
super::opt_server_sync(super::FileOrText::Text(DFLT_FILE.to_string()),
"http://solve.mosek.com:30080".to_string(),
None).unwrap();
}
}
parallel.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : parallel.rs
//!
//! Purpose: Demonstrates parallel optimization using optimizebatch()
//!
extern crate mosek;
extern crate itertools;
use mosek::{Task};
use std::env;
use itertools::{izip};
/// Example of how to use env.optimize_batch().
/// Optimizes tasks whose names were read from command line.
enum FileOrText {
File(String),
Text(String)
}
fn main() -> Result<(),String> {
let mut args = env::args();
if args.len() < 3 {
println!("Syntax: parallel FILENAME FILENAME [ FILENAME ... ]");
Err("Invalid argument list".to_string())
}
else {
let _ = args.next();
parallel(args.map(|s| FileOrText::File(s)).collect())
}
}
fn parallel(files : Vec<FileOrText>) -> Result<(),String> {
// Create an example list of tasks to optimize
let mut tasks : Vec<(String,Task)> = files.iter().filter_map(|fname| {
let mut t = Task::new().unwrap();
match fname {
FileOrText::File(fname) => {
if let Err(_) = t.read_data(fname.as_str()) { None }
else {
t.put_int_param(mosek::Iparam::NUM_THREADS, 2).unwrap();
Some((fname.as_str().to_string(),t))
}
},
FileOrText::Text(data) => {
if let Err(_) = t.read_ptf_string(data.as_str()) { None }
else {
t.put_int_param(mosek::Iparam::NUM_THREADS, 2).unwrap();
Some(("-".to_string(),t))
}
}
}
}).collect();
let mut res = vec![0i32; tasks.len()];
let mut trm = vec![0i32; tasks.len()];
{
let taskrs : Vec<& mut Task> = tasks.iter_mut().map(|(_a,b)| b).collect();
// Size of thread pool available for all tasks
let threadpoolsize : i32 = 6;
// Optimize all the given tasks in parallel
mosek::optimize_batch(false, // No race
-1.0, // No time limit
threadpoolsize,
taskrs.as_slice(), // Array of tasks to optimize
trm.as_mut_slice(),
res.as_mut_slice())?;
}
for (resi,trmi,(fname,ti)) in izip!(res,trm,tasks.iter()) {
println!("Task {} res {} trm {} obj_val {} time {}",
fname,
resi,
trmi,
ti.get_dou_inf(mosek::Dinfitem::INTPNT_PRIMAL_OBJ)?,
ti.get_dou_inf(mosek::Dinfitem::OPTIMIZER_TIME)?);
}
Ok(())
}
#[cfg(test)]
mod tests {
const DFLT_FILE1 : &str = "Task
Objective
Maximize + 2 @x0 + 3 @x1 - @x2
Constraints
@c0 [1] + @x0 + @x1 + @x2
@C0 [QUAD(3)]
@ac1: + 0.03
@ac2: + 1.5 @x0 + 0.1 @x1
@ac3: + 0.3 @x0 + 2.1 @x2 + 0.1
Variables
@x0
@x1
@x2
";
const DFLT_FILE2 : &str = "Task
Objective
Maximize + @x0 + 0.64 @x1
Constraints
@c0 [-inf;250] + 50 @x0 + 31 @x1
@c1 [-4;+inf] + 3 @x0 - 2 @x1
Variables
@x0[0;+inf]
@x1[0;+inf]
Integers
@x0 @x1
";
#[test]
fn test() {
super::parallel(vec![super::FileOrText::Text(DFLT_FILE1.to_string()),
super::FileOrText::Text(DFLT_FILE2.to_string())]).unwrap();
}
}
parameters.rs
//
// Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File : parameters.rs
//
// Purpose : Demonstrates a very simple example about how to get/set
// parameters with MOSEK Julia API
//
extern crate mosek;
use mosek::{Task,Iparam,Dparam,Dinfitem,Iinfitem,Optimizertype,Basindtype};
fn main() -> Result<(),String> {
let mut task = Task::new().unwrap();
println!("Test MOSEK parameter get/set functions");
// Set log level (integer parameter)
task.put_int_param(Iparam::LOG, 1)?;
// Select interior-point optimizer... (integer parameter)
task.put_int_param(Iparam::OPTIMIZER, Optimizertype::INTPNT)?;
// ... without basis identification (integer parameter)
task.put_int_param(Iparam::INTPNT_BASIS,Basindtype::NEVER)?;
// Set relative gap tolerance (double parameter)
task.put_dou_param(Dparam::INTPNT_CO_TOL_REL_GAP, 1.0e-7)?;
// The same using explicit string names
task.put_param("MSK_DPAR_INTPNT_CO_TOL_REL_GAP", "1.0e-7")?;
task.put_na_dou_param("MSK_DPAR_INTPNT_CO_TOL_REL_GAP", 1.0e-7 )?;
// Incorrect value
if let Err(_) = task.put_dou_param(Dparam::INTPNT_CO_TOL_REL_GAP, -1.0) {
println!("Wrong parameter value");
}
let param = task.get_dou_param(Dparam::INTPNT_CO_TOL_REL_GAP)?;
println!("Current value for parameter intpnt_co_tol_rel_gap = $param");
// Define and solve an optimization problem here
// optimize(task,)
// After optimization:
println!("Get MOSEK information items");
let tm = task.get_dou_inf(Dinfitem::OPTIMIZER_TIME)?;
let iter = task.get_int_inf(Iinfitem::INTPNT_ITER)?;
println!("Time: {}",tm);
println!("Iterations: {}",iter);
Ok(())
}
pinfeas.rs
//!
//! File : pinfeas.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! Purpose: Demonstrates how to fetch a primal infeasibility certificate
//! for a linear problem
//!
extern crate mosek;
use mosek::{Streamtype,Boundkey,Soltype,Prosta};
const INF : f64 = 0.0;
fn test_problem() -> Result<mosek::Task,String> {
let mut task = mosek::Task::new().unwrap();
task.append_vars(7)?;
task.append_cons(7)?;
task.put_c_list(&[0,1,2,3,4,5,6],
&[1.0,2.0,5.0,2.0,1.0,2.0,1.0])?;
task.put_aij_list(&[0,0,1,1,2,2,2,3,3,4,5,5,6,6],
&[0,1,2,3,4,5,6,0,4,1,2,5,3,6],
&[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0])?;
task.put_con_bound_slice(0, 7,
&[Boundkey::UP,Boundkey::UP,Boundkey::UP,Boundkey::FX,Boundkey::FX,Boundkey::FX,Boundkey::FX],
&[-INF, -INF, -INF, 1100.0, 200.0, 500.0, 500.0],
&[200.0, 1000.0, 1000.0, 1100.0, 200.0, 500.0, 500.0])?;
task.put_var_bound_slice_const(0, 7, Boundkey::UP, 0.0, INF)?;
Ok(task)
}
// Analyzes and prints infeasibility contributing elements
// sl - dual values for lower bounds
// su - dual values for upper bounds
// eps - tolerance for when a nunzero dual value is significant
fn analyze_certificate(sl : &[f64], su : &[f64], eps : f64) {
for (i,(&sli,&sui)) in sl.iter().zip(su.iter()).enumerate() {
if sli > eps {
println!("#{}, lower, dual = {:e}\n", i, sli);
}
if sui > eps {
println!("#{}, upper, dual = {:e}\n", i, sui);
}
}
}
fn main() -> Result<(),String> {
// In this example we set up a simple problem
// One could use any task or a task read from a file
let mut task = test_problem()?.with_callbacks();
let n = task.get_num_var()?;
let m = task.get_num_con()?;
// Useful for debugging
task.write_data("pinfeas.ptf")?; // Write file in human-readable format
// Attach a log stream printer to the task
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
// Perform the optimization.
task.optimize()?;
task.solution_summary(Streamtype::LOG)?;
// Check problem status, we use the interior point solution
if task.get_pro_sta(Soltype::ITR)? == Prosta::PRIM_INFEAS {
// Set the tolerance at which we consider a dual value as essential
let eps = 1e-7;
println!("Variable bounds important for infeasibility: ");
let mut slx = vec![0.0; n as usize]; task.get_slx(Soltype::ITR, slx.as_mut_slice())?;
let mut sux = vec![0.0; n as usize]; task.get_sux(Soltype::ITR, sux.as_mut_slice())?;
analyze_certificate(slx.as_slice(), sux.as_slice(), eps);
println!("Constraint bounds important for infeasibility: ");
let mut slc = vec![0.0; m as usize]; task.get_slc(Soltype::ITR, slc.as_mut_slice())?;
let mut suc = vec![0.0; m as usize]; task.get_suc(Soltype::ITR, suc.as_mut_slice())?;
analyze_certificate(slc.as_mut_slice(), suc.as_mut_slice(), eps);
}
else {
println!("The problem is not primal infeasible, no certificate to show");
}
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
portfolio_1_basic.rs
//!
//! File : portfolio_1_basic.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! Description : Implements a basic portfolio optimization model.
//!
//! More details can be found at <https://docs.mosek.com/latest/capi/case-portfolio.html#doc-optimizer-case-portfolio>
//!
extern crate mosek;
extern crate itertools;
use mosek::{Task,Objsense,Streamtype,Soltype,Solsta};
use itertools::{iproduct};
/// Solve basic Markowitz portfolio problem: Maximize expected return
/// while bounded the estimated risk.
///
/// ```
/// Maximize mu'x
/// Subject to
/// budget : sum(x) = sum(x0)+w
/// risk: gamma^2 >= || G'x ||^2
/// x >= 0
/// ```
///
/// # Arguments
///
/// - `n` number of assets
/// - `gamma` risk bound (bound on the standard deviation)
/// - `mu` vector of expected returns
/// - `GT` Covariance matrix factor
/// - `x0` vector if initial investment
/// - `w` initial uninvested wealth
#[allow(non_snake_case)]
fn portfolio(n : i32, // number of assets
gamma : f64, // risk bound: maximum stddev
mu : &[f64], // vector of expected returns
GT : &[f64], // covariance matrix factor
x0 : &[f64], // initial investment
w : f64) // initial wealth
-> Result<(Vec<f64>,f64),String> {
let k = (GT.len() / n as usize) as i32;
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* Total budget */
let total_budget = w + x0.iter().sum::<f64>();
/* Constraints. */
task.append_cons(1i32)?;
/* Variables. */
task.append_vars(n)?;
let x : Vec<i32> = (0i32..n).collect();
/* Total budget constraint - set bounds l^c = u^c */
task.put_con_bound(0i32, mosek::Boundkey::FX, total_budget, total_budget)?;
task.put_con_name(0i32,"budget")?;
task.put_obj_sense(Objsense::MAXIMIZE)?;
task.put_c_slice(0,n,mu)?;
/* x variables. */
for (j,xj) in x.iter().enumerate() {
/* Coefficients in the first row of A */
task.put_aij(0, *xj, 1.0)?;
/* No short-selling - x^l = 0, x^u = inf */
task.put_var_bound(*xj, mosek::Boundkey::LO, 0.0, 0.0)?;
task.put_var_name(*xj, format!("x[{}]",j+1).as_str())?;
}
// risk bound
{
let acci = task.get_num_acc()?;
let afei = task.get_num_afe()?;
task.append_afes(k as i64 + 1)?;
let dom = task.append_quadratic_cone_domain(k as i64+1)?;
task.append_acc_seq(dom,
afei,
vec![0.0; k as usize + 1].as_slice())?;
task.put_acc_name(acci,"risk")?;
task.put_afe_g(afei,gamma)?;
for ((i,j),v) in iproduct!(0..n,0..n).zip(GT).filter(|(_,v)| **v != 0.0) {
task.put_afe_f_entry(afei + i as i64 + 1, j as i32, *v)?;
}
}
/* Dump the problem to a human readable OPF file. */
// task.write_data("portfolio_1_basic.ptf")?;
let _trm = task.optimize()?;
// Check if the interior point solution is an optimal point
if task.get_sol_sta(Soltype::ITR)? != Solsta::OPTIMAL {
// See https://docs.mosek.com/latest/rustapi/accessing-solution.html about handling solution statuses.
eprintln!("Solution not optimal!");
std::process::exit(1);
}
/* Display the solution summary for quick inspection of results. */
task.solution_summary(Streamtype::MSG)?;
task.write_data("portfolio_1_basic.ptf")?;
/* Read the x variables one by one and compute expected return. */
/* Can also be obtained as value of the objective. */
let mut level = vec![0.0;n as usize];
task.get_xx_slice(Soltype::ITR,0,n,level.as_mut_slice())?;
let expret = task.get_primal_obj(Soltype::ITR)?;
Ok((level.to_vec(),expret))
}
#[allow(non_snake_case)]
fn main() -> Result<(),String> {
let n = 8i32;
let w = 59.0;
let mu = &[0.07197349, 0.15518171, 0.17535435, 0.0898094 , 0.42895777, 0.39291844, 0.32170722, 0.18378628];
let x0 = &[8.0, 5.0, 3.0, 5.0, 2.0, 9.0, 3.0, 6.0];
let gamma = 36.0;
let 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 ];
let (level,expret) = portfolio(n,gamma,mu,GT,x0,w)?;
println!("Solution x = {:?}",level);
println!("Expected return {:.4e} for gamma = {:.4e}", expret, gamma);
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
portfolio_2_frontier.rs
//!
//! File : portfolio_2_frontier.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! Description : Implements a basic portfolio optimization model.
//! Determines points on the efficient frontier.
//!
extern crate mosek;
extern crate itertools;
use mosek::{Task,Objsense,Solsta,Soltype};
use itertools::{iproduct};
/// Solve basic Markowitz portfolio problem for different risk bounds
/// to produce a list of points on the efficient frontier.
///
/// ```
/// Maximize mu'x - alpha * s
/// Subject to
/// budget : sum(x) = sum(x0)+w
/// risk: s >= || G'x ||
/// x >= 0
/// ```
///
/// # Arguments
///
/// - `n` number of assets
/// - `alphas` list of risk bound (bound on the standard deviation)
/// - `mu` vector of expected returns
/// - `GT` Covariance matrix factor
/// - `x0` vector if initial investment
/// - `w` initial uninvested wealth
#[allow(non_snake_case)]
fn portfolio(n : i32,
mu : &[f64],
GT : &[f64],
x0 : &[f64],
alphas : &[f64],
w : f64) -> Result<Vec<(f64,f64)>,String> {
let k = (GT.len() / n as usize) as i32;
/* Create the optimization task. */
let mut task = match Task::new() {
Some(t) => t,
None => return Err("Failed to create task".to_string()),
};
//task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
task.append_vars(n+1)?;
task.append_cons(1)?;
/* Objective */
task.put_obj_sense(Objsense::MAXIMIZE)?;
let x : Vec<i32> = (0i32..n).collect();
let s = n;
/* Total budget */
let total_budget = w + x0.iter().sum::<f64>();
/* Total budget constraint - set bounds l^c = u^c */
task.put_con_bound(0i32, mosek::Boundkey::FX, total_budget, total_budget)?;
task.put_con_name(0i32,"budget")?;
task.put_c_slice(0,n,mu)?;
/* x variables. */
for (j,xj) in x.iter().enumerate() {
/* Coefficients in the first row of A */
task.put_aij(0, *xj, 1.0)?;
/* No short-selling - x^l = 0, x^u = inf */
task.put_var_bound(*xj, mosek::Boundkey::LO, 0.0, 0.0)?;
task.put_var_name(*xj, format!("x[{}]",j+1).as_str())?;
}
task.put_var_name(s, "s")?;
task.put_var_bound(s, mosek::Boundkey::FR, 0.0, 0.0)?;
// risk bound
// (s,0.5,GT * x) in Q_r
{
let acci = task.get_num_acc()?;
let afei = task.get_num_afe()?;
task.append_afes(k as i64 + 2)?;
let dom = task.append_r_quadratic_cone_domain(k as i64+2)?;
task.append_acc_seq(dom,
afei,
vec![0.0; k as usize + 2].as_slice())?;
task.put_acc_name(acci,"risk")?;
task.put_afe_f_entry(afei,s,1.0)?;
task.put_afe_g(afei+1,0.5)?;
for ((i,j),v) in iproduct!(0..n,0..n).zip(GT).filter(|(_,v)| **v != 0.0) {
task.put_afe_f_entry(afei + i as i64 + 2, j as i32, *v)?;
}
}
let frontier : Vec<(f64,f64)> = alphas.iter().filter_map(|alpha| {
/* Sets the objective function coefficient for s. */
if let Err(_) = task.put_c_j(s, - *alpha) { None }
else if let Err(_) = task.optimize() { None }
else if let Err(_) = task.write_data(format!("portfolio_2_frontier-{}.ptf",alpha).as_str()) { None }
else if let Ok(solsta) = task.get_sol_sta(Soltype::ITR) {
// See https://docs.mosek.com/latest/rustapi/accessing-solution.html about handling solution statuses.
match solsta {
Solsta::OPTIMAL => {
let mut xx = vec![0.0; n as usize+1];
if let Err(_) = task.get_xx(Soltype::ITR,xx.as_mut_slice()) { None }
else {
Some((*alpha,mu.iter().zip(xx.iter()).map(|(m,x)| m * x).sum::<f64>()))
}
}
_ => None
}
}
else {
None
}
}).collect();
Ok(frontier)
}
#[allow(non_snake_case)]
fn main() -> Result<(),String> {
let n : i32 = 8;
let w = 1.0;
let mu = &[0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379];
let x0 = &[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
let 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 ];
let 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];
println!("{:10} {:10}","alpha","exp.ret.");
for (alpha,expret) in portfolio(n,
mu,
GT,
x0,
alphas,
w)? {
println!("{:10.3e} : {:10.e}",alpha,expret);
}
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
portfolio_3_impact.rs
//!
//! File : portfolio_3_impact.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! Description : Implements a basic portfolio optimization model with transaction costs of order x^(3/2).
//!
//! More details can be found at <https://docs.mosek.com/latest/capi/case-portfolio.html#doc-optimizer-case-portfolio>
//!
extern crate mosek;
extern crate itertools;
use mosek::{Task,Objsense,Streamtype,Solsta,Soltype,Boundkey};
use itertools::{izip,iproduct};
const INF : f64 = 0.0;
/// Solve portfolio with market impact terms.
///
/// ```
/// Maximize mu'x
/// Subject to
/// budget : sum(x)+m'c = sum(x0)+w
/// risk : (gamma,G'x) in Q^{k+1}
/// MI : (c_j,1,|x_j-x0_j|) in P^3(2/3,1/3), j = 1..
/// x >= 0
/// ```
///
/// Where
///
/// - m_i is the transaction cost associated with asset i
/// - gamma is the bound on the standard deviation if the portfolio
/// - mu_i is the expected return on asset i
/// - w is the initial wealth held in cash
/// - x0_i is the initial investment in asset i
/// - G'G is the covariance matrix for assets
///
/// The MI constraint is not convex due tot he |.| term, so we relax it:
/// ```
/// MI : (c_j,1,z_j) in P^3(2/3,1/3), j = 1..
/// z_j >= |x_j-x0_j|
/// implemented as
/// z_j >= x_j-x0_j
/// z_j >= x0_j-x_j
/// ```
///
/// # Arguments
///
/// - `n` number of assets
/// - `mu` vector of expected returns
/// - `m` vector of market impact estimates
/// - `GT` factored covariance matrix
/// - `x0` vector if initial investment
/// - `gamma` bound on risk
/// - `w` initial uninvested wealth
///
/// # Returns
///
/// Returns `(solution,objval)`
///
/// - `solution` is the primal investment solution vector
/// - `objval` is the solution expected return
#[allow(non_snake_case)]
pub fn portfolio(n : i32,
mu : &[f64],
m : &[f64],
GT : &[f64],
x0 : &[f64],
gamma : f64,
w : f64) -> Result<(Vec<f64>,f64),String> {
let k = (GT.len() / n as usize) as i32;
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
task.append_vars(3*n)?;
let allvars : Vec<i32> = (0i32..3*n).collect();
let var_x = &allvars[0..n as usize];
let var_c = &allvars[n as usize..2*n as usize];
let var_xc = &allvars[0..2*n as usize];
let var_z = &allvars[2*n as usize..3*n as usize];
for (i,j) in var_x.iter().enumerate() {
task.put_var_bound(*j,mosek::Boundkey::LO, 0.0, 0.0)?;
task.put_var_name(*j,format!("x[{}]",i+1).as_str())?; }
for (i,j) in var_c.iter().enumerate() {
task.put_var_bound(*j,mosek::Boundkey::FR, 0.0, 0.0)?;
task.put_var_name(*j,format!("c[{}]",i+1).as_str())?;
}
for (i,j) in var_z.iter().enumerate() {
task.put_var_bound(*j,mosek::Boundkey::FR, 0.0, 0.0)?;
task.put_var_name(*j,format!("z[{}]",i+1).as_str())?;
}
task.put_obj_sense(Objsense::MAXIMIZE)?;
for i in var_x {
task.put_c_j(*i,mu[*i as usize])?;
}
task.append_cons(1)?;
let con_budget = 0i32;
// budget
task.put_con_name(0,"budget")?;
let wealth = w + x0.iter().sum::<f64>();
task.put_a_row(con_budget,
&var_xc,
(0..n).map(|_| 1.0).chain(m.iter().map(|v| *v)).collect::<Vec<f64>>().as_slice())?;
task.put_con_bound(con_budget,mosek::Boundkey::FX, wealth,wealth)?;
// |x-x0| <= z
{
let coni = task.get_num_con()?;
task.append_cons(2 * n)?;
for i in 0..n {
task.put_con_name(coni+i, format!("zabs1[{}]",1 + i).as_str())?;
task.put_con_name(coni+n+i, format!("zabs2[{}]",1 + i).as_str())?;
}
let ones = vec![1.0; n as usize];
let minusones = vec![-1.0; n as usize];
let con_abs1 : Vec<i32> = (coni..coni+n).collect();
let con_abs2 : Vec<i32> = (coni+n..coni+2*n).collect();
task.put_aij_list(con_abs1.as_slice(), var_x, minusones.as_slice())?;
task.put_aij_list(con_abs1.as_slice(), var_z, ones.as_slice())?;
task.put_con_bound_slice(coni,coni+n, vec![Boundkey::LO; n as usize].as_slice(), x0.iter().map(|&v| -v).collect::<Vec<f64>>().as_slice(), vec![INF; n as usize].as_slice())?;
task.put_aij_list(con_abs2.as_slice(), var_x, ones.as_slice())?;
task.put_aij_list(con_abs2.as_slice(), var_z, ones.as_slice())?;
task.put_con_bound_slice(coni+n,coni+n*2, vec![Boundkey::LO; n as usize].as_slice(), x0, vec![INF; n as usize].as_slice())?;
}
// GT
{
let acci = task.get_num_acc()?;
let afei = task.get_num_afe()?;
task.append_afes(k as i64 + 1)?;
let dom = task.append_quadratic_cone_domain(k as i64+1)?;
task.append_acc_seq(dom,
afei,
vec![0.0; k as usize + 1].as_slice())?;
task.put_acc_name(acci,"risk")?;
task.put_afe_g(afei,gamma)?;
for ((i,j),v) in iproduct!(0..n,0..n).zip(GT).filter(|(_,v)| **v != 0.0) {
task.put_afe_f_entry(afei + i as i64 + 1, j as i32, *v)?;
}
}
// MI
{
let mut acci = task.get_num_acc()?;
let mut afei = task.get_num_afe()?;
let afe0 = afei;
task.append_afes(n as i64 * 2+1)?;
let dom = task.append_primal_power_cone_domain(3,&[2.0, 1.0])?;
task.put_afe_g(afe0,1.0)?;
afei += 1;
for (i,&cj,&zj,&x0j) in izip!(0..n,var_c,var_z,x0) {
task.put_afe_f_entry(afei,cj,1.0)?;
task.put_afe_f_entry(afei+1,zj,1.0)?;
task.put_afe_g(afei+1, - x0j)?;
task.append_acc(dom,
&[afei,afe0,afei+1],
&[0.0, 0.0, 0.0])?;
task.put_acc_name(acci,format!("market_impact[{}]",i+1).as_str())?;
afei += 2;
acci += 1;
}
}
let _ = task.optimize()?;
task.write_data("portfolio_3_impact.ptf")?;
/* Display the solution summary for quick inspection of results. */
task.solution_summary(Streamtype::MSG)?;
if ! task.solution_def(Soltype::ITR)? {
return Err("No solultion defined".to_string());
}
// See https://docs.mosek.com/latest/rustapi/accessing-solution.html about handling solution statuses.
let solsta = task.get_sol_sta(Soltype::ITR)?;
if solsta != Solsta::OPTIMAL {
return Err("Unexpected solution status".to_string());
}
let mut level = vec![0.0;n as usize];
task.get_xx_slice(Soltype::ITR,0,n,level.as_mut_slice())?;
let obj = task.get_primal_obj(Soltype::ITR)?;
Ok((level,obj))
}
#[allow(non_snake_case)]
fn main() -> Result<(),String> {
let n : i32 = 8;
let w = 1.0;
let mu = &[0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379];
let x0 = &[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
let 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 ];
let gamma = 0.36;
let m = vec![0.01; n as usize];
let (level,obj) = portfolio(n,
mu,
m.as_slice(),
GT,
x0,
gamma,
w)?;
println!("Solution x = {:?}",level);
println!("Objective value x = {:?}",obj);
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
portfolio_4_transcost.rs
//! File : portfolio_4_transcost.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! Description : Implements a basic portfolio optimization model
//! with fixed setup costs and transaction costs
//! as a mixed-integer problem.
//!
extern crate mosek;
use mosek::{Task,Objsense,Streamtype,Soltype,Variabletype,Boundkey,Solsta};
extern crate itertools;
use itertools::{iproduct};
const INF : f64 = 0.0;
/// Optimize expected return on an investment with transaction cost.
///
/// ```
/// Maximize mu'x
/// Such That
/// budget: sum(x) + f'y + g'z = w0+sum(x0)
/// risk: gamma > || G'x ||
/// ACC: z_j > |x0_j-x_j|
/// DJC: [ y_j == 0 AND x0_j == x_j ] OR y_j == 1
/// // z_j < U_j y_j, y in {0,1}
/// y free
/// x > 0
/// Where f_i is the fixed cost of a transaction in asset i,
/// g_i is the cost per unit of a transaction in asset i
/// ```
/// # Arguments
///
/// - `n` number of assets
/// - `mu` vector of expected returns
/// - `f` vector of fixed transaction costs
/// - 'g' vector of continuous proportional transaction costs
/// - `GT` Covariance matrix factor
/// - `x0` vector if initial investment
/// - `gamma` risk bound (bound on the standard deviation)
/// - `w` initial uninvested wealth
#[allow(non_snake_case)]
fn portfolio(n : i32,
mu : &[f64],
f : &[f64],
g : &[f64],
GT : &[f64],
x0 : &[f64],
gamma : f64,
w : f64) -> Result<(Vec<f64>,f64),String> {
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
let k = (GT.len() / n as usize) as i32;
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* Compute total wealth */
let w0 = w + x0.iter().sum::<f64>();
task.append_cons(1i32)?;
task.append_vars(3*n)?;
for i in 0i32..n {
task.put_var_name(i, format!("x[{}]",i+1).as_str())?;
task.put_var_name(i+n, format!("y[{}]",i+1).as_str())?;
task.put_var_name(i+2*n,format!("z[{}]",i+1).as_str())?;
task.put_var_type(i+n, Variabletype::TYPE_INT)?;
}
let all_vars : Vec<i32> = (0i32..3*n).collect();
let x = &all_vars[0..n as usize];
let y = &all_vars[n as usize..2*n as usize];
let z = &all_vars[2*n as usize..3*n as usize];
task.put_var_bound_slice_const(0i32,n, mosek::Boundkey::LO, 0.0,0.0)?;
task.put_var_bound_slice_const(n,2*n, mosek::Boundkey::RA, 0.0,1.0)?;
task.put_var_bound_slice_const(2*n,3*n, mosek::Boundkey::FR, 0.0,0.0)?;
/* Constraints. */
task.put_con_name(0,"budget")?;
{
let zeros = vec![0i32; n as usize];
let fones = vec![1.0; n as usize];
task.put_aij_list(zeros.as_slice(), x, fones.as_slice())?;
task.put_aij_list(zeros.as_slice(), y, f)?;
task.put_aij_list(zeros.as_slice(), z, g)?;
task.put_con_bound(0i32,mosek::Boundkey::FX,w0,w0)?;
}
// objective
task.put_obj_sense(Objsense::MAXIMIZE)?;
for (xi,mui) in x.iter().zip(mu.iter()) {
task.put_c_j(*xi, *mui)?;
}
// risk bound
{
let acci = task.get_num_acc()?;
let afei = task.get_num_afe()?;
task.append_afes(k as i64 + 1)?;
let dom = task.append_quadratic_cone_domain(k as i64+1)?;
task.append_acc_seq(dom,
afei,
vec![0.0; k as usize + 1].as_slice())?;
task.put_acc_name(acci,"risk")?;
task.put_afe_g(afei,gamma)?;
for ((i,j),v) in iproduct!(0..n,0..n).zip(GT).filter(|(_,v)| **v != 0.0) {
task.put_afe_f_entry(afei + i as i64 + 1, j as i32, *v)?;
}
}
// |x-x0| <= z
{
let coni = task.get_num_con()?;
task.append_cons(2 * n)?;
for i in 0..n {
task.put_con_name(coni+i, format!("zabs1[{}]",1 + i).as_str())?;
task.put_con_name(coni+n+i, format!("zabs2[{}]",1 + i).as_str())?;
}
let ones = vec![1.0; n as usize];
let minusones = vec![-1.0; n as usize];
let con_abs1 : Vec<i32> = (coni..coni+n).collect();
let con_abs2 : Vec<i32> = (coni+n..coni+2*n).collect();
task.put_aij_list(con_abs1.as_slice(), x, minusones.as_slice())?;
task.put_aij_list(con_abs1.as_slice(), z, ones.as_slice())?;
task.put_con_bound_slice(coni,coni+n, vec![Boundkey::LO; n as usize].as_slice(), x0.iter().map(|&v| -v).collect::<Vec<f64>>().as_slice(), vec![INF; n as usize].as_slice())?;
task.put_aij_list(con_abs2.as_slice(), x, ones.as_slice())?;
task.put_aij_list(con_abs2.as_slice(), z, ones.as_slice())?;
task.put_con_bound_slice(coni+n,coni+n*2, vec![Boundkey::LO; n as usize].as_slice(), x0, vec![INF; n as usize].as_slice())?;
}
// Switch
{
let coni = task.get_num_con()?;
task.append_cons(n)?;
for i in 0..n {
task.put_con_name(coni + i, format!("switch[{}]",i+1).as_str())?;
}
let conlist : Vec<i32> = (coni..coni+n).collect();
task.put_aij_list(conlist.as_slice(), z, vec![1.0; n as usize].as_slice())?;
task.put_aij_list(conlist.as_slice(), y, vec![-w0; n as usize].as_slice())?;
task.put_con_bound_slice_const(coni,coni+n, Boundkey::UP, 0.0,0.0)?;
}
let _ = task.optimize()?;
task.write_data("portfolio_4_transcost.ptf")?;
/* Display the solution summary for quick inspection of results. */
task.solution_summary(Streamtype::MSG)?;
// Check if the integer solution is an optimal point
if task.get_sol_sta(Soltype::ITG)? != Solsta::INTEGER_OPTIMAL {
// See https://docs.mosek.com/latest/rustapi/accessing-solution.html about handling solution statuses.
eprintln!("Solution not optimal!");
std::process::exit(1);
}
let mut xx = vec![0.0;n as usize];
task.get_xx_slice(Soltype::ITG, 0,n, xx.as_mut_slice())?;
let expret = xx[0..n as usize].iter().zip(mu.iter()).map(|(a,b)| a*b).sum::<f64>();
Ok((xx,expret))
}
#[allow(non_snake_case)]
fn main() -> Result<(),String> {
let n = 8i32;
let w = 1.0;
let mu = &[0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379];
let x0 = &[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
let 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 ];
let f = vec![0.01; n as usize];
let g = vec![0.001; n as usize];
let gamma = 0.36;
let (level,expret) = portfolio(n,
mu,
f.as_slice(),
g.as_slice(),
GT,
x0,
gamma,
w)?;
println!("Expected return {:.4e} for gamma {:.4e}\n", expret, gamma);
println!("Solution vector = {:?}\n", level);
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
portfolio_5_card.rs
//!
//! File : portfolio_5_card.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! Description : Implements a basic portfolio optimization model
//! with cardinality constraints on number of assets traded.
//!
extern crate mosek;
use mosek::{Task,Objsense,Soltype,Variabletype,Boundkey,Solsta};
extern crate itertools;
use itertools::{izip,iproduct};
const INF : f64 = 0.0;
/// Solve the Markowitz portfolio problem with cardinality constraints.
///
/// ```
/// Maximize mu'x
/// Such That
/// budget: 𝐞'x = w0+sum(x0)
/// risk: γ^2 > ||G'x||^2
/// cardinality: 𝐞'y < p
/// [ y_j == 0 AND x0_j == x_j ] OR y_j == 1
/// y ∈ R^n
/// y ∈ R+^n
/// ```
///
/// Where
///
/// - `y_j` ∈ `{0,1}` is an indicator of whether we change the
/// investment in asset j, that is if `y_j=0`, then we do not change
/// investment in assert `j`.
///
/// # Arguments
///
/// - `n` number of assets
/// - `mu` vector of expected returns
/// - `GT` Covariance matrix factor
/// - `x0` vector if initial investment
/// - `gamma` risk bound (bound on the standard deviation)
/// - `w` initial uninvested wealth
/// - `p` maximum number of assets to invest in
#[allow(non_snake_case)]
fn portfolio(n : i32,
mu : &[f64],
GT : &[f64],
x0 : &[f64],
gamma : f64,
p : i32,
w : f64) -> Result<(Vec<f64>,f64),String> {
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
};
let k = (GT.len() / n as usize) as i32;
// task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* Compute total wealth */
let w0 = w + x0.iter().sum::<f64>();
task.append_vars(3*n)?;
let all_vars : Vec<i32> = (0..3*n).collect();
let x = &all_vars[0..n as usize];
let y = &all_vars[n as usize..2*n as usize];
let z = &all_vars[2*n as usize..3*n as usize];
task.put_var_bound_slice_const(0,n,mosek::Boundkey::LO,0.0,INF)?;
task.put_var_bound_slice_const(n,2*n,mosek::Boundkey::RA,0.0,1.0)?;
task.put_var_bound_slice_const(2*n,3*n, mosek::Boundkey::FR, -INF,INF)?;
for (i,xj,yj,zj) in izip!(0..n,x,y,z) {
task.put_var_name(*xj,format!("x[{}]",i+1).as_str())?;
task.put_var_name(*yj,format!("y[{}]",i+1).as_str())?;
task.put_var_name(*zj,format!("z[{}]",i+1).as_str())?;
task.put_var_type(*yj, Variabletype::TYPE_INT)?;
}
// objective
task.put_obj_sense(Objsense::MAXIMIZE)?;
for (j,mui) in x.iter().zip(mu.iter()) {
task.put_c_j(*j, *mui)?;
}
let n_ones = vec![1.0; n as usize];
// budget constraint
{
let coni = task.get_num_con()?;
task.append_cons(1)?;
task.put_con_name(coni,"budget")?;
task.put_a_row(coni,
x,
n_ones.as_slice())?;
task.put_con_bound(coni,mosek::Boundkey::FX,w0,w0)?;
}
// |x-x0| <= z
{
let coni = task.get_num_con()?;
task.append_cons(2 * n)?;
for i in 0..n {
task.put_con_name(coni+i, format!("zabs1[{}]",1 + i).as_str())?;
task.put_con_name(coni+n+i, format!("zabs2[{}]",1 + i).as_str())?;
}
let ones = vec![1.0; n as usize];
let minusones = vec![-1.0; n as usize];
let con_abs1 : Vec<i32> = (coni..coni+n).collect();
let con_abs2 : Vec<i32> = (coni+n..coni+2*n).collect();
task.put_aij_list(con_abs1.as_slice(), x, minusones.as_slice())?;
task.put_aij_list(con_abs1.as_slice(), z, ones.as_slice())?;
task.put_con_bound_slice(coni,coni+n, vec![Boundkey::LO; n as usize].as_slice(), x0.iter().map(|&v| -v).collect::<Vec<f64>>().as_slice(), vec![INF; n as usize].as_slice())?;
task.put_aij_list(con_abs2.as_slice(), x, ones.as_slice())?;
task.put_aij_list(con_abs2.as_slice(), z, ones.as_slice())?;
task.put_con_bound_slice(coni+n,coni+n*2, vec![Boundkey::LO; n as usize].as_slice(), x0, vec![INF; n as usize].as_slice())?;
}
// cardinality constraint
{
let coni = task.get_num_con()?;
task.append_cons(1)?;
task.put_con_name(coni,"cardinality")?;
task.put_a_row(coni, y, n_ones.as_slice())?;
task.put_con_bound(coni,mosek::Boundkey::UP,p as f64,p as f64)?;
}
// (gamma,G'x) in Q
{
let afei = task.get_num_afe()?;
let acci = task.get_num_acc()?;
task.append_afes(k as i64+1)?;
let dom = task.append_quadratic_cone_domain(k as i64+1)?;
task.append_acc_seq(dom,
afei,
vec![0.0; k as usize + 1].as_slice())?;
task.put_acc_name(acci,"risk")?;
task.put_afe_g(afei,gamma)?;
for ((i,j),v) in iproduct!(0..n,0..n).zip(GT).filter(|(_,v)| **v != 0.0) {
task.put_afe_f_entry(afei + i as i64 + 1, j as i32, *v)?;
}
}
// Switch
{
let coni = task.get_num_con()?;
task.append_cons(n)?;
for i in 0..n {
task.put_con_name(coni + i, format!("switch[{}]",i+1).as_str())?;
}
let conlist : Vec<i32> = (coni..coni+n).collect();
task.put_aij_list(conlist.as_slice(), z, vec![1.0; n as usize].as_slice())?;
task.put_aij_list(conlist.as_slice(), y, vec![-w0; n as usize].as_slice())?;
task.put_con_bound_slice_const(coni,coni+n, Boundkey::UP, 0.0,0.0)?;
}
let _ = task.optimize()?;
task.write_data(format!("portfolio_5_card-{}.ptf",p).as_str())?;
// Check if the integer solution is an optimal point
if task.get_sol_sta(Soltype::ITG)? != Solsta::INTEGER_OPTIMAL {
// See https://docs.mosek.com/latest/rustapi/accessing-solution.html about handling solution statuses.
eprintln!("Solution not optimal!");
std::process::exit(1);
}
let mut xx = vec![0.0;n as usize];
task.get_xx_slice(Soltype::ITG, 0,n,xx.as_mut_slice())?;
Ok((xx[0..n as usize].to_vec(),task.get_primal_obj(Soltype::ITG)?))
}
#[allow(non_snake_case)]
fn main() -> Result<(),String> {
let n : i32 = 8;
let w = 1.0;
let mu = &[0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379];
let x0 = &[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
let 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 ];
let gamma = 0.25;
for K in 1..n+1 {
let (x,expret) = portfolio(n,
mu,
GT,
x0,
gamma,
K,
w)?;
println!("Bound {}: x = {:?}", K,x);
println!(" Return: {:.5e}\n", expret);
}
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
portfolio_6_factor.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : portfolio_6_factor.rs
//!
//! Purpose : Implements a portfolio optimization model using factor model.
//!
extern crate mosek;
extern crate itertools;
use mosek::{Task,Boundkey,Objsense,Streamtype,Soltype,Transpose,Solsta};
use itertools::{iproduct};
use std::convert::TryInto;
const INF : f64 = 0.0;
#[allow(non_snake_case)]
fn portfolio(w : f64,
mu : &[f64],
x0 : &[f64],
gammas : &[f64],
theta : &[f64],
GT : &Matrix) -> Result<Vec<(f64,f64)>,String> {
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
let (kx,nx) = GT.size();
if mu.len() != nx { panic!("Mismatching data"); }
let k : i32 = kx.try_into().unwrap();
let n : i32 = nx.try_into().unwrap();
let total_budget : f64 = w + x0.iter().sum::<f64>();
//Offset of variables into the API variable.
let numvar = n;
let voff_x : i32 = 0;
// Constraint offset
let coff_bud : i32 = 0;
// Holding variable x of length n
// No other auxiliary variables are needed in this formulation
task.append_vars(numvar)?;
// Setting up variable x
for j in 0..n {
task.put_var_name(voff_x+j,format!("x[{}]",j+1).as_str())?;
}
task.put_var_bound_slice_const(voff_x,voff_x+n, Boundkey::LO,0.0, INF)?;
// One linear constraint: total budget
task.append_cons(1)?;
task.put_con_name(coff_bud, "budget")?;
/* Coefficients in the first row of A */
for j in 0..n {
task.put_aij(coff_bud,voff_x + j, 1.0)?;
}
task.put_con_bound(coff_bud, Boundkey::FX, total_budget, total_budget)?;
// Input (gamma, G_factor_T x, diag(sqrt(theta))*x) in the AFE (affine expression) storage
// We need k+n+1 rows and we fill them in in three parts
task.append_afes((k+n) as i64 + 1)?;
// 1. The first affine expression = gamma, will be specified later
// 2. The next k expressions comprise G_factor_T*x, we add them row by row
// transposing the matrix G_factor on the fly
task.put_afe_f_row_list((1..1+k as i64).collect::<Vec<i64>>().as_slice(), // f row idxs
vec![n; k as usize].as_slice(), // row lengths
(0..GT.len() as i64).step_by(n as usize).collect::<Vec<i64>>().as_slice(), // row ptr
iproduct!(0..k,0..n).map(|(_,b)| b).collect::<Vec<i32>>().as_slice(), // varidx, 0..n repeated k times
GT.data_by_row().as_slice())?;
// 3. The remaining n rows contain sqrt(theta) on the diagonal
for (i,thetai) in (0..n).zip(theta.iter()) {
task.put_afe_f_entry(i as i64 + 1 + k as i64, voff_x + i, thetai.sqrt())?;
}
// Input the affine conic constraint (gamma, GT*x) \in QCone
// Add the quadratic domain of dimension k+1
let qdom = task.append_quadratic_cone_domain(k as i64+ 1)?;
// Add the constraint
task.append_acc_seq(qdom, 0, vec![0.0; k as usize+1].as_slice())?;
task.put_acc_name(0, "risk")?;
// Objective: maximize expected return mu^T x
for (j,&muj) in (0..n).zip(mu.iter()) {
task.put_c_j(voff_x + j, muj)?;
}
task.put_obj_sense(Objsense::MAXIMIZE)?;
Ok(gammas.iter().filter_map(|&gamma| {
// Specify gamma in ACC
task.put_afe_g(0, gamma).ok()?;
task.optimize().ok()?;
/* Display solution summary for quick inspection of results */
let _ = task.solution_summary(Streamtype::LOG);
let _ = task.write_data(format!("portfolio_6_factor-{}.ptf",gamma).as_str());
// Check if the interior point solution is an optimal point
if task.get_sol_sta(Soltype::ITR).ok()? != Solsta::OPTIMAL {
// See https://docs.mosek.com/latest/rustapi/accessing-solution.html about handling solution statuses.
eprintln!("Solution not optimal!");
std::process::exit(1);
}
/* Read the results */
let mut xx = vec![0.0; n as usize];
task.get_xx_slice(Soltype::ITR, voff_x, voff_x + n, xx.as_mut_slice()).ok()?;
Some((gamma,xx.iter().zip(mu.iter()).map(|(&xj,&muj)| xj*muj).sum::<f64>()))
}).collect::<Vec<(f64,f64)>>())
}
#[allow(non_snake_case)]
fn main() -> Result<(),String> {
// Since the value infinity is never used, we define
// 'infinity' for symbolic purposes only
let n : i32 = 8;
let w = 1.0;
let mu = &[0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379];
let x0 = &[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
// Factor exposure matrix, n x 2
let B = Matrix::new_by_row(n as usize, 2,
[ 0.4256, 0.1869,
0.2413, 0.3877,
0.2235, 0.3697,
0.1503, 0.4612,
1.5325, -0.2633,
1.2741, -0.2613,
0.6939, 0.2372,
0.5425, 0.2116 ].to_vec()).unwrap();
// Factor covariance matrix, 2x2
let S_F = Matrix::new_by_row(2,2,
[ 0.0620, 0.0577,
0.0577, 0.0908 ].to_vec()).unwrap();
// Specific risk components
let theta : &[f64] = &[0.0720, 0.0508, 0.0377, 0.0394, 0.0663, 0.0224, 0.0417, 0.0459];
let S_sqrt_theta = Matrix::diag_matrix(theta.iter().map(|&v| v.sqrt()).collect());
let P = S_F.cholesky().unwrap();
let BP = B.mul(&P).unwrap();
//let GT = BP.concat_h(&S_theta.sqrt_element().unwrap()).unwrap().transpose();
let GT = BP.concat_h(&S_sqrt_theta).unwrap().transpose();
let gammas = &[0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48];
for (gamma,expret) in portfolio(w,
mu,
x0,
gammas,
theta,
>)? {
println!("Expected return {:.4e} for gamma {:.2}", expret, gamma);
}
Ok(())
}
// Matrix with data stored in colunn format
#[derive(Copy,Clone)]
pub enum MatrixOrder {
ByRow,
ByCol
}
pub struct Matrix {
fmt : MatrixOrder,
dimi : usize,
dimj : usize,
data : Vec<f64>
}
impl Matrix {
pub fn new(fmt : MatrixOrder, dimi : usize, dimj : usize, data : Vec<f64>) -> Option<Matrix> {
if dimi*dimj == data.len() {
Some(Matrix{fmt,dimi,dimj,data})
}
else {
None
}
}
pub fn new_by_col(dimi : usize, dimj : usize, data : Vec<f64>) -> Option<Matrix> {
Matrix::new(MatrixOrder::ByCol,dimi,dimj,data)
}
pub fn new_by_row(dimi : usize, dimj : usize, data : Vec<f64>) -> Option<Matrix> {
Matrix::new(MatrixOrder::ByRow,dimi,dimj,data)
}
pub fn size(&self) -> (usize,usize) { (self.dimi,self.dimj) }
pub fn len(&self) -> usize { self.data.len() }
pub fn diag_matrix(data : Vec<f64>) -> Matrix {
let mut rdata = vec![0.0; data.len()*data.len()];
for (r,&v) in rdata.iter_mut().step_by(data.len()+1).zip(data.iter()) { *r = v; }
Matrix{fmt:MatrixOrder::ByCol,dimi:data.len(),dimj:data.len(),data:rdata}
}
pub fn transpose(&self) -> Matrix {
Matrix {
fmt : match self.fmt { MatrixOrder::ByRow => MatrixOrder::ByCol, MatrixOrder::ByCol => MatrixOrder::ByRow },
dimi : self.dimj,
dimj : self.dimi,
data : self.data.as_slice().to_vec()
}
}
pub fn data_by_row(&self) -> Vec<f64> {
match self.fmt {
MatrixOrder::ByRow => self.data.clone(),
MatrixOrder::ByCol => iproduct!(0..self.dimi,0..self.dimj).map(|(i,j)| unsafe { *self.data.get_unchecked(self.dimi*j+i) }).collect()
}
}
pub fn data_by_col(&self) -> Vec<f64> {
match self.fmt {
MatrixOrder::ByCol => self.data.clone(),
MatrixOrder::ByRow => iproduct!(0..self.dimj,0..self.dimi).map(|(j,i)| unsafe { *self.data.get_unchecked(self.dimj*i+j) }).collect()
}
}
pub fn data_by_row_to(&self,r : &mut[f64]) {
match self.fmt {
MatrixOrder::ByRow => { r.clone_from_slice(self.data.as_slice()); },
MatrixOrder::ByCol => {
if r.len() != self.data.len() { panic!("Incorrect result array length"); }
for (r,(i,j)) in r.iter_mut().zip(iproduct!(0..self.dimi,0..self.dimj)) {
*r = unsafe { *self.data.get_unchecked(self.dimi*j+i) }
}
}
}
}
pub fn data_by_col_to(&self,r : &mut[f64]) {
match self.fmt {
MatrixOrder::ByCol => { r.clone_from_slice(self.data.as_slice()); },
MatrixOrder::ByRow => {
for (r,(j,i)) in r.iter_mut().zip(iproduct!(0..self.dimj,0..self.dimi)) {
*r = unsafe { *self.data.get_unchecked(self.dimj*i+j) };
}
}
}
}
pub fn cholesky(&self) -> Option<Matrix> {
if self.dimi != self.dimj { None }
else {
let mut resdata = self.data_by_col();
// zero data in the non-tril part
for ((j,i),d) in iproduct!(0..self.dimj,0..self.dimi).zip(resdata.iter_mut()) { if i < j {*d = 0.0}; }
if let Ok(_) = mosek::potrf(mosek::Uplo::LO,
self.dimi.try_into().unwrap(),
resdata.as_mut_slice()) {
Some(Matrix{fmt:MatrixOrder::ByCol,dimi:self.dimi,dimj:self.dimj,data:resdata})
}
else { None }
}
}
pub fn mul(&self,other : &Matrix) -> Option<Matrix> {
if self.dimj != other.dimi { None }
else {
let mut resdata = vec![0.0; self.dimi * other.dimj];
if let Ok(_) = mosek::gemm(match self.fmt { MatrixOrder::ByCol => Transpose::NO, MatrixOrder::ByRow => Transpose::YES },
match other.fmt { MatrixOrder::ByCol => Transpose::NO, MatrixOrder::ByRow => Transpose::YES },
self.dimi.try_into().unwrap(),
other.dimj.try_into().unwrap(),
self.dimj.try_into().unwrap(),
1.0,
self.data.as_slice(),
other.data.as_slice(),
1.0,
resdata.as_mut_slice()) {
Matrix::new_by_col(self.dimi,other.dimj,resdata)
}
else {
None
}
}
}
pub fn concat_h(&self, other : & Matrix) -> Option<Matrix> {
if self.dimi != other.dimi {
None
}
else {
let mut resdata = vec![0.0; self.data.len() + other.data.len()];
self.data_by_col_to(& mut resdata[..self.len()]);
other.data_by_col_to(& mut resdata[self.len()..]);
Matrix::new_by_col(self.dimi,self.dimj+other.dimj,resdata)
}
}
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
pow1.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : pow1.rs
//!
//! 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
//!
extern crate mosek;
use mosek::*;
const INF : f64 = 0.0;
fn main() -> Result<(),String> {
let numcon : i32 = 1;
let numvar : i32 = 5;
// Since the value infinity is never used, we define
// 'infinity' symbolic purposes only
let cval = vec![ 1.0, 1.0, -1.0 ];
let csub = vec![ 3, 4, 0 ];
let aval = vec![ 1.0, 1.0, 0.5 ];
let asub = vec![ 0, 1, 2 ];
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
// Directs the log task stream to the user specified
// method msgclass.streamCB
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* Append 'numcon' empty constraints.
The constraints will initially have no bounds. */
task.append_cons(numcon)?;
/* Append 'numvar' variables.
The variables will initially be fixed at zero (x=0). */
task.append_vars(numvar)?;
/* Set up the linear part of the problem */
task.put_c_list(&csub, &cval)?;
task.put_a_row(0, &asub, &aval)?;
task.put_con_bound(0, Boundkey::FX, 2.0, 2.0)?;
task.put_var_bound_slice_const(0, numvar, Boundkey::FR, -INF, INF)?;
/* Add a conic constraint */
let pc1 = task.append_primal_power_cone_domain(3, &[0.2, 0.8])?;
let pc2 = task.append_primal_power_cone_domain(3, &[4.0, 6.0])?;
// Create data structures F,g so that
//
// F * x + g = (x(0), x(1), x(3), x(2), 1.0, x(4))
//
task.append_afes(6)?;
task.put_afe_f_entry_list(&[0, 1, 2, 3, 5], // Rows
&[0, 1, 3, 2, 4], // Columns
&[1.0, 1.0, 1.0, 1.0, 1.0])?;
task.put_afe_g(4, 1.0)?;
// Append the two conic constraints
task.append_acc(pc1, // Domain
&[0, 1, 2], // Rows from F
&[0.0,0.0,0.0])?; // Unused
task.append_acc(pc2, // Domain
&[3, 4, 5], // Rows from F
&[0.0,0.0,0.0])?; // Unused
task.put_obj_sense(Objsense::MAXIMIZE)?;
task.optimize()?;
task.write_data("pow1.ptf")?;
// Print a summary containing information
// about the solution for debugging purposes
task.solution_summary(Streamtype::LOG)?;
/* Get status information about the solution */
let solsta = task.get_sol_sta(Soltype::ITR)?;
assert!(solsta == Solsta::OPTIMAL);
let mut xx = vec![0.0; numvar as usize];
task.get_xx(Soltype::ITR,
xx.as_mut_slice())?;
println!("Optimal primal solution");
for (j,&xj) in xx[0..3].iter().enumerate() {
println!("x[{}]: {}",j+1,xj);
}
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
qcqo1.rs
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : qcqo1.rs
//!
//! Purpose : Demonstrate how to solve a quadratic
//! optimization problem using the MOSEK API.
//!
//! minimize x0^2 + 0.1 x1^2 + x2^2 - x0 x2 - x1
//! s.t
//! 1 <= x0 + x1 + x2 - x0^2 - x1^2 - 0.1 x2^2 + 0.2 x0 x2
//! x >= 0
extern crate mosek;
extern crate itertools;
use mosek::{Task,Boundkey,Objsense,Streamtype,Solsta,Soltype};
use itertools::izip;
const INF : f64 = 0.0;
fn main() -> Result<(),String> {
const NUMCON : i32 = 1; /* Number of constraints. */
const NUMVAR : i32 = 3; /* Number of variables. */
let c = [0.0, -1.0, 0.0];
let bkc = [Boundkey::LO];
let blc = [1.0];
let buc = [INF];
let bkx = [Boundkey::LO,
Boundkey::LO,
Boundkey::LO ];
let blx = [0.0,
0.0,
0.0 ];
let bux = [INF,
INF,
INF ];
let asub = [ &[0i32], &[0i32], &[0i32] ];
let aval = [ &[1.0], &[1.0], &[1.0] ];
let mut task = Task::new().unwrap().with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
// Give MOSEK an estimate of the size of the input data.
// This is done to increase the speed of inputting data.
// However, it is optional.
// Append 'numcon' empty constraints.
// The constraints will initially have no bounds.
task.append_cons(NUMCON)?;
// Append 'numvar' variables.
// The variables will initially be fixed at zero (x=0).
task.append_vars(NUMVAR)?;
for (j,cj,bkj,blj,buj) in izip!(0..NUMVAR,c,bkx,blx,bux) {
// Set the linear term c_j in the objective.
task.put_c_j(j, cj)?;
// Set the bounds on variable j.
// blx[j] <= x_j <= bux[j]
task.put_var_bound(j, bkj, blj, buj)?;
}
for (j,asubj,avalj) in izip!(0..NUMVAR, asub, aval) {
/* Input column j of A */
task.put_a_col(j, /* Variable (column) index.*/
asubj, /* Row index of non-zeros in column j.*/
avalj)?; /* Non-zero Values of column j. */
}
// Set the bounds on constraints.
// for i=1, ...,numcon : blc[i] <= constraint i <= buc[i]
for (i,bki,bli,bui) in izip!(0..NUMCON,bkc,blc,buc) {
task.put_con_bound(i, bki, bli, bui)?;
}
{
// The lower triangular part of the Q
// matrix in the objective is specified.
let qosubi = &[ 0, 1, 2, 2 ];
let qosubj = &[ 0, 1, 0, 2 ];
let qoval = &[ 2.0, 0.2, -1.0, 2.0 ];
// Input the Q for the objective.
task.put_q_obj(qosubi, qosubj, qoval)?;
}
// The lower triangular part of the Q^0
// matrix in the first constraint is specified.
// This corresponds to adding the term
// x0^2 - x1^2 - 0.1 x2^2 + 0.2 x0 x2
{
let qsubi = &[0, 1, 2, 2 ];
let qsubj = &[0, 1, 2, 0 ];
let qval = &[-2.0, -2.0, -0.2, 0.2];
/* put Q^0 in constraint with index 0. */
task.put_q_con_k(0,
qsubi,
qsubj,
qval)?;
}
task.put_obj_sense(Objsense::MINIMIZE)?;
/* Solve the problem */
let _trm = task.optimize()?;
// Print a summary containing information
// about the solution for debugging purposes
task.solution_summary(Streamtype::MSG)?;
/* Get status information about the solution */
let solsta = task.get_sol_sta(Soltype::ITR)?;
let mut xx = vec![0.0; NUMVAR as usize];
task.get_xx(Soltype::ITR, // Interior solution.
xx.as_mut_slice())?;
match solsta {
Solsta::OPTIMAL => {
println!("Optimal primal solution");
for (j,xj) in izip!(0..NUMVAR,xx) {
println!("x[{}]: {}",j,xj);
}
},
Solsta::DUAL_INFEAS_CER|
Solsta::PRIM_INFEAS_CER =>
println!("Primal or dual infeasibility.\n"),
Solsta::UNKNOWN =>
println!("Unknown solution status.\n"),
_ =>
println!("Other solution status")
}
Ok(())
} /* Main */
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
qo1.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : qo1.rs
//!
//! Purpose: To demonstrate how to solve a quadratic optimization
//! problem using the MOSEK API.
//!
extern crate mosek;
use mosek::{Task,Boundkey,Streamtype,Solsta,Soltype};
const INF : f64 = 0.0;
const NUMCON : usize = 1; /* Number of constraints. */
const NUMVAR : usize = 3; /* Number of variables. */
fn main() -> Result<(),String> {
let c = vec![ 0.0,-1.0,0.0 ];
let bkc = vec![ mosek::Boundkey::LO ];
let blc = vec![ 1.0 ];
let buc = vec![ INF ];
let bkx = vec![ Boundkey::LO,
Boundkey::LO,
Boundkey::LO ];
let blx = vec![ 0.0,
0.0,
0.0 ];
let bux = vec![ INF,
INF,
INF ];
let aptrb = vec![ 0, 1, 2 ];
let aptre = vec![ 1, 2, 3 ];
let asub = vec![ 0, 0, 0 ];
let aval = vec![ 1.0, 1.0, 1.0 ];
let qsubi = vec![ 0, 1, 2, 2 ];
let qsubj = vec![ 0, 1, 0, 2 ];
let qval = vec![ 2.0,0.2,-1.0,2.0 ];
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
//r = MSK_linkfunctotaskstream(task,MSK_STREAM_LOG,NULL,printstr);
task.append_cons(NUMCON as i32)?;
/* Append 'NUMVAR' variables.
* The variables will initially be fixed at zero (x=0). */
task.append_vars(NUMVAR as i32)?;
/* Optionally add a constant term to the objective. */
task.put_cfix(0.0)?;
for j in 0..NUMVAR
{
/* Set the linear term c_j in the objective.*/
task.put_c_j(j as i32,c[j])?;
/* Set the bounds on variable j.
* blx[j] <= x_j <= bux[j] */
task.put_var_bound(j as i32, /* Index of variable.*/
bkx[j], /* Bound key.*/
blx[j], /* Numerical value of lower bound.*/
bux[j])?; /* Numerical value of upper bound.*/
/* Input column j of A */
task.put_a_col(j as i32, /* Variable (column) index.*/
&asub[aptrb[j]..aptre[j]], /* Pointer to row indexes of column j.*/
&aval[aptrb[j]..aptre[j]])?; /* Pointer to Values of column j.*/
}
/* Set the bounds on constraints.
* for i=1, ...,NUMCON : blc[i] <= constraint i <= buc[i] */
for i in 0..NUMCON
{
task.put_con_bound(i as i32, /* Index of constraint.*/
bkc[i], /* Bound key.*/
blc[i], /* Numerical value of lower bound.*/
buc[i])?; /* Numerical value of upper bound.*/
/*
* The lower triangular part of the Q
* matrix in the objective is specified.
*/
/* Input the Q for the objective. */
task.put_q_obj(&qsubi,&qsubj,&qval)?;
}
let _trmcode = task.optimize()?;
/* Run optimizer */
/* Print a summary containing information
about the solution for debugging purposes*/
task.solution_summary(Streamtype::MSG)?;
let solsta = task.get_sol_sta(Soltype::ITR)?;
match solsta
{
Solsta::OPTIMAL =>
{
let mut xx = vec![0.0, 0.0, 0.0];
task.get_xx(Soltype::ITR, /* Request the interior solution. */
& mut xx[..])?;
println!("Optimal primal solution");
for j in 0..NUMVAR
{
println!("x[{}]: {}",j,xx[j]);
}
}
Solsta::DUAL_INFEAS_CER |
Solsta::PRIM_INFEAS_CER =>
{
println!("Primal or dual infeasibility certificate found.");
}
Solsta::UNKNOWN =>
{
println!("The status of the solution could not be determined.");
}
_ =>
{
println!("Other solution status.");
}
}
return Ok(());
} /* main */
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
reoptimization.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! Purpose : Demonstrates how to solve a linear
//! optimization problem using the MOSEK API
//! and modify and re-optimize the problem.
extern crate mosek;
extern crate itertools;
use mosek::{Task,Boundkey,Objsense,Soltype};
use itertools::izip;
const INF : f64 = 0.0;
fn main() -> Result<(),String> {
let numcon = 3;
let numvar = 3;
let c = &[1.5, 2.5, 3.0 ];
let bkc = &[ Boundkey::UP,
Boundkey::UP,
Boundkey::UP ];
let blc = &[ -INF,
-INF,
-INF ];
let buc = &[ 100000.0,
50000.0,
60000.0 ];
let bkx = &[ Boundkey::LO,
Boundkey::LO,
Boundkey::LO
];
let blx = &[ 0.0, 0.0, 0.0 ];
let bux = &[ INF,
INF,
INF ];
let asub = &[
&[ 0i32, 1, 2 ],
&[ 0i32, 1, 2 ],
&[ 0i32, 1, 2 ] ];
let aval = &[
&[ 2.0, 3.0, 2.0 ],
&[ 4.0, 2.0, 3.0 ],
&[ 3.0, 3.0, 2.0 ] ];
let mut task = Task::new().unwrap();
/* Append the constraints. */
task.append_cons(numcon)?;
/* Append the variables. */
task.append_vars(numvar)?;
/* Put C. */
for (j,&cj) in (0..numvar).zip(c.iter()) {
task.put_c_j(j,cj)?;
}
/* Put constraint bounds. */
for (i,&bki,&bli,&bui) in izip!(0..numcon,bkc,blc,buc) {
task.put_con_bound(i, bki, bli, bui)?;
}
/* Put variable bounds. */
for (j,&bki,&bli,&bui) in izip!(0..numvar,bkx,blx,bux) {
task.put_var_bound(j, bki, bli, bui)?;
}
/* Put A. */
if numcon > 0 {
for (j,&asubj,&avalj) in izip!(0..numvar,asub,aval) {
task.put_a_col(j,
asubj,
avalj)?;
}
}
/* A maximization problem */
task.put_obj_sense(Objsense::MAXIMIZE)?;
/* Solve the problem */
let _trm = task.optimize()?;
let mut xx = vec![0.0; task.get_num_var()? as usize];
task.get_xx(Soltype::BAS, // Request the basic solution.
xx.as_mut_slice())?;
for (j,xj) in xx.iter().enumerate() {
println!("x[{}]: {}",j,xj);
}
/****************** Make a change to the A matrix ******************/
task.put_aij(0, 0, 3.0)?;
let _trm = task.optimize()?;
task.get_xx(Soltype::BAS, // Request the basic solution.
xx.as_mut_slice())?;
for (j,xj) in xx.iter().enumerate() {
println!("x[{}]: {}",j,xj);
}
/***************** Add a new variable ******************************/
/* Get index of new variable. */
let varidx = task.get_num_var()?;
/* Append a new variable x_3 to the problem */
task.append_vars(1)?;
let numvar = numvar + 1;
/* Set bounds on new varaible */
task.put_var_bound(varidx, Boundkey::LO, 0.0, INF)?;
/* Change objective */
task.put_c_j(varidx, 1.0)?;
/* Put new values in the A matrix */
let acolsub = &[0i32, 2];
let acolval = &[4.0, 1.0];
task.put_a_col(varidx, /* column index */
acolsub,
acolval)?;
/* Change optimizer to simplex free and reoptimize */
task.put_int_param(mosek::Iparam::OPTIMIZER, mosek::Optimizertype::FREE_SIMPLEX)?;
let _trm = task.optimize()?;
let mut xx = vec![0.0; task.get_num_var()? as usize];
task.get_xx(Soltype::BAS, xx.as_mut_slice())?;
for (j,xj) in (0..numvar).zip(xx.iter()) {
println!("x[{}]: {}",j,xj);
}
/********************** Add a new constraint ***************************/
/* Get index of new constraint. */
let conidx = task.get_num_con()?;
/* Append a new constraint */
task.append_cons(1)?;
let numcon = numcon + 1;
/* Set bounds on new constraint */
task.put_con_bound(conidx,
Boundkey::UP,
-INF,
30000.0)?;
/* Put new values in the A matrix */
let arowsub = &[0i32, 1, 2, 3 ];
let arowval = &[1.0, 2.0, 1.0, 1.0 ];
task.put_a_row(conidx, /* row index */
arowsub,
arowval)?;
let _trm = task.optimize()?;
task.get_xx(Soltype::BAS, // Request the basic solution.
xx.as_mut_slice())?;
for (j,xj) in (0..numvar).zip(xx.iter()) {
println!("x[{}]: {}",j,xj);
}
/********************** Change constraint bounds ********************/
let newbkc = &[Boundkey::UP,
Boundkey::UP,
Boundkey::UP,
Boundkey::UP];
let newblc = &[-INF,
-INF,
-INF,
-INF];
let newbuc = &[ 80000.0, 40000.0, 50000.0, 22000.0 ];
task.put_con_bound_slice(0, numcon, newbkc, newblc, newbuc)?;
let _ = task.optimize()?;
task.get_xx(Soltype::BAS, // Request the basic solution.
xx.as_mut_slice())?;
for (j,xj) in (0..numvar).zip(xx.iter()) {
println!("x[{}]: {}",j,xj);
}
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
response.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : response.rs
//!
//! Purpose : This example demonstrates proper response handling
//! for problems solved with the interior-point optimizers.
//!
extern crate mosek;
use mosek::{Task,Streamtype,Solsta,Soltype};
use std::env;
const CQO1_PTF : &str = "Task 'CQO1 EXAMPLE'
Objective obj
Minimize + x4 + x5 + x6
Constraints
c1 [1] + x1 + x2 + 2 x3
Variables
k1 [QUAD(3)]
x4
x1 [0;+inf]
x2 [0;+inf]
k2 [RQUAD(3)]
x5
x6
x3 [0;+inf]
";
fn main() -> Result<(),String> {
let args: Vec<String> = env::args().collect();
let mut task = Task::new().unwrap().with_callbacks();
if args.len() < 2 {
task.read_ptf_string(CQO1_PTF)?;
}
else {
task.read_data(args[1].as_str())?;
}
// Perform optimization.
let trm = task.optimize()?;
task.solution_summary(Streamtype::LOG)?;
// Handle solution status. We expect Optimal
let solsta = task.get_sol_sta(Soltype::ITR)?;
match solsta {
Solsta::OPTIMAL => {
// Fetch and print the solution
println!("An optimal interior point solution is located.");
let numvar = task.get_num_var()?;
let mut xx = vec![0.0; numvar as usize];
task.get_xx(Soltype::ITR, xx.as_mut_slice())?;
println!("xx = {:?}",xx)
},
Solsta::DUAL_INFEAS_CER =>
println!("Dual infeasibility certificate found."),
Solsta::PRIM_INFEAS_CER =>
println!("Primal infeasibility certificate found."),
Solsta::UNKNOWN => {
// The solutions status is unknown. The termination code
// indicates why the optimizer terminated prematurely.
println!("The solution status is unknown.");
let (symname,desc) = mosek::get_code_desc(trm)?;
println!(" Termination code: {} {}\n", symname, desc)
},
_ =>
println!("Unexpected solution status {}\n",solsta)
}
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
sdo1.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : sdo1.rs
//!
//! Purpose: Solves the following small semidefinite optimization problem
//! using the MOSEK API.
//!
//! minimize Tr [2, 1, 0; 1, 2, 1; 0, 1, 2]*X + x0
//!
//! subject to Tr [1, 0, 0; 0, 1, 0; 0, 0, 1]*X + x0 = 1
//! Tr [1, 1, 1; 1, 1, 1; 1, 1, 1]*X + x1 + x2 = 0.5
//! (x0,x1,x2) \in Q, X \in PSD
//!
extern crate mosek;
use mosek::{Task,Streamtype,Solsta,Soltype};
const INF : f64 = 0.0;
const NUMCON : usize = 2; /* Number of constraints. */
const NUMVAR : usize = 3; /* Number of conic quadratic variables */
fn main() -> Result<(),String>
{
let dimbarvar = vec![3]; /* Dimension of semidefinite cone */
let bkc = &[ mosek::Boundkey::FX, mosek::Boundkey::FX ];
let blc = &[ 1.0, 0.5 ];
let buc = &[ 1.0, 0.5 ];
let barc_i = &[0, 1, 1, 2, 2];
let barc_j = &[0, 0, 1, 1, 2];
let barc_v = &[2.0, 1.0, 2.0, 1.0, 2.0];
let aptrb = &[0, 1];
let aptre = &[1, 3];
let asub = &[0, 1, 2]; /* column subscripts of A */
let aval = &[1.0, 1.0, 1.0];
let bara_i = &[0, 1, 2, 0, 1, 2, 1, 2, 2];
let bara_j = &[0, 1, 2, 0, 0, 0, 1, 1, 2];
let bara_v = &[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
let falpha = 1.0;
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* Append 'NUMCON' empty constraints.
* The constraints will initially have no bounds. */
task.append_cons(NUMCON as i32)?;
/* Append 'NUMVAR' variables.
* The variables will initially be fixed at zero (x=0). */
task.append_vars(NUMVAR as i32)?;
/* Append 'NUMBARVAR' semidefinite variables. */
task.append_barvars(&dimbarvar[..])?;
/* Optionally add a constant term to the objective. */
task.put_cfix(0.0)?;
/* Set the linear term c_j in the objective.*/
task.put_c_j(0,1.0)?;
for j in 0..NUMVAR {
task.put_var_bound(j as i32,
mosek::Boundkey::FR,
-INF,
INF)?;
}
/* Set the linear term barc_j in the objective.*/
let c_symmat_idx = task.append_sparse_sym_mat(dimbarvar[0],
barc_i,
barc_j,
barc_v)?;
task.put_barc_j(0, &[c_symmat_idx], &[falpha])?;
for i in 0..NUMCON
{
/* Input A row by row */
task.put_a_row(i as i32,
& asub[aptrb[i]..aptre[i]],
& aval[aptrb[i]..aptre[i]])?;
/* Set the bounds on constraints.
* for i=1, ...,NUMCON : blc[i] <= constraint i <= buc[i] */
task.put_con_bound(i as i32, /* Index of constraint.*/
bkc[i], /* Bound key.*/
blc[i], /* Numerical value of lower bound.*/
buc[i])?; /* Numerical value of upper bound.*/
}
{
/* Append the conic quadratic cone */
let afei = task.get_num_afe()?;
task.append_afes(3)?;
task.put_afe_f_entry_list(&[0,1,2],
&[0,1,2],
&[1.0,1.0,1.0])?;
let dom = task.append_quadratic_cone_domain(3)?;
task.append_acc_seq(dom,afei,&[0.0,0.0,0.0])?;
}
/* Add the first row of barA */
let a_symmat_idx1 =
task.append_sparse_sym_mat(dimbarvar[0],
& bara_i[..3],
& bara_j[..3],
& bara_v[..3])?;
task.put_bara_ij(0, 0, &[a_symmat_idx1][..], &[falpha][..])?;
/* Add the second row of barA */
let a_symmat_idx2 =
task.append_sparse_sym_mat(dimbarvar[0],
& bara_i[3..9],
& bara_j[3..9],
& bara_v[3..9])?;
task.put_bara_ij(1, 0, &[a_symmat_idx2][..], &[falpha][..])?;
let _trmcode = task.optimize()?;
task.write_data("sdo1.ptf")?;
/* Print a summary containing information
* about the solution for debugging purposes*/
task.solution_summary (Streamtype::MSG)?;
let solsta = task.get_sol_sta(Soltype::ITR)?;
match solsta
{
Solsta::OPTIMAL =>
{
let mut xx = vec![0.0,0.0,0.0];
task.get_xx(Soltype::ITR, /* Request the basic solution. */
& mut xx[..])?;
let mut barx = vec![0.0,0.0,0.0,0.0,0.0,0.0];
task.get_barx_j(Soltype::ITR, /* Request the interior solution. */
0,
& mut barx[..])?;
println!("Optimal primal solution");
for j in 0..NUMVAR as usize
{
println!("x[{}]: {}",j,xx[j]);
}
let n = dimbarvar[0] as usize;
for j in 0..n
{
for i in j..n
{
println!("barx[{},{}]: {}",i,j,barx[j*n+i-j*(j+1)/2]);
}
}
}
Solsta::DUAL_INFEAS_CER |
Solsta::PRIM_INFEAS_CER =>
{
println!("Primal or dual infeasibility certificate found.");
}
Solsta::UNKNOWN =>
{
/* If the solutions status is unknown, print the termination code
* indicating why the optimizer terminated prematurely. */
println!("The solution status is unknown.");
println!("The optimizer terminitated with code: {}",solsta);
}
_ =>
{
println!("Other solution status.");
}
}
Ok(())
} /* main */
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
sdo2.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : sdo2.rs
//!
//! 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.
//!
extern crate mosek;
use mosek::{Task,Streamtype,Solsta,Soltype};
#[allow(non_snake_case)]
fn main() -> Result<(),String> {
/* Input data */
let numcon : i32 = 2; /* Number of constraints. */
let dimbarvar : &[i32] = &[3, 4]; /* Dimension of semidefinite variables */
/* Objective coefficients concatenated */
let Cj : &[i32] = &[ 0, 0, 1, 1, 1, 1 ]; /* Which symmetric variable (j) */
let Ck : &[i32] = &[ 0, 2, 0, 1, 1, 2 ]; /* Which entry (k,l)->v */
let Cl : &[i32] = &[ 0, 2, 0, 0, 1, 2 ];
let Cv : &[f64] = &[ 1.0, 6.0, 1.0, -3.0, 2.0, 1.0 ];
/* Equality constraints coefficients concatenated */
let Ai : &[i32] = &[ 0, 0, 0, 0, 0, 0 ]; /* Which constraint (i = 0) */
let Aj : &[i32] = &[ 0, 0, 0, 1, 1, 1 ]; /* Which symmetric variable (j) */
let Ak : &[i32] = &[ 0, 2, 2, 1, 1, 3 ]; /* Which entry (k,l)->v */
let Al : &[i32] = &[ 0, 0, 2, 0, 1, 3 ];
let Av : &[f64] = &[ 1.0, 1.0, 2.0, 1.0, -1.0, -3.0 ];
/* The second constraint - one-term inequality */
let A2i : &[i32] = &[ 1 ]; /* Which constraint (i = 1) */
let A2j : &[i32] = &[ 1 ]; /* Which symmetric variable (j = 1) */
let A2k : &[i32] = &[ 1 ]; /* Which entry A(1,0) = A(0,1) = 0.5 */
let A2l : &[i32] = &[ 0 ];
let A2v : &[f64] = &[ 0.5 ];
let bkc = &[ mosek::Boundkey::FX,
mosek::Boundkey::UP ];
let blc = &[ 23.0, 0.0 ];
let buc = &[ 23.0, -3.0 ];
/* Create the optimization task. */
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* Append numcon empty constraints.
The constraints will initially have no bounds. */
task.append_cons(numcon)?;
/* Append numbarvar semidefinite variables. */
task.append_barvars(dimbarvar)?;
/* Set objective (6 nonzeros).*/
task.put_barc_block_triplet(Cj, Ck, Cl, Cv)?;
/* Set the equality constraint (6 nonzeros).*/
task.put_bara_block_triplet(Ai, Aj, Ak, Al, Av)?;
/* Set the inequality constraint (1 nonzero).*/
task.put_bara_block_triplet(A2i, A2j, A2k, A2l, A2v)?;
/* Set constraint bounds */
task.put_con_bound_slice(0, 2, bkc, blc, buc)?;
/* Run optimizer */
task.optimize()?;
task.solution_summary(Streamtype::MSG)?;
//mosek.solsta[] solsta = new mosek.solsta[1];
let solsta = task.get_sol_sta (Soltype::ITR)?;
match solsta {
Solsta::OPTIMAL => {
/* Retrieve the soution for all symmetric variables */
println!("Solution (lower triangular part vectorized):");
for (i,dimbarvari) in dimbarvar.iter().enumerate() {
//let dim = dimbarvar[i] * (dimbarvar[i] + 1) / 2;
let dim = dimbarvari * (dimbarvari+1)/2;
//double[] barx = new double[dim];
let mut barx : Vec<f64> = vec![0.0; dim as usize];
task.get_barx_j(Soltype::ITR, i as i32, barx.as_mut_slice())?;
println!("X{}: {:?}",i+1,barx);
// for (int j = 0; j < dim; ++j)
// System.out.print(barx[j] + " ");
// System.out.println();
}
},
Solsta::DUAL_INFEAS_CER|Solsta::PRIM_INFEAS_CER =>
println!("Primal or dual infeasibility certificate found."),
Solsta::UNKNOWN =>
println!("The status of the solution could not be determined."),
_ => println!("Other solution status.")
}
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
sdo_lmi.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : sdo_lmi.rs
//!
//! Purpose : To solve a problem with an LMI and an affine conic constrained problem with a PSD term
//!
//! minimize Tr [1, 0; 0, 1]*X + x(1) + x(2) + 1
//!
//! subject to Tr [0, 1; 1, 0]*X - x(1) - x(2) >= 0
//! x(1) [0, 1; 1, 3] + x(2) [3, 1; 1, 0] - [1, 0; 0, 1] >> 0
//! X >> 0
//!
extern crate mosek;
use mosek::{Task,Boundkey,Objsense,Streamtype,Solsta,Soltype};
const INF : f64 = 0.0;
fn main() -> Result<(),String> {
let numafe : i64 = 4; /* Number of affine expressions. */
let numvar : i32 = 2; /* Number of scalar variables */
let dimbarvar = &[2]; /* Dimension of semidefinite cone */
let lenbarvar = &[2 * (2 + 1) / 2]; /* Number of scalar SD variables */
let barc_j = &[0, 0];
let barc_k = &[0, 1];
let barc_l = &[0, 1];
let barc_v = &[1.0, 1.0];
let afeidx = &[0, 0, 1, 2, 2, 3];
let varidx = &[0, 1, 1, 0, 1, 0];
let f_val = &[-1.0, -1.0, 3.0, 2.0f64.sqrt(), 2.0f64.sqrt(), 3.0];
let g = &[0.0, -1.0, 0.0, -1.0];
let barf_i = &[0, 0];
let barf_j = &[0, 0];
let barf_k = &[0, 1];
let barf_l = &[0, 0];
let barf_v = &[0.0, 1.0];
let mut task = Task::new().unwrap().with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* Append 'NUMAFE' empty affine expressions. */
task.append_afes(numafe)?;
/* Append 'NUMVAR' variables.
The variables will initially be fixed at zero (x=0). */
task.append_vars(numvar)?;
/* Append 'NUMBARVAR' semidefinite variables. */
task.append_barvars(dimbarvar)?;
task.put_obj_sense(Objsense::MINIMIZE)?;
/* Optionally add a constant term to the objective. */
task.put_cfix(1.0)?;
/* Set the linear term c_j in the objective.*/
task.put_c_j(0, 1.0)?;
task.put_c_j(1, 1.0)?;
task.put_var_bound_slice_const(0,numvar, Boundkey::FR, -INF,INF)?;
/* Set the linear term barc_j in the objective.*/
task.put_barc_block_triplet(barc_j, barc_k, barc_l, barc_v)?;
/* Set up the affine conic constraints */
/* Construct the affine expressions */
/* F matrix */
task.put_afe_f_entry_list(afeidx, varidx, f_val)?;
/* g vector */
task.put_afe_g_slice(0, 4, g)?;
/* barF block triplets */
task.put_afe_barf_block_triplet(barf_i, barf_j, barf_k, barf_l, barf_v)?;
/* Append R+ domain and the corresponding ACC */
{
let dom = task.append_rplus_domain(1)?;
task.append_acc(dom, &[0], &[0.0])?;
}
/* Append SVEC_PSD domain and the corresponding ACC */
{
let dom = task.append_svec_psd_cone_domain(3)?;
task.append_acc(dom, &[1,2,3], &[0.0,0.0,0.0])?;
}
/* Run optimizer */
let _ = task.optimize()?;
/* Print a summary containing information
about the solution for debugging purposes*/
task.solution_summary (mosek::Streamtype::MSG)?;
let solsta = task.get_sol_sta(mosek::Soltype::ITR)?;
match solsta {
Solsta::OPTIMAL => {
let mut xx = vec![0.0; numvar as usize];
task.get_xx(Soltype::ITR,xx.as_mut_slice())?;
let mut barx = vec![0.0; lenbarvar[0]];
task.get_barx_j(Soltype::ITR, 0, barx.as_mut_slice())?; /* Request the interior solution. */
println!("Optimal primal solution");
println!(" x = {:?}",xx);
println!(" barx = {:?}",barx);
},
Solsta::DUAL_INFEAS_CER|Solsta::PRIM_INFEAS_CER =>
println!("Primal or dual infeasibility certificate found."),
Solsta::UNKNOWN =>
println!("The status of the solution could not be determined."),
_ =>
println!("Other solution status.")
}
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
sensitivity.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : sensitivity.rs
//!
//! Purpose : To demonstrate how to perform sensitivity
//! analysis from the API on a small problem:
//!
//! minimize
//!
//! obj: +1 x11 + 2 x12 + 5 x23 + 2 x24 + 1 x31 + 2 x33 + 1 x34
//! st
//! c1: + x11 + x12 <= 400
//! c2: + x23 + x24 <= 1200
//! c3: + x31 + x33 + x34 <= 1000
//! c4: + x11 + x31 = 800
//! c5: + x12 = 100
//! c6: + x23 + x33 = 500
//! c7: + x24 + x34 = 500
//!
//! The example uses basis type sensitivity analysis.
//!
extern crate mosek;
use mosek::{Task,Boundkey,Streamtype,Mark,Objsense};
const INFINITY : f64 = 0.0;
fn main() -> Result<(),String> {
let bkc = vec![
Boundkey::UP, Boundkey::UP,
Boundkey::UP, Boundkey::FX,
Boundkey::FX, Boundkey::FX,
Boundkey::FX ];
let bkx = vec![
Boundkey::LO, Boundkey::LO,
Boundkey::LO, Boundkey::LO,
Boundkey::LO, Boundkey::LO,
Boundkey::LO ];
let ptr = [0i64, 2, 4, 6, 8, 10, 12, 14];
let sub = [0i32, 3, 0, 4, 1, 5, 1, 6, 2, 3, 2, 5, 2, 6];
let blc = [ -INFINITY, -INFINITY, -INFINITY, 800.0, 100.0, 500.0, 500.0 ];
let buc = [400.0, 1200.0, 1000.0, 800.0, 100.0, 500.0, 500.0];
let c = [1.0, 2.0, 5.0, 2.0, 1.0, 2.0, 1.0];
let blx = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
let bux = [INFINITY, INFINITY,
INFINITY, INFINITY,
INFINITY, INFINITY,
INFINITY];
let val = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
let numcon = 7; /* Number of constraints. */
let numvar = 7; /* Number of variables. */
/* Create the optimization task. */
let mut task = match Task::new() {
Some(t) => t,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
/* Directs the log task stream to the 'printstr' function. */
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
task.input_data(numcon as i32, numvar as i32,
&c,
0.0,
&ptr[0..numvar as usize],
&ptr[1..numvar as usize+1],
&sub,
&val,
&bkc,
&blc,
&buc,
&bkx,
&blx,
&bux)?;
/* A maximization problem */
task.put_obj_sense(Objsense::MINIMIZE)?;
task.optimize()?;
/* Analyze upper bound on c1 and the equality constraint on c4 */
let mut subi = vec![0i32, 3i32];
let mut marki = vec![Mark::UP, Mark::UP];
/* Analyze lower bound on the variables x12 and x31 */
let mut subj = vec![1i32, 4];
let mut markj = vec![Mark::LO, Mark::LO];
let mut leftpricei = vec![0.0; 2];
let mut rightpricei = vec![0.0; 2];
let mut leftrangei = vec![0.0; 2];
let mut rightrangei = vec![0.0; 2];
let mut leftpricej = vec![0.0; 2];
let mut rightpricej = vec![0.0; 2];
let mut leftrangej = vec![0.0; 2];
let mut rightrangej = vec![0.0; 2];
task.primal_sensitivity(subi.as_mut_slice(),
marki.as_mut_slice(),
subj.as_mut_slice(),
markj.as_mut_slice(),
leftpricei.as_mut_slice(),
rightpricei.as_mut_slice(),
leftrangei.as_mut_slice(),
rightrangei.as_mut_slice(),
leftpricej.as_mut_slice(),
rightpricej.as_mut_slice(),
leftrangej.as_mut_slice(),
rightrangej.as_mut_slice())?;
println!("Results from sensitivity analysis on bounds:");
println!("For constraints:");
for i in 0..2 {
println!("leftprice = {:.5e}, rightprice = {:.5e}, leftrange = {:.5e}, rightrange = {:.5e}",
leftpricei[i], rightpricei[i], leftrangei[i], rightrangei[i]);
}
println!("For variables:\n");
for i in 0..2 {
println!("leftprice = {:.5e}, rightprice = {:.5e}, leftrange = {:.5e}, rightrange = {:.5e}",
leftpricej[i], rightpricej[i], leftrangej[i], rightrangej[i]);
}
let mut leftprice = vec![0.0; 2];
let mut rightprice = vec![0.0; 2];
let mut leftrange = vec![0.0; 2];
let mut rightrange = vec![0.0; 2];
let subc = [2i32, 5i32];
task.dual_sensitivity(&subc,
leftprice.as_mut_slice(),
rightprice.as_mut_slice(),
leftrange.as_mut_slice(),
rightrange.as_mut_slice())?;
println!("Results from sensitivity analysis on objective coefficients:");
for i in 0..2 {
println!("leftprice = {:.5e}, rightprice = {:.5e}, leftrange = {:.5e}, rightrange = {:.5e}",
leftprice[i], rightprice[i], leftrange[i], rightrange[i]);
}
return Result::Ok(());
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
simple.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : simple.rs
//!
//! Purpose : Demonstrates a very simple example using MOSEK by
//! reading a problem file, solving the problem and
//! writing the problem+solution to a file.
extern crate mosek;
use mosek::{Task,Streamtype,Iparam,Onoffkey};
use std::env;
enum FileOrText {
File(String),
Text(String)
}
fn main() -> Result<(),String> {
let mut args = env::args();
if args.len() < 2 {
println!("Syntax: simple FILENAME [ OUTFILE ]");
return Err("Invalid argument list".to_string())
}
let _ = args.next();
simple(FileOrText::File(args.next().unwrap()),
args.next())
}
fn simple(filename : FileOrText, outfile : Option<String>) -> Result<(),String> {
let mut task = Task::new().unwrap().with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
// We assume that a problem file was given as the first command
// line argument (received in `args')
match filename {
FileOrText::File(fname) => {
task.read_data(fname.as_str())?
},
FileOrText::Text(data) => {
task.read_ptf_string(data.as_str())?
}
}
// Solve the problem
let _ = task.optimize()?;
// Print a summary of the solution
task.solution_summary(Streamtype::LOG)?;
// If an output file was specified, save problem to file
if let Some(outfile) = outfile {
// If using OPF format, these parameters will specify what to include in output
task.put_int_param(Iparam::OPF_WRITE_SOLUTIONS, Onoffkey::ON)?;
task.put_int_param(Iparam::OPF_WRITE_PROBLEM, Onoffkey::ON)?;
task.put_int_param(Iparam::OPF_WRITE_HINTS, Onoffkey::OFF)?;
task.put_int_param(Iparam::OPF_WRITE_PARAMETERS, Onoffkey::OFF)?;
task.put_int_param(Iparam::PTF_WRITE_SOLUTIONS, Onoffkey::ON)?;
task.write_data(outfile.as_str())?;
}
Ok(())
}
#[cfg(test)]
mod tests {
const DFLT_FILE : &str = "Task
Objective
Maximize + 2 @x0 + 3 @x1 - @x2
Constraints
@c0 [1] + @x0 + @x1 + @x2
@C0 [QUAD(3)]
@ac1: + 0.03
@ac2: + 1.5 @x0 + 0.1 @x1
@ac3: + 0.3 @x0 + 2.1 @x2 + 0.1
Variables
@x0
@x1
@x2
";
#[test]
fn test() {
super::simple(super::FileOrText::Text(DFLT_FILE.to_string()),None).unwrap();
}
}
solutionquality.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : solutionquality.rs
//!
//! Purpose : To demonstrate how to examine the quality of a solution.
extern crate mosek;
use mosek::{Task,Streamtype,Solsta,Soltype};
use std::env;
enum FileOrText {
File(String),
Text(String)
}
fn main() -> Result<(),String> {
let mut args = env::args();
if args.len() < 2 {
println!("Syntax: solutionquality FILENAME");
return Err("Invalid argument list".to_string())
}
let _ = args.next();
let filename = args.next().unwrap();
solutionquality(FileOrText::File(filename))
}
fn solutionquality(filename : FileOrText) -> Result<(),String> {
let mut task = Task::new().unwrap().with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
// We assume that a problem file was given as the first command
// line argument (received in `args')
match filename {
FileOrText::File(filename) => task.read_data (filename.as_str())?,
FileOrText::Text(data) => task.read_ptf_string(data.as_str())?
}
// Solve the problem
let _ = task.optimize()?;
task.solution_summary(Streamtype::LOG)?;
let solsta = task.get_sol_sta(Soltype::BAS)?;
let mut pobj : f64 = 0.0;
let mut pviolcon : f64 = 0.0;
let mut pviolvar : f64 = 0.0;
let mut pviolbarvar : f64 = 0.0;
let mut pviolcones : f64 = 0.0;
let mut pviolitg : f64 = 0.0;
let mut dobj : f64 = 0.0;
let mut dviolcon : f64 = 0.0;
let mut dviolvar : f64 = 0.0;
let mut dviolbarvar : f64 = 0.0;
let mut dviolcones : f64 = 0.0;
task.get_solution_info(Soltype::BAS,
& mut pobj, & mut pviolcon, & mut pviolvar, & mut pviolbarvar, & mut pviolcones, & mut pviolitg,
& mut dobj, & mut dviolcon, & mut dviolvar, & mut dviolbarvar, & mut dviolcones)?;
match solsta {
Solsta::OPTIMAL => {
let abs_obj_gap = (dobj-pobj).abs();
let rel_obj_gap = abs_obj_gap / (1.0 + f64::min(pobj.abs(), dobj.abs()));
let max_primal_viol = f64::max(pviolcon, pviolvar);
let max_primal_viol = f64::max(max_primal_viol , pviolbarvar);
let max_primal_viol = f64::max(max_primal_viol , pviolcones);
let max_dual_viol = f64::max(dviolcon, dviolvar);
let max_dual_viol = f64::max(max_dual_viol, dviolbarvar);
let max_dual_viol = f64::max(max_dual_viol, dviolcones);
// Assume the application needs the solution to be within
// 1e-6 ofoptimality in an absolute sense. Another approach
// would be looking at the relative objective gap
println!("Customized solution information.");
println!(" Absolute objective gap: {:.3e}", abs_obj_gap);
println!(" Relative objective gap: {:.3e}", rel_obj_gap);
println!(" Max primal violation : {:.3e}", max_primal_viol);
println!(" Max dual violation : {:.3e}", max_dual_viol);
let mut accepted = true;
if rel_obj_gap > 1e-6 {
println!("Warning: The relative objective gap is LARGE.");
accepted = false;
}
// We will accept a primal infeasibility of 1e-8 and
// dual infeasibility of 1e-6. These number should chosen problem
// dependent.
if max_primal_viol > 1e-8 {
println!("Warning: Primal violation is too LARGE");
accepted = false;
}
if max_dual_viol > 1e-6 {
println!("Warning: Dual violation is too LARGE.");
accepted = false;
}
if accepted {
let numvar = task.get_num_var()?;
println!("Optimal primal solution");
let mut xx = vec![0.0; numvar as usize];
task.get_xx(Soltype::BAS,xx.as_mut_slice())?;
for (j,&xj) in (0..numvar).zip(xx.iter()) {
println!("x[{}]: {}",j,xj);
}
} else {
// print etailed information about the solution
task.analyze_solution(Streamtype::LOG, Soltype::BAS)?;
}
},
Solsta::DUAL_INFEAS_CER => println!("Primal or dual infeasibility certificate found."),
Solsta::UNKNOWN => println!("The status of the solution is unknown."),
_ => println!("Other solution status"),
}
Ok(())
}
#[cfg(test)]
mod tests {
const DFLT_FILE : &str = "Task
Objective
Maximize + 3 @x0 + @x1 + 5 @x2 + @x3
Constraints
@c0 [30] + 3 @x0 + @x1 + 2 @x2
@c1 [15;+inf] + 2 @x0 + @x1 + 3 @x2 + @x3
@c2 [-inf;25] + 2 @x1 + 3 @x3
Variables
@x0 [0;+inf]
@x1 [0;10]
@x2 [0;+inf]
@x3 [0;+inf]
";
#[test]
fn test() {
super::solutionquality(super::FileOrText::Text(DFLT_FILE.to_string())).unwrap();
}
}
solvebasis.rs
//!
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : solvebasis.rs
//!
extern crate mosek;
use mosek::{Task,Boundkey,Objsense,Streamtype};
/// Demonstrate the usage of solve_with_basis on the problem:
///
/// ```
/// maximize x0 + x1
/// st.
/// x0 + 2.0 x1 <= 2
/// x0 + x1 <= 6
/// x0 >= 0, x1>= 0
/// ```
///
/// The problem has the slack variables `xc0`, `xc1` on the
/// constraints and the variabels `x0` and `x1`.
///
/// ```
/// maximize x0 + x1
/// st.
/// x0 + 2.0 x1 -xc1 = 2
/// x0 + x1 -xc2 = 6
/// x0 >= 0, x1>= 0,
/// xc1 <= 0 , xc2 <= 0
/// ```
///
fn solve() -> Result<(),String> {
let mut task = match Task::new() {
Some(e) => e,
None => return Err("Failed to create task".to_string()),
}.with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
task.put_obj_name("solvebasis")?;
let numcon : i32 = 2;
let numvar : i32 = 2;
let mut w1 = vec![2.0, 6.0];
let mut w2 = vec![1.0, 0.0];
task.input_data(numcon, numvar,
&[1.0, 1.0], // c
0.0, // cfix
&[0,2], // ptrb
&[2,3], // ptre
&[0,1,
0,1], // sub
&[1.0, 1.0,
2.0, 1.0], // val
&[Boundkey::UP,
Boundkey::UP], // bkc
&[0.0,0.0], // blc
&[2.0,6.0], // buc
&[Boundkey::LO,
Boundkey::LO], // bkx
&[0.0, 0.0], // blx
&[0.0,0.0])?; // bux;
task.put_obj_sense(Objsense::MAXIMIZE)?;
let _ = task.optimize()?;
let mut basis = vec![0i32; numcon as usize];
task.init_basis_solve(basis.as_mut_slice())?;
// List basis variables corresponding to columns of B
let mut varsub = vec![0i32, 1];
for i in 0..numcon {
if basis[varsub[i as usize] as usize] < numcon {
println!("Basis variable no {} is xc{}",i,basis[varsub[i as usize] as usize]);
}
else {
println!("Basis variable no {} is x{}",i,basis[i as usize]-numcon);
// solve Bx = w1
// varsub contains index of non-zeros in b.
// On return b contains the solution x and
// varsub the index of the non-zeros in x.
{
let nz = task.solve_with_basis(false, 2, varsub.as_mut_slice(), w1.as_mut_slice())?;
println!("nz = {}",nz);
println!("Solution to Bx = {:?}",w1);
for vsubi in &varsub[0..nz as usize] {
if basis[*vsubi as usize] < numcon {
println!("xc {} = {}",basis[*vsubi as usize],w1[*vsubi as usize]);
}
else {
println!("x{} = {}",basis[*vsubi as usize] - numcon,w1[*vsubi as usize])
}
}
}
// Solve B^Tx = w2
{
varsub[0] = 1;
let nz = task.solve_with_basis(true,1,varsub.as_mut_slice(),w2.as_mut_slice())?;
println!("nz = {}",nz);
println!("Solution to B^Tx = {:?}",w2);
for vsubi in &varsub[0..nz as usize] {
if basis[*vsubi as usize] < numcon {
print!("xc{} = {}",basis[*vsubi as usize],w2[*vsubi as usize]);
}
else {
print!("x{} = {}",basis[*vsubi as usize] - numcon,w2[*vsubi as usize]);
}
}
}
}
}
Ok(())
}
fn main() -> Result<(),String> {
solve()
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
solvelinear.rs
//! Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//!
//! File : solvelinear.rs
//!
//! Purpose : To demonstrate the usage of MSK_solvewithbasis
//! when solving the linear system:
//!
//! 1.0 x1 = b1
//! -1.0 x0 + 1.0 x1 = b2
//!
//! with two different right hand sides
//!
//! b = (1.0, -2.0)
//!
//! and
//!
//! b = (7.0, 0.0)
extern crate mosek;
use mosek::{Task,Boundkey,Streamtype,Soltype,Stakey};
const INF : f64 = 0.0;
fn setup(task : & mut mosek::Task,
aval : &[f64],
asub : &[i32],
ptrb : &[i64],
ptre : &[i64],
numvar : i32) -> Result<Vec<i32>,String> {
// mosek.stakey[] skx = new mosek.stakey [numvar];
// mosek.stakey[] skc = new mosek.stakey [numvar];
// for (int i = 0; i < numvar ; ++i) {
// skx[i] = mosek.stakey.bas;
// skc[i] = mosek.stakey.fix;
// }
task.append_vars(numvar)?;
task.append_cons(numvar)?;
task.put_a_col_slice(0,numvar,ptrb,ptre,asub,aval)?;
task.put_con_bound_slice_const(0,numvar,Boundkey::FX,0.0,0.0)?;
task.put_var_bound_slice_const(0,numvar,Boundkey::FR,-INF,INF)?;
/* Define a basic solution by specifying
status keys for variables & constraints. */
task.delete_solution(Soltype::BAS)?;
task.put_skc_slice(Soltype::BAS, 0, numvar, vec![Stakey::FIX; numvar as usize].as_slice())?;
task.put_skx_slice(Soltype::BAS, 0, numvar, vec![Stakey::BAS; numvar as usize].as_slice())?;
let mut basis = vec![0; numvar as usize];
task.init_basis_solve(basis.as_mut_slice())?;
Ok(basis)
}
fn main() -> Result<(),String> {
let numcon : i32 = 2;
let numvar : i32 = 2;
let aval = &[ -1.0 ,
1.0, 1.0 ];
let asub = &[ 1,
0, 1 ];
let ptrb = &[0, 1];
let ptre = &[1, 3];
// int[] bsub = new int[numvar];
// double[] b = new double[numvar];
// int[] basis = new int[numvar];
let mut task = Task::new().unwrap();
// Put A matrix and factor A. Call this function only once for a
// given task.
let basis = setup(& mut task,
aval,
asub,
ptrb,
ptre,
numvar)?;
let mut task = task.with_callbacks();
task.put_stream_callback(Streamtype::LOG, |msg| print!("{}",msg))?;
/* now solve rhs */
let mut b = vec![0.0; numvar as usize];
let mut bsub = vec![0; numvar as usize];
b[0] = 1.0; b[1] = -2.0;
bsub[0] = 0; bsub[1] = 1;
let nz = task.solve_with_basis(false, 2, bsub.as_mut_slice(), b.as_mut_slice())?;
println!("\nSolution to Bx = b:\n");
// Print solution and show correspondents to original variables in
// the problem
for &bsubi in bsub[..nz as usize].iter() {
if basis[bsubi as usize] < numcon {
println!("This should never happen");
}
else {
println!("x{} = {}",basis[bsubi as usize] - numcon, b[bsubi as usize]);
}
}
b[0] = 7.0;
bsub[0] = 0;
let nz = task.solve_with_basis(false, 1, bsub.as_mut_slice(), b.as_mut_slice())?;
println!("Solution to Bx = b:");
// Print solution and show correspondents to original variables in
// the problem
for &bsubi in bsub[..nz as usize].iter() {
if basis[bsubi as usize] < numcon {
println!("This should never happen");
}
else {
println!("x{} = {}",basis[bsubi as usize] - numcon,b[bsubi as usize]);
}
}
Ok(())
}
#[cfg(test)]
mod tests {
#[test]
fn test() {
super::main().unwrap();
}
}
sparsecholesky.rs
/*
Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
File: sparsecholesky.rs
Purpose: Demonstrate the sparse Cholesky factorization.
*/
extern crate mosek;
extern crate itertools;
use mosek::Transpose;
use itertools::{iproduct,izip};
fn print_sparse(n : usize,
perm : &[i32],
diag : &[f64],
lnzc : &[i32],
lptrc : &[i64],
lsubc : &[i32],
lvalc : &[f64]) {
println!("P = {:?}",perm);
println!("diag(D) = {:?}",diag);
let mut l = vec![0.0; n*n];
for (j,(&ptr,&sz)) in lptrc.iter().zip(lnzc.iter()).enumerate() {
let pfrom = ptr as usize;
let pto = ptr as usize + sz as usize;
for (&subci,&valci) in lsubc[pfrom..pto].iter().zip(lvalc[pfrom..pto].iter()) {
l[subci as usize * n + j] = valci;
}
};
println!("L = [");
(0..n*n).step_by(n).for_each(|i| {
print!(" {:?}",&l[i..i+n]);
});
println!(" ]");
}
fn main() -> Result<(),String> {
/* Create the mosek environment. */
//Example from the manual
//Observe that anzc, aptrc, asubc and avalc only specify the lower triangular part.
let n = 4;
let anzc = [4, 1, 1, 1];
let asubc = [0, 1, 2, 3, 1, 2, 3];
let aptrc = [0, 4, 5, 6];
let avalc = [4.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
let b = [13.0, 3.0, 4.0, 5.0];
let mut perm = Vec::new();
let mut lnzc = Vec::new();
let mut lsubc = Vec::new();
let mut diag = Vec::new();
let mut lptrc = Vec::new();
let mut lvalc = Vec::new();
let mut lensubnval = 0;
mosek::compute_sparse_cholesky(0, //Mosek chooses number of threads
1, //Apply reordering heuristic
1.0e-14, //Singularity tolerance
&anzc, &aptrc, &asubc, &avalc,
& mut perm, & mut diag,
& mut lnzc, & mut lptrc, & mut lensubnval, & mut lsubc, & mut lvalc)?;
print_sparse(n, perm.as_slice(), diag.as_slice(), lnzc.as_slice(), lptrc.as_slice(), lsubc.as_slice(), lvalc.as_slice());
/* Permuted b is stored as x. */
let mut x : Vec<f64> = perm.iter().map(|&i| b[i as usize]).collect();
/*Compute inv(L)*x.*/
mosek::sparse_triangular_solve_dense(Transpose::NO, lnzc.as_slice(), lptrc.as_slice(), lsubc.as_slice(), lvalc.as_slice(), x.as_mut_slice())?;
/*Compute inv(L^T)*x.*/
mosek::sparse_triangular_solve_dense(Transpose::YES, lnzc.as_slice(), lptrc.as_slice(), lsubc.as_slice(), lvalc.as_slice(), x.as_mut_slice())?;
print!("\nSolution A x = b, x = [ {:?} ]",
iproduct!(0..n,izip!(perm.iter(),x.iter()))
.filter_map(|(i,(&pj,&xj))|if pj as usize == i { Some(xj) } else { None })
.collect::<Vec<f64>>());
let n = 3;
let anzc = [3, 2, 1];
let asubc = [0, 1, 2, 1, 2, 2];
let aptrc = [0, 3, 5, ];
let avalc = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
let mut perm = Vec::new();
let mut lnzc = Vec::new();
let mut lsubc = Vec::new();
let mut diag = Vec::new();
let mut lptrc = Vec::new();
let mut lvalc = Vec::new();
let mut lensubnval = 0;
mosek::compute_sparse_cholesky(0, //Mosek chooses number of threads
1, //Apply reordering heuristic
1.0e-14, //Singularity tolerance
&anzc, &aptrc, &asubc, &avalc,
& mut perm, & mut diag,
& mut lnzc, & mut lptrc, & mut lensubnval, & mut lsubc, & mut lvalc)?;
print_sparse(n, perm.as_slice(), diag.as_slice(), lnzc.as_slice(), lptrc.as_slice(), lsubc.as_slice(), lvalc.as_slice());
Ok(())
}
