Several fixes

DESCRIPTION

@@ -16,7 +16,7 @@ Version: 1.2.0
 Date: 2024-02-29
 Maintainer: Vissarion Fisikopoulos <[email protected]>
 Depends: Rcpp (>= 0.12.17)
-Imports: methods, stats
+Imports: methods, stats, Matrix
 LinkingTo: Rcpp, RcppEigen, BH
 Suggests: testthat
 Encoding: UTF-8

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)

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,

R/HpolytopeSparseClass.R

@@ -0,0 +1,64 @@
+#' An R class to represent an H-polytope defined by a sparse matrix
+#'
+#' 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 \eqn{A} is sparse.
+#'
+#' \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)
+#'
+#' @name HpolytopeSparse-class
+#' @rdname HpolytopeSparse-class
+#' @exportClass HpolytopeSparse
+#' @importFrom "Matrix" "CsparseMatrix"
+HpolytopeSparse <- setClass (
+  # Class name
+  "HpolytopeSparse",
+
+  # Defining slot type
+  representation (
+    Aineq = "sparseMatrix",
+    bineq = "numeric",
+    Aeq = "sparseMatrix",
+    beq = "numeric",
+    lb = "numeric",
+    ub = "numeric",
+    type = "character"
+    ),
+
+  # Initializing slots
+  prototype = list(
+    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"
+  )
+)

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.
@@ -311,7 +287,17 @@ rounding <- function(P, method = NULL, seed = NULL) {
 #' @param n The number of points that the function is going to sample from the convex polytope.
 #' @param 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()}.}
@@ -360,7 +346,7 @@ rounding <- function(P, method = NULL, seed = NULL) {
 #' # 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
@@ -374,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.

README.md

@@ -5,58 +5,68 @@ The `volesti` package provides [R](https://www.r-project.org/) with functions fo
 
 `volesti` computes approximations of volume of polytopes given as a set of points or linear inequalities or as a Minkowski sum of segments (zonotopes). There are algorithms for volume approximation as well as algorithms for sampling, rounding and rotating polytopes. Last but not least, `volesti` provides implementations of geometric algorithms to compute the score of a portfolio given asset returns and to detect financial crises in stock markets.
 
+
+The latest **stable** version is available from CRAN.
+
+[![CRAN status](https://www.r-pkg.org/badges/version/volesti)](https://cran.r-project.org/package=volesti)
+[![CRAN check](https://badges.cranchecks.info/worst/volesti.svg)](https://cran.r-project.org/web/checks/check_results_volesti.html)
+[![CRAN downloads](https://cranlogs.r-pkg.org/badges/volesti)](https://cran.r-project.org/package=volesti)
+![CRAN/METACRAN](https://img.shields.io/cran/l/volesti)
+
+The latest **development** version is available in this repository.
+
+[![R-CMD-ubuntu](https://github.com/GeomScale/Rvolesti/workflows/R-CMD-check-ubuntu/badge.svg)](https://github.com/GeomScale/Rvolesti/actions?query=workflow%3AR-CMD-ubuntu)
+[![R-CMD-macOS](https://github.com/GeomScale/Rvolesti/workflows/R-CMD-check-macOS/badge.svg)](https://github.com/GeomScale/Rvolesti/actions?query=workflow%3AR-CMD-macOS)
+[![R-CMD-windows](https://github.com/GeomScale/Rvolesti/workflows/R-CMD-check-windows/badge.svg)](https://github.com/GeomScale/Rvolesti/actions?query=workflow%3AR-CMD-windows)
+
+
 ##  Download and install
 
-* The latest stable version is available from CRAN.
-* The latest development version is available on Github `https://github.com/GeomScale/Rvolesti`
+To use the development version you need to follow these steps:
 
-To use the development version of `volesti` that includes the C++ code, you need to clone the repository and fetch the submodule. Follow these steps:
+* Clone the `Rvolesti` repository.
 
-* Clone the main `Rvolesti` repository.
-* Fetch the submodule from the [volesti](https://github.com/GeomScale/volesti) repository:
-```
-git submodule update --recursive --init --remote
-```
-* Build package by running:
-```
+* Build and install the package by running (from an `R` terminal):
+```R
 Rcpp::compileAttributes()
 devtools::build()
-```
-* Install Rvolesti by running:
-```
 install.packages("volesti")
 ```
+or from a bash terminal:
+```bash
+Rscript -e 'Rcpp::compileAttributes()'
+R CMD INSTALL --no-multiarch --with-keep.source .
+```
 
-The package-dependencies are: `Rcpp`, `RcppEigen`, `BH`.
+The following packages should be installed: `Rcpp`, `RcppEigen`, `BH`.
 
 ## Documentation
 
 * [Using the R Interface](https://github.com/GeomScale/volesti/blob/v1.1.1/doc/r_interface.md)
 * [Wikipage with Tutorials and Demos](https://github.com/GeomScale/volesti/wiki)
 * [Tutorial given to PyData meetup](https://vissarion.github.io/tutorials/volesti_tutorial_pydata.html)
 
-## How to update the volesti R package?
-
-The C++ source code is retrieved from [volesti](https://github.com/GeomScale/volesti) package and placed in `src/include` of the current repository. To update the current C++ code we have to follow two steps:
+The user can generate or update the documentation:
 
-- Update the `cran_include` branch in [volesti](https://github.com/GeomScale/volesti)
-    - Clone the main `volesti` repository
-    - Checkout to the `cran_include` branch
-    - Update your code and open a PR similar to https://github.com/GeomScale/volesti/pull/277
-- Retrieve the new `include` directory using submodule
-```
-git submodule update --recursive --init --remote
+```R
+Rcpp::compileAttributes() # updates the Rcpp layer from C++ to R
+roxygen2::roxygenize(roclets="rd") # updates the docs based on roxygen comments
 ```
 
-*Note:* it is possible the this update will brake the R interface, thus this operation should be processed with care. 
+and generate a pdf with the docs:
+
+```R
+pack = "volesti"
+path = find.package(pack)
+system(paste(shQuote(file.path(R.home("bin"), "R")), "CMD", "Rd2pdf", shQuote(path)))
+```
 
 ## Credits
 
 * [Contributors and Package History](https://github.com/GeomScale/volesti/blob/v1.1.1/doc/credits.md)
 * [List of Publications](https://github.com/GeomScale/volesti/blob/v1.1.1/doc/publications.md)
 
-Copyright (c) 2012-2023 Vissarion Fisikopoulos
-
-Copyright (c) 2018-2023 Apostolos Chalkis
+Copyright (c) 2012-2024 Vissarion Fisikopoulos
+Copyright (c) 2018-2024 Apostolos Chalkis
 
 You may redistribute or modify the software under the GNU Lesser General Public License as published by Free Software Foundation, either version 3 of the License, or (at your option) any later version. It is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY.

man/examples/nuts_rand_poly.R

@@ -25,8 +25,8 @@ f <- function(x) (norm_vec(x)^2 + sum(x))
 # Negative log-probability gradient oracle
 grad_f <- function(x) (2 * x + 1)
 
-dimension <- 50
-facets <- 200
+dimension <- 5
+facets <- 20
 
 # Create domain of truncation
 H <- gen_rand_hpoly(dimension, facets)
@@ -39,10 +39,11 @@ P <- Hpolytope(A = Tr$Mat[1:nrow(Tr$Mat), 2:ncol(Tr$Mat)], b = Tr$Mat[,1])
 x_min = matrix(0, dimension, 1)
 
 # Warm start point from truncated Gaussian
-warm_start <- sample_points(P, n = 1, random_walk = list("nburns" = 5000), distribution = list("density" = "gaussian", "variance" = 1/2, "mode" = x_min))
+warm_start <- sample_points(P, n = 1, random_walk = list("nburns" = 500),
+                            distribution = list("density" = "gaussian", "variance" = 1/2, "mode" = x_min))
 
 # Sample points
-n_samples <- 20000
+n_samples <- 1000
 
 samples <- sample_points(P, n = n_samples, random_walk = list("walk" = "NUTS", "solver" = "leapfrog", "starting_point" = warm_start[,1]),
                          distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f))
@@ -52,4 +53,6 @@ hist(samples[1,], probability=TRUE, breaks = 100)
 
 psrfs <- psrf_univariate(samples)
 n_ess <- ess(samples)
-
+cat("ESS=", n_ess, "\n")
+cat("\n")
+cat("PSRF=", psrfs, "\n")