Merge pull request #116 from oxfordcontrol/develop

v0.9.0

.github/workflows/ci.yml

@@ -28,7 +28,7 @@ jobs:
 
     - name: Generate code coverage 
       run: |
-        cargo tarpaulin --out xml --features sdp-accelerate --exclude-files "src/python/*,src/julia/*"
+        cargo tarpaulin --out xml --features sdp-accelerate,serde,faer-sparse --exclude-files "src/python/*,src/julia/*"
 
     - name: Upload to codecov.io
       uses: codecov/codecov-action@v4

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "clarabel"
-version = "0.8.1"
+version = "0.9.0"
 authors = ["Paul Goulart <[email protected]>"]
 edition = "2021"
 rust-version = "1.66"
@@ -29,6 +29,10 @@ itertools      = "0.11"
 # -------------------------------
 
 [features]
+default = ["serde"]
+
+# enable reading / writing of problems from json files 
+serde = ["dep:serde", "dep:serde_json"]
 
 # enables blas/lapack for SDP support, with blas/lapack src unspecified
 # also enable packages required for chordal decomposition 
@@ -41,15 +45,20 @@ sdp-openblas   = ["sdp", "blas-src/openblas", "lapack-src/openblas"]
 sdp-mkl        = ["sdp", "blas-src/intel-mkl", "lapack-src/intel-mkl"]
 sdp-r          = ["sdp", "blas-src/r", "lapack-src/r"]
 
+
 # build as the julia interface 
-julia = ["sdp", "dep:libc", "dep:num-derive", "dep:serde", "dep:serde_json"] 
+julia = ["sdp", "dep:libc", "dep:num-derive",  "serde", "faer-sparse"] 
 
 # build as the python interface via maturin.
 # NB: python builds use scipy shared libraries
 # for blas/lapack, and should *not* explicitly 
 # enable a blas/lapack source package 
-python = ["sdp", "dep:libc", "dep:pyo3", "dep:num-derive"]
+python = ["sdp", "dep:libc", "dep:pyo3", "dep:num-derive", "serde", "faer-sparse"]
+
+wasm = ["dep:web-time"]
 
+#compile with faer supernodal solver option
+faer-sparse = ["dep:faer", "dep:faer-entity"]
 
 # -------------------------------
 # SDP configuration
@@ -75,6 +84,14 @@ optional = true
 version = "2.2"
 optional = true
 
+[dependencies.faer]
+version = "0.19"
+optional = true 
+
+[dependencies.faer-entity]
+version = "0.19"
+optional = true 
+
 
 # -------------------------------
 # examples
@@ -109,6 +126,18 @@ name = "sdp"
 path = "examples/rust/example_sdp.rs"
 required-features = ["sdp"]
 
+[[example]]
+name = "box_faer"
+path = "examples/rust/example_box_faer.rs"
+required-features = ["faer-sparse"]
+
+[[example]]
+name = "json"
+path = "examples/rust/example_json.rs"
+required-features = ["serde"]
+
+
+
 
 # -------------------------------
 # custom build profiles 
@@ -162,3 +191,17 @@ crate-type = ["lib","cdylib"]
 rustdoc-args = [ "--html-in-header", "./html/rustdocs-header.html" ]
 features = ["sdp","sdp-mkl"]
 
+# ------------------------------
+# wasm compatibility
+# ------------------------------
+
+[dependencies.web-time]
+optional = true
+version = "0.2.3"
+
+# ------------------------------
+# testing, benchmarking etc 
+# ------------------------------
+[dev-dependencies]
+tempfile = "3"
+

README.md

@@ -12,7 +12,7 @@ Interior Point Conic Optimization for Rust and Python
   <a href="https://codecov.io/gh/oxfordcontrol/Clarabel.rs"><img src="https://codecov.io/gh/oxfordcontrol/Clarabel.rs/branch/main/graph/badge.svg"></a>
   <a href="https://clarabel.org"><img src="https://img.shields.io/badge/Documentation-stable-purple.svg"></a>
   <a href="https://opensource.org/licenses/Apache-2.0"><img src="https://img.shields.io/badge/License-Apache%202.0-blue.svg"></a>
-  <a href="https://github.com/oxfordcontrol/Clarabel.rs/releases"><img src="https://img.shields.io/badge/Release-v0.8.1-blue.svg"></a>
+  <a href="https://github.com/oxfordcontrol/Clarabel.rs/releases"><img src="https://img.shields.io/badge/Release-v0.9.0-blue.svg"></a>
 </p>
 
 <p align="center">

examples/data/hs35.json

@@ -0,0 +1 @@
+{"P":{"m":3,"n":3,"colptr":[0,1,3,5],"rowval":[0,0,1,0,2],"nzval":[4.000000000000001,2.0000000000000004,4.000000000000001,2.0,2.0]},"q":[-8.0,-6.0,-4.0],"A":{"m":4,"n":3,"colptr":[0,2,4,6],"rowval":[0,1,0,2,0,3],"nzval":[1.0,-1.0,1.0,-1.0,2.0,-1.0]},"b":[3.0,0.0,0.0,0.0],"cones":[{"NonnegativeConeT":4}],"settings":{"max_iter":200,"time_limit":1.7976931348623157e308,"verbose":true,"max_step_fraction":0.99,"tol_gap_abs":1e-8,"tol_gap_rel":1e-8,"tol_feas":1e-8,"tol_infeas_abs":1e-8,"tol_infeas_rel":1e-8,"tol_ktratio":1e-6,"reduced_tol_gap_abs":0.00005,"reduced_tol_gap_rel":0.00005,"reduced_tol_feas":0.0001,"reduced_tol_infeas_abs":0.00005,"reduced_tol_infeas_rel":0.00005,"reduced_tol_ktratio":0.0001,"equilibrate_enable":true,"equilibrate_max_iter":10,"equilibrate_min_scaling":0.0001,"equilibrate_max_scaling":10000.0,"linesearch_backtrack_step":0.8,"min_switch_step_length":0.1,"min_terminate_step_length":0.0001,"direct_kkt_solver":true,"direct_solve_method":"qdldl","static_regularization_enable":true,"static_regularization_constant":1e-8,"static_regularization_proportional":4.930380657631324e-32,"dynamic_regularization_enable":true,"dynamic_regularization_eps":1e-13,"dynamic_regularization_delta":2e-7,"iterative_refinement_enable":true,"iterative_refinement_reltol":1e-13,"iterative_refinement_abstol":1e-12,"iterative_refinement_max_iter":10,"iterative_refinement_stop_ratio":5.0,"presolve_enable":true,"chordal_decomposition_enable":true,"chordal_decomposition_merge_method":"clique_graph","chordal_decomposition_compact":true,"chordal_decomposition_complete_dual":true}}

examples/python/example_json.py

@@ -0,0 +1,14 @@
+import clarabel
+import os
+
+thisdir = os.path.dirname(__file__)
+filename = os.path.join(thisdir, "../data/hs35.json")
+print(filename)
+
+# Load problem data from JSON file
+solver = clarabel.read_from_file(filename)
+solution = solver.solve()
+
+# export problem data to JSON file
+# filename = os.path.join(thisdir, "../data/out.json")
+# solver.write_to_file(filename)

examples/rust/example_box.rs

@@ -3,7 +3,7 @@ use clarabel::algebra::*;
 use clarabel::solver::*;
 
 fn problem_data() -> (CscMatrix<f64>, Vec<f64>, CscMatrix<f64>, Vec<f64>) {
-    let n = 20000;
+    let n = 2000000;
 
     let P = CscMatrix::identity(n);
 

examples/rust/example_box_faer.rs

@@ -0,0 +1,39 @@
+#![allow(non_snake_case)]
+use clarabel::algebra::*;
+use clarabel::solver::*;
+
+fn problem_data() -> (CscMatrix<f64>, Vec<f64>, CscMatrix<f64>, Vec<f64>) {
+    let n = 2000000;
+
+    let P = CscMatrix::identity(n);
+
+    // construct A = [I; -I]
+    let I1 = CscMatrix::<f64>::identity(n);
+    let mut I2 = CscMatrix::<f64>::identity(n);
+    I2.negate();
+
+    let A = CscMatrix::vcat(&I1, &I2);
+
+    let q = vec![1.; n];
+    let b = vec![1.; 2 * n];
+
+    (P, q, A, b)
+}
+
+fn main() {
+    let (P, q, A, b) = problem_data();
+
+    let cones = [NonnegativeConeT(b.len())];
+
+    let settings = DefaultSettingsBuilder::default()
+        .direct_solve_method("faer".to_owned())
+        .equilibrate_enable(true)
+        .max_iter(50)
+        .verbose(true)
+        .build()
+        .unwrap();
+
+    let mut solver = DefaultSolver::new(&P, &q, &A, &b, &cones, settings);
+
+    solver.solve();
+}

examples/rust/example_json.rs

@@ -0,0 +1,19 @@
+#![allow(non_snake_case)]
+use clarabel::solver::*;
+use std::fs::File;
+
+fn main() {
+    // HS35 is a small problem QP problem
+    // from the Maros-Meszaros test set
+
+    let filename = "./examples/data/hs35.json";
+    let mut file = File::open(filename).unwrap();
+    let mut solver = DefaultSolver::<f64>::read_from_file(&mut file).unwrap();
+    solver.solve();
+
+    // to write the back to a new file
+
+    // let outfile = "./examples/data/output.json";
+    // let mut file = File::create(outfile).unwrap();
+    // solver.write_to_file(&mut file).unwrap();
+}

src/algebra/csc/core.rs

@@ -6,6 +6,9 @@ use crate::algebra::{Adjoint, MatrixShape, ShapedMatrix, SparseFormatError, Symm
 use num_traits::Num;
 use std::iter::{repeat, zip};
 
+#[cfg(feature = "serde")]
+use serde::{de::DeserializeOwned, Deserialize, Serialize};
+
 /// Sparse matrix in standard Compressed Sparse Column (CSC) format
 ///
 /// __Example usage__ : To construct the 3 x 3 matrix
@@ -38,6 +41,8 @@ use std::iter::{repeat, zip};
 /// ```
 ///
 
+#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
+#[serde(bound = "T: Serialize + DeserializeOwned")]
 #[derive(Debug, Clone, PartialEq, Eq)]
 pub struct CscMatrix<T = f64> {
     /// number of rows
@@ -269,7 +274,7 @@ where
 
     /// Return matrix data in triplet format.
     ///
-    #[cfg_attr(not(sdp), allow(dead_code))]
+    #[cfg_attr(not(feature = "sdp"), allow(dead_code))]
     pub(crate) fn findnz(&self) -> (Vec<usize>, Vec<usize>, Vec<T>) {
         let I = self.rowval.clone();
         let mut J = Vec::with_capacity(self.nnz());

src/algebra/floats.rs

@@ -40,53 +40,50 @@ impl<T> CoreFloatT for T where
 {
 }
 
+// additional traits that add bounds to CoreFloatT
+// when optional dependencies are enabled
+
 // if "sdp" is enabled, we must add an additional trait
 // trait bound to restrict compilation for f32/f64 types
 // since there is no BLAS support otherwise
 
-// Define the documentation string as a macro so that FloatT gets documentation
-// regardless of whether the "sdp" feature is enabled.
-macro_rules! floatT_doc_header {
-    () => {
-        r#"Main trait for floating point types used in the Clarabel solver."#
-    };
-}
-macro_rules! floatT_doc_long {
-    () => {
-        r#"All floating point calculations in Clarabel are represented internally on values 
-         implementing the `FloatT` trait, with implementations provided only for f32 and f64 
-         native types when compiled with BLAS/LAPACK support for SDPs. If SDP support is not 
-         enabled then it should be possible to compile Clarabel to support any any other 
-         floating point type provided that it satisfies the trait bounds of `CoreFloatT`. 
-        \
-        \
-         `FloatT` relies on [`num_traits`](num_traits) for most of its constituent trait bounds."#
-    };
-}
 cfg_if::cfg_if! {
-    if #[cfg(not(feature="sdp"))] {
-    #[doc = floatT_doc_header!()]
-    ///
-    #[doc = floatT_doc_long!()]
-        pub trait FloatT: CoreFloatT {}
-    } else{
-        #[doc = floatT_doc_header!()]
-        ///
-        #[doc = floatT_doc_long!()]
-        ///
-        /// The trait bound `BlasFloatT` is only enforced when compiling with the `sdp` feature.
-        pub trait FloatT: CoreFloatT + BlasFloatT {}
+    if #[cfg(feature="sdp")] {
+        pub trait MaybeBlasFloatT : BlasFloatT {}
+        impl<T> MaybeBlasFloatT for T where T: BlasFloatT {}
+    }
+    else {
+        pub trait MaybeBlasFloatT {}
+        impl<T> MaybeBlasFloatT for T {}
     }
 }
 
+// if "faer" is enabled, we must add an additional bound
+// to restrict compilation to Reals implementing SimpleEntity
+
 cfg_if::cfg_if! {
-    if #[cfg(feature="sdp")] {
-        impl<T> FloatT for T where T: CoreFloatT + BlasFloatT {}
-    } else{
-        impl<T> FloatT for T where T: CoreFloatT {}
+    if #[cfg(feature="faer-sparse")] {
+        pub trait MaybeFaerFloatT : faer::SimpleEntity + faer::ComplexField<Real=Self> {}
+        impl<T> MaybeFaerFloatT for T where T: faer::SimpleEntity + faer::ComplexField<Real=Self> {}
+    }
+    else {
+        pub trait MaybeFaerFloatT {}
+        impl<T> MaybeFaerFloatT for T {}
     }
 }
 
+/// Main trait for floating point types used in the Clarabel solver.
+///
+/// All floating point calculations in Clarabel are represented internally on values
+/// implementing the `FloatT` trait, with implementations provided only for f32 and f64
+/// native types when compiled with BLAS/LAPACK support for SDPs. If SDP support is not
+/// enabled then it should be possible to compile Clarabel to support any any other
+/// floating point type provided that it satisfies the trait bounds of `CoreFloatT`.
+///
+/// `FloatT` relies on [`num_traits`](num_traits) for most of its constituent trait bounds.
+pub trait FloatT: CoreFloatT + MaybeBlasFloatT + MaybeFaerFloatT {}
+impl<T> FloatT for T where T: CoreFloatT + MaybeBlasFloatT + MaybeFaerFloatT {}
+
 /// Trait for convering Rust primitives to [`FloatT`](crate::algebra::FloatT)
 ///
 /// This convenience trait is implemented on f32/64 and u32/64.  This trait

src/algebra/matrix_traits.rs

@@ -3,7 +3,7 @@
 use crate::algebra::MatrixConcatenationError;
 use crate::algebra::MatrixShape;
 
-#[cfg_attr(not(sdp), allow(dead_code))]
+#[cfg_attr(not(feature = "sdp"), allow(dead_code))]
 pub(crate) trait ShapedMatrix {
     fn shape(&self) -> MatrixShape;
     fn size(&self) -> (usize, usize);

src/algebra/scalarmath.rs

@@ -24,7 +24,7 @@ pub(crate) fn triangular_number(k: usize) -> usize {
     (k * (k + 1)) >> 1
 }
 
-#[cfg_attr(not(sdp), allow(dead_code))]
+#[cfg_attr(not(feature = "sdp"), allow(dead_code))]
 pub(crate) fn triangular_index(k: usize) -> usize {
     // 0-based index into a packed triangle. Same as:
     // triangular number(k+1) - 1 = (((k+1) * (k+2)) >> 1) - 1
@@ -33,7 +33,7 @@ pub(crate) fn triangular_index(k: usize) -> usize {
 
 // given an index into the upper triangular part of a matrix, return
 // its row and column position
-#[cfg_attr(not(sdp), allow(dead_code))]
+#[cfg_attr(not(feature = "sdp"), allow(dead_code))]
 pub(crate) fn upper_triangular_index_to_coord(linearidx: usize) -> (usize, usize) {
     if linearidx == 0 {
         return (0, 0);
@@ -47,7 +47,7 @@ pub(crate) fn upper_triangular_index_to_coord(linearidx: usize) -> (usize, usize
 
 // given a row and column position, return the index into the upper
 // triangular part of the matrix
-#[cfg_attr(not(sdp), allow(dead_code))]
+#[cfg_attr(not(feature = "sdp"), allow(dead_code))]
 pub(crate) fn coord_to_upper_triangular_index(coord: (usize, usize)) -> usize {
     if coord == (0, 0) {
         return 0;
@@ -65,7 +65,7 @@ pub(crate) fn coord_to_upper_triangular_index(coord: (usize, usize)) -> usize {
 // here: https://github.com/rust-lang/rust/issues/116226
 // For now, implement a simple truncation method, which works
 // for inputs < 2^53 or so (otherwise possibly off-by-one).
-#[cfg_attr(not(sdp), allow(dead_code))]
+#[cfg_attr(not(feature = "sdp"), allow(dead_code))]
 fn isqrt(v: usize) -> usize {
     (v as f64).sqrt() as usize
 }

src/algebra/tests/vector.rs

@@ -140,6 +140,9 @@ fn test_norm_scaled() {
 fn test_norm_inf() {
     let x = [-3., 4., -12.];
     assert_eq!(x.norm_inf(), 12.);
+
+    let x = [-3., f64::NAN, -12.];
+    assert!(x.norm_inf().is_nan());
 }
 
 #[test]