From b033674516a3fa9ff1f2347c93f882fcb7a42668 Mon Sep 17 00:00:00 2001 From: goulart-paul Date: Tue, 28 May 2024 19:25:10 +0100 Subject: [PATCH 01/13] bump version --- Cargo.toml | 2 +- README.md | 2 +- src/julia/ClarabelRs/Project.toml | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 45d9630..09bf562 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "clarabel" -version = "0.8.1" +version = "0.9" authors = ["Paul Goulart "] edition = "2021" rust-version = "1.66" diff --git a/README.md b/README.md index 1de745b..9ba3270 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ Interior Point Conic Optimization for Rust and Python - +

diff --git a/src/julia/ClarabelRs/Project.toml b/src/julia/ClarabelRs/Project.toml index 450a8b0..c42eaf4 100644 --- a/src/julia/ClarabelRs/Project.toml +++ b/src/julia/ClarabelRs/Project.toml @@ -1,7 +1,7 @@ name = "ClarabelRs" uuid = "a0c58a0a-712c-48b7-9fd4-64369ecb2011" authors = ["Paul Goulart "] -version = "0.8.1" +version = "0.9" [deps] Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" @@ -11,4 +11,4 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] -Clarabel = "0.8.1" +Clarabel = "0.9" From 0e2c26990b3792c0962b7cbf5a36f2d38f438ad5 Mon Sep 17 00:00:00 2001 From: Paul Goulart Date: Tue, 28 May 2024 19:47:39 +0100 Subject: [PATCH 02/13] validation tools for solver settings (#113) --- src/julia/ClarabelRs/src/types.jl | 5 + src/julia/interface.rs | 11 +- src/python/impl_default_py.rs | 15 ++- .../implementations/default/settings.rs | 110 ++++++++++++++++++ 4 files changed, 136 insertions(+), 5 deletions(-) diff --git a/src/julia/ClarabelRs/src/types.jl b/src/julia/ClarabelRs/src/types.jl index e188666..5819f76 100644 --- a/src/julia/ClarabelRs/src/types.jl +++ b/src/julia/ClarabelRs/src/types.jl @@ -107,6 +107,11 @@ mutable struct Solver{T <: Float64} <: Clarabel.AbstractSolver{Float64} ptr:: Ptr{Cvoid} function Solver{T}(ptr) where T + + if ptr == C_NULL + throw(ErrorException("Solver constructor failed")) + end + obj = new(ptr) finalizer(solver_drop_jlrs,obj) return obj diff --git a/src/julia/interface.rs b/src/julia/interface.rs index 2c4a2d2..e2cde05 100644 --- a/src/julia/interface.rs +++ b/src/julia/interface.rs @@ -62,11 +62,18 @@ pub(crate) extern "C" fn solver_new_jlrs( let b = Vec::from(b); let cones = ccall_arrays_to_cones(jlcones); - let settings = settings_from_json(json_settings); - let solver = DefaultSolver::new(&P, &q, &A, &b, &cones, settings); + // manually validate settings from Julia side + match settings.validate() { + Ok(_) => (), + Err(e) => { + println!("Invalid settings: {}", e); + return std::ptr::null_mut(); + } + }; + let solver = DefaultSolver::new(&P, &q, &A, &b, &cones, settings); to_ptr(Box::new(solver)) } diff --git a/src/python/impl_default_py.rs b/src/python/impl_default_py.rs index 198ed2a..72a13fb 100644 --- a/src/python/impl_default_py.rs +++ b/src/python/impl_default_py.rs @@ -13,6 +13,7 @@ use crate::solver::{ }; use num_derive::ToPrimitive; use num_traits::ToPrimitive; +use pyo3::exceptions::PyException; use pyo3::prelude::*; use std::fmt::Write; @@ -401,12 +402,20 @@ impl PyDefaultSolver { b: Vec, cones: Vec, settings: PyDefaultSettings, - ) -> Self { + ) -> PyResult { let cones = _py_to_native_cones(cones); let settings = settings.to_internal(); - let solver = DefaultSolver::new(&P, &q, &A, &b, &cones, settings); - Self { inner: solver } + //manually validate settings from Python side + match settings.validate() { + Ok(_) => (), + Err(e) => { + return Err(PyException::new_err(format!("Invalid settings: {}", e))); + } + } + + let solver = DefaultSolver::new(&P, &q, &A, &b, &cones, settings); + Ok(Self { inner: solver }) } fn solve(&mut self) -> PyDefaultSolution { diff --git a/src/solver/implementations/default/settings.rs b/src/solver/implementations/default/settings.rs index 17a0528..fc6518b 100644 --- a/src/solver/implementations/default/settings.rs +++ b/src/solver/implementations/default/settings.rs @@ -8,6 +8,7 @@ use serde::{Deserialize, Serialize}; /// Standard-form solver type implementing the [`Settings`](crate::solver::core::traits::Settings) trait #[derive(Builder, Debug, Clone)] +#[builder(build_fn(validate = "Self::validate"))] #[cfg_attr(feature = "julia", derive(Serialize, Deserialize))] pub struct DefaultSettings { ///maximum number of iterations @@ -204,3 +205,112 @@ where self } } + +// pre build checker (for auto-validation when using the builder) + +/// Automatic pre-build settings validation +impl DefaultSettingsBuilder +where + T: FloatT, +{ + pub fn validate(&self) -> Result<(), String> { + // check that the direct solve method is valid + if let Some(ref direct_solve_method) = self.direct_solve_method { + validate_direct_solve_method(direct_solve_method.as_str())?; + } + + // check that the chordal decomposition merge method is valid + #[cfg(feature = "sdp")] + if let Some(ref chordal_decomposition_merge_method) = + self.chordal_decomposition_merge_method + { + validate_chordal_decomposition_merge_method( + chordal_decomposition_merge_method.as_str(), + )?; + } + + Ok(()) + } +} + +// post build checker (for ad-hoc validation, e.g. when passing from python/Julia) +// this is not used directly in the solver, but can be called manually by the user + +/// Manual post-build settings validation +impl DefaultSettings +where + T: FloatT, +{ + pub fn validate(&self) -> Result<(), String> { + validate_direct_solve_method(&self.direct_solve_method)?; + + // check that the chordal decomposition merge method is valid + #[cfg(feature = "sdp")] + validate_chordal_decomposition_merge_method(&self.chordal_decomposition_merge_method)?; + + Ok(()) + } +} + +// --------------------------------------------------------- +// individual validation functions go here +// --------------------------------------------------------- + +fn validate_direct_solve_method(direct_solve_method: &str) -> Result<(), String> { + match direct_solve_method { + "qdldl" => Ok(()), + #[cfg(feature = "faer-sparse")] + "faer" => Ok(()), + _ => Err(format!( + "Invalid direct_solve_method: {:?}", + direct_solve_method + )), + } +} + +#[cfg(feature = "sdp")] +fn validate_chordal_decomposition_merge_method( + chordal_decomposition_merge_method: &str, +) -> Result<(), String> { + match chordal_decomposition_merge_method { + "none" => Ok(()), + "parent_child" => Ok(()), + "clique_graph" => Ok(()), + _ => Err(format!( + "Invalid chordal_decomposition_merge_method: {}", + chordal_decomposition_merge_method + )), + } +} + +#[test] +fn test_settings_validate() { + // all standard settings + DefaultSettingsBuilder::::default().build().unwrap(); + + // fail on unknown direct solve method + assert!(DefaultSettingsBuilder::::default() + .direct_solve_method("foo".to_string()) + .build() + .is_err()); + + // fail on solve options in disabled feature + let builder = DefaultSettingsBuilder::::default() + .direct_solve_method("faer".to_string()) + .build(); + cfg_if::cfg_if! { + if #[cfg(feature = "faer-sparse")] { + assert!(build.is_ok()); + } + else { + assert!(builder.is_err()); + } + } + + #[cfg(feature = "sdp")] + // fail on unknown chordal decomposition merge method + assert!(DefaultSettingsBuilder::::default() + .chordal_decomposition_merge_method("foo".to_string()) + .build() + .is_err()); +} From da5546c38773fe8fd92835274455acb564fa670b Mon Sep 17 00:00:00 2001 From: Paul Goulart Date: Tue, 28 May 2024 19:50:03 +0100 Subject: [PATCH 03/13] Read/write problems to JSON files (#111) --- .github/workflows/ci.yml | 2 +- Cargo.toml | 21 ++- examples/data/hs35.json | 1 + examples/python/example_json.py | 14 ++ examples/rust/example_json.rs | 19 +++ src/algebra/csc/core.rs | 5 + src/algebra/utils.rs | 4 +- src/julia/ClarabelRs/src/interface.jl | 47 ++++++ src/julia/interface.rs | 71 ++++++++- src/python/impl_default_py.rs | 13 ++ src/python/module_py.rs | 2 + src/solver/core/cones/supportedcone.rs | 5 +- src/solver/implementations/default/json.rs | 135 ++++++++++++++++++ src/solver/implementations/default/mod.rs | 3 + .../implementations/default/settings.rs | 7 +- 15 files changed, 338 insertions(+), 11 deletions(-) create mode 100644 examples/data/hs35.json create mode 100644 examples/python/example_json.py create mode 100644 examples/rust/example_json.rs create mode 100644 src/solver/implementations/default/json.rs diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e3ac262..e5a7285 100644 --- a/.github/workflows/ci.yml +++ b/.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 --exclude-files "src/python/*,src/julia/*" - name: Upload to codecov.io uses: codecov/codecov-action@v4 diff --git a/Cargo.toml b/Cargo.toml index 09bf562..a4f509c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -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,14 +45,16 @@ 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"] # 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"] + # ------------------------------- @@ -109,6 +115,12 @@ name = "sdp" path = "examples/rust/example_sdp.rs" required-features = ["sdp"] +[[example]] +name = "json" +path = "examples/rust/example_json.rs" +required-features = ["serde"] + + # ------------------------------- # custom build profiles @@ -162,3 +174,8 @@ crate-type = ["lib","cdylib"] rustdoc-args = [ "--html-in-header", "./html/rustdocs-header.html" ] features = ["sdp","sdp-mkl"] + +# ------------------------------ +# testing, benchmarking etc +[dev-dependencies] +tempfile = "3" \ No newline at end of file diff --git a/examples/data/hs35.json b/examples/data/hs35.json new file mode 100644 index 0000000..99e22d7 --- /dev/null +++ b/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}} \ No newline at end of file diff --git a/examples/python/example_json.py b/examples/python/example_json.py new file mode 100644 index 0000000..310f9d4 --- /dev/null +++ b/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) diff --git a/examples/rust/example_json.rs b/examples/rust/example_json.rs new file mode 100644 index 0000000..76f339c --- /dev/null +++ b/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::::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(); +} diff --git a/src/algebra/csc/core.rs b/src/algebra/csc/core.rs index 822fb24..ee87d28 100644 --- a/src/algebra/csc/core.rs +++ b/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 { /// number of rows diff --git a/src/algebra/utils.rs b/src/algebra/utils.rs index 62e5e82..94e63e8 100644 --- a/src/algebra/utils.rs +++ b/src/algebra/utils.rs @@ -106,8 +106,8 @@ fn test_position_all() { let idx = test.iter().position_all(|&v| *v > 2); assert_eq!(idx, vec![0, 3, 4]); - let idx = test.iter().position_all(|&v| *v == 2); - assert_eq!(idx, vec![]); + let idx: Vec = test.iter().position_all(|&v| *v == 2); + assert_eq!(idx, Vec::::new()); } #[test] diff --git a/src/julia/ClarabelRs/src/interface.jl b/src/julia/ClarabelRs/src/interface.jl index 54a1f40..f3fec6d 100644 --- a/src/julia/ClarabelRs/src/interface.jl +++ b/src/julia/ClarabelRs/src/interface.jl @@ -39,6 +39,19 @@ function print_timers(solver::Solver) end +function write_to_file(solver::Solver, filename::String) + + solver_write_to_file_jlrs(solver::Solver, filename::String) + +end + +function read_from_file(filename::String) + + solver_read_from_file_jlrs(filename::String) + +end + + # ------------------------------------- # Wrappers for rust-side interface #-------------------------------------- @@ -77,6 +90,7 @@ function solver_solve_jlrs(solver::Solver) end + function solver_get_info_jlrs(solver::Solver) ccall(Libdl.dlsym(librust,:solver_get_info_jlrs),Clarabel.DefaultInfo{Float64}, @@ -92,6 +106,39 @@ function solver_print_timers_jlrs(solver::Solver) end +function solver_write_to_file_jlrs(solver::Solver, filename::String) + + status = ccall(Libdl.dlsym(librust,:solver_write_to_file_jlrs),Cint, + ( + Ptr{Cvoid}, + Cstring + ), + solver.ptr, + filename, + ) + + if status != 0 + error("Error writing to file $filename") + end + +end + +function solver_read_from_file_jlrs(filename::String) + + ptr = ccall(Libdl.dlsym(librust,:solver_read_from_file_jlrs),Ptr{Cvoid}, + ( + Cstring, + ), + filename, + ) + + if ptr == C_NULL + error("Error reading from file $filename") + end + return Solver{Float64}(ptr) + +end + function solver_drop_jlrs(solver::Solver) ccall(Libdl.dlsym(librust,:solver_drop_jlrs),Cvoid, (Ptr{Cvoid},), solver.ptr) diff --git a/src/julia/interface.rs b/src/julia/interface.rs index e2cde05..1bf2f0e 100644 --- a/src/julia/interface.rs +++ b/src/julia/interface.rs @@ -10,7 +10,11 @@ use crate::solver::{ }; use num_traits::FromPrimitive; use serde_json::*; -use std::{ffi::CStr, os::raw::c_void}; +use std::fs::File; +use std::{ + ffi::CStr, + os::raw::{c_char, c_int, c_void}, +}; // functions for converting solver to / from c void pointers @@ -54,7 +58,7 @@ pub(crate) extern "C" fn solver_new_jlrs( A: &CscMatrixJLRS, b: &VectorJLRS, jlcones: &VectorJLRS, - json_settings: *const std::os::raw::c_char, + json_settings: *const c_char, ) -> *mut c_void { let P = P.to_CscMatrix(); let A = A.to_CscMatrix(); @@ -118,6 +122,69 @@ pub(crate) extern "C" fn solver_print_timers_jlrs(ptr: *mut c_void) { std::mem::forget(solver); } +// dump problem data to a file +// returns -1 on failure, 0 on success +#[no_mangle] +pub(crate) extern "C" fn solver_write_to_file_jlrs( + ptr: *mut c_void, + filename: *const std::os::raw::c_char, +) -> c_int { + let slice = unsafe { CStr::from_ptr(filename) }; + + let filename = match slice.to_str() { + Ok(s) => s, + Err(_) => { + return -1; + } + }; + + let mut file = match File::create(&filename) { + Ok(f) => f, + Err(_) => { + return -1; + } + }; + + let solver = from_ptr(ptr); + let status = solver.write_to_file(&mut file).is_ok(); + let status = if status { 0 } else { -1 } as c_int; + + // don't drop, since the memory is owned by Julia + std::mem::forget(solver); + + return status; +} + +// dump problem data to a file +// returns NULL on failure, pointer to solver on success +#[no_mangle] +pub(crate) extern "C" fn solver_read_from_file_jlrs( + filename: *const std::os::raw::c_char, +) -> *const c_void { + let slice = unsafe { CStr::from_ptr(filename) }; + + let filename = match slice.to_str() { + Ok(s) => s, + Err(_) => { + return std::ptr::null(); + } + }; + + let mut file = match File::open(&filename) { + Ok(f) => f, + Err(_) => { + return std::ptr::null(); + } + }; + + let solver = DefaultSolver::read_from_file(&mut file); + + match solver { + Ok(solver) => to_ptr(Box::new(solver)), + Err(_) => std::ptr::null(), + } +} + // safely drop a solver object through its pointer. // called by the Julia side finalizer when a solver // is out of scope diff --git a/src/python/impl_default_py.rs b/src/python/impl_default_py.rs index 72a13fb..db09381 100644 --- a/src/python/impl_default_py.rs +++ b/src/python/impl_default_py.rs @@ -448,4 +448,17 @@ impl PyDefaultSolver { None => println!("no timers enabled"), }; } + + fn write_to_file(&self, filename: &str) -> PyResult<()> { + let mut file = std::fs::File::create(filename)?; + self.inner.write_to_file(&mut file)?; + Ok(()) + } +} + +#[pyfunction(name = "read_from_file")] +pub fn read_from_file_py(filename: &str) -> PyResult { + let mut file = std::fs::File::open(filename)?; + let solver = DefaultSolver::::read_from_file(&mut file)?; + Ok(PyDefaultSolver { inner: solver }) } diff --git a/src/python/module_py.rs b/src/python/module_py.rs index a0fcc62..339a163 100644 --- a/src/python/module_py.rs +++ b/src/python/module_py.rs @@ -39,6 +39,8 @@ fn clarabel(_py: Python, m: &PyModule) -> PyResult<()> { .unwrap(); m.add_function(wrap_pyfunction!(default_infinity_py, m)?) .unwrap(); + m.add_function(wrap_pyfunction!(read_from_file_py, m)?) + .unwrap(); // API Cone types m.add_class::()?; diff --git a/src/solver/core/cones/supportedcone.rs b/src/solver/core/cones/supportedcone.rs index f95020e..5ea1b39 100644 --- a/src/solver/core/cones/supportedcone.rs +++ b/src/solver/core/cones/supportedcone.rs @@ -2,6 +2,8 @@ use super::*; #[cfg(feature = "sdp")] use crate::algebra::triangular_number; +#[cfg(feature = "serde")] +use serde::{Deserialize, Serialize}; // --------------------------------------------------- // We define some machinery here for enumerating the @@ -9,7 +11,8 @@ use crate::algebra::triangular_number; // --------------------------------------------------- /// API type describing the type of a conic constraint. -/// +/// +#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] #[derive(Debug, Clone)] pub enum SupportedConeT { /// The zero cone (used for equality constraints). diff --git a/src/solver/implementations/default/json.rs b/src/solver/implementations/default/json.rs new file mode 100644 index 0000000..67eecc4 --- /dev/null +++ b/src/solver/implementations/default/json.rs @@ -0,0 +1,135 @@ +use crate::{ + algebra::*, + solver::{DefaultSettings, DefaultSolver, SupportedConeT}, +}; + +use serde::{de::DeserializeOwned, Deserialize, Serialize}; +use std::io::Write; +use std::{fs::File, io, io::Read}; + +// A struct very similar to the problem data, but containing only +// the data types provided by the user (i.e. no internal types). + +#[derive(Serialize, Deserialize)] +#[serde(bound = "T: Serialize + DeserializeOwned")] +struct JsonProblemData { + pub P: CscMatrix, + pub q: Vec, + pub A: CscMatrix, + pub b: Vec, + pub cones: Vec>, + pub settings: DefaultSettings, +} + +impl DefaultSolver +where + T: FloatT + DeserializeOwned + Serialize, +{ + pub fn write_to_file(&self, file: &mut File) -> Result<(), io::Error> { + let mut json_data = JsonProblemData { + P: self.data.P.clone(), + q: self.data.q.clone(), + A: self.data.A.clone(), + b: self.data.b.clone(), + cones: self.data.cones.clone(), + settings: self.settings.clone(), + }; + + // restore scaling to original + let dinv = &self.data.equilibration.dinv; + let einv = &self.data.equilibration.einv; + let c = &self.data.equilibration.c; + + json_data.P.lrscale(dinv, dinv); + json_data.q.hadamard(dinv); + json_data.P.scale(c.recip()); + json_data.q.scale(c.recip()); + + json_data.A.lrscale(einv, dinv); + json_data.b.hadamard(einv); + + // sanitize settings to remove values that + // can't be serialized, i.e. infs + sanitize_settings(&mut json_data.settings); + + // write to file + let json = serde_json::to_string(&json_data)?; + file.write_all(json.as_bytes())?; + + Ok(()) + } + + pub fn read_from_file(file: &mut File) -> Result { + // read file + let mut buffer = String::new(); + file.read_to_string(&mut buffer)?; + let mut json_data: JsonProblemData = serde_json::from_str(&buffer)?; + + // restore sanitized settings to their (likely) original values + desanitize_settings(&mut json_data.settings); + + // create a solver object + let P = json_data.P; + let q = json_data.q; + let A = json_data.A; + let b = json_data.b; + let cones = json_data.cones; + let settings = json_data.settings; + let solver = Self::new(&P, &q, &A, &b, &cones, settings); + + Ok(solver) + } +} + +fn sanitize_settings(settings: &mut DefaultSettings) { + if settings.time_limit == f64::INFINITY { + settings.time_limit = f64::MAX; + } +} + +fn desanitize_settings(settings: &mut DefaultSettings) { + if settings.time_limit == f64::MAX { + settings.time_limit = f64::INFINITY; + } +} + +#[test] +fn test_json_io() { + use crate::solver::IPSolver; + use std::io::{Seek, SeekFrom}; + + let P = CscMatrix { + m: 1, + n: 1, + colptr: vec![0, 1], + rowval: vec![0], + nzval: vec![2.0], + }; + let q = [1.0]; + let A = CscMatrix { + m: 1, + n: 1, + colptr: vec![0, 1], + rowval: vec![0], + nzval: vec![-1.0], + }; + let b = [-2.0]; + let cones = vec![crate::solver::SupportedConeT::NonnegativeConeT(1)]; + + let settings = crate::solver::DefaultSettingsBuilder::default() + .build() + .unwrap(); + + let mut solver = crate::solver::DefaultSolver::::new(&P, &q, &A, &b, &cones, settings); + solver.solve(); + + // write the problem to a file + let mut file = tempfile::tempfile().unwrap(); + solver.write_to_file(&mut file).unwrap(); + + // read the problem from the file + file.seek(SeekFrom::Start(0)).unwrap(); + let mut solver2 = crate::solver::DefaultSolver::::read_from_file(&mut file).unwrap(); + solver2.solve(); + assert_eq!(solver.solution.x, solver2.solution.x); +} diff --git a/src/solver/implementations/default/mod.rs b/src/solver/implementations/default/mod.rs index 1618d5c..f0cea8c 100644 --- a/src/solver/implementations/default/mod.rs +++ b/src/solver/implementations/default/mod.rs @@ -28,3 +28,6 @@ pub use settings::*; pub use solution::*; pub use solver::*; pub use variables::*; + +#[cfg(feature = "serde")] +mod json; diff --git a/src/solver/implementations/default/settings.rs b/src/solver/implementations/default/settings.rs index fc6518b..d7aaee0 100644 --- a/src/solver/implementations/default/settings.rs +++ b/src/solver/implementations/default/settings.rs @@ -2,14 +2,15 @@ use crate::algebra::*; use crate::solver::core::traits::Settings; use derive_builder::Builder; -#[cfg(feature = "julia")] -use serde::{Deserialize, Serialize}; +#[cfg(feature = "serde")] +use serde::{de::DeserializeOwned, Deserialize, Serialize}; /// Standard-form solver type implementing the [`Settings`](crate::solver::core::traits::Settings) trait #[derive(Builder, Debug, Clone)] #[builder(build_fn(validate = "Self::validate"))] -#[cfg_attr(feature = "julia", derive(Serialize, Deserialize))] +#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] +#[serde(bound = "T: Serialize + DeserializeOwned")] pub struct DefaultSettings { ///maximum number of iterations #[builder(default = "200")] From dc0c97fcdbe92a6d008305fc15d5215ed93a7349 Mon Sep 17 00:00:00 2001 From: Paul Goulart Date: Tue, 28 May 2024 19:51:53 +0100 Subject: [PATCH 04/13] supernodal LDL wrapper (#112) --- .github/workflows/ci.yml | 2 +- Cargo.toml | 16 + examples/rust/example_box.rs | 2 +- examples/rust/example_box_faer.rs | 39 ++ src/algebra/floats.rs | 67 ++- src/algebra/tests/vector.rs | 3 + src/algebra/utils.rs | 6 +- src/algebra/vecmath.rs | 4 +- src/qdldl/qdldl.rs | 68 ++- src/qdldl/test.rs | 10 +- .../direct/quasidef/directldlkktsolver.rs | 9 +- .../direct/quasidef/ldlsolvers/faer_ldl.rs | 390 ++++++++++++++++++ .../direct/quasidef/ldlsolvers/mod.rs | 2 + tests/equilibration_bounds.rs | 5 +- 14 files changed, 555 insertions(+), 68 deletions(-) create mode 100644 examples/rust/example_box_faer.rs create mode 100644 src/solver/core/kktsolvers/direct/quasidef/ldlsolvers/faer_ldl.rs diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e5a7285..e2c352a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -28,7 +28,7 @@ jobs: - name: Generate code coverage run: | - cargo tarpaulin --out xml --features sdp-accelerate,serde --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 diff --git a/Cargo.toml b/Cargo.toml index a4f509c..49af1a1 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -56,6 +56,8 @@ julia = ["sdp", "dep:libc", "dep:num-derive", "serde"] python = ["sdp", "dep:libc", "dep:pyo3", "dep:num-derive", "serde"] +#compile with faer supernodal solver option +faer-sparse = ["dep:faer", "dep:faer-entity"] # ------------------------------- # SDP configuration @@ -81,6 +83,14 @@ optional = true version = "2.2" optional = true +[dependencies.faer] +version = "0.19" +optional = true + +[dependencies.faer-entity] +version = "0.19" +optional = true + # ------------------------------- # examples @@ -115,6 +125,11 @@ 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" @@ -122,6 +137,7 @@ required-features = ["serde"] + # ------------------------------- # custom build profiles # ------------------------------- diff --git a/examples/rust/example_box.rs b/examples/rust/example_box.rs index d41b47c..7f9dbad 100644 --- a/examples/rust/example_box.rs +++ b/examples/rust/example_box.rs @@ -3,7 +3,7 @@ use clarabel::algebra::*; use clarabel::solver::*; fn problem_data() -> (CscMatrix, Vec, CscMatrix, Vec) { - let n = 20000; + let n = 2000000; let P = CscMatrix::identity(n); diff --git a/examples/rust/example_box_faer.rs b/examples/rust/example_box_faer.rs new file mode 100644 index 0000000..d0b2588 --- /dev/null +++ b/examples/rust/example_box_faer.rs @@ -0,0 +1,39 @@ +#![allow(non_snake_case)] +use clarabel::algebra::*; +use clarabel::solver::*; + +fn problem_data() -> (CscMatrix, Vec, CscMatrix, Vec) { + let n = 2000000; + + let P = CscMatrix::identity(n); + + // construct A = [I; -I] + let I1 = CscMatrix::::identity(n); + let mut I2 = CscMatrix::::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(); +} diff --git a/src/algebra/floats.rs b/src/algebra/floats.rs index d715992..c337c1c 100644 --- a/src/algebra/floats.rs +++ b/src/algebra/floats.rs @@ -40,53 +40,50 @@ impl 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 MaybeBlasFloatT for T where T: BlasFloatT {} + } + else { + pub trait MaybeBlasFloatT {} + impl 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 FloatT for T where T: CoreFloatT + BlasFloatT {} - } else{ - impl FloatT for T where T: CoreFloatT {} + if #[cfg(feature="faer-sparse")] { + pub trait MaybeFaerFloatT : faer::SimpleEntity + faer::ComplexField {} + impl MaybeFaerFloatT for T where T: faer::SimpleEntity + faer::ComplexField {} + } + else { + pub trait MaybeFaerFloatT {} + impl 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 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 diff --git a/src/algebra/tests/vector.rs b/src/algebra/tests/vector.rs index 4d37433..c8cb4ea 100644 --- a/src/algebra/tests/vector.rs +++ b/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] diff --git a/src/algebra/utils.rs b/src/algebra/utils.rs index 94e63e8..f623a3b 100644 --- a/src/algebra/utils.rs +++ b/src/algebra/utils.rs @@ -6,9 +6,9 @@ // which serves as a vectorized version of the std::iter::position // returning indices of *all* elements satisfying a predicate +use crate::qdldl; use num_traits::Num; use std::cmp::Ordering; -use std::iter::zip; #[cfg_attr(not(sdp), allow(dead_code))] pub(crate) trait PositionAll: Iterator { @@ -34,12 +34,12 @@ where // permutation and inverse permutation pub(crate) fn permute(x: &mut [T], b: &[T], p: &[usize]) { - zip(p, x).for_each(|(p, x)| *x = b[*p]); + qdldl::permute(x, b, p); } #[allow(dead_code)] pub(crate) fn ipermute(x: &mut [T], b: &[T], p: &[usize]) { - zip(p, b).for_each(|(p, b)| x[*p] = *b); + qdldl::ipermute(x, b, p); } // Construct an inverse permutation from a permutation diff --git a/src/algebra/vecmath.rs b/src/algebra/vecmath.rs index e7402c1..84e2b7f 100644 --- a/src/algebra/vecmath.rs +++ b/src/algebra/vecmath.rs @@ -117,11 +117,11 @@ impl VectorMath for [T] { // Returns infinity norm fn norm_inf(&self) -> T { let mut out = T::zero(); - for v in self.iter().map(|v| v.abs()) { + for &v in self { if v.is_nan() { return T::nan(); } - out = if v > out { v } else { out }; + out = T::max(out, v.abs()); } out } diff --git a/src/qdldl/qdldl.rs b/src/qdldl/qdldl.rs index c93c189..053ecaa 100644 --- a/src/qdldl/qdldl.rs +++ b/src/qdldl/qdldl.rs @@ -101,7 +101,7 @@ where // permute b let tmp = &mut self.workspace.fwork; - _permute(tmp, b, &self.perm); + permute(tmp, b, &self.perm); //solve in place with tmp as permuted RHS _solve( @@ -113,7 +113,7 @@ where ); // inverse permutation to put unpermuted soln in b - _ipermute(b, tmp, &self.perm); + ipermute(b, tmp, &self.perm); } pub fn update_values(&mut self, indices: &[usize], values: &[T]) { @@ -195,12 +195,12 @@ fn _qdldl_new( iperm = _invperm(&_perm)?; perm = _perm; } else { - (perm, iperm) = _get_amd_ordering(Ain, opts.amd_dense_scale); + (perm, iperm) = get_amd_ordering(Ain, opts.amd_dense_scale); } //permute to (another) upper triangular matrix and store the //index mapping the input's entries to the permutation's entries - let (A, AtoPAPt) = _permute_symmetric(Ain, &iperm); + let (A, AtoPAPt) = permute_symmetric(Ain, &iperm); // handle the (possibly permuted) vector of // diagonal D signs if one was specified. Otherwise @@ -208,7 +208,7 @@ fn _qdldl_new( let mut Dsigns = vec![1_i8; n]; if let Some(ds) = opts.Dsigns { Dsigns = vec![1_i8; n]; - _permute(&mut Dsigns, &ds, &perm); + permute(&mut Dsigns, &ds, &perm); } // allocate workspace @@ -452,8 +452,8 @@ fn _factor_inner( Lp[0] = 0; let mut acc = 0; for (Lp, Lnz) in zip(&mut Lp[1..], Lnz) { - *Lp = acc + Lnz; - acc = *Lp; + acc += Lnz; + *Lp = acc; } // set all y_idx to be 'unused' initially @@ -681,6 +681,24 @@ fn _ltsolve_unsafe(Lp: &[usize], Li: &[usize], Lx: &[T], x: &mut [T]) } } +// Solves D(L+I)'x = b, with x replacing b. Unchecked version. +fn _dltsolve_unsafe(Lp: &[usize], Li: &[usize], Lx: &[T], Dinv: &[T], x: &mut [T]) { + unsafe { + for i in (0..x.len()).rev() { + let mut s = T::zero(); + let f = *Lp.get_unchecked(i); + let l = *Lp.get_unchecked(i + 1); + for (&Lxj, &Lij) in zip(&Lx[f..l], &Li[f..l]) { + s += Lxj * (*x.get_unchecked(Lij)); + } + + let xi = x.get_unchecked_mut(i); + *xi *= *Dinv.get_unchecked(i); + *xi -= s; + } + } +} + // Solves Ax = b where A has given LDL factors, with x replacing b fn _solve(Lp: &[usize], Li: &[usize], Lx: &[T], Dinv: &[T], b: &mut [T]) { // We call the `unsafe`d version of the forward and backward substitution @@ -688,8 +706,13 @@ fn _solve(Lp: &[usize], Li: &[usize], Lx: &[T], Dinv: &[T], b: &mut [ // compatible dimensions. For super safety or debugging purposes, there // are also `safe` versions implemented above. _lsolve_unsafe(Lp, Li, Lx, b); - zip(b.iter_mut(), Dinv).for_each(|(b, d)| *b *= *d); - _ltsolve_unsafe(Lp, Li, Lx, b); + + // combined D and L^T solve in one pass + _dltsolve_unsafe(Lp, Li, Lx, Dinv, b); + + // in separate passes for reference + //zip(b.iter_mut(), Dinv).for_each(|(b, d)| *b *= *d); + //_ltsolve_unsafe(Lp, Li, Lx, b); } // Construct an inverse permutation from a permutation @@ -706,21 +729,32 @@ fn _invperm(p: &[usize]) -> Result, QDLDLError> { Ok(b) } -// internal permutation and inverse permutation -// functions that require no memory allocations +// permutation and inverse permutation +// functions that require no allocation +// p must be a valid permutation vector +// in both cases for safety -fn _permute(x: &mut [T], b: &[T], p: &[usize]) { - zip(p, x).for_each(|(p, x)| *x = b[*p]); +pub(crate) fn permute(x: &mut [T], b: &[T], p: &[usize]) { + debug_assert!(*p.iter().max().unwrap_or(&0) < x.len()); + unsafe { + zip(p, x).for_each(|(p, x)| *x = *b.get_unchecked(*p)); + } } -fn _ipermute(x: &mut [T], b: &[T], p: &[usize]) { - zip(p, b).for_each(|(p, b)| x[*p] = *b); +pub(crate) fn ipermute(x: &mut [T], b: &[T], p: &[usize]) { + debug_assert!(*p.iter().max().unwrap_or(&0) < x.len()); + unsafe { + zip(p, b).for_each(|(p, b)| *x.get_unchecked_mut(*p) = *b); + } } // Given a sparse symmetric matrix `A` (with only upper triangular entries), return // permuted sparse symmetric matrix `P` (also only upper triangular) given the // inverse permutation vector `iperm`." -fn _permute_symmetric(A: &CscMatrix, iperm: &[usize]) -> (CscMatrix, Vec) { +pub(crate) fn permute_symmetric( + A: &CscMatrix, + iperm: &[usize], +) -> (CscMatrix, Vec) { // perform a number of argument checks let (_m, n) = A.size(); let mut P = CscMatrix::::spalloc((n, n), A.nnz()); @@ -816,7 +850,7 @@ fn _permute_symmetric_inner( } } -fn _get_amd_ordering( +pub(crate) fn get_amd_ordering( A: &CscMatrix, amd_dense_scale: f64, ) -> (Vec, Vec) { diff --git a/src/qdldl/test.rs b/src/qdldl/test.rs index 7b6aafe..c7d2ebc 100644 --- a/src/qdldl/test.rs +++ b/src/qdldl/test.rs @@ -54,10 +54,10 @@ fn test_permute() { let mut x = vec![0.; 4]; let mut y = vec![0.; 4]; - _permute(&mut x, &b, &perm); + permute(&mut x, &b, &perm); assert_eq!(x, vec![4., 1., 3., 2.]); - _ipermute(&mut y, &x, &perm); + ipermute(&mut y, &x, &perm); assert_eq!(y, b); } @@ -124,7 +124,7 @@ fn test_etree() { #[test] fn test_amd() { let A = test_matrix_4x4(); - let (perm, iperm) = _get_amd_ordering(&A, 1.5); + let (perm, iperm) = get_amd_ordering(&A, 1.5); assert_eq!(perm, [3, 0, 1, 2]); assert_eq!(iperm, [1, 2, 3, 0]); } @@ -134,7 +134,7 @@ fn test_permute_symmetric() { //no permutation at all let A = test_matrix_4x4(); let iperm: Vec = vec![0, 1, 2, 3]; - let (P, AtoPAPt) = _permute_symmetric(&A, &iperm); + let (P, AtoPAPt) = permute_symmetric(&A, &iperm); assert_eq!(&A.colptr, &P.colptr); assert_eq!(&A.rowval, &P.rowval); @@ -157,7 +157,7 @@ fn test_permute_symmetric() { let perm: Vec = vec![2, 3, 0, 1]; let iperm = _invperm(&perm).unwrap(); - let (P, _) = _permute_symmetric(&A, &iperm); + let (P, _) = permute_symmetric(&A, &iperm); assert_eq!(&P.colptr, &vec![0, 1, 3, 5, 8]); assert_eq!(&P.rowval, &vec![0, 0, 1, 2, 0, 2, 3, 0]); diff --git a/src/solver/core/kktsolvers/direct/quasidef/directldlkktsolver.rs b/src/solver/core/kktsolvers/direct/quasidef/directldlkktsolver.rs index 8ee0687..9773d3b 100644 --- a/src/solver/core/kktsolvers/direct/quasidef/directldlkktsolver.rs +++ b/src/solver/core/kktsolvers/direct/quasidef/directldlkktsolver.rs @@ -1,5 +1,8 @@ #![allow(non_snake_case)] +#[cfg(feature = "faer-sparse")] +use super::ldlsolvers::faer_ldl::*; + use super::ldlsolvers::qdldl::*; use super::*; use crate::solver::core::kktsolvers::KKTSolver; @@ -345,8 +348,10 @@ where kktshape = QDLDLDirectLDLSolver::::required_matrix_shape(); ldlptr = |M, D, S| Box::new(QDLDLDirectLDLSolver::::new(M, D, S)); } - "custom" => { - unimplemented!(); + #[cfg(feature = "faer-sparse")] + "faer" => { + kktshape = FaerDirectLDLSolver::::required_matrix_shape(); + ldlptr = |M, D, S| Box::new(FaerDirectLDLSolver::::new(M, D, S)); } _ => { panic! {"Unrecognized LDL solver type"}; diff --git a/src/solver/core/kktsolvers/direct/quasidef/ldlsolvers/faer_ldl.rs b/src/solver/core/kktsolvers/direct/quasidef/ldlsolvers/faer_ldl.rs new file mode 100644 index 0000000..980ec55 --- /dev/null +++ b/src/solver/core/kktsolvers/direct/quasidef/ldlsolvers/faer_ldl.rs @@ -0,0 +1,390 @@ +#![allow(non_snake_case)] + +use faer::{ + dyn_stack::{GlobalPodBuffer, PodStack, StackReq}, + linalg::cholesky::ldlt_diagonal::compute::LdltRegularization, + sparse::{ + linalg::{amd::Control, cholesky::*, SupernodalThreshold}, + SparseColMatRef, SymbolicSparseColMatRef, + }, + Conj, Parallelism, Side, +}; + +use crate::algebra::*; +use crate::solver::core::kktsolvers::direct::DirectLDLSolver; +use crate::solver::core::CoreSettings; +use std::iter::zip; + +#[derive(Debug)] + +struct FaerLDLRegularizerParams { + regularize: bool, + eps: T, + delta: T, +} + +impl<'a, T> FaerLDLRegularizerParams +where + T: FloatT, +{ + fn to_faer(&self, Dsigns: &'a [i8]) -> LdltRegularization<'a, T> { + let option = if self.regularize { Some(Dsigns) } else { None }; + LdltRegularization { + dynamic_regularization_signs: option, + dynamic_regularization_delta: self.delta, + dynamic_regularization_epsilon: self.eps, + } + } +} + +pub struct FaerDirectLDLSolver { + // we will pre-permute the KKT matrix and factor it inside faer with + // an identity permutation. This means we will also need to store + // the mapping from the original ordering to the permuted ordering + // for both the KKT matrix and the Dsigns, and will need to permute + // the rhs and solution vectors accordingly + perm: Vec, + iperm: Vec, + + // permuted KKT matrix + perm_kkt: CscMatrix, + perm_map: Vec, + perm_dsigns: Vec, + + // space for permuted LHS/RHS when solving + bperm: Vec, + + // regularizer parameters captured from settings + regularizer_params: FaerLDLRegularizerParams, + + // symbolic + numeric cholesky data + symbolic_cholesky: SymbolicCholesky, + ld_vals: Vec, + + // workspace for faer factor/solve calls + work: GlobalPodBuffer, + + parallelism: Parallelism<'static>, +} + +impl FaerDirectLDLSolver +where + T: FloatT, +{ + pub fn new(KKT: &CscMatrix, Dsigns: &[i8], settings: &CoreSettings) -> Self { + assert!(KKT.is_square(), "KKT matrix is not square"); + + // ----------------------------- + + // Rayon(0) here is equivalent to rayon::current_num_threads() + // PJG: this should be a user-settable option + let parallelism = Parallelism::Rayon(0); + + // manually compute an AMD ordering for the KKT matrix + // and permute it to match the ordering used in QDLDL + let amd_dense_scale = 1.5; // magic number from QDLDL + let (perm, iperm) = crate::qdldl::get_amd_ordering(KKT, amd_dense_scale); + + let (mut perm_kkt, mut perm_map) = crate::qdldl::permute_symmetric(KKT, &iperm); + + //Permute the Dsigns to match the ordering to be used internally + let mut perm_dsigns = vec![1_i8; Dsigns.len()]; + permute(&mut perm_dsigns, Dsigns, &perm); + + // the amd_params will not be used by faer even though + // we put them into CholeskySymbolicParams, provided + // that my ordering is SymmetricOrdering::Identity + + let amd_params = Control { + ..Default::default() + }; + + let supernodal_flop_ratio_threshold = SupernodalThreshold::AUTO; + let cholesky_params = CholeskySymbolicParams { + supernodal_flop_ratio_threshold, + amd_params, + ..Default::default() + }; + // ----------------------------- + + sort_csc_columns_with_map(&mut perm_kkt, &mut perm_map); + + let symbKKT = SymbolicSparseColMatRef::new_checked( + perm_kkt.n, + perm_kkt.n, + &perm_kkt.colptr, + None, + &perm_kkt.rowval, + ); + + let symbolic_cholesky = factorize_symbolic_cholesky( + symbKKT, + Side::Upper, + SymmetricOrdering::Identity, + cholesky_params, + ) + .unwrap(); + + let ld_vals = vec![T::zero(); symbolic_cholesky.len_values()]; + + let regularizer_params = FaerLDLRegularizerParams { + regularize: settings.dynamic_regularization_enable, + delta: settings.dynamic_regularization_delta, + eps: settings.dynamic_regularization_eps, + }; + + // Required workspace for faer factor and solve + let req_factor = symbolic_cholesky + .factorize_numeric_ldlt_req::(true, parallelism) + .unwrap(); + let req_solve = symbolic_cholesky.solve_in_place_req::(1).unwrap(); + let req = StackReq::any_of([req_factor, req_solve]); + let work = GlobalPodBuffer::new(req); + + let bperm = vec![T::zero(); perm_kkt.n]; + + Self { + perm, + iperm, + perm_kkt, + perm_map, + perm_dsigns, + bperm, + regularizer_params, + symbolic_cholesky, + ld_vals, + work, + parallelism, + } + } +} + +impl DirectLDLSolver for FaerDirectLDLSolver +where + T: FloatT + faer::Entity + faer::SimpleEntity, +{ + fn update_values(&mut self, index: &[usize], values: &[T]) { + // PJG: this is replicating the update_values function in qdldl + let nzval = &mut self.perm_kkt.nzval; // post perm internal data + let AtoPAPt = &self.perm_map; //mapping from input matrix entries + + for (i, &idx) in index.iter().enumerate() { + nzval[AtoPAPt[idx]] = values[i]; + } + } + + fn scale_values(&mut self, index: &[usize], scale: T) { + // PJG: this is replicating the scale_values function in qdldl + let nzval = &mut self.perm_kkt.nzval; // post perm internal data + let AtoPAPt = &self.perm_map; //mapping from input matrix entries + + for &idx in index.iter() { + nzval[AtoPAPt[idx]] *= scale; + } + } + + fn offset_values(&mut self, index: &[usize], offset: T, signs: &[i8]) { + // PJG: this is replicating the offset_values function in qdldl + let nzval = &mut self.perm_kkt.nzval; // post perm internal data + let AtoPAPt = &self.perm_map; //mapping from input matrix entries + + for (&idx, &sign) in zip(index, signs) { + let sign: T = T::from_i8(sign).unwrap(); + nzval[AtoPAPt[idx]] += offset * sign; + } + } + + fn solve(&mut self, _kkt: &CscMatrix, x: &mut [T], b: &[T]) { + // NB: faer solves in place. Permute b to match the ordering used internally + permute(&mut self.bperm, b, &self.perm); + + let rhs = faer::mat::from_column_major_slice_mut::(&mut self.bperm[0..], b.len(), 1); + let ldlt = LdltRef::new(&self.symbolic_cholesky, self.ld_vals.as_slice()); + + ldlt.solve_in_place_with_conj( + Conj::No, + rhs, + self.parallelism, + PodStack::new(&mut self.work), + ); + + // bperm is now the solution, permute it back to the original ordering + permute(x, &self.bperm, &self.iperm); + } + + fn refactor(&mut self, _kkt: &CscMatrix) -> bool { + let symbKKT = SymbolicSparseColMatRef::new_checked( + self.perm_kkt.n, + self.perm_kkt.n, + &self.perm_kkt.colptr, + None, + &self.perm_kkt.rowval, + ); + + let a: SparseColMatRef = + SparseColMatRef::new(symbKKT, self.perm_kkt.nzval.as_slice()); + + let regularizer = self.regularizer_params.to_faer(&self.perm_dsigns); + + self.symbolic_cholesky.factorize_numeric_ldlt( + self.ld_vals.as_mut_slice(), + a, + Side::Upper, + regularizer, + self.parallelism, + PodStack::new(&mut self.work), + ); + + true // assume success for experimental purposes + } + + fn required_matrix_shape() -> MatrixTriangle { + MatrixTriangle::Triu + } +} + +// --------------------------------------------------------------------- +// utility functions +// --------------------------------------------------------------------- + +fn sort_csc_columns_with_map(M: &mut CscMatrix, map: &mut [usize]) +where + T: FloatT, +{ + // when we permute the KKT matrix into an upper triangular form, the + // data within each column is no long ordered by increasing row number. + + // Fix that here, and update the map to reflect the new ordering + + // we need to sort on the inverted map, then invert back + let mut imap = invperm(map); + + // it is not obvious to me how to sort all three of the vectors + // M.nzvals, Mrowvals, and map together without allocating. + // Allocate a permutation vector for the whole thing, get sort + // indices for the rowvals columns, then permute everything + let mut perm = (0..M.nnz()).collect::>(); + + for col in 0..M.n { + let start = M.colptr[col]; + let end = M.colptr[col + 1]; + perm[start..end].sort_by_key(|&k| M.rowval[k]); + } + + let mut tmpT = vec![T::zero(); M.nzval.len()]; + let mut tmpint = vec![0usize; M.nzval.len()]; + + tmpT.copy_from_slice(&M.nzval); + permute(&mut M.nzval, &tmpT, &perm); + + tmpint.copy_from_slice(&M.rowval); + permute(&mut M.rowval, &tmpint, &perm); + + tmpint.copy_from_slice(&imap); + permute(&mut imap, &tmpint, &perm); + + // map = invperm(imap), but in place + for (i, j) in imap.iter().enumerate() { + map[*j] = i; + } + + M.check_format().unwrap(); +} + +// --------------------------------------------------------------------- +// tests +// --------------------------------------------------------------------- + +#[test] +fn test_faer_ldl() { + let KKT = CscMatrix { + m: 6, + n: 6, + colptr: vec![0, 1, 2, 4, 6, 8, 10], + rowval: vec![0, 1, 0, 2, 1, 3, 0, 4, 1, 5], + nzval: vec![1.0, 2.0, 1.0, -1.0, 1.0, -2.0, -1.0, -3.0, -1.0, -4.0], + }; + + let Dsigns = vec![1, 1, -1, -1, -1, -1]; + + let mut solver = FaerDirectLDLSolver::new(&KKT, &Dsigns, &CoreSettings::default()); + + let mut x = vec![0.0; 6]; + let b = vec![1.0, 2.0, 3.0, 4., 5., 6.]; + + // check that the permuted values are in what should be natural order + let nzval = &solver.perm_kkt.nzval; // post perm internal data + let AtoPAPt = &solver.perm_map; //mapping from input matrix entries + + for i in 0..nzval.len() { + assert_eq!(KKT.nzval[i], nzval[AtoPAPt[i]]); + } + + solver.refactor(&KKT); + solver.solve(&KKT, &mut x, &b); + + let xsol = vec![ + 1.0, + 0.9090909090909091, + -2.0, + -1.5454545454545454, + -2.0, + -1.7272727272727275, + ]; + assert!(x.norm_inf_diff(&xsol) < 1e-10); + + // assign a new entry at the end of the KKT matrix and resolve + solver.update_values(&[9], &[-10.0]); + + solver.refactor(&KKT); + solver.solve(&KKT, &mut x, &b); + + let xsol = vec![ + 1.0, + 1.3076923076923077, + -2.0, + -1.346153846153846, + -2.0, + -0.7307692307692306, + ]; + assert!(x.norm_inf_diff(&xsol) < 1e-10); + + // scale and update everything for codecov. + solver.offset_values(&[1, 2], 3., &[1, -1]); + solver.scale_values(&[1, 2], 2.); +} + +#[test] +fn test_faer_qp() { + use crate::solver::core::solver::IPSolver; + use crate::solver::SolverStatus; + + let P = CscMatrix { + m: 1, + n: 1, + colptr: vec![0, 1], + rowval: vec![0], + nzval: vec![2.0], + }; + let q = [1.0]; + let A = CscMatrix { + m: 1, + n: 1, + colptr: vec![0, 1], + rowval: vec![0], + nzval: vec![-1.0], + }; + let b = [-2.0]; + let cones = vec![crate::solver::SupportedConeT::NonnegativeConeT(1)]; + + let settings = crate::solver::DefaultSettingsBuilder::default() + .direct_solve_method("faer".to_owned()) + .build() + .unwrap(); + + let mut solver = crate::solver::DefaultSolver::new(&P, &q, &A, &b, &cones, settings); + + solver.solve(); + + assert_eq!(solver.solution.status, SolverStatus::Solved); +} diff --git a/src/solver/core/kktsolvers/direct/quasidef/ldlsolvers/mod.rs b/src/solver/core/kktsolvers/direct/quasidef/ldlsolvers/mod.rs index 7c05652..93ac154 100644 --- a/src/solver/core/kktsolvers/direct/quasidef/ldlsolvers/mod.rs +++ b/src/solver/core/kktsolvers/direct/quasidef/ldlsolvers/mod.rs @@ -1 +1,3 @@ +#[cfg(feature = "faer-sparse")] +pub mod faer_ldl; pub mod qdldl; diff --git a/tests/equilibration_bounds.rs b/tests/equilibration_bounds.rs index ec918bb..7f143a0 100644 --- a/tests/equilibration_bounds.rs +++ b/tests/equilibration_bounds.rs @@ -69,8 +69,6 @@ fn test_equilibrate_upper_bound() { let mut solver = DefaultSolver::new(&P, &c, &A, &b, &cones, settings.clone()); - solver.solve(); - let d = &solver.data.equilibration.d; let e = &solver.data.equilibration.e; @@ -78,6 +76,9 @@ fn test_equilibrate_upper_bound() { assert!(e.minimum() >= settings.equilibrate_min_scaling); assert!(d.maximum() <= settings.equilibrate_max_scaling); assert!(e.maximum() <= settings.equilibrate_max_scaling); + + // forces poorly converging test for codecov + solver.solve(); } #[test] From 7d3b2ffadb13d53d84b04105eadf697f80305c76 Mon Sep 17 00:00:00 2001 From: goulart-paul Date: Tue, 28 May 2024 19:54:08 +0100 Subject: [PATCH 05/13] bump version --- Cargo.toml | 2 +- README.md | 2 +- src/julia/ClarabelRs/Project.toml | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 49af1a1..e70e5d4 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "clarabel" -version = "0.9" +version = "0.9.0" authors = ["Paul Goulart "] edition = "2021" rust-version = "1.66" diff --git a/README.md b/README.md index 9ba3270..8bb417a 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ Interior Point Conic Optimization for Rust and Python - +

diff --git a/src/julia/ClarabelRs/Project.toml b/src/julia/ClarabelRs/Project.toml index c42eaf4..a4a0c7c 100644 --- a/src/julia/ClarabelRs/Project.toml +++ b/src/julia/ClarabelRs/Project.toml @@ -1,7 +1,7 @@ name = "ClarabelRs" uuid = "a0c58a0a-712c-48b7-9fd4-64369ecb2011" authors = ["Paul Goulart "] -version = "0.9" +version = "0.9.0" [deps] Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" @@ -11,4 +11,4 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] -Clarabel = "0.9" +Clarabel = "0.9.0" From 78d5768c4e01431af9cc589307551e5405f048de Mon Sep 17 00:00:00 2001 From: goulart-paul Date: Tue, 28 May 2024 19:57:08 +0100 Subject: [PATCH 06/13] fix unit test --- src/solver/implementations/default/settings.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solver/implementations/default/settings.rs b/src/solver/implementations/default/settings.rs index d7aaee0..7f3ed7e 100644 --- a/src/solver/implementations/default/settings.rs +++ b/src/solver/implementations/default/settings.rs @@ -301,7 +301,7 @@ fn test_settings_validate() { .build(); cfg_if::cfg_if! { if #[cfg(feature = "faer-sparse")] { - assert!(build.is_ok()); + assert!(builder.is_ok()); } else { assert!(builder.is_err()); From e572999c9cfca74a02082ac8bbeeb8a9742f81f9 Mon Sep 17 00:00:00 2001 From: goulart-paul Date: Tue, 28 May 2024 20:23:09 +0100 Subject: [PATCH 07/13] update julia/python opts --- Cargo.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index e70e5d4..676e789 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -47,13 +47,13 @@ sdp-r = ["sdp", "blas-src/r", "lapack-src/r"] # build as the julia interface -julia = ["sdp", "dep:libc", "dep:num-derive", "serde"] +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", "serde"] +python = ["sdp", "dep:libc", "dep:pyo3", "dep:num-derive", "serde", "faer-sparse"] #compile with faer supernodal solver option From eac9d5c4a972c5b352f090e2fc07c425e1c888cc Mon Sep 17 00:00:00 2001 From: goulart-paul Date: Wed, 29 May 2024 09:46:46 +0100 Subject: [PATCH 08/13] skip regularization for symbolic factor --- src/qdldl/qdldl.rs | 38 ++++++++++++++++++++------------------ 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/src/qdldl/qdldl.rs b/src/qdldl/qdldl.rs index 053ecaa..011a204 100644 --- a/src/qdldl/qdldl.rs +++ b/src/qdldl/qdldl.rs @@ -590,27 +590,29 @@ fn _factor_inner( y_markers[cidx] = QDLDL_UNUSED; } - // apply dynamic regularization - if regularize_enable { - let sign = T::from_i8(Dsigns[k]).unwrap(); - if D[k] * sign < regularize_eps { - D[k] = regularize_delta * sign; - *regularize_count += 1; + if !logical_factor { + // apply dynamic regularization + if regularize_enable { + let sign = T::from_i8(Dsigns[k]).unwrap(); + if D[k] * sign < regularize_eps { + D[k] = regularize_delta * sign; + *regularize_count += 1; + } } - } - // Maintain a count of the positive entries - // in D. If we hit a zero, we can't factor - // this matrix, so abort - if D[k] == T::zero() { - return Err(QDLDLError::ZeroPivot); - } - if D[k] > T::zero() { - positiveValuesInD += 1; - } + // Maintain a count of the positive entries + // in D. If we hit a zero, we can't factor + // this matrix, so abort + if D[k] == T::zero() { + return Err(QDLDLError::ZeroPivot); + } + if D[k] > T::zero() { + positiveValuesInD += 1; + } - // compute the inverse of the diagonal - Dinv[k] = T::recip(D[k]); + // compute the inverse of the diagonal + Dinv[k] = T::recip(D[k]); + } } //end for k Ok(positiveValuesInD) From eb6c5d9a69af522628d9a0ce1969b2eb1fc261a1 Mon Sep 17 00:00:00 2001 From: goulart-paul Date: Wed, 29 May 2024 13:13:09 +0100 Subject: [PATCH 09/13] enlarge unsafe block in factor --- src/qdldl/qdldl.rs | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/qdldl/qdldl.rs b/src/qdldl/qdldl.rs index 011a204..d3fba18 100644 --- a/src/qdldl/qdldl.rs +++ b/src/qdldl/qdldl.rs @@ -571,13 +571,14 @@ fn _factor_inner( for (&Lxj, &Lij) in zip(&Lx[f..l], &Li[f..l]) { *(y_vals.get_unchecked_mut(Lij)) -= Lxj * y_vals_cidx; } - } - // Now I have the cidx^th element of y = L\b. - // so compute the corresponding element of - // this row of L and put it into the right place - Lx[tmp_idx] = y_vals_cidx * Dinv[cidx]; - D[k] -= y_vals_cidx * Lx[tmp_idx]; + // Now I have the cidx^th element of y = L\b. + // so compute the corresponding element of + // this row of L and put it into the right place + let Lx_tmp_idx = y_vals_cidx * *Dinv.get_unchecked(cidx); + *Lx.get_unchecked_mut(tmp_idx) = Lx_tmp_idx; + *D.get_unchecked_mut(k) -= y_vals_cidx * Lx_tmp_idx; + } } // record which row it went into From 7069d7831462c34773ba56d447be969e9435d1ad Mon Sep 17 00:00:00 2001 From: Alex Rice Date: Fri, 31 May 2024 21:15:24 +0100 Subject: [PATCH 10/13] Add wasm feature (#114) * Add wasm feature std::time panics on wasm, so we introduce a feature flag that swaps it out with the web-time crate. --- Cargo.toml | 12 +++++++++++- src/timers/timers.rs | 11 ++++++++++- 2 files changed, 21 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 676e789..741628c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -55,6 +55,7 @@ julia = ["sdp", "dep:libc", "dep:num-derive", "serde", "faer-sparse"] # enable a blas/lapack source package 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"] @@ -190,8 +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" \ No newline at end of file +tempfile = "3" + diff --git a/src/timers/timers.rs b/src/timers/timers.rs index 671111c..741d03f 100644 --- a/src/timers/timers.rs +++ b/src/timers/timers.rs @@ -1,6 +1,15 @@ use std::collections::HashMap; use std::ops::{Deref, DerefMut}; -use std::time::{Duration, Instant}; + +cfg_if::cfg_if! { + if #[cfg(feature="wasm")] { + use web_time::{Duration, Instant}; + } + else { + use std::time::{Duration, Instant}; + } +} + #[derive(Debug, Default)] struct InnerTimer { From bb84a5ebf39d5f0943daec9b14399c3a8b120291 Mon Sep 17 00:00:00 2001 From: goulart-paul Date: Sat, 1 Jun 2024 16:27:51 +0100 Subject: [PATCH 11/13] update for docs --- src/solver/core/solver.rs | 8 ++++++++ src/solver/implementations/default/json.rs | 8 ++++---- src/solver/implementations/default/problemdata.rs | 2 +- src/solver/mod.rs | 6 +++++- 4 files changed, 18 insertions(+), 6 deletions(-) diff --git a/src/solver/core/solver.rs b/src/solver/core/solver.rs index d60028c..3c5ae7e 100644 --- a/src/solver/core/solver.rs +++ b/src/solver/core/solver.rs @@ -90,6 +90,14 @@ impl std::fmt::Display for SolverStatus { } } +/// JSON file read/write trait for solver data. +/// Only available with the "serde" feature enabled. +#[cfg(feature = "serde")] +pub trait SolverJSONReadWrite: Sized { + fn write_to_file(&self, file: &mut std::fs::File) -> Result<(), std::io::Error>; + fn read_from_file(file: &mut std::fs::File) -> Result; +} + // --------------------------------- // top level solver container type // --------------------------------- diff --git a/src/solver/implementations/default/json.rs b/src/solver/implementations/default/json.rs index 67eecc4..a01147e 100644 --- a/src/solver/implementations/default/json.rs +++ b/src/solver/implementations/default/json.rs @@ -1,6 +1,6 @@ use crate::{ algebra::*, - solver::{DefaultSettings, DefaultSolver, SupportedConeT}, + solver::{core::SolverJSONReadWrite, DefaultSettings, DefaultSolver, SupportedConeT}, }; use serde::{de::DeserializeOwned, Deserialize, Serialize}; @@ -21,11 +21,11 @@ struct JsonProblemData { pub settings: DefaultSettings, } -impl DefaultSolver +impl SolverJSONReadWrite for DefaultSolver where T: FloatT + DeserializeOwned + Serialize, { - pub fn write_to_file(&self, file: &mut File) -> Result<(), io::Error> { + fn write_to_file(&self, file: &mut File) -> Result<(), io::Error> { let mut json_data = JsonProblemData { P: self.data.P.clone(), q: self.data.q.clone(), @@ -59,7 +59,7 @@ where Ok(()) } - pub fn read_from_file(file: &mut File) -> Result { + fn read_from_file(file: &mut File) -> Result { // read file let mut buffer = String::new(); file.read_to_string(&mut buffer)?; diff --git a/src/solver/implementations/default/problemdata.rs b/src/solver/implementations/default/problemdata.rs index 043ed1d..d071592 100644 --- a/src/solver/implementations/default/problemdata.rs +++ b/src/solver/implementations/default/problemdata.rs @@ -359,7 +359,7 @@ where // an alternative. Necessary since the Options for q and b are &Vec, but // the user supplied data is a slice &[T] -pub fn unwrap_and_slice_or_else<'a, T, F>(opt: &'a Option>, f: F) -> &'a [T] +pub(crate) fn unwrap_and_slice_or_else<'a, T, F>(opt: &'a Option>, f: F) -> &'a [T] where F: FnOnce() -> &'a [T], T: FloatT, diff --git a/src/solver/mod.rs b/src/solver/mod.rs index ee216ba..1259001 100644 --- a/src/solver/mod.rs +++ b/src/solver/mod.rs @@ -38,7 +38,11 @@ pub use crate::solver::core::{IPSolver, SolverStatus}; pub use crate::solver::core::traits; pub use crate::solver::core::CoreSettings; -//If we had implementations for multple alternative +// read/write types if enabled +#[cfg(feature = "serde")] +pub use crate::solver::core::SolverJSONReadWrite; + +//If we had implementations for multiple alternative //problem formats, they would live here. Since we //only have default, it is exposed at the top level //in the use statements directly below instead. From 424bd10a6113fd109a45fca59b2a12fb81b49a71 Mon Sep 17 00:00:00 2001 From: goulart-paul Date: Sat, 1 Jun 2024 16:49:43 +0100 Subject: [PATCH 12/13] resolve lints for rustc >= 1.80 --- src/algebra/csc/core.rs | 2 +- src/algebra/matrix_traits.rs | 2 +- src/algebra/scalarmath.rs | 8 ++++---- src/algebra/utils.rs | 8 ++++---- src/solver/core/cones/supportedcone.rs | 8 ++++---- src/solver/implementations/default/problemdata.rs | 2 +- src/solver/implementations/default/variables.rs | 2 +- 7 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/algebra/csc/core.rs b/src/algebra/csc/core.rs index ee87d28..72ef724 100644 --- a/src/algebra/csc/core.rs +++ b/src/algebra/csc/core.rs @@ -274,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, Vec, Vec) { let I = self.rowval.clone(); let mut J = Vec::with_capacity(self.nnz()); diff --git a/src/algebra/matrix_traits.rs b/src/algebra/matrix_traits.rs index 4224437..c00bb71 100644 --- a/src/algebra/matrix_traits.rs +++ b/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); diff --git a/src/algebra/scalarmath.rs b/src/algebra/scalarmath.rs index 932712f..57d64af 100644 --- a/src/algebra/scalarmath.rs +++ b/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 } diff --git a/src/algebra/utils.rs b/src/algebra/utils.rs index f623a3b..0a5c8f5 100644 --- a/src/algebra/utils.rs +++ b/src/algebra/utils.rs @@ -10,7 +10,7 @@ use crate::qdldl; use num_traits::Num; use std::cmp::Ordering; -#[cfg_attr(not(sdp), allow(dead_code))] +#[cfg_attr(not(feature = "sdp"), allow(dead_code))] pub(crate) trait PositionAll: Iterator { fn position_all(&mut self, predicate: F) -> Vec where @@ -43,7 +43,7 @@ pub(crate) fn ipermute(x: &mut [T], b: &[T], p: &[usize]) { } // Construct an inverse permutation from a permutation -#[cfg_attr(not(sdp), allow(dead_code))] +#[cfg_attr(not(feature = "sdp"), allow(dead_code))] pub(crate) fn invperm(p: &[usize]) -> Vec { let mut b = vec![0; p.len()]; for (i, j) in p.iter().enumerate() { @@ -63,7 +63,7 @@ where p.sort_by_key(|&k| v[k]); } -#[cfg_attr(not(sdp), allow(dead_code))] +#[cfg_attr(not(feature = "sdp"), allow(dead_code))] pub(crate) fn sortperm_rev(p: &mut [usize], v: &[T]) where T: Sized + Ord + Copy, @@ -86,7 +86,7 @@ where // non-float types (e.g. usize). Would require partition of the // vector math traits into those that require FloatT and those // that only require Num + Ord. -#[cfg_attr(not(sdp), allow(dead_code))] +#[cfg_attr(not(feature = "sdp"), allow(dead_code))] pub(crate) fn findmax(v: &[T]) -> Option where T: Num + Copy + Ord, diff --git a/src/solver/core/cones/supportedcone.rs b/src/solver/core/cones/supportedcone.rs index 5ea1b39..52644ef 100644 --- a/src/solver/core/cones/supportedcone.rs +++ b/src/solver/core/cones/supportedcone.rs @@ -203,14 +203,14 @@ impl SupportedConeTag { //PJG: type names are not satisfactory. Try to combine //with the internal cone generators. -#[cfg_attr(not(sdp), allow(dead_code))] +#[cfg_attr(not(feature = "sdp"), allow(dead_code))] pub(crate) struct RangeSupportedConesIterator<'a, T> { cones: &'a [SupportedConeT], index: usize, start: usize, } -#[cfg_attr(not(sdp), allow(dead_code))] +#[cfg_attr(not(feature = "sdp"), allow(dead_code))] impl<'a, T> Iterator for RangeSupportedConesIterator<'a, T> { type Item = std::ops::Range; @@ -227,12 +227,12 @@ impl<'a, T> Iterator for RangeSupportedConesIterator<'a, T> { } } } -#[cfg_attr(not(sdp), allow(dead_code))] +#[cfg_attr(not(feature = "sdp"), allow(dead_code))] pub(crate) trait ConeRanges<'a, T> { fn rng_cones_iter(&'a self) -> RangeSupportedConesIterator<'a, T>; } -#[cfg_attr(not(sdp), allow(dead_code))] +#[cfg_attr(not(feature = "sdp"), allow(dead_code))] impl<'a, T> ConeRanges<'a, T> for [SupportedConeT] { fn rng_cones_iter(&'a self) -> RangeSupportedConesIterator<'a, T> { RangeSupportedConesIterator::<'a, T> { diff --git a/src/solver/implementations/default/problemdata.rs b/src/solver/implementations/default/problemdata.rs index d071592..c22f3ba 100644 --- a/src/solver/implementations/default/problemdata.rs +++ b/src/solver/implementations/default/problemdata.rs @@ -358,7 +358,7 @@ where // -- utility function that tries to unwrap and slice a vector, or return // an alternative. Necessary since the Options for q and b are &Vec, but // the user supplied data is a slice &[T] - +#[cfg(feature = "sdp")] pub(crate) fn unwrap_and_slice_or_else<'a, T, F>(opt: &'a Option>, f: F) -> &'a [T] where F: FnOnce() -> &'a [T], diff --git a/src/solver/implementations/default/variables.rs b/src/solver/implementations/default/variables.rs index dc58106..178e9fd 100644 --- a/src/solver/implementations/default/variables.rs +++ b/src/solver/implementations/default/variables.rs @@ -283,7 +283,7 @@ where self.κ *= scaleinv; } - #[cfg_attr(not(sdp), allow(dead_code))] + #[cfg_attr(not(feature = "sdp"), allow(dead_code))] pub(crate) fn dims(&self) -> (usize, usize) { (self.x.len(), self.s.len()) } From 665177c1e2248da7e55c627c14de92890ba51c73 Mon Sep 17 00:00:00 2001 From: goulart-paul Date: Sat, 1 Jun 2024 17:08:35 +0100 Subject: [PATCH 13/13] update for docs --- src/python/impl_default_py.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/src/python/impl_default_py.rs b/src/python/impl_default_py.rs index db09381..691326a 100644 --- a/src/python/impl_default_py.rs +++ b/src/python/impl_default_py.rs @@ -10,6 +10,7 @@ use crate::solver::{ IPSolver, SolverStatus, }, implementations::default::*, + SolverJSONReadWrite, }; use num_derive::ToPrimitive; use num_traits::ToPrimitive;