diff --git a/.Rbuildignore b/.Rbuildignore index 21d9d939..936d869d 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -10,6 +10,10 @@ ^src/rust/vendor.sh ^src/rust/fix_hash.R ^src/rust/Cargo.lock +^src/rust/vendor$ ^src/*.o ^src/*.so ^update_authors.R$ + +^src/rust/.cargo$ +^src/Makevars$ diff --git a/.gitignore b/.gitignore index 1458c96a..24de750f 100644 --- a/.gitignore +++ b/.gitignore @@ -43,3 +43,6 @@ vignettes/*.pdf # My stuff inst/stuff + +# Makevars is generated +src/Makevars \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 74322caf..64cab7ee 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: clarabel Type: Package Title: Interior Point Conic Optimization Solver -Version: 0.5.1 +Version: 0.9.0 Authors@R: c(person("Balasubramanian", "Narasimhan", role=c("aut", "cre"), email = "naras@stanford.edu"), person("Paul", "Goulart", role=c("aut", "cph")), @@ -10,13 +10,12 @@ Authors@R: c(person("Balasubramanian", "Narasimhan", role=c("aut", "cre"), comment = "For vendoring/Makefile hints/R scripts for generating crate authors/licenses"), person(given = "The authors of the dependency Rust crates", role = c("ctb"), comment = "see inst/AUTHORS file for details")) -Description: A versatile interior point solver that solves linear programs (LPs), quadratic programs (QPs), second-order cone programs (SOCPs), semidefinite programs (SDPs), and problems with exponential and power cone constraints (). For quadratic objectives, unlike interior point solvers based on the standard homogeneous self-dual embedding (HSDE) model, 'Clarabel' handles quadratic objective without requiring any epigraphical reformulation of its objective function. It can therefore be significantly faster than other HSDE-based solvers for problems with quadratic objective functions. Infeasible problems are detected using using a homogeneous embedding technique. +Description: A versatile interior point solver that solves linear programs (LPs), quadratic programs (QPs), second-order cone programs (SOCPs), semidefinite programs (SDPs), and problems with exponential and power cone constraints (). For quadratic objectives, unlike interior point solvers based on the standard homogeneous self-dual embedding (HSDE) model, 'Clarabel' handles quadratic objective without requiring any epigraphical reformulation of its objective function. It can therefore be significantly faster than other HSDE-based solvers for problems with quadratic objective functions. Infeasible problems are detected using using a homogeneous embedding technique. License: Apache License (== 2.0) -SystemRequirements: Cargo (rustc package manager) and GNU make Encoding: UTF-8 Roxygen: list(markdown = TRUE) Config/rextendr/version: 0.3.1 -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.1 URL: https://oxfordcontrol.github.io/clarabel-r/ BugReports: https://github.com/oxfordcontrol/clarabel-r/issues Suggests: @@ -25,3 +24,6 @@ Suggests: rmarkdown, tinytest VignetteBuilder: knitr +SystemRequirements: Cargo (Rust package manager), rustc and GNU Make +Imports: + methods diff --git a/LICENSE.note b/LICENSE.note index a8bb0f92..fd7e598c 100644 --- a/LICENSE.note +++ b/LICENSE.note @@ -1,20 +1,16 @@ This package contains the Rust source code of the dependencies in src/rust/vendor.tar.xz The authorships and the licenses are listed below. In summary, all libraries are distributed either under the MIT license or under MIT/Apache-2.0 dual license [1]. - Note that, when Cargo (Rust’s build system and package manager) is not installed on the machine, the pre-compiled binary will be downloaded on building this package. The binary is compiled using the same Rust code, so the authorships and the licenses are the same as listed here. - [1]: The unicode-indent library shows 'Unicode-DFS-2016' license because it contains some test data generated by using the Unicode Character Database. So, this license is not applied to the actual sources that get compiled. Please refer to the License section of the library's README (https://crates.io/crates/unicode-ident) for the details. - =============================== - Name: aho-corasick Files: vendor/aho-corasick/* Authors: Andrew Gallant @@ -38,7 +34,7 @@ License: Apache-2.0 OR MIT Name: blas-src Files: vendor/blas-src/* -Authors: Balasubramanian Narasimhan, Ivan Ukhov, Jed Brown, Stefan Kroboth, Toshiki Teramura, bluss +Authors: Balasubramanian Narasimhan, Ivan Ukhov, Jed Brown, Michael Zietz, Stefan Kroboth, Toshiki Teramura, bluss License: Apache-2.0/MIT ------------------------------ @@ -57,6 +53,13 @@ License: Apache-2.0/MIT ------------------------------ +Name: cc +Files: vendor/cc/* +Authors: Alex Crichton +License: MIT OR Apache-2.0 + +------------------------------ + Name: cfg-if Files: vendor/cfg-if/* Authors: Alex Crichton @@ -113,31 +116,24 @@ License: MIT/Apache-2.0 ------------------------------ -Name: enum_dispatch -Files: vendor/enum_dispatch/* -Authors: Anton Lazarev +Name: either +Files: vendor/either/* +Authors: bluss License: MIT OR Apache-2.0 ------------------------------ -Name: extendr-api -Files: vendor/extendr-api/* -Authors: andy-thomason, Thomas Down, Mossa Merhi Reimert, Claus O. Wilke, Hiroaki Yutani, Ilia A. Kosenkov, Michael Milton -License: MIT - ------------------------------- - -Name: extendr-engine -Files: vendor/extendr-engine/* -Authors: andy-thomason, Thomas Down, Mossa Merhi Reimert, Claus O. Wilke, Hiroaki Yutani, Ilia A. Kosenkov -License: MIT +Name: enum_dispatch +Files: vendor/enum_dispatch/* +Authors: Anton Lazarev +License: MIT OR Apache-2.0 ------------------------------ -Name: extendr-macros -Files: vendor/extendr-macros/* -Authors: andy-thomason, Thomas Down, Mossa Merhi Reimert, Claus O. Wilke, Hiroaki Yutani, Ilia A. Kosenkov -License: MIT +Name: equivalent +Files: vendor/equivalent/* +Authors: +License: Apache-2.0 OR MIT ------------------------------ @@ -148,6 +144,13 @@ License: Apache-2.0 / MIT ------------------------------ +Name: hashbrown +Files: vendor/hashbrown/* +Authors: Amanieu d'Antras +License: MIT OR Apache-2.0 + +------------------------------ + Name: ident_case Files: vendor/ident_case/* Authors: Ted Driggs @@ -155,9 +158,30 @@ License: MIT/Apache-2.0 ------------------------------ +Name: indexmap +Files: vendor/indexmap/* +Authors: +License: Apache-2.0 OR MIT + +------------------------------ + +Name: itertools +Files: vendor/itertools/* +Authors: bluss +License: MIT OR Apache-2.0 + +------------------------------ + +Name: itoa +Files: vendor/itoa/* +Authors: David Tolnay +License: MIT OR Apache-2.0 + +------------------------------ + Name: lapack-src Files: vendor/lapack-src/* -Authors: Balasubramanian Narasimhan, Ivan Ukhov, Mitsutoshi Aoe, Stefan Kroboth, Toshiki Teramura +Authors: Balasubramanian Narasimhan, Ivan Ukhov, Michael Zietz, Mitsutoshi Aoe, Stefan Kroboth, Toshiki Teramura License: Apache-2.0/MIT ------------------------------ @@ -190,17 +214,10 @@ License: MIT OR Apache-2.0 ------------------------------ -Name: libR-sys -Files: vendor/libR-sys/* -Authors: andy-thomason, Thomas Down, Mossa Merhi Reimert, Claus O. Wilke, Ilia A. Kosenkov, Hiroaki Yutani -License: MIT - ------------------------------- - Name: memchr Files: vendor/memchr/* Authors: Andrew Gallant, bluss -License: Unlicense/MIT +License: Unlicense OR MIT ------------------------------ @@ -225,13 +242,6 @@ License: MIT OR Apache-2.0 ------------------------------ -Name: paste -Files: vendor/paste/* -Authors: David Tolnay -License: MIT OR Apache-2.0 - ------------------------------- - Name: proc-macro2 Files: vendor/proc-macro2/* Authors: David Tolnay, Alex Crichton @@ -253,16 +263,79 @@ License: MIT ------------------------------ +Name: regex-automata +Files: vendor/regex-automata/* +Authors: The Rust Project Developers, Andrew Gallant +License: MIT OR Apache-2.0 + +------------------------------ + Name: regex-syntax Files: vendor/regex-syntax/* -Authors: The Rust Project Developers +Authors: The Rust Project Developers, Andrew Gallant License: MIT OR Apache-2.0 ------------------------------ Name: regex Files: vendor/regex/* -Authors: The Rust Project Developers +Authors: The Rust Project Developers, Andrew Gallant +License: MIT OR Apache-2.0 + +------------------------------ + +Name: ryu +Files: vendor/ryu/* +Authors: David Tolnay +License: Apache-2.0 OR BSL-1.0 + +------------------------------ + +Name: savvy-bindgen +Files: vendor/savvy-bindgen/* +Authors: Hiroaki Yutani +License: MIT + +------------------------------ + +Name: savvy-ffi +Files: vendor/savvy-ffi/* +Authors: Hiroaki Yutani +License: MIT + +------------------------------ + +Name: savvy-macro +Files: vendor/savvy-macro/* +Authors: Hiroaki Yutani +License: MIT + +------------------------------ + +Name: savvy +Files: vendor/savvy/* +Authors: Hiroaki Yutani +License: MIT + +------------------------------ + +Name: serde_derive +Files: vendor/serde_derive/* +Authors: Erick Tryzelaar, David Tolnay +License: MIT OR Apache-2.0 + +------------------------------ + +Name: serde_json +Files: vendor/serde_json/* +Authors: Erick Tryzelaar, David Tolnay +License: MIT OR Apache-2.0 + +------------------------------ + +Name: serde +Files: vendor/serde/* +Authors: Erick Tryzelaar, David Tolnay License: MIT OR Apache-2.0 ------------------------------ diff --git a/NAMESPACE b/NAMESPACE index a9020475..15269377 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,5 +3,6 @@ export(clarabel) export(clarabel_control) export(solver_status_descriptions) +importFrom(methods,as) useDynLib(clarabel) useDynLib(clarabel, .registration = TRUE) diff --git a/NEWS.md b/NEWS.md index d7fa8f7f..121a6e97 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,11 @@ +# clarabel 0.9.0 + +- Synced up to version 0.9.0 of Clarabel.rs +- Added all applicable tests from Clarabel.rs +- Updated documentation of cone specification +- Updated vignette +- Switched to savvy from rextendr + # clarabel 0.5.1 - Clarabel now supports semidefinite programs (syncing up to version 0.5.1 of Clarabel.rs) diff --git a/R/000-wrappers.R b/R/000-wrappers.R new file mode 100644 index 00000000..d8ecd860 --- /dev/null +++ b/R/000-wrappers.R @@ -0,0 +1,31 @@ +# Generated by savvy: do not edit by hand +# +# Note: +# This wrapper file is named as `000-wrappers.R` so that this file is loaded +# first, which allows users to override the functions defined here (e.g., a +# print() method for an enum). + +#' @useDynLib clarabel, .registration = TRUE +#' @keywords internal +NULL + +# Check class and extract the external pointer embedded in the environment +.savvy_extract_ptr <- function(e, class) { + if(is.null(e)) { + return(NULL) + } + + if(inherits(e, class)) { + e$.ptr + } else { + msg <- paste0("Expected ", class, ", got ", class(e)[1]) + stop(msg, call. = FALSE) + } +} + + +clarabel_solve <- function(m, n, Ai, Ap, Ax, b, q, Pi, Pp, Px, cone_spec, r_settings) { + .Call(savvy_clarabel_solve__impl, m, n, Ai, Ap, Ax, b, q, Pi, Pp, Px, cone_spec, r_settings) +} + + diff --git a/R/clarabel-package.R b/R/clarabel-package.R index 5da8ed44..eecd81bb 100644 --- a/R/clarabel-package.R +++ b/R/clarabel-package.R @@ -1,10 +1,9 @@ #' Interface to Clarabel solver implemented in Rust. #' -#' @description Clarabel is a versatile interior point solver for convex programs using a new homogeneous embedding. It solves solves linear programs (LPs), quadratic programs (QPs), second-order cone programs (SOCPs), and problems with exponential and power cone constraints. For quadratic objectives, unlike interior point solvers based on the standard homogeneous self-dual embedding (HSDE) model, Clarabel handles quadratic objective without requiring any epigraphical reformulation of its objective function. It can therefore be significantly faster than other HSDE-based solvers for problems with quadratic objective functions. Infeasible problems are detected using a homogeneous embedding technique. See . +#' @description Clarabel is a versatile interior point solver for convex programs using a new homogeneous embedding. It solves solves linear programs (LPs), quadratic programs (QPs), second-order cone programs (SOCPs), and problems with exponential and power cone constraints. For quadratic objectives, unlike interior point solvers based on the standard homogeneous self-dual embedding (HSDE) model, Clarabel handles quadratic objective without requiring any epigraphical reformulation of its objective function. It can therefore be significantly faster than other HSDE-based solvers for problems with quadratic objective functions. Infeasible problems are detected using a homogeneous embedding technique. See . #' #' @name clarabel-package #' @useDynLib clarabel -#' @docType package #' @author Balasubramanian Narasimhan, Paul Goulart, Yuwen Chen -#' @keywords package -NULL +#' @keywords internal +"_PACKAGE" diff --git a/R/clarabel.R b/R/clarabel.R index 8fcdbc6c..bf90a822 100644 --- a/R/clarabel.R +++ b/R/clarabel.R @@ -56,16 +56,17 @@ #' listed in \dQuote{Cone Parameters} below } } #' #' \subsection{Cone Parameters}{ -#' The table below shows the cone parameter specifications +#' The table below shows the cone parameter specifications. Mathematical definitions are in the vignette. #' \tabular{rllll}{ #' \tab \bold{Parameter} \tab \bold{Type} \tab \bold{Length} \tab \bold{Description} \cr -#' \tab \code{z} \tab integer \tab \eqn{1} \tab number of primal zero cones (dual free cones), \cr -#' \tab \tab \tab \tab which corresponds to the primal equality constraints \cr -#' \tab \code{l} \tab integer \tab \eqn{1} \tab number of linear cones (non-negative cones) \cr -#' \tab \code{q} \tab integer \tab \eqn{\geq1} \tab vector of second-order cone sizes \cr -#' \tab \code{s} \tab integer \tab \eqn{\geq1} \tab vector of positive semidefinite cone sizes \cr -#' \tab \code{ep} \tab integer \tab \eqn{1} \tab number of primal exponential cones \cr -#' \tab \code{p} \tab numeric \tab \eqn{\geq1} \tab vector of primal power cone parameters +#' \tab \code{z} \tab integer \tab \eqn{1} \tab number of primal zero cones (dual free cones), \cr +#' \tab \tab \tab \tab which corresponds to the primal equality constraints \cr +#' \tab \code{l} \tab integer \tab \eqn{1} \tab number of linear cones (non-negative cones) \cr +#' \tab \code{q} \tab integer \tab \eqn{\ge 1} \tab vector of second-order cone sizes \cr +#' \tab \code{s} \tab integer \tab \eqn{\ge 1} \tab vector of positive semidefinite cone sizes \cr +#' \tab \code{ep} \tab integer \tab \eqn{1} \tab number of primal exponential cones \cr +#' \tab \code{p} \tab numeric \tab \eqn{\ge 1} \tab vector of primal power cone parameters \cr +#' \tab \code{gp} \tab list \tab \eqn{\ge 1} \tab list of named lists of two items, `a` : a numeric vector of at least 2 exponent terms in the product summing to 1, and `n` : an integer dimension of generalized power cone parameters #' } } #' #' When the parameter `strict_cone_order` is `FALSE`, one can specify @@ -77,7 +78,9 @@ #' 2L, l1 = 2L, q1 = 2L, zb = 3L, qx = 3L)`, indicating a zero #' cone of size 2, followed by a linear cone of size 2, followed by a second-order #' cone of size 2, followed by a zero cone of size 3, and finally a second-order -#' cone of size 3. +#' cone of size 3. Generalized power cones parameters have to specified as named lists, e.g., `list(z = 2L, gp1 = list(a = c(0.3, 0.7), n = 3L), gp2 = list(a = c(0.5, 0.5), n = 1L))`. +#' +#' _Note that when `strict_cone_order = FALSE`, types of cone parameters such as integers, reals have to be explicit since the parameters are directly passed to the Rust interface with no sanity checks.!_ #' #' @examples #' A <- matrix(c(1, 1), ncol = 1) @@ -91,14 +94,29 @@ clarabel <- function(A, b, q, P = NULL, cones, control = list(), strict_cone_order = TRUE) { - ## TBD check box cone parameters, bsize > 0 & bl, bu have lengths bsize - 1 - - n_variables <- NCOL(A) - n_constraints <- NROW(A) + m <- length(b); n <- length(q); + n_variables <- ncol(A) + n_constraints <- nrow(A) + + if (m != n_constraints) stop("A and b incompatible dimensions.") + if (n != n_variables) stop("A and q incompatible dimensions.") + if (!is.null(P)) { + if (n != ncol(P)) stop("P and q incompatible dimensions.") + if (ncol(P) != nrow(P)) stop("P not square.") + } - ## if (n_variables <= 0) { - ## stop("No variables") - ## } + # Sanitize control parameters + control <- do.call(clarabel_control, control) + if (strict_cone_order) { + cones_and_nvars <- sanitize_cone_spec(cones) + cones <- cones_and_nvars[["cones"]] + nvars <- cones_and_nvars[["nvars"]] + } else { + nvars <- nvars(cones) + } + if (sum(nvars) != m) stop("Constraint dimensions inconsistent with size of cones.") + + ## TBD check box cone parameters, bsize > 0 & bl, bu have lengths bsize - 1 if ( inherits(A, "dgCMatrix") ) { Ai <- A@i @@ -111,9 +129,6 @@ clarabel <- function(A, b, q, P = NULL, cones, control = list(), Ax <- csc[["values"]] } - ## data <- list(m = n_constraints, n = n_variables, c = obj, - ## Ai = Ai, Ap = Ap, Ax = Ax, b = b) - if (!is.null(P)) { if (inherits(P, "dsCMatrix") ) { Pi <- P@i @@ -131,10 +146,9 @@ clarabel <- function(A, b, q, P = NULL, cones, control = list(), Px <- numeric(0) } - # Sanitize control parameters - control <- do.call(clarabel_control, control) - if (strict_cone_order) cones <- sanitize_cone_spec(cones) - clarabel_solve(n_constraints, n_variables, Ai, Ap, Ax, b, q, Pi, Pp, Px, cones, control) + .Call(savvy_clarabel_solve__impl, n_constraints, n_variables, Ai, Ap, Ax, b, q, Pi, Pp, Px, cones, control, PACKAGE = "clarabel") + + ## clarabel_solve(n_constraints, n_variables, Ai, Ap, Ax, b, q, Pi, Pp, Px, cones, control) } #' Control parameters with default values and types in parenthesis @@ -176,6 +190,10 @@ clarabel <- function(A, b, q, P = NULL, cones, control = list(), #' @param iterative_refinement_max_iter iterative refinement maximum iterations (`10L`) #' @param iterative_refinement_stop_ratio iterative refinement stalling tolerance (`5.0`) #' @param presolve_enable whether to enable presolvle (`TRUE`) +#' @param chordal_decomposition_enable whether to enable chordal decomposition for SDPs (`FALSE`) +#' @param chordal_decomposition_merge_method chordal decomposition merge method, one of `'none'`, `'parent_child'` or `'clique_graph'`, for SDPs (`'none'`) +#' @param chordal_decomposition_compact a boolean flag for SDPs indicating whether to assemble decomposed system in _compact_ form for SDPs (`FALSE`) +#' @param chordal_decomposition_complete_dual a boolean flag indicating complete PSD dual variables after decomposition for SDPs #' @return a list containing the control parameters. #' @export clarabel_control clarabel_control <- function( @@ -223,20 +241,30 @@ clarabel_control <- function( iterative_refinement_abstol = 1e-12, iterative_refinement_max_iter = 10L, iterative_refinement_stop_ratio = 5.0, - presolve_enable = TRUE) { + presolve_enable = TRUE, + chordal_decomposition_enable = FALSE, + chordal_decomposition_merge_method = c('none', 'parent_child', 'clique_graph'), + chordal_decomposition_compact = FALSE, + chordal_decomposition_complete_dual = FALSE + ) { params <- as.list(environment()) - ## Match string arg to make sure it is kosher + + ## Match string args to make sure it is kosher direct_solve_method <- match.arg(direct_solve_method) params$direct_solve_method <- direct_solve_method + chordal_decomposition_merge_method <- match.arg(chordal_decomposition_merge_method) + params$chordal_decomposition_merge_method <- chordal_decomposition_merge_method ## Rust has strict type and length checks, so try to avoid panics bool_params <- c("verbose", "equilibrate_enable", "direct_kkt_solver", "static_regularization_enable", - "dynamic_regularization_enable", "iterative_refinement_enable", "presolve_enable") + "dynamic_regularization_enable", "iterative_refinement_enable", "presolve_enable", + "chordal_decomposition_enable", "chordal_decomposition_compact", + "chordal_decomposition_complete_dual") int_params <- c("max_iter", "equilibrate_max_iter", "iterative_refinement_max_iter") - string_params <- "direct_solve_method" # Might need to uncomment character coercion below, if length > 1 + string_params <- c("direct_solve_method", "chordal_decomposition_merge_method") # Might need to uncomment character coercion below, if length > 1 if (any(sapply(params, length) != 1L)) stop("clarabel_control: arguments should be scalars!") if (any(unlist(params[int_params]) < 0)) stop("clarabel_control: integer arguments should be >= 0!") @@ -283,52 +311,128 @@ solver_status_descriptions <- function() { ### Sanitize cone specifications ### @param cone_spec a list of cone specifications -### @return sanitized cone specifications +### @return a named list of sanitized cone specifications and the number of variables for each cone (`nvars`) sanitize_cone_spec <- function(cone_spec) { cone_names <- names(cone_spec) - + ## Simple sanity checks if ((nc <- length(cone_names)) == 0L) { stop("sanitize_cone_spec: no cone parameters specified") } if (length(intersect(cone_names, c("z", "l", "q", "s", "ep", "p"))) != nc) { - stop("sanitize_cone_spec: unknown cone parameters specified") + stop("sanitize_cone_spec: repeated cone parameters or unknown cone parameters specified") } - + ## Check lengths as noted cone parameters table for ?clarabel ## First, scalars - z <- as.integer(cone_spec[["z"]]); zl <- length(z) - l <- as.integer(cone_spec[["l"]]); ll <- length(l) - ep <- as.integer(cone_spec[["ep"]]); epl <- length(ep) + z <- as.integer(cone_spec[["z"]]); zl <- length(z); nvar_z <- sum(z); + l <- as.integer(cone_spec[["l"]]); ll <- length(l); nvar_l <- sum(l); + ep <- as.integer(cone_spec[["ep"]]); epl <- length(ep); nvar_ep <- sum(ep) * 3 ## 3 variables per exp cone if (zl > 1 || ll > 1 || epl > 1) { stop("sanitize_cone_spec: z, l, ep should be scalars") } if (any(c(z, l, ep) < 0L)) { stop("sanitize_cone_spec: z, l, ep should be scalars > 0") } - + ## Now the others ## SOC - q <- as.integer(cone_spec[["q"]]); ql <- length(q) + q <- as.integer(cone_spec[["q"]]); ql <- length(q); nvar_q <- sum(q); if (any(q <= 0L)) stop("sanitize_cone_spec: SOC dimensions should be > 0") q <- as.list(q); names(q) <- rep("q", ql); - + ## PSD - s <- as.integer(cone_spec[["s"]]); sl <- length(s) + s <- as.integer(cone_spec[["s"]]); sl <- length(s); nvar_s <- if (sl > 0) sum(sapply(s, triangular_number)) else 0; ## triangular number if (any(s <= 0L)) stop("sanitize_cone_spec: PSD dimensions should be > 0") s <- as.list(s); names(s) <- rep("s", sl); - + ## Power Cone - p <- as.numeric(cone_spec[["p"]]); pl <- length(p) + p <- as.numeric(cone_spec[["p"]]); pl <- length(p); nvar_p <- pl * 3; ## 3 variables per power cone if (any(p <= 0)) stop("sanitize_cone_spec: Power cone parameter should be > 0") p <- as.list(p); names(p) <- rep("p", pl); + + ## Power Cone + gp <- cone_spec[["gp"]]; gpl <- length(gp); result <- sanitize_gp_params_and_get_nvars(gp); + gp <- result$sanitized_gp_params; names(gp) <- rep("gp", gpl); nvar_gp <- result$nvar; - c( + cones <- c( if (zl > 0) list(z = z), if (ll > 0) list(l = l), if (ql > 0) q, if (sl > 0) s, if (epl > 0) list(ep = ep), - if (pl > 0) p + if (pl > 0) p, + if (gpl > 0) gp ) + + nvars <- c(z = nvar_z, + l = nvar_l, + q = nvar_q, + s = nvar_s, + ep = nvar_ep, + p = nvar_p, + gp = nvar_gp) + list(cones = cones, nvars = nvars) } + +### Return the number of variables used for each type of cone based on cone specification +### @param cones the cone specifications, expecting strict_cone_order to be FALSE +### @return a named vector of number of variables per each type of cone +nvars <- function(cones) { + cone_names <- names(cones) + ## separate out gp cones from others. + gp_indices <- grep("^gp", cone_names) + if (length(gp_indices) > 0) { + gp_cones <- cones[gp_indices] + nvar_gp <- sum(sapply(gp_cones, function(x) length(x[["a"]]) + x[["n"]])) + } else { + nvar_gp <- 0L + } + other_indices <- grep("^gp", cone_names, invert = TRUE) + cones <- unlist(cones[other_indices]) + cone_names <- names(cones) + nvar_z <- if (length(matched <- grep("^z", cone_names)) > 0) sum(cones[matched]) else 0L + nvar_l <- if (length(matched <- grep("^l", cone_names)) > 0) sum(cones[matched]) else 0L + nvar_q <- if (length(matched <- grep("^q", cone_names)) > 0) sum(cones[matched]) else 0L + nvar_s <- if (length(matched <- grep("^s", cone_names)) > 0) sum(sapply(cones[matched], triangular_number)) else 0L + nvar_ep <- if (length(matched <- grep("^ep", cone_names)) > 0) 3 * sum(cones[matched]) else 0L + nvar_p <- if (length(matched <- grep("^p", cone_names)) > 0) 3 * length(cones[matched]) else 0L + c(z = nvar_z, + l = nvar_l, + q = nvar_q, + s = nvar_s, + ep = nvar_ep, + p = nvar_p, + gp = nvar_gp) +} + +## Return the n-th triangular number +triangular_number <- function(n) { + n * (n + 1) / 2 +} + +## Check Generalized Power Cone Params +sanitize_gp_params_and_get_nvars <- function(gp) { + # Generalized power cone α.len() + dim + if (length(gp) == 0) { + list(sanitized_gp_params = list(), nvar = 0L) + } else { + sanitized_gp_params <- + lapply(gp, function(x) { + par_names <- sort(names(x)) + if (length(x) != 2L || !identical(par_names, c("a", "n"))) { + stop("Generalized power cone param list should be a list of two elements named 'a' and 'n'") + } + exps <- x[["a"]] + if (length(exps) < 2L || any(exps <= 0) || any(exps >= 1) || abs(sum(exps) - 1.0) > 0.0) { + stop("Improper Generalized power cone exponents!") + } + n <- x[["n"]] + if (length(n) != 1L || n <= 0) stop("Improper Generalized power cone dimension!") + list(a = as.numeric(exps), n = as.integer(n)) + }) + list(sanitized_gp_params = sanitized_gp_params, nvar = sum(sapply(sanitized_gp_params, function(x) length(x[["a"]]) + x[["n"]]))) + } +} + + diff --git a/R/extendr-wrappers.R b/R/extendr-wrappers.R deleted file mode 100644 index b6610fd9..00000000 --- a/R/extendr-wrappers.R +++ /dev/null @@ -1,17 +0,0 @@ -# Generated by extendr: Do not edit by hand - -# nolint start - -# -# This file was created with the following call: -# .Call("wrap__make_clarabel_wrappers", use_symbols = TRUE, package_name = "clarabel") - -#' @docType package -#' @usage NULL -#' @useDynLib clarabel, .registration = TRUE -NULL - -clarabel_solve <- function(m, n, Ai, Ap, Ax, b, q, Pi, Pp, Px, cone_spec, r_settings) .Call(wrap__clarabel_solve, m, n, Ai, Ap, Ax, b, q, Pi, Pp, Px, cone_spec, r_settings) - - -# nolint end diff --git a/R/sparse.R b/R/sparse.R index 66586612..5c1fd9a4 100644 --- a/R/sparse.R +++ b/R/sparse.R @@ -13,6 +13,7 @@ make_csc_matrix <- function(x) UseMethod("make_csc_matrix") +#' @method make_csc_matrix matrix make_csc_matrix.matrix <- function(x) { if(!is.matrix(x)) stop("Argument 'x' must be a matrix.") @@ -23,6 +24,7 @@ make_csc_matrix.matrix <- function(x) { values = x[ind]) } +#' @method make_csc_matrix simple_triplet_matrix make_csc_matrix.simple_triplet_matrix <- function(x) { if(!inherits(x, "simple_triplet_matrix")) stop("Argument 'x' must be of class 'simple_triplet_matrix'.") @@ -40,6 +42,7 @@ make_csc_matrix.simple_triplet_matrix <- function(x) { make_csc_symm_matrix <- function(m) UseMethod("make_csc_symm_matrix") +#' @method make_csc_symm_matrix matrix make_csc_symm_matrix.matrix <- function(m) { ind <- which(m!=0 ,arr.ind = TRUE) ind <- ind[ind[, 1] <= ind[, 2], ] ## keep upper part only @@ -51,6 +54,7 @@ make_csc_symm_matrix.matrix <- function(m) { values = values) } +#' @method make_csc_symm_matrix simple_triplet_matrix make_csc_symm_matrix.simple_triplet_matrix <- function(m) { ind <- which(m$i <= m$j) x <- list(i = m$i[ind] + 1L, j = m$j[ind] + 1L) ##make it 1-based @@ -61,3 +65,12 @@ make_csc_symm_matrix.simple_triplet_matrix <- function(m) { values = values) } +#' @method make_csc_symm_matrix dgCMatrix +#' @importFrom methods as +make_csc_symm_matrix.dgCMatrix <- function(m) { + x <- as(m, "symmetricMatrix") + list(matbeg = x@p, + matind = x@i, + values = x@x) +} + diff --git a/README.md b/README.md index 55a20fad..ff57338e 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ [![](https://cranlogs.r-pkg.org/badges/clarabel)](https://CRAN.R-project.org/package=clarabel) R interface to the -[Clarabel](https://oxfordcontrol.github.io/ClarabelDocs/stable/) +[Clarabel](https://clarabel.org/stable/) interior point conic optimization solver from the [Oxford Control Group](https://github.com/oxfordcontrol). diff --git a/configure b/configure new file mode 100755 index 00000000..2add568b --- /dev/null +++ b/configure @@ -0,0 +1,5 @@ +if [ "$(uname)" = "Emscripten" ]; then + TARGET="wasm32-unknown-emscripten" +fi + +sed -e "s/@TARGET@/${TARGET}/" src/Makevars.in > src/Makevars diff --git a/docs/404.html b/docs/404.html index e3f50325..a11e1f49 100644 --- a/docs/404.html +++ b/docs/404.html @@ -12,7 +12,7 @@ - + @@ -39,7 +39,7 @@ clarabel - 0.5.1 + 0.9.0 @@ -98,7 +98,7 @@

Page not found (404)

-

Site built with pkgdown 2.0.7.

+

Site built with pkgdown 2.0.9.

diff --git a/docs/articles/clarabel.html b/docs/articles/clarabel.html index b258baa0..379f1a4a 100644 --- a/docs/articles/clarabel.html +++ b/docs/articles/clarabel.html @@ -12,7 +12,7 @@ - + @@ -40,7 +40,7 @@ clarabel - 0.5.1 + 0.9.0 @@ -89,8 +89,8 @@

Clarabel Solver Examples

Introduction

-

The first two examples are from the original Clarabel -documentation and the third is from SCS.

+

The first two examples are from the original Clarabel documentation and the third is +from SCS.

1. Basic Quadratic Program Example @@ -154,8 +154,9 @@

1.3. Solutions <- clarabel(A = A, b = b, q = q, P = P, cones = cones) cat(sprintf("Solution status, description: = (%d, %s)\n", s$status, solver_status_descriptions()[s$status])) -#> Solution status, description: = (2, Solver terminated with a solution.) -cat(sprintf("Solution: (x1, x2) = (%f, %f)\n", s$x[1], s$x[2])) +#> Solution status, description: = (2, Solver terminated with a solution.)

+
+cat(sprintf("Solution: (x1, x2) = (%f, %f)\n", s$x[1], s$x[2]))
 #> Solution: (x1, x2) = (0.428571, 0.214286)
@@ -202,7 +203,7 @@

2.2. Constraints

2.3. Solution

-
+
 P <- Matrix::Matrix(2 * c(0, 0, 0, 1), nrow = 2, ncol = 2, sparse = TRUE)
 P <- as(P, "symmetricMatrix") # P needs to be a symmetric matrix
 q <- c(0, 0)
@@ -212,8 +213,9 @@ 

2.3. Solutions <- clarabel(A = A, b = b, q = q, P = P, cones = cones) cat(sprintf("Solution status, description: = (%d, %s)\n", s$status, solver_status_descriptions()[s$status])) -#> Solution status, description: = (2, Solver terminated with a solution.) -cat(sprintf("Solution (x1, x2) = (%f, %f)\n", s$x[1], s$x[2])) +#> Solution status, description: = (2, Solver terminated with a solution.)

+
+cat(sprintf("Solution (x1, x2) = (%f, %f)\n", s$x[1], s$x[2]))
 #> Solution (x1, x2) = (1.000000, -1.000000)
@@ -269,8 +271,8 @@

3. Semidefinite Cone Programming\(1/\sqrt{2}\) and the upper triangular entries are filled in by copying the values of lower triangular entries. -Explicitly, the inverse operation takes vector \(s \in \mathbb{R}^{k(k+1)/2}\) and produces -the matrix

+Explicitly, the inverse operation takes vector \(s \in +\mathbb{R}^{k(k+1)/2}\) and produces the matrix

\[ \begin{aligned} \text{mat}(s) = \begin{bmatrix} @@ -291,7 +293,7 @@

3. Semidefinite Cone Programming

Below are two functions to implement both \(\text{vec}\) and \(\text{mat}\).

-
+
 #' Return an vectorization of symmetric matrix using the upper triangular part,
 #' still in column order.
 #' @param S a symmetric matrix
@@ -353,13 +355,16 @@ 

3.1. Example\(x \in \mathbb{R}^n\) and -\(S \in \mathbb{R}^{k \times k}\)

+\(S \in +\mathbb{R}^{k \times k}\)

\[ B - \sum_{i=1}^n \mathcal{A}_i x_i = S \succeq 0 \]

where data \(B, \mathcal{A}_1, \ldots, -\mathcal{A}_n \in \mathbb{R}^{k \times k}\) are symmetric. We can -write this in the canonical form over a new variable \(s \in \mathcal{S}_+^k\):

+\mathcal{A}_n \in \mathbb{R}^{k +\times k}\) are symmetric. We can write this in the canonical +form over a new variable \(s \in +\mathcal{S}_+^k\):

\[ \begin{aligned} \begin{align} @@ -371,8 +376,8 @@

3.1. Example\(\text{vec}\) -is linear, where \(b = \text{vec}(B)\) -and

+is linear, where \(b = +\text{vec}(B)\) and

\[ A = \begin{bmatrix} @@ -388,7 +393,7 @@

3.1. Example
+
 q <- c(1, -1, 1) # objective: x_1 - x2 + x_3
 A11 <- matrix(c(-7, -11, -11, 3), nrow = 2)
 A12 <- matrix(c(7, -18, -18, 8), nrow = 2)
@@ -407,22 +412,122 @@ 

3.1. Example) b <- c(vec(B1), vec(B2)) # stack both psd constraints cones <- list(s = c(2, 3)) # cone dimensions -s <- clarabel(A = A, b = b, q = q, cone = cones) +s <- clarabel(A = A, b = b, q = q, cones = cones) cat(sprintf("Solution status, description: = (%d, %s)\n", s$status, solver_status_descriptions()[s$status])) -#> Solution status, description: = (2, Solver terminated with a solution.) -cat(sprintf("Solution (x1, x2, x3) = (%f, %f, %f)\n", s$x[1], s$x[2], s$x[3])) -#> Solution (x1, x2, x3) = (-0.367752, 1.898333, -0.887460)

+#> Solution status, description: = (2, Solver terminated with a solution.)

+
+cat(sprintf("Solution (x1, x2, x3) = (%f, %f, %f)\n", s$x[1], s$x[2], s$x[3]))
+#> Solution (x1, x2, x3) = (-0.367746, 1.898333, -0.887466)
-

4. Control parameters +

4. Cone Specifications +

+

The following cones can be specified in Clarabel.

+ +++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
ParameterTypeLengthDescriptionDefinition (per parameter element)
zinteger1Number of primal zero cones (dual free cones), which +corresponds to the primal equality constraints\(\{ 0 \}^{z}\)
linteger1Number of linear cones (non-negative cones)\(\{ x \in \mathbb{R}^{l} : +x_i \ge 0, \forall i=1,\dots,l \}\)
qinteger>= 1Vector of second-order cone sizes\(\{ (t,x) \in +\mathbb{R}^{q} : \lVert x\rVert_2 \leq t \}\)
sinteger>= 1Vector of positive semidefinite cone sizesUpper triangular part of the positive semidefinite cone +\(S^s_+\). The elements \(x\) of this cone represent the columnwise +stacking of the upper triangular part of a positive semidefinite matrix +\(X \in S^s_+\), so that \(x \in R^d\) with \(d = s(s+1)/2\) +
epinteger1Number of primal exponential cones\(\{(x, y, z) : y > 0,~~ +ye^{x/y} \le z \}\)
pnumeric>= 1Vector of primal power cone parameters +\(\{(x, y, z) : x^p y^{(1-p)} +\ge \lVert z\rVert,~ (x,y) \ge 0 \}\) with \(p \in (0,1)\) +
gplist>= 1List of named lists of two items, a : the +numeric vector of at least 2 exponent terms, and n : an +integer dimension of generalized power cone parameters +\(\{(x, y) \in R^{len(a)} +\times R^n : \prod\limits_{a_i \in a} x_i^{a_i} \ge \lVert y\rVert_2,~ x +\ge 0 \}\) with \(a_i \in +(0,1)\) and \(\sum a_i = +1\) +
+

Generalized power cone parameters are specified as list of two-item +lists, with component named \(a\) +denoting the exponents and the named component \(n\) denoting the dimension.

+

One can specify cones in any order if strict_cone_order +is set to FALSE in the call to clarabel() but +one has to ensure that parameter types are strictly specified for the +values, e.g. 5L for integers, 0. for reals +etc.

+
+
+

5. Control parameters

Clarabel has a number of parameters that control its behavior, including verbosity, time limits, and tolerances; see help on clarabel_control(). As an example, in the last problem, we can reduce the number of iterations.

-
+
 P <- Matrix::Matrix(2 * c(0, 0, 0, 1), nrow = 2, ncol = 2, sparse = TRUE)
 P <- as(P, "symmetricMatrix") # P needs to be a symmetric matrix
 q <- c(0, 0)
@@ -433,8 +538,9 @@ 

4. Control parameters= list(max_iter = 3)) ## Reduced number of iterations cat(sprintf("Solution status, description: = (%d, %s)\n", s$status, solver_status_descriptions()[s$status])) -#> Solution status, description: = (5, Solver terminated with a solution (reduced accuracy)) -cat(sprintf("Solution (x1, x2) = (%f, %f)\n", s$x[1], s$x[2])) +#> Solution status, description: = (5, Solver terminated with a solution (reduced accuracy))

+
+cat(sprintf("Solution (x1, x2) = (%f, %f)\n", s$x[1], s$x[2]))
 #> Solution (x1, x2) = (1.000000, -0.999998)

Note the different status, which should always be checked in code.

@@ -458,7 +564,7 @@

4. Control parameters

-

Site built with pkgdown 2.0.7.

+

Site built with pkgdown 2.0.9.

diff --git a/docs/articles/index.html b/docs/articles/index.html index 28bc6291..d764454e 100644 --- a/docs/articles/index.html +++ b/docs/articles/index.html @@ -1,5 +1,5 @@ -Articles • clarabelArticles • clarabel @@ -17,7 +17,7 @@ clarabel - 0.5.1 + 0.9.0
@@ -65,7 +65,7 @@

All vignettes

-

Site built with pkgdown 2.0.7.

+

Site built with pkgdown 2.0.9.

diff --git a/docs/authors.html b/docs/authors.html index cadd3396..dbea3b7e 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -1,5 +1,5 @@ -Authors and Citation • clarabelAuthors and Citation • clarabel @@ -17,7 +17,7 @@ clarabel - 0.5.1 + 0.9.0 @@ -47,7 +47,7 @@
@@ -80,15 +80,49 @@

Citation

-

Goulart P, Chen Y (2021). -“Clarabel: A Library for Optimization and Control.” -https://oxfordcontrol.github.io/ClarabelDocs/stable/. +

Goulart PJ, Chen Y (2024). +“Clarabel: An interior-point solver for conic programs with quadratic objectives.” +doi:10.48550/arXiv.2405.12762, 2405.12762.

-
@Misc{,
-  title = {Clarabel: A Library for Optimization and Control},
-  author = {Paul Goulart and Yuwen Chen},
-  year = {2021},
-  url = {https://oxfordcontrol.github.io/ClarabelDocs/stable/},
+    
@Misc{Clarabel_2024,
+  title = {Clarabel: An interior-point solver for conic programs with quadratic objectives},
+  author = {Paul J. Goulart and Yuwen Chen},
+  year = {2024},
+  eprint = {2405.12762},
+  archiveprefix = {arXiv},
+  doi = {10.48550/arXiv.2405.12762},
+  primaryclass = {math.OC},
+}
+

Garstka M, Cannon M, Goulart P (2020). +“A clique graph based merging strategy for decomposable SDPs.” +In IFAC-PapersOnLine, volume 53 number 2, 7355-7361. +doi:10.1016/j.ifacol.2020.12.1255, 21th IFAC World Congress. +

+
@InProceedings{Garstka_2020,
+  author = {Michael Garstka and Mark Cannon and Paul Goulart},
+  title = {A clique graph based merging strategy for decomposable {SDPs}},
+  booktitle = {IFAC-PapersOnLine},
+  year = {2020},
+  note = {21th IFAC World Congress},
+  number = {2},
+  pages = {7355-7361},
+  volume = {53},
+  doi = {10.1016/j.ifacol.2020.12.1255},
+  issn = {2405-8963},
+  journal = {IFAC-PapersOnLine},
+}
+

Chen Y, Goulart P (2023). +“An Efficient IPM Implementation for A Class of Nonsymmetric Cones.” +doi:10.48550/arXiv.2305.12275, 2305.12275. +

+
@Misc{chen2023efficient,
+  title = {An Efficient IPM Implementation for A Class of Nonsymmetric Cones},
+  author = {Yuwen Chen and Paul Goulart},
+  year = {2023},
+  eprint = {2305.12275},
+  archiveprefix = {arXiv},
+  doi = {10.48550/arXiv.2305.12275},
+  primaryclass = {math.OC},
 }
@@ -102,7 +136,7 @@

Citation

-

Site built with pkgdown 2.0.7.

+

Site built with pkgdown 2.0.9.

diff --git a/docs/index.html b/docs/index.html index 0ccd660d..6cd4db8a 100644 --- a/docs/index.html +++ b/docs/index.html @@ -12,13 +12,13 @@ - + - + Changelog • clarabelChangelog • clarabel @@ -17,7 +17,7 @@ clarabel - 0.5.1 + 0.9.0 @@ -51,7 +51,15 @@

Changelog

- + +
  • Synced up to version 0.9.0 of Clarabel.rs
  • +
  • Added all applicable tests from Clarabel.rs
  • +
  • Updated documentation of cone specification
  • +
  • Updated vignette
  • +
  • Switched to savvy from rextendr
  • +
+
+
  • Clarabel now supports semidefinite programs (syncing up to version 0.5.1 of Clarabel.rs)
  • Added tests
  • Updated documentation of cone specification
  • @@ -74,7 +82,7 @@
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 6d318038..22e1fe88 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -1,7 +1,7 @@ -pandoc: 2.17.1.1 -pkgdown: 2.0.7 +pandoc: '3.2' +pkgdown: 2.0.9 pkgdown_sha: ~ articles: clarabel: clarabel.html -last_built: 2023-06-24T16:30Z +last_built: 2024-06-22T18:01Z diff --git a/docs/reference/clarabel-package.html b/docs/reference/clarabel-package.html index 00aef776..164bf3a5 100644 --- a/docs/reference/clarabel-package.html +++ b/docs/reference/clarabel-package.html @@ -1,5 +1,5 @@ -Interface to Clarabel solver implemented in Rust. — clarabel-package • clarabelInterface to Clarabel solver implemented in Rust. — clarabel-package • clarabel @@ -17,7 +17,7 @@ clarabel - 0.5.1 + 0.9.0 @@ -52,10 +52,16 @@

Interface to Clarabel solver implemented in Rust.

-

Clarabel is a versatile interior point solver for convex programs using a new homogeneous embedding. It solves solves linear programs (LPs), quadratic programs (QPs), second-order cone programs (SOCPs), and problems with exponential and power cone constraints. For quadratic objectives, unlike interior point solvers based on the standard homogeneous self-dual embedding (HSDE) model, Clarabel handles quadratic objective without requiring any epigraphical reformulation of its objective function. It can therefore be significantly faster than other HSDE-based solvers for problems with quadratic objective functions. Infeasible problems are detected using a homogeneous embedding technique. See https://oxfordcontrol.github.io/ClarabelDocs/stable/.

+

Clarabel is a versatile interior point solver for convex programs using a new homogeneous embedding. It solves solves linear programs (LPs), quadratic programs (QPs), second-order cone programs (SOCPs), and problems with exponential and power cone constraints. For quadratic objectives, unlike interior point solvers based on the standard homogeneous self-dual embedding (HSDE) model, Clarabel handles quadratic objective without requiring any epigraphical reformulation of its objective function. It can therefore be significantly faster than other HSDE-based solvers for problems with quadratic objective functions. Infeasible problems are detected using a homogeneous embedding technique. See https://clarabel.org/stable/.

+

Author

Balasubramanian Narasimhan, Paul Goulart, Yuwen Chen

@@ -73,7 +79,7 @@

Author

-

Site built with pkgdown 2.0.7.

+

Site built with pkgdown 2.0.9.

diff --git a/docs/reference/clarabel.html b/docs/reference/clarabel.html index ecd01bb7..4eb8646b 100644 --- a/docs/reference/clarabel.html +++ b/docs/reference/clarabel.html @@ -1,5 +1,5 @@ -Interface to 'Clarabel', an interior point conic solver — clarabel • clarabelInterface to 'Clarabel', an interior point conic solver — clarabel • clarabel clarabel - 0.5.1 + 0.9.0 @@ -149,7 +149,7 @@

Clarabel can solve

Cone Parameters

-

The table below shows the cone parameter specifications

ParameterTypeLengthDescription
zinteger\(1\)number of primal zero cones (dual free cones),
which corresponds to the primal equality constraints
linteger\(1\)number of linear cones (non-negative cones)
qinteger\(\geq1\)vector of second-order cone sizes
sinteger\(\geq1\)vector of positive semidefinite cone sizes
epinteger\(1\)number of primal exponential cones
pnumeric\(\geq1\)vector of primal power cone parameters

+

The table below shows the cone parameter specifications. Mathematical definitions are in the vignette.

ParameterTypeLengthDescription
zinteger\(1\)number of primal zero cones (dual free cones),
which corresponds to the primal equality constraints
linteger\(1\)number of linear cones (non-negative cones)
qinteger\(\ge 1\)vector of second-order cone sizes
sinteger\(\ge 1\)vector of positive semidefinite cone sizes
epinteger\(1\)number of primal exponential cones
pnumeric\(\ge 1\)vector of primal power cone parameters
gplist\(\ge 1\)list of named lists of two items, a : a numeric vector of at least 2 exponent terms in the product summing to 1, and n : an integer dimension of generalized power cone parameters

When the parameter strict_cone_order is FALSE, one can specify @@ -159,7 +159,8 @@

Cone Parameterslist(z = 2L, l = 2L, q = 2L, z = 3L, q = 3L), or, list(z1 = 2L, l1 = 2L, q1 = 2L, zb = 3L, qx = 3L), indicating a zero cone of size 2, followed by a linear cone of size 2, followed by a second-order cone of size 2, followed by a zero cone of size 3, and finally a second-order -cone of size 3.

+cone of size 3. Generalized power cones parameters have to specified as named lists, e.g., list(z = 2L, gp1 = list(a = c(0.3, 0.7), n = 3L), gp2 = list(a = c(0.5, 0.5), n = 1L)).

+

Note that when strict_cone_order = FALSE, types of cone parameters such as integers, reals have to be explicit since the parameters are directly passed to the Rust interface with no sanity checks.!

See also

@@ -186,14 +187,17 @@

Examples

#> $obj_val #> [1] 1 #> +#> $obj_val_dual +#> [1] 1 +#> #> $status #> [1] 2 #> #> $solve_time -#> [1] 4.5391e-05 +#> [1] 0.000111267 #> #> $iterations -#> [1] 3 +#> [1] 0 #> #> $r_prim #> [1] 0 @@ -201,6 +205,53 @@

Examples

#> $r_dual #> [1] 0 #> +#> $info +#> $info$μ +#> [1] 1 +#> +#> $info$sigma +#> [1] 1 +#> +#> $info$step_length +#> [1] 0 +#> +#> $info$cost_primal +#> [1] 1 +#> +#> $info$cost_dual +#> [1] 1 +#> +#> $info$res_primal +#> [1] 0 +#> +#> $info$res_dual +#> [1] 0 +#> +#> $info$res_primal_inf +#> [1] 1 +#> +#> $info$res_dual_inf +#> [1] 1.414214 +#> +#> $info$gap_abs +#> [1] 0 +#> +#> $info$gap_rel +#> [1] 0 +#> +#> $info$ktratio +#> [1] 1 +#> +#> $info$solve_time +#> [1] 0.000111267 +#> +#> $info$iterations +#> [1] 0 +#> +#> $info$status +#> [1] 2 +#> +#>
@@ -216,7 +267,7 @@

Examples

-

Site built with pkgdown 2.0.7.

+

Site built with pkgdown 2.0.9.

diff --git a/docs/reference/clarabel_control.html b/docs/reference/clarabel_control.html index 58d19b0c..e2fe59a7 100644 --- a/docs/reference/clarabel_control.html +++ b/docs/reference/clarabel_control.html @@ -1,5 +1,5 @@ -Control parameters with default values and types in parenthesis — clarabel_control • clarabelControl parameters with default values and types in parenthesis — clarabel_control • clarabel @@ -17,7 +17,7 @@ clarabel - 0.5.1 + 0.9.0 @@ -93,7 +93,11 @@

Control parameters with default values and types in parenthesis

iterative_refinement_abstol = 1e-12, iterative_refinement_max_iter = 10L, iterative_refinement_stop_ratio = 5, - presolve_enable = TRUE + presolve_enable = TRUE, + chordal_decomposition_enable = FALSE, + chordal_decomposition_merge_method = c("none", "parent_child", "clique_graph"), + chordal_decomposition_compact = FALSE, + chordal_decomposition_complete_dual = FALSE ) @@ -246,6 +250,22 @@

Arguments

presolve_enable

whether to enable presolvle (TRUE)

+ +
chordal_decomposition_enable
+

whether to enable chordal decomposition for SDPs (FALSE)

+ + +
chordal_decomposition_merge_method
+

chordal decomposition merge method, one of 'none', 'parent_child' or 'clique_graph', for SDPs ('none')

+ + +
chordal_decomposition_compact
+

a boolean flag for SDPs indicating whether to assemble decomposed system in compact form for SDPs (FALSE)

+ + +
chordal_decomposition_complete_dual
+

a boolean flag indicating complete PSD dual variables after decomposition for SDPs

+

Value

@@ -266,7 +286,7 @@

Value

-

Site built with pkgdown 2.0.7.

+

Site built with pkgdown 2.0.9.

diff --git a/docs/reference/index.html b/docs/reference/index.html index e3f4e88e..92e4da2d 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -1,5 +1,5 @@ -Function reference • clarabelFunction reference • clarabel @@ -17,7 +17,7 @@ clarabel - 0.5.1 + 0.9.0 @@ -54,10 +54,6 @@

All functions

-

clarabel-package

- -

Interface to Clarabel solver implemented in Rust.

-

clarabel()

Interface to 'Clarabel', an interior point conic solver

@@ -82,7 +78,7 @@

All functions
-

Site built with pkgdown 2.0.7.

+

Site built with pkgdown 2.0.9.

diff --git a/docs/reference/solver_status_descriptions.html b/docs/reference/solver_status_descriptions.html index 97ae90df..e35bdc98 100644 --- a/docs/reference/solver_status_descriptions.html +++ b/docs/reference/solver_status_descriptions.html @@ -1,5 +1,5 @@ -Return the solver status description as a named character vector — solver_status_descriptions • clarabelReturn the solver status description as a named character vector — solver_status_descriptions • clarabel @@ -17,7 +17,7 @@ clarabel - 0.5.1 + 0.9.0 @@ -88,7 +88,7 @@

Examples

-

Site built with pkgdown 2.0.7.

+

Site built with pkgdown 2.0.9.

diff --git a/inst/AUTHORS b/inst/AUTHORS index 74b66dbc..98306e93 100644 --- a/inst/AUTHORS +++ b/inst/AUTHORS @@ -1,18 +1,18 @@ The authors of the dependency Rust crates: - -aho-corasick (version 1.0.2): +aho-corasick (version 1.1.3): Andrew Gallant amd (version 0.2.2): -autocfg (version 1.1.0): +autocfg (version 1.3.0): Josh Stone -blas-src (version 0.9.0): +blas-src (version 0.10.0): Balasubramanian Narasimhan Ivan Ukhov Jed Brown + Michael Zietz Stefan Kroboth Toshiki Teramura bluss @@ -30,10 +30,13 @@ blas (version 0.22.0): Mason Smith Toshiki Teramura +cc (version 1.0.99): + Alex Crichton + cfg-if (version 1.0.0): Alex Crichton -clarabel (version 0.5.1): +clarabel (version 0.9.0): Paul Goulart darling_core (version 0.14.4): @@ -63,43 +66,37 @@ derive_builder (version 0.11.2): Jan-Erik Rediger Ted Driggs -enum_dispatch (version 0.3.11): - Anton Lazarev +either (version 1.12.0): + bluss -extendr-api (version 0.4.0): - andy-thomason - Thomas Down - Mossa Merhi Reimert - Claus O. Wilke - Hiroaki Yutani - Ilia A. Kosenkov - Michael Milton - -extendr-engine (version 0.4.0): - andy-thomason - Thomas Down - Mossa Merhi Reimert - Claus O. Wilke - Hiroaki Yutani - Ilia A. Kosenkov +enum_dispatch (version 0.3.13): + Anton Lazarev -extendr-macros (version 0.4.0): - andy-thomason - Thomas Down - Mossa Merhi Reimert - Claus O. Wilke - Hiroaki Yutani - Ilia A. Kosenkov +equivalent (version 1.0.1): + fnv (version 1.0.7): Alex Crichton +hashbrown (version 0.14.5): + Amanieu d'Antras + ident_case (version 1.0.1): Ted Driggs -lapack-src (version 0.9.0): +indexmap (version 2.2.6): + + +itertools (version 0.11.0): + bluss + +itoa (version 1.0.11): + David Tolnay + +lapack-src (version 0.10.0): Balasubramanian Narasimhan Ivan Ukhov + Michael Zietz Mitsutoshi Aoe Stefan Kroboth Toshiki Teramura @@ -122,48 +119,70 @@ lapack (version 0.19.0): lazy_static (version 1.4.0): Marvin Löbel -libc (version 0.2.146): +libc (version 0.2.155): The Rust Project Developers -libR-sys (version 0.4.0): - andy-thomason - Thomas Down - Mossa Merhi Reimert - Claus O. Wilke - Ilia A. Kosenkov - Hiroaki Yutani - -memchr (version 2.5.0): +memchr (version 2.7.4): Andrew Gallant bluss -num-complex (version 0.4.3): +num-complex (version 0.4.6): The Rust Project Developers -num-traits (version 0.2.15): +num-traits (version 0.2.19): The Rust Project Developers -once_cell (version 1.18.0): +once_cell (version 1.19.0): Aleksey Kladov -paste (version 1.0.12): - David Tolnay - -proc-macro2 (version 1.0.60): +proc-macro2 (version 1.0.85): David Tolnay Alex Crichton -quote (version 1.0.28): +quote (version 1.0.36): David Tolnay r-src (version 0.1.0): Balasubramanian Narasimhan -regex-syntax (version 0.7.2): +regex-automata (version 0.4.7): The Rust Project Developers + Andrew Gallant + +regex-syntax (version 0.8.4): + The Rust Project Developers + Andrew Gallant -regex (version 1.8.4): +regex (version 1.10.5): The Rust Project Developers + Andrew Gallant + +ryu (version 1.0.18): + David Tolnay + +savvy-bindgen (version 0.6.4): + Hiroaki Yutani + +savvy-ffi (version 0.6.4): + Hiroaki Yutani + +savvy-macro (version 0.6.4): + Hiroaki Yutani + +savvy (version 0.6.4): + Hiroaki Yutani + +serde_derive (version 1.0.203): + Erick Tryzelaar + David Tolnay + +serde_json (version 1.0.117): + Erick Tryzelaar + David Tolnay + +serde (version 1.0.203): + Erick Tryzelaar + David Tolnay strsim (version 0.10.0): Danny Guo @@ -171,14 +190,14 @@ strsim (version 0.10.0): syn (version 1.0.109): David Tolnay -syn (version 2.0.18): +syn (version 2.0.66): David Tolnay -thiserror-impl (version 1.0.40): +thiserror-impl (version 1.0.61): David Tolnay -thiserror (version 1.0.40): +thiserror (version 1.0.61): David Tolnay -unicode-ident (version 1.0.9): +unicode-ident (version 1.0.12): David Tolnay diff --git a/inst/CITATION b/inst/CITATION index 6be4a0c6..e63c16e2 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -1,8 +1,51 @@ bibentry(bibtype = "Misc", - title = "Clarabel: A Library for Optimization and Control", - author = c(person(given = "Paul", family = "Goulart"), - person(given = "Yuwen", family = "Chen")), - year = "2021", - url = "https://oxfordcontrol.github.io/ClarabelDocs/stable/", - header = "To cite Clarabel in publications use:" -) + key = "Clarabel_2024", + title = "Clarabel: An interior-point solver for conic programs with quadratic objectives", + author = c(person(given = c("Paul", "J."), + family = "Goulart"), + person(given = "Yuwen", + family = "Chen")), + year = "2024", + eprint = "2405.12762", + archiveprefix = "arXiv", + doi = "10.48550/arXiv.2405.12762", + primaryclass = "math.OC", + header = "To cite Clarabel in publications use:" + ) + +bibentry(bibtype = "InProceedings", + key = "Garstka_2020", + author = c(person(given = "Michael", + family = "Garstka"), + person(given = "Mark", + family = "Cannon"), + person(given = "Paul", + family = "Goulart")), + title = "A clique graph based merging strategy for decomposable {SDPs}", + booktitle = "IFAC-PapersOnLine", + year = "2020", + note = "21th IFAC World Congress", + number = "2", + pages = "7355-7361", + volume = "53", + doi = "10.1016/j.ifacol.2020.12.1255", + issn = "2405-8963", + journal = "IFAC-PapersOnLine", + header = "To cite Clarabel to solve SDPs with decomposable structure use:" + ) + +bibentry(bibtype = "Misc", + key = "chen2023efficient", + title = "An Efficient IPM Implementation for A Class of Nonsymmetric Cones", + author = c(person(given = "Yuwen", + family = "Chen"), + person(given = "Paul", + family = "Goulart")), + year = "2023", + eprint = "2305.12275", + archiveprefix = "arXiv", + doi = "10.48550/arXiv.2305.12275", + primaryclass = "math.OC", + header = "To cite Clarabel to solve problems with generalized power cones use:" + ) + diff --git a/inst/misc/code_gen.R b/inst/misc/code_gen.R new file mode 100644 index 00000000..f2f99412 --- /dev/null +++ b/inst/misc/code_gen.R @@ -0,0 +1,24 @@ +library(clarabel) +# Function to match on types using savvy +f <- function(name, type = c("integer", "double", "logical", "character")) { + type <- match.arg(type) + c( + sprintf('"%s" => ', name), + sprintf(' match typed_value {'), + switch(type, + integer = sprintf(' TypedSexp::Integer(i) => settings.%s = i.as_slice()[0],', name), + double = sprintf(' TypedSexp::Real(f) => settings.%s = f.as_slice()[0],', name), + logical = sprintf(' TypedSexp::Logical(b) => settings.%s = b.as_slice_raw()[0] != 0,', name), + character = sprintf(' TypedSexp::String(s) => if let Some(result) = s.to_vec().get(0) {\\n + settings.%s = result.to_string();\\n },\\n _ => (),', name) + ), + "_ => ()," + "}," + ) +} + +s <- clarabel_control() +sn <- names(s) +st <- sapply(s, typeof) +out <- lapply(seq_along(st), function(i) f(sn[i], st[i])) +writeLines(unlist(out), "update_settings_snippet.rs") diff --git a/inst/tinytest/test_api_dimension.R b/inst/tinytest/test_api_dimension.R new file mode 100644 index 00000000..df8c8fcb --- /dev/null +++ b/inst/tinytest/test_api_dimension.R @@ -0,0 +1,51 @@ + +api_dim_check_data <- function() { + P <- matrix(0.0, nrow = 4, ncol = 4) + q <- numeric(4) + A <- matrix(0.0, nrow = 6, ncol = 4) + b <- numeric(6) + cones = list(z = 1L, l = 5L) + list(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +} + +## This example should work because dimensions are +## all compatible. All following checks vary one +## of these sizes to test dimension checks +problem_data <- api_dim_check_data() +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) + +## Check that bad P dimension is detected +problem_data$P <- matrix(0.0, nrow = 3, ncol = 3) +expect_error(solution <- do.call(clarabel, problem_data)) + +## Check that bad A rows are detected +problem_data <- api_dim_check_data() +problem_data$A <- matrix(0.0, nrow = 5, ncol = 4) +expect_error(solution <- do.call(clarabel, problem_data)) + +## Check that bad A columns are detected +problem_data <- api_dim_check_data() +problem_data$A <- matrix(0.0, nrow = 6, ncol = 3) +expect_error(solution <- do.call(clarabel, problem_data)) + +## Check that P not square is detected +problem_data <- api_dim_check_data() +problem_data$P <- matrix(0.0, nrow = 5, ncol = 4) +expect_error(solution <- do.call(clarabel, problem_data)) + +## Check that repeated specification of cone types with strict_cone_order TRUE fails +problem_data <- api_dim_check_data() +problem_data$cones <- list(z = 1L, l = 2L, l = 4L) +expect_error(solution <- do.call(clarabel, problem_data)) + +## Check that the cone dimensions mismatch with matrices is detected +problem_data <- api_dim_check_data() +problem_data$cones <- list(z = 1L, l = 6L) +expect_error(solution <- do.call(clarabel, problem_data)) + +## Check that repeated specification of cone types succeeds with strict_cone_order FALSE +problem_data <- api_dim_check_data() +problem_data$cones <- list(z = 1L, l = 2L, l = 3L) +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) diff --git a/inst/tinytest/test_api_dimension.R~ b/inst/tinytest/test_api_dimension.R~ new file mode 100644 index 00000000..17c0b11d --- /dev/null +++ b/inst/tinytest/test_api_dimension.R~ @@ -0,0 +1,57 @@ + +api_dim_check_data <- function() { + P <- matrix(0.0, nrow = 4, ncol = 4) + q <- numeric(4) + A <- matrix(0.0, nrow = 6, ncol = 4) + b <- numeric(6) + cones = list(z = 1L, l = 5L) + list(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +} + +## This example should work because dimensions are +## all compatible. All following checks vary one +## of these sizes to test dimension checks +problem_data <- api_dim_check_data() +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) + +## Check that bad P dimension is detected +problem_data$P <- matrix(0.0, nrow = 3, ncol = 3) +solution <- do.call(clarabel, problem_data) +expect_error(clarabel(P = P, q = q, A = A, b = b, cones = cones, control = settings)) + +## Check that bad A rows are detected +problem_data <- api_dim_check_data() +problem_data$A <- matrix(0.0, nrow = 5, ncol = 4) +solution <- do.call(clarabel, problem_data) +expect_error(clarabel(P = P, q = q, A = A, b = b, cones = cones, control = settings)) + +## Check that bad A columns are detected +problem_data <- api_dim_check_data() +problem_data$A <- matrix(0.0, nrow = 6, ncol = 3) +solution <- do.call(clarabel, problem_data) +expect_error(clarabel(P = P, q = q, A = A, b = b, cones = cones, control = settings)) + +## Check that P not square is detected +problem_data <- api_dim_check_data() +problem_data$P <- matrix(0.0, nrow = 5, ncol = 4) +solution <- do.call(clarabel, problem_data) +expect_error(clarabel(P = P, q = q, A = A, b = b, cones = cones, control = settings)) + +## Check that repeated specification of cone types with strict_cone_order TRUE fails +problem_data <- api_dim_check_data() +problem_data$cones <- list(z = 1L, l = 2L, l = 4L) +solution <- do.call(clarabel, problem_data) +expect_error(clarabel(P = P, q = q, A = A, b = b, cones = cones, control = settings)) + +## Check that the cone dimensions mismatch with matrices is detected +problem_data <- api_dim_check_data() +problem_data$cones <- list(z = 1L, l = 6L) +solution <- do.call(clarabel, problem_data) +expect_error(clarabel(P = P, q = q, A = A, b = b, cones = cones, control = settings)) + +## Check that repeated specification of cone types succeeds with strict_cone_order FALSE +problem_data <- api_dim_check_data() +problem_data$cones <- list(z = 1L, l = 2L, l = 3L) +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) diff --git a/inst/tinytest/test_basic_eq_constrained.R b/inst/tinytest/test_basic_eq_constrained.R new file mode 100644 index 00000000..7734dd5c --- /dev/null +++ b/inst/tinytest/test_basic_eq_constrained.R @@ -0,0 +1,59 @@ +eq_constrained_A1 <- function() { + ## A = + ##[ 0. 1. 1.; + ## 0. 1. -1.] + Matrix::sparseMatrix( + dims = c(2, 3), ## m x n + p = c(0, 0, 2, 4), ##colptr + i = c(0, 1, 0, 1), ##rowval + x = c(1., 1., 1., -1.), ##nzva; + index1 = FALSE + ) +} + +eq_constrained_A2 <- function() { + ## A = [ + ## 0 1.0 1.0; + ## 0 1.0 -1.0; + ##1.0 2.0 -1.0l + ##2.0 -1.0 3.0l + ##] + Matrix::sparseMatrix( + dims = c(4, 3), ## m x n + p = c(0, 2, 6, 10), ##colptr + i = c(2, 3, 0, 1, 2, 3, 0, 1, 2, 3), ##rowval + x = c(1., 2., 1., 1., 2., -1., 1., -1., -1., 3.), ##nzva; + index1 = FALSE + ) +} + +## test_eq_constrained_feasible +P <- diag(3) +q <- numeric(3) +A <- eq_constrained_A1() ## <- two constraints +b <- c(2., 0.) +cones <- list(z = 2) +solution <- clarabel(P = P, q = q, A = A, b = b, cones = cones) +refsol <- c(0., 1., 1.) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +expect_true(l2_dist(solution$x, refsol) <= 1e-6) + + +## test_eq_constrained_primal_infeasible +P <- diag(3) +q <- numeric(3) +A <- eq_constrained_A2() ## <- 4 constraints, 3 vars +b <- rep(1.0, 4) +cones <- list(z = 4) +solution <- clarabel(P = P, q = q, A = A, b = b, cones = cones) +expect_equal(status_codes[[solution$status]], status_codes[["PrimalInfeasible"]]) + + +## test_eq_constrained_dual_infeasible +P <- diag(3); P[1, 1] <- 0; +q <- rep(1.0, 3) +A <- eq_constrained_A1() ## <- two constraints +b <- c(2., 0.) +cones <- list(z = 2L) +solution <- clarabel(P = P, q = q, A = A, b = b, cones = cones) +expect_equal(status_codes[[solution$status]], status_codes[["DualInfeasible"]]) diff --git a/inst/tinytest/test_basic_expcone.R b/inst/tinytest/test_basic_expcone.R new file mode 100644 index 00000000..99363f09 --- /dev/null +++ b/inst/tinytest/test_basic_expcone.R @@ -0,0 +1,66 @@ +basic_expcone_data <- function() { + ## produces data for the following exponential cone problem + ## max x + ## s.t. y * exp(x / y) <= z + ## y == 1, z == exp(5) + + P <- matrix(0, nrow = 3, ncol = 3) + q <- c(-1., 0., 0.) + A1 <- -Matrix::Diagonal(3) + b1 <- rep(0., 3) + A2 <- Matrix::sparseMatrix( + dims = c(2, 3), ## m x n + p = c(0, 0, 1, 2), ##colptr + i = c(0, 1), ##rowval + x = c(1., 1.), ##nzval + index1 = FALSE + ) + b2 <- c(1.0, exp(5.0)) + A <- rbind(A1, A2) + b <- c(b1, b2) + + cones <- list(ep = 1L, z = 2L) + list(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +} + +## test_expcone_feasible +## solve the following exponential cone problem +## max x +## s.t. y * exp(x / y) <= z +## y == 1, z == exp(5) +## +## This is just the default problem data above + +problem_data <- basic_expcone_data() +refsol <- c(5.0, 1., exp(5.0)) +refobj <- -5.0 +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +expect_true(abs(solution$info$cost_primal - refobj) <= 1e-6) +expect_true(l2_dist(solution$x, refsol) <= 1e-6) + +## test_expcone_primal_infeasible +## solve the following exponential cone problem +## max x +## s.t. y * exp(x / y) <= z +## y == 1, z == -1 +## +## Same as default, but last element of b is different +problem_data <- basic_expcone_data() +problem_data$b[5] = -1. +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["PrimalInfeasible"]]) + +## test_expcone_dual_infeasible +## solve the following exponential cone problem +## max x +## s.t. y * exp(x / y) <= z +## +## Same as default, but no equality constraint +P <- NULL +q <- c(-1., 0., 0.) +A <- -diag(3) +b <- numeric(3) +cones <- list(ep = 1L) +solution <- clarabel(P = P, q = q, A = A, b = b, cones = cones) +expect_equal(status_codes[[solution$status]], status_codes[["DualInfeasible"]]) diff --git a/inst/tinytest/test_basic_genpowcone.R b/inst/tinytest/test_basic_genpowcone.R new file mode 100644 index 00000000..6bc1f93c --- /dev/null +++ b/inst/tinytest/test_basic_genpowcone.R @@ -0,0 +1,41 @@ +## test_powcone +## solve the following power cone problem +## max x1^0.6 y^0.4 + x2^0.1 +## s.t. x1, y, x2 >= 0 +## x1 + 2y + 3x2 == 3 +## which is equivalent to +## max z1 + z2 +## s.t. (x1, y, z1) in K_pow(0.6) +## (x2, 1, z2) in K_pow(0.1) +## x1 + 2y + 3x2 == 3 + +## x = (x1, y, z1, x2, y2, z2) + +n <- 6; +P <- Matrix::sparseMatrix(i = integer(0), j = integer(0), x = numeric(0), dims = c(n, n)) +q <- c(0., 0., -1., 0., 0., -1.) + +## (x1, y, z1) in K_pow(0.6) +## (x2, y2, z2) in K_pow(0.1) + +A1 <- -diag(n) +b1 <- numeric(n) +cones1 <- list(gp1 = list(a = c(0.6, 0.4), n = 1L), + gp2 = list(a = c(0.1, 0.9), n = 1L) + ) + +## x1 + 2y + 3x2 == 3 +## y2 == 1 + +A2 <- Matrix::sparseMatrix(i = c(rep(1, 3), 2), j = c(1, 2, 4, 5), x = c(1.0, 2.0, 3.0, 1.0), + dims = c(2, 6)) +b2 <- c(3., 1.) +cones2 <- list(z = 2L) +A <- rbind(A1, A2) +b <- c(b1, b2) +cones <- c(cones1, cones2) + +solution <- clarabel(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +refobj <- -1.8458; +expect_true(abs(solution$info$cost_primal - refobj) <= 1e-3) diff --git a/inst/tinytest/test_basic_lp.R b/inst/tinytest/test_basic_lp.R new file mode 100644 index 00000000..b6b9b156 --- /dev/null +++ b/inst/tinytest/test_basic_lp.R @@ -0,0 +1,53 @@ + +#source("inst/tinytest/preamble.R") + +basic_lp_data <- function() { + P <- Matrix::sparseMatrix(i = integer(0), j = integer(0), x = numeric(0), dims = c(3, 3)) + #P <- NULL + I1 <- diag(3) + I2 <- -diag(3) + A <- rbind(I1, I2) + A <- A * 2 + q <- c(3., -2., 1.) + b <- rep(1.0, 6) + cones <- list(l1 = 3L, l2 = 3L) + list(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +} + +## test_lp_feasible +problem_data <- basic_lp_data(); +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +refsol <- c(-0.5, 0.5, -0.5) +expect_true(l2_dist(solution$x, refsol) <= 1e-8) +refobj <- -3.0 +expect_true(abs(solution$obj_val - refobj) <= 1e-8) +expect_true(abs(solution$obj_val_dual- refobj) <= 1e-8) + + +## test_lp_primal_infeasible +problem_data <- basic_lp_data(); +problem_data$b[1] <- -1; problem_data$b[4] <- -1; +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["PrimalInfeasible"]]) +expect_true(is.nan(solution$obj_val)) +expect_true(is.nan(solution$obj_val_dual)) + +## test_lp_dual_infeasible +problem_data <- basic_lp_data(); +problem_data$A[4, 1] <- 1. ##swap lower bound on first variable to redundant upper bound +problem_data$q <- c(1., 0., 0.) +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["DualInfeasible"]]) +expect_true(is.nan(solution$obj_val)) +expect_true(is.nan(solution$obj_val_dual)) + +## test_lp_dual_infeasible_ill_cond() { +problem_data <- basic_lp_data(); +problem_data$A[1, 1] <- .Machine$double.eps +problem_data$A[4, 1] <- 0 +problem_data$q <- c(1., 0., 0.) +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["DualInfeasible"]]) +expect_true(is.nan(solution$obj_val)) +expect_true(is.nan(solution$obj_val_dual)) diff --git a/inst/tinytest/test_basic_powcone.R b/inst/tinytest/test_basic_powcone.R new file mode 100644 index 00000000..9d196755 --- /dev/null +++ b/inst/tinytest/test_basic_powcone.R @@ -0,0 +1,38 @@ + +#source("inst/tinytest/preamble.R") + +## test_powcone +## solve the following power cone problem +## max x1^0.6 y^0.4 + x2^0.1 +## s.t. x1, y, x2 >= 0 +## x1 + 2y + 3x2 == 3 +## which is equivalent to +## max z1 + z2 +## s.t. (x1, y, z1) in K_pow(0.6) +## (x2, 1, z2) in K_pow(0.1) +## x1 + 2y + 3x2 == 3 + +## x = (x1, y, z1, x2, y2, z2) +n <- 6; +P <- Matrix::sparseMatrix(i = integer(0), j = integer(0), x = numeric(0), dims = c(n, n)) +q <- c(0., 0., -1., 0., 0., -1.) + +## (x1, y, z1) in K_pow(0.6) +## (x2, y2, z2) in K_pow(0.1) +A1 <- -diag(n) +b1 <- numeric(n) +cones1 <- list(p1 = 0.6, p2 = 0.1) + +## x1 + 2y + 3x2 == 3 +## y2 == 1 +A2 <- Matrix::sparseMatrix(i = c(rep(1, 3), 2), j = c(1, 2, 4, 5), x = c(1.0, 2.0, 3.0, 1.0), + dims = c(2, 6)) +b2 <- c(3., 1.) +cones2 <- list(z = 2L) +A <- rbind(A1, A2) +b <- c(b1, b2) +cones <- c(cones1, cones2) +solution <- clarabel(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +refobj <- -1.8458; +expect_true(abs(solution$info$cost_primal - refobj) <= 1e-3) diff --git a/inst/tinytest/test_basic_psd.R b/inst/tinytest/test_basic_psd.R new file mode 100644 index 00000000..ad05c5d1 --- /dev/null +++ b/inst/tinytest/test_basic_psd.R @@ -0,0 +1,97 @@ +## Semidefinite cone example + +## Basic example + +basic_sdp_data <- function() { + ## problem will be 3x3, so upper triangle + ## of problem data has 6 entries + P <- as(Matrix::Diagonal(6), "symmetricMatrix") + A <- diag(6) + q <- numeric(6) + b <- c(-3., 1., 4., 1., 2., 5.) + cones <- list(s = 3L) + list(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +} + +## test_sdp_feasible +problem_data <- basic_sdp_data() +refsol <- c(-3.0729833267361095, 0.3696004167288786, -0.022226685581313674, + 0.31441213129613066, -0.026739700851545107, -0.016084530571308823) +refobj <- 4.840076866013861 +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +expect_true(abs(solution$info$cost_primal - refobj) <= 1e-6) +expect_true(l2_dist(solution$x, refsol) <= 1e-6) + +## empty_sdp_cone +problem_data <- basic_sdp_data() +problem_data$cones <- c(problem_data$cones, list(s = 0L)) +refsol <- c(-3.0729833267361095, 0.3696004167288786, -0.022226685581313674, + 0.31441213129613066, -0.026739700851545107, -0.016084530571308823) +refobj <- 4.840076866013861 +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +expect_true(abs(solution$info$cost_primal - refobj) <= 1e-6) +expect_true(l2_dist(solution$x, refsol) <= 1e-6) + +## sdp_primal_infeasible +problem_data <- basic_sdp_data() +problem_data$A <- rbind(problem_data$A, -problem_data$A) +problem_data$b <- c(problem_data$b, numeric(length(problem_data$b))) +problem_data$cones <- c(problem_data$cones, problem_data$cones) +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["PrimalInfeasible"]]) + +## Example from scs + +#' Return an vectorization of symmetric matrix using the upper triangular part, +#' still in column order. +#' @param S a symmetric matrix +#' @return vector of values +vec <- function(S) { + n <- nrow(S) + sqrt2 <- sqrt(2.0) + upper_tri <- upper.tri(S, diag = FALSE) + S[upper_tri] <- S[upper_tri] * sqrt2 + S[upper.tri(S, diag = TRUE)] +} + +#' Return the symmetric matrix from the [vec] vectorization +#' @param v a vector +#' @return a symmetric matrix +mat <- function(v) { + n <- (sqrt(8 * length(v) + 1) - 1) / 2 + sqrt2 <- sqrt(2.0) + S <- matrix(0, n, n) + upper_tri <- upper.tri(S, diag = TRUE) + S[upper_tri] <- v / sqrt2 + S <- S + t(S) + diag(S) <- diag(S) / sqrt(2) + S +} + +q <- c(1, -1, 1) # objective: x_1 - x2 + x_3 +A11 <- matrix(c(-7, -11, -11, 3), nrow = 2) +A12 <- matrix(c(7, -18, -18, 8), nrow = 2) +A13 <- matrix(c(-2, -8, -8, 1), nrow = 2) + +A21 <- matrix(c(-21, -11, 0, -11, 10, 8, 0, 8, 5), nrow = 3) +A22 <- matrix(c(0, 10, 16, 10, -10, -10, 16, -10, 3), nrow = 3) +A23 <- matrix(c(-5, 2, -17, 2, -6, 8, -17, 8, 6), nrow = 3) + +B1 <- matrix(c(33, -9, -9, 26), nrow = 2) +B2 <- matrix(c(14, 9, 40, 9, 91, 10, 40, 10, 15), nrow = 3) + +A <- rbind( + cbind(vec(A11), vec(A12), vec(A13)), # first psd constraint + cbind(vec(A21), vec(A22), vec(A23)) # second psd constraint +) +b <- c(vec(B1), vec(B2)) # stack both psd constraints +cones <- list(s = c(2, 3)) # cone dimensions +sol <- clarabel(A = A, b = b, q = q, cone = cones) + +expect_equal(sol$status, 2L) +expect_equal(sol$x, + c(-0.367751526880478, 1.89833324679172, -0.887460228128519), + tolerance = 1e-5) + diff --git a/inst/tinytest/test_basic_qp.R b/inst/tinytest/test_basic_qp.R new file mode 100644 index 00000000..bb6dc6d0 --- /dev/null +++ b/inst/tinytest/test_basic_qp.R @@ -0,0 +1,82 @@ + +## Simple QP test (old remnant) + +P <- Matrix::Matrix(2 * c(3, 0, 0, 2), nrow = 2, ncol = 2, sparse = TRUE) +P <- as(P, "symmetricMatrix") # P needs to be a symmetric matrix +q <- c(-1, -4) +A <- Matrix::Matrix(c(1, 1, 0, -1, 0, -2, 0, 1, 0, -1), ncol = 2, sparse = TRUE) +b <- c(0, 1, 1, 1, 1) +cones <- list(z = 1L, l = 4L) ## 1 equality and 4 inequalities, in order +sol <- clarabel(A = A, b = b, q = q, P = P, cones = cones) +expect_equal(sol$status, 2L) +expect_equal(sol$x, + c(0.428571428198843, 0.214285714099422), + tolerance = 1e-7) + +basic_qp_data <- function() { + P <- matrix(c(4, 1, 1, 2), byrow = TRUE, ncol = 2) + A <- matrix(c(1., 1, 1, 0, 0, 1), byrow = TRUE, ncol = 2) + A <- rbind(-A, A) + q <- c(1., 1.) + b <- c(-1., 0., 0., 1., 0.7, 0.7) + cones = list(l1 = 3L, l2 = 3L) + list(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +} + +basic_qp_data_dual_inf <- function() { + P <- matrix(1, nrow = 2, ncol = 2) + A <- matrix(c(1., 1, 1, 0), byrow = TRUE, ncol = 2) + q <- c(1., -1.) + b <- c(1., 1.) + cones = list(l = 2L) + list(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +} + +## test_qp_univariate NEEDS FIXING because of univariate matrices!! The S3 methods for +## coercion don't work. +## P <- as.matrix(1.0) +## q <- 0.0 +## A <- P +## b <- 1.0 +## cones <- list(l = 1L) +## solution <- clarabel(P = P, q = q, A = A, b = b, cones = cones) +## expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +## expect_true(abs(solution$x) < 1e-6) +## expect_true(abs(solution$obj_val) < 1e-6) +## expect_true(abs(solution$obj_val_dual) < 1e-6) + +## test_qp_feasible +problem_data <- basic_qp_data(); +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +refsol <- c(0.3, 0.7) +expect_true(l2_dist(solution$x, refsol) <= 1e-6) +refobj <- 1.8800000298331538; +expect_true(abs(solution$obj_val - refobj) <= 1e-6) +expect_true(abs(solution$obj_val_dual- refobj) <= 1e-6) + +## test_qp_primal_infeasible +problem_data <- basic_qp_data(); +problem_data$b[1] <- -1; problem_data$b[4] <- -1; +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["PrimalInfeasible"]]) +expect_true(is.nan(solution$obj_val)) +expect_true(is.nan(solution$obj_val_dual)) + +## test_qp_dual_infeasible +problem_data <- basic_qp_data_dual_inf() +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["DualInfeasible"]]) +expect_true(is.nan(solution$obj_val)) +expect_true(is.nan(solution$obj_val_dual)) + +## test_qp_dual_infeasible_ill_cond +problem_data <- basic_qp_data_dual_inf() +problem_data$A <- Matrix::sparseMatrix(i = rep(0, 2), p = c(0, 1, 2), x = rep(1.0, 2), index1 = FALSE, dims = c(1, 2)) +problem_data$cones <- list(l = 1L) +problem_data$b <- 1.0 +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["DualInfeasible"]]) +expect_true(is.nan(solution$obj_val)) +expect_true(is.nan(solution$obj_val_dual)) + diff --git a/inst/tinytest/test_basic_socp.R b/inst/tinytest/test_basic_socp.R new file mode 100644 index 00000000..3e117829 --- /dev/null +++ b/inst/tinytest/test_basic_socp.R @@ -0,0 +1,64 @@ +#source("inst/tinytest/preamble.R") + +basic_socp_data <- function() { + ## P matrix data taken from corresponding Julia unit test. + ## These nzvals form a 3x3 positive definite matrix + + P <- Matrix::sparseMatrix( + p = c(0, 3, 6, 9), ## colptr + i = c(0, 1, 2, 0, 1, 2, 0, 1, 2), ## rowval + x = c( ## nzval + 1.4652521089139698, + 0.6137176286085666, + -1.1527861771130112, + 0.6137176286085666, + 2.219109946678485, + -1.4400420548730628, + -1.1527861771130112, + -1.4400420548730628, + 1.6014483534926371 + ), + index1 = FALSE, + dims = c(3, 3) + ) + ## A = [2I;-2I;I] + A <- rbind(2*diag(3), -2*diag(3), diag(3)) + q <- c(0.1, -2.0, 1.0) + b <- c(1., 1., 1., 1., 1., 1., 0., 0., 0.) + cones <- list(l1 = 3L, l2 = 3L, q = 3L) + list(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +} + +## test_socp_feasible +problem_data <- basic_socp_data() +refsol <- c(-0.5, 0.435603, -0.245459) +refobj <- -8.4590e-01 +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +expect_true(l2_dist(solution$x, refsol) <= 1e-4) +expect_true(abs(solution$obj_val - refobj) <= 1e-4) +expect_true(abs(solution$obj_val_dual - refobj) <= 1e-4) + +## test_socp_infeasible +problem_data <- basic_socp_data() +## make the cone constraint unsatisfiable +problem_data$b[7] <- -10. +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["PrimalInfeasible"]]) +expect_true(is.nan(solution$obj_val)) +expect_true(is.nan(solution$obj_val_dual)) + + +## Basic SOCP + +P <- Matrix::Matrix(2 * c(0, 0, 0, 1), nrow = 2, ncol = 2, sparse = TRUE) +P <- as(P, "symmetricMatrix") # P needs to be a symmetric matrix +q <- c(0, 0) +A <- Matrix::Matrix(c(0, -2.0, 0, 0, 0, 1.0), nrow = 3, ncol = 2, sparse = TRUE) +b <- c(1, -2, -2) +cones <- list(q = 3L) +sol <- clarabel(A = A, b = b, q = q, P = P, cones = cones) +expect_equal(sol$status, 2L) +expect_equal(sol$x, + c(1, -1), + tolerance = 1e-7) diff --git a/inst/tinytest/test_basic_unconstrained.R b/inst/tinytest/test_basic_unconstrained.R new file mode 100644 index 00000000..aa358693 --- /dev/null +++ b/inst/tinytest/test_basic_unconstrained.R @@ -0,0 +1,19 @@ +## test_unconstrained_feasible +P <- diag(3) +q <- c(1.0, 2.0, -3.0) +A <- matrix(0, nrow = 0, ncol = 3) ## <- no constraints +b <- numeric(0) +cones <- list() +solution <- clarabel(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +refsol <- -q +expect_true(l2_dist(solution$x, refsol) <= 1e-6) + +## test_unconstrained_dual_infeasible +P <- matrix(0, nrow = 3, ncol = 3) +q <- c(1.0, 0, 0) +A <- matrix(0, nrow = 0, ncol = 3) ## <- no constraints +b <- numeric(0) +cones <- list() +solution <- clarabel(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +expect_equal(status_codes[[solution$status]], status_codes[["DualInfeasible"]]) diff --git a/inst/tinytest/test_mixed_conic.R b/inst/tinytest/test_mixed_conic.R new file mode 100644 index 00000000..6ac1d690 --- /dev/null +++ b/inst/tinytest/test_mixed_conic.R @@ -0,0 +1,23 @@ +## test_mixed_conic_feasible +## solves a problem with a mix of symmetric and asymmetric +## cones. This exercises the barrier methods and unit +## initializations of the symmetric cones + +n <- 3 +P <- diag(3) +q <- c(1., 1., 1.) +## put a 3 dimensional vector into the composition of multiple +## cones, all with b = 0 on the RHS +cones <- list(z = 3L, + l = 3L, + q = 3L, + p = 0.5, + ep = 1L + ) +A <- rbind(diag(3), diag(3)) +A <- rbind(A, A, diag(3)) +b <- numeric(5 * n) +solution <- clarabel(P = P, q = q, A = A, b = b, cones = cones) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +expect_true(abs(solution$info$cost_primal - 0.) <= 1e-8) + diff --git a/inst/tinytest/test_presolve.R b/inst/tinytest/test_presolve.R new file mode 100644 index 00000000..88515e16 --- /dev/null +++ b/inst/tinytest/test_presolve.R @@ -0,0 +1,35 @@ +## This is not the complete suite since we don't expose the solver object! +## TBD +presolve_test_data <- function() { + n <- 3 + P <- diag(n) + A <- 2 * rbind(P, -diag(n)) + q <- c(3., -2., 1.) + b <- rep(1.0, 2 * n) + cones <- list(l1 = 3L, l2 = 3L) + list(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +} + +## test_presolve_single_unbounded +problem_data <- presolve_test_data() +problem_data$b[4] <- 1e30 +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +expect_equal(solution$z[4], 0.) + +## test_presolve_completely_redundant_cone +problem_data <- presolve_test_data() +problem_data$b[1:3] <- 1e30 +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +expect_equal(solution$z[1:3], numeric(3)) +refsol <- c(-0.5, 2., -0.5) +expect_true(l2_dist(solution$x, refsol) <= 1e-6) + +##test_presolve_every_constraint_redundant +problem_data <- presolve_test_data() +problem_data$b <- rep(1e30, length(problem_data$b)) +solution <- do.call(clarabel, problem_data) +expect_equal(status_codes[[solution$status]], status_codes[["Solved"]]) +expect_true(l2_dist(solution$x, -problem_data$q) <= 1e-6) + diff --git a/inst/tinytest/test_psd.R b/inst/tinytest/test_psd.R deleted file mode 100644 index be409338..00000000 --- a/inst/tinytest/test_psd.R +++ /dev/null @@ -1,72 +0,0 @@ -## Semidefinite cone example - -## Basic example - -## problem will be 3x3, so upper triangle -## of problem data has 6 entries - -P <- as(Matrix::Diagonal(6), "symmetricMatrix") -A <- diag(6) -c <- numeric(6) -b <- c(-3., 1., 4., 1., 2., 5.) -cones <- list(s = 3L) -sol <- clarabel(P = P, A = A, b = b, q = c, cone = cones) -expect_equal(sol$status, 2L) -expect_equal(sol$x, - c(-3.0729833267361095, 0.3696004167288786, -0.022226685581313674, - 0.31441213129613066, -0.026739700851545107, -0.016084530571308823), - tolerance = 1e-7) - -## Example from scs - -#' Return an vectorization of symmetric matrix using the upper triangular part, -#' still in column order. -#' @param S a symmetric matrix -#' @return vector of values -vec <- function(S) { - n <- nrow(S) - sqrt2 <- sqrt(2.0) - upper_tri <- upper.tri(S, diag = FALSE) - S[upper_tri] <- S[upper_tri] * sqrt2 - S[upper.tri(S, diag = TRUE)] -} - -#' Return the symmetric matrix from the [vec] vectorization -#' @param v a vector -#' @return a symmetric matrix -mat <- function(v) { - n <- (sqrt(8 * length(v) + 1) - 1) / 2 - sqrt2 <- sqrt(2.0) - S <- matrix(0, n, n) - upper_tri <- upper.tri(S, diag = TRUE) - S[upper_tri] <- v / sqrt2 - S <- S + t(S) - diag(S) <- diag(S) / sqrt(2) - S -} - -q <- c(1, -1, 1) # objective: x_1 - x2 + x_3 -A11 <- matrix(c(-7, -11, -11, 3), nrow = 2) -A12 <- matrix(c(7, -18, -18, 8), nrow = 2) -A13 <- matrix(c(-2, -8, -8, 1), nrow = 2) - -A21 <- matrix(c(-21, -11, 0, -11, 10, 8, 0, 8, 5), nrow = 3) -A22 <- matrix(c(0, 10, 16, 10, -10, -10, 16, -10, 3), nrow = 3) -A23 <- matrix(c(-5, 2, -17, 2, -6, 8, -17, 8, 6), nrow = 3) - -B1 <- matrix(c(33, -9, -9, 26), nrow = 2) -B2 <- matrix(c(14, 9, 40, 9, 91, 10, 40, 10, 15), nrow = 3) - -A <- rbind( - cbind(vec(A11), vec(A12), vec(A13)), # first psd constraint - cbind(vec(A21), vec(A22), vec(A23)) # second psd constraint -) -b <- c(vec(B1), vec(B2)) # stack both psd constraints -cones <- list(s = c(2, 3)) # cone dimensions -sol <- clarabel(A = A, b = b, q = q, cone = cones) - -expect_equal(sol$status, 2L) -expect_equal(sol$x, - c(-0.367751526880478, 1.89833324679172, -0.887460228128519), - tolerance = 1e-7) - diff --git a/inst/tinytest/test_qp.R b/inst/tinytest/test_qp.R deleted file mode 100644 index 76fe482b..00000000 --- a/inst/tinytest/test_qp.R +++ /dev/null @@ -1,15 +0,0 @@ -## Simple QP test - -P <- Matrix::Matrix(2 * c(3, 0, 0, 2), nrow = 2, ncol = 2, sparse = TRUE) -P <- as(P, "symmetricMatrix") # P needs to be a symmetric matrix -q <- c(-1, -4) -A <- Matrix::Matrix(c(1, 1, 0, -1, 0, -2, 0, 1, 0, -1), ncol = 2, sparse = TRUE) -b <- c(0, 1, 1, 1, 1) -cones <- list(z = 1L, l = 4L) ## 1 equality and 4 inequalities, in order -sol <- clarabel(A = A, b = b, q = q, P = P, cones = cones) -expect_equal(sol$status, 2L) -expect_equal(sol$x, - c(0.428571428198843, 0.214285714099422), - tolerance = 1e-7) - - diff --git a/inst/tinytest/test_sdp_chordal.R b/inst/tinytest/test_sdp_chordal.R new file mode 100644 index 00000000..9fcc425c --- /dev/null +++ b/inst/tinytest/test_sdp_chordal.R @@ -0,0 +1,87 @@ +SQRT_2 <- sqrt(2.0) + +sdp_chordal_data <- function() { + P <- Matrix::sparseMatrix(i = integer(0), j = integer(0), x = numeric(0), dims = c(8, 8)) + q <- c(-1.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0) + A <- Matrix::sparseMatrix( + dims = c(28, 8), + p = c(0, 1, 4, 5, 8, 9, 10, 13, 16), + i = c(24, 7, 10, 22, 8, 12, 15, 25, 9, 13, 18, 21, 26, 0, 23, 27), + x = c( + -1.0, + -SQRT_2, + -1.0, + -1.0, + -SQRT_2, + -SQRT_2, + -1.0, + -1.0, + -SQRT_2, + -SQRT_2, + -SQRT_2, + -1.0, + -1.0, + -1.0, + -1.0, + -1.0 + ), + index1 = FALSE + ) + b <- c( + 0.0, + 3.0, + 2. * SQRT_2, + 2.0, + SQRT_2, + SQRT_2, + 3.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0 + ) + cones <- list(l = 1L, s = 6L, p1 = 1.0 / 3.0, p2 = 0.5) + list(P = P, q = q, A = A, b = b, cones = cones, strict_cone_order = FALSE) +} + +##test_sdp_chordal() { +problem_data <- sdp_chordal_data() +problem_data$control <- clarabel_control( + chordal_decomposition_enable = TRUE, + chordal_decomposition_compact = TRUE, + chordal_decomposition_complete_dual = TRUE, + chordal_decomposition_merge_method = "clique_graph", + max_iter = 50L) + +for (compact in c(FALSE, TRUE)) { + for (complete_dual in c(FALSE, TRUE)) { + for (merge_method in c("clique_graph", "parent_child", "none")) { + problem_data$control <- clarabel_control( + chordal_decomposition_enable = TRUE, + chordal_decomposition_compact = compact, + chordal_decomposition_complete_dual = complete_dual, + chordal_decomposition_merge_method = merge_method, + max_iter = 50L) + solution <- do.call(clarabel, problem_data) + print(expect_equal(status_codes[[solution$status]], status_codes[["Solved"]])) + } + } +} + diff --git a/inst/tinytest/test_socp.R b/inst/tinytest/test_socp.R deleted file mode 100644 index 215518d5..00000000 --- a/inst/tinytest/test_socp.R +++ /dev/null @@ -1,13 +0,0 @@ -## Basic SOCP - -P <- Matrix::Matrix(2 * c(0, 0, 0, 1), nrow = 2, ncol = 2, sparse = TRUE) -P <- as(P, "symmetricMatrix") # P needs to be a symmetric matrix -q <- c(0, 0) -A <- Matrix::Matrix(c(0, -2.0, 0, 0, 0, 1.0), nrow = 3, ncol = 2, sparse = TRUE) -b <- c(1, -2, -2) -cones <- list(q = 3L) -sol <- clarabel(A = A, b = b, q = q, P = P, cones = cones) -expect_equal(sol$status, 2L) -expect_equal(sol$x, - c(1, -1), - tolerance = 1e-7) diff --git a/man/clarabel-package.Rd b/man/clarabel-package.Rd index 57c366ff..4e965952 100644 --- a/man/clarabel-package.Rd +++ b/man/clarabel-package.Rd @@ -5,9 +5,17 @@ \alias{clarabel-package} \title{Interface to Clarabel solver implemented in Rust.} \description{ -Clarabel is a versatile interior point solver for convex programs using a new homogeneous embedding. It solves solves linear programs (LPs), quadratic programs (QPs), second-order cone programs (SOCPs), and problems with exponential and power cone constraints. For quadratic objectives, unlike interior point solvers based on the standard homogeneous self-dual embedding (HSDE) model, Clarabel handles quadratic objective without requiring any epigraphical reformulation of its objective function. It can therefore be significantly faster than other HSDE-based solvers for problems with quadratic objective functions. Infeasible problems are detected using a homogeneous embedding technique. See \url{https://oxfordcontrol.github.io/ClarabelDocs/stable/}. +Clarabel is a versatile interior point solver for convex programs using a new homogeneous embedding. It solves solves linear programs (LPs), quadratic programs (QPs), second-order cone programs (SOCPs), and problems with exponential and power cone constraints. For quadratic objectives, unlike interior point solvers based on the standard homogeneous self-dual embedding (HSDE) model, Clarabel handles quadratic objective without requiring any epigraphical reformulation of its objective function. It can therefore be significantly faster than other HSDE-based solvers for problems with quadratic objective functions. Infeasible problems are detected using a homogeneous embedding technique. See \url{https://clarabel.org/stable/}. +} +\seealso{ +Useful links: +\itemize{ + \item \url{https://oxfordcontrol.github.io/clarabel-r/} + \item Report bugs at \url{https://github.com/oxfordcontrol/clarabel-r/issues} +} + } \author{ Balasubramanian Narasimhan, Paul Goulart, Yuwen Chen } -\keyword{package} +\keyword{internal} diff --git a/man/clarabel.Rd b/man/clarabel.Rd index 88372bec..45ce49a0 100644 --- a/man/clarabel.Rd +++ b/man/clarabel.Rd @@ -67,16 +67,17 @@ problems with any combination of cones, defined by the parameters listed in \dQuote{Cone Parameters} below } } \subsection{Cone Parameters}{ -The table below shows the cone parameter specifications +The table below shows the cone parameter specifications. Mathematical definitions are in the vignette. \tabular{rllll}{ \tab \bold{Parameter} \tab \bold{Type} \tab \bold{Length} \tab \bold{Description} \cr -\tab \code{z} \tab integer \tab \eqn{1} \tab number of primal zero cones (dual free cones), \cr -\tab \tab \tab \tab which corresponds to the primal equality constraints \cr -\tab \code{l} \tab integer \tab \eqn{1} \tab number of linear cones (non-negative cones) \cr -\tab \code{q} \tab integer \tab \eqn{\geq1} \tab vector of second-order cone sizes \cr -\tab \code{s} \tab integer \tab \eqn{\geq1} \tab vector of positive semidefinite cone sizes \cr -\tab \code{ep} \tab integer \tab \eqn{1} \tab number of primal exponential cones \cr -\tab \code{p} \tab numeric \tab \eqn{\geq1} \tab vector of primal power cone parameters +\tab \code{z} \tab integer \tab \eqn{1} \tab number of primal zero cones (dual free cones), \cr +\tab \tab \tab \tab which corresponds to the primal equality constraints \cr +\tab \code{l} \tab integer \tab \eqn{1} \tab number of linear cones (non-negative cones) \cr +\tab \code{q} \tab integer \tab \eqn{\ge 1} \tab vector of second-order cone sizes \cr +\tab \code{s} \tab integer \tab \eqn{\ge 1} \tab vector of positive semidefinite cone sizes \cr +\tab \code{ep} \tab integer \tab \eqn{1} \tab number of primal exponential cones \cr +\tab \code{p} \tab numeric \tab \eqn{\ge 1} \tab vector of primal power cone parameters \cr +\tab \code{gp} \tab list \tab \eqn{\ge 1} \tab list of named lists of two items, \code{a} : a numeric vector of at least 2 exponent terms in the product summing to 1, and \code{n} : an integer dimension of generalized power cone parameters } } When the parameter \code{strict_cone_order} is \code{FALSE}, one can specify @@ -86,7 +87,9 @@ argument in such a case should be a named list with names matching and so on. For example, either of the following would be valid: \code{list(z = 2L, l = 2L, q = 2L, z = 3L, q = 3L)}, or, \code{list(z1 = 2L, l1 = 2L, q1 = 2L, zb = 3L, qx = 3L)}, indicating a zero cone of size 2, followed by a linear cone of size 2, followed by a second-order cone of size 2, followed by a zero cone of size 3, and finally a second-order -cone of size 3. +cone of size 3. Generalized power cones parameters have to specified as named lists, e.g., \code{list(z = 2L, gp1 = list(a = c(0.3, 0.7), n = 3L), gp2 = list(a = c(0.5, 0.5), n = 1L))}. + +\emph{Note that when \code{strict_cone_order = FALSE}, types of cone parameters such as integers, reals have to be explicit since the parameters are directly passed to the Rust interface with no sanity checks.!} } \examples{ A <- matrix(c(1, 1), ncol = 1) diff --git a/man/clarabel_control.Rd b/man/clarabel_control.Rd index 433ebb0b..2ce8cff8 100644 --- a/man/clarabel_control.Rd +++ b/man/clarabel_control.Rd @@ -41,7 +41,11 @@ clarabel_control( iterative_refinement_abstol = 1e-12, iterative_refinement_max_iter = 10L, iterative_refinement_stop_ratio = 5, - presolve_enable = TRUE + presolve_enable = TRUE, + chordal_decomposition_enable = FALSE, + chordal_decomposition_merge_method = c("none", "parent_child", "clique_graph"), + chordal_decomposition_compact = FALSE, + chordal_decomposition_complete_dual = FALSE ) } \arguments{ @@ -118,6 +122,14 @@ clarabel_control( \item{iterative_refinement_stop_ratio}{iterative refinement stalling tolerance (\code{5.0})} \item{presolve_enable}{whether to enable presolvle (\code{TRUE})} + +\item{chordal_decomposition_enable}{whether to enable chordal decomposition for SDPs (\code{FALSE})} + +\item{chordal_decomposition_merge_method}{chordal decomposition merge method, one of \code{'none'}, \code{'parent_child'} or \code{'clique_graph'}, for SDPs (\code{'none'})} + +\item{chordal_decomposition_compact}{a boolean flag for SDPs indicating whether to assemble decomposed system in \emph{compact} form for SDPs (\code{FALSE})} + +\item{chordal_decomposition_complete_dual}{a boolean flag indicating complete PSD dual variables after decomposition for SDPs} } \value{ a list containing the control parameters. diff --git a/src/.gitignore b/src/.gitignore index 0acc02d8..780abd8e 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -2,3 +2,5 @@ *.so *.dll target +.cargo +Makevars \ No newline at end of file diff --git a/src/Makevars b/src/Makevars deleted file mode 100644 index 16f956ef..00000000 --- a/src/Makevars +++ /dev/null @@ -1,36 +0,0 @@ -# FLAG to indicate use of bundled source; use yes/no -VENDORING = yes -VENDOR_SRC = ./rust/vendor.tar.xz -TARGET_DIR = ./rust/target -LIBDIR = $(TARGET_DIR)/release -STATLIB = $(LIBDIR)/libclarabel.a -PKG_LIBS = -L$(LIBDIR) -lclarabel $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) - -all: C_clean - -$(SHLIB): $(STATLIB) - -CARGOTMP = $(CURDIR)/rust/.cargo - -$(STATLIB): - # vendoring (Note: to avoid NOTE of "Found the following hidden files and - # directories", .cargo needs to be created here) - if [ "$(VENDORING)" = "yes" ]; then \ - mkdir -p $(CARGOTMP); \ - cp ./rust/cargo_vendor_config.toml $(CARGOTMP)/config.toml; \ - $(TAR) --extract --xz -f $(VENDOR_SRC) -C ./rust ; \ - export CARGO_HOME=$(CARGOTMP); \ - export PATH="$(PATH):$(HOME)/.cargo/bin"; \ - cargo build --lib --release --offline --manifest-path=./rust/Cargo.toml --target-dir $(TARGET_DIR); \ - fi && \ - if [ "$(VENDORING)" = "no" ]; then \ - cargo build --lib --verbose --release --manifest-path=./rust/Cargo.toml --target-dir $(TARGET_DIR); \ - fi - rm -Rf $(CARGOTMP) - rm -Rf $(LIBDIR)/build - -C_clean: - rm -Rf $(SHLIB) $(STATLIB) $(OBJECTS) - -clean: - rm -Rf $(SHLIB) $(STATLIB) $(OBJECTS) $(CARGOTMP) $(TARGET_DIR) rust/vendor rust/Cargo.lock diff --git a/src/Makevars.in b/src/Makevars.in new file mode 100644 index 00000000..e23da092 --- /dev/null +++ b/src/Makevars.in @@ -0,0 +1,53 @@ +TARGET = @TARGET@ +VENDORING = yes +VENDOR_SRC = ./rust/vendor.tar.xz +VENDOR_DIR = ./rust/vendor + +# catch DEBUG envvar, which is passed from pkgbuild::compile_dll() +PROFILE = $(subst x,release,$(subst truex,dev,$(DEBUG)x)) + +TARGET_DIR = $(CURDIR)/rust/target +LIBDIR = $(TARGET_DIR)/release +STATLIB = $(LIBDIR)/libclarabel.a +PKG_LIBS = -L$(LIBDIR) -lclarabel $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) + +CARGOTMP = $(CURDIR)/rust/.cargo +CARGO_BUILD_ARGS = --jobs 2 --verbose --lib --manifest-path=./rust/Cargo.toml --target-dir $(TARGET_DIR) + +all: C_clean + +$(SHLIB): $(STATLIB) + +$(STATLIB): + # In some environments, ~/.cargo/bin might not be included in PATH, so we need + # to set it here to ensure cargo can be invoked. It is appended to PATH and + # therefore is only used if cargo is absent from the user's PATH. + export PATH="$(PATH):$(HOME)/.cargo/bin" && \ + if [ "$(TARGET)" != "wasm32-unknown-emscripten" ]; then \ + if [ "$(VENDORING)" = "yes" ]; then \ + mkdir -p $(CARGOTMP); \ + cp ./rust/cargo_vendor_config.toml $(CARGOTMP)/config.toml; \ + $(TAR) --extract --xz -f $(VENDOR_SRC) -C ./rust ; \ + export CARGO_HOME=$(CARGOTMP); \ + export PATH="$(PATH):$(HOME)/.cargo/bin"; \ + cargo build $(CARGO_BUILD_ARGS) --release --offline; \ + else \ + cargo build $(CARGO_BUILD_ARGS); \ + fi \ + else \ + export CC="$(CC)" && \ + export CFLAGS="$(CFLAGS)" && \ + export CARGO_PROFILE_DEV_PANIC="abort" && \ + export CARGO_PROFILE_RELEASE_PANIC="abort" && \ + cargo +nightly build $(CARGO_BUILD_ARGS) -Zbuild-std=panic_abort,std; \ + fi + rm -Rf $(CARGOTMP) + rm -Rf $(VENDOR_DIR) + rm -Rf $(LIBDIR)/build + +C_clean: + rm -Rf $(SHLIB) $(STATLIB) $(OBJECTS) + +clean: + rm -Rf $(SHLIB) $(STATLIB) $(OBJECTS) rust/target + diff --git a/src/Makevars.ucrt b/src/Makevars.ucrt deleted file mode 100644 index 17b153e3..00000000 --- a/src/Makevars.ucrt +++ /dev/null @@ -1,5 +0,0 @@ -# Rtools42 doesn't have the linker in the location that cargo expects, so we -# need to overwrite it via configuration. -CARGO_LINKER = x86_64-w64-mingw32.static.posix-gcc.exe - -include Makevars.win diff --git a/src/Makevars.win b/src/Makevars.win index d4b44bbd..98ca41b1 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -1,37 +1,43 @@ -# For Windows, we use vendoring by default! -TARGET = $(subst 64,x86_64,$(subst 32,i686,$(WIN)))-pc-windows-gnu +TARGET = x86_64-pc-windows-gnu -TARGET_DIR = ./rust/target +VENDOR_DIR = ./rust/vendor + +# catch DEBUG envvar, which is passed from pkgbuild::compile_dll() +PROFILE = $(subst x,release,$(subst truex,dev,$(DEBUG)x)) + +TARGET_DIR = $(CURDIR)/rust/target LIBDIR = $(TARGET_DIR)/$(TARGET)/release STATLIB = $(LIBDIR)/libclarabel.a -PKG_LIBS = -L$(LIBDIR) -lclarabel -lntdll -lws2_32 -ladvapi32 -luserenv -lbcrypt $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) +PKG_LIBS = -L$(LIBDIR) -lclarabel -lws2_32 -ladvapi32 -luserenv -lbcrypt -lntdll $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) + +CARGOTMP = $(CURDIR)/rust/.cargo + +# Rtools doesn't have the linker in the location that cargo expects, so we need +# to overwrite it via configuration. +CARGO_LINKER = x86_64-w64-mingw32.static.posix-gcc.exe all: C_clean $(SHLIB): $(STATLIB) -CARGOTMP = $(CURDIR)/rust/.cargo - $(STATLIB): + # When the GNU toolchain is used (i.e. on CRAN), -lgcc_eh is specified for + # building proc-macro2, but Rtools doesn't contain libgcc_eh. This isn't used + # in actual, but we need this tweak to please the compiler. + mkdir -p $(LIBDIR)/libgcc_mock && touch $(LIBDIR)/libgcc_mock/libgcc_eh.a mkdir -p $(CARGOTMP) && cp ./rust/cargo_vendor_config.toml $(CARGOTMP)/config.toml - # `rustc` adds `-lgcc_eh` flags to the compiler, but Rtools' GCC doesn't have - # `libgcc_eh` due to the compilation settings. So, in order to please the - # compiler, we need to add empty `libgcc_eh` to the library search paths. - # - # For more details, please refer to - # https://github.com/r-windows/rtools-packages/blob/2407b23f1e0925bbb20a4162c963600105236318/mingw-w64-gcc/PKGBUILD#L313-L316 - mkdir -p $(TARGET_DIR)/libgcc_mock - touch $(TARGET_DIR)/libgcc_mock/libgcc_eh.a $(TAR) --extract --xz -f ./rust/vendor.tar.xz -C ./rust - export CARGO_HOME="$(CARGOTMP)"; \ - export CARGO_TARGET_X86_64_PC_WINDOWS_GNU_LINKER="$(CARGO_LINKER)"; \ - export LIBRARY_PATH="$${LIBRARY_PATH};$(CURDIR)/$(TARGET_DIR)/libgcc_mock"; \ - cargo build --target=$(TARGET) --lib --release --offline --manifest-path=./rust/Cargo.toml --target-dir $(TARGET_DIR) + + export CARGO_HOME="$(CARGOTMP)" && \ + export CARGO_TARGET_X86_64_PC_WINDOWS_GNU_LINKER="$(CARGO_LINKER)" && \ + export LIBRARY_PATH="$${LIBRARY_PATH};$(LIBDIR)/libgcc_mock" && \ + cargo build --jobs 2 --target $(TARGET) --lib --offline --release --manifest-path ./rust/Cargo.toml --target-dir $(TARGET_DIR) rm -Rf $(CARGOTMP) + rm -Rf $(VENDOR_DIR) rm -Rf $(LIBDIR)/build C_clean: rm -Rf $(SHLIB) $(STATLIB) $(OBJECTS) clean: - rm -Rf $(SHLIB) $(STATLIB) $(OBJECTS) $(TARGET_DIR) rust/vendor rust/Cargo.lock + rm -Rf $(SHLIB) $(STATLIB) $(OBJECTS) rust/target diff --git a/src/entrypoint.c b/src/entrypoint.c deleted file mode 100644 index d2f6767d..00000000 --- a/src/entrypoint.c +++ /dev/null @@ -1,8 +0,0 @@ -// We need to forward routine registration from C to Rust -// to avoid the linker removing the static library. - -void R_init_clarabel_extendr(void *dll); - -void R_init_clarabel(void *dll) { - R_init_clarabel_extendr(dll); -} diff --git a/src/init.c b/src/init.c new file mode 100644 index 00000000..49573d88 --- /dev/null +++ b/src/init.c @@ -0,0 +1,54 @@ + +#include +#include +#include + +#include "rust/api.h" + +static uintptr_t TAGGED_POINTER_MASK = (uintptr_t)1; + +SEXP handle_result(SEXP res_) { + uintptr_t res = (uintptr_t)res_; + + // An error is indicated by tag. + if ((res & TAGGED_POINTER_MASK) == 1) { + // Remove tag + SEXP res_aligned = (SEXP)(res & ~TAGGED_POINTER_MASK); + + // Currently, there are two types of error cases: + // + // 1. Error from Rust code + // 2. Error from R's C API, which is caught by R_UnwindProtect() + // + if (TYPEOF(res_aligned) == CHARSXP) { + // In case 1, the result is an error message that can be passed to + // Rf_errorcall() directly. + Rf_errorcall(R_NilValue, "%s", CHAR(res_aligned)); + } else { + // In case 2, the result is the token to restart the + // cleanup process on R's side. + R_ContinueUnwind(res_aligned); + } + } + + return (SEXP)res; +} + +SEXP savvy_clarabel_solve__impl(SEXP m, SEXP n, SEXP Ai, SEXP Ap, SEXP Ax, SEXP b, SEXP q, SEXP Pi, SEXP Pp, SEXP Px, SEXP cone_spec, SEXP r_settings) { + SEXP res = savvy_clarabel_solve__ffi(m, n, Ai, Ap, Ax, b, q, Pi, Pp, Px, cone_spec, r_settings); + return handle_result(res); +} + + +static const R_CallMethodDef CallEntries[] = { + {"savvy_clarabel_solve__impl", (DL_FUNC) &savvy_clarabel_solve__impl, 12}, + {NULL, NULL, 0} +}; + +void R_init_clarabel(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); + + // Functions for initialzation, if any. + +} diff --git a/src/rust/Cargo.lock b/src/rust/Cargo.lock index e782b4ef..aa441dd2 100644 --- a/src/rust/Cargo.lock +++ b/src/rust/Cargo.lock @@ -4,9 +4,9 @@ version = 3 [[package]] name = "aho-corasick" -version = "1.0.2" +version = "1.1.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "43f6cb1bf222025340178f382c426f13757b2960e89779dfcb319c32542a5a41" +checksum = "8e60d3430d3a69478ad0993f19238d2df97c507009a52b3c10addcd7f6bcb916" dependencies = [ "memchr", ] @@ -22,9 +22,9 @@ dependencies = [ [[package]] name = "autocfg" -version = "1.1.0" +version = "1.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" +checksum = "0c4b4d0bd25bd0b74681c0ad21497610ce1b7c91b1022cd21c80c6fbdd9476b0" [[package]] name = "blas" @@ -39,9 +39,9 @@ dependencies = [ [[package]] name = "blas-src" -version = "0.9.0" +version = "0.10.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "aa443ee19b4cde6cdbd49043eb8964f9dd367b6d98d67f04395958ebfa28f39d" +checksum = "b95e83dc868db96e69795c0213143095f03de9dd3252f205d4ac716e4076a7e0" dependencies = [ "r-src", ] @@ -55,6 +55,12 @@ dependencies = [ "libc", ] +[[package]] +name = "cc" +version = "1.0.99" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "96c51067fd44124faa7f870b4b1c969379ad32b2ba805aa959430ceaa384f695" + [[package]] name = "cfg-if" version = "1.0.0" @@ -63,19 +69,19 @@ checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" [[package]] name = "clarabel" -version = "0.5.1" +version = "0.9.0" dependencies = [ - "clarabel 0.5.1 (registry+https://github.com/rust-lang/crates.io-index)", - "extendr-api", + "clarabel 0.9.0 (registry+https://github.com/rust-lang/crates.io-index)", "lazy_static", "regex", + "savvy", ] [[package]] name = "clarabel" -version = "0.5.1" +version = "0.9.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "68d8445d4b01aae626d2530d19cfba03ef9e4f062f3aa7a353069a58a616eaee" +checksum = "83e62eacd93b899251364a22bd4dca0f293f8175d05311e42bbf6dbb5edcc762" dependencies = [ "amd", "blas", @@ -83,10 +89,14 @@ dependencies = [ "cfg-if", "derive_builder", "enum_dispatch", + "indexmap", + "itertools", "lapack", "lapack-src", "lazy_static", "num-traits", + "serde", + "serde_json", "thiserror", ] @@ -156,62 +166,72 @@ dependencies = [ "syn 1.0.109", ] +[[package]] +name = "either" +version = "1.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3dca9240753cf90908d7e4aac30f630662b02aebaa1b58a3cadabdb23385b58b" + [[package]] name = "enum_dispatch" -version = "0.3.11" +version = "0.3.13" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "11f36e95862220b211a6e2aa5eca09b4fa391b13cd52ceb8035a24bf65a79de2" +checksum = "aa18ce2bc66555b3218614519ac839ddb759a7d6720732f979ef8d13be147ecd" dependencies = [ "once_cell", "proc-macro2", "quote", - "syn 1.0.109", + "syn 2.0.66", ] [[package]] -name = "extendr-api" -version = "0.4.0" +name = "equivalent" +version = "1.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1e36d66fa307948c291a6fc5b09d8295dd58e88ab5e8d782d30e23670113e9ab" -dependencies = [ - "extendr-engine", - "extendr-macros", - "lazy_static", - "libR-sys", - "paste", -] +checksum = "5443807d6dff69373d433ab9ef5378ad8df50ca6298caf15de6e52e24aaf54d5" [[package]] -name = "extendr-engine" -version = "0.4.0" +name = "fnv" +version = "1.0.7" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8298d5a2e38bb91820b92bbd7e5aaf1d3b95ed9f096fc66393c50af38ff8155d" -dependencies = [ - "libR-sys", -] +checksum = "3f9eec918d3f24069decb9af1554cad7c880e2da24a9afd88aca000531ab82c1" [[package]] -name = "extendr-macros" -version = "0.4.0" +name = "hashbrown" +version = "0.14.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "09bf0849f0d48209be8163378248137fed5ccb5f464d171cf93a19f31a9e6c67" +checksum = "e5274423e17b7c9fc20b6e7e208532f9b19825d82dfd615708b70edd83df41f1" + +[[package]] +name = "ident_case" +version = "1.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b9e0384b61958566e926dc50660321d12159025e767c18e043daf26b70104c39" + +[[package]] +name = "indexmap" +version = "2.2.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "168fb715dda47215e360912c096649d23d58bf392ac62f73919e831745e40f26" dependencies = [ - "proc-macro2", - "quote", - "syn 1.0.109", + "equivalent", + "hashbrown", ] [[package]] -name = "fnv" -version = "1.0.7" +name = "itertools" +version = "0.11.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "3f9eec918d3f24069decb9af1554cad7c880e2da24a9afd88aca000531ab82c1" +checksum = "b1c173a5686ce8bfa551b3563d0c2170bf24ca44da99c7ca4bfdab5418c3fe57" +dependencies = [ + "either", +] [[package]] -name = "ident_case" -version = "1.0.1" +name = "itoa" +version = "1.0.11" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b9e0384b61958566e926dc50660321d12159025e767c18e043daf26b70104c39" +checksum = "49f1f14873335454500d59611f1cf4a4b0f786f9ac11f4312a78e4cf2566695b" [[package]] name = "lapack" @@ -226,9 +246,9 @@ dependencies = [ [[package]] name = "lapack-src" -version = "0.9.0" +version = "0.10.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "24c81fcc728418323178fd40407619d0ed26dbbbd1a553693c6290ef5d6698c6" +checksum = "a960a9d20df9abedf2e136c087c8fcc29f693ed985349bb03fe00816de8b00d2" dependencies = [ "r-src", ] @@ -248,68 +268,56 @@ version = "1.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" -[[package]] -name = "libR-sys" -version = "0.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cd728a97b9b0975f546bc865a7413e0ce6f98a8f6cea52e77dc5ee0bcea00adf" - [[package]] name = "libc" -version = "0.2.146" +version = "0.2.155" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f92be4933c13fd498862a9e02a3055f8a8d9c039ce33db97306fd5a6caa7f29b" +checksum = "97b3888a4aecf77e811145cadf6eef5901f4782c53886191b2f693f24761847c" [[package]] name = "memchr" -version = "2.5.0" +version = "2.7.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2dffe52ecf27772e601905b7522cb4ef790d2cc203488bbd0e2fe85fcb74566d" +checksum = "78ca9ab1a0babb1e7d5695e3530886289c18cf2f87ec19a575a0abdce112e3a3" [[package]] name = "num-complex" -version = "0.4.3" +version = "0.4.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "02e0d21255c828d6f128a1e41534206671e8c3ea0c62f32291e808dc82cff17d" +checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" dependencies = [ "num-traits", ] [[package]] name = "num-traits" -version = "0.2.15" +version = "0.2.19" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "578ede34cf02f8924ab9447f50c28075b4d3e5b269972345e7e0372b38c6cdcd" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" dependencies = [ "autocfg", ] [[package]] name = "once_cell" -version = "1.18.0" +version = "1.19.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dd8b5dd2ae5ed71462c540258bedcb51965123ad7e7ccf4b9a8cafaa4a63576d" - -[[package]] -name = "paste" -version = "1.0.12" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9f746c4065a8fa3fe23974dd82f15431cc8d40779821001404d10d2e79ca7d79" +checksum = "3fdb12b2476b595f9358c5161aa467c2438859caa136dec86c26fdd2efe17b92" [[package]] name = "proc-macro2" -version = "1.0.60" +version = "1.0.85" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dec2b086b7a862cf4de201096214fa870344cf922b2b30c167badb3af3195406" +checksum = "22244ce15aa966053a896d1accb3a6e68469b97c7f33f284b99f0d576879fc23" dependencies = [ "unicode-ident", ] [[package]] name = "quote" -version = "1.0.28" +version = "1.0.36" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1b9ab9c7eadfd8df19006f1cf1a4aed13540ed5cbc047010ece5826e10825488" +checksum = "0fa76aaf39101c457836aec0ce2316dbdc3ab723cdda1c6bd4e6ad4208acaca7" dependencies = [ "proc-macro2", ] @@ -322,9 +330,21 @@ checksum = "ea397956e1043a8d947ea84b13971d9cb30fce65ca66a921081755ff2e899b6a" [[package]] name = "regex" -version = "1.8.4" +version = "1.10.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d0ab3ca65655bb1e41f2a8c8cd662eb4fb035e67c3f78da1d61dffe89d07300f" +checksum = "b91213439dad192326a0d7c6ee3955910425f441d7038e0d6933b0aec5c4517f" +dependencies = [ + "aho-corasick", + "memchr", + "regex-automata", + "regex-syntax", +] + +[[package]] +name = "regex-automata" +version = "0.4.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "38caf58cc5ef2fed281f89292ef23f6365465ed9a41b7a7754eb4e26496c92df" dependencies = [ "aho-corasick", "memchr", @@ -333,9 +353,87 @@ dependencies = [ [[package]] name = "regex-syntax" -version = "0.7.2" +version = "0.8.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7a66a03ae7c801facd77a29370b4faec201768915ac14a721ba36f20bc9c209b" + +[[package]] +name = "ryu" +version = "1.0.18" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f3cb5ba0dc43242ce17de99c180e96db90b235b8a9fdc9543c96d2209116bd9f" + +[[package]] +name = "savvy" +version = "0.6.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a0ee0755ee6b80da65ee795ff2bd3484d53bd53f6a833433cea6f3248b7e501d" +dependencies = [ + "cc", + "once_cell", + "savvy-ffi", + "savvy-macro", +] + +[[package]] +name = "savvy-bindgen" +version = "0.6.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "958ef28111b7e18f3b056345d9cdbde9565c65f4358abdc95bf45cebe104e772" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.66", +] + +[[package]] +name = "savvy-ffi" +version = "0.6.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "436b050e76ed2903236f032a59761c1eb99e1b0aead2c257922771dab1fc8c78" +checksum = "fcd85a5d36724531674b53f6ac0ef03b250992d78758cf7ef7850566f4e5e7f4" + +[[package]] +name = "savvy-macro" +version = "0.6.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e922d5f96f21e6253d3566d09683edb519a161f9916d1c992835b68ea55b7fad" +dependencies = [ + "proc-macro2", + "quote", + "savvy-bindgen", + "syn 2.0.66", +] + +[[package]] +name = "serde" +version = "1.0.203" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7253ab4de971e72fb7be983802300c30b5a7f0c2e56fab8abfc6a214307c0094" +dependencies = [ + "serde_derive", +] + +[[package]] +name = "serde_derive" +version = "1.0.203" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "500cbc0ebeb6f46627f50f3f5811ccf6bf00643be300b4c3eabc0ef55dc5b5ba" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.66", +] + +[[package]] +name = "serde_json" +version = "1.0.117" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "455182ea6142b14f93f4bc5320a2b31c1f266b66a4a5c858b013302a5d8cbfc3" +dependencies = [ + "itoa", + "ryu", + "serde", +] [[package]] name = "strsim" @@ -356,9 +454,9 @@ dependencies = [ [[package]] name = "syn" -version = "2.0.18" +version = "2.0.66" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "32d41677bcbe24c20c52e7c70b0d8db04134c5d1066bf98662e2871ad200ea3e" +checksum = "c42f3f41a2de00b01c0aaad383c5a45241efc8b2d1eda5661812fda5f3cdcff5" dependencies = [ "proc-macro2", "quote", @@ -367,26 +465,26 @@ dependencies = [ [[package]] name = "thiserror" -version = "1.0.40" +version = "1.0.61" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "978c9a314bd8dc99be594bc3c175faaa9794be04a5a5e153caba6915336cebac" +checksum = "c546c80d6be4bc6a00c0f01730c08df82eaa7a7a61f11d656526506112cc1709" dependencies = [ "thiserror-impl", ] [[package]] name = "thiserror-impl" -version = "1.0.40" +version = "1.0.61" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f9456a42c5b0d803c8cd86e73dd7cc9edd429499f37a3550d286d5e86720569f" +checksum = "46c3384250002a6d5af4d114f2845d37b57521033f30d5c3f46c4d70e1197533" dependencies = [ "proc-macro2", "quote", - "syn 2.0.18", + "syn 2.0.66", ] [[package]] name = "unicode-ident" -version = "1.0.9" +version = "1.0.12" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b15811caf2415fb889178633e7724bad2509101cde276048e013b9def5e51fa0" +checksum = "3354b9ac3fae1ff6755cb6db53683adb661634f67557942dea4facebec0fee4b" diff --git a/src/rust/Cargo.toml b/src/rust/Cargo.toml index 8c2b9054..4bd34db9 100644 --- a/src/rust/Cargo.toml +++ b/src/rust/Cargo.toml @@ -1,14 +1,25 @@ [package] -name = 'clarabel' -version = '0.5.1' -edition = '2021' +name = "clarabel" +version = "0.9.0" +edition = "2021" [lib] -crate-type = [ 'staticlib' ] -name = 'clarabel' +crate-type = ["staticlib", "lib"] [dependencies] -extendr-api = '0.4.0' -clarabel = {version = '0.5.1', features = [ 'sdp-r' ]} +savvy = "0.6.4" +clarabel = {version = '0.9.0', features = [ 'sdp-r' ]} lazy_static = "1.4.0" -regex = "1.8.4" +regex = "1.10.4" + +[profile.release] +# By default, on release build, savvy terminates the R session when a panic +# occurs. This is the right behavior in that a panic means such a fatal event +# where we can have no hope of recovery. +# +# cf. https://doc.rust-lang.org/book/ch09-03-to-panic-or-not-to-panic.html +# +# However, it's possible that the panic is thrown by some of the dependency +# crate and there's little you can do. In such cases, you can change the +# following line to `panic = "unwind"` to always catch a panic. +panic = "abort" diff --git a/src/rust/api.h b/src/rust/api.h new file mode 100644 index 00000000..58ad2989 --- /dev/null +++ b/src/rust/api.h @@ -0,0 +1 @@ +SEXP savvy_clarabel_solve__ffi(SEXP m, SEXP n, SEXP Ai, SEXP Ap, SEXP Ax, SEXP b, SEXP q, SEXP Pi, SEXP Pp, SEXP Px, SEXP cone_spec, SEXP r_settings); diff --git a/src/rust/src/lib.rs b/src/rust/src/lib.rs index 66db4c98..2a2025d8 100644 --- a/src/rust/src/lib.rs +++ b/src/rust/src/lib.rs @@ -1,31 +1,25 @@ -use extendr_api::prelude::*; +#![allow(non_snake_case)] + +// Example functions +use savvy::savvy; +use savvy::{IntegerSexp, OwnedIntegerSexp, RealSexp, OwnedRealSexp, ListSexp, OwnedListSexp, TypedSexp}; use clarabel::algebra::*; use clarabel::solver::*; use lazy_static::lazy_static; use regex::Regex; + // Solve using Clarabel Rust solver. -#[extendr] -fn clarabel_solve(m: i32, n: i32, Ai: &[i32], Ap: &[i32], Ax: &[f64], b: &[f64], q: &[f64], - Pi: &[i32], Pp: &[i32], Px: &[f64], cone_spec: List, r_settings: List) -> Robj { - #![allow(non_snake_case)] +#[savvy] +fn clarabel_solve(m: i32, n: i32, Ai: IntegerSexp, Ap: IntegerSexp, Ax: RealSexp, b: RealSexp, q: RealSexp, Pi: IntegerSexp, Pp: IntegerSexp, Px: RealSexp, cone_spec: ListSexp, r_settings: ListSexp) -> savvy::Result { // QP Example let P = if Px.len() > 0 { - let mut SPp = Vec::::with_capacity(Pp.len()); - for i in 0..Pp.len() { - SPp.push(Pp[i] as usize); - } - let mut SPi = Vec::::with_capacity(Pi.len()); - for i in 0..Pi.len() { - SPi.push(Pi[i] as usize); - } - CscMatrix::new( n as usize, // m n as usize, // n - SPp, // colptr - SPi, // rowval - Px.to_vec(), // nzval + Pp.as_slice().iter().map(|&x| x as usize).collect(), // colptr + Pi.as_slice().iter().map(|&x| x as usize).collect(), // rowval + Px.as_slice().to_vec(), // nzval ) } else { CscMatrix::zeros((n as usize, n as usize)) // For P = 0 @@ -36,29 +30,22 @@ fn clarabel_solve(m: i32, n: i32, Ai: &[i32], Ap: &[i32], Ax: &[f64], b: &[f64], // assert!(P.check_format().is_ok()); let A = { - let mut SAp = Vec::::with_capacity(Ap.len()); - for i in 0..Ap.len() { - SAp.push(Ap[i] as usize); - } - let mut SAi = Vec::::with_capacity(Ai.len()); - for i in 0..Ai.len() { - SAi.push(Ai[i] as usize); - } - CscMatrix::new( m as usize, // m n as usize, // n - SAp, // colptr - SAi, // rowval - Ax.to_vec(), // nzval + Ap.as_slice().iter().map(|&x| x as usize).collect(), // colptr + Ai.as_slice().iter().map(|&x| x as usize).collect(), // rowval + Ax.as_slice().to_vec(), // nzval ) }; // println!("A: {:?}", A); // assert!(A.check_format().is_ok()); - // println!("b: {:?}", b); - // println!("q: {:?}", q); + let bb = b.as_slice().to_vec(); + let qq = q.as_slice().to_vec(); + // println!("b: {:?}", bb); + // println!("q: {:?}", qq); // Handle cones lazy_static! { @@ -68,101 +55,364 @@ fn clarabel_solve(m: i32, n: i32, Ai: &[i32], Ap: &[i32], Ax: &[f64], b: &[f64], static ref EPC: Regex = Regex::new("^ep").unwrap(); // Exponential Cone static ref PC: Regex = Regex::new("^p").unwrap(); // Power Cone static ref PSDTC: Regex = Regex::new("^s").unwrap(); // PSD Triangle Cone + static ref GPC: Regex = Regex::new("^gp").unwrap(); // Generalized Power Cone } let mut cones: Vec::> = Vec::new(); for (key, value) in cone_spec.iter() { + let typed_value = value.into_typed(); if ZC.is_match(key.as_ref()) { - cones.push(ZeroConeT(value.as_integer().expect("Positive integer expected") as usize)); + match typed_value { + TypedSexp::Integer(i) => cones.push(ZeroConeT(i.as_slice()[0] as usize)), + _ => (), + } } else if NNC.is_match(key.as_ref()) { - cones.push(NonnegativeConeT(value.as_integer().expect("Positive integer expected") as usize)); + match typed_value { + TypedSexp::Integer(i) => cones.push(NonnegativeConeT(i.as_slice()[0] as usize)), + _ => (), + } } else if SOC.is_match(key.as_ref()) { - cones.push(SecondOrderConeT(value.as_integer().expect("Positive integer expected") as usize)); + match typed_value { + TypedSexp::Integer(i) => cones.push(SecondOrderConeT(i.as_slice()[0] as usize)), + _ => (), + } } else if EPC.is_match(key.as_ref()) { - for _i in 0..value.as_integer().expect("Positive integer expected") { - cones.push(ExponentialConeT()); + match typed_value { + TypedSexp::Integer(i) => for _i in 0..(i.as_slice()[0] as usize) { + cones.push(ExponentialConeT()); + } + _ => (), } } else if PSDTC.is_match(key.as_ref()) { - cones.push(PSDTriangleConeT(value.as_integer().expect("Positive integer expected") as usize)); + match typed_value { + TypedSexp::Integer(i) => cones.push(PSDTriangleConeT(i.as_slice()[0] as usize)), + _ => (), + } } else if PC.is_match(key.as_ref()) { - cones.push(PowerConeT(value.as_real().expect("Positive real value expected"))); + match typed_value { + TypedSexp::Real(f) => cones.push(PowerConeT(f.as_slice()[0] as f64)), + _ => (), + } + } else if GPC.is_match(key.as_ref()) { + match typed_value { + TypedSexp::List(l) => { + match l.get("a").expect("GenPowerCone exponents vector").into_typed() { + TypedSexp::Real(f) => { + match l.get("n").expect("GenPowerCone dimension").into_typed() { + TypedSexp::Integer(i) => { + cones.push(GenPowerConeT(f.as_slice().to_vec(), i.as_slice()[0] as usize)) + }, + _ => (), + } + + }, + _ => (), + } + }, + _ => (), + } } else { - println!("Unknown cone {}; ignoring", &key); + let msg = format!("Ignoring unknown cone: {}", key); + let _ = savvy::io::r_warn(&msg); } } - // println!("cones: {:?}", cones); + println!("cones: {:?}", cones); // Update default settings with specified R settings for use below let settings = update_settings(r_settings); - let mut solver = DefaultSolver::new(&P, &q, &A, &b, &cones, settings); + let mut solver = DefaultSolver::new(&P, &qq, &A, &bb, &cones, settings); solver.solve(); - r!(list!(x = solver.solution.x.iter().collect_robj(), - z = solver.solution.z.iter().collect_robj(), - s = solver.solution.s.iter().collect_robj(), - obj_val = solver.solution.obj_val, - status = solver.solution.status as i32 + 1, // R's 1-based index - solve_time = solver.solution.solve_time, - iterations = solver.solution.iterations, - r_prim = solver.solution.r_prim, - r_dual = solver.solution.r_dual)) + + let mut obj_val = OwnedRealSexp::new(1)?; + obj_val[0] = solver.solution.obj_val; + + let mut obj_val_dual = OwnedRealSexp::new(1)?; + obj_val_dual[0] = solver.solution.obj_val_dual; + + let mut status = OwnedIntegerSexp::new(1)?; + status[0] = solver.solution.status as i32 + 1; // R's one-based index + + let mut solve_time = OwnedRealSexp::new(1)?; + solve_time[0] = solver.solution.solve_time; + + let mut iterations = OwnedIntegerSexp::new(1)?; + iterations[0] = solver.solution.iterations as i32; + + let mut r_prim = OwnedRealSexp::new(1)?; + r_prim[0] = solver.solution.r_prim; + + let mut r_dual = OwnedRealSexp::new(1)?; + r_dual[0] = solver.solution.r_dual; + + let mut out = OwnedListSexp::new(11, true)?; + out.set_name_and_value(0, "x", OwnedRealSexp::try_from_slice(solver.solution.x)?)?; + out.set_name_and_value(1, "z", OwnedRealSexp::try_from_slice(solver.solution.z)?)?; + out.set_name_and_value(2, "s", OwnedRealSexp::try_from_slice(solver.solution.s)?)?; + out.set_name_and_value(3, "obj_val", obj_val)?; + out.set_name_and_value(4, "obj_val_dual", obj_val_dual)?; + out.set_name_and_value(5, "status", status)?; + out.set_name_and_value(6, "solve_time", solve_time)?; + out.set_name_and_value(7, "iterations", iterations)?; + out.set_name_and_value(8, "r_prim", r_prim)?; + out.set_name_and_value(9, "r_dual", r_dual)?; + + // Handle Info + let fields = [ + ("μ", solver.info.μ), + ("sigma", solver.info.sigma), + ("step_length", solver.info.step_length), + ("cost_primal", solver.info.cost_primal), + ("cost_dual", solver.info.cost_dual), + ("res_primal", solver.info.res_primal), + ("res_dual", solver.info.res_dual), + ("res_primal_inf", solver.info.res_primal_inf), + ("res_dual_inf", solver.info.res_dual_inf), + ("gap_abs", solver.info.gap_abs), + ("gap_rel", solver.info.gap_rel), + ("ktratio", solver.info.ktratio), + ("solve_time", solver.info.solve_time), + ]; + + let mut info = OwnedListSexp::new(15, true)?; + for (i, (name, value)) in fields.iter().enumerate() { + let mut tmp = OwnedRealSexp::new(1)?; + tmp[0] = *value; + info.set_name_and_value(i, name, tmp)?; + } + let mut tmp = OwnedIntegerSexp::new(1)?; + tmp[0] = solver.info.iterations as i32; + info.set_name_and_value(13, "iterations", tmp)?; + let mut tmp = OwnedIntegerSexp::new(1)?; + tmp[0] = solver.info.status as i32 + 1; + info.set_name_and_value(14, "status", tmp)?; + + out.set_name_and_value(10, "info", info)?; + out.into() } -fn update_settings(r_settings: List) -> DefaultSettings { + +fn update_settings(r_settings: ListSexp) -> DefaultSettings { let mut settings = DefaultSettings::default(); // Should really be using something like structmap (https://github.com/ex0dus-0x/structmap) but // unsuccessful so far, so a directly generated implementation for (key, value) in r_settings.iter() { + let typed_value = value.into_typed(); + // auto-generated, see: inst/misc/code-gen.R match key.as_ref() { - "max_iter" => settings.max_iter = value.as_integer().expect("max_iter should be scalar integer") as u32, - "time_limit" => settings.time_limit = value.as_real().expect("should be scalar double"), - "verbose" => settings.verbose = value.as_bool().expect("verbose should be scalar boolean"), - "max_step_fraction" => settings.max_step_fraction = value.as_real().expect("should be scalar double"), - "tol_gap_abs" => settings.tol_gap_abs = value.as_real().expect("should be scalar double"), - "tol_gap_rel" => settings.tol_gap_rel = value.as_real().expect("should be scalar double"), - "tol_feas" => settings.tol_feas = value.as_real().expect("should be scalar double"), - "tol_infeas_abs" => settings.tol_infeas_abs = value.as_real().expect("should be scalar double"), - "tol_infeas_rel" => settings.tol_infeas_rel = value.as_real().expect("should be scalar double"), - "tol_ktratio" => settings.tol_ktratio = value.as_real().expect("should be scalar double"), - "reduced_tol_gap_abs" => settings.reduced_tol_gap_abs = value.as_real().expect("should be scalar double"), - "reduced_tol_gap_rel" => settings.reduced_tol_gap_rel = value.as_real().expect("should be scalar double"), - "reduced_tol_feas" => settings.reduced_tol_feas = value.as_real().expect("should be scalar double"), - "reduced_tol_infeas_abs" => settings.reduced_tol_infeas_abs = value.as_real().expect("should be scalar double"), - "reduced_tol_infeas_rel" => settings.reduced_tol_infeas_rel = value.as_real().expect("should be scalar double"), - "reduced_tol_ktratio" => settings.reduced_tol_ktratio = value.as_real().expect("should be scalar double"), - "equilibrate_enable" => settings.equilibrate_enable = value.as_bool().expect("equilibrate_enable should be scalar boolean"), - "equilibrate_max_iter" => settings.equilibrate_max_iter = value.as_integer().expect("Should be scalar integer") as u32, - "equilibrate_min_scaling" => settings.equilibrate_min_scaling = value.as_real().expect("should be scalar double"), - "equilibrate_max_scaling" => settings.equilibrate_max_scaling = value.as_real().expect("should be scalar double"), - "linesearch_backtrack_step" => settings.linesearch_backtrack_step = value.as_real().expect("should be scalar double"), - "min_switch_step_length" => settings.min_switch_step_length = value.as_real().expect("should be scalar double"), - "min_terminate_step_length" => settings.min_terminate_step_length = value.as_real().expect("should be scalar double"), - "direct_kkt_solver" => settings.direct_kkt_solver = value.as_bool().expect("direct_kkt_solver should be scalar boolean"), - "direct_solve_method" => settings.direct_solve_method = String::from(value.as_str().expect("direct_solve_method should be qdldl")), - "static_regularization_enable" => settings.static_regularization_enable = value.as_bool().expect("static_regularization_enable should be scalar boolean"), - "static_regularization_constant" => settings.static_regularization_constant = value.as_real().expect("should be scalar double"), - "static_regularization_proportional" => settings.static_regularization_proportional = value.as_real().expect("should be scalar double"), - "dynamic_regularization_enable" => settings.dynamic_regularization_enable = value.as_bool().expect("dynamic_regularization_enable should be scalar boolean"), - "dynamic_regularization_eps" => settings.dynamic_regularization_eps = value.as_real().expect("should be scalar double"), - "dynamic_regularization_delta" => settings.dynamic_regularization_delta = value.as_real().expect("should be scalar double"), - "iterative_refinement_enable" => settings.iterative_refinement_enable = value.as_bool().expect("iterative_refinement_enable should be scalar boolean"), - "iterative_refinement_reltol" => settings.iterative_refinement_reltol = value.as_real().expect("should be scalar double"), - "iterative_refinement_abstol" => settings.iterative_refinement_abstol = value.as_real().expect("should be scalar double"), - "iterative_refinement_max_iter" => settings.iterative_refinement_max_iter = value.as_integer().expect("Should be scalar integer") as u32, - "iterative_refinement_stop_ratio" => settings.iterative_refinement_stop_ratio = value.as_real().expect("should be scalar double"), - "presolve_enable" => settings.presolve_enable = value.as_bool().expect("presolve_enable should be scalar boolean"), - _ => (), - + "max_iter" => + match typed_value { + TypedSexp::Integer(i) => settings.max_iter = i.as_slice()[0] as u32, + _ => (), + }, + "time_limit" => + match typed_value { + TypedSexp::Real(f) => settings.time_limit = f.as_slice()[0], + _ => (), + }, + "verbose" => + match typed_value { + TypedSexp::Logical(b) => settings.verbose = b.as_slice_raw()[0] != 0, + _ => (), + }, + "max_step_fraction" => + match typed_value { + TypedSexp::Real(f) => settings.max_step_fraction = f.as_slice()[0], + _ => (), + }, + "tol_gap_abs" => + match typed_value { + TypedSexp::Real(f) => settings.tol_gap_abs = f.as_slice()[0], + _ => (), + }, + "tol_gap_rel" => + match typed_value { + TypedSexp::Real(f) => settings.tol_gap_rel = f.as_slice()[0], + _ => (), + }, + "tol_feas" => + match typed_value { + TypedSexp::Real(f) => settings.tol_feas = f.as_slice()[0], + _ => (), + }, + "tol_infeas_abs" => + match typed_value { + TypedSexp::Real(f) => settings.tol_infeas_abs = f.as_slice()[0], + _ => (), + }, + "tol_infeas_rel" => + match typed_value { + TypedSexp::Real(f) => settings.tol_infeas_rel = f.as_slice()[0], + _ => (), + }, + "tol_ktratio" => + match typed_value { + TypedSexp::Real(f) => settings.tol_ktratio = f.as_slice()[0], + _ => (), + }, + "reduced_tol_gap_abs" => + match typed_value { + TypedSexp::Real(f) => settings.reduced_tol_gap_abs = f.as_slice()[0], + _ => (), + }, + "reduced_tol_gap_rel" => + match typed_value { + TypedSexp::Real(f) => settings.reduced_tol_gap_rel = f.as_slice()[0], + _ => (), + }, + "reduced_tol_feas" => + match typed_value { + TypedSexp::Real(f) => settings.reduced_tol_feas = f.as_slice()[0], + _ => (), + }, + "reduced_tol_infeas_abs" => + match typed_value { + TypedSexp::Real(f) => settings.reduced_tol_infeas_abs = f.as_slice()[0], + _ => (), + }, + "reduced_tol_infeas_rel" => + match typed_value { + TypedSexp::Real(f) => settings.reduced_tol_infeas_rel = f.as_slice()[0], + _ => (), + }, + "reduced_tol_ktratio" => + match typed_value { + TypedSexp::Real(f) => settings.reduced_tol_ktratio = f.as_slice()[0], + _ => (), + }, + "equilibrate_enable" => + match typed_value { + TypedSexp::Logical(b) => settings.equilibrate_enable = b.as_slice_raw()[0] != 0, + _ => (), + }, + "equilibrate_max_iter" => + match typed_value { + TypedSexp::Integer(i) => settings.equilibrate_max_iter = i.as_slice()[0] as u32, + _ => (), + }, + "equilibrate_min_scaling" => + match typed_value { + TypedSexp::Real(f) => settings.equilibrate_min_scaling = f.as_slice()[0], + _ => (), + }, + "equilibrate_max_scaling" => + match typed_value { + TypedSexp::Real(f) => settings.equilibrate_max_scaling = f.as_slice()[0], + _ => (), + }, + "linesearch_backtrack_step" => + match typed_value { + TypedSexp::Real(f) => settings.linesearch_backtrack_step = f.as_slice()[0], + _ => (), + }, + "min_switch_step_length" => + match typed_value { + TypedSexp::Real(f) => settings.min_switch_step_length = f.as_slice()[0], + _ => (), + }, + "min_terminate_step_length" => + match typed_value { + TypedSexp::Real(f) => settings.min_terminate_step_length = f.as_slice()[0], + _ => (), + }, + "direct_kkt_solver" => + match typed_value { + TypedSexp::Logical(b) => settings.direct_kkt_solver = b.as_slice_raw()[0] != 0, + _ => (), + }, + "direct_solve_method" => + match typed_value { + TypedSexp::String(s) => if let Some(result) = s.to_vec().get(0) { + settings.direct_solve_method = result.to_string(); + }, + _ => (), + }, + "static_regularization_enable" => + match typed_value { + TypedSexp::Logical(b) => settings.static_regularization_enable = b.as_slice_raw()[0] != 0, + _ => (), + }, + "static_regularization_constant" => + match typed_value { + TypedSexp::Real(f) => settings.static_regularization_constant = f.as_slice()[0], + _ => (), + }, + "static_regularization_proportional" => + match typed_value { + TypedSexp::Real(f) => settings.static_regularization_proportional = f.as_slice()[0], + _ => (), + }, + "dynamic_regularization_enable" => + match typed_value { + TypedSexp::Logical(b) => settings.dynamic_regularization_enable = b.as_slice_raw()[0] != 0, + _ => (), + }, + "dynamic_regularization_eps" => + match typed_value { + TypedSexp::Real(f) => settings.dynamic_regularization_eps = f.as_slice()[0], + _ => (), + }, + "dynamic_regularization_delta" => + match typed_value { + TypedSexp::Real(f) => settings.dynamic_regularization_delta = f.as_slice()[0], + _ => (), + }, + "iterative_refinement_enable" => + match typed_value { + TypedSexp::Logical(b) => settings.iterative_refinement_enable = b.as_slice_raw()[0] != 0, + _ => (), + }, + "iterative_refinement_reltol" => + match typed_value { + TypedSexp::Real(f) => settings.iterative_refinement_reltol = f.as_slice()[0], + _ => (), + }, + "iterative_refinement_abstol" => + match typed_value { + TypedSexp::Real(f) => settings.iterative_refinement_abstol = f.as_slice()[0], + _ => (), + }, + "iterative_refinement_max_iter" => + match typed_value { + TypedSexp::Integer(i) => settings.iterative_refinement_max_iter = i.as_slice()[0] as u32, + _ => (), + }, + "iterative_refinement_stop_ratio" => + match typed_value { + TypedSexp::Real(f) => settings.iterative_refinement_stop_ratio = f.as_slice()[0], + _ => (), + }, + "presolve_enable" => + match typed_value { + TypedSexp::Logical(b) => settings.presolve_enable = b.as_slice_raw()[0] != 0, + _ => (), + }, + "chordal_decomposition_enable" => + match typed_value { + TypedSexp::Logical(b) => settings.chordal_decomposition_enable = b.as_slice_raw()[0] != 0, + _ => (), + }, + "chordal_decomposition_merge_method" => + match typed_value { + TypedSexp::String(s) => if let Some(result) = s.to_vec().get(0) { + settings.chordal_decomposition_merge_method = result.to_string(); + }, + _ => (), + }, + "chordal_decomposition_compact" => + match typed_value { + TypedSexp::Logical(b) => settings.chordal_decomposition_compact = b.as_slice_raw()[0] != 0, + _ => (), + }, + "chordal_decomposition_complete_dual" => + match typed_value { + TypedSexp::Logical(b) => settings.chordal_decomposition_complete_dual = b.as_slice_raw()[0] != 0, + _ => (), + }, + _ => (), } } settings } -// Macro to generate exports. -// This ensures exported functions are registered with R. -// See corresponding C code in `entrypoint.c`. -extendr_module! { - mod clarabel; - fn clarabel_solve; -} diff --git a/src/rust/vendor.sh b/src/rust/vendor.sh index 33c6d73c..7ea5d0a6 100644 --- a/src/rust/vendor.sh +++ b/src/rust/vendor.sh @@ -5,6 +5,7 @@ cargo vendor +## Clarabel has a CITATION.bib file that CRAN flags # c.f. https://reproducible-builds.org/docs/archives/ echo echo "Tarring up files" diff --git a/src/rust/vendor.tar.xz b/src/rust/vendor.tar.xz index 2349619e..37d5cc4e 100644 Binary files a/src/rust/vendor.tar.xz and b/src/rust/vendor.tar.xz differ diff --git a/tests/tinytest.R b/tests/tinytest.R index adf43dc0..619aa692 100644 --- a/tests/tinytest.R +++ b/tests/tinytest.R @@ -1,3 +1,7 @@ +l2_dist <- function(x, y) sqrt(sum((x-y)^2)) + +status_codes <- clarabel::solver_status_descriptions() + if ( requireNamespace("tinytest", quietly = TRUE) ){ tinytest::test_package("clarabel") } diff --git a/update_authors.R b/update_authors.R index 67e02d59..36af54b2 100644 --- a/update_authors.R +++ b/update_authors.R @@ -9,7 +9,8 @@ library(RcppTOML) VENDOR_PATH <- "src/rust/vendor" manifests <- list.files(VENDOR_PATH, pattern = "Cargo.toml", recursive = TRUE) - +manifests <- stringr::str_extract(manifests, "^[^/]+/Cargo\\.toml$") +manifests <- manifests[!is.na(manifests)] l <- lapply(manifests, \(x) RcppTOML::parseTOML(file.path(VENDOR_PATH, x))$package) names <- vapply(l, \(x) x[["name"]], FUN.VALUE = character(1L)) diff --git a/vignettes/clarabel.Rmd b/vignettes/clarabel.Rmd index 8f4b8b1e..69f981a6 100644 --- a/vignettes/clarabel.Rmd +++ b/vignettes/clarabel.Rmd @@ -21,8 +21,8 @@ library(clarabel) ## Introduction The first two examples are from the original [Clarabel -documentation](https://github.com/oxfordcontrol/ClarabelDocs) and the -third is from [SCS](https://www.cvxgrp.org/scs/). +documentation](https://clarabel.org) and the third is from +[SCS](https://www.cvxgrp.org/scs/). ## 1. Basic Quadratic Program Example @@ -349,13 +349,52 @@ A <- rbind( ) b <- c(vec(B1), vec(B2)) # stack both psd constraints cones <- list(s = c(2, 3)) # cone dimensions -s <- clarabel(A = A, b = b, q = q, cone = cones) +s <- clarabel(A = A, b = b, q = q, cones = cones) cat(sprintf("Solution status, description: = (%d, %s)\n", s$status, solver_status_descriptions()[s$status])) cat(sprintf("Solution (x1, x2, x3) = (%f, %f, %f)\n", s$x[1], s$x[2], s$x[3])) ``` -## 4. Control parameters +## 4. Cone Specifications + +The following cones can be specified in Clarabel. + +```{r, echo = FALSE} +parameter_df <- data.frame( + Parameter = c("z", "l", "q", "s", "ep", "p", "gp"), + Type = c("integer", "integer", "integer", "integer", "integer", "numeric", "list"), + Length = c("1", "1", ">= 1", ">= 1", "1", ">= 1", ">= 1"), + Description = c( + "Number of primal zero cones (dual free cones), which corresponds to the primal equality constraints", + "Number of linear cones (non-negative cones)", + "Vector of second-order cone sizes", + "Vector of positive semidefinite cone sizes", + "Number of primal exponential cones", + "Vector of primal power cone parameters", + "List of named lists of two items, `a` : the numeric vector of at least 2 exponent terms, and `n` : an integer dimension of generalized power cone parameters" + ), + Definition = c("$\\{ 0 \\}^{z}$", + "$\\{ x \\in \\mathbb{R}^{l} : x_i \\ge 0, \\forall i=1,\\dots,l \\}$", + "$\\{ (t,x) \\in \\mathbb{R}^{q} : \\lVert x\\rVert_2 \\leq t \\}$", + "Upper triangular part of the positive semidefinite cone $S^s_+$. The elements $x$ of this cone represent the columnwise stacking of the upper triangular part of a positive semidefinite matrix $X \\in S^s_+$, so that $x \\in R^d$ with $d = s(s+1)/2$", + "$\\{(x, y, z) : y > 0,~~ ye^{x/y} \\le z \\}$", + "$\\{(x, y, z) : x^p y^{(1-p)} \\ge \\lVert z\\rVert,~ (x,y) \\ge 0 \\}$ with $p \\in (0,1)$", + "$\\{(x, y) \\in R^{len(a)} \\times R^n : \\prod\\limits_{a_i \\in a} x_i^{a_i} \\ge \\lVert y\\rVert_2,~ x \\ge 0 \\}$ with $a_i \\in (0,1)$ and $\\sum a_i = 1$" + ) +) +names(parameter_df)[5] <- "Definition (per parameter element)" +knitr::kable(parameter_df) +``` +Generalized power cone parameters are specified as list of two-item +lists, with component named $a$ denoting the exponents and the named +component $n$ denoting the dimension. + +One can specify cones in any order if `strict_cone_order` is set to +`FALSE` in the call to `clarabel()` but one has to ensure that +parameter types are strictly specified for the values, e.g. `5L` for integers, `0.` +for reals etc. + +## 5. Control parameters Clarabel has a number of parameters that control its behavior, including verbosity, time limits, and tolerances; see help on