Skip to content

Commit 05a3b6a

Browse files
authored
Merge pull request #16 from vfisikop/several_fixes
Several fixes
2 parents 001005f + 87cb409 commit 05a3b6a

22 files changed

+358
-705
lines changed

DESCRIPTION

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ Version: 1.2.0
1616
Date: 2024-02-29
1717
Maintainer: Vissarion Fisikopoulos <[email protected]>
1818
Depends: Rcpp (>= 0.12.17)
19-
Imports: methods, stats
19+
Imports: methods, stats, Matrix
2020
LinkingTo: Rcpp, RcppEigen, BH
2121
Suggests: testthat
2222
Encoding: UTF-8

NAMESPACE

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -19,8 +19,6 @@ export(gen_simplex)
1919
export(gen_skinny_cube)
2020
export(geweke)
2121
export(inner_ball)
22-
export(loadSdpaFormatFile)
23-
export(ode_solve)
2422
export(pinvweibull_with_loc)
2523
export(psrf_multivariate)
2624
export(psrf_univariate)
@@ -30,10 +28,10 @@ export(rotate_polytope)
3028
export(round_polytope)
3129
export(sample_points)
3230
export(volume)
33-
export(writeSdpaFormatFile)
3431
export(write_sdpa_format_file)
3532
export(zonotope_approximation)
3633
exportClasses(Hpolytope)
34+
exportClasses(HpolytopeSparse)
3735
exportClasses(Spectrahedron)
3836
exportClasses(Vpolytope)
3937
exportClasses(VpolytopeIntersection)

NEWS.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -54,7 +54,7 @@
5454
# volesti 1.2.0
5555

5656
- New functions: dinvweibull_with_loc, ess, estimtate_lipschitz_constant, gen_birkhoff, geweke
57-
ode_solve, pinvweibull_with_loc, psrf_multivariate, psrf_univariate, raftery
57+
pinvweibull_with_loc, psrf_multivariate, psrf_univariate, raftery
5858

5959
- New features in sample_points function:
6060
a) new walks: i) Dikin walk, ii) Vaidya walk, iii) John walk,

R/HpolytopeSparseClass.R

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
#' An R class to represent an H-polytope defined by a sparse matrix
2+
#'
3+
#' A sparse H-polytope is a convex polytope defined by a set of linear inequalities or equivalently a
4+
#' \eqn{d}-dimensional H-polytope with \eqn{m} facets is defined by a \eqn{m\times d} matrix A and
5+
#' a \eqn{m}-dimensional vector b, s.t.: \eqn{Ax\leq b}, where \eqn{A} is sparse.
6+
#'
7+
#' \describe{
8+
#' \item{Aineq}{An \eqn{m\times d} sparse numerical matrix for the inequalities.}
9+
#'
10+
#' \item{bineq}{An \eqn{m}-dimensional vector for the ineqaulities.}
11+
#'
12+
#' \item{Aineq}{An \eqn{m'\times d} sparse numerical matrix for the equalities.}
13+
#'
14+
#' \item{bineq}{An \eqn{m'}-dimensional vector for the eqaulities.}
15+
#'
16+
#' \item{lb}{\eqn{d}-dimensional vector lb.}
17+
#'
18+
#' \item{ub}{\eqn{d}-dimensional vector ub.}
19+
#'
20+
#' \item{type}{A character with default value 'HpolytopeSparse', to declare the representation of the polytope.}
21+
#' }
22+
#'
23+
#' @examples
24+
#' library(Matrix)
25+
#' bineq=c(10,10,10,10,10)
26+
#' Aineq = matrix(c(1,0,-0.25,-1,2.5,1,0.4,-1,-0.9,0.5), nrow=5, ncol=2, byrow = TRUE)
27+
#' Aineq = as( Aineq, 'dgCMatrix' )
28+
#' beq=(0)
29+
#' Aeq = matrix(, nrow=0, ncol=2, byrow = TRUE)
30+
#' Aeq=as( Aeq, 'dgCMatrix' )
31+
#' lb=-100000*c(1,1);
32+
#' ub=100000*c(1,1);
33+
#' P <- HpolytopeSparse(Aineq = Aineq, bineq = bineq, Aeq = Aeq, beq = beq, lb = lb, ub = ub)
34+
#'
35+
#' @name HpolytopeSparse-class
36+
#' @rdname HpolytopeSparse-class
37+
#' @exportClass HpolytopeSparse
38+
#' @importFrom "Matrix" "CsparseMatrix"
39+
HpolytopeSparse <- setClass (
40+
# Class name
41+
"HpolytopeSparse",
42+
43+
# Defining slot type
44+
representation (
45+
Aineq = "sparseMatrix",
46+
bineq = "numeric",
47+
Aeq = "sparseMatrix",
48+
beq = "numeric",
49+
lb = "numeric",
50+
ub = "numeric",
51+
type = "character"
52+
),
53+
54+
# Initializing slots
55+
prototype = list(
56+
Aineq = Matrix::sparseMatrix(i = integer(0), j = integer(0), x = double(0), dims = c(0, 0)),
57+
bineq = numeric(0),
58+
Aeq = Matrix::sparseMatrix(i = integer(0), j = integer(0), x = double(0), dims = c(0, 0)),
59+
beq = numeric(0),
60+
lb = numeric(0),
61+
ub = numeric(0),
62+
type = "HpolytopeSparse"
63+
)
64+
)

R/RcppExports.R

Lines changed: 12 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -187,30 +187,6 @@ load_sdpa_format_file <- function(input_file = NULL) {
187187
.Call(`_volesti_load_sdpa_format_file`, input_file)
188188
}
189189

190-
#' Solve an ODE of the form dx^n / dt^n = F(x, t)
191-
#'
192-
#' @param n The number of steps.
193-
#' @param step_size The step size.
194-
#' @param order The ODE order (default is n = 1)
195-
#' @param dimension The dimension of each derivative
196-
#' @param initial_time The initial time
197-
#' @param F The function oracle F(x, t) in the ODE.
198-
#' @param method The method to be used
199-
#' @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
200-
#' @param domains A list of n H-polytopes with keys "P_1", "P_2", ..., "P_n" that correspond to each derivative's domain
201-
#'
202-
#' @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.
203-
#'
204-
#' @examples
205-
#' F <- function (x) (-x)
206-
#' initial_conditions <- list("x_1" = c(0), "x_2" = c(1))
207-
#' 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())
208-
#'
209-
#' @export
210-
ode_solve <- function(n, step_size, order, dimension, initial_time, F, method, domains = NULL, initial_conditions = NULL) {
211-
.Call(`_volesti_ode_solve`, n, step_size, order, dimension, initial_time, F, method, domains, initial_conditions)
212-
}
213-
214190
#' An internal Rccp function as a polytope generator
215191
#'
216192
#' @param kind_gen An integer to declare the type of the polytope.
@@ -311,7 +287,17 @@ rounding <- function(P, method = NULL, seed = NULL) {
311287
#' @param n The number of points that the function is going to sample from the convex polytope.
312288
#' @param random_walk Optional. A list that declares the random walk and some related parameters as follows:
313289
#' \itemize{
314-
#' \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.}
290+
#' \item{\code{walk} }{ A string to declare the random walk: i) \code{'CDHR'} for Coordinate Directions Hit-and-Run,
291+
#' ii) \code{'RDHR'} for Random Directions Hit-and-Run, iii) \code{'BaW'} for Ball Walk, iv) \code{'BiW'} for Billiard walk,
292+
#' v) \code{'dikin'} for dikin walk, vi) \code{'vaidya'} for vaidya walk, vii) \code{'john'} for john walk,
293+
#' 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,
294+
#' x) \code{'NUTS'} for NUTS Hamiltonian Monte Carlo sampler (logconcave densities), xi) \code{'HMC'} for Hamiltonian Monte Carlo (logconcave densities),
295+
#' xii) CRHMC for Riemannian HMC with H-polytope constraints (uniform and general logconcave densities),
296+
#' xiii) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method (logconcave densities),
297+
#' xiii) \code{'ExactHMC'} for exact Hamiltonian Monte Carlo with reflections (spherical Gaussian or exponential distribution).
298+
#' The default walk is \code{'aBiW'} for the uniform distribution, \code{'CDHR'} for the Gaussian distribution and H-polytopes and
299+
#' \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'}
300+
#' for logconcave densities with H-polytope and sparse constrainted problems.}
315301
#' \item{\code{walk_length} }{ The number of the steps per generated point for the random walk. The default value is \eqn{1}.}
316302
#' \item{\code{nburns} }{ The number of points to burn before start sampling. The default value is \eqn{1}.}
317303
#' \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) {
360346
#' # gaussian distribution from the 2d unit simplex in H-representation with variance = 2
361347
#' A = matrix(c(-1,0,0,-1,1,1), ncol=2, nrow=3, byrow=TRUE)
362348
#' b = c(0,0,1)
363-
#' P = Hpolytope(A = A, b = b)
349+
#' P = Hpolytope(A=A,b=b)
364350
#' points = sample_points(P, n = 100, distribution = list("density" = "gaussian", "variance" = 2))
365351
#'
366352
#' # 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 =
374360
.Call(`_volesti_sample_points`, P, n, random_walk, distribution, seed)
375361
}
376362

377-
#' Write a SDPA format file
378-
#'
379-
#' Outputs a spectrahedron (the matrices defining a linear matrix inequality) and a vector (the objective function)
380-
#' to a SDPA format file.
381-
#'
382-
#' @param spectrahedron A spectrahedron in n dimensions; must be an object of class Spectrahedron
383-
#' @param objectiveFunction A numerical vector of length n
384-
#' @param outputFile Name of the output file
385-
#'
386-
#' @examples
387-
#' \dontrun{
388-
#' A0 = matrix(c(-1,0,0,0,-2,1,0,1,-2), nrow=3, ncol=3, byrow = TRUE)
389-
#' A1 = matrix(c(-1,0,0,0,0,1,0,1,0), nrow=3, ncol=3, byrow = TRUE)
390-
#' A2 = matrix(c(0,0,-1,0,0,0,-1,0,0), nrow=3, ncol=3, byrow = TRUE)
391-
#' lmi = list(A0, A1, A2)
392-
#' S = Spectrahedron(matrices = lmi)
393-
#' objFunction = c(1,1)
394-
#' writeSdpaFormatFile(S, objFunction, "output.txt")
395-
#' }
396-
#' @export
397-
writeSdpaFormatFile <- function(spectrahedron = NULL, objectiveFunction = NULL, outputFile = NULL) {
398-
invisible(.Call(`_volesti_writeSdpaFormatFile`, spectrahedron, objectiveFunction, outputFile))
399-
}
400-
401-
#' Read a SDPA format file
402-
#'
403-
#' @param inputFile Name of the input file
404-
#'
405-
#' @return A list with two named items: an item "matrices" which is a list of the matrices and an vector "objFunction"
406-
#'
407-
#' @examples
408-
#' path = system.file('extdata', package = 'volesti')
409-
#' l = loadSdpaFormatFile(paste0(path,'/sdpa_n2m3.txt'))
410-
#' @export
411-
loadSdpaFormatFile <- function(inputFile = NULL) {
412-
.Call(`_volesti_loadSdpaFormatFile`, inputFile)
413-
}
414-
415363
#' 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
416364
#'
417365
#' 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.

R/zzz.R

Lines changed: 0 additions & 11 deletions
This file was deleted.

README.md

Lines changed: 38 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -5,58 +5,68 @@ The `volesti` package provides [R](https://www.r-project.org/) with functions fo
55

66
`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.
77

8+
9+
The latest **stable** version is available from CRAN.
10+
11+
[![CRAN status](https://www.r-pkg.org/badges/version/volesti)](https://cran.r-project.org/package=volesti)
12+
[![CRAN check](https://badges.cranchecks.info/worst/volesti.svg)](https://cran.r-project.org/web/checks/check_results_volesti.html)
13+
[![CRAN downloads](https://cranlogs.r-pkg.org/badges/volesti)](https://cran.r-project.org/package=volesti)
14+
![CRAN/METACRAN](https://img.shields.io/cran/l/volesti)
15+
16+
The latest **development** version is available in this repository.
17+
18+
[![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)
19+
[![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)
20+
[![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)
21+
22+
823
## Download and install
924

10-
* The latest stable version is available from CRAN.
11-
* The latest development version is available on Github `https://github.com/GeomScale/Rvolesti`
25+
To use the development version you need to follow these steps:
1226

13-
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:
27+
* Clone the `Rvolesti` repository.
1428

15-
* Clone the main `Rvolesti` repository.
16-
* Fetch the submodule from the [volesti](https://github.com/GeomScale/volesti) repository:
17-
```
18-
git submodule update --recursive --init --remote
19-
```
20-
* Build package by running:
21-
```
29+
* Build and install the package by running (from an `R` terminal):
30+
```R
2231
Rcpp::compileAttributes()
2332
devtools::build()
24-
```
25-
* Install Rvolesti by running:
26-
```
2733
install.packages("volesti")
2834
```
35+
or from a bash terminal:
36+
```bash
37+
Rscript -e 'Rcpp::compileAttributes()'
38+
R CMD INSTALL --no-multiarch --with-keep.source .
39+
```
2940

30-
The package-dependencies are: `Rcpp`, `RcppEigen`, `BH`.
41+
The following packages should be installed: `Rcpp`, `RcppEigen`, `BH`.
3142

3243
## Documentation
3344

3445
* [Using the R Interface](https://github.com/GeomScale/volesti/blob/v1.1.1/doc/r_interface.md)
3546
* [Wikipage with Tutorials and Demos](https://github.com/GeomScale/volesti/wiki)
3647
* [Tutorial given to PyData meetup](https://vissarion.github.io/tutorials/volesti_tutorial_pydata.html)
3748

38-
## How to update the volesti R package?
39-
40-
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:
49+
The user can generate or update the documentation:
4150

42-
- Update the `cran_include` branch in [volesti](https://github.com/GeomScale/volesti)
43-
- Clone the main `volesti` repository
44-
- Checkout to the `cran_include` branch
45-
- Update your code and open a PR similar to https://github.com/GeomScale/volesti/pull/277
46-
- Retrieve the new `include` directory using submodule
47-
```
48-
git submodule update --recursive --init --remote
51+
```R
52+
Rcpp::compileAttributes() # updates the Rcpp layer from C++ to R
53+
roxygen2::roxygenize(roclets="rd") # updates the docs based on roxygen comments
4954
```
5055

51-
*Note:* it is possible the this update will brake the R interface, thus this operation should be processed with care.
56+
and generate a pdf with the docs:
57+
58+
```R
59+
pack = "volesti"
60+
path = find.package(pack)
61+
system(paste(shQuote(file.path(R.home("bin"), "R")), "CMD", "Rd2pdf", shQuote(path)))
62+
```
5263

5364
## Credits
5465

5566
* [Contributors and Package History](https://github.com/GeomScale/volesti/blob/v1.1.1/doc/credits.md)
5667
* [List of Publications](https://github.com/GeomScale/volesti/blob/v1.1.1/doc/publications.md)
5768

58-
Copyright (c) 2012-2023 Vissarion Fisikopoulos
59-
60-
Copyright (c) 2018-2023 Apostolos Chalkis
69+
Copyright (c) 2012-2024 Vissarion Fisikopoulos
70+
Copyright (c) 2018-2024 Apostolos Chalkis
6171

6272
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/HpolytopeSparse-class.Rd

Lines changed: 42 additions & 0 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

man/examples/nuts_rand_poly.R

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -25,8 +25,8 @@ f <- function(x) (norm_vec(x)^2 + sum(x))
2525
# Negative log-probability gradient oracle
2626
grad_f <- function(x) (2 * x + 1)
2727

28-
dimension <- 50
29-
facets <- 200
28+
dimension <- 5
29+
facets <- 20
3030

3131
# Create domain of truncation
3232
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])
3939
x_min = matrix(0, dimension, 1)
4040

4141
# Warm start point from truncated Gaussian
42-
warm_start <- sample_points(P, n = 1, random_walk = list("nburns" = 5000), distribution = list("density" = "gaussian", "variance" = 1/2, "mode" = x_min))
42+
warm_start <- sample_points(P, n = 1, random_walk = list("nburns" = 500),
43+
distribution = list("density" = "gaussian", "variance" = 1/2, "mode" = x_min))
4344

4445
# Sample points
45-
n_samples <- 20000
46+
n_samples <- 1000
4647

4748
samples <- sample_points(P, n = n_samples, random_walk = list("walk" = "NUTS", "solver" = "leapfrog", "starting_point" = warm_start[,1]),
4849
distribution = list("density" = "logconcave", "negative_logprob" = f, "negative_logprob_gradient" = grad_f))
@@ -52,4 +53,6 @@ hist(samples[1,], probability=TRUE, breaks = 100)
5253

5354
psrfs <- psrf_univariate(samples)
5455
n_ess <- ess(samples)
55-
56+
cat("ESS=", n_ess, "\n")
57+
cat("\n")
58+
cat("PSRF=", psrfs, "\n")

0 commit comments

Comments
 (0)