diff --git a/DESCRIPTION b/DESCRIPTION index 5af65f37..453cd81e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -16,7 +16,7 @@ Version: 1.2.0 Date: 2024-02-29 Maintainer: Vissarion Fisikopoulos Depends: Rcpp (>= 0.12.17) -Imports: methods, stats +Imports: methods, stats, Matrix LinkingTo: Rcpp, RcppEigen, BH Suggests: testthat Encoding: UTF-8 diff --git a/NAMESPACE b/NAMESPACE index 0d88effd..b1fba2f0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -19,8 +19,6 @@ export(gen_simplex) export(gen_skinny_cube) export(geweke) export(inner_ball) -export(loadSdpaFormatFile) -export(ode_solve) export(pinvweibull_with_loc) export(psrf_multivariate) export(psrf_univariate) @@ -30,10 +28,10 @@ export(rotate_polytope) export(round_polytope) export(sample_points) export(volume) -export(writeSdpaFormatFile) export(write_sdpa_format_file) export(zonotope_approximation) exportClasses(Hpolytope) +exportClasses(HpolytopeSparse) exportClasses(Spectrahedron) exportClasses(Vpolytope) exportClasses(VpolytopeIntersection) diff --git a/NEWS.md b/NEWS.md index cbc58db0..e3e98b6a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -54,7 +54,7 @@ # volesti 1.2.0 - New functions: dinvweibull_with_loc, ess, estimtate_lipschitz_constant, gen_birkhoff, geweke -ode_solve, pinvweibull_with_loc, psrf_multivariate, psrf_univariate, raftery +pinvweibull_with_loc, psrf_multivariate, psrf_univariate, raftery - New features in sample_points function: a) new walks: i) Dikin walk, ii) Vaidya walk, iii) John walk, diff --git a/R/HpolytopeSparseClass.R b/R/HpolytopeSparseClass.R index 394396f4..6bde225c 100644 --- a/R/HpolytopeSparseClass.R +++ b/R/HpolytopeSparseClass.R @@ -1,26 +1,41 @@ -#' An R class to represent an H-polytope +#' An R class to represent an H-polytope defined by a sparse matrix #' -#' An H-polytope is a convex polytope defined by a set of linear inequalities or equivalently a \eqn{d}-dimensional H-polytope with \eqn{m} facets is defined by a \eqn{m\times d} matrix A and a \eqn{m}-dimensional vector b, s.t.: \eqn{Ax\leq b}. +#' A sparse H-polytope is a convex polytope defined by a set of linear inequalities or equivalently a +#' \eqn{d}-dimensional H-polytope with \eqn{m} facets is defined by a \eqn{m\times d} matrix A and +#' a \eqn{m}-dimensional vector b, s.t.: \eqn{Ax\leq b}, where $A$ is sparse. #' #' \describe{ -#' \item{A}{An \eqn{m\times d} numerical matrix.} +#' \item{Aineq}{An \eqn{m\times d} sparse numerical matrix for the inequalities.} #' -#' \item{b}{An \eqn{m}-dimensional vector b.} +#' \item{bineq}{An \eqn{m}-dimensional vector for the ineqaulities.} #' -#' \item{volume}{The volume of the polytope if it is known, \eqn{NaN} otherwise by default.} +#' \item{Aineq}{An \eqn{m'\times d} sparse numerical matrix for the equalities.} #' -#' \item{type}{A character with default value 'Hpolytope', to declare the representation of the polytope.} +#' \item{bineq}{An \eqn{m'}-dimensional vector for the eqaulities.} +#' +#' \item{lb}{\eqn{d}-dimensional vector lb.} +#' +#' \item{ub}{\eqn{d}-dimensional vector ub.} +#' +#' \item{type}{A character with default value 'HpolytopeSparse', to declare the representation of the polytope.} #' } #' #' @examples -#' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) -#' b = c(0,0,1) -#' P = Hpolytope(A = A, b = b) +#' library(Matrix) +#' bineq=c(10,10,10,10,10) +#' Aineq = matrix(c(1,0,-0.25,-1,2.5,1,0.4,-1,-0.9,0.5), nrow=5, ncol=2, byrow = TRUE) +#' Aineq = as( Aineq, 'dgCMatrix' ) +#' beq=(0) +#' Aeq = matrix(, nrow=0, ncol=2, byrow = TRUE) +#' Aeq=as( Aeq, 'dgCMatrix' ) +#' lb=-100000*c(1,1); +#' ub=100000*c(1,1); +#' P <- HpolytopeSparse(Aineq = Aineq, bineq = bineq, Aeq = Aeq, beq = beq, lb = lb, ub = ub) #' #' @name HpolytopeSparse-class #' @rdname HpolytopeSparse-class #' @exportClass HpolytopeSparse -library(Matrix) +#' @importFrom "Matrix" "CsparseMatrix" HpolytopeSparse <- setClass ( # Class name "HpolytopeSparse", @@ -33,19 +48,17 @@ HpolytopeSparse <- setClass ( beq = "numeric", lb = "numeric", ub = "numeric", - dimension = "numeric", type = "character" ), # Initializing slots prototype = list( - Aineq = as(0, "CsparseMatrix"), - bineq = as.numeric(NULL), - Aeq = as(0, "CsparseMatrix"), - beq = as.numeric(NULL), - lb = as.numeric(NULL), - ub = as.numeric(NULL), - dimension = as.numeric(NaN), + Aineq = Matrix::sparseMatrix(i = integer(0), j = integer(0), x = double(0), dims = c(0, 0)), + bineq = numeric(0), + Aeq = Matrix::sparseMatrix(i = integer(0), j = integer(0), x = double(0), dims = c(0, 0)), + beq = numeric(0), + lb = numeric(0), + ub = numeric(0), type = "HpolytopeSparse" ) ) diff --git a/R/RcppExports.R b/R/RcppExports.R index 2ada6786..c4e93dd8 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -187,30 +187,6 @@ load_sdpa_format_file <- function(input_file = NULL) { .Call(`_volesti_load_sdpa_format_file`, input_file) } -#' Solve an ODE of the form dx^n / dt^n = F(x, t) -#' -#' @param n The number of steps. -#' @param step_size The step size. -#' @param order The ODE order (default is n = 1) -#' @param dimension The dimension of each derivative -#' @param initial_time The initial time -#' @param F The function oracle F(x, t) in the ODE. -#' @param method The method to be used -#' @param initial_conditions The initial conditions provided to the solver. Must be provided in a list with keys "x_1", ..., "x_n" and column vectors as values. The state "x_n" represents the (n-1)-th order derivative with respect to time -#' @param domains A list of n H-polytopes with keys "P_1", "P_2", ..., "P_n" that correspond to each derivative's domain -#' -#' @return A list which contains elements "x_1", ..., "x_n" representing each derivative results. Each "x_i" corresponds to a d x n matrix where each column represents a certain timestep of the solver. -#' -#' @examples -#' F <- function (x) (-x) -#' initial_conditions <- list("x_1" = c(0), "x_2" = c(1)) -#' states <- ode_solve(dimension=1, n=1000, F=F, initial_time=0, step_size=0.01, order=2, method="leapfrog", initial_conditions=initial_conditions, domains = list()) -#' -#' @export -ode_solve <- function(n, step_size, order, dimension, initial_time, F, method, domains = NULL, initial_conditions = NULL) { - .Call(`_volesti_ode_solve`, n, step_size, order, dimension, initial_time, F, method, domains, initial_conditions) -} - #' An internal Rccp function as a polytope generator #' #' @param kind_gen An integer to declare the type of the polytope. @@ -384,44 +360,6 @@ sample_points <- function(P, n, random_walk = NULL, distribution = NULL, seed = .Call(`_volesti_sample_points`, P, n, random_walk, distribution, seed) } -#' Write a SDPA format file -#' -#' Outputs a spectrahedron (the matrices defining a linear matrix inequality) and a vector (the objective function) -#' to a SDPA format file. -#' -#' @param spectrahedron A spectrahedron in n dimensions; must be an object of class Spectrahedron -#' @param objectiveFunction A numerical vector of length n -#' @param outputFile Name of the output file -#' -#' @examples -#' \dontrun{ -#' A0 = matrix(c(-1,0,0,0,-2,1,0,1,-2), nrow=3, ncol=3, byrow = TRUE) -#' A1 = matrix(c(-1,0,0,0,0,1,0,1,0), nrow=3, ncol=3, byrow = TRUE) -#' A2 = matrix(c(0,0,-1,0,0,0,-1,0,0), nrow=3, ncol=3, byrow = TRUE) -#' lmi = list(A0, A1, A2) -#' S = Spectrahedron(matrices = lmi) -#' objFunction = c(1,1) -#' writeSdpaFormatFile(S, objFunction, "output.txt") -#' } -#' @export -writeSdpaFormatFile <- function(spectrahedron = NULL, objectiveFunction = NULL, outputFile = NULL) { - invisible(.Call(`_volesti_writeSdpaFormatFile`, spectrahedron, objectiveFunction, outputFile)) -} - -#' Read a SDPA format file -#' -#' @param inputFile Name of the input file -#' -#' @return A list with two named items: an item "matrices" which is a list of the matrices and an vector "objFunction" -#' -#' @examples -#' path = system.file('extdata', package = 'volesti') -#' l = loadSdpaFormatFile(paste0(path,'/sdpa_n2m3.txt')) -#' @export -loadSdpaFormatFile <- function(inputFile = NULL) { - .Call(`_volesti_loadSdpaFormatFile`, inputFile) -} - #' The main function for volume approximation of a convex Polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes). It returns a list with two elements: (a) the logarithm of the estimated volume and (b) the estimated volume #' #' For the volume approximation can be used three algorithms. Either CoolingBodies (CB) or SequenceOfBalls (SOB) or CoolingGaussian (CG). An H-polytope with \eqn{m} facets is described by a \eqn{m\times d} matrix \eqn{A} and a \eqn{m}-dimensional vector \eqn{b}, s.t.: \eqn{P=\{x\ |\ Ax\leq b\} }. A V-polytope is defined as the convex hull of \eqn{m} \eqn{d}-dimensional points which correspond to the vertices of P. A zonotope is desrcibed by the Minkowski sum of \eqn{m} \eqn{d}-dimensional segments. diff --git a/man/HpolytopeSparse-class.Rd b/man/HpolytopeSparse-class.Rd new file mode 100644 index 00000000..c3e1f4bf --- /dev/null +++ b/man/HpolytopeSparse-class.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/HpolytopeSparseClass.R +\docType{class} +\name{HpolytopeSparse-class} +\alias{HpolytopeSparse-class} +\alias{HpolytopeSparse} +\title{An R class to represent an H-polytope defined by a sparse matrix} +\description{ +A sparse H-polytope is a convex polytope defined by a set of linear inequalities or equivalently a +\eqn{d}-dimensional H-polytope with \eqn{m} facets is defined by a \eqn{m\times d} matrix A and +a \eqn{m}-dimensional vector b, s.t.: \eqn{Ax\leq b}, where $A$ is sparse. +} +\details{ +\describe{ + \item{Aineq}{An \eqn{m\times d} sparse numerical matrix for the inequalities.} + + \item{bineq}{An \eqn{m}-dimensional vector for the ineqaulities.} + + \item{Aineq}{An \eqn{m'\times d} sparse numerical matrix for the equalities.} + + \item{bineq}{An \eqn{m'}-dimensional vector for the eqaulities.} + + \item{lb}{\eqn{d}-dimensional vector lb.} + + \item{ub}{\eqn{d}-dimensional vector ub.} + + \item{type}{A character with default value 'HpolytopeSparse', to declare the representation of the polytope.} +} +} +\examples{ +library(Matrix) +bineq=c(10,10,10,10,10) +Aineq = matrix(c(1,0,-0.25,-1,2.5,1,0.4,-1,-0.9,0.5), nrow=5, ncol=2, byrow = TRUE) +Aineq = as( Aineq, 'dgCMatrix' ) +beq=(0) +Aeq = matrix(, nrow=0, ncol=2, byrow = TRUE) +Aeq=as( Aeq, 'dgCMatrix' ) +lb=-100000*c(1,1); +ub=100000*c(1,1); +P <- HpolytopeSparse(Aineq = Aineq, bineq = bineq, Aeq = Aeq, beq = beq, lb = lb, ub = ub) + +} diff --git a/man/examples/simple_ode.R b/man/examples/simple_ode.R deleted file mode 100644 index 10d0c6d8..00000000 --- a/man/examples/simple_ode.R +++ /dev/null @@ -1,31 +0,0 @@ -# VolEsti (volume computation and sampling library) - -# Copyright (c) 2012-2020 Vissarion Fisikopoulos -# Copyright (c) 2018-2020 Apostolos Chalkis -# Copyright (c) 2020-2020 Marios Papachristou - -# Contributed and/or modified by Marios Papachristou, as part of Google Summer of Code 2020 program. - -# Licensed under GNU LGPL.3, see LICENCE file - -# Example script for ODE solvers -library(volesti) - -F <- function (x) (-x) -order <- 2 -step_size <- 0.01 -n <- 1000 -initial_conditions <- list("x_1" = c(0), "x_2" = c(1)) -initial_time <- 0 - -# Do not impose constraints -domains <- list() - -# Call the ode solver -states <- volesti::ode_solve(dimension=1, n=n, F=F, initial_time=initial_time, step_size=step_size, order=order, method="leapfrog", initial_conditions=initial_conditions, domains = list()) - -x <- states[["x_1"]] -v <- states[["x_2"]] - -plot(x, v) - diff --git a/man/examples/simple_ode_truncated.R b/man/examples/simple_ode_truncated.R deleted file mode 100644 index 17470a3c..00000000 --- a/man/examples/simple_ode_truncated.R +++ /dev/null @@ -1,34 +0,0 @@ -# VolEsti (volume computation and sampling library) - -# Copyright (c) 2012-2020 Vissarion Fisikopoulos -# Copyright (c) 2018-2020 Apostolos Chalkis -# Copyright (c) 2020-2020 Marios Papachristou - -# Contributed and/or modified by Marios Papachristou, as part of Google Summer of Code 2020 program. - -# Licensed under GNU LGPL.3, see LICENCE file - -# Example script for truncated ODE solvers -library(volesti) - -F <- function (x) (x) -order <- 1 -step_size <- 0.01 -n <- 1000 -initial_conditions <- list("x_1" = c(0.1)) -initial_time <- 0 - -A <- matrix(c(1, -1), ncol=1, nrow=2, byrow=TRUE) -b <- c(1, 0) - -# Create domain of truncation -P_1 <- volesti::Hpolytope(A = A, b = b) -domains <- list("P_1" = P_1) - -# Call the ode solver -states <- ode_solve(dimension=1, n=n, F=F, initial_time=initial_time, step_size=step_size, order=order, method="euler", initial_conditions=initial_conditions, domains = domains) - -x <- states[["x_1"]] -t <- step_size * seq(0, n - 1) - -plot(t, x) diff --git a/man/loadSdpaFormatFile.Rd b/man/loadSdpaFormatFile.Rd deleted file mode 100644 index bb6aa755..00000000 --- a/man/loadSdpaFormatFile.Rd +++ /dev/null @@ -1,21 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{loadSdpaFormatFile} -\alias{loadSdpaFormatFile} -\title{Read a SDPA format file} -\usage{ -loadSdpaFormatFile(inputFile = NULL) -} -\arguments{ -\item{inputFile}{Name of the input file} -} -\value{ -A list with two named items: an item "matrices" which is a list of the matrices and an vector "objFunction" -} -\description{ -Read a SDPA format file -} -\examples{ -path = system.file('extdata', package = 'volesti') -l = loadSdpaFormatFile(paste0(path,'/sdpa_n2m3.txt')) -} diff --git a/man/ode_solve.Rd b/man/ode_solve.Rd deleted file mode 100644 index 2836f342..00000000 --- a/man/ode_solve.Rd +++ /dev/null @@ -1,49 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{ode_solve} -\alias{ode_solve} -\title{Solve an ODE of the form dx^n / dt^n = F(x, t)} -\usage{ -ode_solve( - n, - step_size, - order, - dimension, - initial_time, - F, - method, - domains = NULL, - initial_conditions = NULL -) -} -\arguments{ -\item{n}{The number of steps.} - -\item{step_size}{The step size.} - -\item{order}{The ODE order (default is n = 1)} - -\item{dimension}{The dimension of each derivative} - -\item{initial_time}{The initial time} - -\item{F}{The function oracle F(x, t) in the ODE.} - -\item{method}{The method to be used} - -\item{domains}{A list of n H-polytopes with keys "P_1", "P_2", ..., "P_n" that correspond to each derivative's domain} - -\item{initial_conditions}{The initial conditions provided to the solver. Must be provided in a list with keys "x_1", ..., "x_n" and column vectors as values. The state "x_n" represents the (n-1)-th order derivative with respect to time} -} -\value{ -A list which contains elements "x_1", ..., "x_n" representing each derivative results. Each "x_i" corresponds to a d x n matrix where each column represents a certain timestep of the solver. -} -\description{ -Solve an ODE of the form dx^n / dt^n = F(x, t) -} -\examples{ -F <- function (x) (-x) -initial_conditions <- list("x_1" = c(0), "x_2" = c(1)) -states <- ode_solve(dimension=1, n=1000, F=F, initial_time=0, step_size=0.01, order=2, method="leapfrog", initial_conditions=initial_conditions, domains = list()) - -} diff --git a/man/sample_points.Rd b/man/sample_points.Rd index ed0cc961..ab633930 100644 --- a/man/sample_points.Rd +++ b/man/sample_points.Rd @@ -13,7 +13,17 @@ sample_points(P, n, random_walk = NULL, distribution = NULL, seed = NULL) \item{random_walk}{Optional. A list that declares the random walk and some related parameters as follows: \itemize{ -\item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, v) \code{'dikin'} for dikin walk, vi) \code{'vaidya'} for vaidya walk, vii) \code{'john'} for john walk, viii) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or ix) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR x) \code{'HMC'} for Hamiltonian Monte Carlo (logconcave densities) xi) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method xii) \code{'ExactHMC'} for exact Hamiltonian Monte Carlo with reflections (spherical Gaussian or exponential distribution). The default walk is \code{'aBiW'} for the uniform distribution or \code{'CDHR'} for the Gaussian distribution and H-polytopes and \code{'BiW'} or \code{'RDHR'} for the same distributions and V-polytopes and zonotopes.} +\item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run, +ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk, +v) \code{'dikin'} for dikin walk, vi) \code{'vaidya'} for vaidya walk, vii) \code{'john'} for john walk, +viii) \code{'BCDHR'} boundary sampling by keeping the extreme points of CDHR or ix) \code{'BRDHR'} boundary sampling by keeping the extreme points of RDHR, +x) \code{'NUTS'} for NUTS Hamiltonian Monte Carlo sampler (logconcave densities), xi) \code{'HMC'} for Hamiltonian Monte Carlo (logconcave densities), +xii) CRHMC for Riemannian HMC with H-polytope constraints (uniform and general logconcave densities), +xiii) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method (logconcave densities), +xiii) \code{'ExactHMC'} for exact Hamiltonian Monte Carlo with reflections (spherical Gaussian or exponential distribution). +The default walk is \code{'aBiW'} for the uniform distribution, \code{'CDHR'} for the Gaussian distribution and H-polytopes and +\code{'BiW'} or \code{'RDHR'} for the same distributions and V-polytopes and zonotopes. \code{'NUTS'} is the default sampler for logconcave densities and \code{'CRHMC'} +for logconcave densities with H-polytope and sparse constrainted problems.} \item{\code{walk_length} }{ The number of the steps per generated point for the random walk. The default value is \eqn{1}.} \item{\code{nburns} }{ The number of points to burn before start sampling. The default value is \eqn{1}.} \item{\code{starting_point} }{ A \eqn{d}-dimensional numerical vector that declares a starting point in the interior of the polytope for the random walk. The default choice is the center of the ball as that one computed by the function \code{inner_ball()}.} @@ -51,7 +61,7 @@ points = sample_points(P, n = 100, random_walk = list("walk" = "BaW", "walk_leng # gaussian distribution from the 2d unit simplex in H-representation with variance = 2 A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE) b = c(0,0,1) -P = Hpolytope(A = A, b = b) +P = Hpolytope(A=A,b=b) points = sample_points(P, n = 100, distribution = list("density" = "gaussian", "variance" = 2)) # uniform points from the boundary of a 2-dimensional random H-polytope diff --git a/man/writeSdpaFormatFile.Rd b/man/writeSdpaFormatFile.Rd deleted file mode 100644 index f70cac91..00000000 --- a/man/writeSdpaFormatFile.Rd +++ /dev/null @@ -1,34 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/RcppExports.R -\name{writeSdpaFormatFile} -\alias{writeSdpaFormatFile} -\title{Write a SDPA format file} -\usage{ -writeSdpaFormatFile( - spectrahedron = NULL, - objectiveFunction = NULL, - outputFile = NULL -) -} -\arguments{ -\item{spectrahedron}{A spectrahedron in n dimensions; must be an object of class Spectrahedron} - -\item{objectiveFunction}{A numerical vector of length n} - -\item{outputFile}{Name of the output file} -} -\description{ -Outputs a spectrahedron (the matrices defining a linear matrix inequality) and a vector (the objective function) -to a SDPA format file. -} -\examples{ -\dontrun{ -A0 = matrix(c(-1,0,0,0,-2,1,0,1,-2), nrow=3, ncol=3, byrow = TRUE) -A1 = matrix(c(-1,0,0,0,0,1,0,1,0), nrow=3, ncol=3, byrow = TRUE) -A2 = matrix(c(0,0,-1,0,0,0,-1,0,0), nrow=3, ncol=3, byrow = TRUE) -lmi = list(A0, A1, A2) -S = Spectrahedron(matrices = lmi) -objFunction = c(1,1) -writeSdpaFormatFile(S, objFunction, "output.txt") -} -} diff --git a/src/Makevars b/src/Makevars index 173e17f8..eb27e4d5 100644 --- a/src/Makevars +++ b/src/Makevars @@ -1,5 +1,5 @@ PKG_CPPFLAGS=-Iexternal -Iexternal/lpSolve/src -Iexternal/minimum_ellipsoid -Ivolesti/include -Ivolesti/include/convex_bodies/spectrahedra -PKG_CXXFLAGS= -DBOOST_NO_AUTO_PTR -DDISABLE_NLP_ORACLES -Wno-deprecated-declarations -lm -ldl -Wno-ignored-attributes +PKG_CXXFLAGS= -DBOOST_NO_AUTO_PTR -DDISABLE_NLP_ORACLES PKG_LIBS=-Lexternal/lpSolve/src -llp_solve -Lexternal/PackedCSparse/qd -lqd $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 93e873c4..594a5cbb 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -110,25 +110,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// ode_solve -Rcpp::List ode_solve(Rcpp::Nullable n, Rcpp::Nullable step_size, Rcpp::Nullable order, Rcpp::Nullable dimension, Rcpp::Nullable initial_time, Rcpp::Function F, Rcpp::Nullable method, Rcpp::Nullable domains, Rcpp::Nullable initial_conditions); -RcppExport SEXP _volesti_ode_solve(SEXP nSEXP, SEXP step_sizeSEXP, SEXP orderSEXP, SEXP dimensionSEXP, SEXP initial_timeSEXP, SEXP FSEXP, SEXP methodSEXP, SEXP domainsSEXP, SEXP initial_conditionsSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::Nullable >::type n(nSEXP); - Rcpp::traits::input_parameter< Rcpp::Nullable >::type step_size(step_sizeSEXP); - Rcpp::traits::input_parameter< Rcpp::Nullable >::type order(orderSEXP); - Rcpp::traits::input_parameter< Rcpp::Nullable >::type dimension(dimensionSEXP); - Rcpp::traits::input_parameter< Rcpp::Nullable >::type initial_time(initial_timeSEXP); - Rcpp::traits::input_parameter< Rcpp::Function >::type F(FSEXP); - Rcpp::traits::input_parameter< Rcpp::Nullable >::type method(methodSEXP); - Rcpp::traits::input_parameter< Rcpp::Nullable >::type domains(domainsSEXP); - Rcpp::traits::input_parameter< Rcpp::Nullable >::type initial_conditions(initial_conditionsSEXP); - rcpp_result_gen = Rcpp::wrap(ode_solve(n, step_size, order, dimension, initial_time, F, method, domains, initial_conditions)); - return rcpp_result_gen; -END_RCPP -} // poly_gen Rcpp::NumericMatrix poly_gen(int kind_gen, bool Vpoly_gen, bool Zono_gen, int dim_gen, int m_gen, Rcpp::Nullable seed); RcppExport SEXP _volesti_poly_gen(SEXP kind_genSEXP, SEXP Vpoly_genSEXP, SEXP Zono_genSEXP, SEXP dim_genSEXP, SEXP m_genSEXP, SEXP seedSEXP) { @@ -223,29 +204,6 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// writeSdpaFormatFile -void writeSdpaFormatFile(Rcpp::Reference spectrahedron, Rcpp::Nullable objectiveFunction, Rcpp::Nullable outputFile); -RcppExport SEXP _volesti_writeSdpaFormatFile(SEXP spectrahedronSEXP, SEXP objectiveFunctionSEXP, SEXP outputFileSEXP) { -BEGIN_RCPP - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::Reference >::type spectrahedron(spectrahedronSEXP); - Rcpp::traits::input_parameter< Rcpp::Nullable >::type objectiveFunction(objectiveFunctionSEXP); - Rcpp::traits::input_parameter< Rcpp::Nullable >::type outputFile(outputFileSEXP); - writeSdpaFormatFile(spectrahedron, objectiveFunction, outputFile); - return R_NilValue; -END_RCPP -} -// loadSdpaFormatFile -Rcpp::List loadSdpaFormatFile(Rcpp::Nullable inputFile); -RcppExport SEXP _volesti_loadSdpaFormatFile(SEXP inputFileSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::Nullable >::type inputFile(inputFileSEXP); - rcpp_result_gen = Rcpp::wrap(loadSdpaFormatFile(inputFile)); - return rcpp_result_gen; -END_RCPP -} // volume Rcpp::List volume(Rcpp::Reference P, Rcpp::Nullable settings, Rcpp::Nullable rounding); RcppExport SEXP _volesti_volume(SEXP PSEXP, SEXP settingsSEXP, SEXP roundingSEXP) { @@ -295,7 +253,6 @@ static const R_CallMethodDef CallEntries[] = { {"_volesti_geweke", (DL_FUNC) &_volesti_geweke, 3}, {"_volesti_inner_ball", (DL_FUNC) &_volesti_inner_ball, 2}, {"_volesti_load_sdpa_format_file", (DL_FUNC) &_volesti_load_sdpa_format_file, 1}, - {"_volesti_ode_solve", (DL_FUNC) &_volesti_ode_solve, 9}, {"_volesti_poly_gen", (DL_FUNC) &_volesti_poly_gen, 6}, {"_volesti_psrf_multivariate", (DL_FUNC) &_volesti_psrf_multivariate, 1}, {"_volesti_psrf_univariate", (DL_FUNC) &_volesti_psrf_univariate, 2}, @@ -303,8 +260,6 @@ static const R_CallMethodDef CallEntries[] = { {"_volesti_rotating", (DL_FUNC) &_volesti_rotating, 3}, {"_volesti_rounding", (DL_FUNC) &_volesti_rounding, 3}, {"_volesti_sample_points", (DL_FUNC) &_volesti_sample_points, 5}, - {"_volesti_writeSdpaFormatFile", (DL_FUNC) &_volesti_writeSdpaFormatFile, 3}, - {"_volesti_loadSdpaFormatFile", (DL_FUNC) &_volesti_loadSdpaFormatFile, 1}, {"_volesti_volume", (DL_FUNC) &_volesti_volume, 3}, {"_volesti_write_sdpa_format_file", (DL_FUNC) &_volesti_write_sdpa_format_file, 3}, {"_volesti_zono_approx", (DL_FUNC) &_volesti_zono_approx, 4}, diff --git a/src/ode_solve.cpp b/src/ode_solve.cpp deleted file mode 100644 index eb062342..00000000 --- a/src/ode_solve.cpp +++ /dev/null @@ -1,207 +0,0 @@ -// [[Rcpp::depends(BH)]] - -// VolEsti (volume computation and sampling library) - -// Copyright (c) 2012-2020 Vissarion Fisikopoulos -// Copyright (c) 2018-2020 Apostolos Chalkis -// Copyright (c) 2020-2020 Marios Papachristou - -//Contributed and/or modified by Marios Papachristou, as part of Google Summer of Code 2018 and 2019 program. - - -#include -#include -#include -#include -#include -#include -#include -#include "random_walks/random_walks.hpp" -#include "volume/volume_sequence_of_balls.hpp" -#include "volume/volume_cooling_gaussians.hpp" -#include "sampling/sampling.hpp" -#include "ode_solvers/ode_solvers.hpp" -#include "oracle_functors_rcpp.h" - -template < - typename Solver, - typename MT -> -void run_ode_solver( - Solver &solver, - unsigned int &order, - unsigned int &num_steps, - unsigned int &dimension, - Rcpp::List &results) { - - std::vector results_temp; - - for (unsigned int i = 0; i < order; i++) { - MT temp_result; - temp_result.resize(dimension, num_steps); - results_temp.push_back(temp_result); - } - - for (unsigned int i = 0; i < num_steps; i++) { - for (unsigned int j = 0; j < order; j++) { - results_temp[j].col(i) = solver.xs[j].getCoefficients(); - } - solver.step(i, true); - } - - for (unsigned int i = 0; i < order; i++) { - std::ostringstream stringStream; - stringStream << "x_" << i + 1; - std::string state_name = stringStream.str(); - - results.push_back(Rcpp::wrap(results_temp[i]), state_name.c_str()); - - } - -} - -//' Solve an ODE of the form dx^n / dt^n = F(x, t) -//' -//' @param n The number of steps. -//' @param step_size The step size. -//' @param order The ODE order (default is n = 1) -//' @param dimension The dimension of each derivative -//' @param initial_time The initial time -//' @param F The function oracle F(x, t) in the ODE. -//' @param method The method to be used -//' @param initial_conditions The initial conditions provided to the solver. Must be provided in a list with keys "x_1", ..., "x_n" and column vectors as values. The state "x_n" represents the (n-1)-th order derivative with respect to time -//' @param domains A list of n H-polytopes with keys "P_1", "P_2", ..., "P_n" that correspond to each derivative's domain -//' -//' @return A list which contains elements "x_1", ..., "x_n" representing each derivative results. Each "x_i" corresponds to a d x n matrix where each column represents a certain timestep of the solver. -//' -//' @examples -//' F <- function (x) (-x) -//' initial_conditions <- list("x_1" = c(0), "x_2" = c(1)) -//' states <- ode_solve(dimension=1, n=1000, F=F, initial_time=0, step_size=0.01, order=2, method="leapfrog", initial_conditions=initial_conditions, domains = list()) -//' -//' @export -// [[Rcpp::export]] -Rcpp::List ode_solve(Rcpp::Nullable n, - Rcpp::Nullable step_size, - Rcpp::Nullable order, - Rcpp::Nullable dimension, - Rcpp::Nullable initial_time, - Rcpp::Function F, - Rcpp::Nullable method, - Rcpp::Nullable domains = R_NilValue, - Rcpp::Nullable initial_conditions = R_NilValue) -{ - - typedef double NT; - typedef Cartesian Kernel; - typedef BoostRandomNumberGenerator RNGType; - typedef typename Kernel::Point Point; - typedef HPolytope Hpolytope; - typedef Eigen::Matrix VT; - typedef Eigen::Matrix MT; - typedef std::vector pts; - typedef std::vector hp_bounds; - typedef RcppFunctor::GradientFunctor func; - - unsigned int n_ = Rcpp::as(n); - NT step_size_ = Rcpp::as(step_size); - unsigned int order_ = Rcpp::as(order); - std::string method_ = Rcpp::as(method); - unsigned int dim = Rcpp::as(dimension); - NT initial_time_ = Rcpp::as(initial_time); - - // Create functors - RcppFunctor::parameters rcpp_functor_params(1, 1, order_); - - // Create C++ functor - func F_(rcpp_functor_params, F, false); - - // Initialize initial conditions - pts initial_conditions_; - VT temp_initial_condition; - - for (unsigned int i = 1; i <= order_; i++) { - std::ostringstream stringStream; - stringStream << "x_" << i; - std::string state_name = stringStream.str(); - - if (Rcpp::as(initial_conditions).containsElementNamed(state_name.c_str())) { - temp_initial_condition = Rcpp::as(Rcpp::as(initial_conditions)[state_name.c_str()]); - Point p(temp_initial_condition); - initial_conditions_.push_back(p); - } else { - Point p(dim); - initial_conditions_.push_back(p); - } - - } - - unsigned int i = 0; - unsigned int t = 0; - F_(i, initial_conditions_, t); - - hp_bounds domains_; - - // Initialize domains - // TODO Test and add other polytope types - for (unsigned int i = 0; i < order_; i++) { - std::ostringstream stringStream; - stringStream << "P_" << i + 1; - std::string domain_name = stringStream.str(); - - if (Rcpp::as(domains).containsElementNamed(domain_name.c_str())) { - - Hpolytope HP(dim, Rcpp::as( - Rcpp::as(Rcpp::as(domains) - [domain_name.c_str()]).slot("A")), - Rcpp::as( - Rcpp::as(Rcpp::as(domains) - [domain_name.c_str()]).slot("b")) - ); - - HP.normalize(); - - Hpolytope *HP_ref = &HP; - - if (!HP_ref->is_in(initial_conditions_[i])) { - std::ostringstream errStream; - errStream << "Initial condition out of bounds for state " << i + 1 << ".\n"; - throw Rcpp::exception(errStream.str().c_str()); - } - domains_.push_back(HP_ref); - } else { - domains_.push_back(NULL); - } - } - - Rcpp::List results; - - if (method_ == "euler") { - EulerODESolver euler_solver = - EulerODESolver - (initial_time_, step_size_, initial_conditions_, F_, domains_); - run_ode_solver, MT>(euler_solver, order_, n_, dim, results); - } else if (method_ == "leapfrog") { - if (order_ % 2 == 1) { - throw Rcpp::exception("Leapfrog is an even order solver."); - } - LeapfrogODESolver leapfrog_solver = - LeapfrogODESolver - (initial_time_, step_size_, initial_conditions_, F_, domains_); - run_ode_solver, MT>(leapfrog_solver, order_, n_, dim, results); - } else if (method_ == "runge_kutta") { - RKODESolver rk_solver = - RKODESolver - (initial_time_, step_size_, initial_conditions_, F_, domains_); - run_ode_solver, MT>(rk_solver, order_, n_, dim, results); - } else if (method_ == "richardson") { - RichardsonExtrapolationODESolver r_solver = - RichardsonExtrapolationODESolver - (initial_time_, step_size_, initial_conditions_, F_, domains_); - run_ode_solver, MT>(r_solver, order_, n_, dim, results); - } else { - throw Rcpp::exception("Unrecognized solver. Aborting."); - } - - return results; -} diff --git a/src/spectrahedron.cpp b/src/spectrahedron.cpp deleted file mode 100644 index f084fd65..00000000 --- a/src/spectrahedron.cpp +++ /dev/null @@ -1,117 +0,0 @@ -// VolEsti (volume computation and sampling library) - -// Copyright (c) 2012-2020 Vissarion Fisikopoulos -// Copyright (c) 2020 Apostolos Chalkis - -//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program. - -// Licensed under GNU LGPL.3, see LICENCE file - - -#include -#include -#include -#include "cartesian_geom/cartesian_kernel.h" -#include -#include -#include -#include -#include "convex_bodies/spectrahedra/LMI.h" -#include "convex_bodies/spectrahedra/spectrahedron.h" -#include "SDPAFormatManager.h" - -//' Write a SDPA format file -//' -//' Outputs a spectrahedron (the matrices defining a linear matrix inequality) and a vector (the objective function) -//' to a SDPA format file. -//' -//' @param spectrahedron A spectrahedron in n dimensions; must be an object of class Spectrahedron -//' @param objectiveFunction A numerical vector of length n -//' @param outputFile Name of the output file -//' -//' @examples -//' \dontrun{ -//' A0 = matrix(c(-1,0,0,0,-2,1,0,1,-2), nrow=3, ncol=3, byrow = TRUE) -//' A1 = matrix(c(-1,0,0,0,0,1,0,1,0), nrow=3, ncol=3, byrow = TRUE) -//' A2 = matrix(c(0,0,-1,0,0,0,-1,0,0), nrow=3, ncol=3, byrow = TRUE) -//' lmi = list(A0, A1, A2) -//' S = Spectrahedron(matrices = lmi) -//' objFunction = c(1,1) -//' writeSdpaFormatFile(S, objFunction, "output.txt") -//' } -//' @export -// [[Rcpp::export]] -void writeSdpaFormatFile(Rcpp::Reference spectrahedron = R_NilValue, - Rcpp::Nullable objectiveFunction = R_NilValue, - Rcpp::Nullable outputFile = R_NilValue) { - - typedef double NT; - typedef Eigen::Matrix VT; - typedef Eigen::Matrix MT; - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef LMI LMI; - typedef Spectrahedron Spectrahedron; - - std::vector matrices = Rcpp::as > (Rcpp::as (spectrahedron).slot("matrices")); - LMI lmi(matrices); - Spectrahedron _spectrahedron(lmi); - Point c(Rcpp::as (objectiveFunction)); - - std::filebuf fb; - fb.open(Rcpp::as (outputFile), std::ios::out); - std::ostream os(&fb); - - SdpaFormatManager sdpaFormatManager; - sdpaFormatManager.writeSDPAFormatFile(os, _spectrahedron, c); - - return; -} - - - -//' Read a SDPA format file -//' -//' @param inputFile Name of the input file -//' -//' @return A list with two named items: an item "matrices" which is a list of the matrices and an vector "objFunction" -//' -//' @examples -//' path = system.file('extdata', package = 'volesti') -//' l = loadSdpaFormatFile(paste0(path,'/sdpa_n2m3.txt')) -//' @export -// [[Rcpp::export]] -Rcpp::List loadSdpaFormatFile(Rcpp::Nullable inputFile = R_NilValue) { - - typedef double NT; - typedef Eigen::Matrix VT; - typedef Eigen::Matrix MT; - typedef Cartesian Kernel; - typedef typename Kernel::Point Point; - typedef boost::mt19937 RNGType; - typedef LMI LMI; - typedef Spectrahedron Spectrahedron; - - Spectrahedron _spectrahedron; - Point c; - - // open stream - std::ifstream os; - os.open(Rcpp::as (inputFile),std::ifstream::in); - - // read file - SdpaFormatManager sdpaFormatManager; - sdpaFormatManager.loadSDPAFormatFile(os, _spectrahedron, c); - - std::vector matrices = _spectrahedron.getLMI().getMatrices(); - - // return spectrahedron and objective function - Rcpp::List _matrices; - - for (auto matrix : matrices) - _matrices.push_back(Rcpp::wrap(matrix)); - - Rcpp::List retList = Rcpp::List::create(Rcpp::Named("matrices") = _matrices , Rcpp::_["objFunction"] = Rcpp::wrap(c.getCoefficients())); - return Rcpp::wrap(retList); -}