forked from GeomScale/volesti
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathRcppExports.R
475 lines (454 loc) · 28.7 KB
/
RcppExports.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
#' Construct a copula using uniform sampling from the unit simplex
#'
#' Given two families of parallel hyperplanes or a family of parallel hyperplanes and a family of concentric ellispoids centered at the origin intersecting the canonical simplex, this function uniformly samples from the canonical simplex and construct an approximation of the bivariate probability distribution, called copula (see \url{https://en.wikipedia.org/wiki/Copula_(probability_theory)}).
#' At least two families of hyperplanes or one family of hyperplanes and one family of ellipsoids have to be given as input.
#'
#' @param r1 The \eqn{d}-dimensional normal vector of the first family of parallel hyperplanes.
#' @param r2 Optional. The \eqn{d}-dimensional normal vector of the second family of parallel hyperplanes.
#' @param sigma Optional. The \eqn{d\times d} symmetric positive semidefine matrix that describes the family of concentric ellipsoids centered at the origin.
#' @param m The number of the slices for the copula. The default value is 100.
#' @param n The number of points to sample. The default value is \eqn{5\cdot 10^5}.
#' @param seed Optional. A fixed seed for the number generator.
#'
#' @references \cite{L. Cales, A. Chalkis, I.Z. Emiris, V. Fisikopoulos,
#' \dQuote{Practical volume computation of structured convex bodies, and an application to modeling portfolio dependencies and financial crises,} \emph{Proc. of Symposium on Computational Geometry, Budapest, Hungary,} 2018.}
#'
#' @return A \eqn{m\times m} numerical matrix that corresponds to a copula.
#' @examples
#' # compute a copula for two random families of parallel hyperplanes
#' h1 = runif(n = 10, min = 1, max = 1000)
#' h1 = h1 / 1000
#' h2=runif(n = 10, min = 1, max = 1000)
#' h2 = h2 / 1000
#' cop = copula(r1 = h1, r2 = h2, m = 10, n = 100000)
#'
#' # compute a copula for a family of parallel hyperplanes and a family of conentric ellipsoids
#' h = runif(n = 10, min = 1, max = 1000)
#' h = h / 1000
#' E = replicate(10, rnorm(20))
#' E = cov(E)
#' cop = copula(r1 = h, sigma = E, m = 10, n = 100000)
#'
#' @export
copula <- function(r1, r2 = NULL, sigma = NULL, m = NULL, n = NULL, seed = NULL) {
.Call(`_volesti_copula`, r1, r2, sigma, m, n, seed)
}
#' Sample perfect uniformly distributed points from well known convex bodies: (a) the unit simplex, (b) the canonical simplex, (c) the boundary of a hypersphere or (d) the interior of a hypersphere.
#'
#' The \eqn{d}-dimensional unit simplex is the set of points \eqn{\vec{x}\in \R^d}, s.t.: \eqn{\sum_i x_i\leq 1}, \eqn{x_i\geq 0}. The \eqn{d}-dimensional canonical simplex is the set of points \eqn{\vec{x}\in \R^d}, s.t.: \eqn{\sum_i x_i = 1}, \eqn{x_i\geq 0}.
#'
#' @param body A list to request exact uniform sampling from special well known convex bodies through the following input parameters:
#' \itemize{
#' \item{\code{type} }{ A string that declares the type of the body for the exact sampling: a) \code{'unit_simplex'} for the unit simplex, b) \code{'canonical_simplex'} for the canonical simplex, c) \code{'hypersphere'} for the boundary of a hypersphere centered at the origin, d) \code{'ball'} for the interior of a hypersphere centered at the origin.}
#' \item{\code{dimension} }{ An integer that declares the dimension when exact sampling is enabled for a simplex or a hypersphere.}
#' \item{\code{radius} }{ The radius of the \eqn{d}-dimensional hypersphere. The default value is \eqn{1}.}
#' }
#' @param n The number of points that the function is going to sample.
#' @param seed Optional. A fixed seed for the number generator.
#'
#' @references \cite{R.Y. Rubinstein and B. Melamed,
#' \dQuote{Modern simulation and modeling} \emph{ Wiley Series in Probability and Statistics,} 1998.}
#' @references \cite{A Smith, Noah and W Tromble, Roy,
#' \dQuote{Sampling Uniformly from the Unit Simplex,} \emph{ Center for Language and Speech Processing Johns Hopkins University,} 2004.}
#'
#' @return A \eqn{d\times n} matrix that contains, column-wise, the sampled points from the convex polytope P.
#' @examples
#' # 100 uniform points from the 2-d unit ball
#' points = direct_sampling(n = 100, body = list("type" = "ball", "dimension" = 2))
#' @export
direct_sampling <- function(body, n, seed = NULL) {
.Call(`_volesti_direct_sampling`, body, n, seed)
}
#' Gelman-Rubin and Brooks-Gelman Potential Scale Reduction Factor (PSRF) for each marginal
#'
#' @param samples A matrix that contans column-wise the sampled points from a geometric random walk.
#' @param method A string to reauest diagnostic: (i) \code{'normal'} for psrf of Gelman-Rubin and (ii) \code{'interval'} for psrf of Brooks-Gelman.
#'
#' @references \cite{Gelman, A. and Rubin, D. B.,
#' \dQuote{Inference from iterative simulation using multiple sequences,} \emph{Statistical Science,} 1992.}
#'
#' @references \cite{Brooks, S. and Gelman, A.,
#' \dQuote{General Methods for Monitoring Convergence of Iterative Simulations,} \emph{Journal of Computational and Graphical Statistics,} 1998.}
#'
#' @return A vector that contains the values of PSRF for each coordinate
#'
#' @export
ess <- function(samples) {
.Call(`_volesti_ess`, samples)
}
#' Compute the exact volume of (a) a zonotope (b) an arbitrary simplex in V-representation or (c) if the volume is known and declared by the input object.
#'
#' Given a zonotope (as an object of class Zonotope), this function computes the sum of the absolute values of the determinants of all the \eqn{d \times d} submatrices of the \eqn{m\times d} matrix \eqn{G} that contains row-wise the \eqn{m} \eqn{d}-dimensional segments that define the zonotope.
#' For an arbitrary simplex that is given in V-representation this function computes the absolute value of the determinant formed by the simplex's points assuming it is shifted to the origin.
#'
#' @param P A polytope
#'
#' @references \cite{E. Gover and N. Krikorian,
#' \dQuote{Determinants and the Volumes of Parallelotopes and Zonotopes,} \emph{Linear Algebra and its Applications, 433(1), 28 - 40,} 2010.}
#'
#' @return The exact volume of the input polytope, for zonotopes, simplices in V-representation and polytopes with known exact volume
#' @examples
#'
#' # compute the exact volume of a 5-dimensional zonotope defined by the Minkowski sum of 10 segments
#' Z = gen_rand_zonotope(2, 5)
#' vol = exact_vol(Z)
#'
#' \donttest{# compute the exact volume of a 2-d arbitrary simplex
#' V = matrix(c(2,3,-1,7,0,0),ncol = 2, nrow = 3, byrow = TRUE)
#' P = Vpolytope$new(V)
#' vol = exact_vol(P)
#' }
#'
#' # compute the exact volume the 10-dimensional cross polytope
#' P = gen_cross(10,'V')
#' vol = exact_vol(P)
#' @export
exact_vol <- function(P) {
.Call(`_volesti_exact_vol`, P)
}
#' Compute the percentage of the volume of the simplex that is contained in the intersection of a half-space and the simplex.
#'
#' A half-space \eqn{H} is given as a pair of a vector \eqn{a\in R^d} and a scalar \eqn{z0\in R} s.t.: \eqn{a^Tx\leq z0}. This function calls the Ali's version of the Varsi formula to compute a frustum of the simplex.
#'
#' @param a A \eqn{d}-dimensional vector that defines the direction of the hyperplane.
#' @param z0 The scalar that defines the half-space.
#'
#' @references \cite{Varsi, Giulio,
#' \dQuote{The multidimensional content of the frustum of the simplex,} \emph{Pacific J. Math. 46, no. 1, 303--314,} 1973.}
#'
#' @references \cite{Ali, Mir M.,
#' \dQuote{Content of the frustum of a simplex,} \emph{ Pacific J. Math. 48, no. 2, 313--322,} 1973.}
#'
#' @return The percentage of the volume of the simplex that is contained in the intersection of a given half-space and the simplex.
#'
#' @examples
#' # compute the frustum of H: -x1+x2<=0
#' a=c(-1,1)
#' z0=0
#' frustum = frustum_of_simplex(a, z0)
#' @export
frustum_of_simplex <- function(a, z0) {
.Call(`_volesti_frustum_of_simplex`, a, z0)
}
#' Internal rcpp function to compute the full dimensional polytope when a low dimensional is given
#'
#' @param P A low dimensional convex polytope in H-representation.
#'
#' @keywords internal
#'
#' @return A numerical matrix that describes the full dimensional polytope, a numerical matrix of the inverse
#' linear transformation that is applied on the input polytope, the numerical vector - point that the
#' input polytope is shifted and the product of the singular values of the matrix of the linear map
#' applied on the input polytope.
#'
full_dimensional_polytope <- function(P) {
.Call(`_volesti_full_dimensional_polytope`, P)
}
#' Geweke's MCMC diagnostic
#'
#' @param samples A matrix that contans column-wise the sampled points from a geometric random walk.
#' @param frac_first Optional. The portion of the first in order points in matrix samples.
#' @param frac_last Optional. The portion of the last in order points in matrix samples.
#'
#' @references \cite{Geweke, J.,
#' \dQuote{Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments,} \emph{ In Bayesian Statistics 4. Proceedings of the Fourth Valencia International Meeting,} 1992.}
#'
#' @return A boolean to denote if the result of Geweke diagnostic: (i) false if the null hypothesis is rejected, (ii) true if the null hypothesis is not rejected.
#'
#' @export
geweke <- function(samples, frac_first = NULL, frac_last = NULL) {
.Call(`_volesti_geweke`, samples, frac_first, frac_last)
}
#' Compute an inscribed ball of a convex polytope
#'
#' For a H-polytope 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\} }, this function computes the largest inscribed ball (Chebychev ball) by solving the corresponding linear program.
#' For both zonotopes and V-polytopes the function computes the minimum \eqn{r} s.t.: \eqn{ r e_i \in P} for all \eqn{i=1, \dots ,d}. Then the ball centered at the origin with radius \eqn{r/ \sqrt{d}} is an inscribed ball.
#'
#' @param P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope or (d) VpolytopeIntersection.
#' @param lpsolve Optional. A boolean variable to compute the Chebychev ball of an H-polytope using the lpsolve library.
#'
#' @return A \eqn{(d+1)}-dimensional vector that describes the inscribed ball. The first \eqn{d} coordinates corresponds to the center of the ball and the last one to the radius.
#'
#' @examples
#' # compute the Chebychev ball of the 2d unit simplex
#' P = gen_cube(10,'H')
#' ball_vec = inner_ball(P)
#'
#' # compute an inscribed ball of the 3-dimensional unit cube in V-representation
#' P = gen_cube(3, 'V')
#' ball_vec = inner_ball(P, lpsolve = TRUE)
#' @export
inner_ball <- function(P, lpsolve = NULL) {
.Call(`_volesti_inner_ball`, P, lpsolve)
}
#' 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
#' # Please visit the examples directory on examples demonstrating usage of the ODE solvers.
#'
#' @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.
#' @param Vpoly_gen A boolean parameter to declare if the requested polytope has to be in V-representation.
#' @param Zono_gen A boolean parameter to declare if the requested polytope has to be a zonotope.
#' @param dim_gen An integer to declare the dimension of the requested polytope.
#' @param m_gen An integer to declare the number of generators for the requested random zonotope or the number of vertices for a V-polytope.
#' @param seed Optional. A fixed seed for the random polytope generator.
#'
#' @keywords internal
#'
#' @return A numerical matrix describing the requested polytope
poly_gen <- function(kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen, seed = NULL) {
.Call(`_volesti_poly_gen`, kind_gen, Vpoly_gen, Zono_gen, dim_gen, m_gen, seed)
}
#' Gelman-Rubin Potential Scale Reduction Factor (PSRF)
#'
#' @param samples A matrix that contans column-wise the sampled points from a geometric random walk.
#'
#' @references \cite{Gelman, A. and Rubin, D. B.,
#' \dQuote{Inference from iterative simulation using multiple sequences,} \emph{Statistical Science,} 1992.}
#'
#' @references \cite{Brooks, S. and Gelman, A.,
#' \dQuote{General Methods for Monitoring Convergence of Iterative Simulations,} \emph{Journal of Computational and Graphical Statistics,} 1998.}
#'
#' @return The value of multivariate PSRF by S. Brooks and A. Gelman.
#'
#' @export
psrf_multivariate <- function(samples) {
.Call(`_volesti_psrf_multivariate`, samples)
}
#' Gelman-Rubin and Brooks-Gelman Potential Scale Reduction Factor (PSRF) for each marginal
#'
#' @param samples A matrix that contans column-wise the sampled points from a geometric random walk.
#' @param method A string to reauest diagnostic: (i) \code{'normal'} for psrf of Gelman-Rubin and (ii) \code{'interval'} for psrf of Brooks-Gelman.
#'
#' @references \cite{Gelman, A. and Rubin, D. B.,
#' \dQuote{Inference from iterative simulation using multiple sequences,} \emph{Statistical Science,} 1992.}
#'
#' @references \cite{Brooks, S. and Gelman, A.,
#' \dQuote{General Methods for Monitoring Convergence of Iterative Simulations,} \emph{Journal of Computational and Graphical Statistics,} 1998.}
#'
#' @return A vector that contains the values of PSRF for each coordinate
#'
#' @export
psrf_univariate <- function(samples, method = NULL) {
.Call(`_volesti_psrf_univariate`, samples, method)
}
#' Raftery and Lewis MCMC diagnostic
#'
#' @param samples A matrix that contans column-wise the sampled points from a geometric random walk.
#' @param q Optional. The quantile of the quantity of interest. The default value is 0.025.
#' @param r Optional. The level of precision desired. The default value is 0.01.
#' @param s Optional. The probability associated with r. The default value is 0.95.
#'
#' @references \cite{Raftery, A. E. and Lewis, S. M.,
#' \dQuote{How many iterations in the Gibbs sampler?,} \emph{Bayesian Statistics 4. Proceedings of the Fourth Valencia International Meeting,} 1992.}
#'
#' @return (i) The number of draws required for burn-in, (ii) the skip parameter for 1st-order Markov chain, (iii) the skip parameter sufficient to get independence chain, (iv) the number of draws required to achieve r precision, (v) the number of draws if the chain is white noise, (vi) the I-statistic from Raftery and Lewis (1992).
#'
#' @export
raftery <- function(samples, q = NULL, r = NULL, s = NULL) {
.Call(`_volesti_raftery`, samples, q, r, s)
}
#' An internal Rccp function for the random rotation of a convex polytope
#'
#' @param P A convex polytope (H-, V-polytope or a zonotope).
#' @param T Optional. A rotation matrix.
#' @param seed Optional. A fixed seed for the random linear map generator.
#'
#' @keywords internal
#'
#' @return A matrix that describes the rotated polytope
rotating <- function(P, T = NULL, seed = NULL) {
.Call(`_volesti_rotating`, P, T, seed)
}
#' Internal rcpp function for the rounding of a convex polytope
#'
#' @param P A convex polytope (H- or V-representation or zonotope).
#' @param method Optional. The method to use for rounding, a) \code{'min_ellipsoid'} for the method based on mimimmum volume enclosing ellipsoid of a uniform sample from P, b) \code{'max_ellipsoid'} for the method based on maximum volume enclosed ellipsoid in P, (c) \code{'isotropy'} for the method based on isotropy. The default method is \code{'min_ellipsoid'} for all the representations.
#' @param seed Optional. A fixed seed for the number generator.
#'
#' @keywords internal
#'
#' @return A numerical matrix that describes the rounded polytope, a numerical matrix of the inverse linear transofmation that is applied on the input polytope, the numerical vector the the input polytope is shifted and the determinant of the matrix of the linear transformation that is applied on the input polytope.
rounding <- function(P, method = NULL, seed = NULL) {
.Call(`_volesti_rounding`, P, method, seed)
}
#' Sample uniformly, normally distributed, or logconcave distributed points from a convex Polytope (H-polytope, V-polytope, zonotope or intersection of two V-polytopes).
#'
#' @param P A convex polytope. It is an object from class (a) Hpolytope or (b) Vpolytope or (c) Zonotope or (d) VpolytopeIntersection.
#' @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_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()}.}
#' \item{\code{BaW_rad} }{ The radius for the ball walk.}
#' \item{\code{L} }{ The maximum length of the billiard trajectory or the radius for the step of dikin, vaidya or john walk.}
#' \item{\code{solver} }{ Specify ODE solver for logconcave sampling. Options are i) leapfrog, ii) euler iii) runge-kutta iv) richardson}
#' \item{\code{step_size }{ Optionally chosen step size for logconcave sampling. Defaults to a theoretical value if not provided.}
#' }
#' @param distribution Optional. A list that declares the target density and some related parameters as follows:
#' \itemize{
#' \item{\code{density} }{ A string: (a) \code{'uniform'} for the uniform distribution or b) \code{'gaussian'} for the multidimensional spherical distribution c) \code{logconcave} with form proportional to exp(-f(x)) where f(x) is L-smooth and m-strongly-convex d) \code{'exponential'} for the exponential distribution. The default target distribution is the uniform distribution.}
#' \item{\code{variance} }{ The variance of the multidimensional spherical gaussian or the exponential distribution. The default value is 1.}
#' \item{\code{mode} }{ A \eqn{d}-dimensional numerical vector that declares the mode of the Gaussian distribution. The default choice is the center of the as that one computed by the function \code{inner_ball()}.}
#' \item{\code{bias} }{ The bias vector for the exponential distribution. The default vector is \eqn{c_1 = 1} and \eqn{c_i = 0} for \eqn{i \neq 1}.}
#' \item{\code{L_} }{ Smoothness constant (for logconcave). }
#' \item{\code{m} }{ Strong-convexity constant (for logconcave). }
#' \item{\code{negative_logprob} }{ Negative log-probability (for logconcave). }
#' \item{\code{negative_logprob_gradient} }{ Negative log-probability gradient (for logconcave). }
#' }
#' @param seed Optional. A fixed seed for the number generator.
#'
#' @references \cite{Robert L. Smith,
#' \dQuote{Efficient Monte Carlo Procedures for Generating Points Uniformly Distributed Over Bounded Regions,} \emph{Operations Research,} 1984.},
#'
#' @references \cite{B.T. Polyak, E.N. Gryazina,
#' \dQuote{Billiard walk - a new sampling algorithm for control and optimization,} \emph{IFAC Proceedings Volumes,} 2014.},
#'
#' @references \cite{Y. Chen, R. Dwivedi, M. J. Wainwright and B. Yu,
#' \dQuote{Fast MCMC Sampling Algorithms on Polytopes,} \emph{Journal of Machine Learning Research,} 2018.}
#'
#' @references \cite{Lee, Yin Tat, Ruoqi Shen, and Kevin Tian,
#' \dQuote{"Logsmooth Gradient Concentration and Tighter Runtimes for Metropolized Hamiltonian Monte Carlo,"} \emph{arXiv preprint arXiv:2002.04121}, 2020.}
#'
#' @references \cite{Shen, Ruoqi, and Yin Tat Lee,
#' \dQuote{"The randomized midpoint method for log-concave sampling.",} \emph{Advances in Neural Information Processing Systems}, 2019.}
#'
#' @references \cite{Augustin Chevallier, Sylvain Pion, Frederic Cazals,
#' \dQuote{"Hamiltonian Monte Carlo with boundary reflections, and application to polytope volume calculations,"} \emph{Research Report preprint hal-01919855}, 2018.}
#'
#' @return A \eqn{d\times n} matrix that contains, column-wise, the sampled points from the convex polytope P.
#' @examples
#' # uniform distribution from the 3d unit cube in H-representation using ball walk
#' P = gen_cube(3, 'H')
#' points = sample_points(P, n = 100, random_walk = list("walk" = "BaW", "walk_length" = 5))
#'
#' # 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$new(A,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
#' P = gen_rand_hpoly(2,20)
#' points = sample_points(P, n = 100, random_walk = list("walk" = "BRDHR"))
#'
#' # For sampling from logconcave densities see the examples directory
#'
#' @export
sample_points <- function(P, n, random_walk = NULL, distribution = NULL, seed = NULL) {
.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$new(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.
#'
#' @param P A convex polytope. It is an object from class a) Hpolytope or b) Vpolytope or c) Zonotope or d) VpolytopeIntersection.
#' @param settings Optional. A list that declares which algorithm, random walk and values of parameters to use, as follows:
#' \itemize{
#' \item{\code{algorithm} }{ A string to set the algorithm to use: a) \code{'CB'} for CB algorithm, b) \code{'SoB'} for SOB algorithm or b) \code{'CG'} for CG algorithm. The defalut algorithm is \code{'CB'}.}
#' \item{\code{error} }{ A numeric value to set the upper bound for the approximation error. The default value is \eqn{1} for SOB algorithm and \eqn{0.1} otherwise.}
#' \item{\code{random_walk} }{ A string that declares the random walk method: a) \code{'CDHR'} for Coordinate Directions Hit-and-Run, b) \code{'RDHR'} for Random Directions Hit-and-Run, c) \code{'BaW'} for Ball Walk, or \code{'BiW'} for Billiard walk. For CB algorithm the default walk is \code{'BiW'}. For CG and SOB algorithms the default walk is \code{'CDHR'} for H-polytopes and \code{'RDHR'} for the other representations.}
#' \item{\code{walk_length} }{ An integer to set the number of the steps for the random walk. The default value is \eqn{\lfloor 10 + d/10\rfloor} for \code{'SOB'} and \eqn{1} otherwise.}
#' \item{\code{win_len} }{ The length of the sliding window for CB or CG algorithm. The default value is \eqn{250} for CB with BiW and \eqn{400+3d^2} for CB and any other random walk and \eqn{500+4d^2} for CG.}
#' \item{\code{hpoly} }{ A boolean parameter to use H-polytopes in MMC of CB algorithm when the input polytope is a zonotope. The default value is \code{TRUE} when the order of the zonotope is \eqn{<5}, otherwise it is \code{FALSE}.}
#' }
#' @param rounding Optional. A string parameter to request a rounding method to be applied in the input polytope before volume computation: a) \code{'min_ellipsoid'}, b) \code{'svd'}, c) \code{'max_ellipsoid'} and d) \code{'none'} for no rounding.
#' @param seed Optional. A fixed seed for the number generator.
#'
#' @references \cite{I.Z.Emiris and V. Fisikopoulos,
#' \dQuote{Practical polytope volume approximation,} \emph{ACM Trans. Math. Soft.,} 2018.},
#' @references \cite{A. Chalkis and I.Z.Emiris and V. Fisikopoulos,
#' \dQuote{Practical Volume Estimation by a New Annealing Schedule for Cooling Convex Bodies,} \emph{CoRR, abs/1905.05494,} 2019.},
#' @references \cite{B. Cousins and S. Vempala, \dQuote{A practical volume algorithm,} \emph{Springer-Verlag Berlin Heidelberg and The Mathematical Programming Society,} 2015.}
#'
#'
#' @return The approximation of the volume of a convex polytope.
#' @examples
#'
#' # calling SOB algorithm for a H-polytope (5d unit simplex)
#' HP = gen_cube(5,'H')
#' pair_vol = volume(HP)
#'
#' # calling CG algorithm for a V-polytope (3d simplex)
#' VP = gen_simplex(3,'V')
#' pair_vol = volume(VP, settings = list("algorithm" = "CG"))
#'
#' # calling CG algorithm for a 2-dimensional zonotope defined as the Minkowski sum of 4 segments
#' Z = gen_rand_zonotope(2, 4)
#' pair_vol = volume(Z, settings = list("random_walk" = "RDHR", "walk_length" = 2))
#'
#' @export
volume <- function(P, settings = NULL, rounding = NULL, seed = NULL) {
.Call(`_volesti_volume`, P, settings, rounding, seed)
}
#' An internal Rccp function for the over-approximation of a zonotope
#'
#' @param Z A zonotope.
#' @param fit_ratio Optional. A boolean parameter to request the computation of the ratio of fitness.
#' @param settings Optional. A list that declares the values of the parameters of CB algorithm.
#' @param seed Optional. A fixed seed for the number generator.
#'
#' @keywords internal
#'
#' @return A List that contains a numerical matrix that describes the PCA approximation as a H-polytope and the ratio of fitness.
zono_approx <- function(Z, fit_ratio = NULL, settings = NULL, seed = NULL) {
.Call(`_volesti_zono_approx`, Z, fit_ratio, settings, seed)
}