diff --git a/.github/workflows/R-CMD-check-ubuntu.yml b/.github/workflows/R-CMD-check-ubuntu.yml index 645396120..9fe81de99 100644 --- a/.github/workflows/R-CMD-check-ubuntu.yml +++ b/.github/workflows/R-CMD-check-ubuntu.yml @@ -22,8 +22,8 @@ jobs: config: - {os: ubuntu-latest, r: 'devel'} - {os: ubuntu-latest, r: 'release'} - - {os: ubuntu-18.04, r: 'devel'} - - {os: ubuntu-18.04, r: 'release'} + - {os: ubuntu-20.04, r: 'devel'} + - {os: ubuntu-20.04, r: 'release'} env: R_REMOTES_NO_ERRORS_FROM_WARNINGS: true diff --git a/.github/workflows/cmake-clang.yml b/.github/workflows/cmake-clang.yml index 5386c3963..caa22fff2 100644 --- a/.github/workflows/cmake-clang.yml +++ b/.github/workflows/cmake-clang.yml @@ -1,5 +1,5 @@ ############################################################################## -# GitHub Actions Workflow for volesti to build tests with GCC +# GitHub Actions Workflow for volesti to build tests with Clang # # Copyright (c) 2020-2022 Vissarion Fisikopoulos # @@ -36,4 +36,4 @@ jobs: cd build; cmake -D CMAKE_CXX_COMPILER=${{ matrix.config.compiler }} -D CMAKE_CXX_FLAGS=-fsanitize=memory -D CMAKE_CXX_FLAGS=-fsanitize=undefined -D CMAKE_CXX_FLAGS=-g -D DISABLE_NLP_ORACLES=ON -D USE_MKL=OFF ../test; make; - ctest --verbose; + ctest --verbose; \ No newline at end of file diff --git a/.github/workflows/cmake-examples.yml b/.github/workflows/cmake-examples.yml new file mode 100644 index 000000000..b7784fa30 --- /dev/null +++ b/.github/workflows/cmake-examples.yml @@ -0,0 +1,40 @@ +############################################################################## +# GitHub Actions Workflow for volesti to build tests with GCC +# +# Copyright (c) 2020-2022 Vissarion Fisikopoulos +# +# Licensed under GNU LGPL.3, see LICENCE file +############################################################################## +name: cmake-examples + +on: [push, pull_request] + +jobs: + build: + name: ${{ matrix.config.os }} - ${{ matrix.config.compiler }} + strategy: + fail-fast: false + matrix: + config: + - {os: ubuntu-22.04, compiler_pkg: clang-11, compiler: clang++-11} + - {os: ubuntu-22.04, compiler_pkg: g++-11, compiler: g++-11} + runs-on: ${{ matrix.config.os }} + steps: + - uses: actions/checkout@v1 + - run: sudo apt-get update || true; + sudo apt-get install ${{ matrix.config.compiler_pkg }} lp-solve libomp-dev libopenblas-dev libarpack2-dev; + - name: Build examples + run: | + cd examples + for dir in */; do + if [ "$dir" != "EnvelopeProblemSOS/" ] && [ "$dir" != "python_utilities/" ]; then + echo + echo "Building examples in $dir ....................." + cd "$dir" + mkdir build && cd build + cmake -DCMAKE_CXX_COMPILER=${{ matrix.config.compiler }} -DUSE_MKL=OFF .. + make + cd ../.. + fi + done + diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 5d2414da9..3af68ab7b 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -2,42 +2,42 @@ :+1::tada: First off, thanks for taking the time to contribute! :tada::+1: -The following is a set of guidelines for contributing to volesti, -which are hosted in the [GeomScale Organization](https://github.com/GeomScale) on GitHub. -These are mostly guidelines, not rules. +The following is a set of guidelines for contributing to volesti, +which are hosted in the [GeomScale Organization](https://github.com/GeomScale) on GitHub. +These are mostly guidelines, not rules. Use your best judgment, and feel free to propose changes to this document in a pull request. ## Table of Contents - * [Prerequisites (how to start)](#prerequisites--how-to-start-) - * [Testing the development branch of volesti (get the tools ready)](#testing-the-development-branch-of-volesti--get-the-tools-ready-) - * [Fork volesti repository (this is your repo now!)](#fork--volesti-repository--this-is-your-repo-now--) - + [Verify if your fork works (optional)](#verify-if-your-fork-works--optional-) - * [Working with volesti (get ready to contribute)](#working-with-volesti--get-ready-to-contribute-) - + [GitFlow workflow](#gitflow-workflow) - + [Create new branch for your work](#create-new-branch-for-your-work) - + [Verify your new branch (optional)](#verify-your-new-branch--optional-) - * [Modify the branch (implement, implement, implement)](#modify-the-branch--implement--implement--implement-) - + [Tests](#tests) - + [Push](#push) - * [Pull request (the joy of sharing)](#pull-request--the-joy-of-sharing-) - * [Review (ok this is not an exam)](#review--ok-this-is-not-an-exam-) - +- [Prerequisites (how to start)](#prerequisites--how-to-start-) +- [Testing the development branch of volesti (get the tools ready)](#testing-the-development-branch-of-volesti--get-the-tools-ready-) +- [Fork volesti repository (this is your repo now!)](#fork--volesti-repository--this-is-your-repo-now--) + - [Verify if your fork works (optional)](#verify-if-your-fork-works--optional-) +- [Working with volesti (get ready to contribute)](#working-with-volesti--get-ready-to-contribute-) + - [GitFlow workflow](#gitflow-workflow) + - [Create new branch for your work](#create-new-branch-for-your-work) + - [Verify your new branch (optional)](#verify-your-new-branch--optional-) +- [Modify the branch (implement, implement, implement)](#modify-the-branch--implement--implement--implement-) + - [Tests](#tests) + - [Push](#push) +- [Pull request (the joy of sharing)](#pull-request--the-joy-of-sharing-) +- [Review (ok this is not an exam)](#review--ok-this-is-not-an-exam-) + ## Prerequisites (how to start) -* git (see [Getting Started with Git](https://help.github.com/en/github/using-git/getting-started-with-git-and-github)) -* a compiler to run tests - gcc, clang, etc. -* configured GitHub account - +- git (see [Getting Started with Git](https://help.github.com/en/github/using-git/getting-started-with-git-and-github)) +- a compiler to run tests - gcc, clang, etc. +- configured GitHub account + Other helpful links: -* http://git-scm.com/documentation -* https://help.github.com/articles/set-up-git -* https://opensource.com/article/18/1/step-step-guide-git +- http://git-scm.com/documentation +- https://help.github.com/articles/set-up-git +- https://opensource.com/article/18/1/step-step-guide-git ## Testing the development branch of volesti (get the tools ready) -Clone the repository, +Clone the repository, git clone git@github.com:GeomScale/volume_approximation.git volesti cd volesti @@ -52,12 +52,14 @@ To compile the `C++` code you have to specify the path to external library `libl # e.g. on linux: cmake -DLP_SOLVE=/usr/lib/lp_solve/liblpsolve55.so .. make -Run the tests, +You can check [here](/docs/getting_started/install.md) to see more installation guide. - ctest -jK +Run the tests, -where `K` is the number of CPU threads. By adding the option `--verbose` to `ctest` you get more information about the tests, -*e.g.* time per test, volume computed and the name of the polytope or convex body. + ctest -jK + +where `K` is the number of CPU threads. By adding the option `--verbose` to `ctest` you get more information about the tests, +_e.g._ time per test, volume computed and the name of the polytope or convex body. ![test_cube](https://user-images.githubusercontent.com/3660366/72348403-0524df00-36e3-11ea-9b6d-288a2bddc22c.png) @@ -65,8 +67,8 @@ If everything works for you, you may move forward. ## Fork volesti repository (this is your repo now!) -You can't work directly in the original volesti repository, therefore you should create your fork of this library. -This way you can modify the code and when the job is done send a pull request to merge your changes with the original +You can't work directly in the original volesti repository, therefore you should create your fork of this library. +This way you can modify the code and when the job is done send a pull request to merge your changes with the original repository. ![fork](https://user-images.githubusercontent.com/3660366/72348562-57fe9680-36e3-11ea-9746-385ff61c752a.png) @@ -90,9 +92,10 @@ clone your repository and checkout develop branch git clone git@github.com:vissarion/volume_approximation.git volesti_fork cd volesti_fork - git checkout develop + git remote add upstream git@github.com:GeomScale/volesti.git + git fetch upstream + git checkout upstream/develop git branch -vv - git pull see commits @@ -105,17 +108,17 @@ For now you should see exactly the same commits as in `volesti` repository. ### GitFlow workflow -Volesit is using the [GitFlow](http://nvie.com/posts/a-successful-git-branching-model/) workflow. -It's because it is very well suited to collaboration and scaling the development team. +Volesit is using the [GitFlow](http://nvie.com/posts/a-successful-git-branching-model/) workflow. +It's because it is very well suited to collaboration and scaling the development team. Each repository using this model should contain two main branches: -* master - release-ready version of the library -* develop - development version of the library - -and could contain various supporting branches for new features and hotfixes. +- master - release-ready version of the library +- develop - development version of the library + +and could contain various supporting branches for new features and hotfixes. -As a contributor you'll most likely be adding new features or fixing bugs in the development version of the library. -This means that for each contribution you should create a new branch originating from the develop branch, +As a contributor you'll most likely be adding new features or fixing bugs in the development version of the library. +This means that for each contribution you should create a new branch originating from the develop branch, modify it and send a pull request in order to merge it, again with the develop branch. ### Create new branch for your work @@ -128,7 +131,7 @@ you should see ![branch -vv](https://user-images.githubusercontent.com/3660366/72348696-a1e77c80-36e3-11ea-93ec-70f5622c0675.png) -Now you should pick a name for your new branch that doesn't already exist. +Now you should pick a name for your new branch that doesn't already exist. The following checks for existing remote branches git branch -a @@ -136,7 +139,7 @@ The following checks for existing remote branches ![List of branches](https://user-images.githubusercontent.com/3660366/72348763-c5aac280-36e3-11ea-8f2c-c66e2c107929.png) Alternatively, you can check them on `GitHub`. -Assume you want to add some new functionality (i.e. a new feature) for example a new sampling algorithm. Then you have +Assume you want to add some new functionality (i.e. a new feature) for example a new sampling algorithm. Then you have to create a new branch e.g. `feature/the_fastest_sampling_algo_ever` Create new local branch @@ -168,12 +171,12 @@ Alternatively, your newly created remote branch is also available on GitHub ## Modify the branch (implement, implement, implement) -Before contributiong to a library by adding a new feature, or a bugfix, or improving documentation, +Before contributiong to a library by adding a new feature, or a bugfix, or improving documentation, it is always wise to interact with the community of developers, for example by opening an issue. ### Tests -Tests are placed in the `test` directory and use the [doctest](https://github.com/onqtam/doctest) library. +Tests are placed in the `test` directory and use the [doctest](https://github.com/onqtam/doctest) library. It is recommended to add new test whenever you contribute a new functionality/feature. Also if your contribution is a bugfix then consider adding this case to the test-suite. @@ -202,16 +205,16 @@ and click the "Create pull request" button. ## Review (ok this is not an exam) -After creating a pull request your code will be reviewed. You can propose one or more reviewers +After creating a pull request your code will be reviewed. You can propose one or more reviewers by clicking on the "Reviewers" button ![reviewer](https://user-images.githubusercontent.com/3660366/72349476-44ecc600-36e5-11ea-81cd-d0938d923529.png) -If there are no objections your changes will be merged. +If there are no objections your changes will be merged. Otherwise you'll see some comments under the pull request and/or under specific lines of your code. -Then you have to make the required changes, commit them and push to your branch. -Those changes will automatically be a part of the same pull request. This procedure will be repeated until the code +Then you have to make the required changes, commit them and push to your branch. +Those changes will automatically be a part of the same pull request. This procedure will be repeated until the code is ready for merging. -If you're curious how it looks like you may see one of the open or closed +If you're curious how it looks like you may see one of the open or closed [pull requests](https://github.com/GeomScale/volume_approximation/pulls). diff --git a/R-proj/R/RcppExports.R b/R-proj/R/RcppExports.R index 5713d11a3..c3f192aee 100644 --- a/R-proj/R/RcppExports.R +++ b/R-proj/R/RcppExports.R @@ -361,11 +361,6 @@ sample_points <- function(P, n, random_walk = NULL, distribution = NULL, seed = .Call(`_volesti_sample_points`, P, n, random_walk, distribution, seed) } -#' @export -sample_spectra <- function(file = NULL, N = NULL, walk_length = NULL) { - .Call(`_volesti_sample_spectra`, file, N, walk_length) -} - #' Write a SDPA format file #' #' Outputs a spectrahedron (the matrices defining a linear matrix inequality) and a vector (the objective function) diff --git a/R-proj/examples/logconcave/nuts_rand_poly.R b/R-proj/examples/logconcave/nuts_rand_poly.R new file mode 100644 index 000000000..00400ba72 --- /dev/null +++ b/R-proj/examples/logconcave/nuts_rand_poly.R @@ -0,0 +1,55 @@ +# 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 using the logconcave sampling methods + +# Import required libraries +library(ggplot2) +library(volesti) + +# Sampling from logconcave density example + +# Helper function +norm_vec <- function(x) sqrt(sum(x^2)) + +# Negative log-probability oracle +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 + +# Create domain of truncation +H <- gen_rand_hpoly(dimension, facets, seed = 15) + +# Rounding +Tr <- rounding(H, seed = 127) + +P <- Hpolytope$new(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)) + +# Sample points +n_samples <- 20000 + +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)) + +# Plot histogram +hist(samples[1,], probability=TRUE, breaks = 100) + +psrfs <- psrf_univariate(samples) +n_ess <- ess(samples) + diff --git a/R-proj/examples/logconcave/simple_hmc_rand_poly.R b/R-proj/examples/logconcave/simple_hmc_rand_poly.R index 9f785aaa8..b352c9a16 100644 --- a/R-proj/examples/logconcave/simple_hmc_rand_poly.R +++ b/R-proj/examples/logconcave/simple_hmc_rand_poly.R @@ -29,10 +29,10 @@ dimension <- 50 facets <- 200 # Create domain of truncation -H <- gen_rand_hpoly(dimension, facets) +H <- gen_rand_hpoly(dimension, facets, seed = 15) # Rounding -Tr <- rounding(H) +Tr <- rounding(H, seed = 127) P <- Hpolytope$new(A = Tr$Mat[1:nrow(Tr$Mat), 2:ncol(Tr$Mat)], b = Tr$Mat[,1]) diff --git a/R-proj/man/sample_points.Rd b/R-proj/man/sample_points.Rd index 2f08d6df9..e969e557d 100644 --- a/R-proj/man/sample_points.Rd +++ b/R-proj/man/sample_points.Rd @@ -13,7 +13,7 @@ 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) xi) \code{'ULD'} for Underdamped Langevin Dynamics using the Randomized Midpoint Method. 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 or \code{'CRHMC'} for logconcave distributions and H-polytopes and sparse_constraint_problems.} +\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()}.} @@ -25,13 +25,14 @@ sample_points(P, n, random_walk = NULL, distribution = NULL, seed = NULL) \item{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. The default target distribution is uniform. c) Logconcave with form proportional to exp(-f(x)) where f(x) is L-smooth and m-strongly-convex. } -\item{\code{variance} }{ The variance of the multidimensional spherical gaussian. The default value is 1.} +\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{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). } +\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). } }} \item{seed}{Optional. A fixed seed for the number generator.} @@ -75,4 +76,7 @@ points = sample_points(P, n = 100, random_walk = list("walk" = "BRDHR")) \cite{Shen, Ruoqi, and Yin Tat Lee, \dQuote{"The randomized midpoint method for log-concave sampling.",} \emph{Advances in Neural Information Processing Systems}, 2019.} + +\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.} } diff --git a/R-proj/src/sample_points.cpp b/R-proj/src/sample_points.cpp index a0b50eec9..1014e74ea 100644 --- a/R-proj/src/sample_points.cpp +++ b/R-proj/src/sample_points.cpp @@ -36,6 +36,7 @@ enum random_walks { brdhr, bcdhr, hmc, + nuts, gaussian_hmc, exponential_hmc, uld, @@ -219,6 +220,26 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo execute_crhmc(P, rng, randPoints, walkL, numpoints, nburns, F, f, h); break; } + case nuts: + + logconcave_sampling < + PointList, + Polytope, + RNGType, + NutsHamiltonianMonteCarloWalk, + NT, + Point, + NegativeGradientFunctor, + NegativeLogprobFunctor, + LeapfrogODESolver < + Point, + NT, + Polytope, + NegativeGradientFunctor + > + >(randPoints, P, rng, walkL, numpoints, StartingPoint, nburns, *F, *f); + + break; case uld: logconcave_sampling < @@ -250,7 +271,17 @@ void sample_from_polytope(Polytope &P, int type, RNGType &rng, PointList &randPo //' @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()}.} @@ -357,7 +388,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, random_walks walk; ode_solvers solver; // Used only for logconcave sampling - NT eta; + NT eta = 1; std::list randPoints; std::pair InnerBall; Point mode(dim); @@ -421,18 +452,24 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, Rcpp::Function negative_logprob = Rcpp::as(distribution)["negative_logprob"]; Rcpp::Function negative_logprob_gradient = Rcpp::as(distribution)["negative_logprob_gradient"]; - NT L_, m, eta; + NT L_ = 1, m = 1; if (Rcpp::as(distribution).containsElementNamed("L_")) { L_ = Rcpp::as(Rcpp::as(distribution)["L_"]); + if (L_ <= NT(0)) { + throw Rcpp::exception("The smoothness constant must be positive"); + } } else { - throw Rcpp::exception("The smoothness constant is absent"); + L_ = -1; } if (Rcpp::as(distribution).containsElementNamed("m")) { m = Rcpp::as(Rcpp::as(distribution)["m"]); + if (m <= NT(0)) { + throw Rcpp::exception("The strong-convexity constant must be positive"); + } } else { - throw Rcpp::exception("The strong-convexity constant is absent"); + m = -1; } if (Rcpp::as(random_walk).containsElementNamed("step_size")) { @@ -456,7 +493,6 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, throw Rcpp::exception("Invalid ODE solver specified. Aborting."); } } else { - Rcpp::warning("Solver set to leapfrog."); solver = leapfrog; } // Create functors @@ -474,8 +510,6 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, functor_defined = false; - NT eta; - if (Rcpp::as(random_walk).containsElementNamed("step_size")) { eta = NT(Rcpp::as(Rcpp::as(random_walk)["step_size"])); if (eta <= NT(0)) { @@ -513,18 +547,12 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, throw Rcpp::exception("Exponential sampling is supported only for H-polytopes"); } walk = exponential_hmc; + } else if (logconcave) { + walk = (type == 5) ? crhmc : nuts; } else if (gaussian) { - if (type == 1) { - walk = cdhr; - } else { - walk = rdhr; - } + walk = (type == 1) ? cdhr : rdhr; } else { - if (type == 1) { - walk = accelarated_billiard; - } else { - walk = billiard; - } + walk = (type == 1) ? accelarated_billiard : billiard; } } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("CDHR")) == 0) { walk = cdhr; @@ -589,7 +617,12 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, } } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("HMC")) == 0) { if (!logconcave) throw Rcpp::exception("HMC is not supported for non first-order sampling"); + if (F->params.L < 0) throw Rcpp::exception("The smoothness constant is absent"); + if (F->params.m < 0) throw Rcpp::exception("The strong-convexity constant is absent"); walk = hmc; + } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("NUTS")) == 0) { + if (!logconcave) throw Rcpp::exception("NUTS is not supported for non first-order sampling"); + walk = nuts; } else if (Rcpp::as(Rcpp::as(random_walk)["walk"]).compare(std::string("ULD")) == 0) { if (!logconcave) throw Rcpp::exception("ULD is not supported for non first-order sampling"); walk = uld; @@ -749,7 +782,7 @@ Rcpp::NumericMatrix sample_points(Rcpp::Nullable P, if (numpoints % 2 == 1 && (walk == brdhr || walk == bcdhr)) numpoints--; MT RetMat(dim, numpoints); unsigned int jj = 0; - + for (typename std::list::iterator rpit = randPoints.begin(); rpit!=randPoints.end(); rpit++, jj++) { if (gaussian) { RetMat.col(jj) = (*rpit).getCoefficients() + mode.getCoefficients(); diff --git a/R-proj/tests/testthat/test_sampling.R b/R-proj/tests/testthat/test_sampling.R index 47b526c83..d051dc70e 100644 --- a/R-proj/tests/testthat/test_sampling.R +++ b/R-proj/tests/testthat/test_sampling.R @@ -87,15 +87,5 @@ for (i in 1:2) { res = runsample(Z, 'zonotope_4_8', distribution) expect_equal(res, 1) }) - - if (Sys.info()["sysname"] != "Windows") - { - test_that("Sampling test", { - P = gen_simplex(10, 'H') - psrf = logconcave_sample(P,distribution,5000,2000) - expect_lte(psrf, 1.2) - }) - } - - + } diff --git a/cran_gen/genCRANpkg.R b/cran_gen/genCRANpkg.R index 45c41b109..abf51d88a 100644 --- a/cran_gen/genCRANpkg.R +++ b/cran_gen/genCRANpkg.R @@ -42,6 +42,8 @@ dir_lp = paste0(path,"/lpSolveAPI/inst/include") h_files = dir(dir_lp, "*.h", ignore.case = TRUE, all.files = TRUE) lp_dist = paste0(path,"/external/LPsolve_src/include") file.copy(file.path(dir_lp, h_files), lp_dist, recursive=TRUE, overwrite=TRUE) +# replace the lpsolve header file that issues a warning in mac's cran test +file.copy(paste0(path,"/external/cmake-files/lpsolve_modified_header_files/lp_rlp.h"), lp_dist, recursive=TRUE, overwrite=TRUE) dir_lp = paste0(path,"/lpSolve/src") h_files = dir(dir_lp, "*.h", ignore.case = TRUE, all.files = TRUE) lp_dist = paste0(path,"/external/LPsolve_src/run_headers") diff --git a/docs/getting_started/install.md b/docs/getting_started/install.md index 57962f8ff..cfdb773e7 100644 --- a/docs/getting_started/install.md +++ b/docs/getting_started/install.md @@ -17,6 +17,21 @@ make ``` For example: `-DLP_SOLVE=/usr/lib/lpsolve/liblpsolve55.so` +In WSL (Windows Subsystem Linux), you can run the following command to install libc6-dev-i386. This will be required for `ieeefp.h` which is used by `qd` library, + + sudo apt-get install libc6-dev-i386 + +Also to install `mkl` related dependencies, run the following, + + wget https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB + sudo apt-key add GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB + sudo sh -c 'echo deb https://apt.repos.intel.com/mkl all main > /etc/apt/sources.list.d/intel-mkl.list' + sudo apt-get update + sudo apt-get install intel-mkl-2020.4-912 + sudo sh -c "echo /opt/intel/mkl/lib/intel64 > /etc/ld.so.conf.d/intel-mkl.conf" + sudo ldconfig + export CPLUS_INCLUDE_PATH="/opt/intel/mkl/include:$CPLUS_INCLUDE_PATH" + You can run the tests by `cmake test` or `ctest -jK` where `K` the number of `CPU` threads. By adding the option `--verbose` to `ctest` you get more information about the tests, *e.g.* time per test, volume computed and the name of the polytope or convex body. ### Development environment from Docker container diff --git a/examples/hpolytope-volume/hpolytopeVolume.cpp b/examples/hpolytope-volume/hpolytopeVolume.cpp index 487ce98f4..b61190de9 100644 --- a/examples/hpolytope-volume/hpolytopeVolume.cpp +++ b/examples/hpolytope-volume/hpolytopeVolume.cpp @@ -29,7 +29,7 @@ typedef typename Kernel::Point Point; typedef BoostRandomNumberGenerator RNGType; typedef HPolytope HPOLYTOPE; -void calculateVolumes(const HPOLYTOPE &HP) { +void calculateVolumes(HPOLYTOPE &HP) { // Setup parameters for calculating volume int walk_len = 10 + HP.dimension()/10; NT e=0.1; diff --git a/examples/logconcave/CMakeLists.txt b/examples/logconcave/CMakeLists.txt index 0448ce3e3..33725d0c9 100644 --- a/examples/logconcave/CMakeLists.txt +++ b/examples/logconcave/CMakeLists.txt @@ -20,7 +20,7 @@ endif(COMMAND cmake_policy) option(DISABLE_NLP_ORACLES "Disable non-linear oracles (used in collocation)" ON) option(BUILTIN_EIGEN "Use eigen from ../external" OFF) -option(USE_MKL "Use MKL library to build eigen" ON) +option(USE_MKL "Use MKL library to build eigen" OFF) if(DISABLE_NLP_ORACLES) add_definitions(-DDISABLE_NLP_ORACLES) @@ -73,16 +73,24 @@ else () include_directories(BEFORE /usr/include/eigen3) endif(BUILTIN_EIGEN) +find_library(BLAS NAMES libblas.so libblas.dylib PATHS /usr/local/Cellar/lapack/3.9.1_1/lib /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu /usr/local/Cellar/openblas/0.3.15_1/lib /usr/lib) + if (USE_MKL) - find_library(BLAS NAME libblas.so PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu) + find_library(BLAS NAMES libblas.so libblas.dylib PATHS /usr/local/Cellar/lapack/3.9.1_1/lib /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu /usr/local/Cellar/openblas/0.3.15_1/lib /usr/lib) + find_library(GFORTRAN NAME libgfortran.dylib PATHS /usr/local/Cellar/gcc/10.2.0_4/lib/gcc/10) + find_library(LAPACK NAME liblapack.dylib PATHS /usr/lib) + find_library(OPENMP NAME libiomp5.dylib PATHS /opt/intel/oneapi/compiler/2021.1.1/mac/compiler/lib) + include_directories (BEFORE ${MKLROOT}/include) - set(PROJECT_LIBS ${BLAS_LIBRARIES}) + set(PROJECT_LIBS ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${GFORTRAN_LIBRARIES}) set(MKL_LINK "-L${MKLROOT}/lib -Wl,-rpath,${MKLROOT}/lib -lmkl_intel_ilp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl") add_definitions(-DEIGEN_USE_MKL_ALL) +else() + set(MKL_LINK "") endif(USE_MKL) # Find lpsolve library -find_library(LP_SOLVE NAMES liblpsolve55.so PATHS /usr/lib/lp_solve) +find_library(LP_SOLVE NAMES liblpsolve55.so liblpsolve55.dylib PATHS /usr/lib/lp_solve /usr/local/lib) if (NOT LP_SOLVE) message(FATAL_ERROR "This program requires the lp_solve library, and will not be compiled.") @@ -145,6 +153,6 @@ else () TARGET_LINK_LIBRARIES(simple_hmc ${LP_SOLVE} ${BLAS} ${MKL_LINK} ) TARGET_LINK_LIBRARIES(exponential_exact_hmc ${LP_SOLVE} ${BLAS} ${MKL_LINK} ) TARGET_LINK_LIBRARIES(gaussian_exact_hmc ${LP_SOLVE} ${BLAS} ${MKL_LINK} ) - + endif() diff --git a/examples/mmcs_method/CMakeLists.txt b/examples/mmcs_method/CMakeLists.txt index 1ee6ce0a5..01b8a040c 100644 --- a/examples/mmcs_method/CMakeLists.txt +++ b/examples/mmcs_method/CMakeLists.txt @@ -105,7 +105,7 @@ else () add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-DBOOST_NO_AUTO_PTR") #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgslcblas") #add_definitions( "-O3 -lgsl -lm -ldl -lgslcblas" ) - + find_package(OpenMP REQUIRED) add_executable (skinny_cube_10_dim skinny_cube_10_dim.cpp) diff --git a/examples/multithread_sampling/CMakeLists.txt b/examples/multithread_sampling/CMakeLists.txt index fe4d5ee75..ded4b4dbb 100644 --- a/examples/multithread_sampling/CMakeLists.txt +++ b/examples/multithread_sampling/CMakeLists.txt @@ -97,7 +97,7 @@ else () endif () - add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11") # enable C++11 standard + #add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11") # enable C++11 standard add_definitions(${CMAKE_CXX_FLAGS} "-O3") # optimization of the compiler #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm") @@ -105,11 +105,11 @@ else () add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-DBOOST_NO_AUTO_PTR") #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgslcblas") #add_definitions( "-O3 -lgsl -lm -ldl -lgslcblas" ) - + find_package(OpenMP REQUIRED) add_executable (sampling_multithread_walks sampling_multithread_walks.cpp) - + TARGET_LINK_LIBRARIES(sampling_multithread_walks ${LP_SOLVE} OpenMP::OpenMP_CXX) diff --git a/examples/multithread_sampling/sampling_multithread_walks.cpp b/examples/multithread_sampling/sampling_multithread_walks.cpp index 191b8e6c5..cb416c72c 100644 --- a/examples/multithread_sampling/sampling_multithread_walks.cpp +++ b/examples/multithread_sampling/sampling_multithread_walks.cpp @@ -48,12 +48,12 @@ MT get_uniform_samples(Polytope &P, std::list randPoints; typedef RandomPointGeneratorMultiThread _RandomPointGeneratorMultiThread; - + _RandomPointGeneratorMultiThread::apply(P, p, N, walk_len, num_threads, randPoints, push_back_policy, rng); unsigned int d = P.dimension(); - + MT samples(d, N); unsigned int jj = 0; @@ -92,7 +92,7 @@ MT get_gaussian_samples(Polytope &P, std::list randPoints; typedef GaussianPointGeneratorMultiThread _GaussianPointGeneratorMultiThread; - + _GaussianPointGeneratorMultiThread::apply(P, p, a_i, N, walk_len, num_threads, randPoints, push_back_policy, rng); @@ -118,11 +118,10 @@ void test_uniform_random_walk(std::string random_walk, unsigned int const& num_t typedef Eigen::Matrix MT; typedef Eigen::Matrix VT; typedef BoostRandomNumberGenerator RNGType; - Hpolytope P; unsigned int d = 10, walk_len = 5, N = 5000; std::cout << "--- Sampling with " + random_walk + " walk from H-cube10" << std::endl; - P = generate_cube(d, false); + Hpolytope P = generate_cube(d, false); RNGType rng(P.dimension()); MT samples = get_uniform_samples(P, rng, walk_len, N, num_threads); @@ -155,7 +154,7 @@ void test_gaussian_random_walk(std::string random_walk, unsigned int const& num_ int main() { - + unsigned int num_threads = 2; test_uniform_random_walk("BRDHR", num_threads); diff --git a/examples/optimization_spectrahedra/CMakeLists.txt b/examples/optimization_spectrahedra/CMakeLists.txt index 356a3e6ac..a35aee4c7 100644 --- a/examples/optimization_spectrahedra/CMakeLists.txt +++ b/examples/optimization_spectrahedra/CMakeLists.txt @@ -78,11 +78,12 @@ IF (OPENBLAS_LIB) IF (GFORTRAN_LIB) message( STATUS "GFORTRAN_LIB found: ${GFORTRAN_LIB}" ) - add_executable (boltzmann_hmc_walk boltzmann_hmc_walk.cpp) - TARGET_LINK_LIBRARIES(boltzmann_hmc_walk ${ARPACK_LIB} ${LAPACK_LIBRARIES} ${GFORTRAN_LIB}) + # These are not compiling, see issue https://github.com/GeomScale/volesti/issues/264 + #add_executable (boltzmann_hmc_walk boltzmann_hmc_walk.cpp) + #TARGET_LINK_LIBRARIES(boltzmann_hmc_walk ${ARPACK_LIB} ${LAPACK_LIBRARIES} ${GFORTRAN_LIB}) - add_executable (solve_sdp solve_sdp.cpp) - TARGET_LINK_LIBRARIES(solve_sdp ${GFORTRAN_LIB} ${LAPACK_LIBRARIES} ${ARPACK_LIB}) + #add_executable (solve_sdp solve_sdp.cpp) + #TARGET_LINK_LIBRARIES(solve_sdp ${GFORTRAN_LIB} ${LAPACK_LIBRARIES} ${ARPACK_LIB}) ELSE() MESSAGE(STATUS "gfortran is required but it could not be found") diff --git a/examples/sampling-hpolytope-with-billiard-walks/sampler.cpp b/examples/sampling-hpolytope-with-billiard-walks/sampler.cpp index 16b6f2367..22dd6d9e5 100644 --- a/examples/sampling-hpolytope-with-billiard-walks/sampler.cpp +++ b/examples/sampling-hpolytope-with-billiard-walks/sampler.cpp @@ -37,7 +37,7 @@ void write_to_file(std::string filename, std::vector const& randPoints) { } -void sample_using_uniform_billiard_walk(HPOLYTOPE const& HP, RNGType& rng, unsigned int walk_len, unsigned int num_points) { +void sample_using_uniform_billiard_walk(HPOLYTOPE& HP, RNGType& rng, unsigned int walk_len, unsigned int num_points) { std::string filename = "uniform_billiard_walk.txt"; typedef RandomPointGenerator Generator; @@ -49,7 +49,7 @@ void sample_using_uniform_billiard_walk(HPOLYTOPE const& HP, RNGType& rng, unsig } -void sample_using_uniform_accelerated_billiard_walk(HPOLYTOPE const& HP, RNGType& rng, unsigned int walk_len, unsigned int num_points) { +void sample_using_uniform_accelerated_billiard_walk(HPOLYTOPE& HP, RNGType& rng, unsigned int walk_len, unsigned int num_points) { std::string filename = "uniform_accelerated_billiard_walk.txt"; typedef RandomPointGenerator Generator; @@ -61,7 +61,7 @@ void sample_using_uniform_accelerated_billiard_walk(HPOLYTOPE const& HP, RNGType } -void sample_using_gaussian_billiard_walk(HPOLYTOPE const& HP, RNGType& rng, unsigned int walk_len, unsigned int num_points) { +void sample_using_gaussian_billiard_walk(HPOLYTOPE& HP, RNGType& rng, unsigned int walk_len, unsigned int num_points) { std::string filename = "gaussian_billiard_walk.txt"; typedef MultivariateGaussianRandomPointGenerator Generator; diff --git a/examples/vpolytope-volume/CMakeLists.txt b/examples/vpolytope-volume/CMakeLists.txt index fdcb3095f..5ffd4784b 100644 --- a/examples/vpolytope-volume/CMakeLists.txt +++ b/examples/vpolytope-volume/CMakeLists.txt @@ -106,7 +106,7 @@ else () add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-ldl") add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-DBOOST_NO_AUTO_PTR") - add_executable (vpolytopeVolume vpolytopeVolume.cpp) - TARGET_LINK_LIBRARIES(vpolytopeVolume ${LP_SOLVE}) + add_executable (vpolytopevolume vpolytopevolume.cpp) + TARGET_LINK_LIBRARIES(vpolytopevolume ${LP_SOLVE}) endif() diff --git a/examples/vpolytope-volume/vpolytopevolume.cpp b/examples/vpolytope-volume/vpolytopevolume.cpp index e814b5875..007964232 100644 --- a/examples/vpolytope-volume/vpolytopevolume.cpp +++ b/examples/vpolytope-volume/vpolytopevolume.cpp @@ -28,7 +28,7 @@ typedef typename Kernel::Point Point; typedef BoostRandomNumberGenerator RNGType; typedef VPolytope VPOLYTOPE; -void calculateVolumes(const VPOLYTOPE &VP) { +void calculateVolumes(VPOLYTOPE &VP) { // Setup parameters for calculating volume int walk_len = 10 + VP.dimension()/10; NT e=0.1; diff --git a/external/cmake-files/LPSolve.cmake b/external/cmake-files/LPSolve.cmake index b3f3cf98b..c5a173837 100644 --- a/external/cmake-files/LPSolve.cmake +++ b/external/cmake-files/LPSolve.cmake @@ -26,7 +26,7 @@ function(GetLPSolve) message(STATUS "lp_solve library found: ${LP_SOLVE_DIR}") endif() - + include_directories(${LP_SOLVE_DIR}) endfunction() diff --git a/external/cmake-files/lpsolve_modified_header_files/lp_rlp.h b/external/cmake-files/lpsolve_modified_header_files/lp_rlp.h new file mode 100644 index 000000000..058b45232 --- /dev/null +++ b/external/cmake-files/lpsolve_modified_header_files/lp_rlp.h @@ -0,0 +1,2460 @@ +// Copyright(c) 2016-2018 Kjell Konis . +// Version: 5.5.2.0-17 +// Description: The lpSolveAPI package provides an R interface to 'lp_solve', +// a Mixed Integer Linear Programming (MILP) solver with support for pure +// linear, (mixed) integer/binary, semi-continuous and special ordered sets +// (SOS) models. +// License: LGPL-2 +// Repository: CRAN + + +#define YY_INT_ALIGNED short int + +/* A lexical scanner generated by flex */ + +#define FLEX_SCANNER +#define YY_FLEX_MAJOR_VERSION 2 +#define YY_FLEX_MINOR_VERSION 5 +#define YY_FLEX_SUBMINOR_VERSION 35 +#if YY_FLEX_SUBMINOR_VERSION > 0 +#define FLEX_BETA +#endif + +/* First, we deal with platform-specific or compiler-specific issues. */ + +/* begin standard C headers. */ +#include +#include +#include +#include + +/* end standard C headers. */ + +/* flex integer type definitions */ + +#ifndef FLEXINT_H +#define FLEXINT_H + +/* C99 systems have . Non-C99 systems may or may not. */ + +#if defined (__STDC_VERSION__) && __STDC_VERSION__ >= 199901L + +/* C99 says to define __STDC_LIMIT_MACROS before including stdint.h, + * if you want the limit (max/min) macros for int types. + */ +#ifndef __STDC_LIMIT_MACROS +#define __STDC_LIMIT_MACROS 1 +#endif + +#include +typedef int8_t flex_int8_t; +typedef uint8_t flex_uint8_t; +typedef int16_t flex_int16_t; +typedef uint16_t flex_uint16_t; +typedef int32_t flex_int32_t; +typedef uint32_t flex_uint32_t; +#else +typedef signed char flex_int8_t; +typedef short int flex_int16_t; +typedef int flex_int32_t; +typedef unsigned char flex_uint8_t; +typedef unsigned short int flex_uint16_t; +typedef unsigned int flex_uint32_t; +#endif /* ! C99 */ + +/* Limits of integral types. */ +#ifndef INT8_MIN +#define INT8_MIN (-128) +#endif +#ifndef INT16_MIN +#define INT16_MIN (-32767-1) +#endif +#ifndef INT32_MIN +#define INT32_MIN (-2147483647-1) +#endif +#ifndef INT8_MAX +#define INT8_MAX (127) +#endif +#ifndef INT16_MAX +#define INT16_MAX (32767) +#endif +#ifndef INT32_MAX +#define INT32_MAX (2147483647) +#endif +#ifndef UINT8_MAX +#define UINT8_MAX (255U) +#endif +#ifndef UINT16_MAX +#define UINT16_MAX (65535U) +#endif +#ifndef UINT32_MAX +#define UINT32_MAX (4294967295U) +#endif + +#endif /* ! FLEXINT_H */ + +#ifdef __cplusplus + +/* The "const" storage-class-modifier is valid. */ +#define YY_USE_CONST + +#else /* ! __cplusplus */ + +/* C99 requires __STDC__ to be defined as 1. */ +#if defined (__STDC__) + +#define YY_USE_CONST + +#endif /* defined (__STDC__) */ +#endif /* ! __cplusplus */ + +#ifdef YY_USE_CONST +#define lp_yyconst const +#else +#define lp_yyconst +#endif + +/* Returned upon end-of-file. */ +#define YY_NULL 0 + +/* Promotes a possibly negative, possibly signed char to an unsigned + * integer for use as an array index. If the signed char is negative, + * we want to instead treat it as an 8-bit unsigned char, hence the + * double cast. + */ +#define YY_SC_TO_UI(c) ((unsigned int) (unsigned char) c) + +/* An opaque pointer. */ +#ifndef YY_TYPEDEF_YY_SCANNER_T +#define YY_TYPEDEF_YY_SCANNER_T +typedef void* lp_yyscan_t; +#endif + +/* For convenience, these vars (plus the bison vars far below) + are macros in the reentrant scanner. */ +#define lp_yyin lp_yyg->lp_yyin_r +#define lp_yyout lp_yyg->lp_yyout_r +#define lp_yyextra lp_yyg->lp_yyextra_r +#define lp_yyleng lp_yyg->lp_yyleng_r +#define lp_yytext lp_yyg->lp_yytext_r +#define lp_yylineno (YY_CURRENT_BUFFER_LVALUE->lp_yy_bs_lineno) +#define lp_yycolumn (YY_CURRENT_BUFFER_LVALUE->lp_yy_bs_column) +#define lp_yy_flex_debug lp_yyg->lp_yy_flex_debug_r + +/* Enter a start condition. This macro really ought to take a parameter, + * but we do it the disgusting crufty way forced on us by the ()-less + * definition of BEGIN. + */ +#define BEGIN lp_yyg->lp_yy_start = 1 + 2 * + +/* Translate the current start state into a value that can be later handed + * to BEGIN to return to the state. The YYSTATE alias is for lex + * compatibility. + */ +#define YY_START ((lp_yyg->lp_yy_start - 1) / 2) +#define YYSTATE YY_START + +/* Action number for EOF rule of a given start state. */ +#define YY_STATE_EOF(state) (YY_END_OF_BUFFER + state + 1) + +/* Special action meaning "start processing a new file". */ +#define YY_NEW_FILE lp_yyrestart(lp_yyin ,lp_yyscanner ) + +#define YY_END_OF_BUFFER_CHAR 0 + +/* Size of default input buffer. */ +#ifndef YY_BUF_SIZE +#define YY_BUF_SIZE 16384 +#endif + +/* The state buf must be large enough to hold one state per character in the main buffer. + */ +#define YY_STATE_BUF_SIZE ((YY_BUF_SIZE + 2) * sizeof(lp_yy_state_type)) + +#ifndef YY_TYPEDEF_YY_BUFFER_STATE +#define YY_TYPEDEF_YY_BUFFER_STATE +typedef struct lp_yy_buffer_state *YY_BUFFER_STATE; +#endif + +#define EOB_ACT_CONTINUE_SCAN 0 +#define EOB_ACT_END_OF_FILE 1 +#define EOB_ACT_LAST_MATCH 2 + + /* Note: We specifically omit the test for lp_yy_rule_can_match_eol because it requires + * access to the local variable lp_yy_act. Since lp_yyless() is a macro, it would break + * existing scanners that call lp_yyless() from OUTSIDE lp_yylex. + * One obvious solution it to make lp_yy_act a global. I tried that, and saw + * a 5% performance hit in a non-lp_yylineno scanner, because lp_yy_act is + * normally declared as a register variable-- so it is not worth it. + */ + #define YY_LESS_LINENO(n) \ + do { \ + int lp_yyl;\ + for ( lp_yyl = n; lp_yyl < lp_yyleng; ++lp_yyl )\ + if ( lp_yytext[lp_yyl] == '\n' )\ + --lp_yylineno;\ + }while(0) + +/* Return all but the first "n" matched characters back to the input stream. */ +#define lp_yyless(n) \ + do \ + { \ + /* Undo effects of setting up lp_yytext. */ \ + int lp_yyless_macro_arg = (n); \ + YY_LESS_LINENO(lp_yyless_macro_arg);\ + *lp_yy_cp = lp_yyg->lp_yy_hold_char; \ + YY_RESTORE_YY_MORE_OFFSET \ + lp_yyg->lp_yy_c_buf_p = lp_yy_cp = lp_yy_bp + lp_yyless_macro_arg - YY_MORE_ADJ; \ + YY_DO_BEFORE_ACTION; /* set up lp_yytext again */ \ + } \ + while ( 0 ) + +#define unput(c) lp_yyunput( c, lp_yyg->lp_yytext_ptr , lp_yyscanner ) + +#ifndef YY_TYPEDEF_YY_SIZE_T +#define YY_TYPEDEF_YY_SIZE_T +typedef size_t lp_yy_size_t; +#endif + +#ifndef YY_STRUCT_YY_BUFFER_STATE +#define YY_STRUCT_YY_BUFFER_STATE +struct lp_yy_buffer_state + { + FILE *lp_yy_input_file; + + char *lp_yy_ch_buf; /* input buffer */ + char *lp_yy_buf_pos; /* current position in input buffer */ + + /* Size of input buffer in bytes, not including room for EOB + * characters. + */ + lp_yy_size_t lp_yy_buf_size; + + /* Number of characters read into lp_yy_ch_buf, not including EOB + * characters. + */ + int lp_yy_n_chars; + + /* Whether we "own" the buffer - i.e., we know we created it, + * and can realloc() it to grow it, and should free() it to + * delete it. + */ + int lp_yy_is_our_buffer; + + /* Whether this is an "interactive" input source; if so, and + * if we're using stdio for input, then we want to use getc() + * instead of fread(), to make sure we stop fetching input after + * each newline. + */ + int lp_yy_is_interactive; + + /* Whether we're considered to be at the beginning of a line. + * If so, '^' rules will be active on the next match, otherwise + * not. + */ + int lp_yy_at_bol; + + int lp_yy_bs_lineno; /**< The line count. */ + int lp_yy_bs_column; /**< The column count. */ + + /* Whether to try to fill the input buffer when we reach the + * end of it. + */ + int lp_yy_fill_buffer; + + int lp_yy_buffer_status; + +#define YY_BUFFER_NEW 0 +#define YY_BUFFER_NORMAL 1 + /* When an EOF's been seen but there's still some text to process + * then we mark the buffer as YY_EOF_PENDING, to indicate that we + * shouldn't try reading from the input source any more. We might + * still have a bunch of tokens to match, though, because of + * possible backing-up. + * + * When we actually see the EOF, we change the status to "new" + * (via lp_yyrestart()), so that the user can continue scanning by + * just pointing lp_yyin at a new input file. + */ +#define YY_BUFFER_EOF_PENDING 2 + + }; +#endif /* !YY_STRUCT_YY_BUFFER_STATE */ + +/* We provide macros for accessing buffer states in case in the + * future we want to put the buffer states in a more general + * "scanner state". + * + * Returns the top of the stack, or NULL. + */ +#define YY_CURRENT_BUFFER ( lp_yyg->lp_yy_buffer_stack \ + ? lp_yyg->lp_yy_buffer_stack[lp_yyg->lp_yy_buffer_stack_top] \ + : NULL) + +/* Same as previous macro, but useful when we know that the buffer stack is not + * NULL or when we need an lvalue. For internal use only. + */ +#define YY_CURRENT_BUFFER_LVALUE lp_yyg->lp_yy_buffer_stack[lp_yyg->lp_yy_buffer_stack_top] + +void lp_yyrestart (FILE *input_file ,lp_yyscan_t lp_yyscanner ); +void lp_yy_switch_to_buffer (YY_BUFFER_STATE new_buffer ,lp_yyscan_t lp_yyscanner ); +YY_BUFFER_STATE lp_yy_create_buffer (FILE *file,int size ,lp_yyscan_t lp_yyscanner ); +void lp_yy_delete_buffer (YY_BUFFER_STATE b ,lp_yyscan_t lp_yyscanner ); +void lp_yy_flush_buffer (YY_BUFFER_STATE b ,lp_yyscan_t lp_yyscanner ); +void lp_yypush_buffer_state (YY_BUFFER_STATE new_buffer ,lp_yyscan_t lp_yyscanner ); +void lp_yypop_buffer_state (lp_yyscan_t lp_yyscanner ); + +static void lp_yyensure_buffer_stack (lp_yyscan_t lp_yyscanner ); +static void lp_yy_load_buffer_state (lp_yyscan_t lp_yyscanner ); +static void lp_yy_init_buffer (YY_BUFFER_STATE b,FILE *file ,lp_yyscan_t lp_yyscanner ); + +#define YY_FLUSH_BUFFER lp_yy_flush_buffer(YY_CURRENT_BUFFER ,lp_yyscanner) + +YY_BUFFER_STATE lp_yy_scan_buffer (char *base,lp_yy_size_t size ,lp_yyscan_t lp_yyscanner ); +YY_BUFFER_STATE lp_yy_scan_string (lp_yyconst char *lp_yy_str ,lp_yyscan_t lp_yyscanner ); +YY_BUFFER_STATE lp_yy_scan_bytes (lp_yyconst char *bytes,int len ,lp_yyscan_t lp_yyscanner ); + +void *lp_yyalloc (lp_yy_size_t ,lp_yyscan_t lp_yyscanner ); +void *lp_yyrealloc (void *,lp_yy_size_t ,lp_yyscan_t lp_yyscanner ); +void lp_yyfree (void * ,lp_yyscan_t lp_yyscanner ); + +#define lp_yy_new_buffer lp_yy_create_buffer + +#define lp_yy_set_interactive(is_interactive) \ + { \ + if ( ! YY_CURRENT_BUFFER ){ \ + lp_yyensure_buffer_stack (lp_yyscanner); \ + YY_CURRENT_BUFFER_LVALUE = \ + lp_yy_create_buffer(lp_yyin,YY_BUF_SIZE ,lp_yyscanner); \ + } \ + YY_CURRENT_BUFFER_LVALUE->lp_yy_is_interactive = is_interactive; \ + } + +#define lp_yy_set_bol(at_bol) \ + { \ + if ( ! YY_CURRENT_BUFFER ){\ + lp_yyensure_buffer_stack (lp_yyscanner); \ + YY_CURRENT_BUFFER_LVALUE = \ + lp_yy_create_buffer(lp_yyin,YY_BUF_SIZE ,lp_yyscanner); \ + } \ + YY_CURRENT_BUFFER_LVALUE->lp_yy_at_bol = at_bol; \ + } + +#define YY_AT_BOL() (YY_CURRENT_BUFFER_LVALUE->lp_yy_at_bol) + +/* Begin user sect3 */ + +#define lp_yywrap(n) 1 +#define YY_SKIP_YYWRAP + +typedef unsigned char YY_CHAR; + +typedef int lp_yy_state_type; + +#define lp_yytext_ptr lp_yytext_r + +static lp_yy_state_type lp_yy_get_previous_state (lp_yyscan_t lp_yyscanner ); +static lp_yy_state_type lp_yy_try_NUL_trans (lp_yy_state_type current_state ,lp_yyscan_t lp_yyscanner); +static int lp_yy_get_next_buffer (lp_yyscan_t lp_yyscanner ); +static void lp_yy_fatal_error (lp_yyconst char msg[] ,lp_yyscan_t lp_yyscanner ); + +/* Done after the current pattern has been matched and before the + * corresponding action - sets up lp_yytext. + */ +#define YY_DO_BEFORE_ACTION \ + lp_yyg->lp_yytext_ptr = lp_yy_bp; \ + lp_yyleng = (size_t) (lp_yy_cp - lp_yy_bp); \ + lp_yyg->lp_yy_hold_char = *lp_yy_cp; \ + *lp_yy_cp = '\0'; \ + lp_yyg->lp_yy_c_buf_p = lp_yy_cp; + +#define YY_NUM_RULES 33 +#define YY_END_OF_BUFFER 34 +/* This struct is not used in this scanner, + but its presence is necessary. */ +struct lp_yy_trans_info + { + flex_int32_t lp_yy_verify; + flex_int32_t lp_yy_nxt; + }; +static lp_yyconst flex_int16_t lp_yy_accept[144] = + { 0, + 0, 0, 0, 0, 0, 0, 34, 32, 10, 10, + 27, 17, 11, 32, 32, 14, 26, 31, 29, 28, + 30, 25, 25, 10, 25, 25, 25, 25, 3, 4, + 3, 3, 9, 7, 8, 10, 17, 17, 0, 15, + 1, 6, 15, 14, 0, 29, 30, 0, 25, 24, + 0, 25, 25, 10, 0, 0, 0, 0, 25, 25, + 25, 25, 25, 2, 0, 15, 0, 15, 22, 0, + 25, 25, 0, 0, 0, 0, 0, 19, 25, 18, + 20, 25, 25, 21, 0, 25, 0, 13, 25, 0, + 12, 25, 19, 0, 18, 20, 21, 25, 23, 25, + + 20, 21, 21, 16, 16, 0, 25, 25, 0, 23, + 0, 21, 25, 25, 0, 0, 25, 25, 0, 0, + 19, 25, 0, 0, 25, 25, 19, 0, 18, 0, + 0, 25, 25, 18, 0, 0, 0, 0, 0, 0, + 0, 0, 0 + } ; + +static lp_yyconst flex_int32_t lp_yy_ec[256] = + { 0, + 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, + 1, 1, 4, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 2, 1, 1, 5, 6, 5, 5, 5, 1, + 1, 7, 8, 9, 10, 11, 12, 13, 14, 14, + 13, 13, 13, 13, 13, 13, 13, 15, 16, 17, + 18, 19, 1, 5, 20, 21, 22, 23, 24, 25, + 26, 23, 27, 23, 23, 23, 28, 29, 30, 23, + 23, 31, 32, 33, 34, 23, 23, 35, 36, 37, + 5, 1, 5, 5, 5, 1, 20, 21, 22, 23, + + 24, 25, 26, 23, 27, 23, 23, 23, 28, 29, + 30, 23, 23, 31, 32, 33, 34, 23, 23, 35, + 36, 37, 5, 1, 5, 5, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1 + } ; + +static lp_yyconst flex_int32_t lp_yy_meta[38] = + { 0, + 1, 2, 3, 3, 4, 5, 6, 3, 6, 3, + 5, 5, 5, 5, 7, 6, 7, 6, 6, 4, + 4, 4, 4, 4, 4, 4, 8, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4 + } ; + +static lp_yyconst flex_int16_t lp_yy_base[150] = + { 0, + 0, 36, 36, 38, 43, 45, 366, 388, 48, 62, + 388, 338, 388, 40, 48, 60, 388, 388, 346, 388, + 326, 60, 65, 91, 81, 74, 85, 102, 388, 388, + 388, 330, 388, 388, 388, 125, 313, 134, 308, 96, + 388, 388, 117, 132, 139, 388, 388, 88, 146, 320, + 0, 149, 152, 0, 307, 301, 294, 83, 153, 156, + 157, 189, 160, 388, 286, 126, 65, 108, 388, 289, + 181, 185, 272, 273, 250, 249, 220, 199, 203, 208, + 192, 211, 219, 227, 243, 109, 163, 225, 202, 174, + 224, 215, 213, 207, 191, 388, 189, 227, 228, 231, + + 244, 240, 253, 276, 388, 170, 260, 262, 166, 388, + 169, 179, 263, 241, 166, 159, 270, 272, 149, 155, + 284, 288, 130, 124, 296, 303, 388, 103, 300, 96, + 45, 324, 328, 388, 82, 311, 79, 68, 54, 56, + 25, 12, 388, 345, 353, 360, 367, 372, 379 + } ; + +static lp_yyconst flex_int16_t lp_yy_def[150] = + { 0, + 143, 1, 144, 144, 145, 145, 143, 143, 143, 143, + 143, 146, 143, 143, 143, 143, 143, 143, 143, 143, + 143, 147, 147, 143, 147, 147, 147, 147, 143, 143, + 143, 143, 143, 143, 143, 143, 146, 143, 143, 143, + 143, 143, 143, 143, 143, 143, 143, 143, 147, 143, + 148, 147, 147, 24, 143, 143, 143, 143, 147, 147, + 147, 147, 147, 143, 143, 143, 143, 143, 143, 148, + 147, 147, 143, 143, 143, 143, 143, 147, 147, 147, + 147, 147, 147, 147, 149, 143, 143, 143, 62, 143, + 143, 62, 143, 143, 143, 143, 143, 62, 62, 62, + + 62, 62, 62, 143, 143, 143, 62, 62, 143, 143, + 143, 143, 62, 62, 143, 143, 62, 62, 143, 143, + 62, 62, 143, 143, 62, 62, 143, 143, 62, 143, + 143, 147, 147, 143, 143, 149, 143, 143, 143, 143, + 143, 143, 0, 143, 143, 143, 143, 143, 143 + } ; + +static lp_yyconst flex_int16_t lp_yy_nxt[426] = + { 0, + 8, 9, 10, 9, 8, 8, 11, 12, 13, 12, + 14, 15, 16, 16, 17, 18, 19, 20, 21, 22, + 22, 22, 22, 22, 22, 22, 22, 23, 22, 22, + 22, 22, 22, 22, 22, 22, 22, 24, 30, 31, + 30, 31, 32, 96, 32, 34, 35, 34, 35, 36, + 36, 36, 40, 40, 41, 37, 25, 37, 142, 42, + 26, 48, 27, 36, 36, 36, 48, 28, 136, 37, + 43, 37, 44, 44, 50, 48, 51, 68, 68, 50, + 136, 51, 48, 45, 52, 141, 48, 140, 50, 48, + 51, 53, 54, 36, 36, 50, 139, 51, 37, 50, + + 37, 51, 50, 48, 60, 138, 76, 59, 40, 40, + 48, 55, 77, 61, 137, 56, 50, 57, 51, 45, + 68, 68, 58, 50, 135, 62, 36, 36, 36, 66, + 66, 63, 37, 134, 37, 38, 38, 38, 66, 66, + 45, 38, 43, 38, 44, 44, 67, 48, 67, 45, + 48, 68, 68, 48, 48, 45, 131, 48, 48, 130, + 50, 48, 51, 50, 87, 51, 50, 50, 51, 51, + 50, 50, 51, 51, 50, 90, 51, 88, 128, 79, + 72, 78, 87, 71, 127, 124, 90, 123, 91, 80, + 48, 84, 112, 48, 120, 88, 119, 51, 116, 91, + + 48, 51, 112, 50, 48, 51, 50, 89, 51, 48, + 81, 92, 48, 50, 111, 51, 82, 50, 98, 51, + 48, 83, 50, 49, 51, 50, 99, 51, 48, 107, + 110, 100, 109, 50, 49, 51, 49, 101, 69, 69, + 103, 50, 108, 51, 104, 104, 104, 49, 49, 49, + 102, 97, 49, 115, 49, 49, 114, 113, 49, 49, + 49, 49, 49, 49, 122, 49, 103, 49, 49, 106, + 96, 49, 49, 49, 49, 81, 49, 104, 104, 104, + 49, 49, 95, 49, 49, 49, 117, 49, 118, 49, + 49, 49, 49, 49, 49, 49, 94, 49, 121, 49, + + 93, 125, 49, 126, 49, 49, 125, 86, 126, 49, + 85, 49, 104, 104, 104, 49, 49, 49, 129, 132, + 49, 49, 75, 49, 49, 87, 133, 49, 49, 90, + 49, 74, 49, 73, 69, 49, 65, 143, 88, 39, + 51, 64, 91, 47, 51, 29, 29, 29, 29, 29, + 29, 29, 29, 33, 33, 33, 33, 33, 33, 33, + 33, 38, 38, 46, 39, 143, 143, 38, 49, 143, + 49, 49, 143, 49, 49, 70, 70, 143, 143, 70, + 105, 105, 143, 105, 105, 105, 105, 7, 143, 143, + 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, + + 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, + 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, + 143, 143, 143, 143, 143 + } ; + +static lp_yyconst flex_int16_t lp_yy_chk[426] = + { 0, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 2, 3, 3, + 4, 4, 3, 142, 4, 5, 5, 6, 6, 9, + 9, 9, 14, 14, 15, 9, 2, 9, 141, 15, + 2, 22, 2, 10, 10, 10, 23, 2, 131, 10, + 16, 10, 16, 16, 22, 26, 22, 67, 67, 23, + 131, 23, 25, 16, 23, 140, 27, 139, 26, 48, + 26, 23, 24, 24, 24, 25, 138, 25, 24, 27, + + 24, 27, 48, 28, 26, 137, 58, 25, 40, 40, + 86, 24, 58, 27, 135, 24, 28, 24, 28, 40, + 68, 68, 24, 86, 130, 28, 36, 36, 36, 43, + 43, 28, 36, 128, 36, 38, 38, 38, 66, 66, + 43, 38, 44, 38, 44, 44, 45, 49, 45, 66, + 52, 45, 45, 53, 59, 44, 124, 60, 61, 123, + 49, 63, 49, 52, 87, 52, 53, 59, 53, 59, + 60, 61, 60, 61, 63, 90, 63, 87, 120, 60, + 53, 59, 71, 52, 119, 116, 72, 115, 90, 61, + 62, 63, 112, 81, 111, 71, 109, 71, 106, 72, + + 78, 72, 97, 62, 79, 62, 81, 71, 81, 80, + 62, 72, 82, 78, 95, 78, 62, 79, 78, 79, + 83, 62, 80, 89, 80, 82, 79, 82, 84, 89, + 94, 80, 93, 83, 89, 83, 92, 82, 91, 88, + 84, 84, 92, 84, 85, 85, 85, 92, 98, 99, + 83, 77, 100, 101, 98, 99, 100, 98, 100, 98, + 99, 102, 114, 100, 114, 101, 103, 102, 114, 85, + 76, 101, 102, 114, 103, 101, 101, 104, 104, 104, + 103, 107, 75, 108, 113, 103, 107, 107, 108, 108, + 113, 117, 107, 118, 108, 113, 74, 117, 113, 118, + + 73, 117, 117, 118, 118, 121, 117, 70, 118, 122, + 65, 121, 136, 136, 136, 122, 121, 125, 122, 125, + 122, 129, 57, 125, 126, 132, 126, 129, 125, 133, + 126, 56, 129, 55, 50, 126, 39, 136, 132, 37, + 132, 32, 133, 21, 133, 144, 144, 144, 144, 144, + 144, 144, 144, 145, 145, 145, 145, 145, 145, 145, + 145, 146, 146, 19, 12, 7, 0, 146, 147, 0, + 147, 147, 0, 147, 147, 148, 148, 0, 0, 148, + 149, 149, 0, 149, 149, 149, 149, 143, 143, 143, + 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, + + 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, + 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, + 143, 143, 143, 143, 143 + } ; + +/* Table of booleans, true if rule could match eol. */ +static lp_yyconst flex_int32_t lp_yy_rule_can_match_eol[34] = + { 0, +0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, }; + +/* The intent behind this definition is that it'll catch + * any uses of REJECT which flex missed. + */ +#define REJECT reject_used_but_not_detected +#define lp_yymore() lp_yymore_used_but_not_detected +#define YY_MORE_ADJ 0 +#define YY_RESTORE_YY_MORE_OFFSET +/* + made reentrant with help of + http://www.usualcoding.eu/post/2007/09/03/Building-a-reentrant-parser-in-C-with-Flex/Bison +*/ +/* + Note that a minimum version of flex is needed to be able to compile this. + Older version don't know the reentrant code. + Version 2.5.4 is not enough. Probably at least v2.5.31 is needed. Tested with v2.5.35 +*/ +/* +** We want the scanner to be reentrant, therefore generate no global variables. +** That what the 'reentrant' option is for. +** 'bison-bridge' is used to create a bison compatible scanner and share lp_yylval +*/ + +#define INITIAL 0 +#define COMMENT 1 +#define LINECOMMENT 2 + +#ifndef YY_NO_UNISTD_H +/* Special case for "unistd.h", since it is non-ANSI. We include it way + * down here because we want the user's section 1 to have been scanned first. + * The user has a chance to override it with an option. + */ +#include +#endif + +#ifndef YY_EXTRA_TYPE +#define YY_EXTRA_TYPE void * +#endif + +/* Holds the entire state of the reentrant scanner. */ +struct lp_yyguts_t + { + + /* User-defined. Not touched by flex. */ + YY_EXTRA_TYPE lp_yyextra_r; + + /* The rest are the same as the globals declared in the non-reentrant scanner. */ + FILE *lp_yyin_r, *lp_yyout_r; + size_t lp_yy_buffer_stack_top; /**< index of top of stack. */ + size_t lp_yy_buffer_stack_max; /**< capacity of stack. */ + YY_BUFFER_STATE * lp_yy_buffer_stack; /**< Stack as an array. */ + char lp_yy_hold_char; + int lp_yy_n_chars; + int lp_yyleng_r; + char *lp_yy_c_buf_p; + int lp_yy_init; + int lp_yy_start; + int lp_yy_did_buffer_switch_on_eof; + int lp_yy_start_stack_ptr; + int lp_yy_start_stack_depth; + int *lp_yy_start_stack; + lp_yy_state_type lp_yy_last_accepting_state; + char* lp_yy_last_accepting_cpos; + + int lp_yylineno_r; + int lp_yy_flex_debug_r; + + char *lp_yytext_r; + int lp_yy_more_flag; + int lp_yy_more_len; + + YYSTYPE * lp_yylval_r; + + }; /* end struct lp_yyguts_t */ + +static int lp_yy_init_globals (lp_yyscan_t lp_yyscanner ); + + /* This must go here because YYSTYPE and YYLTYPE are included + * from bison output in section 1.*/ + # define lp_yylval lp_yyg->lp_yylval_r + +int lp_yylex_init (lp_yyscan_t* scanner); + +int lp_yylex_init_extra (YY_EXTRA_TYPE user_defined,lp_yyscan_t* scanner); + +/* Accessor methods to globals. + These are made visible to non-reentrant scanners for convenience. */ + +int lp_yylex_destroy (lp_yyscan_t lp_yyscanner ); + +int lp_yyget_debug (lp_yyscan_t lp_yyscanner ); + +void lp_yyset_debug (int debug_flag ,lp_yyscan_t lp_yyscanner ); + +YY_EXTRA_TYPE lp_yyget_extra (lp_yyscan_t lp_yyscanner ); + +void lp_yyset_extra (YY_EXTRA_TYPE user_defined ,lp_yyscan_t lp_yyscanner ); + +FILE *lp_yyget_in (lp_yyscan_t lp_yyscanner ); + +void lp_yyset_in (FILE * in_str ,lp_yyscan_t lp_yyscanner ); + +FILE *lp_yyget_out (lp_yyscan_t lp_yyscanner ); + +void lp_yyset_out (FILE * out_str ,lp_yyscan_t lp_yyscanner ); + +int lp_yyget_leng (lp_yyscan_t lp_yyscanner ); + +char *lp_yyget_text (lp_yyscan_t lp_yyscanner ); + +int lp_yyget_lineno (lp_yyscan_t lp_yyscanner ); + +void lp_yyset_lineno (int line_number ,lp_yyscan_t lp_yyscanner ); + +YYSTYPE * lp_yyget_lval (lp_yyscan_t lp_yyscanner ); + +void lp_yyset_lval (YYSTYPE * lp_yylval_param ,lp_yyscan_t lp_yyscanner ); + +/* Macros after this point can all be overridden by user definitions in + * section 1. + */ + +#ifndef YY_SKIP_YYWRAP +#ifdef __cplusplus +extern "C" int lp_yywrap (lp_yyscan_t lp_yyscanner ); +#else +extern int lp_yywrap (lp_yyscan_t lp_yyscanner ); +#endif +#endif + + static void lp_yyunput (int c,char *buf_ptr ,lp_yyscan_t lp_yyscanner); + +#ifndef lp_yytext_ptr +static void lp_yy_flex_strncpy (char *,lp_yyconst char *,int ,lp_yyscan_t lp_yyscanner); +#endif + +#ifdef YY_NEED_STRLEN +static int lp_yy_flex_strlen (lp_yyconst char * ,lp_yyscan_t lp_yyscanner); +#endif + +#ifndef YY_NO_INPUT + +#ifdef __cplusplus +static int lp_yyinput (lp_yyscan_t lp_yyscanner ); +#else +static int input (lp_yyscan_t lp_yyscanner ); +#endif + +#endif + +/* Amount of stuff to slurp up with each read. */ +#ifndef YY_READ_BUF_SIZE +#define YY_READ_BUF_SIZE 8192 +#endif + +/* Copy whatever the last rule matched to the standard output. */ +#ifndef ECHO +/* This used to be an fputs(), but since the string might contain NUL's, + * we now use fwrite(). + */ +#define ECHO fwrite( lp_yytext, lp_yyleng, 1, lp_yyout ) +#endif + +/* Gets input and stuffs it into "buf". number of characters read, or YY_NULL, + * is returned in "result". + */ +#ifndef YY_INPUT +#define YY_INPUT(buf,result,max_size) \ + if ( YY_CURRENT_BUFFER_LVALUE->lp_yy_is_interactive ) \ + { \ + int c = '*'; \ + int n; \ + for ( n = 0; n < max_size && \ + (c = getc( lp_yyin )) != EOF && c != '\n'; ++n ) \ + buf[n] = (char) c; \ + if ( c == '\n' ) \ + buf[n++] = (char) c; \ + if ( c == EOF && ferror( lp_yyin ) ) \ + YY_FATAL_ERROR( "input in flex scanner failed" ); \ + result = n; \ + } \ + else \ + { \ + errno=0; \ + while ( (result = fread(buf, 1, max_size, lp_yyin))==0 && ferror(lp_yyin)) \ + { \ + if( errno != EINTR) \ + { \ + YY_FATAL_ERROR( "input in flex scanner failed" ); \ + break; \ + } \ + errno=0; \ + clearerr(lp_yyin); \ + } \ + }\ +\ + +#endif + +/* No semi-colon after return; correct usage is to write "lp_yyterminate();" - + * we don't want an extra ';' after the "return" because that will cause + * some compilers to complain about unreachable statements. + */ +#ifndef lp_yyterminate +#define lp_yyterminate() return YY_NULL +#endif + +/* Number of entries by which start-condition stack grows. */ +#ifndef YY_START_STACK_INCR +#define YY_START_STACK_INCR 25 +#endif + +/* Report a fatal error. */ +#ifndef YY_FATAL_ERROR +#define YY_FATAL_ERROR(msg) lp_yy_fatal_error( msg , lp_yyscanner) +#endif + +/* end tables serialization structures and prototypes */ + +/* Default declaration of generated scanner - a define so the user can + * easily add parameters. + */ +#ifndef YY_DECL +#define YY_DECL_IS_OURS 1 + +extern int lp_yylex \ + (YYSTYPE * lp_yylval_param ,lp_yyscan_t lp_yyscanner); + +#define YY_DECL int lp_yylex \ + (YYSTYPE * lp_yylval_param , lp_yyscan_t lp_yyscanner) +#endif /* !YY_DECL */ + +/* Code executed at the beginning of each rule, after lp_yytext and lp_yyleng + * have been set up. + */ +#ifndef YY_USER_ACTION +#define YY_USER_ACTION +#endif + +/* Code executed at the end of each rule. */ +#ifndef YY_BREAK +#define YY_BREAK break; +#endif + +#define YY_RULE_SETUP \ + if ( lp_yyleng > 0 ) \ + YY_CURRENT_BUFFER_LVALUE->lp_yy_at_bol = \ + (lp_yytext[lp_yyleng - 1] == '\n'); \ + YY_USER_ACTION + +/** The main scanner function which does all the work. + */ +YY_DECL +{ + register lp_yy_state_type lp_yy_current_state; + register char *lp_yy_cp, *lp_yy_bp; + register int lp_yy_act; + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + lp_yylval = lp_yylval_param; + + if ( !lp_yyg->lp_yy_init ) + { + lp_yyg->lp_yy_init = 1; + +#ifdef YY_USER_INIT + YY_USER_INIT; +#endif + + if ( ! lp_yyg->lp_yy_start ) + lp_yyg->lp_yy_start = 1; /* first start state */ + + if ( ! lp_yyin ) + lp_yyin = stdin; + + if ( ! lp_yyout ) + lp_yyout = stdout; + + if ( ! YY_CURRENT_BUFFER ) { + lp_yyensure_buffer_stack (lp_yyscanner); + YY_CURRENT_BUFFER_LVALUE = + lp_yy_create_buffer(lp_yyin,YY_BUF_SIZE ,lp_yyscanner); + } + + lp_yy_load_buffer_state(lp_yyscanner ); + } + + while ( 1 ) /* loops until end-of-file is reached */ + { + lp_yy_cp = lp_yyg->lp_yy_c_buf_p; + + /* Support of lp_yytext. */ + *lp_yy_cp = lp_yyg->lp_yy_hold_char; + + /* lp_yy_bp points to the position in lp_yy_ch_buf of the start of + * the current run. + */ + lp_yy_bp = lp_yy_cp; + + lp_yy_current_state = lp_yyg->lp_yy_start; + lp_yy_current_state += YY_AT_BOL(); +lp_yy_match: + do + { + register YY_CHAR lp_yy_c = lp_yy_ec[YY_SC_TO_UI(*lp_yy_cp)]; + if ( lp_yy_accept[lp_yy_current_state] ) + { + lp_yyg->lp_yy_last_accepting_state = lp_yy_current_state; + lp_yyg->lp_yy_last_accepting_cpos = lp_yy_cp; + } + while ( lp_yy_chk[lp_yy_base[lp_yy_current_state] + lp_yy_c] != lp_yy_current_state ) + { + lp_yy_current_state = (int) lp_yy_def[lp_yy_current_state]; + if ( lp_yy_current_state >= 144 ) + lp_yy_c = lp_yy_meta[(unsigned int) lp_yy_c]; + } + lp_yy_current_state = lp_yy_nxt[lp_yy_base[lp_yy_current_state] + (unsigned int) lp_yy_c]; + ++lp_yy_cp; + } + while ( lp_yy_base[lp_yy_current_state] != 388 ); + +lp_yy_find_action: + lp_yy_act = lp_yy_accept[lp_yy_current_state]; + if ( lp_yy_act == 0 ) + { /* have to back up */ + lp_yy_cp = lp_yyg->lp_yy_last_accepting_cpos; + lp_yy_current_state = lp_yyg->lp_yy_last_accepting_state; + lp_yy_act = lp_yy_accept[lp_yy_current_state]; + } + + YY_DO_BEFORE_ACTION; + + if ( lp_yy_act != YY_END_OF_BUFFER && lp_yy_rule_can_match_eol[lp_yy_act] ) + { + int lp_yyl; + for ( lp_yyl = 0; lp_yyl < lp_yyleng; ++lp_yyl ) + if ( lp_yytext[lp_yyl] == '\n' ) + + do{ lp_yylineno++; + lp_yycolumn=0; + }while(0) +; + } + +do_action: /* This label is used only to access EOF actions. */ + + switch ( lp_yy_act ) + { /* beginning of action switch */ + case 0: /* must back up */ + /* undo the effects of YY_DO_BEFORE_ACTION */ + *lp_yy_cp = lp_yyg->lp_yy_hold_char; + lp_yy_cp = lp_yyg->lp_yy_last_accepting_cpos; + lp_yy_current_state = lp_yyg->lp_yy_last_accepting_state; + goto lp_yy_find_action; + +case 1: +YY_RULE_SETUP +{ + BEGIN COMMENT; +} /* begin skip comment */ + YY_BREAK +case 2: +YY_RULE_SETUP +{ + BEGIN INITIAL; +} /* end skip comment */ + YY_BREAK +case 3: +YY_RULE_SETUP +{ +} + YY_BREAK +case 4: +/* rule 4 can match eol */ +YY_RULE_SETUP +{ +} + YY_BREAK +case 5: +YY_RULE_SETUP +{ +} + YY_BREAK +case 6: +YY_RULE_SETUP +{ + BEGIN LINECOMMENT; +} /* begin skip LINECOMMENT */ + YY_BREAK +case 7: +/* rule 7 can match eol */ +YY_RULE_SETUP +{ + BEGIN INITIAL; +} /* end skip LINECOMMENT */ + YY_BREAK +case 8: +YY_RULE_SETUP +{ + BEGIN INITIAL; +} /* end skip LINECOMMENT */ + YY_BREAK +case 9: +YY_RULE_SETUP +{ +} + YY_BREAK +case 10: +/* rule 10 can match eol */ +YY_RULE_SETUP +{ +} + YY_BREAK +case 11: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + + pp->lineno = lp_yylineno; + return(COMMA); +} + YY_BREAK +case 12: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + + pp->lineno = lp_yylineno; + return(MINIMISE); +} + YY_BREAK +case 13: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + + pp->lineno = lp_yylineno; + return(MAXIMISE); +} + YY_BREAK +case 14: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + + pp->lineno = lp_yylineno; + pv->f = atof((char *)lp_yytext); + return(INTCONS); +} /* f contains the last float */ + YY_BREAK +case 15: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + + pp->lineno = lp_yylineno; + pv->f = atof((char *)lp_yytext); + return(CONS); +} /* f contains the last float */ + YY_BREAK +case 16: +/* rule 16 can match eol */ +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + char *ptr, c; + + pp->lineno = lp_yylineno; + pv->f = DEF_INFINITE; + pv->Sign = 0; + ptr = (char *)lp_yytext; + while (isspace(*ptr)) ptr++; + if(*ptr == '-') + pv->Sign = 1; + if(lp_yyleng > 0) { + c = lp_yytext[lp_yyleng - 1]; + if(!isalnum(c)) + unput(c); + } + return(INF); +} /* f contains the last float */ + YY_BREAK +case 17: +/* rule 17 can match eol */ +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + int x; + + pp->lineno = lp_yylineno; + pv->Sign = 0; + for(x = 0; x < lp_yyleng; x++) + if(lp_yytext[x] == '-' || lp_yytext[x] == '+') + pv->Sign = (pv->Sign == (lp_yytext[x] == '+')); + return (TOK_SIGN); + /* Sign is TRUE if the sign-string + represents a '-'. Otherwise Sign + is FALSE */ +} + YY_BREAK +case 18: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + + pp->lineno = lp_yylineno; + if((!pv->Within_int_decl) && (!pv->Within_sec_decl) && (!pv->Within_sos_decl) && (!pv->Within_free_decl)) { + pv->Within_int_decl = 1; + pv->Within_sos_decl1 = FALSE; + } + return(SEC_INT); +} + YY_BREAK +case 19: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + + pp->lineno = lp_yylineno; + if((!pv->Within_int_decl) && (!pv->Within_sec_decl) && (!pv->Within_sos_decl) && (!pv->Within_free_decl)) { + pv->Within_int_decl = 2; + pv->Within_sos_decl1 = FALSE; + } + return(SEC_BIN); +} + YY_BREAK +case 20: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + + pp->lineno = lp_yylineno; + if((!pv->Within_int_decl) && (!pv->Within_sec_decl) && (!pv->Within_sos_decl) && (!pv->Within_free_decl)) { + pv->Within_sec_decl = TRUE; + pv->Within_sos_decl1 = FALSE; + } + return(SEC_SEC); +} + YY_BREAK +case 21: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + + pp->lineno = lp_yylineno; + if(!pv->Within_sos_decl) + pv->SOStype0 = (short)atoi(((char *)lp_yytext) + 3); + if((!pv->Within_int_decl) && (!pv->Within_sec_decl) && (!pv->Within_sos_decl) && (!pv->Within_free_decl)) + pv->Within_sos_decl = TRUE; + return(SEC_SOS); +} + YY_BREAK +case 22: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + + pp->lineno = lp_yylineno; + FREE(pv->Last_var); + pv->Last_var = strdup((char *)lp_yytext); + pv->Last_var[strlen(pv->Last_var) - 2] = 0; + return(SOSDESCR); +} + YY_BREAK +case 23: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + + pp->lineno = lp_yylineno; + if((!pv->Within_int_decl) && (!pv->Within_sec_decl) && (!pv->Within_sos_decl) && (!pv->Within_free_decl)) { + pv->Within_free_decl = TRUE; + pv->Within_sos_decl1 = FALSE; + } + return(SEC_FREE); +} + YY_BREAK +case 24: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + char *ptr; + + pp->lineno = lp_yylineno; + FREE(pv->Last_var); + pv->Last_var = strdup((char *)lp_yytext); + ptr = pv->Last_var + strlen(pv->Last_var); + ptr[-1] = ' '; + while ((--ptr >= pv->Last_var) && (isspace(*ptr))) + *ptr = 0; + return(VARIABLECOLON); +} + YY_BREAK +case 25: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + + pp->lineno = lp_yylineno; + FREE(pv->Last_var); + pv->Last_var = strdup((char *)lp_yytext); + return(VAR); +} + YY_BREAK +case 26: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + + pp->lineno = lp_yylineno; + return (COLON); +} + YY_BREAK +case 27: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + + pp->lineno = lp_yylineno; + return(AR_M_OP); +} + YY_BREAK +case 28: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + + pp->lineno = lp_yylineno; + pv->OP = *lp_yytext; + return(RE_OPEQ); +} + YY_BREAK +case 29: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + + pp->lineno = lp_yylineno; + pv->OP = *lp_yytext; + return(RE_OPLE); +} + YY_BREAK +case 30: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + + pp->lineno = lp_yylineno; + pv->OP = *lp_yytext; + return(RE_OPGE); +} + YY_BREAK +case 31: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + parse_vars *pv = (parse_vars *) pp->parse_vars; + + pp->lineno = lp_yylineno; + pv->Within_int_decl = pv->Within_sec_decl = pv->Within_sos_decl = pv->Within_free_decl = FALSE; + check_int_sec_sos_free_decl(pp, (int) pv->Within_int_decl, (int) pv->Within_sec_decl, (int) pv->Within_sos_decl, (int) pv->Within_free_decl); + return(END_C); +} + YY_BREAK +case 32: +YY_RULE_SETUP +{ + parse_parm *pp = PARM; + + pp->lineno = lp_yylineno; + report(NULL, CRITICAL, "LEX ERROR : %s lineno %d\n", lp_yytext, lp_yylineno); + return(UNDEFINED); +} + YY_BREAK +case 33: +YY_RULE_SETUP +ECHO; + YY_BREAK +case YY_STATE_EOF(INITIAL): +case YY_STATE_EOF(COMMENT): +case YY_STATE_EOF(LINECOMMENT): + lp_yyterminate(); + + case YY_END_OF_BUFFER: + { + /* Amount of text matched not including the EOB char. */ + int lp_yy_amount_of_matched_text = (int) (lp_yy_cp - lp_yyg->lp_yytext_ptr) - 1; + + /* Undo the effects of YY_DO_BEFORE_ACTION. */ + *lp_yy_cp = lp_yyg->lp_yy_hold_char; + YY_RESTORE_YY_MORE_OFFSET + + if ( YY_CURRENT_BUFFER_LVALUE->lp_yy_buffer_status == YY_BUFFER_NEW ) + { + /* We're scanning a new file or input source. It's + * possible that this happened because the user + * just pointed lp_yyin at a new source and called + * lp_yylex(). If so, then we have to assure + * consistency between YY_CURRENT_BUFFER and our + * globals. Here is the right place to do so, because + * this is the first action (other than possibly a + * back-up) that will match for the new input source. + */ + lp_yyg->lp_yy_n_chars = YY_CURRENT_BUFFER_LVALUE->lp_yy_n_chars; + YY_CURRENT_BUFFER_LVALUE->lp_yy_input_file = lp_yyin; + YY_CURRENT_BUFFER_LVALUE->lp_yy_buffer_status = YY_BUFFER_NORMAL; + } + + /* Note that here we test for lp_yy_c_buf_p "<=" to the position + * of the first EOB in the buffer, since lp_yy_c_buf_p will + * already have been incremented past the NUL character + * (since all states make transitions on EOB to the + * end-of-buffer state). Contrast this with the test + * in input(). + */ + if ( lp_yyg->lp_yy_c_buf_p <= &YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf[lp_yyg->lp_yy_n_chars] ) + { /* This was really a NUL. */ + lp_yy_state_type lp_yy_next_state; + + lp_yyg->lp_yy_c_buf_p = lp_yyg->lp_yytext_ptr + lp_yy_amount_of_matched_text; + + lp_yy_current_state = lp_yy_get_previous_state( lp_yyscanner ); + + /* Okay, we're now positioned to make the NUL + * transition. We couldn't have + * lp_yy_get_previous_state() go ahead and do it + * for us because it doesn't know how to deal + * with the possibility of jamming (and we don't + * want to build jamming into it because then it + * will run more slowly). + */ + + lp_yy_next_state = lp_yy_try_NUL_trans( lp_yy_current_state , lp_yyscanner); + + lp_yy_bp = lp_yyg->lp_yytext_ptr + YY_MORE_ADJ; + + if ( lp_yy_next_state ) + { + /* Consume the NUL. */ + lp_yy_cp = ++lp_yyg->lp_yy_c_buf_p; + lp_yy_current_state = lp_yy_next_state; + goto lp_yy_match; + } + + else + { + lp_yy_cp = lp_yyg->lp_yy_c_buf_p; + goto lp_yy_find_action; + } + } + + else switch ( lp_yy_get_next_buffer( lp_yyscanner ) ) + { + case EOB_ACT_END_OF_FILE: + { + lp_yyg->lp_yy_did_buffer_switch_on_eof = 0; + + if ( lp_yywrap(lp_yyscanner ) ) + { + /* Note: because we've taken care in + * lp_yy_get_next_buffer() to have set up + * lp_yytext, we can now set up + * lp_yy_c_buf_p so that if some total + * hoser (like flex itself) wants to + * call the scanner after we return the + * YY_NULL, it'll still work - another + * YY_NULL will get returned. + */ + lp_yyg->lp_yy_c_buf_p = lp_yyg->lp_yytext_ptr + YY_MORE_ADJ; + + lp_yy_act = YY_STATE_EOF(YY_START); + goto do_action; + } + + else + { + if ( ! lp_yyg->lp_yy_did_buffer_switch_on_eof ) + YY_NEW_FILE; + } + break; + } + + case EOB_ACT_CONTINUE_SCAN: + lp_yyg->lp_yy_c_buf_p = + lp_yyg->lp_yytext_ptr + lp_yy_amount_of_matched_text; + + lp_yy_current_state = lp_yy_get_previous_state( lp_yyscanner ); + + lp_yy_cp = lp_yyg->lp_yy_c_buf_p; + lp_yy_bp = lp_yyg->lp_yytext_ptr + YY_MORE_ADJ; + goto lp_yy_match; + + case EOB_ACT_LAST_MATCH: + lp_yyg->lp_yy_c_buf_p = + &YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf[lp_yyg->lp_yy_n_chars]; + + lp_yy_current_state = lp_yy_get_previous_state( lp_yyscanner ); + + lp_yy_cp = lp_yyg->lp_yy_c_buf_p; + lp_yy_bp = lp_yyg->lp_yytext_ptr + YY_MORE_ADJ; + goto lp_yy_find_action; + } + break; + } + + default: + YY_FATAL_ERROR( + "fatal flex scanner internal error--no action found" ); + } /* end of action switch */ + } /* end of scanning one token */ +} /* end of lp_yylex */ + +/* lp_yy_get_next_buffer - try to read in a new buffer + * + * Returns a code representing an action: + * EOB_ACT_LAST_MATCH - + * EOB_ACT_CONTINUE_SCAN - continue scanning from current position + * EOB_ACT_END_OF_FILE - end of file + */ +static int lp_yy_get_next_buffer (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + register char *dest = YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf; + register char *source = lp_yyg->lp_yytext_ptr; + register int number_to_move, i; + int ret_val; + + if ( lp_yyg->lp_yy_c_buf_p > &YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf[lp_yyg->lp_yy_n_chars + 1] ) + YY_FATAL_ERROR( + "fatal flex scanner internal error--end of buffer missed" ); + + if ( YY_CURRENT_BUFFER_LVALUE->lp_yy_fill_buffer == 0 ) + { /* Don't try to fill the buffer, so this is an EOF. */ + if ( lp_yyg->lp_yy_c_buf_p - lp_yyg->lp_yytext_ptr - YY_MORE_ADJ == 1 ) + { + /* We matched a single character, the EOB, so + * treat this as a final EOF. + */ + return EOB_ACT_END_OF_FILE; + } + + else + { + /* We matched some text prior to the EOB, first + * process it. + */ + return EOB_ACT_LAST_MATCH; + } + } + + /* Try to read more data. */ + + /* First move last chars to start of buffer. */ + number_to_move = (int) (lp_yyg->lp_yy_c_buf_p - lp_yyg->lp_yytext_ptr) - 1; + + for ( i = 0; i < number_to_move; ++i ) + *(dest++) = *(source++); + + if ( YY_CURRENT_BUFFER_LVALUE->lp_yy_buffer_status == YY_BUFFER_EOF_PENDING ) + /* don't do the read, it's not guaranteed to return an EOF, + * just force an EOF + */ + YY_CURRENT_BUFFER_LVALUE->lp_yy_n_chars = lp_yyg->lp_yy_n_chars = 0; + + else + { + int num_to_read = + YY_CURRENT_BUFFER_LVALUE->lp_yy_buf_size - number_to_move - 1; + + while ( num_to_read <= 0 ) + { /* Not enough room in the buffer - grow it. */ + + /* just a shorter name for the current buffer */ + YY_BUFFER_STATE b = YY_CURRENT_BUFFER; + + int lp_yy_c_buf_p_offset = + (int) (lp_yyg->lp_yy_c_buf_p - b->lp_yy_ch_buf); + + if ( b->lp_yy_is_our_buffer ) + { + int new_size = b->lp_yy_buf_size * 2; + + if ( new_size <= 0 ) + b->lp_yy_buf_size += b->lp_yy_buf_size / 8; + else + b->lp_yy_buf_size *= 2; + + b->lp_yy_ch_buf = (char *) + /* Include room in for 2 EOB chars. */ + lp_yyrealloc((void *) b->lp_yy_ch_buf,b->lp_yy_buf_size + 2 ,lp_yyscanner ); + } + else + /* Can't grow it, we don't own it. */ + b->lp_yy_ch_buf = 0; + + if ( ! b->lp_yy_ch_buf ) + YY_FATAL_ERROR( + "fatal error - scanner input buffer overflow" ); + + lp_yyg->lp_yy_c_buf_p = &b->lp_yy_ch_buf[lp_yy_c_buf_p_offset]; + + num_to_read = YY_CURRENT_BUFFER_LVALUE->lp_yy_buf_size - + number_to_move - 1; + + } + + if ( num_to_read > YY_READ_BUF_SIZE ) + num_to_read = YY_READ_BUF_SIZE; + + /* Read in more data. */ + YY_INPUT( (&YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf[number_to_move]), + lp_yyg->lp_yy_n_chars, (size_t) num_to_read ); + + YY_CURRENT_BUFFER_LVALUE->lp_yy_n_chars = lp_yyg->lp_yy_n_chars; + } + + if ( lp_yyg->lp_yy_n_chars == 0 ) + { + if ( number_to_move == YY_MORE_ADJ ) + { + ret_val = EOB_ACT_END_OF_FILE; + lp_yyrestart(lp_yyin ,lp_yyscanner); + } + + else + { + ret_val = EOB_ACT_LAST_MATCH; + YY_CURRENT_BUFFER_LVALUE->lp_yy_buffer_status = + YY_BUFFER_EOF_PENDING; + } + } + + else + ret_val = EOB_ACT_CONTINUE_SCAN; + + if ((lp_yy_size_t) (lp_yyg->lp_yy_n_chars + number_to_move) > YY_CURRENT_BUFFER_LVALUE->lp_yy_buf_size) { + /* Extend the array by 50%, plus the number we really need. */ + lp_yy_size_t new_size = lp_yyg->lp_yy_n_chars + number_to_move + (lp_yyg->lp_yy_n_chars >> 1); + YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf = (char *) lp_yyrealloc((void *) YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf,new_size ,lp_yyscanner ); + if ( ! YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf ) + YY_FATAL_ERROR( "out of dynamic memory in lp_yy_get_next_buffer()" ); + } + + lp_yyg->lp_yy_n_chars += number_to_move; + YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf[lp_yyg->lp_yy_n_chars] = YY_END_OF_BUFFER_CHAR; + YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf[lp_yyg->lp_yy_n_chars + 1] = YY_END_OF_BUFFER_CHAR; + + lp_yyg->lp_yytext_ptr = &YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf[0]; + + return ret_val; +} + +/* lp_yy_get_previous_state - get the state just before the EOB char was reached */ + + static lp_yy_state_type lp_yy_get_previous_state (lp_yyscan_t lp_yyscanner) +{ + register lp_yy_state_type lp_yy_current_state; + register char *lp_yy_cp; + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + lp_yy_current_state = lp_yyg->lp_yy_start; + lp_yy_current_state += YY_AT_BOL(); + + for ( lp_yy_cp = lp_yyg->lp_yytext_ptr + YY_MORE_ADJ; lp_yy_cp < lp_yyg->lp_yy_c_buf_p; ++lp_yy_cp ) + { + register YY_CHAR lp_yy_c = (*lp_yy_cp ? lp_yy_ec[YY_SC_TO_UI(*lp_yy_cp)] : 1); + if ( lp_yy_accept[lp_yy_current_state] ) + { + lp_yyg->lp_yy_last_accepting_state = lp_yy_current_state; + lp_yyg->lp_yy_last_accepting_cpos = lp_yy_cp; + } + while ( lp_yy_chk[lp_yy_base[lp_yy_current_state] + lp_yy_c] != lp_yy_current_state ) + { + lp_yy_current_state = (int) lp_yy_def[lp_yy_current_state]; + if ( lp_yy_current_state >= 144 ) + lp_yy_c = lp_yy_meta[(unsigned int) lp_yy_c]; + } + lp_yy_current_state = lp_yy_nxt[lp_yy_base[lp_yy_current_state] + (unsigned int) lp_yy_c]; + } + + return lp_yy_current_state; +} + +/* lp_yy_try_NUL_trans - try to make a transition on the NUL character + * + * synopsis + * next_state = lp_yy_try_NUL_trans( current_state ); + */ + static lp_yy_state_type lp_yy_try_NUL_trans (lp_yy_state_type lp_yy_current_state , lp_yyscan_t lp_yyscanner) +{ + register int lp_yy_is_jam; + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; /* This var may be unused depending upon options. */ + register char *lp_yy_cp = lp_yyg->lp_yy_c_buf_p; + + register YY_CHAR lp_yy_c = 1; + if ( lp_yy_accept[lp_yy_current_state] ) + { + lp_yyg->lp_yy_last_accepting_state = lp_yy_current_state; + lp_yyg->lp_yy_last_accepting_cpos = lp_yy_cp; + } + while ( lp_yy_chk[lp_yy_base[lp_yy_current_state] + lp_yy_c] != lp_yy_current_state ) + { + lp_yy_current_state = (int) lp_yy_def[lp_yy_current_state]; + if ( lp_yy_current_state >= 144 ) + lp_yy_c = lp_yy_meta[(unsigned int) lp_yy_c]; + } + lp_yy_current_state = lp_yy_nxt[lp_yy_base[lp_yy_current_state] + (unsigned int) lp_yy_c]; + lp_yy_is_jam = (lp_yy_current_state == 143); + + return lp_yy_is_jam ? 0 : lp_yy_current_state; +} + + static void lp_yyunput (int c, register char * lp_yy_bp , lp_yyscan_t lp_yyscanner) +{ + register char *lp_yy_cp; + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + lp_yy_cp = lp_yyg->lp_yy_c_buf_p; + + /* undo effects of setting up lp_yytext */ + *lp_yy_cp = lp_yyg->lp_yy_hold_char; + + if ( lp_yy_cp < YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf + 2 ) + { /* need to shift things up to make room */ + /* +2 for EOB chars. */ + register int number_to_move = lp_yyg->lp_yy_n_chars + 2; + register char *dest = &YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf[ + YY_CURRENT_BUFFER_LVALUE->lp_yy_buf_size + 2]; + register char *source = + &YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf[number_to_move]; + + while ( source > YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf ) + *--dest = *--source; + + lp_yy_cp += (int) (dest - source); + lp_yy_bp += (int) (dest - source); + YY_CURRENT_BUFFER_LVALUE->lp_yy_n_chars = + lp_yyg->lp_yy_n_chars = YY_CURRENT_BUFFER_LVALUE->lp_yy_buf_size; + + if ( lp_yy_cp < YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf + 2 ) + YY_FATAL_ERROR( "flex scanner push-back overflow" ); + } + + *--lp_yy_cp = (char) c; + + if ( c == '\n' ){ + --lp_yylineno; + } + + lp_yyg->lp_yytext_ptr = lp_yy_bp; + lp_yyg->lp_yy_hold_char = *lp_yy_cp; + lp_yyg->lp_yy_c_buf_p = lp_yy_cp; +} + +#ifndef YY_NO_INPUT +#ifdef __cplusplus + static inline int lp_yyinput (lp_yyscan_t lp_yyscanner) +#else + static inline int input (lp_yyscan_t lp_yyscanner) +#endif + +{ + int c; + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + *lp_yyg->lp_yy_c_buf_p = lp_yyg->lp_yy_hold_char; + + if ( *lp_yyg->lp_yy_c_buf_p == YY_END_OF_BUFFER_CHAR ) + { + /* lp_yy_c_buf_p now points to the character we want to return. + * If this occurs *before* the EOB characters, then it's a + * valid NUL; if not, then we've hit the end of the buffer. + */ + if ( lp_yyg->lp_yy_c_buf_p < &YY_CURRENT_BUFFER_LVALUE->lp_yy_ch_buf[lp_yyg->lp_yy_n_chars] ) + /* This was really a NUL. */ + *lp_yyg->lp_yy_c_buf_p = '\0'; + + else + { /* need more input */ + int offset = lp_yyg->lp_yy_c_buf_p - lp_yyg->lp_yytext_ptr; + ++lp_yyg->lp_yy_c_buf_p; + + switch ( lp_yy_get_next_buffer( lp_yyscanner ) ) + { + case EOB_ACT_LAST_MATCH: + /* This happens because lp_yy_g_n_b() + * sees that we've accumulated a + * token and flags that we need to + * try matching the token before + * proceeding. But for input(), + * there's no matching to consider. + * So convert the EOB_ACT_LAST_MATCH + * to EOB_ACT_END_OF_FILE. + */ + + /* Reset buffer status. */ + lp_yyrestart(lp_yyin ,lp_yyscanner); + + /*FALLTHROUGH*/ + + case EOB_ACT_END_OF_FILE: + { + if ( lp_yywrap(lp_yyscanner ) ) + return EOF; + + if ( ! lp_yyg->lp_yy_did_buffer_switch_on_eof ) + YY_NEW_FILE; +#ifdef __cplusplus + return lp_yyinput(lp_yyscanner); +#else + return input(lp_yyscanner); +#endif + } + + case EOB_ACT_CONTINUE_SCAN: + lp_yyg->lp_yy_c_buf_p = lp_yyg->lp_yytext_ptr + offset; + break; + } + } + } + + c = *(unsigned char *) lp_yyg->lp_yy_c_buf_p; /* cast for 8-bit char's */ + *lp_yyg->lp_yy_c_buf_p = '\0'; /* preserve lp_yytext */ + lp_yyg->lp_yy_hold_char = *++lp_yyg->lp_yy_c_buf_p; + + YY_CURRENT_BUFFER_LVALUE->lp_yy_at_bol = (c == '\n'); + if ( YY_CURRENT_BUFFER_LVALUE->lp_yy_at_bol ) + + do{ lp_yylineno++; + lp_yycolumn=0; + }while(0) +; + + return c; +} +#endif /* ifndef YY_NO_INPUT */ + +/** Immediately switch to a different input stream. + * @param input_file A readable stream. + * @param lp_yyscanner The scanner object. + * @note This function does not reset the start condition to @c INITIAL . + */ + void lp_yyrestart (FILE * input_file , lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + if ( ! YY_CURRENT_BUFFER ){ + lp_yyensure_buffer_stack (lp_yyscanner); + YY_CURRENT_BUFFER_LVALUE = + lp_yy_create_buffer(lp_yyin,YY_BUF_SIZE ,lp_yyscanner); + } + + lp_yy_init_buffer(YY_CURRENT_BUFFER,input_file ,lp_yyscanner); + lp_yy_load_buffer_state(lp_yyscanner ); +} + +/** Switch to a different input buffer. + * @param new_buffer The new input buffer. + * @param lp_yyscanner The scanner object. + */ + void lp_yy_switch_to_buffer (YY_BUFFER_STATE new_buffer , lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + /* TODO. We should be able to replace this entire function body + * with + * lp_yypop_buffer_state(); + * lp_yypush_buffer_state(new_buffer); + */ + lp_yyensure_buffer_stack (lp_yyscanner); + if ( YY_CURRENT_BUFFER == new_buffer ) + return; + + if ( YY_CURRENT_BUFFER ) + { + /* Flush out information for old buffer. */ + *lp_yyg->lp_yy_c_buf_p = lp_yyg->lp_yy_hold_char; + YY_CURRENT_BUFFER_LVALUE->lp_yy_buf_pos = lp_yyg->lp_yy_c_buf_p; + YY_CURRENT_BUFFER_LVALUE->lp_yy_n_chars = lp_yyg->lp_yy_n_chars; + } + + YY_CURRENT_BUFFER_LVALUE = new_buffer; + lp_yy_load_buffer_state(lp_yyscanner ); + + /* We don't actually know whether we did this switch during + * EOF (lp_yywrap()) processing, but the only time this flag + * is looked at is after lp_yywrap() is called, so it's safe + * to go ahead and always set it. + */ + lp_yyg->lp_yy_did_buffer_switch_on_eof = 1; +} + +static void lp_yy_load_buffer_state (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + lp_yyg->lp_yy_n_chars = YY_CURRENT_BUFFER_LVALUE->lp_yy_n_chars; + lp_yyg->lp_yytext_ptr = lp_yyg->lp_yy_c_buf_p = YY_CURRENT_BUFFER_LVALUE->lp_yy_buf_pos; + lp_yyin = YY_CURRENT_BUFFER_LVALUE->lp_yy_input_file; + lp_yyg->lp_yy_hold_char = *lp_yyg->lp_yy_c_buf_p; +} + +/** Allocate and initialize an input buffer state. + * @param file A readable stream. + * @param size The character buffer size in bytes. When in doubt, use @c YY_BUF_SIZE. + * @param lp_yyscanner The scanner object. + * @return the allocated buffer state. + */ + YY_BUFFER_STATE lp_yy_create_buffer (FILE * file, int size , lp_yyscan_t lp_yyscanner) +{ + YY_BUFFER_STATE b; + + b = (YY_BUFFER_STATE) lp_yyalloc(sizeof( struct lp_yy_buffer_state ) ,lp_yyscanner ); + if ( ! b ) + YY_FATAL_ERROR( "out of dynamic memory in lp_yy_create_buffer()" ); + + b->lp_yy_buf_size = size; + + /* lp_yy_ch_buf has to be 2 characters longer than the size given because + * we need to put in 2 end-of-buffer characters. + */ + b->lp_yy_ch_buf = (char *) lp_yyalloc(b->lp_yy_buf_size + 2 ,lp_yyscanner ); + if ( ! b->lp_yy_ch_buf ) + YY_FATAL_ERROR( "out of dynamic memory in lp_yy_create_buffer()" ); + + b->lp_yy_is_our_buffer = 1; + + lp_yy_init_buffer(b,file ,lp_yyscanner); + + return b; +} + +/** Destroy the buffer. + * @param b a buffer created with lp_yy_create_buffer() + * @param lp_yyscanner The scanner object. + */ + void lp_yy_delete_buffer (YY_BUFFER_STATE b , lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + if ( ! b ) + return; + + if ( b == YY_CURRENT_BUFFER ) /* Not sure if we should pop here. */ + YY_CURRENT_BUFFER_LVALUE = (YY_BUFFER_STATE) 0; + + if ( b->lp_yy_is_our_buffer ) + lp_yyfree((void *) b->lp_yy_ch_buf ,lp_yyscanner ); + + lp_yyfree((void *) b ,lp_yyscanner ); +} + +#ifndef __cplusplus +extern int isatty (int ); +#endif /* __cplusplus */ + +/* Initializes or reinitializes a buffer. + * This function is sometimes called more than once on the same buffer, + * such as during a lp_yyrestart() or at EOF. + */ + static void lp_yy_init_buffer (YY_BUFFER_STATE b, FILE * file , lp_yyscan_t lp_yyscanner) + +{ + int oerrno = errno; + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + lp_yy_flush_buffer(b ,lp_yyscanner); + + b->lp_yy_input_file = file; + b->lp_yy_fill_buffer = 1; + + /* If b is the current buffer, then lp_yy_init_buffer was _probably_ + * called from lp_yyrestart() or through lp_yy_get_next_buffer. + * In that case, we don't want to reset the lineno or column. + */ + if (b != YY_CURRENT_BUFFER){ + b->lp_yy_bs_lineno = 1; + b->lp_yy_bs_column = 0; + } + + b->lp_yy_is_interactive = file ? (isatty( fileno(file) ) > 0) : 0; + + errno = oerrno; +} + +/** Discard all buffered characters. On the next scan, YY_INPUT will be called. + * @param b the buffer state to be flushed, usually @c YY_CURRENT_BUFFER. + * @param lp_yyscanner The scanner object. + */ + void lp_yy_flush_buffer (YY_BUFFER_STATE b , lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + if ( ! b ) + return; + + b->lp_yy_n_chars = 0; + + /* We always need two end-of-buffer characters. The first causes + * a transition to the end-of-buffer state. The second causes + * a jam in that state. + */ + b->lp_yy_ch_buf[0] = YY_END_OF_BUFFER_CHAR; + b->lp_yy_ch_buf[1] = YY_END_OF_BUFFER_CHAR; + + b->lp_yy_buf_pos = &b->lp_yy_ch_buf[0]; + + b->lp_yy_at_bol = 1; + b->lp_yy_buffer_status = YY_BUFFER_NEW; + + if ( b == YY_CURRENT_BUFFER ) + lp_yy_load_buffer_state(lp_yyscanner ); +} + +/** Pushes the new state onto the stack. The new state becomes + * the current state. This function will allocate the stack + * if necessary. + * @param new_buffer The new state. + * @param lp_yyscanner The scanner object. + */ +void lp_yypush_buffer_state (YY_BUFFER_STATE new_buffer , lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + if (new_buffer == NULL) + return; + + lp_yyensure_buffer_stack(lp_yyscanner); + + /* This block is copied from lp_yy_switch_to_buffer. */ + if ( YY_CURRENT_BUFFER ) + { + /* Flush out information for old buffer. */ + *lp_yyg->lp_yy_c_buf_p = lp_yyg->lp_yy_hold_char; + YY_CURRENT_BUFFER_LVALUE->lp_yy_buf_pos = lp_yyg->lp_yy_c_buf_p; + YY_CURRENT_BUFFER_LVALUE->lp_yy_n_chars = lp_yyg->lp_yy_n_chars; + } + + /* Only push if top exists. Otherwise, replace top. */ + if (YY_CURRENT_BUFFER) + lp_yyg->lp_yy_buffer_stack_top++; + YY_CURRENT_BUFFER_LVALUE = new_buffer; + + /* copied from lp_yy_switch_to_buffer. */ + lp_yy_load_buffer_state(lp_yyscanner ); + lp_yyg->lp_yy_did_buffer_switch_on_eof = 1; +} + +/** Removes and deletes the top of the stack, if present. + * The next element becomes the new top. + * @param lp_yyscanner The scanner object. + */ +void lp_yypop_buffer_state (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + if (!YY_CURRENT_BUFFER) + return; + + lp_yy_delete_buffer(YY_CURRENT_BUFFER ,lp_yyscanner); + YY_CURRENT_BUFFER_LVALUE = NULL; + if (lp_yyg->lp_yy_buffer_stack_top > 0) + --lp_yyg->lp_yy_buffer_stack_top; + + if (YY_CURRENT_BUFFER) { + lp_yy_load_buffer_state(lp_yyscanner ); + lp_yyg->lp_yy_did_buffer_switch_on_eof = 1; + } +} + +/* Allocates the stack if it does not exist. + * Guarantees space for at least one push. + */ +static void lp_yyensure_buffer_stack (lp_yyscan_t lp_yyscanner) +{ + int num_to_alloc; + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + if (!lp_yyg->lp_yy_buffer_stack) { + + /* First allocation is just for 2 elements, since we don't know if this + * scanner will even need a stack. We use 2 instead of 1 to avoid an + * immediate realloc on the next call. + */ + num_to_alloc = 1; + lp_yyg->lp_yy_buffer_stack = (struct lp_yy_buffer_state**)lp_yyalloc + (num_to_alloc * sizeof(struct lp_yy_buffer_state*) + , lp_yyscanner); + if ( ! lp_yyg->lp_yy_buffer_stack ) + YY_FATAL_ERROR( "out of dynamic memory in lp_yyensure_buffer_stack()" ); + + memset(lp_yyg->lp_yy_buffer_stack, 0, num_to_alloc * sizeof(struct lp_yy_buffer_state*)); + + lp_yyg->lp_yy_buffer_stack_max = num_to_alloc; + lp_yyg->lp_yy_buffer_stack_top = 0; + return; + } + + if (lp_yyg->lp_yy_buffer_stack_top >= (lp_yyg->lp_yy_buffer_stack_max) - 1){ + + /* Increase the buffer to prepare for a possible push. */ + int grow_size = 8 /* arbitrary grow size */; + + num_to_alloc = lp_yyg->lp_yy_buffer_stack_max + grow_size; + lp_yyg->lp_yy_buffer_stack = (struct lp_yy_buffer_state**)lp_yyrealloc + (lp_yyg->lp_yy_buffer_stack, + num_to_alloc * sizeof(struct lp_yy_buffer_state*) + , lp_yyscanner); + if ( ! lp_yyg->lp_yy_buffer_stack ) + YY_FATAL_ERROR( "out of dynamic memory in lp_yyensure_buffer_stack()" ); + + /* zero only the new slots.*/ + memset(lp_yyg->lp_yy_buffer_stack + lp_yyg->lp_yy_buffer_stack_max, 0, grow_size * sizeof(struct lp_yy_buffer_state*)); + lp_yyg->lp_yy_buffer_stack_max = num_to_alloc; + } +} + +/** Setup the input buffer state to scan directly from a user-specified character buffer. + * @param base the character buffer + * @param size the size in bytes of the character buffer + * @param lp_yyscanner The scanner object. + * @return the newly allocated buffer state object. + */ +YY_BUFFER_STATE lp_yy_scan_buffer (char * base, lp_yy_size_t size , lp_yyscan_t lp_yyscanner) +{ + YY_BUFFER_STATE b; + + if ( size < 2 || + base[size-2] != YY_END_OF_BUFFER_CHAR || + base[size-1] != YY_END_OF_BUFFER_CHAR ) + /* They forgot to leave room for the EOB's. */ + return 0; + + b = (YY_BUFFER_STATE) lp_yyalloc(sizeof( struct lp_yy_buffer_state ) ,lp_yyscanner ); + if ( ! b ) + YY_FATAL_ERROR( "out of dynamic memory in lp_yy_scan_buffer()" ); + + b->lp_yy_buf_size = size - 2; /* "- 2" to take care of EOB's */ + b->lp_yy_buf_pos = b->lp_yy_ch_buf = base; + b->lp_yy_is_our_buffer = 0; + b->lp_yy_input_file = 0; + b->lp_yy_n_chars = b->lp_yy_buf_size; + b->lp_yy_is_interactive = 0; + b->lp_yy_at_bol = 1; + b->lp_yy_fill_buffer = 0; + b->lp_yy_buffer_status = YY_BUFFER_NEW; + + lp_yy_switch_to_buffer(b ,lp_yyscanner ); + + return b; +} + +/** Setup the input buffer state to scan a string. The next call to lp_yylex() will + * scan from a @e copy of @a str. + * @param lp_yystr a NUL-terminated string to scan + * @param lp_yyscanner The scanner object. + * @return the newly allocated buffer state object. + * @note If you want to scan bytes that may contain NUL values, then use + * lp_yy_scan_bytes() instead. + */ +YY_BUFFER_STATE lp_yy_scan_string (lp_yyconst char * lp_yystr , lp_yyscan_t lp_yyscanner) +{ + + return lp_yy_scan_bytes(lp_yystr,strlen(lp_yystr) ,lp_yyscanner); +} + +/** Setup the input buffer state to scan the given bytes. The next call to lp_yylex() will + * scan from a @e copy of @a bytes. + * @param bytes the byte buffer to scan + * @param len the number of bytes in the buffer pointed to by @a bytes. + * @param lp_yyscanner The scanner object. + * @return the newly allocated buffer state object. + */ +YY_BUFFER_STATE lp_yy_scan_bytes (lp_yyconst char * lp_yybytes, int _lp_yybytes_len , lp_yyscan_t lp_yyscanner) +{ + YY_BUFFER_STATE b; + char *buf; + lp_yy_size_t n; + int i; + + /* Get memory for full buffer, including space for trailing EOB's. */ + n = _lp_yybytes_len + 2; + buf = (char *) lp_yyalloc(n ,lp_yyscanner ); + if ( ! buf ) + YY_FATAL_ERROR( "out of dynamic memory in lp_yy_scan_bytes()" ); + + for ( i = 0; i < _lp_yybytes_len; ++i ) + buf[i] = lp_yybytes[i]; + + buf[_lp_yybytes_len] = buf[_lp_yybytes_len+1] = YY_END_OF_BUFFER_CHAR; + + b = lp_yy_scan_buffer(buf,n ,lp_yyscanner); + if ( ! b ) + YY_FATAL_ERROR( "bad buffer in lp_yy_scan_bytes()" ); + + /* It's okay to grow etc. this buffer, and we should throw it + * away when we're done. + */ + b->lp_yy_is_our_buffer = 1; + + return b; +} + +#ifndef YY_EXIT_FAILURE +#define YY_EXIT_FAILURE 2 +#endif + +static void lp_yy_fatal_error (lp_yyconst char* msg , lp_yyscan_t lp_yyscanner) +{ + (void) fprintf( stderr, "%s\n", msg ); + exit( YY_EXIT_FAILURE ); +} + +/* Redefine lp_yyless() so it works in section 3 code. */ + +#undef lp_yyless +#define lp_yyless(n) \ + do \ + { \ + /* Undo effects of setting up lp_yytext. */ \ + int lp_yyless_macro_arg = (n); \ + YY_LESS_LINENO(lp_yyless_macro_arg);\ + lp_yytext[lp_yyleng] = lp_yyg->lp_yy_hold_char; \ + lp_yyg->lp_yy_c_buf_p = lp_yytext + lp_yyless_macro_arg; \ + lp_yyg->lp_yy_hold_char = *lp_yyg->lp_yy_c_buf_p; \ + *lp_yyg->lp_yy_c_buf_p = '\0'; \ + lp_yyleng = lp_yyless_macro_arg; \ + } \ + while ( 0 ) + +/* Accessor methods (get/set functions) to struct members. */ + +/** Get the user-defined data for this scanner. + * @param lp_yyscanner The scanner object. + */ +YY_EXTRA_TYPE lp_yyget_extra (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + return lp_yyextra; +} + +/** Get the current line number. + * @param lp_yyscanner The scanner object. + */ +int lp_yyget_lineno (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + if (! YY_CURRENT_BUFFER) + return 0; + + return lp_yylineno; +} + +/** Get the current column number. + * @param lp_yyscanner The scanner object. + */ +int lp_yyget_column (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + if (! YY_CURRENT_BUFFER) + return 0; + + return lp_yycolumn; +} + +/** Get the input stream. + * @param lp_yyscanner The scanner object. + */ +FILE *lp_yyget_in (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + return lp_yyin; +} + +/** Get the output stream. + * @param lp_yyscanner The scanner object. + */ +FILE *lp_yyget_out (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + return lp_yyout; +} + +/** Get the length of the current token. + * @param lp_yyscanner The scanner object. + */ +int lp_yyget_leng (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + return lp_yyleng; +} + +/** Get the current token. + * @param lp_yyscanner The scanner object. + */ + +char *lp_yyget_text (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + return lp_yytext; +} + +/** Set the user-defined data. This data is never touched by the scanner. + * @param user_defined The data to be associated with this scanner. + * @param lp_yyscanner The scanner object. + */ +void lp_yyset_extra (YY_EXTRA_TYPE user_defined , lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + lp_yyextra = user_defined ; +} + +/** Set the current line number. + * @param line_number + * @param lp_yyscanner The scanner object. + */ +void lp_yyset_lineno (int line_number , lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + /* lineno is only valid if an input buffer exists. */ + if (! YY_CURRENT_BUFFER ) + lp_yy_fatal_error( "lp_yyset_lineno called with no buffer" , lp_yyscanner); + + lp_yylineno = line_number; +} + +/** Set the current column. + * @param line_number + * @param lp_yyscanner The scanner object. + */ +void lp_yyset_column (int column_no , lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + /* column is only valid if an input buffer exists. */ + if (! YY_CURRENT_BUFFER ) + lp_yy_fatal_error( "lp_yyset_column called with no buffer" , lp_yyscanner); + + lp_yycolumn = column_no; +} + +/** Set the input stream. This does not discard the current + * input buffer. + * @param in_str A readable stream. + * @param lp_yyscanner The scanner object. + * @see lp_yy_switch_to_buffer + */ +void lp_yyset_in (FILE * in_str , lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + lp_yyin = in_str ; +} + +void lp_yyset_out (FILE * out_str , lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + lp_yyout = out_str ; +} + +int lp_yyget_debug (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + return lp_yy_flex_debug; +} + +void lp_yyset_debug (int bdebug , lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + lp_yy_flex_debug = bdebug ; +} + +/* Accessor methods for lp_yylval and lp_yylloc */ + +YYSTYPE * lp_yyget_lval (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + return lp_yylval; +} + +void lp_yyset_lval (YYSTYPE * lp_yylval_param , lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + lp_yylval = lp_yylval_param; +} + +/* User-visible API */ + +/* lp_yylex_init is special because it creates the scanner itself, so it is + * the ONLY reentrant function that doesn't take the scanner as the last argument. + * That's why we explicitly handle the declaration, instead of using our macros. + */ + +int lp_yylex_init(lp_yyscan_t* ptr_lp_yy_globals) + +{ + if (ptr_lp_yy_globals == NULL){ + errno = EINVAL; + return 1; + } + + *ptr_lp_yy_globals = (lp_yyscan_t) lp_yyalloc ( sizeof( struct lp_yyguts_t ), NULL ); + + if (*ptr_lp_yy_globals == NULL){ + errno = ENOMEM; + return 1; + } + + /* By setting to 0xAA, we expose bugs in lp_yy_init_globals. Leave at 0x00 for releases. */ + memset(*ptr_lp_yy_globals,0x00,sizeof(struct lp_yyguts_t)); + + return lp_yy_init_globals ( *ptr_lp_yy_globals ); +} + +/* lp_yylex_init_extra has the same functionality as lp_yylex_init, but follows the + * convention of taking the scanner as the last argument. Note however, that + * this is a *pointer* to a scanner, as it will be allocated by this call (and + * is the reason, too, why this function also must handle its own declaration). + * The user defined value in the first argument will be available to lp_yyalloc in + * the lp_yyextra field. + */ + +int lp_yylex_init_extra(YY_EXTRA_TYPE lp_yy_user_defined,lp_yyscan_t* ptr_lp_yy_globals ) + +{ + struct lp_yyguts_t dummy_lp_yyguts; + + lp_yyset_extra (lp_yy_user_defined, &dummy_lp_yyguts); + + if (ptr_lp_yy_globals == NULL){ + errno = EINVAL; + return 1; + } + + *ptr_lp_yy_globals = (lp_yyscan_t) lp_yyalloc ( sizeof( struct lp_yyguts_t ), &dummy_lp_yyguts ); + + if (*ptr_lp_yy_globals == NULL){ + errno = ENOMEM; + return 1; + } + + /* By setting to 0xAA, we expose bugs in + lp_yy_init_globals. Leave at 0x00 for releases. */ + memset(*ptr_lp_yy_globals,0x00,sizeof(struct lp_yyguts_t)); + + lp_yyset_extra (lp_yy_user_defined, *ptr_lp_yy_globals); + + return lp_yy_init_globals ( *ptr_lp_yy_globals ); +} + +static int lp_yy_init_globals (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + /* Initialization is the same as for the non-reentrant scanner. + * This function is called from lp_yylex_destroy(), so don't allocate here. + */ + + lp_yyg->lp_yy_buffer_stack = 0; + lp_yyg->lp_yy_buffer_stack_top = 0; + lp_yyg->lp_yy_buffer_stack_max = 0; + lp_yyg->lp_yy_c_buf_p = (char *) 0; + lp_yyg->lp_yy_init = 0; + lp_yyg->lp_yy_start = 0; + + lp_yyg->lp_yy_start_stack_ptr = 0; + lp_yyg->lp_yy_start_stack_depth = 0; + lp_yyg->lp_yy_start_stack = NULL; + +/* Defined in main.c */ +#ifdef YY_STDINIT + lp_yyin = stdin; + lp_yyout = stdout; +#else + lp_yyin = (FILE *) 0; + lp_yyout = (FILE *) 0; +#endif + + /* For future reference: Set errno on error, since we are called by + * lp_yylex_init() + */ + return 0; +} + +/* lp_yylex_destroy is for both reentrant and non-reentrant scanners. */ +int lp_yylex_destroy (lp_yyscan_t lp_yyscanner) +{ + struct lp_yyguts_t * lp_yyg = (struct lp_yyguts_t*)lp_yyscanner; + + /* Pop the buffer stack, destroying each element. */ + while(YY_CURRENT_BUFFER){ + lp_yy_delete_buffer(YY_CURRENT_BUFFER ,lp_yyscanner ); + YY_CURRENT_BUFFER_LVALUE = NULL; + lp_yypop_buffer_state(lp_yyscanner); + } + + /* Destroy the stack itself. */ + lp_yyfree(lp_yyg->lp_yy_buffer_stack ,lp_yyscanner); + lp_yyg->lp_yy_buffer_stack = NULL; + + /* Destroy the start condition stack. */ + lp_yyfree(lp_yyg->lp_yy_start_stack ,lp_yyscanner ); + lp_yyg->lp_yy_start_stack = NULL; + + /* Reset the globals. This is important in a non-reentrant scanner so the next time + * lp_yylex() is called, initialization will occur. */ + lp_yy_init_globals( lp_yyscanner); + + /* Destroy the main struct (reentrant only). */ + lp_yyfree ( lp_yyscanner , lp_yyscanner ); + lp_yyscanner = NULL; + return 0; +} + +/* + * Internal utility routines. + */ + +#ifndef lp_yytext_ptr +static void lp_yy_flex_strncpy (char* s1, lp_yyconst char * s2, int n , lp_yyscan_t lp_yyscanner) +{ + register int i; + for ( i = 0; i < n; ++i ) + s1[i] = s2[i]; +} +#endif + +#ifdef YY_NEED_STRLEN +static int lp_yy_flex_strlen (lp_yyconst char * s , lp_yyscan_t lp_yyscanner) +{ + register int n; + for ( n = 0; s[n]; ++n ) + ; + + return n; +} +#endif + +void *lp_yyalloc (lp_yy_size_t size , lp_yyscan_t lp_yyscanner) +{ + return (void *) malloc( size ); +} + +void *lp_yyrealloc (void * ptr, lp_yy_size_t size , lp_yyscan_t lp_yyscanner) +{ + /* The cast to (char *) in the following accommodates both + * implementations that use char* generic pointers, and those + * that use void* generic pointers. It works with the latter + * because both ANSI C and C++ allow castless assignment from + * any pointer type to void*, and deal with argument conversions + * as though doing an assignment. + */ + return (void *) realloc( (char *) ptr, size ); +} + +void lp_yyfree (void * ptr , lp_yyscan_t lp_yyscanner) +{ + free( (char *) ptr ); /* see lp_yyrealloc() for (char *) cast */ +} + +#define YYTABLES_NAME "lp_yytables" + diff --git a/include/cartesian_geom/point.h b/include/cartesian_geom/point.h index a2209f4f5..b8c8e32b4 100644 --- a/include/cartesian_geom/point.h +++ b/include/cartesian_geom/point.h @@ -81,6 +81,7 @@ class point void set_dimension(const unsigned int dim) { d = dim; + coeffs.setZero(d); } void set_coord(const unsigned int i, FT coord) diff --git a/include/convex_bodies/ballintersectconvex.h b/include/convex_bodies/ballintersectconvex.h index 12965d08c..b61ae5324 100644 --- a/include/convex_bodies/ballintersectconvex.h +++ b/include/convex_bodies/ballintersectconvex.h @@ -28,7 +28,7 @@ class BallIntersectPolytope { BallIntersectPolytope() {} - BallIntersectPolytope(Polytope const& PP, CBall &BB) : P(PP), B(BB) {}; + BallIntersectPolytope(Polytope& PP, CBall &BB) : P(PP), B(BB) {}; Polytope first() const { return P; } CBall second() const { return B; } diff --git a/include/diagnostics/raftery.hpp b/include/diagnostics/raftery.hpp index de7928202..1fa087a5d 100644 --- a/include/diagnostics/raftery.hpp +++ b/include/diagnostics/raftery.hpp @@ -29,10 +29,9 @@ #define DIAGNOSTICS_RAFTERY_HPP template -NT fix(NT x) +NT round_to_zero(NT x) { - if (x > 0.0) return std::floor(x); - return std::ceil(x); + return (x > 0.0) ? std::floor(x) : std::ceil(x); } #include "raftery_subroutines/empquant.hpp" @@ -121,13 +120,13 @@ MT perform_raftery(MT const& samples, NT const& q, NT const& r, NT const& s) NT psum = alpha + beta; NT tmp1 = std::log(psum * epss / std::max(alpha, beta)) / std::log(std::abs(1.0 - psum)); - NT nburn = fix((tmp1 + 1.0) * NT(kthin)); + NT nburn = round_to_zero((tmp1 + 1.0) * NT(kthin)); NT phi = ppnd((s + 1.0) / 2.0); NT tmp2 = (2.0 - psum) * alpha * beta * (phi * phi) / (psum * psum * psum * (r * r)); - NT nprec = fix(tmp2 + 1.0) * kthin; - NT nmin = fix(((1.0 - q) * q * (phi * phi) / (r * r)) + 1.0); + NT nprec = round_to_zero(tmp2 + 1.0) * kthin; + NT nmin = round_to_zero(((1.0 - q) * q * (phi * phi) / (r * r)) + 1.0); NT irl = (nburn + nprec) / nmin; - NT kind = std::max(fix(irl + 1.0), NT(kmind)); + NT kind = std::max(round_to_zero(irl + 1.0), NT(kmind)); results(i, 0) = NT(kthin); results(i, 1) = NT(nburn); diff --git a/include/diagnostics/raftery_subroutines/empquant.hpp b/include/diagnostics/raftery_subroutines/empquant.hpp index cd305372e..6604f0ad7 100644 --- a/include/diagnostics/raftery_subroutines/empquant.hpp +++ b/include/diagnostics/raftery_subroutines/empquant.hpp @@ -17,10 +17,10 @@ template NT empquant(VT const& sorted_samples, NT const& q) { unsigned int n = sorted_samples.rows(); - + NT order = (n - 1) * q + 1.0; NT fract = order - NT(int(order)); - int low = std::max(fix(order), 1.0); + int low = std::max(round_to_zero(order), 1.0); int high = std::min(low + 1.0, NT(n)); NT y = (1.0 - fract) * sorted_samples(low - 1) + fract * sorted_samples(high - 1); diff --git a/include/ode_solvers/leapfrog.hpp b/include/ode_solvers/leapfrog.hpp index 4c1f299b0..b7f2724b4 100644 --- a/include/ode_solvers/leapfrog.hpp +++ b/include/ode_solvers/leapfrog.hpp @@ -56,6 +56,8 @@ struct LeapfrogODESolver { pts xs; pts xs_prev; + Point grad_x; + MT _AA; std::pair pbpair; @@ -69,11 +71,11 @@ struct LeapfrogODESolver { eta(step), eta0(step), t(initial_time), F(oracle), Ks(boundaries), xs(initial_state), adaptive(adaptive_) { dim = xs[0].dimension(); _update_parameters = update_parameters(); + grad_x.set_dimension(dim); initialize(); }; - void initialize() { for (unsigned int i = 0; i < xs.size(); i++) { VT ar, av; @@ -88,7 +90,6 @@ struct LeapfrogODESolver { Av.push_back(av); lambda_prev.push_back(NT(0)); } - //step(); } void disable_adaptive() { @@ -101,27 +102,28 @@ struct LeapfrogODESolver { void step(int k, bool accepted) { num_steps++; - if (adaptive) eta = (eta0 * num_steps) / (num_steps + num_reflections); xs_prev = xs; unsigned int x_index, v_index, it; t += eta; + Point y; for (unsigned int i = 1; i < xs.size(); i += 2) { - //pbpair.second = -1; + x_index = i - 1; v_index = i; // v' <- v + eta / 2 F(x) - Point z = F(v_index, xs_prev, t); - z = (eta / 2) * z; - xs[v_index] = xs[v_index] + z; + if (k == 0 && !accepted) { + grad_x = F(v_index, xs_prev, t); + } + xs[v_index] += (eta / 2) * grad_x; // x <- x + eta v' - Point y = xs[v_index]; + y = xs[v_index]; if (Ks[x_index] == NULL) { - xs[x_index] = xs_prev[x_index] + eta*y; + xs[x_index] += eta*y; } else { // Find intersection (assuming a line trajectory) between x and y @@ -173,10 +175,8 @@ struct LeapfrogODESolver { } // tilde v <- v + eta / 2 F(tilde x) - z = F(v_index, xs, t); - z = (eta / 2) * z; - xs[v_index] = xs[v_index] + z; - + grad_x = F(v_index, xs, t); + xs[v_index] += (eta / 2) * grad_x; } } @@ -202,6 +202,19 @@ struct LeapfrogODESolver { xs[index] = p; } + NT get_eta() { + return eta; + } + + void set_eta(NT &eta_) { + eta = eta_; + eta0 = eta_; + } + + bounds get_bounds() { + return Ks; + } + }; diff --git a/include/preprocess/estimate_L_smooth_parameter.hpp b/include/preprocess/estimate_L_smooth_parameter.hpp new file mode 100644 index 000000000..8f5a1e89c --- /dev/null +++ b/include/preprocess/estimate_L_smooth_parameter.hpp @@ -0,0 +1,73 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis + +//Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program. + +// Licensed under GNU LGPL.3, see LICENCE file + + +#ifndef ESTIMATE_L_SMOOTH_PARAMETER_HPP +#define ESTIMATE_L_SMOOTH_PARAMETER_HPP + +#include "random_walks/random_walks.hpp" + +template +< + typename WalkTypePolicy = AcceleratedBilliardWalk, + typename Polytope, + typename Point, + typename NegativeGradientFunctor, + typename RandomNumberGenerator +> +double estimate_L_smooth(Polytope &P, Point &p, unsigned int const& walk_length, + NegativeGradientFunctor F, RandomNumberGenerator &rng) +{ + typedef typename Point::FT NT; + typedef typename WalkTypePolicy::template Walk + < + Polytope, + RandomNumberGenerator + > RandomWalk; + + P.ComputeInnerBall(); + + unsigned int d = P.dimension(); + unsigned int rnum = 5 * d; + std::vector randPoints(1); + std::vector vecPoint1; + std::vector vecPoint2; + std::vector< std::vector > listOfPoints; + Point F1; + + RandomWalk walk(P, p, rng); + for (unsigned int i=0; i::lowest(), Ltemp; + + for (auto pit=listOfPoints.begin(); pit!=(listOfPoints.end()-1); ++pit) + { + F1 = F(1, *pit, 0); + + for (auto qit=(pit+1); qit!=listOfPoints.end(); ++qit) + { + vecPoint2 = *qit; + Ltemp = (F1 - F(1, *qit, 0)).length() / ((*pit)[0] - (*qit)[0]).length(); + + if (Ltemp > L) + { + L = Ltemp; + } + } + } + return L; +} + + +#endif diff --git a/include/random_walks/boundary_cdhr_walk.hpp b/include/random_walks/boundary_cdhr_walk.hpp index 40d35df2d..5e6b5224c 100644 --- a/include/random_walks/boundary_cdhr_walk.hpp +++ b/include/random_walks/boundary_cdhr_walk.hpp @@ -26,7 +26,7 @@ struct BCDHRWalk typedef typename Point::FT NT; template - Walk(GenericPolytope const& P, Point const& p, RandomNumberGenerator& rng) + Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator& rng) { initialize(P, p, rng); } diff --git a/include/random_walks/boundary_rdhr_walk.hpp b/include/random_walks/boundary_rdhr_walk.hpp index 0cf582957..29d500b6c 100644 --- a/include/random_walks/boundary_rdhr_walk.hpp +++ b/include/random_walks/boundary_rdhr_walk.hpp @@ -27,7 +27,7 @@ struct BRDHRWalk typedef typename Point::FT NT; template - Walk(GenericPolytope const& P, Point const& p, RandomNumberGenerator& rng) + Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator& rng) { initialize(P, p, rng); } diff --git a/include/random_walks/compute_diameter.hpp b/include/random_walks/compute_diameter.hpp index e3ec0f441..b20aa839a 100644 --- a/include/random_walks/compute_diameter.hpp +++ b/include/random_walks/compute_diameter.hpp @@ -205,7 +205,7 @@ template struct compute_diameter> { template - static NT compute(OrderPolytope const& P) + static NT compute(OrderPolytope& P) { return std::sqrt(NT(P.dimension())); } @@ -215,7 +215,7 @@ template struct compute_diameter, Ellipsoid > > { template - static NT compute(BallIntersectPolytope, Ellipsoid> const& P) + static NT compute(BallIntersectPolytope, Ellipsoid>& P) { NT polytope_diameter = std::sqrt(NT(P.dimension())); return std::min(polytope_diameter, (NT(2) * P.radius())); diff --git a/include/random_walks/exponential_hamiltonian_monte_carlo_exact_walk.hpp b/include/random_walks/exponential_hamiltonian_monte_carlo_exact_walk.hpp index 0d1199410..6e5d2e000 100644 --- a/include/random_walks/exponential_hamiltonian_monte_carlo_exact_walk.hpp +++ b/include/random_walks/exponential_hamiltonian_monte_carlo_exact_walk.hpp @@ -84,7 +84,7 @@ struct Walk < typename GenericPolytope > - inline bool apply(GenericPolytope const& P, + inline bool apply(GenericPolytope& P, Point& p, // a point to start unsigned int const& walk_length, RandomNumberGenerator &rng) @@ -104,7 +104,7 @@ struct Walk } T = rng.sample_urdist() * _Len; _v = GetDirection::apply(n, rng, false); - + it = 0; while (it < _rho) { @@ -122,7 +122,7 @@ struct Walk P.compute_reflection(_v, _p, pbpair.second); it++; } - + } while (P.is_in(_p, INSIDE_BODY_TOLLERANCE) == 0); if (it == _rho) { _p = p0; @@ -137,7 +137,7 @@ struct Walk < typename GenericPolytope > - inline void get_starting_point(GenericPolytope const& P, + inline void get_starting_point(GenericPolytope& P, Point const& center, Point &q, unsigned int const& walk_length, @@ -158,7 +158,7 @@ struct Walk < typename GenericPolytope > - inline void parameters_burnin(GenericPolytope const& P, + inline void parameters_burnin(GenericPolytope& P, Point const& center, unsigned int const& num_points, unsigned int const& walk_length, @@ -171,7 +171,7 @@ struct Walk pointset.push_back(_p); NT rad = NT(0), max_dist, Lmax = NT(0), radius = P.InnerBall().second; - for (int i = 0; i < num_points; i++) + for (int i = 0; i < num_points; i++) { p = GetPointInDsphere::apply(n, radius, rng); p += center; @@ -179,18 +179,18 @@ struct Walk apply(P, p, walk_length, rng); max_dist = get_max_distance(pointset, p, rad); - if (max_dist > Lmax) + if (max_dist > Lmax) { Lmax = max_dist; } - if (2.0 * rad > Lmax) + if (2.0 * rad > Lmax) { Lmax = 2.0 * rad; } pointset.push_back(p); } - if (Lmax > _Len) + if (Lmax > _Len) { update_delta(Lmax); } @@ -210,7 +210,7 @@ private : < typename GenericPolytope > - inline void initialize(GenericPolytope const& P, + inline void initialize(GenericPolytope& P, Point const& p, RandomNumberGenerator &rng) { @@ -219,7 +219,7 @@ private : _lambdas.setZero(P.num_of_hyperplanes()); _Av.setZero(P.num_of_hyperplanes()); _v = GetDirection::apply(n, rng, false); - + do { _p = p; T = rng.sample_urdist() * _Len; @@ -262,11 +262,11 @@ private : } - inline double get_max_distance(std::vector &pointset, Point const& q, double &rad) + inline double get_max_distance(std::vector &pointset, Point const& q, double &rad) { double dis = -1.0, temp_dis; int jj = 0; - for (auto vecit = pointset.begin(); vecit!=pointset.end(); vecit++, jj++) + for (auto vecit = pointset.begin(); vecit!=pointset.end(); vecit++, jj++) { temp_dis = (q.getCoefficients() - (*vecit).getCoefficients()).norm(); if (temp_dis > dis) { diff --git a/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/include/random_walks/gaussian_accelerated_billiard_walk.hpp index a44608fb2..fbf43f443 100644 --- a/include/random_walks/gaussian_accelerated_billiard_walk.hpp +++ b/include/random_walks/gaussian_accelerated_billiard_walk.hpp @@ -68,7 +68,7 @@ struct GaussianAcceleratedBilliardWalk typedef typename Point::FT NT; template - Walk(GenericPolytope const& P, + Walk(GenericPolytope& P, Point const& p, Ellipsoid const& E, // ellipsoid representing the Gaussian distribution RandomNumberGenerator &rng) @@ -83,7 +83,7 @@ struct GaussianAcceleratedBilliardWalk } template - Walk(GenericPolytope const& P, + Walk(GenericPolytope& P, Point const& p, Ellipsoid const& E, // ellipsoid representing the Gaussian distribution RandomNumberGenerator &rng, @@ -104,7 +104,7 @@ struct GaussianAcceleratedBilliardWalk typename GenericPolytope, typename Ellipsoid > - inline void apply(GenericPolytope const& P, + inline void apply(GenericPolytope& P, Point &p, // a point to return the result Ellipsoid const& E, // ellipsoid representing the Gaussian distribution unsigned int const& walk_length, @@ -181,7 +181,7 @@ struct GaussianAcceleratedBilliardWalk typename GenericPolytope, typename Ellipsoid > - inline void initialize(GenericPolytope const& P, + inline void initialize(GenericPolytope& P, Point const& p, // a point to start Ellipsoid const& E, // ellipsoid representing the Gaussian distribution RandomNumberGenerator &rng) diff --git a/include/random_walks/gaussian_ball_walk.hpp b/include/random_walks/gaussian_ball_walk.hpp index c641b7344..8c267e0f9 100644 --- a/include/random_walks/gaussian_ball_walk.hpp +++ b/include/random_walks/gaussian_ball_walk.hpp @@ -46,19 +46,19 @@ struct Walk typedef typename Point::FT NT; template - static inline NT compute_delta(GenericPolytope const& P, NT const& a) + static inline NT compute_delta(GenericPolytope& P, NT const& a) { //return ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); return (NT(4) * (P.InnerBall()).second) / std::sqrt(std::max(NT(1), a) * NT(P.dimension())); } - Walk (Polytope const& P, Point const& p, NT const& a, + Walk (Polytope& P, Point const& p, NT const& a, RandomNumberGenerator &rng) { _delta = compute_delta(P, a); } - Walk (Polytope const& P, + Walk (Polytope& P, Point const& p, NT const& a, RandomNumberGenerator &rng, diff --git a/include/random_walks/gaussian_cdhr_walk.hpp b/include/random_walks/gaussian_cdhr_walk.hpp index 8962ce413..04ebedbc4 100644 --- a/include/random_walks/gaussian_cdhr_walk.hpp +++ b/include/random_walks/gaussian_cdhr_walk.hpp @@ -75,7 +75,7 @@ struct Walk typedef typename Polytope::PointType Point; typedef typename Point::FT NT; - Walk(Polytope const& P, + Walk(Polytope& P, Point const& p, NT const& a_i, RandomNumberGenerator &rng) @@ -83,7 +83,7 @@ struct Walk initialize(P, p, a_i, rng); } - Walk(Polytope const& P, + Walk(Polytope& P, Point const& p, NT const& a_i, RandomNumberGenerator& rng, diff --git a/include/random_walks/gaussian_hamiltonian_monte_carlo_exact_walk.hpp b/include/random_walks/gaussian_hamiltonian_monte_carlo_exact_walk.hpp index 954fda731..8b83b053b 100644 --- a/include/random_walks/gaussian_hamiltonian_monte_carlo_exact_walk.hpp +++ b/include/random_walks/gaussian_hamiltonian_monte_carlo_exact_walk.hpp @@ -19,7 +19,7 @@ struct GaussianHamiltonianMonteCarloExactWalk GaussianHamiltonianMonteCarloExactWalk(double L, unsigned int _rho) : param(L, true, _rho, true) {} - + GaussianHamiltonianMonteCarloExactWalk(double L) : param(L, true, 0, false) {} @@ -28,7 +28,7 @@ struct GaussianHamiltonianMonteCarloExactWalk : param(0, false, 0, false) {} - + struct parameters { parameters(double L, bool set, unsigned int _rho, bool _set_rho) @@ -80,7 +80,7 @@ struct Walk < typename GenericPolytope > - inline void apply(GenericPolytope const& P, + inline void apply(GenericPolytope& P, Point& p, NT const& a_i, unsigned int const& walk_length, @@ -120,7 +120,7 @@ struct Walk < typename GenericPolytope > - inline void get_starting_point(GenericPolytope const& P, + inline void get_starting_point(GenericPolytope& P, Point const& center, Point &q, unsigned int const& walk_length, @@ -141,7 +141,7 @@ struct Walk < typename GenericPolytope > - inline void parameters_burnin(GenericPolytope const& P, + inline void parameters_burnin(GenericPolytope& P, Point const& center, unsigned int const& num_points, unsigned int const& walk_length, @@ -154,7 +154,7 @@ struct Walk pointset.push_back(_p); NT rad = NT(0), max_dist, Lmax = NT(0), radius = P.InnerBall().second; - for (int i = 0; i < num_points; i++) + for (int i = 0; i < num_points; i++) { p = GetPointInDsphere::apply(n, radius, rng); p += center; @@ -162,18 +162,18 @@ struct Walk apply(P, p, walk_length, rng); max_dist = get_max_distance(pointset, p, rad); - if (max_dist > Lmax) + if (max_dist > Lmax) { Lmax = max_dist; } - if (2.0*rad > Lmax) + if (2.0*rad > Lmax) { Lmax = 2.0 * rad; } pointset.push_back(p); } - if (Lmax > _Len) + if (Lmax > _Len) { update_delta(Lmax); } @@ -191,7 +191,7 @@ private : < typename GenericPolytope > - inline void initialize(GenericPolytope const& P, + inline void initialize(GenericPolytope& P, Point const& p, NT const& a_i, RandomNumberGenerator &rng) @@ -239,14 +239,14 @@ private : p.set_coord(i, C * std::cos(omega * T + Phi)); v.set_coord(i, -C * omega * std::sin(omega * T + Phi)); } - + } - inline double get_max_distance(std::vector &pointset, Point const& q, double &rad) + inline double get_max_distance(std::vector &pointset, Point const& q, double &rad) { double dis = -1.0, temp_dis; int jj = 0; - for (auto vecit = pointset.begin(); vecit!=pointset.end(); vecit++, jj++) + for (auto vecit = pointset.begin(); vecit!=pointset.end(); vecit++, jj++) { temp_dis = (q.getCoefficients() - (*vecit).getCoefficients()).norm(); if (temp_dis > dis) { diff --git a/include/random_walks/gaussian_rdhr_walk.hpp b/include/random_walks/gaussian_rdhr_walk.hpp index de7fd7bf9..b7a53abe7 100644 --- a/include/random_walks/gaussian_rdhr_walk.hpp +++ b/include/random_walks/gaussian_rdhr_walk.hpp @@ -80,10 +80,10 @@ struct Walk typedef typename Polytope::PointType Point; typedef typename Point::FT NT; - Walk(Polytope const&, Point const&, NT const&, RandomNumberGenerator&) + Walk(Polytope&, Point const&, NT const&, RandomNumberGenerator&) {} - Walk(Polytope const&, Point const&, NT const&, RandomNumberGenerator&, + Walk(Polytope&, Point const&, NT const&, RandomNumberGenerator&, parameters&) {} diff --git a/include/random_walks/hamiltonian_monte_carlo_walk.hpp b/include/random_walks/hamiltonian_monte_carlo_walk.hpp index 3c50ca5e9..23342ea57 100644 --- a/include/random_walks/hamiltonian_monte_carlo_walk.hpp +++ b/include/random_walks/hamiltonian_monte_carlo_walk.hpp @@ -114,12 +114,10 @@ struct HamiltonianMonteCarloWalk { }; - inline void apply( - RandomNumberGenerator &rng, - int walk_length=1, - bool metropolis_filter=true) + inline void apply(RandomNumberGenerator &rng, + int walk_length=1, + bool metropolis_filter=true) { - num_runs++; // Pick a random velocity @@ -156,6 +154,7 @@ struct HamiltonianMonteCarloWalk { } } else { x = x_tilde; + accepted = true; } discard_ratio = (1.0 * total_discarded_samples) / num_runs; diff --git a/include/random_walks/multithread_walks.hpp b/include/random_walks/multithread_walks.hpp index 45b2b67cb..c42fba1fa 100644 --- a/include/random_walks/multithread_walks.hpp +++ b/include/random_walks/multithread_walks.hpp @@ -58,7 +58,7 @@ struct BCDHRWalk_multithread //typedef thread_params thread_parameters_; template - Walk(GenericPolytope const& P, thread_parameters_ ¶meters, RandomNumberGenerator& rng) + Walk(GenericPolytope& P, thread_parameters_ ¶meters, RandomNumberGenerator& rng) { initialize(P, parameters, rng); } @@ -68,7 +68,7 @@ struct BCDHRWalk_multithread typename BallPolytope > inline void apply(BallPolytope const& P, - thread_parameters_ ¶ms, // parameters + thread_parameters_ ¶ms, // parameters unsigned int const& walk_length, RandomNumberGenerator& rng) { @@ -97,13 +97,13 @@ struct BCDHRWalk_multithread template inline void initialize(GenericBody const& P, - thread_parameters_ ¶ms, // parameters + thread_parameters_ ¶ms, // parameters RandomNumberGenerator& rng) { params.lambdas.setZero(P.num_of_hyperplanes()); params.rand_coord = rng.sample_uidist(); NT kapa = rng.sample_urdist(); - + std::pair bpair = P.line_intersect_coord(params.p, params.rand_coord, params.lambdas); params.p_prev = params.p; @@ -158,7 +158,7 @@ struct BRDHRWalk_multithread //typedef thread_params thread_parameters_; template - Walk(GenericPolytope const& P, thread_parameters_ ¶meters, RandomNumberGenerator& rng) + Walk(GenericPolytope& P, thread_parameters_ ¶meters, RandomNumberGenerator& rng) { initialize(P, parameters, rng); } @@ -168,7 +168,7 @@ struct BRDHRWalk_multithread typename BallPolytope > inline void apply(BallPolytope const& P, - thread_parameters_ ¶ms, // parameters + thread_parameters_ ¶ms, // parameters unsigned int const& walk_length, RandomNumberGenerator& rng) { @@ -191,7 +191,7 @@ struct BRDHRWalk_multithread template inline void initialize(GenericBody const& P, - thread_parameters_ ¶ms, // parameters + thread_parameters_ ¶ms, // parameters RandomNumberGenerator& rng) { params.lambdas.setZero(P.num_of_hyperplanes()); @@ -243,7 +243,7 @@ struct Walk typedef thread_parameters thread_parameters_; //typedef thread_params thread_parameters_; - Walk(Polytope const& P, + Walk(Polytope& P, thread_parameters_ ¶ms, NT const& a_i, RandomNumberGenerator &rng) @@ -252,7 +252,7 @@ struct Walk } template - Walk(Polytope const& P, + Walk(Polytope& P, thread_parameters_ ¶ms, NT const& a_i, RandomNumberGenerator &rng, @@ -260,14 +260,14 @@ struct Walk { initialize(P, params, a_i, rng); } - + template < typename BallPolytope > inline void apply(BallPolytope const& P, - thread_parameters_ ¶ms, // parameters + thread_parameters_ ¶ms, // parameters NT const& a_i, unsigned int const& walk_length, RandomNumberGenerator &rng) @@ -293,7 +293,7 @@ private : template inline void initialize(BallPolytope const& P, - thread_parameters_ ¶ms, // parameters + thread_parameters_ ¶ms, // parameters NT const& a_i, RandomNumberGenerator &rng) { @@ -345,11 +345,11 @@ struct Walk typedef thread_parameters thread_parameters_; //typedef thread_params thread_parameters_; - Walk(Polytope const&, thread_parameters_ &, NT const&, RandomNumberGenerator&) + Walk(Polytope&, thread_parameters_ &, NT const&, RandomNumberGenerator&) {} template - Walk(Polytope const&, thread_parameters_ &, NT const&, RandomNumberGenerator&, + Walk(Polytope&, thread_parameters_ &, NT const&, RandomNumberGenerator&, parameters&) {} @@ -358,7 +358,7 @@ struct Walk typename BallPolytope > inline void apply(BallPolytope const& P, - thread_parameters_ ¶ms, // parameters + thread_parameters_ ¶ms, // parameters NT const& a_i, unsigned int const& walk_length, RandomNumberGenerator &rng) @@ -440,7 +440,7 @@ struct Walk //typedef thread_params thread_parameters_; template - Walk(GenericPolytope const& P, thread_parameters_ ¶meters, RandomNumberGenerator &rng) + Walk(GenericPolytope& P, thread_parameters_ ¶meters, RandomNumberGenerator &rng) { _Len = compute_diameter ::template compute(P); @@ -448,7 +448,7 @@ struct Walk } template - Walk(GenericPolytope const& P, thread_parameters_ ¶meters, RandomNumberGenerator &rng, + Walk(GenericPolytope& P, thread_parameters_ ¶meters, RandomNumberGenerator &rng, parameters_ const& params) { _Len = params.set_L ? params.m_L @@ -461,7 +461,7 @@ struct Walk < typename GenericPolytope > - inline void apply(GenericPolytope const& P, + inline void apply(GenericPolytope& P, thread_parameters_ ¶meters, unsigned int const& walk_length, RandomNumberGenerator &rng) @@ -491,7 +491,7 @@ struct Walk P.compute_reflection(parameters.v, parameters.p, pbpair.second); it++; } - if (it == 50*n) + if (it == 50*n) { parameters.p = parameters.p0; } @@ -509,7 +509,7 @@ private : < typename GenericPolytope > - inline void initialize(GenericPolytope const& P, + inline void initialize(GenericPolytope& P, thread_parameters_ ¶meters, RandomNumberGenerator &rng) { @@ -524,7 +524,7 @@ private : std::pair pbpair = P.line_positive_intersect(parameters.p, parameters.v, parameters.lambdas, parameters.Av); - if (T <= pbpair.first) + if (T <= pbpair.first) { parameters.p += (T * parameters.v); parameters.lambda_prev = T; @@ -538,9 +538,9 @@ private : while (it <= 50*n) { std::pair pbpair - = P.line_positive_intersect(parameters.p, parameters.v, parameters.lambdas, + = P.line_positive_intersect(parameters.p, parameters.v, parameters.lambdas, parameters.Av, parameters.lambda_prev); - if (T <= pbpair.first) + if (T <= pbpair.first) { parameters.p += (T * parameters.v); parameters.lambda_prev = T; @@ -599,7 +599,7 @@ struct Walk //typedef thread_params thread_parameters_; template - Walk(GenericPolytope const& P, thread_parameters_ ¶meters, RandomNumberGenerator& rng) + Walk(GenericPolytope& P, thread_parameters_ ¶meters, RandomNumberGenerator& rng) { initialize(P, parameters, rng); } @@ -610,7 +610,7 @@ struct Walk typename BallPolytope > inline void apply(BallPolytope const& P, - thread_parameters_ ¶ms, // parameters + thread_parameters_ ¶ms, // parameters unsigned int const& walk_length, RandomNumberGenerator &rng) { @@ -634,7 +634,7 @@ private : template inline void initialize(BallPolytope const& P, - thread_parameters_ ¶ms, // parameters + thread_parameters_ ¶ms, // parameters RandomNumberGenerator &rng) { params.lambdas.setZero(P.num_of_hyperplanes()); @@ -690,7 +690,7 @@ struct Walk //typedef thread_params thread_parameters_; template - Walk(GenericPolytope const& P, thread_parameters_ ¶meters, RandomNumberGenerator& rng) + Walk(GenericPolytope& P, thread_parameters_ ¶meters, RandomNumberGenerator& rng) { initialize(P, parameters, rng); } @@ -700,7 +700,7 @@ struct Walk typename BallPolytope > inline void apply(BallPolytope const& P, - thread_parameters_ ¶ms, // parameters + thread_parameters_ ¶ms, // parameters unsigned int const& walk_length, RandomNumberGenerator& rng) { @@ -719,7 +719,7 @@ private : template inline void initialize(BallPolytope const& P, - thread_parameters_ ¶ms, // parameters + thread_parameters_ ¶ms, // parameters RandomNumberGenerator &rng) { params.lambdas.setZero(P.num_of_hyperplanes()); diff --git a/include/random_walks/nuts_hmc_walk.hpp b/include/random_walks/nuts_hmc_walk.hpp new file mode 100644 index 000000000..857064c0e --- /dev/null +++ b/include/random_walks/nuts_hmc_walk.hpp @@ -0,0 +1,380 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2022 Vissarion Fisikopoulos +// Copyright (c) 2018-2022 Apostolos Chalkis +// Copyright (c) 2020-2022 Elias Tsigaridas +// Copyright (c) 2020-2022 Marios Papachristou + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Matthew D. Hoffman, Andrew Gelman. "The No-U-Turn Sampler: +// Adaptively Setting Path Lengths in Hamiltonian Monte Carlo", 2011. + +// Comment: Compared to [Matthew D. Hoffman, Andrew Gelman, 2011] +// we modify the step of Nesterov's algorithm in the burn in phase. + +#ifndef NUTS_HAMILTONIAN_MONTE_CARLO_WALK_HPP +#define NUTS_HAMILTONIAN_MONTE_CARLO_WALK_HPP + + +#include "generators/boost_random_number_generator.hpp" +#include "random_walks/gaussian_helpers.hpp" +#include "ode_solvers/ode_solvers.hpp" +#include "preprocess/estimate_L_smooth_parameter.hpp" + +struct NutsHamiltonianMonteCarloWalk { + + template + < + typename NT, + typename OracleFunctor + > + struct parameters { + NT epsilon; // tolerance in mixing + NT eta; // step size + + parameters( + OracleFunctor const& F, + unsigned int dim, + NT epsilon_=2) + { + epsilon = epsilon_; + eta = F.params.L > 0 ? 10.0 / (dim * sqrt(F.params.L)) : 0.005; + } + }; + + template + < + typename Point, + typename Polytope, + typename RandomNumberGenerator, + typename NegativeGradientFunctor, + typename NegativeLogprobFunctor, + typename Solver + > + struct Walk { + + typedef std::vector pts; + typedef typename Point::FT NT; + typedef std::vector bounds; + + // Hyperparameters of the sampler + parameters ¶ms; + + // Numerical ODE solver + Solver *solver; + + // Dimension + unsigned int dim; + + // Discarded Samples + long num_runs = 0; + long total_acceptance = 0; + + // Average acceptance probability + NT average_acceptance = 0; + + // References to xs + Point x, v; + + // Helper points + Point v_pl, v_min, v_min_j, v_pl_j, X_pl, X_pl_j, X_min, X, X_rnd_j, X_min_j, x_pl_min; + + // Gradient function + NegativeGradientFunctor &F; + + bool accepted; + + // Burnin parameters + NT eps_step, mu, log_tilde_eps, H_tilde, alpha, na; + const NT delta = NT(0.65), Delta_max = NT(1000), gamma = NT(0.05), t0 = NT(10), kk = NT(0.85); + + // Density exponent + NegativeLogprobFunctor &f; + + Walk(Polytope *P, + Point &p, + NegativeGradientFunctor &neg_grad_f, + NegativeLogprobFunctor &neg_logprob_f, + parameters ¶m, + bool burn_in_phase = true) : + params(param), + F(neg_grad_f), + f(neg_logprob_f) + { + dim = p.dimension(); + + v_pl.set_dimension(dim); + v_min.set_dimension(dim); + v_min_j.set_dimension(dim); + v_pl_j.set_dimension(dim); + X_pl.set_dimension(dim); + X_pl_j.set_dimension(dim); + X_min.set_dimension(dim); + X.set_dimension(dim); + X_rnd_j.set_dimension(dim); + X_min_j.set_dimension(dim); + x_pl_min.set_dimension(dim); + + eps_step = params.eta; + mu = std::log(10*eps_step); + log_tilde_eps = NT(0); + H_tilde = NT(0); + alpha = NT(0); + na = NT(0); + + // Starting point is provided from outside + x = p; + + accepted = false; + + // Initialize solver + solver = new Solver(0, params.eta, pts{x, x}, F, bounds{P, NULL}); + disable_adaptive(); + + if (burn_in_phase) + { + RandomNumberGenerator rng(dim); + burnin(rng); + } + }; + + + inline void burnin(RandomNumberGenerator &rng, + unsigned int N = 1000, + unsigned int walk_length=1) + { + reset_num_runs(); + Point p = x; + NT L; + + if ((solver->get_bounds())[0] == NULL) + { + L = (NT(100) / NT(dim)) * (NT(100) / NT(dim)); + } + else + { + Polytope K = *(solver->get_bounds())[0]; + L = estimate_L_smooth(K, p, walk_length, F, rng); + } + + eps_step = NT(5) / (NT(dim) * std::sqrt(L)); + solver->set_eta(eps_step); + + for (int i = 0; i < N; i++) + { + apply(rng, walk_length, true); + solver->set_eta(eps_step); + } + reset_num_runs(); + } + + + inline void apply(RandomNumberGenerator &rng, + unsigned int walk_length=1, + bool burnin = false) + { + num_runs++; + + int x_counting_total = 0; + + // Pick a random velocity + v = GetDirection::apply(dim, rng, false); + + v_pl = v; + v_min = NT(-1) * v; + X_pl = x; + X_min = x; + + NT h1 = hamiltonian(x, v); + + NT uu = std::log(rng.sample_urdist()) - h1; + int j = -1; + bool s = true; + bool updated = false; + bool pos_state_single = false; + + if (burnin) + { + alpha = NT(0); + } + + while (s) + { + j++; + + if (burnin) + { + na = std::pow(NT(2), NT(j)); + } + + NT dir = rng.sample_urdist(); + + if (dir > 0.5) + { + v = v_pl; + X = X_pl; + } + else + { + v = v_min; + X = X_min; + } + X_rnd_j = X; + + int x_counting = 0; + int num_samples = int(std::pow(NT(2), NT(j))); + accepted = false; + + for (int k = 1; k <= num_samples; k++) + { + if (!accepted) + { + solver->set_state(0, X); + solver->set_state(1, v); + } + + // Get proposals + solver->steps(walk_length, accepted); + accepted = true; + + X = solver->get_state(0); + v = solver->get_state(1); + + NT hj = hamiltonian(X, v); + + if (burnin) + { + alpha += std::min(NT(1), std::exp(-hj + h1)); + } + + if (uu > Delta_max - hj) + { + s = false; + break; + } + + bool pos_state = false; + if (uu < -hj) + { + pos_state = true; + pos_state_single = true; + x_counting = x_counting + 1; + x_counting_total = x_counting_total + 1; + } + + if (k == 1) + { + if (dir > 0.5) + { + X_min_j = X; + v_min_j = v; + } + else + { + X_pl_j = X; + v_pl_j = v; + } + } + if (k == num_samples) + { + if (dir > 0.5) + { + x_pl_min = X - X_min_j; + if ((x_pl_min.dot(v) < 0) || (x_pl_min.dot(v_min_j) < 0)) + { + s = false; + } + } + else + { + x_pl_min = X_pl_j - X; + if ((x_pl_min.dot(v) < 0) || (x_pl_min.dot(v_pl_j) < 0)) + { + s = false; + } + } + } + if ((rng.sample_urdist() < (1/NT(x_counting))) && pos_state) + { + X_rnd_j = X; + } + } + + if (dir > 0.5) + { + X_pl = X; + v_pl = v; + } + else + { + X_min = X; + v_min = v; + } + + if (s && (rng.sample_urdist() < (NT(x_counting) / NT(x_counting_total)))) + { + x = X_rnd_j; + if (pos_state_single) + { + updated = true; + } + } + + if (s) + { + x_pl_min = X_pl - X_min; + if ((x_pl_min.dot(v_min) < 0) || (x_pl_min.dot(v_pl) < 0)) + { + s = false; + } + } + } + + if (updated) + { + total_acceptance++; + } + average_acceptance = NT(total_acceptance) / NT(num_runs); + + if (burnin) + { + H_tilde = (NT(1) - NT(1) / (NT(num_runs) + t0)) * H_tilde + (NT(1) / (NT(num_runs) + t0)) * (delta - alpha / na); + NT log_eps = mu - (std::sqrt(NT(num_runs)) / gamma) * H_tilde; + + // TODO: use the following to generalize Nesterov's algorithm + //log_tilde_eps = std::pow(mu,-kk) * log_eps + (NT(1) - std::pow(mu,-kk))*log_tilde_eps; + + eps_step = std::exp(log_eps); + } + } + + inline NT hamiltonian(Point &pos, Point &vel) const { + return f(pos) + 0.5 * vel.dot(vel); + } + + inline NT get_eta_solver() { + return solver->get_eta(); + } + + void disable_adaptive() { + solver->disable_adaptive(); + } + + void enable_adaptive() { + solver->enable_adaptive(); + } + + void reset_num_runs() { + num_runs = 0; + total_acceptance = 0; + } + + NT get_ratio_acceptance() { + return average_acceptance; + } + }; +}; + +#endif // HAMILTONIAN_MONTE_CARLO_WALK_HPP diff --git a/include/random_walks/random_walks.hpp b/include/random_walks/random_walks.hpp index 4affee74d..6c36457c1 100644 --- a/include/random_walks/random_walks.hpp +++ b/include/random_walks/random_walks.hpp @@ -28,6 +28,7 @@ #include "random_walks/exponential_hamiltonian_monte_carlo_exact_walk.hpp" #include "random_walks/uniform_accelerated_billiard_walk_parallel.hpp" #include "random_walks/hamiltonian_monte_carlo_walk.hpp" +#include "random_walks/nuts_hmc_walk.hpp" #include "random_walks/langevin_walk.hpp" #include "random_walks/crhmc/crhmc_walk.hpp" #endif // RANDOM_WALKS_RANDOM_WALKS_HPP diff --git a/include/random_walks/uniform_ball_walk.hpp b/include/random_walks/uniform_ball_walk.hpp index aeee06537..77c608b72 100644 --- a/include/random_walks/uniform_ball_walk.hpp +++ b/include/random_walks/uniform_ball_walk.hpp @@ -46,14 +46,14 @@ struct BallWalk typedef typename Point::FT NT; template - Walk(GenericPolytope const& P, Point const& /*p*/, + Walk(GenericPolytope& P, Point const& /*p*/, RandomNumberGenerator& /*rng*/) { _delta = compute_delta(P); } template - Walk(GenericPolytope const& P, Point const& /*p*/, + Walk(GenericPolytope& P, Point const& /*p*/, RandomNumberGenerator& /*rng*/, parameters const& params) { _delta = params.set_delta ? params.m_L @@ -61,7 +61,7 @@ struct BallWalk } template - static inline NT compute_delta(GenericPolytope const& P) + static inline NT compute_delta(GenericPolytope& P) { //return ((P.InnerBall()).second * NT(4)) / NT(P.dimension()); return (NT(4) * (P.InnerBall()).second) / std::sqrt(NT(P.dimension())); diff --git a/include/random_walks/uniform_cdhr_walk.hpp b/include/random_walks/uniform_cdhr_walk.hpp index b1a1a5386..12d2ff78e 100644 --- a/include/random_walks/uniform_cdhr_walk.hpp +++ b/include/random_walks/uniform_cdhr_walk.hpp @@ -30,13 +30,13 @@ struct Walk typedef typename Point::FT NT; template - Walk(GenericPolytope const& P, Point const& p, RandomNumberGenerator& rng) + Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator& rng) { initialize(P, p, rng); } template - Walk(GenericPolytope const& P, Point const& p, + Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator& rng, parameters const& params) { initialize(P, p, rng); diff --git a/include/random_walks/uniform_rdhr_walk.hpp b/include/random_walks/uniform_rdhr_walk.hpp index a322fe12c..5e0eb4750 100644 --- a/include/random_walks/uniform_rdhr_walk.hpp +++ b/include/random_walks/uniform_rdhr_walk.hpp @@ -31,13 +31,13 @@ struct Walk typedef typename Point::FT NT; template - Walk(GenericPolytope const& P, Point const& p, RandomNumberGenerator& rng) + Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator& rng) { initialize(P, p, rng); } template - Walk(GenericPolytope const& P, Point const& p, + Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator& rng, parameters const& params) { initialize(P, p, rng); @@ -47,7 +47,7 @@ struct Walk < typename BallPolytope > - inline void apply(BallPolytope const& P, + inline void apply(BallPolytope& P, Point& p, // a point to start unsigned int const& walk_length, RandomNumberGenerator& rng) @@ -67,7 +67,7 @@ struct Walk private : template - inline void initialize(BallPolytope const& P, + inline void initialize(BallPolytope& P, Point const& p, RandomNumberGenerator &rng) { diff --git a/include/sampling/random_point_generators.hpp b/include/sampling/random_point_generators.hpp index 89ebdfa1f..e70e14617 100644 --- a/include/sampling/random_point_generators.hpp +++ b/include/sampling/random_point_generators.hpp @@ -143,7 +143,7 @@ struct GaussianRandomPointGenerator typename WalkPolicy, typename RandomNumberGenerator > - static void apply(Polytope const& P, + static void apply(Polytope& P, Point &p, // a point to start NT const& a_i, unsigned int const& rnum, @@ -170,7 +170,7 @@ struct GaussianRandomPointGenerator typename RandomNumberGenerator, typename Parameters > - static void apply(Polytope const& P, + static void apply(Polytope& P, Point &p, // a point to start NT const& a_i, unsigned int const& rnum, @@ -203,7 +203,7 @@ struct BoundaryRandomPointGenerator typename WalkPolicy, typename RandomNumberGenerator > - static void apply(Polytope const& P, + static void apply(Polytope& P, Point &p, // a point to start unsigned int const& rnum, unsigned int const& walk_length, @@ -231,26 +231,15 @@ struct LogconcaveRandomPointGenerator { template - < - typename Polytope, - typename Point, - typename PointList, - typename WalkPolicy, - typename RandomNumberGenerator, - typename NegativeGradientFunctor, - typename NegativeLogprobFunctor, - typename Parameters + < typename PointList, + typename WalkPolicy, + typename RandomNumberGenerator > - static void apply(Polytope &P, - Point &p, // a point to start - unsigned int const& rnum, + static void apply(unsigned int const& rnum, unsigned int const& walk_length, PointList &randPoints, WalkPolicy &policy, RandomNumberGenerator &rng, - NegativeGradientFunctor &F, - NegativeLogprobFunctor &f, - Parameters ¶meters, Walk &walk) { typedef double NT; @@ -341,7 +330,7 @@ struct ExponentialRandomPointGenerator typename WalkPolicy, typename RandomNumberGenerator > - static void apply(Polytope const& P, + static void apply(Polytope& P, Point &p, // a point to start Point const& c, // bias function NT const& T, // temperature/variance @@ -374,7 +363,7 @@ struct ExponentialRandomPointGenerator typename RandomNumberGenerator, typename Parameters > - static void apply(Polytope const& P, + static void apply(Polytope& P, Point &p, // a point to start Point const& c, // bias function NT const& T, // temperature/variance diff --git a/include/sampling/sampling.hpp b/include/sampling/sampling.hpp index aa432990b..696186098 100644 --- a/include/sampling/sampling.hpp +++ b/include/sampling/sampling.hpp @@ -23,10 +23,10 @@ #define SAMPLE_ONLY_H template void uniform_sampling(PointList &randPoints, Polytope &P, @@ -232,14 +232,14 @@ template < typename Solver > void logconcave_sampling(PointList &randPoints, - Polytope &P, - RandomNumberGenerator &rng, - const unsigned int &walk_len, - const unsigned int &rnum, - const Point &starting_point, - unsigned int const& nburns, - NegativeGradientFunctor &F, - NegativeLogprobFunctor &f) + Polytope &P, + RandomNumberGenerator &rng, + const unsigned int &walk_len, + const unsigned int &rnum, + const Point &starting_point, + unsigned int const& nburns, + NegativeGradientFunctor &F, + NegativeLogprobFunctor &f) { typedef typename WalkTypePolicy::template Walk < @@ -272,16 +272,16 @@ void logconcave_sampling(PointList &randPoints, walk logconcave_walk(&P, p, F, f, params); typedef LogconcaveRandomPointGenerator RandomPointGenerator; - - RandomPointGenerator::apply(P, p, nburns, walk_len, randPoints, - push_back_policy, rng, F, f, params, logconcave_walk); - + + if (nburns > 0) { + RandomPointGenerator::apply(nburns, walk_len, randPoints, + push_back_policy, rng, logconcave_walk); + } logconcave_walk.disable_adaptive(); randPoints.clear(); - RandomPointGenerator::apply(P, p, rnum, walk_len, randPoints, - push_back_policy, rng, F, f, params, logconcave_walk); - + RandomPointGenerator::apply(rnum, walk_len, randPoints, + push_back_policy, rng, logconcave_walk); } #include "preprocess/crhmc/crhmc_input.h" #include "preprocess/crhmc/crhmc_problem.h" diff --git a/include/volume/volume_cooling_balls.hpp b/include/volume/volume_cooling_balls.hpp index 29731a2e9..d90604fa4 100644 --- a/include/volume/volume_cooling_balls.hpp +++ b/include/volume/volume_cooling_balls.hpp @@ -694,7 +694,7 @@ template typename RandomNumberGenerator = BoostRandomNumberGenerator > -std::pair volume_cooling_balls(Polytope const& Pin, +std::pair volume_cooling_balls(Polytope& Pin, RandomNumberGenerator &rng, double const& error = 0.1, unsigned int const& walk_length = 1, diff --git a/include/volume/volume_cooling_gaussians.hpp b/include/volume/volume_cooling_gaussians.hpp index ef7ab0f70..ca06ce315 100644 --- a/include/volume/volume_cooling_gaussians.hpp +++ b/include/volume/volume_cooling_gaussians.hpp @@ -57,7 +57,7 @@ struct update_delta> // Compute the first variance a_0 for the starting gaussian template -void get_first_gaussian(Polytope const& P, +void get_first_gaussian(Polytope& P, NT const& frac, NT const& chebychev_radius, NT const& error, @@ -125,7 +125,7 @@ template typename NT, typename RandomNumberGenerator > -NT get_next_gaussian(Polytope const& P, +NT get_next_gaussian(Polytope& P, Point &p, NT const& a, const unsigned int &N, @@ -186,7 +186,7 @@ template typename NT, typename RandomNumberGenerator > -void compute_annealing_schedule(Polytope const& P, +void compute_annealing_schedule(Polytope& P, NT const& ratio, NT const& C, NT const& frac, @@ -289,7 +289,7 @@ template typename RandomNumberGenerator > -double volume_cooling_gaussians(Polytope const& Pin, +double volume_cooling_gaussians(Polytope& Pin, RandomNumberGenerator& rng, double const& error = 0.1, unsigned int const& walk_length = 1) @@ -483,7 +483,7 @@ double volume_cooling_gaussians(Polytope &Pin, { RandomNumberGenerator rng(Pin.dimension()); Pin.set_interior_point(interior_point); - + return volume_cooling_gaussians(Pin, rng, error, walk_length); } diff --git a/include/volume/volume_cooling_hpoly.hpp b/include/volume/volume_cooling_hpoly.hpp index 1175d9b64..c8183384c 100644 --- a/include/volume/volume_cooling_hpoly.hpp +++ b/include/volume/volume_cooling_hpoly.hpp @@ -135,7 +135,7 @@ bool get_next_zonoball(std::vector &HPolySet, return false; } -template +template < typename RandomPointGenerator, typename ZonoHP, @@ -374,7 +374,7 @@ template typename HPolytope, typename Polytope > -double volume_cooling_hpoly(Polytope const& Pin, +double volume_cooling_hpoly(Polytope& Pin, double const& error = 0.1, unsigned int const& walk_length = 1) { diff --git a/include/volume/volume_sequence_of_balls.hpp b/include/volume/volume_sequence_of_balls.hpp index fc1342cee..b79358a6e 100644 --- a/include/volume/volume_sequence_of_balls.hpp +++ b/include/volume/volume_sequence_of_balls.hpp @@ -44,7 +44,7 @@ template typename RandomNumberGenerator > -double volume_sequence_of_balls(Polytope const& Pin, +double volume_sequence_of_balls(Polytope& Pin, RandomNumberGenerator &rng, double const& error = 1.0, unsigned int const& walk_length = 1, @@ -258,7 +258,7 @@ double volume_sequence_of_balls(Polytope &Pin, { RandomNumberGenerator rng(Pin.dimension()); Pin.set_interior_point(interior_point); - + return volume_sequence_of_balls(Pin, rng, error, walk_length, n_threads); } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 076298929..3adbe4681 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -318,6 +318,11 @@ else () COMMAND logconcave_sampling_test -tc=uld) add_test(NAME logconcave_sampling_test_exponential_biomass_sampling COMMAND logconcave_sampling_test -tc=exponential_biomass_sampling) + add_test(NAME logconcave_sampling_test_nuts_hmc_truncated + COMMAND logconcave_sampling_test -tc=benchmark_nuts_hmc_truncated) + add_test(NAME logconcave_sampling_test_nuts_hmc + COMMAND logconcave_sampling_test -tc=benchmark_nuts_hmc) + add_executable (crhmc_sampling_test crhmc_sampling_test.cpp $) diff --git a/test/logconcave_sampling_test.cpp b/test/logconcave_sampling_test.cpp index f8f6901cd..8e53f489b 100644 --- a/test/logconcave_sampling_test.cpp +++ b/test/logconcave_sampling_test.cpp @@ -39,6 +39,8 @@ #include "generators/h_polytopes_generator.h" #include "generators/convex_bodies_generator.h" +#include "diagnostics/univariate_psrf.hpp" + #include "preprocess/svd_rounding.hpp" #include "misc/misc.h" @@ -308,6 +310,80 @@ void benchmark_hmc(bool truncated) { } +template +void benchmark_nuts_hmc(bool truncated) { + typedef Cartesian Kernel; + typedef typename Kernel::Point Point; + typedef std::vector pts; + typedef HPolytope Hpolytope; + typedef BoostRandomNumberGenerator RandomNumberGenerator; + typedef CustomFunctor::GradientFunctor NegativeGradientFunctor; + typedef CustomFunctor::FunctionFunctor NegativeLogprobFunctor; + typedef LeapfrogODESolver Solver; + + typedef Eigen::Matrix MT; + typedef Eigen::Matrix VT; + + NegativeGradientFunctor F; + NegativeLogprobFunctor f; + RandomNumberGenerator rng(1); + unsigned int dim_min = 10; + unsigned int dim_max = 100; + int n_samples = 1000; + long ETA; + bool automatic_burnin = false; + std::chrono::time_point start, stop; + + for (unsigned int dim = dim_min; dim <= dim_max; dim+=10) + { + MT samples(dim, n_samples); + Point x0(dim); + NutsHamiltonianMonteCarloWalk::parameters hmc_params(F, dim); + if (truncated) + { + Hpolytope P = generate_cube(dim, false); + + std::cout << "eta0: " << hmc_params.eta << std::endl; + + NutsHamiltonianMonteCarloWalk::Walk + hmc(&P, x0, F, f, hmc_params, automatic_burnin); + + hmc.burnin(rng); + std::cout << "eta: " << hmc.get_eta_solver() << std::endl; + + start = std::chrono::high_resolution_clock::now(); + for (int i = 0; i < n_samples; i++) { + hmc.apply(rng); + samples.col(i) = hmc.x.getCoefficients(); + } + stop = std::chrono::high_resolution_clock::now(); + std::cout << "proportion of sucussfull steps: " << hmc.get_ratio_acceptance() << std::endl; + } + else + { + NutsHamiltonianMonteCarloWalk::Walk + hmc(NULL, x0, F, f, hmc_params, automatic_burnin); + + hmc.burnin(rng); + std::cout << "eta: " << hmc.get_eta_solver() << std::endl; + + start = std::chrono::high_resolution_clock::now(); + for (int i = 0; i < n_samples; i++) { + hmc.apply(rng); + samples.col(i) = hmc.x.getCoefficients(); + } + stop = std::chrono::high_resolution_clock::now(); + std::cout << "proportion of sucussfull steps: " << hmc.get_ratio_acceptance() << std::endl; + } + + std::cout << "PSRF: " << univariate_psrf(samples).maxCoeff() << std::endl; + + ETA = (long) std::chrono::duration_cast(stop - start).count(); + std::cout<< "time: " << ETA << "\n" << std::endl; + } + +} + template void test_hmc() { typedef Cartesian Kernel; @@ -878,6 +954,11 @@ void call_test_benchmark_hmc(bool truncated) { benchmark_hmc(truncated); } +template +void call_test_benchmark_nuts_hmc(bool truncated) { + benchmark_nuts_hmc(truncated); +} + template void call_test_benchmark_polytopes_grid_search() { typedef Cartesian Kernel; @@ -1196,6 +1277,14 @@ TEST_CASE("exponential_biomass_sampling") { call_test_exp_sampling(); } +TEST_CASE("benchmark_nuts_hmc_truncated") { + call_test_benchmark_nuts_hmc(true); +} + +TEST_CASE("benchmark_nuts_hmc") { + call_test_benchmark_nuts_hmc(false); +} + TEST_CASE("benchmark_hmc") { call_test_benchmark_hmc(false); }