René Lozi

A Dieudonné

COMPLEXITY LEADS TO RANDOMNESS IN CHAOTIC SYSTEMS

Complexity of a particular coordinated system is the degree of difficulty in predicting the properties of the system if the properties of the system's correlated parts are given. The coordinated system manifests properties not carried by individual parts. The subject system can be said to emerge without any "guiding hand". In systems theory and science, emergence is the way complex systems and patterns arise out of a multiplicity of relatively simple interactions. Emergence is central to the theories of integrative levels and of complex systems. The emergent property of the ultra weak multidimensional coupling of p 1-dimensional dynamical chaotic systems for which complexity leads from chaos to randomness has been recently pointed out. Pseudorandom or chaotic numbers are nowadays used in many areas of contemporary technology such as modern communication systems and engineering applications. Efficient Chaotic Pseudo Random Number Generators (CPRNG) have been recently introduced. They use the ultra weak multidimensional coupling of p 1-dimensional dynamical systems which preserves the chaotic properties of the continuous models in numerical experiments. Together with chaotic sampling and mixing processes, the complexity of ultra weak coupling leads to families of CPRNG which are noteworthy. In this paper we improve again these families using a double threshold chaotic sampling instead of a single one. A window of emergence of randomness for some parameter value is numerically displayed. Moreover we emphasize that a determining property of such improved CPRNG is the high number of parameters used and the high sensitivity to the parameters value which allows choosing it as cipher-keys.

Introduction

Characterizing complexity is not easy and there are in science a number of approaches to do it. Many definitions tend to postulate or assume that complexity expresses a condition of numerous elements in a system and numerous forms of relationships among the elements. Some others definitions relate to the algorithmic basis for the expression of a complex phenomenon or model or mathematical expression. Warren Weaver [1] has posited that the (organized) complexity of a particular system is the degree of difficulty in predicting the properties of the system if the properties of the system's parts are given. In Weaver's view organized complexity, resides in nothing else than the non-random, or correlated, interaction between the parts. These correlated relationships create a differentiated structure which can, as a system, interact with other systems. The coordinated system manifests properties not carried by individual parts. The organized aspect of this form of complexity versus other systems than the subject system can be said to emerge without any "guiding hand". The number of parts does not have to be very large for a particular system to have emergent properties.

In systems theory and science, emergence is the way complex systems and patterns arise out of a multiplicity of relatively simple interactions. Emergence is central to the theories of integrative levels and of complex systems (M. A. Aziz-Alaoui et al. [START_REF] Aziz-Alaoui | From System Complexity to Emergent Properties (Understanding Complex Systems[END_REF]).

In this paper we use the emergent property of the ultra weak multidimensional coupling of p 1-dimensional dynamical chaotic systems for which complexity leads from chaos to randomness. Efficient Chaotic Pseudo Random Number Generators (CPRNG) have been recently introduced (Lozi [3,[START_REF] Lozi | Proceedings of 6th EUROMECH Non Linear Dynamics Conference[END_REF][START_REF] Lozi | Conference Proceedings of ICCSA 2009 20[END_REF][START_REF] Lozi | [END_REF]) and their properties analyzed (Hénaff et al. [START_REF] Hénaff | Conference Proceedings of ICCSA 2009[END_REF][START_REF] Hénaff | [END_REF][START_REF] Hénaff | Proceedings of the Physics and Control Conference[END_REF][START_REF] Hénaff | [END_REF]). The idea of applying discrete chaotic dynamical systems, intrinsically, exploits the property of extreme sensitivity of trajectories to small changes of initial conditions. The ultra weak multidimensional coupling of p 1-dimensional dynamical systems preserves the chaotic properties of the continuous models in numerical experiments. The process of chaotic sampling and mixing of chaotic sequences, which is pivotal for these families, works perfectly in numerical simulation when floating point (or any multi-precision) numbers are handled by a computer.

It is noteworthy that these families of ultra weakly coupled maps are more powerful than the usual formulas used to generate chaotic sequences mainly because only additions and multiplications are used in the computation process; no division being required. Moreover the computations are done using floating point or double precision numbers, allowing the use of the powerful Floating Point Unit (FPU) of the modern microprocessors (built by both Intel and Advanced Micro Devices (AMD)). In addition, a large part of the computations can be parallelized taking advantage of the multicore microprocessors which appear on the market of laptop computers.

In this paper we improve the properties of these families using a double threshold chaotic sampling instead of a single one. The genuine map f used as one-dimensional dynamical systems to generate them is henceforth perfectly hidden. A window of emergence of randomness for some parameter value is numerically displayed.

A determining property of such improved CPRNG is the high number of parameters used ( ( 1) p p ×for p coupled equations) which allows to choose it as cipher-keys due to the high sensitivity to the parameters values. We call these families multi-parameter chaotic pseudo-random number generators (M-p CPRNG). Several applications can be found for these families as for example producing Gaussian noise, computing hash function or in chaotic cryptography.

In Sec. 2 we review some of the most popular chaotic mappings in low dimension in the scope of their use in numerical algorithms and PRNG.

In Sec. 3 we improve the properties of ultra weak multidimensional coupled of p 1-dimensional dynamical chaotic systems using a double threshold chaotic sampling instead of a single one.

In Sec. [START_REF] Lozi | Proceedings of 6th EUROMECH Non Linear Dynamics Conference[END_REF] we describe the emergence of randomness from complexity in a particular window of parameter value. We point out the parameter sensitivity in Sec. 5, with some applications of the M-p CPRNG and we give a conclusion in Sec. 6.

Discrete Dynamical Systems in Low Dimension

Chaotic dynamical systems in low dimension are often used since their discovery in the 70' in order to generate chaotic numbers, because they are very easy to implement in numerical algorithms [START_REF] Sprott | Chaos and Time-Series Analysis[END_REF]. However, as we point out in this section the computation of numerical approximation of their periodic orbits leads to very different results from the theoretical ones. Then they are unable to generate Pseudo Random Numbers (PRN). We review some of the most used maps in dimension from 1 to 3 in this scope.

1-Dimensional Chaotic Dynamical Systems

Logistic map

The very well known logistic map

[ [ [ [ ] ] ] ] [ [ [ [ ] ] ] ] a g : 0,1 0,1 → → → → is simply defined as ( ( ( ( ) ) ) ) ( ( ( ( ) ) ) ) a g x ax 1 x = - = - = - = - (1)
and generally considered for

[ [ [ [ ] ] ] ] a 0,4 ∈ ∈ ∈ ∈
(see Fig. 1). It is associated to the discrete dynamical system [START_REF] Alligood | Chaos. An introduction to dynamical systems Springer[END_REF] ( ( ( ( ) ) )

) This dynamical system which has excellent ergodic properties on the real interval

n 1 a n x g x + + + + = = = = (2)

[ [ [ [ ] ] ] ]

0,1 has been extensively studied especially by R. M. May [13], and J.

Feigenbaum [14] who introduced what is now called the Feigenbaum's constant 4.66920160910299067185320382...

δ = = = =
explaining by a new theory (period doubling bifurcation) the onset of chaos.

For every value of a there exist two fixed points: x 0 = = = = which is always unstable and

a 1 x a - -- - = = = = which is stable for ] ] ] ] [ [ [ [ a 1,3 ∈ ∈ ∈ ∈ and unstable for ] ] ] ] [ [ [ [ ] ] ] ] [ [ [ [ a 0,1 1,4 ∈ ∪ ∈ ∪ ∈ ∪ ∈ ∪ . When a 4 = = = = , the system is chaotic. The set { { { { } } } } 5 5 5 5 , 0.3454915028,0.9045084972 8 8         - + - + - + - +         = = = =                        
is the period-2 orbit. In fact there exist infinity of periodic orbits and infinity of periods. This dynamical system possesses an invariant measure

P( x )

x( 1 x )

π = = = = - -- - (see

Numerical approximation of the logistic map

In order to compute longer periodic orbits the use of computer is required, as it is equivalent to find roots of polynomial equation of degree greater than 4 for which Galois theory teaches that no closed formula is available. However, numerical computation uses ordinarily double precision numbers (IEEE-754) so that the working interval contains roughly 10 16 representable points. Doing such a computation in Eq. ( 2) with 1,000 randomly chosen initial guesses, 596, i.e., the majority, converge to the unstable fixed point x 0 = = = = , and 404 converge to a cycle of period 15,784,521. (see Table 1) [15]. Thus, in this case at least, the very long-term behaviour of numerical orbits is, for a substantial fraction of initial points, in flagrant disagreement with the true behaviour of typical orbits of the original smooth logistic map.

In others numerical experiments we have performed, the computer working with fixed finite precision is able to represent finitely many points in the interval in question. It is probably good, for purposes of orientation, to think of the case where the representable points are uniformly spaced in the interval. The true logistic map is then approximated by a discretized map, sending the finite set of representable points in the interval to itself.

Describing the discretized mapping exactly is usually complicated, but it is roughly the mapping obtained by applying the exact smooth mapping to each of the discrete representable points and "rounding" the result to the nearest representable point. In our experiments [16,17], uniformly spaced points in the interval with several order of discretization (ranging from 9 to 2,001 points) are involved, the results for 2,000 and 2,001 points are displayed in Table 2. In each experiment the questions addressed are: • how many periodic cycles are there and what are their periods ? • how large are their respective basins of attraction, i.e. , for each periodic cycle, how many initial points give orbits with eventually land on the cycle in question ? The existence of very short periodic orbit (see Table 1), the existence of a non constant invariant measure (see Fig. 2) and the easily recognized shape of the function in the phase space ( ( ( ( ) ) ) )

n n 1

x ,x + + + + avoid the use of the logistic map as a PRN generator. However, its very simple implementation in computer program led some authors to use it as a base of cryptosytem [18,19].

[ [ [ [ ] ] ] ] a J 1,1 , f : J J = - → = - → = - → = - → a f ( x ) 1 a x = - = - = - = - (3) ( ( ( ( ) ) )
)

n 1 a n x f x + + + + = = = = (4)
Despite its simple shape (see Fig. 3), it has several interesting properties. First, when the parameter value a 2 = = = = , the system possesses chaotic orbits. Because of its piecewise-linear structure, it is easy to find those orbits explicitly. More, owing to its simple definition, the symmetric tent map's shape under iteration is very well understood. The invariant measure is the Lebesgue measure. Finally, and perhaps the most important, the tent map is conjugate to the logistic map, which in turn is conjugate to the Hénon map for small values of b [START_REF] Alligood | Chaos. An introduction to dynamical systems Springer[END_REF]. However the symmetric tent map is dramatically numerically instable:

Sharkovskiǐ's theorem applies for it [20]. When a 2 = = = = there exists a period three orbit, which implies that there is infinity of periodic orbits. Nevertheless the orbit of almost every point of the interval J of the discretized tent map converges to the (unstable) fixed point x 1 = -= -= -= -(this is due to the binary structure of floating points) and there is no numerical attracting periodic orbit [START_REF] Sprott | Chaos and Time-Series Analysis[END_REF]. The numerical behaviour of iterates with respect to chaos is worse than the numerical behaviour of iterates of the approximated logistic map. This is why the tent map is never used to generate numerically chaotic numbers. However in Sec. 3 we will show that it is possible to preserve its chaotic properties when several logistic maps are ultra weakly coupled.

2-Dimensional System Chaotic Dynamical Systems

Hénon map

In order to study numerically the properties of the Lorenz attractor [START_REF] Sprott | Chaos and Time-Series Analysis[END_REF], M. Hénon in 1976 [21] introduced a simplified model of the Poincaré map [START_REF] Alligood | Chaos. An introduction to dynamical systems Springer[END_REF] of this attractor. The Lorenz attractor being in dimension 3, the corresponding Poincaré map is a map from 2 ℝ to 2 ℝ . The Hénon map is then also defined in dimension 2 as

2 x y 1 ax F : y bx         + - + - + - + -         → → → →                                 (5)
It is associated to the dynamical system

2 n 1 n n n 1 n
x y 1 ax y bx

+ + + + + + + +     = + - = + - = + - = + -     = = = =     (6)
For the parameter value a 1.4

= = = = , b 0.3 = = = =
Hénon pointed out numerically that there exists an attractor with fractal structure (see Fig. 4). This was the first example of strange attractor (previously introduced by D. Ruelle and F. Takens [START_REF] Ruelle | and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM[END_REF]) for a mapping defined by an analytic formula.

Nowadays hundreds of research papers have been published on this prototypical map in order to fully understand its innermost structure. However as in dimension 1, there is a discrepancy between the mathematical properties of this map in the plane 2 ℝ and the numerical computations done using (IEEE-754) double precision numbers.

If we call Megaperiodic orbits [23], those whose length of the period belongs to the interval of natural numbers

6 9 10 ,10               
 and Gigaperiodic orbits, those whose length of the period belongs to the interval i.e. two hundred forty times longer than the longest period of the one-dimensional logistic map (see Table 1). The same digits of these points are bold printed, they are very close.

9 12 10 ,10                 ,
Nevertheless, as displayed in Fig. 4 the orbit are not uniformly distributed on the phase space, then it is not possible to use this map as a PRN generator.

Beside the problem of PRN generator, logistic and Hénon maps are recently used together with a secret key, by N. Pareek et al. [START_REF] Pareek | [END_REF], in order to build a chaotic block cipher which is extremely robust, due to the excellent confusion and diffusion properties of these maps. The results of the statistical analysis show that the chaotic cipher possesses all features needed for a secure system and useable for the security of communication system.

L. dos Santos Coelho et al. [25], introduced a chaotic particle swarm optimisation (PSO) which is a population-based swarm intelligence algorithm driven by the stimulation of a social psychological metaphor instead of the survival of the fittest individual. Based on the chaotic systems theory (using Hénon map sequences which increase its convergence rate and resulting precision) the novel PSO combined with an implicit filtering allows solving economic dispatch problems.

Lozi map

The Lozi map [26] is a linearized version of the Hénon map, built in order to simplify the computations, mainly because it is possible to compute explicitly any periodic orbits solving a linear system. It is defined as

n 1 n n n 1 n
x y 1 a x y bx

+ + + + + + + +     = + - = + - = + - = + -     = = = =     (7) or equivalently n 1 n n 1 x 1 a x bx + - + - + - + - = - + = - + = - + = - + (8)
For a 1.7

= = = = and b 0.5 = = = =
there exists a strange attractor. The particularity of this strange attractor is that it has been rigorously proved by Misiurewicz in 1980 [27].

In the same conditions of computation as for Hénon map, running the computation during nineteen hours, one can find a Gigaperiodic attracting orbit of period 436,170,188,959 more than one hundred times longer than the period of Orbit 1 found for the Hénon map. ). There is a transient regime before the orbit is reached. It seems that there is no periodic orbit with a smaller length. This could be due to the quasihyperbolic nature of the attractor. However, the orbit-shifted shadowing property of Lozi map (and generalized Lozi map), which is the property which ensures that pseudo-orbits of a homeomorphism f can be traceable by actual orbits even if rounding errors in computing are not inevitable has been recently proved [28].

Hence this attractor is very efficient, in order to generate chaotic numbers without repetition for standard simulation using either the first or the second component. However they are not equally distributed on the plane (see Fig. 5). The non constant density forbids its direct use as a PRN generator. Nevertheless there are some U.S. patents for "method of generating pseudo-random numbers in an electronic device, and a method of encrypting and decrypting electronic data" in which the Lozi map is involved [29,30]. 12 founding paper, R. Caponetto et al. [31] propose an experimental analysis on the convergence of EA. The effect of introducing chaotic sequences instead of random ones during all the phases of the evolution process is investigated. The approach is based on the substitution of the PRNG with chaotic sequences. Several numerical examples are reported in order to compare the performance of the EA using random and chaotic generators as regards to both the results and the convergence speed. The results obtained show that chaotic sequences obtained from Lozi map are always able to increase the value of some measured algorithm-performance indexes with respect to random sequences. Several authors following this idea use Lozi map in chaotic optimization in order to avoid local optima stagnation and embed a superior search strategy [32 -40].

3-Dimensional System Chaotic Dynamical Systems

In order to generalize in higher dimension the tent map, G. Manjunath et al. [41] introduce a three-dimensional map

3 3 F : I I → → → → where [ [ [ [ ] ] ] ] I 0,1 = = = = (see Fig. 6)
which is continuous in the Euclidian topology and prove its chaotic properties:

y z 1 2x 1 2 x z F( x, y,z ) 1 2 y 1 2 x y 1 2z 1 2         + + + + - + - - + - - + - - + -                         + + + +         = - + - = - + - = - + - = - + -                 + + + +         - + - - + - - + - - + -                 (9)
The related dynamical system is

( ( ( ( ) ) ) ) ( ( ( ( ) ) )
)

n 1 n 1 n 1 n n n
x , y ,z , F x , y ,z ,

+ + + + + + + + + + + + = = = = (10)
They emphasize that most of the well known examples of higher dimensional chaotic dynamical systems belong to the class of hyperbolic diffeomorphisms on a n-torus. These higher dimensional maps on the torus are not continuous on the standard topology of the Euclidean space since they exhibit jump discontinuities. The realization of such jump discontinuities in an electronic circuit implementation is not reliable. They prove the following theorem:

δ         ∈ ∈ ∈ ∈                        
as sensitivity constant.

Once again the sequence of iterated points ( ( ( ( ) ) )

)

n n n
x , y ,z obtained from the dynamical system ( 10) is not equally distributed on the volume To allow the generation of PRN using complexity and emergence theory we consider in the next section how to generate these numbers with uniform repartition on a given interval, or on a given square of the plane or more generally in a given hypercube of

Multi-parameter Chaotic Pseudo-Random Number Generator (M-p CPRNG)

As previously seen, when a dynamical system is realized on a computer using floating point or double precision numbers, the computation is of a discretization, where finite machine arithmetic replaces continuum state space. For chaotic dynamical systems, the discretization often has collapsing effects to a fixed point or to short cycles [15,42]. In order to preserve the chaotic properties of the continuous models in numerical experiments we consider an ultra weak multidimensional coupling of p one-dimensional dynamical systems.

System of p-Coupled Symmetric Tent Map

In order to simplify the presentation of the M-p CPRNG we introduce, we use as an example the symmetric tent map previously defined (3), even though others chaotic map of the interval (as the logistic map, the baker transform, …) can be used for the same purpose (as a matter of course, the invariant measure of the chaotic map chosen is preserved).

The considered system of the p-coupled dynamical systems is described by

)

n 1 n n X F X A f ( X ) + + + + = = ⋅ = = ⋅ = = ⋅ = = ⋅ (11) with 1 1 ( ) , ( ) ( )
n n n n p p n n x f x X f X x f x ⋮ ⋮         = =             (12) and 1,1 1, 2 2,2 2, 1, 2 1 ,
, 1

j p j 1,2 1,p-1 1,p j j p 2,1 j 2,p-1 2,p j j j p p,1 p,p-1 p p p j j =1-ε ε ε ε ε =1- ε ε ε A ε ε =1- ε ε ε ε ⋯ ⋯ ⋮ ⋱ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋯ ⋯ = = = = ≠ = - =             =               ∑ ∑ ∑ ( 13
)
F is a map of [ ] 1,1 p p p J ℝ = - ⊂ into itself. Considering j p i , i i , j j 1, j i = 1- ε ε = = ≠
∑ , the matrix A is always a stochastic matrix iff the coupling constants verify 0 i , j ε > for every i and j.

If i, j i j i,j ε 0 ≠ ∀ = the maps are totally decoupled, instead they are fully crisscross coupled when for example i, j i j

1 ε p 1 ≠ = -
. Generally, researchers do not consider very small values of , i j ε because it seems that the maps are quasi-decoupled with those values and no special effect of the coupling is expected. In fact it is not the case and ultra small coupling constants (as small as 10 -7 for floating point numbers or 10 -16 for double precision numbers), allows the construction of very long periodic orbits, leading to sterling chaotic generators. This is the way in complexity leads to randomness from chaos.

Chaotic Sampling and Mixing

However chaotic numbers are not pseudo-random numbers because the plot of the couples of any component ( (

)

l l n n 1 x , x + + + + of iterated points ( ) 1 , n n X X + in
the corresponding phase plane reveals the map f used as one-dimensional dynamical systems to generate them via Eq. ( 11). Nevertheless we have recently introduced a family of enhanced Chaotic Pseudo Random Number Generators (CPRNG) in order to compute very fast long series of pseudorandom numbers with desktop computer [START_REF] Lozi | [END_REF][START_REF] Lozi | Proceedings of 6th EUROMECH Non Linear Dynamics Conference[END_REF][START_REF] Lozi | Conference Proceedings of ICCSA 2009 20[END_REF]. This family is based on the previous ultra weak coupling which is improved in order to conceal the chaotic genuine function.

In order to hide f in the phase space ( )

l l n n 1
x , x + two mechanisms are used.

The pivotal idea of the first one mechanism is to sample chaotically the sequence

( ( ( ( ) ) ) ) l l l l l 0 1 2 n n 1
x , x , x , , x , x ,

+ + + + … … generated
( 0) (1) ( 2) ( ) ( 1) , , , , , , q q
l l l l l n n n n n x x x x x … … + denoted here ( ( ( ( ) ) )
)

0 1 2 q q 1
x , x , x , ,x , x ,

+ + + + ⋯ ⋯ of the original one, in the following way Given 1 , , l m p l m ≤ ≤ ≠ { } ( ) ( 1)
( ) The sequence ( ( ( ( ) ) ) )

, q l m q n q q r r n x x with n Min r n x T ℕ - - ∈ = -    = = > >   ( 1) 1
0 1 2 q q 1
x , x , x , ,x , x ,

+ + + + ⋯ ⋯ is then the sequence of chaotic pseudo-random numbers.
The mathematical formula ( 14) can be best understood in algorithmic way. The pseudo-code, for computing iterates of ( 14) corresponding to N iterates of ( 11) is:

( )

1 2 p 1 p 0 0 0 0 0 X x , x , , x , x seed … - = = n 0; q 0 = = do { while n < N do { while ( ) m n x T ≤ compute ( ) 1 2 p 1 p n n n n
x , x , , x , x … -

; n + +} compute ( )

1 2 p 1 p n n n n x , x , , x , x … - ; then ( ) q l q n n( q ) n ; x x ; = = n + +; q + +}
This chaotic sampling is possible due to the independence of each component of the iterated points n X vs. the others [START_REF] Lozi | [END_REF].

Remark 1: Albeit the number iter NSampl of pseudo-random numbers q x corresponding to the computation of N iterates is not known a priori, considering that the selecting process is again linked to the uniform distribution of the iterates of the tent map on J , this number is equivalent to 2 1

N T -.

A second mechanism can improve the unpredictability of the pseudo-random sequence generated as above, using synergistically all the components of the vector n X , instead of two. Given p -1 thresholds

1 2 1 p T T T J ⋯ - < < < ∈ (15)
and the corresponding partition of J

1 0 p k k J J ∪ - = = (16)
with

[ ] 0 1 J 1,T = - , ] [ 1 1 2 J T ,T = , [ [ 1 1 1 k k k J T ,T for k p + = < < -and 1 1 1 p p J T , - - - - - - - -         = = = =         ,
⋯ ⋯ by { } { } ( q ) k k ( 1 ) k p q n ( q ) k k ( q 1 ) r k 1 k p 1 r n 1 x x , with n Min s ( q ) Min r n x J ℕ - - ≤ ≤ - ∈ = -    = = = > ∈   (17)
The pseudo-code, for computing the iterates of ( 17) corresponding to N iterates of ( 11) is:

( )

1 2 p 1 p 0 0 0 0 0 X x , x , , x , x seed … - = = n 0; q 0 = = do { while n < N do {while ( ) p n 0 x J ∈ compute ( ) 1 2 p 1 p n n n n
x , x , , x , x … -

; n + +} compute ( )

1 2 p 1 p n n n n x , x , , x , x … - let k be such that p n k x J ∈ then ( ) q k q n n( q ) n ; x x = =
; n + +; q + +} Remark 2: In this case also, iter NSampl is not known a priori, however, considering that the selecting process is linked to the uniform distribution of the iterates of the tent map on J , one has Remark 3: This second mechanism is more or less linked to the whitening process [START_REF] Viega | and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM[END_REF]44].

Remark 4: Actually, one can choose any of the components in order to sample and mix the sequence, not only the last one.

Double Threshold Chaotic Sampling

On can eventually improve the CPRNG previously introduced with respect to the infinity norm instead of the 1

L or 2

L norms because the L ∞ norm is more sensitive than the others ones to reveal the concealed function f [START_REF] Lozi | Conference Proceedings of ICCSA 2009 20[END_REF]. For this purpose we introduce a second kind of threshold

' T ℕ ∈ , together with 1 1
, , p T T J

⋯ -∈ such that the subsequence ( ( ( ( ) ) ) )

0 1 2 q q 1
x , x , x , ,x , x ,

+ + + + ⋯ ⋯ is defined by { } { } ( q ) k k ( 1 ) k p q n ( q ) k k ( q 1 ) r k 1 k p 1 r n 1 x x , with n Min s ( q ) Min r n T ' x J ℕ - - ≤ ≤ - ∈ = -    = = = > + ∈   (18)
In pseudo-code Eq. ( 18) is then:

( )

1 2 p 1 p 0 0 0 0 0 X x , x , , x , x seed … - = = n 0, q 0 = = do { while n < N do {while ( ) p ( q 1 ) n 0 n n T ' and x J - ≤ + ∈ compute ( ) 1 2 p 1 p n n n n
x , x , , x , x … -

; n + +} compute ( ) process is linked to the uniform distribution of the iterates of the tent map on J , and to the second threshold T', it comes that

1 2 p 1 p n n n n x , x , , x , x … - let k be such that p n k x J ∈ then ( ) q k q n n( q ) n ; x x = = ; n + +; q + +} Remark
1 2 , 1 iter N N NSampl Min T T   ≤   ′ -  
.

Remark 6: the second kind of threshold ' T can also be used with only the chaotic sampling, without the chaotic mixing.

Emergence of Randomness

Numerical results about chaotic numbers produced by ( 11) - (17) show that they are equally distributed over the interval J with a very good precision [START_REF] Lozi | [END_REF][START_REF] Lozi | Proceedings of 6th EUROMECH Non Linear Dynamics Conference[END_REF].

In this section we emphasize that when the parameters , i j ε belong to a special window (called the window of emergence) the M-p CPRNG defined above behaves well.

Approximated Invariant Measures

In order to perform numerical computation, we have to define some numerical tools: the approximated invariant measures.

First we define an approximation , ( )

M N P
x of the invariant measure also called the probability distribution function linked to the 1-dimensional map f when computed with floating numbers (or numbers in double precision). In this scope we consider a regular partition of M small intervals (boxes) r i of J defined by:

i 2i s 1 , i 0, M M = -+ = i i i 1 r s , s , i 0, M -2 + = =     [ ] M 1 M 1 r s , 1 - - = 1 0 M i J r ∪ - = the length of each box is i 1 i 2 s s M + + + + -= -= -= -=
(note that this regular partition of J is different from the previous one linked to the threshold values i T (16)).

All iterates f (n)

M N i i 1 M P s r 2 N =
where #r i is the number of iterates belonging to the interval r i and the constant 1/2 allows the normalisation of M ,N P ( x ) on the interval J .

M ,N M ,N i i P ( x ) P ( s ) x r = ∀ ∈ = ∀ ∈ = ∀ ∈ = ∀ ∈
In the case of p-coupled maps, we are more interested by the distribution of each component ( ) , , , ,

1 2 1 p 2 x x x
x … of X rather than the distribution of the variable X itself in p J . We then consider the approximated probability distribution function , ( )

j M N P
x associated to one among several components of F(X) defined by [START_REF] Sprott | Chaos and Time-Series Analysis[END_REF] which are one-dimensional maps. In this paper we use equally disc N for M and iter N for N when they are more explicit.

The discrepancies 1

E (in norm 1 L ), 2 E (in norm 2 L ) and E ∞ (in norm L ∞ ) between disc iter
N , N P ( x ) and the Lebesgue measure which is the invariant measure associated to the symmetric tent map, are defined by

5

disc iter disc iter N N N N L E x P x = - 2 2, , , ( ) ( ) 0.
5

disc iter disc iter N N N N L E x P x = - , , , ( ) ( ) 0.
5

disc iter disc iter N N N N L E x P x ∞ ∞ = -
In the same way an approximation of the correlation distribution function , ( , )

M N C
x y is obtained numerically building a regular partition of 2 M small squares (boxes) of 2 J imbedded in the phase subspace (x l , x m ) , , , ,

i j 2i 2 j s 1 t 1 i j 0 M M M = -+ = -+ = i , j i i 1 j j 1 0, M-2 r
s , s t , t , i, j

+ + + + + + + +         = × = = × = = × = = × =                         [ [ [ ] i ,M 1 i i 1 M 1 r s , s t , 1 , j 0, M-2 - + - = × = Published In
[ ] [ ] M 1,M 1 M 1 M 1 r s , 1 t , 1 - - - - = ×
the measure of the area of each box is

( ( ( ( ) ) ) ) ( ( ( ( ) ) )
)

2 i 1 i i 1 i 2 s s t t M + + + + + + + +         -× -= -× -= -× -= -× -=                 Once N + Q iterated points ( ( ( ( ) ) )
)

l m n n
x , x belonging to these boxes are collected the relative number of iterates with respect to N/M 2 in each box r i,j represents the value C N (s i , t j ). The approximated probability distribution function C N (x, y) defined here is then a 2-dimensional step function, with M 2 steps. As M can take several values in the next sections, we define

( ( ( ( ) )
)

) 2 M ,N i j i , j 1 M C ( s , t ) # r 4 N = = = = (19)
where #r i,j is the number of iterates belonging to the square r i,j and the constant 1/4

allows the normalisation of M ,N C ( x, y ) on the square 2 J M ,N M ,N i j i , j C ( x, y ) C ( s ,t ) ( x, y ) r = ∀ ∈ = ∀ ∈ = ∀ ∈ = ∀ ∈ (20)
The discrepancies

1 C E (in norm 1 L ), 2 C E (in norm 2 L ) and C E ∞ (in norm L ∞ ) between disc iter
N , N C ( x, y ) and the uniform distribution on the square, are defined by

) 0.25

disc iter disc iter C N N N N L E x y C x y = - 2 2 , , , ( , ) ( , ) 0.25
disc iter disc iter C N N N N L E x y C x y = - , , , ( , ) ( , ) 0.25
disc iter disc iter C N N N N L E x y C x y ∞ ∞ = -
)

l l n n 1 x , x + + + +
instead of the phase space ( )

, l m
x x . In order to control that the enhanced chaotic numbers ( (

)

0 1 2 q q 1
x , x , x , ,x , x , ),

+ + + + ⋯ ⋯ are
disc iter disc iter AC N N q q AC N N q q E x x E x x + + , , 1
( , )

disc iter AC N N q q E x x ∞ + vanish).

A window of Emergence of Randomness

In order to point out the usefulness of the double threshold chaotic sampling with simply consider the case of only 4 coupled equations, and such that i , j i

disc iter disc iter disc iter N N N N N N E x E x E x ∞ and , , 1
( , ) ), . We choose, 0.9 T = and ' 20 T = . We display on Fig. 8 the values of the six computed error when

disc iter AC N N q q E x x ∞ + 1 2 , , 1 , , 1 ( , ), ( ,
disc iter disc iter AC N N q q AC N N q q E x x E x x + +
17 1 1 ε 10 ,10 - - - - - - - -         ∈ ∈ ∈ ∈       
 . The seed (initial values) being

- - - - - - - -         ∈ ∈ ∈ ∈       
 if one considers all together the six errors.

The errors , , ( ),

disc iter N N E x ∞ , , 1
( , )

disc iter AC N N q q E x x

The Underneath of Randomness

The double threshold chaotic sampling is very efficient because its aim is mainly to conceal f in the most drastic way. In order to understand the underneath mechanism consider first that in the phase space ( )

l l n n 1
x , x + the graph of the

25

chaotically sampled chaotic numbers is a mix of the graphs of the f (r) for all r ∈ ∈ ∈ ∈ ℕ (see Fig. 9).

Figure 9. Graphs of the symmetric tent map f, f (2) and f (3) on the interval [-1, 1].

It is obvious as showed on Fig. 10 The autocorrelation function is different from zero only if M > 2 (see Fig. 11).

In the same way as displayed on Fig. 12, 13 and 14, M N AC x y is constant for the symmetric tent map f (3) are shaded for M = 16. ( , )

Testing the Randomness

disc iter AC N N q q E x x ∞ + , 1 10

28

Computations are done using double precision numbers (~14-15 digits).

As shown previously [START_REF] Lozi | [END_REF] the errors in 1 L or 2 L norms decrease with the number of chaotic points (as in the law of large numbers) and conversely increase with the number M of boxes used to define , ( , )

M N AC x y . It is the same for the error in L ∞ norm. Fig. 15 shows that when M is greater than 2 5 , the sequence defined by (18) behaves better than the one defined by ( 14) or ( 17) when applied to Eq. ( 21). ( , ) ( , )

disc iter AC N N q q E x x ∞ + 1
disc iter AC N N q q E x x ∞ +
decreases drastically. If for example T ' 100 > > > >

, it is necessary to use a huge grid of 2 100 x2 100 boxes splitting the square 2 J in order to find a trace of the genuine function f. This is numerically impossible with double precision numbers. Then the chaotic numbers emerge as random numbers.

Applications

Generation of random or pseudorandom numbers, nowadays, is a key feature of industrial mathematics. Pseudorandom or chaotic numbers are used in many areas of contemporary technology such as modern communication systems and engineering applications.

More and more European or US patents using discrete mappings for this purpose are obtained by researchers of discrete dynamical systems [29,30].

When an efficient M-p CPRNG is defined, there exists a huge number of applications for the pseudo-random numbers it can generate, as for example chaotic masking, chaotic modulation or chaotic shift keying in the fields of secure communications [START_REF] Hénaff | Conference Proceedings of ICCSA 2009[END_REF][START_REF] Hénaff | [END_REF][START_REF] Hénaff | Proceedings of the Physics and Control Conference[END_REF][START_REF] Hénaff | [END_REF].

Parameter sensitivity

A determining property of the M-p CPRNG we have improved in this paper via Eq. ( 21) and double threshold chaotic sampling (18) is the high number of parameters used ( ( 1) p p ×for p coupled equations) which allows to choose it as cipher-keys however this achievement is possible only if there is a high sensitivity to the parameters values.

In order to point up this sensitivity, it is enough to consider the simplest case of 2-coupled equations with two sets of slightly different parameters ( )

, ε ε and ( ) ⋯ ⋯ differ completely as displayed in Table 3 (In fact all the components ( )

( q ) ( q ) 1 2 n n
x , x and ( )

( q ) ( q ) *1 *2 n n
x , x are different).

( q ) ( q ) 1 *1 n n
x , x and ( )

( q ) ( q ) 2 *2 n n
x , x of Eq. ( 22) and ( 23) with ( ) ( )

1 2 * *1 *2 0 0 0 0 0 0 X x , x X x , x = = =
Then rather than a unique CPRNG which is introduced here, there is a quasi-infinite family of CPRNG that the M-p CPRNG define allowing several possibilities of applications.

Gaussian Noise

As an example of such application, the generation of Gaussian noise from the sequences ( ( ( ( ) ) ) )

0 1 2 q q 1
x , x , x , ,x , x , + + + + ⋯ ⋯ is very easy when a Box-Muller transform is applied. A Box-Muller transform [START_REF] Box | [END_REF] is a method of generating pairs of independent standard normally distributed (zero expectation, unit variance) random numbers, given a source of uniformly distributed random numbers. The polar form [46] of such a transform takes two samples from a different interval, [ ] 1,1 and maps them to two normally distributed samples without the use of sine or cosine functions. This form of the polar transform is widely used, in part due to its inclusion in Numerical Recipes. As the sequences ( ( ( ( ) ) ) )

0 1 2 q q 1
x , x , x , ,x , x ,

+ + + + ⋯ ⋯ are uniformly distributed in [ ] 1,1 J ℝ = -
⊂ , the application is straightforward.

Hash Function

Another example of application could be the computation of hash function. A hash function is any well-defined procedure or mathematical function that converts a large, possibly variable-sized amount of data into a small one. The values returned by a hash function are called hash values, hash codes, hash sums, checksums or simply hashes.

Hash functions are mostly used to speed up table lookup or data comparison tasks -such as finding items in a database, detecting duplicated or similar records in a large file, finding similar stretches in DNA sequences, and so on.

A hash function may map two or more keys to the same hash value. In many applications, it is desirable to minimize the occurrence of such collisions, which means that the hash function must map the keys to the hash values as evenly as possible. Depending on the application, other properties may be required as well. Although the idea was conceived in the 1950s, the design of good hash functions is still a topic of active research.

Although hash function generally involve integers, on can consider that the application which maps the initial seed ( ) ) )

0 1 2 q q 1
x , x , x , ,x , x ,

+ + + +
⋯ ⋯ is a hash function working on floating point numbers. We will explore this application in a forthcoming paper. Others applications show the high-potency of such M-p CPRNG. Due to limitation of this article, they will be published elsewhere.

Conclusion

Using a double threshold in order to sample a chaotic sequence, we have improved with respect to the infinity norm the M-p CPRNG previously introduced. When the value of the second threshold T ' is greater than 100, it is impossible to find the genuine function used to generate the chaotic numbers. The new M-p CPRNG family is robust versus the choice of the weak parameter of the system for 10 -14 < ε < 10 -5 , allowing the use of this family in several applications as

Published

  In Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August 2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher, Singapore (2011) pp. 93-125.

2

 2

Figure 1 .

 1 Figure 1. Graph of the logistic map for a 4 = = = = .

Fig. 2 )Figure 2 .

 22 Fig. 2).

Figure 3 .

 3 Figure 3. Graph of the symmetric tent map on J for a 2 = = = = .

8

 8

  Hénon map possesses Gigaperiodic orbits. On a Dell computer with a Pentium IV microprocessor running at the frequency of 1.5 Gigahertz, using a Borland C compiler and computing with ordinary (IEEE-754) double precision numbers, one can find for a 1of length 3,800,716,788 Published In Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August 2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher, Singapore (2011) pp. 93-125.

9

 9

Figure 4 .

 4 Figure 4. The strange attractor of the Hénon map for a 1.4 = = = = and b 0.3 = = = = .

Figure 5 .

 5 Figure 5. Invariant density of the second component y of Eq. (7) computed using 10 10 iterations.

  Published In Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August 2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher, Singapore (2011) pp. 93-125. 13 Theorem (G. Manjunath et al.) The map defined in (9) is topologically transitive and exhibits sensitive dependence on initial conditions with any real number

3 I

 3 . The invariant density of the first component n x is displayed in Fig.7. The relative discrepancy of this invariant density versus the uniform one lies between 4% and 5%.

Figure 6 . Tessellation of 3 I

 63 Figure 6. Tessellation of 3 I into 3 3 regions by parallel set of critical planes

nℝFigure 7 .

 7 Figure 7. Invariant density of the first component x of Eq. (10) computed using 10 11 iterations.

( 14 )

 14 Published In Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August 2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher, Singapore (2011) pp. 93-125.

17

 17

  Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August 2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher, Singapore (2011) pp. 93-125.

20

 20

  (x) belonging to these boxes are collected (after a transient regime of Q iterations decided a priori, i.e. the first Q iterates are neglected). Once Published In Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August 2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher, Singapore (2011) pp. 93-125. 21 the computation of N + Q iterates is completed, the relative number of iterates with respect to N/M in each box r i represents the value ( ) in this article is then a step function, with M steps. As M may vary,

  y ) of (20) defined in the phase space ( (

Figure 10 .

 10 Figure 10. In shaded regions the autocorrelation distribution

26 Figure 11 .

 2611 Figure 11. Regions where the autocorrelation distribution

Figure 12 .

 12 Figure 12. In shaded regions the autocorrelation distribution

Figure 13 . 27 Figure 14 .

 132714 Figure 13. Regions where the autocorrelation distribution

Figure 15 .

 15 Figure 15. Error of

  0.9 and T' = 20, Niter = 10^9 Published In Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August 2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher, Singapore (2011) pp. 93-125.

Figure 16 .

 16 Figure 16. Error of

  .9 and T'= 20 , Niter = 10^9 Threshold T= 0.9 and T'= 20 , Niter = 10^10 Threshold T= 0.9 and T'= 20 , Niter = 10^11 Published In Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August 2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher, Singapore (2011) pp. 93-125. 29 digits).Fig.16shows that when the number of chaotic points increases the error , , 1

  Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August 2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher, Singapore (2011) pp. 93-125.

32

 32

Table 1

 1

	. Coexisting periodic orbits found using 1,000 random initial
		points for double precision numbers
	Period 1	Orbit = unstable fixed point x 0 = = =	Relative Basin size 596 over 1,000
	15,784,521	Scattered over the interval	404 over 1,000

Table 2 .

 2 Coexisting periodic orbits for the discretization with regular

		meshes of N 2,000 and 2,001 = = = =	points.
	N	Period	Orbit	Relative Basin Size
	2,000	1	{0}		2 over 2,000
	2,000	2	{1,499}		14 over 2,000
	2,000	2	{691;1,808}		138 over 2,000
	2,000	3	{276;1,221;1,900}		6 over 2,000
	2,000	8	{3;11;43;168;615;1,703.1,008.1,998} 1,840 over 2,000
	2,001	1	{0}		5 over 2,001
	2,001	1	{1500}		34 over 2,001
	2,001	2	{91; 1809}		92 over 2,001
	2,001	8	{3;11;43;168;615;1,703;1,011;1,999} 608 over 2,001
	2,001	18	{35;137;510;1,519;1,461;1,574;…}	263 over 2,001
	2,001	25	{27;106;401;1,282;1,840;588;…}	1,262 over 2,001

  Published In Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August 2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher,

	Singapore (2011) pp. 93-125.
	Starting with (x 0 , y 0 ) = (0.88187777591, 0.0000322222356) one obtains
	(x 686,295,403,186 , y 686,295,403,186 ) = (x 250,125,214,227 , y 250,125,214,227 )
	= (1.34348444450739479, -2.07670041497986548. 10 -7
	11

5 :

 5 In this case also,

	Published In Mathematics in Science and Technology: Mathematical Methods, Models and
	Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August
	2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher,
	Singapore (2011) pp. 93-125.
	iter NSampl is not known a priori, it is more
	complicated to give an equivalent to it. However, considering that the selecting

  Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August 2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher, Singapore (2011) pp. 93-125.

	22

  Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August 2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher, Singapore (2011) pp. 93-125.

	Published In
	1,E+06
	1,E+05
	1,E+04
	1,E+03
	1,E+02		Window of emergence of
			randomness
	1,E+01
	1,E+00
	1,E-01
	1,E-02
	1,E-03
	1,E-04
	1,E-02	1,E-05	1,E-08	1,E-11	1,E-14	1,E-17
			ε
		E A C Sup	Esup
		E A C 1		E1
		E A C 2		E2

Table 3 .

 3 Published In Mathematics in Science and Technology: Mathematical Methods, Models and Algorithms in Science and Technology. Proceedings of the Satellite Conference of ICM, 15-17 August 2010 Delhi, India (A.H. Siddiqi, R.C. Singh, P. Manchanda Editors), World Scientific Publisher, Singapore (2011) pp. 93-125. Sequences ( )

	ε	1	0.000, 001	ε	* 1	0.000,001, 000, 000,000, 000,3
	1 0 x	0.330,000, 013, 113,021,851	*1 0 x	0.330,000, 013, 113,021,851
	( 0 ) 1 n x	0.959, 214,817, 207,605,153 -	( 0 ) *1 n x	-0.058,536,729, 173,974,455,5
	( 1 ) 1 n x	0.657,775,688, 600,752,417	( 1 ) *1 n x	0.386,129,403, 866,398,935
	( 2 ) 1 n x	-0.784,600,935, 471,051,031	( 2 ) *1 n x	0.471,824,729, 381,262,631

 ( ( ( ( ) )) )

x , x , x , , x , x ,

) )

x , x , x , , x , x ,

)

x , x , x , , x , x ,

… in order to distribute the iterated points with respect to this given partition defining the subsequence ( )

, , , , , ,

)

x , x , x , ,x , x , 11) becomes

Moreover we assume that i

For the shake of simplicity we consider only the chaotic sampling method (i.e.

we use only one threshold T ), without the chaotic mixing. We then compute 1, , 0.000, 001, 000, 000, 000, 000,3 ε = , and 2 0.000, 002

The double threshold sampling is done using 0.9 T = and ' 20 T = and the same seed is taken ( ) ( )

Despite the fact that the difference between 1 ) )

x , x , x , ,x , x ,

) )

x , x , x , ,x , x ,


Emmanuelle Comets

PhD Karl Brendel

PhD France Mentré

Why should prediction discrepancies be renamed Standardized VPC ?

Keywords: non-linear mixed effect models, model evaluation, visual predictive check, prediction discrepancies, normalised prediction distribution errors

We have read the paper by Wang and Zhang [START_REF] Wang | Standardized Visual Predictive Check versus Visual Predictive Check for model evaluation[END_REF] . We are pleased to see that other authors further illustrate and also find nice properties for an evaluation tool that we first used more than ten years ago [START_REF] Mesnil | Population pharmacokinetics analysis of mizolastine and validation from sparse data on patients using the nonparametric maximum likelihood method[END_REF] , presented at PAGE in 2000 [START_REF] Mentré | Validation methods in population pharmacokinetics: a new approach based on predictive distributions with an evaluation by simulation[END_REF] and published extensively first in 2006 [START_REF] Mentré | Prediction discrepancies for the evaluation of nonlinear mixedeffects models[END_REF] . However we disagree with some of the conclusions made by the authors, and we feel that an important reference should have been included in their manuscript.

As Wang and Zhang point out themselves, what they call SVPC is nothing other than the prediction discrepancies (pd) named that way, after being called pseudo-residuals by Mentré and Escolano in 2006 [START_REF] Mentré | Prediction discrepancies for the evaluation of nonlinear mixedeffects models[END_REF] . It is therefore misleading to present SVPC as something novel when in fact it goes back on something that our group has published and presented in conferences.

Plots of SVPC versus time in the paper by Wang and Zhang 1 are very similar to plots of pd versus time in Mentré and Escolano [START_REF] Mentré | Prediction discrepancies for the evaluation of nonlinear mixedeffects models[END_REF] , therefore we disagree with the statement on page 2 that "Neither pd nor npde was intended/recommended for evaluation of model predictions over a time course". pd and npde were developed for their improved statistical properties over linearisationbased residuals but are used as visual diagnostic tools in a similar way. Also the paper is not complete when it comes to the current state of literature. Wang and Zang 1 missed one important reference, in which we have compared pd, npde, VPC, as well as tests based on prediction intervals with or without decorrelation [START_REF] Brendel | Evaluation of different tests based on observations for external model evaluation of population analyses[END_REF] . We have also proposed tests for covariate models illustrating how npde can be used to evaluate covariate models. Their paper is rather inaccurate when comparing the properties of the pd to those of the npde. The idea behind the npde, as recalled in the paper, is to decorrelate observations within individuals. In Mentré and Escolano [START_REF] Mentré | Prediction discrepancies for the evaluation of nonlinear mixedeffects models[END_REF] , it had been shown that the type I error of the test based on pd was inflated when subjects contribute several observations, because of the within-subject correlation between observations, and it was anticipated that the decorrelation would improve this feature. Indeed, in Brendel et al. [START_REF] Brendel | Evaluation of different tests based on observations for external model evaluation of population analyses[END_REF] we have shown that the type I error of the npde is close to the expected 5% while the type I error for the pd or the test based on prediction intervals of the VPC was larger. In that paper, we performed a full simulation study, i.e. with several replications of the simulated data set, while in the manuscript of Wang and Zhang only one simulated dataset is given, from which it is very hard to draw meaningful conclusions. It is quite possible that for one simulated dataset the pd detects something that the npde does not, indeed, our simulations have shown that on a large number of simulated datasets, the tests with npde maintain type I error while the tests with pd have an increased type I error and hence a higher power. A discussion of the power to detect model misspecification should therefore take into account the increase in type I error under the null hypothesis. Also in their table V, there was a significant departure from 0 of the mean npde for the wrong model (p=0.03 with a Wilcoxon test).

We would also like to take this opportunity to mention that more informative VPC graphs have been proposed by Wilkins 6 et al in 2006, which use prediction intervals around the simulated percentiles, and that we have adapted these graphs to pd and npde [START_REF] Comets | Model evaluation in nonlinear mixed effect models, with applications to pharmacokinetics[END_REF] . In that paper, we made a case of using pd instead of npde to plot diagnostic graphs because the decorrelation tends to blur the relationship with time when used for visual diagnostics.

Finally, we join Wang and Zhang to stress the nice properties of simulation-based model evaluation tools such as pd (SVPC) over VPC and to encourage readers to use them.


Viet-Trung Tran
email: viet-trung.tran@irisa.fr

Gabriel Antoniu
email: gabriel.antoniu@inria.fr

Towards a storage backend optimized for atomic MPI-I/O for parallel scientific applications

I. INTRODUCTION

Scientific applications are becoming increasingly dataintensive: high-resolution simulations of natural phenomenas, climate modeling, large-scale image analysis, etc. Such applications currently manipulate data volumes in the petabyte scale and with the growing trend of data sizes we are rapidly advancing towards the exabyte scale. In this context, I/O performance has been repeatedly pointed out as a source of bottleneck that negatively impacts the performance of the applications. One of problems that is causing the I/O bottleneck are the I/O access patterns generated by such scientific applications that do not match the I/O access interfaces exposed by the file systems used as underlying storage back-ends.

One particularly difficult challenge in this context is the need to efficiently address the I/O needs of scientific applications [START_REF] Smirni | I/O requirements of scientific applications: An evolutionary view[END_REF], [START_REF] Smirni | Lessons from characterizing the I/O behavior of parallel scientific applications[END_REF] that partition multi-dimensional domains into overlapping subdomains. The subdomains need to be processed in parallel and then stored in a globally shared file. Since the file is a flat sequence of bytes, subdomains map to non-contiguous regions in the file. Because the subdomains overlap, under concurrent accesses such non-contiguous regions may interleave in an inconsistent fashion if they are not grouped together as a single atomic transaction. Therefore, atomicity of non-contiguous, overlapping reads and writes of data from a shared file is a crucial issue.

The purpose of this Phd research aims at addressing this shortcoming of existing approaches by optimizing the storage back-end specifically for the access patterns described above.

II. PROBLEM DESCRIPTION

In a large class of scientific applications, especially largescale simulations, input and output data represents huge spatial domains made of billions of cells associated with a set of parameters (e.g., temperature, pressure, etc.). These spatial domains represent the state of the simulated system at a specific moment in time. They are iteratively refined by the simulation in order to obtain an insight into how the system evolves in time. The large size of the spatial domains triggers the need to parallelize the simulation, by splitting them into subdomains that are distributed and processed by a large number of compute elements, typically MPI processes.

In many such simulations, the contents of a cell depends on the contents of the neighboring cells. Thus, the cells at the border of a subdomain (called "ghost cells") need to be shared by more than one MPI process. In order to avoid repeated exchanges of border cells between MPI processes during the simulation, a large class of applications [START_REF] Smirni | I/O requirements of scientific applications: An evolutionary view[END_REF], [START_REF] Smirni | Lessons from characterizing the I/O behavior of parallel scientific applications[END_REF], [START_REF] Crandall | Input/Output characteristics of scalable parallel applications[END_REF], [START_REF] Nieuwejaar | File-access characteristics of parallel scientific workloads[END_REF] partition the spatial domain in such way that the resulting subdomains overlap at their borders.

At each iteration, the MPI processes typically dump their subdomains in parallel in a globally shared file, which is then later used to interpret and/or visualize the results. Since the spatial domain is a multidimensional structure which is stored as a single, flat sequence of bytes, the data corresponding to each subdomain maps to a set of noncontiguous regions within the file. This in turn translates to a write-intensive access pattern where the MPI processes concurrently write a set of non-contiguous regions in the same file. Moreover, because some subdomains overlap, the non-contiguous regions in the file belonging to such subdomains may overlap too.

In order to obtain a globally shared file that represents the whole spatial domain in a consistent way, it is important to write all non-contiguous regions belonging to the same subdomain in an atomic fashion. This atomicity requirement is defined in MPI standardization as a guarantee that in concurrent, overlapping MPI-I/O write operations (which can possibly involve sets of non-contigous regions), the overlapped regions shall contain data only from one of the MPI processes that participates in the concurrent writes.

However, this work is a non-trivial issue, and that has been difficult to be implemented. Currently, MPI-I/O implementations like ROMIO [START_REF] Thakur | On implementing MPI-I/O portably and with high performance[END_REF] poorly support MPI atomic mode because of the poor implementations of storage back-ends.

Either they do not offer guarantees with respect to atomicity, which means concurrent writes need to be coordinated at application level. Or they offer limited guarantees, such as the POSIX atomicity semantics [START_REF]Information technology -Portable Operating System Interface (POSIX)[END_REF], where a write of a contiguous region is guaranteed to end up as a single sequence of bytes, without interleaving with other writes. However, in our context, the POSIX atomicity guarantees are not sufficient, because a single write operation involves a whole set of non-contiguous regions. Under concurrency, this can lead to an interleaved overlapping that generates inconsistent states.

III. RELATED WORK

Previous work has shown that providing MPI atomicity efficiently enough is not a trivial task in practice, especially when dealing with concurrent, non-contiguous I/O, most of which rely on locking-based techniques. Several approaches have been proposed at various levels: at the level of the MPI-I/O layer, at the level of the file system.

A first series of approaches assumes no specific support at the level of the parallel file system This is typically the case of PVFS [START_REF] Haddad | PVFS: A parallel virtual file system for Linux clusters[END_REF], a widely-used parallel file system which makes the choice of enabling high-performance data access for both contiguous and non-contiguous I/O operations without guaranteeing atomicity at all for I/O operations. For applications where MPI atomicity requirement needs to be satisfied, a solution (e.g. illustrated in [START_REF] Ross | Implementing MPI-I/O atomic mode without file system support[END_REF]) consists in guaranteeing MPI atomicity portably at the level of the MPI-I/O layer. The lack of guarantees on the semantics of I/O operations provided by the file system comes however at a high cost introduced by the use of coarse-grain locking at a higher level. Typically, the whole file is locked for each I/O request and thus concurrent accesses are serialized, which is an obvious source of overhead. To avoid this bottleneck in the case of concurrent non-overlapping accesses to the same shared file, an alternative approach [START_REF] Sehrish | Conflict detection algorithm to minimize locking for MPI-I/O atomicity[END_REF] proposes to introduce a mechanism for automatic detection of non-overlapping I/O and thus avoid locking in this case. However, as acknowledged by the authors of this approach, an unnecessary overhead due to the detection mechanism is then introduced for non-overlapping concurrent I/O.

Another class of approches addresses the case where the underlying parallel file system supports POSIX atomicity. Atomic contiguous I/O can then seamlessly be mapped to atomic read/write primitives provided by the POSIX interface. However, POSIX atomicity alone is not enough to provide the necessary atomicity guarantees for applications that exhibit concurrent, non-contiguous I/O operations. It is important to note that, to enable MPI atomicity, both contiguous and non-contiguous I/O requests need to be considered.

Parallel file systems such as GPFS [START_REF] Schmuck | GPFS: A shared-disk file system for large computing clusters[END_REF] and Lustre [START_REF] Schwan | Lustre: Building a file system for 1000-node clusters[END_REF] provide POSIX atomicity semantics using a distributed locking approach: locks are stored and managed on the storage servers hosting the objects they control. Whereas POSIX atomicity can simply and directly be leveraged for contiguous I/O operations using byte range locking, enabling atomic non-contiguous I/O based on POSIX atomicity is not efficient. In the default scheme, if we consider a set of non-contiguous byte ranges to be atomically accessed by an individual I/O request, it is then necessary to lock the smallest contiguous byte range that covers all elements of the set of ranges to be accessed. This leads to unnecessary synchronization and thus to a potential bottleneck, since this contiguous byte range also covers unaccessed data that would not need to be locked.

IV. OUR APPROACH

We propose a general approach to solve the issue of enabling a high throughput under concurrency for writes of non-contiguous, overlapped regions under MPI atomicity guarantees. This approach relies on three key design principles:

Dedicated API at the level of the storage back-end: Traditional approaches address the issue of providing support for MPI atomic, non-contiguous writes by implementing the MPI-I/O access interface as a layer on top of highly standardized consistency semantics models (e.g. POSIX). The rationale behind this approach is to be able to easily plug in a different storage back-end without the need to rewrite the MPI-I/O layer. However, this advantage comes at a high price: the MPI-I/O layer needs to translate the MPI atomicity into a different consistency model, which greatly limits the potential to optimize for the access patterns that are needed in our context. By contrast, we propose to extend the storage back-end with a data access interface that provides native support for non-contiguous, MPI-atomic data accesses. Using this approach circumvents the need to translate to a different consistency model and enables introducing optimizations directly at the level of the storage back-end, enhancing the potential to implement a a better concurrency control scheme.

Data striping: Data striping a well-known technique to increase the throughput of data accesses, by splitting the huge file where the spatial domain is stored into chunks that distributed among multiple storage elements. Using a loadbalancing allocation strategy that redirects write operations to different storage elements in a round robin fashion, the I/O workload is distributed by itself, which effectively increases the overall throughput that can be achieved.

Versioning as a key to enhance data access under concurrency: Most storage back-ends manipulate a single version of the file at a time under concurrency by providing locking-based mechanisms that guarantee exclusive access to overlapped regions, in order to eliminate inconsistencies. However, in our context, such an approach does not scale, even if the storage back-end is able to deliver high throughputs. The problem comes from the fact that a locking-based scheme is expensive, as it enables only one writer at a time to gain exclusive access to a region. Since there are many overlapped regions, this leads to a situation where many writers are idle while waiting for their turn to lock, which greatly limits the potential to achieve a globally high aggregated throughput. In order to avoid this issue, we propose a versioning-based access scheme that avoids the need to lock, effectively eliminating the need of writers to wait for each other.

Our approach is based on shadowing techniques [START_REF] Rodeh | B-trees, shadowing, and clones[END_REF], which offers the illusion of creating a new standalone snapshot of the file for each update to it but to physically store only the differences and manipulate metadata in such way that the aforementioned illusion is upheld. Starting from the principles introduced in [START_REF] Nicolae | BlobSeer: Next-generation data management for large scale infrastructures[END_REF], we propose to enable concurrent MPI processes to write their non-contiguous regions in complete isolation, without having to care about overlappings and synchronization. This is made possible by keeping data immutable: new differences are added instead of modifying an existing snapshot. It is at the metadata level where the ordering is done and the overlappings are resolved in such way as to expose a snapshot of the file that looks as if all differences were applied in an arbitrary sequential order.

V. IMPLEMENTATION

Following the approach proposed in Section IV, we have chosen to build the storage back-end on top of BlobSeer [START_REF] Nicolae | BlobSeer: Next-generation data management for large scale infrastructures[END_REF], a versioning-oriented data sharing service specifically designed to meet the requirements of data-intensive applications that are distributed at large scale: scalable aggregation of storage space from a large number of participating machines with minimal overhead, support to store huge data objects, efficient fine-grain access to data subsets and ability to sustain a high throughput under heavy access concurrency.

The choice of building the storage back-end on top of BlobSeer was motivated by two factors. First, BlobSeer supports transparent striping of BLOBs into chunks and enables fine-grain access to them, which enables to store each spatial domain directly as a BLOB. This in turn avoids the need to perform data striping explicitly. Second, BlobSeer offers out-of-the-box support for shadowing by generating a new BLOB snapshot for each fine-grain update while physically storing only the differences. This provides a solid foundation to introduce versioning as a key principle to support MPI atomicity.

Since BlobSeer supports atomic reads and writes of contiguous regions only, we can not leverage its versioningoriented access interface directly to support MPI atomicity which consists of both atomic contiguous and noncontiguous accesses. As mentioned in Section IV, we want to avoid locking and introduce optimizations directly at the level of the storage back-end. These following tunings were done inside BlobSeer in order to support MPI atomicity:

• We extended the access interface of BlobSeer such that it can describe complex non-contiguous data access in a single call. We introduce a series of versioningoriented primitives that facilitate non-contiguous accesses. These primitives are closely matched the List I/O interface proposal introduced in [START_REF] Ching | Noncontiguous I/O through PVFS[END_REF]. • We added support for atomicity of non-contiguous I/O under concurrency by hacking inside BlobSeer's shadowing technique. We built the metadata corresponding to update made by non-contiguous writes such that only consistent snapshots that obey MPI atomicity are published. A detail implementation of our prototype is described in [START_REF] Tran | Efficient support for MPI-IO atomicity based on versioning[END_REF]. To leverage our storage back-end at the level of the MPI-I/O layer, we integrated our prototype into ROMIO [START_REF] Thakur | On implementing MPI-I/O portably and with high performance[END_REF], a popular MPI-I/O implementation. Since our prototype provides native support for both contiguous and non-contiguous writes under MPI atomicity, ROMIO thus has no need to translate between consistency model, which was source of overhead that impacts the performance of the applications.

VI. PRELIMINARY RESULTS

We targeted to build a storage back-end that explicitly optimizes for non-contiguous, overlapped I/O accesses that obey MPI atomicity guarantees. To evaluate our prototype, we compare our approach to the locking-based approach that leverages a POSIX-compatible file system at the level of the MPI-I/O layer, which is the traditional way of addressing MPI atomicity. The Lustre parallel file system [START_REF] Schwan | Lustre: Building a file system for 1000-node clusters[END_REF] was chosen against our prototype. Both Lustre and our prototype are plugged into ROMIO by using their default ROMIO ADIO modules.

Two series of experiments have been realized. Our first experiment aims at evaluating the scalability of our approach when increasing the number of clients that concurrently write non-contigous regions into the same file. In this scenario, we consider the extreme case where each of the clients writes a large set of non-contiguous regions that are intentionally selected in such way as to generate a large number of overlapping that need to obey MPI atomicity. In the second experiment, we performed an evaluation of the performance of our approach using a standard benchmark, MPI-tile-IO, that closely simulates the access patterns of real scientific applications that split the input data into overlapped subdomains that need to be concurrently written in the same file under MPI atomicity guarantees.

The preliminary experiments are performed on the Grid'5000 [START_REF] Jégou | Grid'5000: a large scale and highly reconfigurable experimental grid testbed[END_REF] testbed, a reconfigurable, controllable and monitor-able experimental Grid platform gathering 9 sites geographically distributed in France. The results have been reported in [START_REF] Tran | Efficient support for MPI-IO atomicity based on versioning[END_REF]. Our approach demonstrated excellent scalability under concurrency when compared to the Lustre-based approach. It achieved an aggregated throughput ranging from 3.5 times to 10 times higher in several experimental setups, including highly standardized MPI benchmarks specifically designed to measure the performance of MPI-I/O for non-contiguous overlapped writes that need to obey MPI-atomicity semantics.

VII. CONCLUSIONS

We have proposed an original versioning-based mechanism that can be leveraged to efficiently address the I/O needs of data-intensive MPI applications involving data-partitioning schemes that exhibit overlapping noncontiguous I/O where MPI-I/O atomicity needs to be guaranteed under concurrency.

Unlike traditional approaches that leverage POSIXcompliant parallel file systems as storage back-ends and employ locking schemes at the level of the MPI-I/O layer, we propose to use versioning techniques as a key principle to achieve high throughput under concurrency while guaranteeing MPI atomicity. We implemented this idea in practice by extending BlobSeer, an existing versioning-oriented, distributed data sharing service, with a non-contiguous data access interface that we directly integrated with ROMIO, a standard MPI-I/O implementation. A number of experiments were performed which demonstrated our approach having excellent scalability under concurrency when compared to locking-based approach.

Such promising results encouraged us to purse interesting future work directions. In particular, one advantage of Blob-Seer that we have not developed is to expose its versioning interface directly at application level. The ability to make use of versioning at application level brings several potential benefits, such as in the case of producer-consumer workloads where for example the output of simulations is concurrently used as the input of visualizations. Using versioning at application level could avoid expensive synchronization schemes, which is an acknowledged problem of current approaches.