Augustin Chevallier
email: chevallier.augustin@gmail.com

Sylvain Pion
email: sylvain.pion@inria.fr

Frédéric Cazals
email: frederic.cazals@inria.fr

Improved polytope volume calculations based on Hamiltonian Monte Carlo with boundary reflections and sweet arithmetics

Keywords: polytope volume, randomized algorithms, multi-phase Monte Carlo, Hamiltonian Monte Carlo, numerics, multi-precision arithmetic, exact predicates

Computing the volume of a high dimensional polytope is a fundamental problem in geometry, also connected to the calculation of densities of states in statistical physics, and a central building block of such algorithms is the method used to sample a target probability distribution.

This paper studies Hamiltonian Monte Carlo (HMC) with reflections on the boundary of a domain, providing an enhanced alternative to Hit-and-run (HAR) to sample a target distribution restricted to the polytope. We make three contributions. First, we provide a convergence bound, paving the way to more precise mixing time analysis. Second, we present a robust implementation based on multi-precision arithmetic, a mandatory ingredient to guarantee exact predicates and robust constructions. We however allow controlled failures to happen, introducing the Sweeten Exact Geometric Computing (SEGC) paradigm. Third, we use our HMC random walk to perform H-polytope volume calculations, using it as an alternative to HAR within the volume algorithm by Cousins and Vempala. The systematic tests conducted up to dimension n = 100 on the cube, the isotropic and the standard simplex show that HMC significantly outperforms HAR both in terms of accuracy and running time. Additional tests show that calculations may be handled up to dimension n = 500. These tests also establish that multiprecision is mandatory to avoid exits from the polytope.

Introduction

Polytope volume calculations and related problems

Volume calculations. Computing the volume of a polytope -a bounded region of R n defined by the intersection of a fixed set of half spaces (H-polytope) or the convex hull of a finite set of vertices (V-polytope), is a classical problem in science and engineering. Complexitywise, the problem is #-P hard irrespective of the representation of the polytope (H-polytope or V-polytope) [START_REF] Dyer | On the complexity of computing the volume of a polyhedron[END_REF]. This observation naturally calls for approximation algorithms [START_REF] Bárány | Computing the volume is difficult[END_REF][START_REF] Levy | Flavors of Geometry[END_REF] delivering (ε, δ) approximations. Over the years, the complexity of volume calculation algorithms, measured by the number of calls to the oracle stating whether a point is inside the polytope, has been lowered from O (n 23 ) [START_REF] Dyer | A random polynomial-time algorithm for approximating the volume of convex bodies[END_REF] to O (n 4 ) [START_REF] László | Simulated annealing in convex bodies and an O * (n 4 ) volume algorithm[END_REF], and O (n 3 ) [START_REF] Cousins | Bypassing KLS: Gaussian cooling and an O * (n 3 ) volume algorithm[END_REF], the latter for well rounded bodies. More recently, the complexity O (hn 2/3 hn ω-1 ) has been established [START_REF] Tat | Convergence rate of riemannian hamiltonian Monte Carlo and faster polytope volume computation[END_REF], with h the number of hyperplanes and ω the matrix multiplication exponent. The reader is referred to [START_REF] Tat | Convergence rate of riemannian hamiltonian Monte Carlo and faster polytope volume computation[END_REF] for the full history. Interestingly, recent volume calculation algorithms are of the multi-phase Monte-Carlo type, an iterative strategy where each step benefits from the progress made at the previous step. Sketchily, the volume calculation boils down to estimating ratios in a telescoping product, each ratio being the integral over the convex of functions carefully chosen according to a cooling schedule. In recent algorithms, exponential functions are used [START_REF] Cousins | Bypassing KLS: Gaussian cooling and an O * (n 3 ) volume algorithm[END_REF][START_REF] Cousins | Gaussian cooling and O * (n 3 ) algorithms for volume and gaussian volume[END_REF]. We also note that recent work has focused on the reduction of the cooling schedule size, using statistical tests to bound the aforementioned successive ratios [START_REF] Chalkis | Practical volume estimation by a new annealing schedule for cooling convex bodies[END_REF]. The complexity of the algorithm therefore depends on (i) the number of functions in the cooling schedule, (ii) the number of points sampled at each step, and (iii) the complexity of generating a sample according to the specified distribution.

Sampling a target distribution in a bounded domain

Volume calculations require an algorithm to sample a target distribution π in a bounded domain. We describe such algorithms, providing mixing times for the uniform distributionfor the sake of conciseness-from a warm start.

General idea: MCMC. In non trivial cases, the default strategy is to build a Markov Chain leaving the target distribution π invariant. Under mild additional conditions, computing integrals of a function with respect to the target distribution can be approximated by averaging the function on the samples [START_REF] Brooks | Handbook of Markov Chain Monte Carlo[END_REF]. This method is called Markov Chain Monte Carlo (MCMC). The convergence of an MCMC schemed is usually assessed by mixing properties, e.g. based on how fast the Markov chains converges to its stationary distribution [START_REF] Levin | Markov chains and mixing times[END_REF].

Hit-and-run and ball walk. In the case of polytopes, two important Markov chains have been introduced in the literature: Hit and Run (HAR) [START_REF] Lovász | Hit-and-run mixes fast[END_REF][START_REF] Lovász | Hit-and-run is fast and fun[END_REF][START_REF] Lovász | Hit-and-run from a corner[END_REF][START_REF] Haraldsdóttir | CHRR: coordinate hit-and-run with rounding for uniform sampling of constraint-based models[END_REF] and Ball Walk [START_REF] Lovász | Faster mixing via average conductance[END_REF]. These two Markov chains use different strategies to stay inside the polytope: ball walk samples within a ball and rejects points outside of the polytope, while Hit and Run only proposes inside the polytope-no rejection step needed. Furthermore, HAR is amenable to several optimizations [START_REF] Emiris | Efficient random-walk methods for approximating polytope volume[END_REF][START_REF] Emiris | Practical polytope volume approximation[END_REF], including the choice of the random line used, and the calculation of the facet of the polytope intersected by a line. In the context of polytope volume computation, upper bound for the mixing times of the uniform distribution from a warm start exists: HAR mixes in O (n 3 ) steps while ball-walk mixes in O (n 2.5 ) steps [START_REF] Tat | Convergence rate of riemannian hamiltonian Monte Carlo and faster polytope volume computation[END_REF][START_REF] Lovász | Hit-and-run from a corner[END_REF][START_REF] Lee | Eldan's stochastic localization and the KLS hyperplane conjecture: An improved lower bound for expansion[END_REF].

Hamiltonian Monte Carlo. An efficient sampler is Hamiltonian Monte Carlo (HMC) [START_REF] Betancourt | A conceptual introduction to Hamiltonian Monte Carlo[END_REF][START_REF] Radford | MCMC using Hamiltonian Dynamics[END_REF]. HMC relies on the measure preserving properties of Hamiltonian flows on phase space to build a Markov Chain leaving π invariant. In a nutshell, a HMC works by adding a velocity/momentum and an HMC step involves three sub-steps which are (i) picking a random velocity p, (ii) following the Hamiltonian flow associated to H(q, p) = ∇ log π(q) + 1 2 p 2 for a fixed time, and (iii) projecting down in position space. As seen from Hamilton's equation, the fact that the gradient of the target density is used to twist the momentum p rather than the position q helps forcing the dynamical system to glide across the typical set [START_REF] Betancourt | A conceptual introduction to Hamiltonian Monte Carlo[END_REF], which is useful to deal with concentration phenomena. HMC generally uses a Metropolis step when a numerical integration is required for the flow of the Hamiltonian. In our case, we bypass this difficulty by using analytical trajectories. We note in passing that our trajectories have also been derived in the context of sampling truncated multivariate Gaussians [START_REF] Pakman | Exact hamiltonian Monte Carlo for truncated multivariate gaussians[END_REF]. Setting aside the analytical expression of trajectories, our processing is different in two respects: the choice of the travel time (based on ESF and ESS in [START_REF] Pakman | Exact hamiltonian Monte Carlo for truncated multivariate gaussians[END_REF] versus problem dependent for us), and the handling of robustness issues (no handling in [START_REF] Pakman | Exact hamiltonian Monte Carlo for truncated multivariate gaussians[END_REF], versus based on the Sweeten Exact Geometric Computing paradigm for us).

Riemannian HMC. RHMC uses a Hamiltonian exploiting a metric defined on the domain of interest [START_REF] Girolami | Riemann manifold langevin and hamiltonian monte carlo methods[END_REF]. For polytopes, a metric based on the Hessian of the barrier function (φ(x) = -i=1,...,d log(a i T x-b i )) can be used to constrain the random walk (RW) within the polytope [START_REF] Alvarez | Hessian riemannian gradient flows in convex programming[END_REF]. This strategy has been used to sample polytopes and compute their volume [START_REF] Tat | Convergence rate of riemannian hamiltonian Monte Carlo and faster polytope volume computation[END_REF], which yields O (hn 2/3 ) as mixing time in the context of polytope volume computation, with h the number of hyperplanes [START_REF] Tat | Convergence rate of riemannian hamiltonian Monte Carlo and faster polytope volume computation[END_REF]. However, the Hessian of φ is ill-conditioned near barriers [START_REF] Margaret | Some properties of the hessian of the logarithmic barrier function[END_REF], and we are not aware of any implementation for polytope volume calculations.

Dynamical billiard. Random walks in polytopes are also connected to billiards. In two dimensions, it has been proven that there is a G δ dense subset of ergodic billiards in the set of all possible billiards [START_REF] Masur | Chapter 13) Rational billiards and flat structures[END_REF]. For those billiards, Birkhoff theorem implies that almost every trajectory is uniformly distributed on the phase space of the polygon. While analogous properties are unknown for billiards in dimension n > 2 [START_REF] Bálint | Geometry of multi-dimensional dispersing billiards[END_REF], it has been conjectured that billiard trajectories provide valuable building blocks for random walks sampling uniform distributions.

Billiard walk and billiard HMC. A step of billiard walk consists of choosing a direction at random, and following the corresponding billiard trajectory for a fixed time. Properties of billiard trajectories, including their ability to escape from corners, have motivated the sampling algorithm from [START_REF] Polyak | Billiard walk-a new sampling algorithm for control and optimization[END_REF] for general n-dimensional domains, with the velocity refresh ensuring ergodicity. This work is however limited to uniform distributions. In the case of polytopes, billiard walk can be seen as a special case of HMC with reflections on boundaries [START_REF] Afshar | Reflection, refraction, and Hamiltonian Monte Carlo[END_REF]. A general extension to HMC termed reflective HMC is indeed introduced in [START_REF] Afshar | Reflection, refraction, and Hamiltonian Monte Carlo[END_REF], incorporating refraction and reflexions on discontinuities of the potential -respectively corresponding to cases where the discontinuity is passed or not. Since Hamiltonian dynamics can be computed exactly in very few cases, reflective HMC resorts to a leapfrog integrator, a modification preserving the correct stationary distribution. In our case, we prove the uniform ergodicity of billiard HMC (Thm 3, Sect. 2), and experimentally assess its efficiency for polytope volume estimation (Sections 3 and 4).

Robustness issues and the SEGC paradigm

It has long been known that geometric algorithms are prone to (almost) degenerate situations, which manifest even on the simplest expressions in 2D [START_REF] Kettner | Classroom examples of robustness problems in geometric computations[END_REF]. The design of robust geometric algorithms can be done in a general way using the Exact Geometric Computation paradigm [START_REF] Yap | The exact computation paradigm[END_REF], as it is done for example in the CGAL [START_REF] Cgal | [END_REF] software library. In order to do so, the CGAL kernel distinguishes between predicates and constructions. However, the EGC paradigm faces limitations of theoretical and practical nature. The former refer to the impossibility to provably determine the sign (positive, negative or zero) of some real functions in finite time in all cases. The latter relate to the added complexity which can cause the use of resources like processing time and memory space to hit limits.

Without compromising on the robustness properties of the EGC paradigm, we simply allow controlled failures by allowing the sign functions to return a fourth possible value meaning I can't compute. This generalization can practically be implemented using (C++) exception handling mechanisms. We name this extended paradigm Sweeten Exact Computing Paradigm (SEGC). And we refer to sweet arithmetics to qualify the kind of arithmetic functions libraries complying with it. In our implementation, we use the iRRAM C++ library [START_REF] Müller | The iRRAM: Exact arithmetic in C++[END_REF], which provides iterative multiple precision computations for many real functions. This library does not fit into the original EGC paradigm definition, as the zero determination cannot be decided using algebraic separation bounds.

Contributions

This paper makes contributions touching upon random walks for MCMC algorithms, polytope volume calculations, and Hamiltonian Monte Carlo. More precisely:

1. In section 2, we present a robust version of billiard HMC, and prove its uniform ergodicity for convex bodies.

2. In section 3, we instantiate our random walk to sample distributions used for Hpolytope volume calculations. Exploiting analytical expressions for HMC trajectories, we provide a robust implementation of the random walk based on multi-precision interval arithmetic.

3. Finally, section 4 reports experiments comparing HAR and our random walk. The first test samples a target distribution. The second one embeds our random walk into the practical polytope volume calculation of Cousins and Vempala [START_REF] Cousins | A practical volume algorithm[END_REF]. In both cases, we show superior performances over HAR, for dimension up to n = 500.

It is important to note that Riemannian HMC [START_REF] Tat | Convergence rate of riemannian hamiltonian Monte Carlo and faster polytope volume computation[END_REF] forces the random walk to stay within the domain of interest. This is achieved via a distorted metric based on the barrier function, which also yields ill conditioned numerics. We instead use billiard HMC, and control that the RW remains in the polytope using exact predicates based on multiprecision.

Notations. To conform with previous work, the following notations are used in this paper: (i) ε: criterion used to assess the quality of the volume approximation [START_REF] Cousins | A practical volume algorithm[END_REF]; (ii) : notation used for the total variation bound on the mixing time. See theorems 3 and 4. The total variation norm is denoted

• T V .
Proofs and pseudo-code. The reader is referred to the appendix.

Billiard Hamiltonian Monte Carlo

Consider a bounded open set Q ⊂ R n with piecewise smooth boundary and a target probability measure with density π : R n → R ≥0 such that π(x) = 0 for all x not in the closure Q of Q.

Billiard HMC

Denoting q (i) and p (i) the i-th coordinates of position and momentum respectively, recall Hamilton's equations which state that the velocity field in phase space is orthogonal to the gradient of the Hamiltonian H:

dq (i) dt = ∂H ∂p (i) , dp (i) dt = - ∂H ∂q (i) . (1)
As with HMC, we define a potential energy U (q) = -log(π(q)) and the Hamiltonian H(q, p) = U (q) + 1 2 p 2 but this time restricted to Γ = Q × R n . We assume that π is the restriction to Q of positive smooth function defined on R n and we use Φ t the Hamiltonian flow. However, the trajectories of this flow are not included in Q even if the initial point (q, p) is in Q × R n . For every q ∈ Q and every p ∈ R n , we define T (q, p) as the largest T such that for all 0 ≤ t < T , Φ t (q, p) ∈ Q. We also define T (q, p) = 0 when q is in the boundary of Q.

Following [START_REF] Afshar | Reflection, refraction, and Hamiltonian Monte Carlo[END_REF] and [START_REF] Gryazinaa | Random sampling: Billiard walk algorithm[END_REF][START_REF] Polyak | Billiard walk-a new sampling algorithm for control and optimization[END_REF], we modify the flow by forcing reflections on the boundary of Q. This flow, illustrated on Fig. 1, is denoted as follows:

Φt : Q × R n → Q × R n (q, p) → Φt (q, p) ≡ ( Φ(q) t (q, p), Φ(p) t (q, p)). (2)
Note that the latter equation defines position and velocity upon applying the flow. However, as noted in [START_REF] Gryazinaa | Random sampling: Billiard walk algorithm[END_REF][START_REF] Polyak | Billiard walk-a new sampling algorithm for control and optimization[END_REF], this new flow might exhibit problematic trajectories: some of them might not be defined for all t because of singularities on the boundary, others might have huge number of reflections or even an infinite number of reflections in finite time. This does not happen in 2D [START_REF] Chernov | Chaotic billiards[END_REF]Section 2.4 ], but we are not aware of any proof when n > 2 or the trajectories are curved. Hence, we introduce the upper bound M > 0 for the maximum number of reflections, and cull problematic trajectories accordingly.

Remark 1. For q ∈ Q and any p, T (q, p) > 0 and ΦT(q,p) (q, p) is in the boundary of Q.

Remark 2. When q is in the boundary of Q and t > 0, Φt (q, p) is only defined for the momenta p such that the open half-line with origin q and direction p, is included in Q in a neighborhood of q.

Algorithm. Let M ∈ N * be the maximum number of reflections allowed. Given a point q (t) ∈ Q, the algorithm is as follows (Algorithm S1):

Algorithm 1 Hamiltonian Monte Carlo with reflections • (Step 1) Choose the traveling time L ∼ unif(0, 1)

• (Step 2) Pick the momentum p ∼ N (0, I n ) • (
Step 3) If the flow ΦL (q (t) , p) is defined and does not involve more than M reflections between t = 0 and L, and if ΦL (q (t) , p) ∈ Q -Take q (t+1) = Φ(q) L (q (t) , p) -Else, take q (t+1) = q (t)

For a fixed L > 0, steps from 2. and 3. define a Markov kernel P π,L . The full algorithm (steps from 1. to 3.) define a Markov kernel P π that can be expressed with P π,L .

For L > 0, let Γ L be the largest subset of the phase space Γ = Q × R n where ΦL is defined; this set admits no more than M reflections, and the trajectory does not finish in a singularity (a point of the boundary where the normal is not defined) at time L. Γ L is open and therefore measurable. Let

Φ(q, p) = ΦL (q, p) if (q, p) ∈ Γ L (q, p) if (q, p) / ∈ Γ L
the application from Γ to Γ which corresponds to steps 3 and 4.

Remark 3. For any point (q, p) ∈ Γ, if L is small enough, (q, p) ∈ Γ L . However, if L is too large, Γ L could be empty.

-n

H i Φ(q) t c i (q (0) , p) p(t + c i ) p(t - c i ) K H i H j H k Φ(q) t c j (q (0) , p) Φ(q) t c k (q (0) , p) Φ(q) t c l (q (0) , p) = q (1)
q (0)

Figure 1: One HMC step starting at q (0) , with reflections of the HMC trajectory on the boundary of the polytope K. The trajectory successively reflects on hyperplanes, before stopping at q (n L ) . The normal to a facet is denoted n H .

Measure invariance via detailed balance

We recall the definition of detailed balance:

Definition 1 (detailed balance). A Markov chain P is said to satisfy detailed balance (or reversibility) with respect to the measure π if for every A and B measurable subsets,

B P (x, A)π(dx) = A P (x, B)π(dx).
To prove detailed balance for P π , we establish the following intermediate result:

Theorem 1. For a fixed time L, we consider the Markov kernel P π,L associated to steps 2 and 3 of Algorithm S1. Then P π,L satisfies detailed balance with respect to π.

From which one derives the invariance of P π :

Theorem 2. The Markov kernel P π associated to Algorithm S1 satisfies detailed balance with respect to π.

It is well known that satisfying detailed balance implies the invariance of the measure.

Convergence result

Detailed balance ensures that the Markov kernel P π leaves π invariant, but it does not imply the convergence of P n π to π by itself, with P n π is n-times iterated of P π . In this section, we prove uniform ergodicity to get such a convergence result. This requires extra assumptions on Q.

Definition 2 is taken from [START_REF] Roberts | General state space Markov chains and MCMC algorithms[END_REF] with P is a Markov kernel. Practically, we will use P = P π .

Definition 2. A subset C ⊂ X is small (or, (n 0 , , ν)-small) if there exists a positive integer n 0 , a real > 0, and a probability measure ν(.) on X such that the following minorisation condition holds:

P n 0 (x, .) ≥ ν(.) x ∈ C
i.e. P n 0 (x, A) ≥ ν(A) for all x ∈ C and all measurable A ⊂ X Intuitively, the previous definition states that whatever the starting point x-whence the adjective small, the iterated kernels P n 0 (x, .) cover a common measure ν(.). Note that the value of n 0 is the number of steps needed to reach ν -and can be equal to one. The following theorem provides a sufficient condition for Q to be small with respect to P π :

Lemma 1. If Q if convex and the gradient of the potential energy ∇U is bounded on Q, then Q is small for P π , with n 0 = 1.
The ergodicity of P π follows from the above Lemma 1, applied to Theorem 8 from [START_REF] Roberts | General state space Markov chains and MCMC algorithms[END_REF]:

Theorem 3. If Q if convex
and the gradient of the potential energy ∇U is bounded on Q, then P π is uniformly ergodic, that is, for all x ∈ Q:

P n π (x, .) -π(.) T V ≤ (1 -) n (3)
3 Application: computing the volume of a polytope

We now specialize Algorithm S1 and use it as a building block for polytope volume calculation from [START_REF] Cousins | A practical volume algorithm[END_REF]. The general domain Q of Sect. 2 is now the polytope K.

Volume algorithm

In the following, we review the volume calculation algorithm, stressing the role of constants related to the random walk -as these play an important role when the dimension increases.

Overview. Consider a polytope defined in matrix form by Ax ≤ b with the origin 0 inside the polytope. The algorithm used in [START_REF] Cousins | A practical volume algorithm[END_REF] computes the volume of such a polytope with target relative error ε. The principle is a multi-phase Monte Carlo computation, which splits the calculation into m steps. Let {f 0 , . . . , f m-1 } be m isotropic Gaussians i.e.

f i (x) = exp(-x 2 /(2σ 2 i )) with σ i = 1/ √ 2a i , or equivalently f i (x) = exp(-a i x 2 ), (4)
such that the first one is highly concentrated around a point deep inside the convex, and the last one is an almost flat distribution. The volume calculation reduces to computing the telescoping product

Vol(K) = K f 0 (x)dx K f 1 (x)dx K f 0 (x)dx . . . K dx K f m-1 (x)dx ≡ K f 0 (x)dx i=1,...,m R i (5)
Each of these ratio R i are estimated using Monte-Carlo integration with respect to the density

π i-1 (x) = f i-1 (x) K f i-1 (y)dy 1 K (x) (6)
via importance sampling:

R i = K f i (x)dx K f i-1 (x)dx = K f i (x) f i-1 (x) π i-1 (x)dx. (7)
The method then uses as core block an algorithm sampling the previous distribution, usually using HAR. In our case, HMC with reflections will be used instead.

Starting Gaussian f 0 (x). The constants involved in the choice of f 0 are not related to the random walki.e. this step does not use the random walk.

Cooling schedule incremental construction. The sequence of Gaussians f i (x) is chosen so as to estimate R i efficiently. One takes a i = a i-1 (1 -1/n) r , trying to maximize r so as to reach the uniform distribution as rapidly as possible. In practice, this is done by running a dichotomy on the values of r -one knows that r = 1 suffices, controlling the variance of the RV Y = exp((a i-1 -a i ))X). The following numbers of samples are used to estimate that the condition

Var [Y ] /E[Y ] 2 ≤ C holds:
• HAR: N cs = 500 * C + n 2 /2 (From [START_REF] Cousins | A practical volume algorithm[END_REF], with C = 2).

• HMC: N cs = 500 * C + n (Our algorithm, also with C = 2).

Cooling schedule termination. To decide whether the cooling schedule needs another Gaussian or the value m has been reached, the ratio R i (Eq. 7) is evaluated with a number N m of samples. If this ratio of close to one (≤ 1.001 in the MATLAB implementation of [START_REF] Cousins | A practical volume algorithm[END_REF]), the cooling schedule is complete. In practice, this evaluation is carried out using N m = 150/ε samples returned by the random walk (MATLAB implementation of [START_REF] Cousins | A practical volume algorithm[END_REF]).

Convergence of ratios R i . For X 1 , ..., X k consecutive points given by the random walk

sampling π i-1 , the Monte-Carlo estimation R (k) i of R i is given by: R (k) i = 1 k k j=1 f i (X j ) f i-1 (X j ) . (8)
The following stop criterion is introduced in [START_REF] Cousins | A practical volume algorithm[END_REF]. Let ε = ε/ √ m; this is the relative ratio error allocated for each ratio R i estimation. Consider a sliding window of size W . When W consecutive estimated ratios R

(k-W +1) i , ..., R (k) i are within ε /2, the convergence for R i is declared.
In this work, we test the following window sizes:

• HAR with window size W HAR = 500 + 4 * n 2 (From [START_REF] Cousins | A practical volume algorithm[END_REF]).

• HMC0: HMC with window size W HM C0 = 250

• HMC1.5: HMC with window size W HM C2 = 250 + n √ n

HMC algorithm

Our HMC specialization has analytical trajectories -see also [START_REF] Pakman | Exact hamiltonian Monte Carlo for truncated multivariate gaussians[END_REF], whose intersection with the polytope have simple expressions.

Analytical trajectories. We build an HMC sampler for π i (x) from Eq. 6. Let U (q) = log(exp(-a i ||q|| 2 )) = -a i ||q|| 2 . U can be interpret as a potential energy whose associated Boltzmann measure is π i . The associated Hamiltonian in phase space is given by H(q, p) = U (q)+1/2 p 2 . Note that the normalization constant 1 K f i (y)dy was discarded because it does not change the trajectory (its gradient is 0). Rewriting the dynamical system associated to this Hamiltonian yields the following differential equations:

d 2 q j dt 2 (t) = -2a i q j (t) for j ≤ n. ( 9
)
Each coordinate is independent and has a solution of the form

q j (t) = C j cos(ωt + φ j ), (10)
with C j , ω, φ j ∈ R. The parameter Φ j satisfies the following 2 equations:

cos(φ j ) = q j (0) C j , sin(φ j ) = -p j (0) ωC j . ( 11
)
Thus we deduce:

       ω = √ 2a i C j = q j (0) 2 + p j (0) 2 /ω 2 φ j = arctan - p j (0) q j (0)ω + 1 {q j (0)<0} π (12)
Note that these equations are the same for any choice of coordinates as long as the basis is orthonormal-an observation exploited below.

Collision with convex boundary. To restrict the sampling algorithm to the polytope K, trajectories should reflect on the boundary of K. The polytope K is defined by a set of hyperplanes. Hence, we need to compute the intersection time of a trajectory with each hyperplane and take the smallest time. Thankfully, there is an analytical expression. The hyperplanes are defined by a matrix A and a vector b, and hyperplane i is defined by the equation (Ax) i = b i .

To compute the collision time with a hyperplane H, we make the following remark: let n H be the normal of the hyperplane. We can complete n H to an orthonormal basis. In this basis, the collision time depends only on what happens for the coordinate on n H . Consider q n H (t) =< q(t), n h > and p n H (t) =< p(t), n H > the coordinates along the normal of the hyperplane. Let ω,C n H and φ n H be the parameters of the trajectory along direction n H . Finding the intersection times is equivalent to solving the equation for t:

C n H cos(ωt + φ n H ) = b i (13)
We deduce that if |C n H | < b i , there is no solution, and else, the following times are solution:

t 1 = (arccos(b i /C n H ) -φ n H ) /ω. t 2 = (-arccos(b i /C n H ) -φ n H ) /ω. ( 14
)
One solution corresponds to the entry into K, while the other corresponds to the exit out of K. We select the exit trajectory via a dot product between the velocity at t 1 and t 2 and the outward normal n H . In the sequel, the corresponding value is denoted t c for time of collision.

Travel time choice with respect to a i

Algorithm S1 can be slightly generalized by choosing a travel time L uniformly in [0, L max ] instead of [0, 1]. We propose here a strategy to chose L max with respect to the parameter a i of the Gaussian sampled (σ i = 1/ √ 2a i ). We distinguish two main cases here. First, if a i is large, then the probability measure is concentrated in a neighborhood of 0, and the trajectories will seldom hit the boundaries and the strategy is largely unaffected by the hyperplanes bounding the polytope. Second, if a i is close to 0, the distribution is close the uniform distribution of the convex.

Case 1.

If a i is large, we consider q j (t) a solution of the dynamical system 9 and we introduce space-time rescaling for the trajectory:

qj (t) = √ a i q j (t/ √ a i ).
It satisfies Eq. 9 with a i = 1:

d 2 qj dt 2 (t) = qj (t) with initial condition dq j dt (0) = dq j dt (0) ∼ N (0, 1). (15)
It follows that in the absence of convex boundaries, the rescaled process q is sampling the distribution associated to a i = 1, i.e. the standard normal distribution. Therefore, any sensible value for L max (1) taken for a i = 1 should be scaled as

L max (a i ) = L max (1)/ √ a i .
Case 2. When a i goes to 0, the previous strategy leads L a i max → ∞. We argue that in this case, the trajectories converges to billiard trajectories with reflections on the boundary and that the optimal L max should be therefore close to the optimal L 0 max for a i = 0. To the best of our knowledge, this optimal value is not known. However as we shall see in section 3.5, if the polytope is a cube, choosing L 0 max such that the number of reflection on the boundary is proportional to the dimension leads to a O(log(n)) mixing time (where n is the dimension). Hence we propose a heuristic where for a given parameter r, L 0 max is chosen such that on average for a i = 0 the trajectory hits the boundary r × n times. Practically, we perform a binary search on the values of L 0 max , and for each candidate value evaluate the average number of reflections over N 0 (= 1000) steps.

Values used in experiments. Mixing the two cases, we chose L max as follows:

L max (a i ) = min(π/ √ 2a i , L 0 max ) (16)
Remark 4. As already noticed, our trajectories are identical to those used in [START_REF] Pakman | Exact hamiltonian Monte Carlo for truncated multivariate gaussians[END_REF]. In this latter work, heuristics based on the Effective Sample Factor and the Effective Sample Size are presented to set the travel time T . But a fixed value (T = π/2) is finally used.

HMC implementation based on interval arithmetic

Robustness issues

The algorithm as stated before can be implemented in the real RAM model (Algo. S1), but is prone to numerical rounding errors. As a particular case, one may consider the situation where rounding errors would be such that the point computed on the HMC trajectory would be outside the convex. More generally, all geometric constructions and geometric predicates on them potentially raise robustness issues -see list in section 3.4.3.

As discussed in Introduction, we guarantee robustness using the SEGC paradigm. More specifically, recall that efficient arithmetic operations usually combine two ingredients: first, an interval representation of the numbers, as non overlapping intervals yield exact predicates; second, an arbitrary precision representation of the interval bounds, as precision can be increased so as to yield exact predicates and constructions of controlled accuracy. In the sequel, we use the iRRAM library which provides these two ingredients [START_REF] Müller | The iRRAM: Exact arithmetic in C++[END_REF].

For the sake of clarity, all functions for which the iRRAM library plays a key role are highlighted in blue.

iRRAM and used features

We represent points as n-dimensional points whose coordinates are of the iRRAM number type. In iRRAM, a real number is represented by two types of data: firstly a symbolic representation memorizing the way it was defined (type of function and pointers to the operands), and secondly a numeric approximation using an interval with rational endpoints guaranteed to enclose the exact real value. The accuracy of the latter interval can be increased if needed, by recomputing it using recursively increased precision of the operands intervals that defines it. Therefore, the iRRAM encoding q iRRAM (t) of the position q(t) of the HMC trajectory is numerically represented by a n-dimensional box certified to contain the exact real position.

Two specific operations of iRRAM may trigger numerical refinement:

• operator x < y: predicate answering the < comparison operator. If the interval representations of x and y overlap, these intervals are automatically refined, a feature called precision refinement thereafter.

• Near_inf_double(iRRAM x) : returns the nearest double < the iRRAM number x.

Robust operations and robust HMC step

Our robust implementation calls for the following operations (i) Computing the trajectory parameters -Eq. 12, (ii) Finding the exit intersection time -Eq. 14, (iii) Finding the smallest exit intersection time amongst all hyperplanes, and (iv) Evolving the trajectory. In the sequel, we detail these operations, and refer the reader to the SI section 7; in particular, Algorithm S2 refines Algorithm S1 based on these robust primitives.

Trajectory parameters. Following Eq. 12, we construct the numbers C j , w and φ j as iRRAM number types. See Algorithm S4.

Exit intersection time t c with one hyperplane.

Intersecting the trajectory with a hyperplane yields two solutions (Eq. 14) respectively exiting and entering the convex. The collision time t c corresponds to the former. Assuming that the normal vector to the hyperplane is oriented outwards, the value t c is such that < p(t c ), n H >> 0. The evaluation of this predicate triggers a precision refinement if needed.

Smallest exit intersection time. Since the boundary of K involves several hyperplanes, the nearest one, which corresponds to the smallest exit time, must be determined. To do so, we first construct the intersection time t c i with respect to each hyperplane. Then we compute t c = min i t c i .

This calculation is tantamount to sorting the individual intersection times, which in turn requires the comparison operator <. iRRAM provides such an operator, which triggers precision refinement if needed. See Algorithm S5.

Remark 5. We note that in case the trajectory would hit a face of dimension < d -1, an equality between exit times occurs. The absence of separation bound in iRRAM does not allow us to handle such cases, and an infinite refinement loop is entered. However, for each starting point, the measure of velocities leading to such sets is null. Practically, such a case was never faced-as expected.

The Is_strictly_in_convex[N T ](P oint q) predicate. To constrain the trajectory within the convex, we resort to a predicate using a number type NT, telling whether a given position q belongs to the interior K o of K. The predicate checks whether < op, n > < b i holds for every hyperplane. (Nb: < •, • > stands for the dot product of two vectors.) The robust implementation Is_strictly_in_convex[iRRAM](iRRAM_point_d p) of this predicate uses iRRAM as number type, triggering the iRRAM refinement if needed.

Robust implementation of one HMC step. The robust implementation of one step hinges on two operations (Algorithm S2):

• Computing the nearest inferior double t < c = Near_inf_double(t c ). The intersection point between the trajectory and a hyperplane is defined analytically by Eq. 14. The iRRAM representation of the collision time is an interval certified to contain the exact solution, and the corresponding n-dimensional point q iRRAM (t < c ) is represented as a box. We note that the box q iRRAM (t < c ) intersects the interior K o of the polytope K. Indeed:

the exact collision point q(t c ) lies on its defining hyperplane i.e. q(t c ) ∈ H i by definition of t < c , the exact embedding q(t < c ) satisfies

q(t < c ) ∈ K o -the iRRAM box q iRRAM (t < c ) corresponding to t < c intersects K o since q iRRAM (t < c ) q(t < c ) ∈ K o (17)
• Calling the predicate Is_strictly_in_convex[iRRAM](q(t < c )). Recall that in iRRAM, a n-dimensional point is represented as a box. This is in particular the case for the collision points with the hyperplanes, and for the final point returned. To ensure that all such points are strictly within the convex, we call the aforementioned Is_strictly_in_convex[iRRAM].

Mixed precision and MATLAB interface

In practice, a robust implementation faces two issues: it is much slower than the corresponding double precision implementation, and it currently has to be interfaced with MATLAB, which cannot handle iRRAM number types.

To solve the first problem, we use a lazy two step strategy (Algorithm S3). First, the HMC step is computed using double precision. Then the programs checks whether the resulting point is inside the polytope K, using the CGAL library. If not, the computation is restarted using the robust iRRAM implementation.

The second issue forces us to convert the resulting point from a representation using the iRRAM number type (when relevant) to a representation using doubles. This conversion may result in a point outside of the convex. We check this property using the Is_strictly_in_convex predicate, using here a lazy approach again: first, Is_strictly_in_convex[CGAL::Interval_nt_advanced] is computed; in case of failure, Is_strictly_in_convex[CGAL::Lazy_exact_nt<CGAL::Quotient<CGAL::MP_Float> >] is evaluated [START_REF] Fabri | A generic lazy evaluation scheme for exact geometric computations[END_REF]. When the point is detected outside the convex, we refuse it, resulting in a move step keeping the initial point invariant. Remark 6. Having started strictly inside the convex, with respect to the SEGC paradigm (Sec. 1.3), we never faced a situation I can't compute.

Cube: HMC mixing time is O(log n)

We make a small digression in the case of the cube. It turns out that in this special case, the mixing time of the HMC random walk is O(log n) with n the dimension, and the average complexity per step (the average number of calls to the oracle per step) is linear with the dimension. We assume without any loss of generality that the cube is [0, 1] n . Rigorously:

Theorem 4. Let P (n) k
the distribution after k steps in dimension n and g n an isotropic Gaussian of parameter a i restricted to [0, 1] n . Then there exists 0 < ρ < 1 such that for every > 0, every n > 1 and every x ∈ R n , we have for k ≥ (log n -log )/(log 1/ρ):

P (n) k (x, .) -g n T V ≤ . ( 18
)
Furthermore, the average number of reflections per step is O(n).

Experiments

Implementation and code availability

Implementation.

The implementation of our algorithm for the particular case of a convex polytope Q = K is provided in the Structural Bioinformatics Library (SBL, http: //sbl.inria.fr), a state-of-the-art environment targeting molecular simulation at large [START_REF] Cazals | The Structural Bioinformatics Library: modeling in biomolecular science and beyond[END_REF]. The corresponding package, Hamiltonian_Monte_Carlo is ascribed to the Core / Geometry-Topology component of the library. The user manual of the package, as well as a jupyter notebook illustrating the main functionalities, can be accessed from https: //sbl.inria.fr/doc/Hamiltonian_Monte_Carlo-user-manual.html. Volume calculations. As described in section 3, we also embed our random walk in the framework of [START_REF] Cousins | A practical volume algorithm[END_REF]. We reuse the MATLAB code provided by [START_REF] Cousins | A practical volume algorithm[END_REF] and adapt it so as to call our HMC random walk instead of the usual HAR random walk. Also following [START_REF] Cousins | A practical volume algorithm[END_REF], we use HAR (rather than ball walk) as contender.

Illustrations of the HMC random walk

Our first illustration features samples generated by HMC for the cube [-1, 1] 3 and the standard simplex

x i ≤ 1 in dimension three. Varying the parameter σ rapidly yields concentrated samples (Fig. 2). Our second illustration shows the ability of the algorithm to escape corners. As opposed to HAR, HMC yields an almost uniform distribution after 10 steps even when the dimension increases (Fig. 3).

Experimental setup

Models. In [START_REF] Cousins | A practical volume algorithm[END_REF], five types of polytopes are used for testing: Birhkoff polytopes, triplystochastic polytopes, the cube, the isotropic simplex, and the standard simplex. To present detailed statistics on the quality of the output, we focus on the latter three, since their exact volume can be computed. We perform a systematic study in dimension n = 10, 30, 50, 70, 100. As target error, we use ε = 0.01 -see Section 3.1. We also present selected results in dimension n = 250 and n = 500, with ε = 0.05.

Algorithms and parameters used. As specified in Section 3.1, volume algorithms based on HAR and HMC require three sets of constants: N cs for the cooling schedule incremental construction, N m for the cooling schedule termination, and the window size W to assess the convergence of ratios R i . In addition, fixing the travel time in HMC requires the parameter r (see the discussion accompanying Eq. 16).

As specified in Section 3.1, we compare three contenders HAR, HMC0, and HMC1.5.Table 1 summarizes all parameters.

To perform statistics, 50 runs of each algorithm on a given polytope of prescribed dimension were obtained. Summarizing, our experiments involve 3 algorithms × 3 polytopes × 5 dimensions × 50 repeats -whence a total of 2250 runs. A nearly flat isotropic Gaussian was used to approach the uniform distribution (σ = √ 500). For n > 3, the plot displays projections onto the first two coordinates. Simulations were started from a corner (q (0) i = 0.9). Samples generated after 10 steps of HAR or HMC, repeated 500 times -whence 500 samples. Remark 7. Our implementation do not handle the degenerated case a f lat = 0, hence we use a f lat = 10 -13 (Eq. 4) to tune the travel time for the flat Gaussian.

N cs N m W r HAR 1000 + n 2 /
Statistics. We collect the following statistics:

• the relative error of the estimated Ṽ :

err r =| Ṽ -Vol(K) | /Vol(K). (19)
• #S: number of sampled points for a volume computation.

• #O: complexity, i.e. the number of calls to the oracle. For HAR, this is equal to the number of sampled points. For HMC, #O takes the number of reflections into account.

• #R: multi-precision refinements triggered in iRRAM (Algo. S3). Nb: MATLAB first attempts a calculation using double floating point numbers. In case of failure, a precision refinement is triggered.

• #E: the number of exits of the polytope (Algo. S3). As discussed in Section 3.4.4, this may happen when converting the coordinates of a point from the iRRAM representation to the double representation. Note that such events are not related to our implementation, but instead to the number type (doubles) used in the MATLAB implementation of the volume algorithm.

• t: the running time, t in seconds.

Computer. Calculations were run on a desktop DELL Precision 7920 Tower (Intel Xeon Silver 4214 CPU at 2.20GHz, 64 Go of RAM), under Linux Fedora core 32.

Results

We first analyze results up to dimension n = 100.

Quality of estimates vs running time. Analyzing the relative error err r and the running time t per se, as expected, a clear correlation between the median values and the stdev appears (Fig. S1, Table S1). In terms of comparisons, the global trend states that HMC0 is faster than HMC1.5, which is faster than HAR, a fact especially visible for large values of the dimension, whatever the polytope.

Next, we use scatter plots to relate the relative error (median, stdev) to the running time.

The first two scatter plots link symbols for a given polytope/dimension, which is instrumental to compare the contenders (Fig. 4(Tow row)). We first notice that compared to the simplices, the cube appears as an easy case. For the two simplices, HMC is faster and more accurate than HAR; also, HMC1.5 is more accurate but slower than HMC0. As a specific illustration, consider the standard simplex in dimension n = 100: HMC0 is ∼ 16 times faster than HAR, achieving a median accuracy 7.34% against14.9% for HAR. HMC1.5 gets an accuracy of 4.7% while remaining ∼ 5.7 times faster than HAR.

The second two scatter plots link symbols for a given algorithm, which is useful to assess the incidence of dimension (Fig. 4(Bottom row)). For HMC, the median and the stdev errors remain in a narrow band, that is, these order one and two moments do not increase as a function of the dimension. We also note that statistics for HMC1.5 are more stable than those of HMC0. For HAR, one needs to distinguish cases. As already noticed, the cube is an easy case, for which a behavior similar to HMC is observed. For the two simplices, both the median and the stdev of the error significantly increases beyond the dimension n = 50.

From these observations, depending on the desired trade-off between accuracy an speed, one may choose HMC0 or HMC1.5.

Other statistics -robustness. The main issue with trajectories computed with doubles is that they may exit the polytope and remain outside: this is observed in practice, with selected runs using double floating point numbers never ending. Our robust implementation fixes this difficulty. As a quantitative measure, we track the number of times #R precision refinement is triggered in iRRAM-see Section 4.3. It clearly appears that multiprecision gets used more often when the dimension increases, say in the range 50-100 for our experiments (Table S2). This is particularly striking for the standard simplex, with a number of calls in thousands in dimension n = 100. We speculate that this is due to local geometry properties near corners of the simplex. For HMC, another interesting statistic is #O/#S, namely the number of reflections per step. As expected, the ratio between the number of reflections and the number of step is a bit smaller than the target r × n (Section 3.3), since for most of the Gaussian the length of the trajectory will be smaller than L 0 max .

Results in dimension n = 250, 500. We now briefly comment on results for higher dimensional cases (Tables S3, S4, S5, S6). For the dimensions tested, only the cube is handled by the implementation of the volume algorithm, due to the format of the floating point numbers used in MATLAB. (Note that the exact volume is not computed either due to this issue.) Using ε = 0.05, we notice that HMC manages to estimate the order of magnitude of the volume of the cube satisfactorily. As for robustness, we note that precision refinements are called often, especially for the standard simplex, an observation already raised for n = 100.

Conclusion

To compute the volume of high dimensional polytopes, this paper exploits a novel strategy based on billiard Hamiltonian Monte Carlo. This strategy exploits the ability of billiard walks to escape corners, and the simple analytical expression of the Hamiltonian makes a robust but fallible implementation based on our SEGC paradigm possible and effective, as shown by experiments on polytope volume calculations using the iRRAM and CGAL libraries.

On the theoretical side, our work leaves stimulating questions open, two of which are of prominent importance. The first one is the choice of the optimal travel time required to sample a convex uniformly, which would ideally be determined at runtime for each convex. The second one is the analysis of the incidence of reflections on the mixing time, so as to quantify the speed at which reflections decorrelate successive points.

Our work also prompts a connection with statistical physics, for the calculation of density of states (DoS), which measure the volume in phase space of the pre-image of an energy stratum. DoS are central to compute partition functions, whence all thermodynamic quantities. We anticipate that billiard on level set surfaces bounding strata will prove beneficial, with a potential impact in statistical physics at large.

Proofs

Lemmas for Thm. 2

Let A and B be measurable subsets of Q. Then let

A B = {(q, p) ∈ Γ L |q ∈ A, Φ(q, p) ∈ B × R n }
be the subset of Γ of all positions in A with momenta that brings them in B after time L. Similarly, we define on Γ Ψ(q, p) = ( Φ(q) (q, p), -Φ(p) (q, p)).

(

) 20
Lemma 2. The maps Φ and Ψ preserve the Lebesgue measure on Γ L . ie for every

A ⊂ Γ L measurable, λ( Φ-1 (A)) = λ(A) Lemma 2.
Clearly it is enough to prove that Φt preserves the Lebesgue measure. A numerical flow based on Leapfrog integration is introduced in [START_REF] Afshar | Reflection, refraction, and Hamiltonian Monte Carlo[END_REF], that we call Φh t for a stepsize h. For a fixed t, letting the step size h going to zero, we see that Φt is the pointwise limit of Φh t . Let A be a bounded open subset. Clearly:

λ( Φ-1 t (A)) = 1 A ( Φt (x))dx and λ(( Φh t ) -1 (A)) = 1 A ( Φh t (x))dx,
where 1 A is the indicator function of the set A. Furthermore, theorem 1 of [START_REF] Afshar | Reflection, refraction, and Hamiltonian Monte Carlo[END_REF] shows that the numerical flow Φh t preserves the Lebesgue measure:

λ(A) = λ(( Φh t ) -1 (A)).
Let h k > 0 be a sequence of step sizes that converges to 0.

Let B = Φ-1 t (A)∪ ∪ k ( Φh k t ) -1 (A)
. Since A is bounded, and using properties of the numerical flow, B is also bounded. Furthermore, since A is an open set, 1 A • Φt is also the pointwise limit of 1 A • Φh k t when k goes to +∞. Using the dominated convergence theorem with dominating function 1 B , we deduce that λ( Φ-1 (A)) = λ(A) for all bounded open subset, which implies that it is true for all measurable subsets. Lemma 3.

1.

Ψ(Γ L ) ⊂ Γ L 2. Ψ • Ψ = I on Γ L 3.
For any measurable sets A and B of R n ,

B A = Ψ(A B )
4. H(Ψ(q, p)) = H(q, p) for all (q, p) ∈ Γ. Thm. 1. Let A and B be measurable sets of R n . The left hand side of the detailed balance equation becomes A P π,L (q, B)dπ(q) = A P π,L (q, B) exp(-U (q))dq = A {p|(q,p)∈Γ L , Φ(q) (q,p)∈B} exp(-p 2 )dp + {p|(q,p) / ∈Γ L , Φ(q) (q,p)∈B} exp(-p 2 )dp exp(-U (q))dq

= A B exp(-H(q, p))dqdp + A {p|(q,p) / ∈Γ L } 1 B (q) exp(-H(q, p))dqdp = A B
exp(-H(q, p))dqdp + A∩B {p|(q,p) / ∈Γ L } exp(-H(q, p))dqdp.

By symmetry:

B P π,L (q, A)dπ(q) = B A exp(-H(q, p))dqdp + A∩B {p|(q,p) / ∈Γ L } exp(-H(q, p))dqdp

Using lemma 3 and the measure conservation of lemma 2 we obtain:

A B exp(-H(q, p))dqdp = Ψ(B A )
exp(-H(q, p))dqdp

= B A exp(-H(Ψ(q, p)))dqdp = B A exp(-H(q, p))dqdp
Which concludes the proof.

Proof of Lemma 1

Lemma 4. Let Q be an open convex subset in R n that contains the ball B(0, 2ρ), let x be a point in Q and let v be in B(0, ρ) -x. Let f : R n → R n be a continuous and bounded map.

Suppose that (x(t), v(t)) ∈ R 2n is the solution of the Cauchy problem

x(0) = x v(0) = v x (t) = v(t) v (t) = af (t)
where a is a real number. If |a| ≤ ρ f ∞ then for all t ∈ [0, 1], x(t) is in the convex hull of x and of the ball B(0, 2ρ) and therefore in Q.

Lemma. 4. Suppose |a| ≤ ρ f ∞ . By the mean value theorem, v(t) = v +w(t) where w(t) ≤ |a| f ∞ t ≤ ρt for t ≥ 0. The derivative of function y(t) = x(t) -vt is w(t), therefore by the mean value theorem, y(t) = y(0) + tz(t) where z(t) ≤ ρt/2. Therefore for all t ∈ [0, 1],

x(t) = y(t) + vt = (1 -t)x + t((v + x) + z(t))
is in the convex hull of x and of the ball B(0, 2ρ).

f (t) = dϕ(x + tu)(u) • u = Id(u) • u + (dϕ(x + tu) -Id)(u) • u ≥ 1 -dϕ(x + tu) -Id u u > 0. It follows that f ( y -x ) > f (0). Now f (0) = ϕ(x) • u and f ( y -x ) = ϕ(y) • u, hence ϕ(x) = ϕ(y).
C i \ E is open. Now C i is connected, hence E or C i \ E is empty. Therefore E = C i .
Since ϕ(B) is compact, E is bounded, and therefore E = C i = C 1 . Moreover, by assumption, ϕ(0) ∈ B o (0, r/2) ∩ E which implies that B o (0, r/2) ⊂ E. Lemma 1. We assume without any loss of generality that 0 ∈ Q.

Let q (1) , q (2) ∈ Q. One way to get a trajectory going from the point q (1) to a point very close to q (2) , is to select a very high momentum p α = 1 α (q (2) -q (1) ) and a short time t α = α with α > 0. When α → 0, the potential energy term U of the Hamiltonian becomes less and less relevant, thus the trajectory converges to a straight line and lim α→0 Φ q tα (q (1) , p α ) = q (2) . This intuition is formalized using the scaled variables q(t) = q(αt) p(t) = αp(αt).

(

) 21
The equation of motion becomes:

dq dt (t) = p(t) (22)
dp dt (t) = α 2 ∇ q U (q(t)). (23)
which defines a flow φ(α, t, q, p). It should be noted that φ(-α, t, q, p) = φ(α, t, q, p) for every α, t, q, p and that φ is correctly defined for α = 0. In this case, it is easy to see that: φ(0, t, q, p) = q + pt.

As π is the restriction of a positive smooth function, φ is defined for every (α, t, q, p) ∈

[-1, 1] × R + × R n × R n .
Furthermore, α 2 ∇ q U is smooth. Hence, using the differentiability of the solutions of differential equations on parameters and initial conditions (see [START_REF] Burke | Continuity and differentiability of solutions[END_REF]), we see that φ is C

2 on ] -1, 1[×Q × R + × R n .
Since Q is open, there exists ρ > 0 such that B(0, 2ρ) ⊂ Q. Let ν be the measure on Q which is the Lebesgue measure on B(0, ρ/2) and zero outside the ball B(0, ρ/2).

Our aim is to show that there exists > 0 such that for every q ∈ Q, P π (q, .) ≥ ν(.). Note that we would only need that P n 0 π (q, .) ≥ ν(.) for some n 0 , but since Q is convex we will be able to take n 0 = 1.

The density of the probability measure associated with momenta is exp(-1/2 p 2 ), and taking into account the rescaling of the momenta, we define the probability density

γ(p) = α exp(- 1 2 1 α p 2 ). (25)
Q is bounded, hence there exists > 0 such that for every (q, p) ∈ {(q, p)|q ∈ Q, p ∈ B(-q, ρ)}, γ(p) >

Mixing time for the cube

Lemma 4. We consider the canonical basis of R n . As shown before in Eq. ( 10), the trajectory coordinates associated to the Gaussian are all independent from each other. Furthermore, when a reflection with a boundary occurs, it means that one of the coordinates reached 0 or 1. The reflection simply switches the sign of the momentum for this coordinate, leaving other coordinates unchanged. Finally, the initial momentum vector is sampled according to p(0) ∼ N (0, I n ), therefore each coordinate of p(0) is sampled from an independent Gaussian N (0, 1) in R.

Hence we conclude that each coordinate has the behavior of a 1-dimensional HMC random walk sampling a 1-dimensional Gaussian, all independent from each other.

Let us consider the 1-dimensional random walk for a given Gaussian. We write P k (x, .) the distribution after k steps starting from x ∈ [0, 1], and g the probability density associated with the restriction of the Gaussian to [0, 1]. Using theorem 3 combined with theorem 1 for the 1-D Gaussian restricted to [0, 1], we deduce that there exists 0 < ρ < 1 such that for all x ∈ [0, 1],

P k (x, .) -g T V ≤ ρ k .
Observe that writing P (n) k the distribution after k steps in dimension n and g n the Gaussian restricted to [0, 1] n , we have

P (n) k (x, y) = P k (x 1 , y 1 )P (n-1) k ((x 2 , ..., x n ), (y 2 , ..., y n )), g n (x) = g(x 1 )g n-1 ((x 2 , ..., x n )) (33)
The total variation distance can be written as

P (n) k (x, .) -g n T V = 1 2 [0,1] n P (n) (x, y) -g n (y) dy = 1 2 [0,1]×[0,1] n-1 P k (x, y 1 )P (n-1) k (x, y 2 ) -g(y 1 )g n-1 (y 2 ) dy 1 dy 2 = 1 2 [0,1]×[0,1] n-1 (P k (x, y 1 ) -g(y 1 ) + g(y 1 )) P (n-1) k (x, y 2 ) -g(y 1 )g n-1 (y 2 ) dy 1 dy 2 ≤ 1 2 [0,1]×[0,1] n-1 (P k (x, y 1 ) -g(y 1 )) P (n-1) k (x, y 2 ) dy 1 dy 2 + 1 2 [0,1]×[0,1] n-1 g(y 1 )P (n-1) k (x, y 2 ) -g(y 1 )g n-1 (y 2 ) dy 1 dy 2 ≤ 1 2 [0,1] |P k (x, y 1 ) -g(y 1 )| dy 1 + 1 2 [0,1] n-1 P (n-1) k (x, y 2 ) -g n-1 (y 2 ) dy 2
We deduce

P (n) k (x, .) -g n T V ≤ n P k (x, .) -g T V ≤ nρ k (34)
Hence for a fixed , if k satisfies nρ k ≤ , then for every x, P

k (x, .) -g n T V ≤ . Thus we take k ≥ (log n -log )/(log 1/ρ), and the mixing time is O(log n).

In addition, as each coordinate is from each other, the total number of reflections is the sum of reflections per coordinate. Hence, the number of reflections is proportional to the dimension.

Supporting information: pseudo-code

In the following, we provide the pseudo-code for the high level description of the HMC algorithm (Algorithm 1).

Algorithm S 1 Hamiltonian Monte Carlo step with real RAM arithmetic model 1: HMC_step(q) 2: Choose a travel time L ∼ U nif (0, L max ). (intersection, t c ) ← Intersect_hyper_planes(q, p) // find intersection with hyperplanes 7:

if intersection = False OR dist < t c then 8:

(q, p) = Update_positions_momenta(q, p, dist) // update traj. with distance dist

9:

Set dist = 0 10: else 11:

(q, p) = Update_positions_momenta(q, p, t c )

12:

Reflext_normal(p(t c ), n c )

13: Set dist = dist -t c
Algorithm S 2 Hamiltonian Monte Carlo step with iRRAM number type 1: HMC_step(q) 2: Choose a travel time L ∼ U nif (0, L max ). (intersection, t c ) ← Intersect_hyper_planes(q, p) with t c an iRRAM REAL.

7:

t < c = Near_inf_double(t c ) 8:
if intersection = False OR dist < t < c then 9:

(q, p) = Update_positions_momenta(q, p, dist) // update trajectory with distance dist (q, p) = Update_positions_momenta(q, p, t < c )

12:

Reflext_normal(p(t < c ), n c )

13:

Set dist = dist -t < c 14:

Is_strictly_in_convex[iRRAM :: Real](q)

15:

Return q

Algorithm S 3 Hamiltonian Monte Carlo step with mixed precision. The statistics updated are as follows: #R: multi-precision refinements triggered in iRRAM; #E: the number of exits of the polytope. See Section 4.3. In this pseudo-code, CGAL refers to the numbers types described in section 3.4.4, which are used in a lazy manner using a try-catch.

1: HMC_step(q) 2: Choose a travel time L ∼ U nif (0, L max ).

3: choose p ∼ N (0, I n ) double precision 4: q proposed = HM C_step < double > (q, p, L) 5: if Is_strictly_in_convex[CGAL](q proposed ) then 6:

Return q proposed 7: else q iRRAM = HM C_step < iRRAM > (q, p, L) 10: q proposed = to_double(q iRRAM ) 11:

if Is_strictly_in_convex[CGAL](q proposed ) then increment #E 15:

Return q

Algorithm S 4 Find trajectory parameters for a given direction.

1: Compute_traj_params(q dir , p dir ) 2: Set ω = √ 2a 3: Compute C = q 2 dir + p 2 dir /w 2 4: Compute φ = arctan(-p dir q dir ω ) 5: if p dir < 0 and φ < 0 then Compute the outward pointing normal n H to the hyper-plane 5:

Compute the dot products q n H =< q, n H > and p n H =< p, n H > 6:

[ω, C, Φ] = Compute_traj_params(q n H , p n H )

7:
if C > b i then if intersection = False then 16:

Set t c = t 17:

Set t c = n H

18:

Set intersection = True [ω, C, Φ] = Compute_traj_params (q i ,p i )

4:

Set q i = C cos(ωt + φ)

5:

Set p i = -ωC sin(ωt + φ) 6: Return (q, p) median(err r ) × stdev(err r ) median(time) × stdev(time)

σ = 1 σ

 1 = 0.01

Figure 2 :Figure 3 :

 23 Figure 2: Using HMC to sample a Gaussian of standard deviation σ restricted to the cube [-1, 1] 3 and the standard simplex. In all cases, N = 10, 000 samples. (Color gradient defined w.r.t. a coordinate value.)

Figure 4 :

 4 Figure 4: Estimating the volume of three polytopes with known volumes (cube, isotropic simplex, standard simplex), in dimension n = 10, 30, 50, 70, 100, with three algorithms. Statistics reported over 50 runs. The error of an estimate Ṽ is defined as err r = | Ṽ -V |/V . Symbol shape vs polytope: square: cube; triangle up: isotropic simplex; triangle down: standard simplex Symbol size vs dimension : smallest for n = 10, largest for n = 100 Color vs algorithm: purple/light blue/light green HAR/HMC0/HMC1.5 (Top row) Comparing the three algorithms on a given polytope (type and dimension fixed). (Bottom row) Assessing the incidence of dimension by running a given algorithm on instances of a fixed polytope type whose dimension decreases.

Lemma 3 .

 3 1. and 2. are simple consequences of the reversibility of the Hamiltonian flow with reflections. 3. Let A and B be measurable sets of R n . Ψ q (A B ) ⊂ B and Ψ(Ψ(A B )) = A B thus Ψ(A B ) ⊂ B A . By symmetry, Ψ(B A ) ⊂ A B . Hence by composing with Ψ: Ψ(Ψ(B A )) ⊂ Ψ(A B ). Using 2., we deduce B A ⊂ Ψ(A B ). 4. is clear.

Lemma 5 .Lemma 5 .

 55 Let B be an open ball in R n and let ϕ : B → R n be a differentiable map. If for each x in B, dϕ(x) -Id < 1, then ϕ is a diffeomorphism. The only thing to prove is that ϕ is injective. Let x = y be two points in B. Consider the map f : t ∈ [0, y -x ] → ϕ(x + tu) • u where u = y-x y-x . Using the Cauchy-Schwarz inequality and the assumption we obtain

Lemma 6 .Lemma 6 .

 66 Let B = B(0, r) be a closed ball in R n of center 0 and radius r > 0 and let ϕ : B → R n be an injective continuous map. If for all x in B, d(x, ϕ(x)) ≤ r/4 then ϕ(B) contains the ball B(0, r/2). By Jordan-Brouwer Theorem, the complement in R n of the image Σ of the sphere S = ∂B has exactly two connected components C 1 and C 2 , one which is bounded, say C 1 , and one which is not. By assumption the open ball• B(O, r/2) doesn't intersect Σ, hence is included in C 1 or C 2 .The image E = ϕ( • B) of the interior of the ball B is included in C i , one of the two connected components of R n \ Σ. On the one hand, by Jordan-Brouwer invariance of the domain theorem, E is open. On the other hand, E = ϕ(B) ∩ C i and since ϕ(B) is compact,

3 : 4 :

 34 choose p ∼ N (0, I n ) Set dist = L 5: while dist > 0 do 6:

3 : 4 :

 34 choose p ∼ N (0, I n ) iRRAM REAL Set dist = L 5: while dist > 0 do 6:

6 :φ = φ + π 7 : 2 :

 672 if p dir > 0 and φ > 0 then 8: φ = φ -π 9: Return [ω, C, φ] Algorithm S 5 Intersecting the trajectory with hyperplanes bounding the polytope 1: Intersect_hyper_planes(q, p) Set intersection = False 3: for each hyperplane H of equation (Ax) i = b i do 4:

8 :

 8 t1 = (arccos(b i /C) -φ) /ω

Figure S 1 : 5 Table S 2 : 1 1 . 3 4 2 2 . 2 2 . 3 8 3 1 31 7 33 8

 15211342223813133 Figure S 1: Estimating the volume of three polytopes with known volumes (cube, isotropic simplex, standard simplex), in dimension n = 10, 30, 50, 70, 100, with three algorithms. Statistics reported over 50 runs. The error of an estimate Ṽ is defined as err r = | Ṽ -V |/V . Symbol shape vs polytope: square: cube; triangle up: isotropic simplex; triangle down: standard simplex Symbol size vs dimension : smallest for n = 10, largest for n = 100 Color vs algorithm: purple/light blue/light green HAR/HMC0/HMC1.5

Table 1 :

 1 Summary of parameters used for all contenders. N cs : number of samples used for the incremental construction of the cooling schedule; N m : number of samples used in the cooling schedule termination condition; W : window size to estimate the R i s; r: parameter used to adjust the travel time. Nb: ε is the target quality.

		2 150/ε 500 + 4 * n 2 NA
	HMC0 HMC1.5	1000 + n 150/ε 1000 + n 150/ε 250 + n 250 0.1 √ n 0.1

  Set n c = n H 23: Return (intersection, t c )

	19:	else
	20:	if t < t c then
	21:	Set t c = t
	22:
	Algorithm S 6 Reflecting the normal
	1: Reflext_normal(p, n)
	2: n = n/ n // unit normal

3: Return p -2 < p, n > n Algorithm S 7 Update trajectory with distance t 1: Update_positions_momenta(q, p, t) 2: for i from 1 to n do 3:

Acknowledgments. We are grateful to B. Cousins and S. Vempala [34] for providing a high quality MATLAB code.

This work has been partially supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.

Let

A = {(α, t, q, p) : |α| ≤ 1/2, t ∈ [1/2, 3/2], q ∈ Q, p ∈ B(-q, ρ)}.

A is compact, hence the first and second order derivatives of φ are bounded on A. Thus there exists Z > 0 such that for every (α, t, q, p) ∈ A, φ q (α, t, q, p) -φ q (0, 1, q, p) ≤ (α + |t -1|)Z [START_REF] Masur | Chapter 13) Rational billiards and flat structures[END_REF] and d p φ q (α, t, q, p) -d v φ q (0, 1, q, p) ≤ (α + |t -1|)Z.

Using equation 24, the two above inequalities are equivalent to

and

The determinant is continuous, so there exists

Finally, using Lemma 4 together with the fact that ∇U is bounded, we deduce that there exists α 0 > 0 such that for every 0 < α ≤ α 0 and every q ∈ Q and p ∈ B(-q, ρ), the trajectory φ q (α, t, q, p) for t ≤ 1 stays in Q. It follows that φ and Φ coincide, for there is no reflection.

and q be in Q. By lemma 5 below, the map f : p → φ q (α, t, q, p) is a C 1 diffeomorphism from B(-q, ρ) to V = φ q (α, t, q, B(-q, ρ)) and by Lemma 6 (applied to f translated by -q), B(0, ρ/2) ⊂ V . Furthermore, equation 30 implies that for every

Hence, the push forward measure ξ of ν by f is non zero on B(0, ρ/2), and has a density

on B(0, ρ/2). Hence, under the condition that the travel time t is in [1 -β 2Z , 1], we get the following transition probability:

For the random walk, t is sampled uniformly in [0, 1], hence

which concludes the proof.

Remark 8. The previous proof uses n 0 = 1. We believe it should be possible to extend the proof to non convex Q by taking n 0 > 1, and some extra regularity assumptions on Q.

8.2 Dimensions n = 250, 500

stdev( Ṽ ol) med(err r ) stdev(err r ) 5 HMC-Wa250-Wb0.0-Wc0.0-r0.1 cube 250 0.05 1.809E+75 1.463E+75 2.285E+75 1.817E+75 3.463E+74 1.886E-01 1.037E-01 5 HMC-Wa250-Wb0.0-Wc0.0-r0.1 isotropic-∆ 250 0.05 NAN 1.306E+108 2.849E+108 2.445E+108 6.030E+107 NAN NAN 5 HMC-Wa250-Wb0.0-Wc0.0-r0.1 standard-∆ 250 0.05 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.000E+00 -1.000E+00 0.000E+00 Table S 3


Andrew Taylor

Scott Mcgrane

Martina Heer

Jonathan Beauchamp

Loic Briand

Do space conditions change flavour perception and decrease food intake by astronauts?

Keywords: space food, odour interaction, flavour perception

Several hundred astronauts have experienced spaceflight for periods of up to 15 months on visits to the International Space Station (ISS). Each astronaut's health is monitored and there is clear evidence that their calorie intake only meets around 80% of the recommended daily allowance. There is concern that this underconsumption may have negative consequences on astronauts' mental and physical health during the planned interplanetary missions which will last about three years. There are many hypotheses for the calorie deficit and we have focused on just one; the concept that chemical contaminants in recycled air and water could change the flavour quality of meals during spaceflight and therefore affect food intake. Using detailed NASA data on the composition of ISS recycled air and water, we applied the Dose-over-Threshold (DoT) concept to identify contaminants that might be odour active and could change the perceived overall food flavour. However, because of the large variation in published odour threshold values, the DoT concept delivered inconclusive results. The limitations of the DoT concept and ways to improve it are discussed.

Introduction

The welfare of astronauts is a key factor in space exploration. Since the first human spaceflight in 1961, and the launch of the International Space Station (ISS) over two decades ago, much has been learnt about the changes that occur to the human body during spaceflight over periods of three months to one year. One factor that is still not explained is why astronauts only consume around 80% of their daily calorie needs [START_REF] Heer | Nutrient supply during recent European missions[END_REF][START_REF] Smith | Evidence report: Risk factor of inadequate nutrition[END_REF]. While this seems to be tolerable for (short) visits to the ISS, the proposed crewed mission to Mars will take around three years and the weight loss could be detrimental to astronauts' physical health and mental wellbeing. Although there could be many reasons for underconsumption, one hypothesis is that the conditions onboard the spacecraft alter the perceived flavour of food [START_REF] Taylor | Factors affecting flavor perception in space: Does the spacecraft environment influence food intake by astronauts? Critical Reviews in Food Science and Food Safety[END_REF] and Earth-based studies have established that flavour quality is linked to palatability of food, as well as food intake [START_REF] Yeomans | Taste, palatability and the control of appetite[END_REF][START_REF] Yin | Effects of aroma and taste, independently or in combination, on appetite sensation and subsequent food intake[END_REF].

Within the ISS, air and water are extensively recycled to make the ISS more self-sufficient. Close monitoring of recycled air and water quality is carried out to ensure these supplies are safe for the astronauts [START_REF]Refinements to the methods for developing spacecraft exposure guidelines[END_REF][START_REF]Spacecraft maximum allowable concentrations for selected airborne contaminants[END_REF]. However, the recycled air and water contain contaminants, some of which could affect the flavour perception of food in space and decrease food intake. For example, NASA has identified certain contaminants that reduce water consumption (RWC) due to poor taste at very high levels, but that are acceptable at the levels reported by the regular analyses [START_REF] Ryder | Spacecraft water exposure guidelines (SWEGs). Lyndon B. Johnson Space Center[END_REF]. Further, ISS air contains a range of volatile compounds and a ten-fold higher carbon dioxide concentration compared to the levels found on Earth. Detailed analyses of ISS air and water samples from different missions are available from the NASA archives [START_REF] Straub | ISS expeditions 16 thru 20: Chemical analysis results for potable water[END_REF][START_REF] Straub | Chemical analysis results for potable water from ISS expeditions 21 through 25[END_REF][START_REF] Straub | International space station potable water quality for expeditions 26 through 29[END_REF][START_REF] Straub | Potable water quality for International Space Station expeditions 30-33[END_REF], where the concentrations of a wide range of contaminants (several hundred) are reported. The data published by NASA has focused primarily on safety but it is welldocumented that some contaminants in potable water are detectable by smell and taste, well below the safety levels set for health reasons (see Table 9 in [START_REF] Young | Taste and odour threshold concentrations of potential potable water contaminants[END_REF]). Despite the detailed chemical analyses of ISS air and water, our search of the space literature found no systematic studies on the potential taste activity, or odour activity, of these contaminants.

Current methods to relate the perceived intensity of taste and aroma compounds with their concentrations are limited to Stevens Law and the Dose-over-Threshold (DoT) concept (concentration of a compound divided by its odour threshold). Stevens law states that perceived intensity is related to concentration through a power law function, with different flavour compounds having different exponents; unfortunately, there are limited data available on the values of these exponents. In contrast, threshold data on a range of taste and odour compounds are available, through experimental [START_REF] Young | Taste and odour threshold concentrations of potential potable water contaminants[END_REF][START_REF] Van Gemert | Compilation of odor threshold values in air and water[END_REF] or calculated values [START_REF] Abraham | An algorithm for 353 odor detection thresholds in humans[END_REF]. Therefore, the DoT concept is really the only way to assess the potential flavour intensity of the compounds listed in the NASA analyses. Experimental values of odour threshold (OT) range over several orders of magnitude, apparently due to variation in odour presentation to panellists, as well as the difference between human assessors [START_REF] Schmidt | Making scents: Dynamic olfactometry for threshold measurement[END_REF], so a mean or median value is usually taken when calculating the DoT value, where values >1 indicate detectable odour or taste intensity while values <1 indicate a lack of detectable odour/taste.

A limitation of the DoT concept is that it considers each flavour compound individually, whereas there is a considerable amount of data that inter-and intra-actions between taste and odour compounds can deliver detectable odour when one or more of the compounds is below its individual taste/odour threshold. Studies in the 1960s measured the odour intensity of sub-threshold mixtures using sensory analysis and showed that there were some "additive" effects. More recently, the published data has been reviewed [START_REF] Ferreira | Revisiting psychophysical work on the quantitative and qualitative odour properties of simple odour mixtures: A flavour chemistry view. Part 1: Intensity and detectability[END_REF][START_REF] Ferreira | Revisiting psychophysical work on the quantitative and qualitative odour properties of simple odour mixtures: A flavour chemistry view. Part 2: Qualitative aspects[END_REF] and it is clear that sub-threshold interactions do occur, not just in fundamental studies using model mixtures, but also in real food and beverage products.

Therefore, assessing the potential odour/taste activity within the ISS air and water supplies is severely limited due to the following factors:

1. The DoT concept relies on the availability and accuracy of odour and taste threshold data.

2. The threshold data in compilations and individual publications cover a few hundred compounds only and time-intensive sensory analysis is the primary method to obtain threshold values. 3. Data are recorded as the detection and/or recognition threshold concentrations. 4. Published odour threshold values vary significantly due to differences in methodology as well as human variation and it is difficult to decide which value to use. 5. The DoT concept treats each flavour compound as a discrete entity and ignores the well-documented interactions between flavour compounds. 6. There is strong evidence of interactions between odour compounds that change perception.

Experimental

Data sources

Literature searches covered the mainstream food and flavour journals as well as journals associated with space travel (for example, Aviation Space and Environmental Medicine, New Space, Advances in Space Research) and publications from the American Institute of Aeronautics and Astronautics. In the absence of science-based information on astronauts' perception of space food, anecdotal information was obtained by a web search. The NASA archives contain reports, presentations and policy documents, which can be found at https://ntrs.nasa.gov/. There are also datasets that can be downloaded from https://nasa.github.io/data-nasa-gov-frontpage/. Data from the European Space Agency can be found at https://www.esa.int/About_Us/ESA_Publications/(archive).

Results and discussion

The RWC compounds identified by NASA as potentially flavour-active in the ISS water supply are shown in Table 1, along with their concentrations in samples of potable ISS water. Values vary across the different samples taken, so high and low levels are quoted in Table 1. [START_REF] Amoore | Odor as an aid to chemical safety: Odor thresholds compared with threshold limit values and volatilities for 214 industrial chemicals in air and water dilution[END_REF] were determined in drinking water as the concentration in water that would elicit a detectable odour in the headspace. c) The odour threshold (water) of 2-mercaptobenzothiazole is taken from Nielsen [START_REF] Nielsen | 2-mercapto-benzothiazole (MBT)[END_REF].

Table 1 uses OT data from a single source [START_REF] Amoore | Odor as an aid to chemical safety: Odor thresholds compared with threshold limit values and volatilities for 214 industrial chemicals in air and water dilution[END_REF] (with the exception of 2-mercaptobenzothiazole, taken from [START_REF] Nielsen | 2-mercapto-benzothiazole (MBT)[END_REF]) to calculate the DoT values, in an attempt to produce DoT values that could be compared with each other. On this basis, only trimethylamine shows a DoT >1, however, the OT value for trimethylamine in Table 1 (0.2 µg/L; [START_REF] Amoore | Odor as an aid to chemical safety: Odor thresholds compared with threshold limit values and volatilities for 214 industrial chemicals in air and water dilution[END_REF]) is much lower than the values quoted elsewhere (500 and 1700 µg/L; [START_REF] Van Gemert | Compilation of odor threshold values in air and water[END_REF]), so it is debatable whether trimethylamine is perceivable at the concentrations present in ISS water. Given the use of chlorine to disinfect drinking water, there are limited data on the flavour thresholds values of some organochlorine compounds [START_REF] Young | Taste and odour threshold concentrations of potential potable water contaminants[END_REF]. Using the Amoore & Hautala OT data, neither of the organochlorine compounds in Table 1 will be perceivable at their concentrations present in ISS water, and even using the lowest reported OTs, only dichloromethane has a DoT >1 if the OT is taken as 4.5 µg/L. A study of other compounds in ISS water depends on the availability of odour threshold data, requiring searches on individual compounds to establish the OT ranges. Our experience is that a web search can identify OT values in various scientific and commercial web sources although the units used vary. It is not always apparent whether the values are quoted as the OT in air or the OT above aqueous solutions and further searching and calculation is required to validate the values.

The examples given above illustrate the difficulty when applying DoT values to evaluate whether a contaminant could change the flavour of ISS water. A similar situation is found when examining the contaminants in air. Air circulation within the ISS lacks the convection element found on Earth (cold air cannot sink to the ground in microgravity!), so fans are used to circulate air and it is acknowledged that this is partially successful, but does lead to higher localised concentrations of contaminants in certain areas of the ISS. The stand out difference is that CO2 levels in ISS air are around 0.4%, 10-fold greater than on Earth. In humans, gaseous CO2 is not detected at concentrations up to 30% [START_REF] Shusterman | Review of the upper airway, including olfaction, as mediator of symptoms[END_REF], but the presence of CO2 in carbonated beverages suggests its detection is through the conversion of CO2 into carbonic acid by the enzyme carbonic anhydrase (Figure 1).

Figure 1: Schematic representation of carbonation perception. (A) CO2 is odourless to human but activates the taste and trigeminal systems. (B) The reaction catalysed by carbonic anhydrase (CA). (C) CA operates in CO2

sensing and provokes a taste response to acidity. Reproduced from [START_REF] Taylor | Factors affecting flavor perception in space: Does the spacecraft environment influence food intake by astronauts? Critical Reviews in Food Science and Food Safety[END_REF].

Carbonic acid has been reported to activate not only the TRPA1 but also the sour taste receptors. Most published work on the effects of CO2 on flavour perception comes from studies on carbonated beverages where decreases in sweetness and increases in bitter after-taste have been noted [START_REF] Taylor | Factors affecting flavor perception in space: Does the spacecraft environment influence food intake by astronauts? Critical Reviews in Food Science and Food Safety[END_REF]. Recirculated ISS air also contains numerous other volatile organic compounds (Figure 2) that are derived from two main sources: the metabolic products from the astronauts and "off-gassing" of volatile compounds from the construction materials used within the ISS [START_REF] Macatangay | Impacts of microbial growth on the air quality of the International Space Station. 40th international conference on environmental systems[END_REF]. Depending on the level of physical activity by the astronauts, and the nature of the payloads that are uploaded to carry out different types of scientific experiments, the concentrations of airborne contaminants can vary from day to day, although CO2 levels are kept fairly constant by the recycling plant because elevated levels are known to adversely affect astronauts [START_REF] Law | Relationship between carbon dioxide levels and reported headaches on the International Space Station[END_REF]. [START_REF] Macatangay | Impacts of microbial growth on the air quality of the International Space Station. 40th international conference on environmental systems[END_REF].

Figure 2: Composition of volatile compounds in ISS air (volume basis). Left hand pie chart shows composition excluding nitrogen and oxygen, right hand chart shows breakdown of non-methane hydrocarbons

In Figure 2, siloxanes are a class of compounds released through off-gassing, while methane, ketones, alcohols and aldehydes are produced by human (or microbial) metabolism [START_REF] Macatangay | Impacts of microbial growth on the air quality of the International Space Station. 40th international conference on environmental systems[END_REF]. DoT analysis of the volatile compounds listed by NASA suggests that only the C2 to C7 aldehydes and n-butanol are potentially odour active [START_REF] Taylor | Factors affecting flavor perception in space: Does the spacecraft environment influence food intake by astronauts? Critical Reviews in Food Science and Food Safety[END_REF], while methane is described in the literature as odourless. Although the data on chemical analysis of ISS water and air is detailed and robust, systematic information on astronauts' perception of the odour inside the ISS or the taste of food in space was not found in our literature search. The anecdotal information from astronauts during debriefing or during public interviews is also scarce and variable (see [START_REF] Taylor | Factors affecting flavor perception in space: Does the spacecraft environment influence food intake by astronauts? Critical Reviews in Food Science and Food Safety[END_REF] for more information). Carrying out conventional sensory tests in orbit is currently restricted, due to the small numbers of participants available, as well as the difficulty in delivering odours or tastes under microgravity conditions.

Conclusion

Our attempts to use the DoT concept to identify compounds that might cause changes in flavour quality during spaceflight have led to inconclusive results and there are many reasons for this. Obviously, the factors described in the Introduction play a role (variation in odour threshold values, potential odour interactions) and this is not helped by the lack of sensory data from space-based experiments. To make further progress, it is possible to design Earth-based experiments to reproduce some of the conditions on the ISS and test whether they affect flavour perception of food samples. For example, the effect of higher CO2 levels on flavour perception of a test food could be assessed by a sensory panel, while exposed to atmospheres containing 0.04 and 0.4% CO2, providing a suitable ventilation system could be designed and the safety of the panellists could be ensured. Designing Earth-based experiments to determine whether the levels of odour-active contaminants found on the ISS could affect the perceived flavour of food are possible, but could not easily reproduce the effects of microgravity nor the composition of the ISS atmosphere. For contaminants found in the ISS potable water supply, it is technically feasible to add selected contaminants to food and drink samples and test whether their addition causes a significant sensory change. The potential sticking point is the sheer number of compounds found in the ISS air and water supplies, as well as the need to obtain permission from the safety and ethical perspectives. The situation may be made easier by the extensive testing that has been carried out to assess the safe levels of contaminants for astronauts over extended exposure periods. Space-based experiments are possible but there are restrictions on the weight, size and time needed to set up and execute experiments. Sensory experiments would have to be longitudinal-type studies since the number of participants during any one mission is limited to two to three astronauts. Studies would involve choosing odour or taste stimuli that could be measured in individual astronauts before, during and after the spaceflight, and repeated over several missions to get sufficient data that could be statistically analysed. The stimulus would have to be light, easy to use and reproducible, so the options are really limited to using something like odour pens as the stimuli. In turn, these experimental design choices restrict the results that could be obtained and only simple measures could be studied, such as determining odour intensity or odour threshold of single compounds. Both outputs are far away from the ultimate goal of determining whether a contaminant can alter the flavour of a space food to such an extent that it is under-consumed.

The other conclusion of this work is to question whether the DoT concept is a useful and valid tool to identify odour-active compounds. The wide variability of odour thresholds reported in the literature means that compounds may, or may not be odour active, depending on which OT value is to be believed. Genetic variability in the human taste and odour receptors is well known [START_REF] Trimmer | Genetic variation across the human olfactory receptor repertoire alters odor perception[END_REF], and accounts for some of the variation, while anatomical and physiological differences can affect the degree to which flavour stimuli are transferred to the receptors [START_REF] Salles | In-mouth mechanisms leading to flavor release and perception[END_REF].

Therefore, it could be argued that any sensory panel will produce a different set of threshold values due to genetic and physiological differences. It is difficult to see how an "average" panel could be convened to provide a set of benchmark threshold values. Could molecular biology help by assembling odour and taste receptors that are representative of the general population and determining threshold values using EC50 data? This is possible, but it still does not address the strong evidence that some flavour compounds can affect flavour perception below their odour thresholds. The evidence for odour-odour interactions is compelling [START_REF] Ferreira | Revisiting psychophysical work on the quantitative and qualitative odour properties of simple odour mixtures: A flavour chemistry view. Part 1: Intensity and detectability[END_REF]. Psychophysical data from the 1960s, showed that subthreshold levels of one odorant could change the perception of another odorant (see, for example, [START_REF] Guadagni | Additive effect of sub-threshold concentrations of some organic compounds associated with food aromas[END_REF]) and more recent work has found similar interactions in a vodka-based beverage [START_REF] Niu | Characterization of odor-active volatiles and odor contribution based on binary interaction effects in mango and vodka cocktail[END_REF]. Besides the psychophysical effects, interactions may occur at other stages in the flavour perception process. For example, competition for enzymes at the peri-receptor level can change the odour profile reaching the odour receptors and therefore change perception [START_REF] Hanser | Odorant-odorant metabolic interaction, a novel actor in olfactory perception and behavioral responsiveness[END_REF]; some compounds appear to modulate the signal from the receptors [START_REF] Xu | Widespread receptor-driven modulation in peripheral olfactory coding[END_REF], while others act as antagonists at the odour receptors [START_REF] Zak | Antagonistic odor interactions in olfactory sensory neurons are widespread in freely breathing mice[END_REF]. There are also taste-odour interactions at subthreshold levels that alter the sensory properties [START_REF] Dalton | The merging of the senses: Integration of subthreshold taste and smell[END_REF] or the neural coding patterns [START_REF] March | Modulation of the combinatorial code of odorant receptor response patterns in odorant mixtures[END_REF]. Given the complexity of the data, Artificial Intelligence (AI) is a potential tool to analyse the data available. AI studies have focused on the link between odour structure and odour quality through the DREAM project [START_REF] Keller | Predicting human olfactory perception from chemical features of odor molecules[END_REF] and similar work on bitterness has also been reported [START_REF] Dagan-Wiener | Bitter or not? Bitterpredict, a tool for predicting taste from chemical structure[END_REF][START_REF] Margulis | Intense bitterness of molecules: Machine learning for expediting drug discovery[END_REF]. These reports do not consider the intensity of the odour or bitter compound, nor do they calculate thresholds. A method to calculate odour thresholds based on physicochemical properties was proposed in 2002 [START_REF] Abraham | A model for odour thresholds[END_REF] and a recent publication described the use of machine learning to predict odour-odour interactions [START_REF] Yan | Visual analysis of odor interaction based on support vector regression method[END_REF]. AI is a powerful tool but, for a successful outcome, high-quality data is still needed and it is unclear how and where this data will come from. Identifying the problems in using DoT values is simple, but finding a solution to provide flavour researchers with a more sophisticated tool to understand how interactions affect flavour perception is a massive task and will require a truly multidisciplinary effort to obtain quality data from sensory, genetic and physicochemical sources, then combine the datasets to investigate the overall effect of odour-odour, taste-taste and taste-odour interactions. Until progress is made in these areas, it is difficult to test the hypothesis that contaminants in the recycled air and water of spacecraft can alter the flavour of food and decrease food intake during spaceflight. Spiking some food samples with known levels of contaminants is theoretically possible but raises ethical and safety issues, and deciding which of the several hundred contaminants to use makes experimentation more trial and error than sharply focused.

Table 1 . Estimation of the sensory effect of RWC compounds in ISS potable water supplies Compound Lowest level a

 1

			Highest	OT in	DoT	DoT
			level a	water b	level/OT	level/OT
		µg/L (ppb)	µg/L (ppb) µg/L (ppb)	Low	High
	Methylethylketone	5	20	8400	0.0006	0.024
	Dichloromethane	5	20	9100	0.0005	0.0022
	Chloroform	5	20	2400	0.0021	0.0083
	Ammonia	2	34	1500	0.0013	0.022
	Acetone	3	34	20000	0.0002	0.00017
	Methylamine	125	125	2400	0.052	0.052
	Trimethylamine	125	125	0.20	625	625
	2-Mercaptobenzothiazole	40	320	1760 c	0.0227	0.182
	Phenol	4	20	790	0.0005	0.0025

a) Data taken from Appendix 4 in Straub et al.

[START_REF] Straub | ISS expeditions 16 thru 20: Chemical analysis results for potable water[END_REF]

. The data in the reference are quoted in increments e.g., <5 or <20 µg/L but, to calculate DoT values they have been quoted in the Table as discrete values. b) Odour threshold values

Acknowledgements

The European Space Agency funded two meetings of the team through a Topical Team award (Contract # 4000125158). Victor Demaria Pesce (INSERM, ESA) and Ines Antunes (ESA) are thanked for their support of the project. We appreciate constructive discussions on this topic with further members of the Topical Team, namely Thomas Hummel, Christian Margot, Paola Pittia, Charles Spence and Serge Pieters. We are grateful to Neil DaCosta (IFF) for organising the 2019 ACS Orlando meeting on "Food for Space Travel & Extreme Environments" and connecting researchers from different backgrounds.