The Coq mailing list

[Fontaine Allyx <fontaine AT labri.fr>]

Le 24/05/2011 13:15, AUGER Cedric a écrit :
> Le Tue, 24 May 2011 10:50:44 +0200,
> Fontaine 
> Allyx<fontaine AT labri.fr>
>   a écrit :
>
> > Le 23/05/2011 11:18, AUGER Cedric a écrit :
> > > Le Sun, 22 May 2011 17:41:50 +0200,
> > > David 
> > > Baelde<david.baelde AT gmail.com>
> > >    a écrit :
> > >
> > > > Hi Allyx,
> > > >
> > > > If you haven't already, you should have a look at the paper by
> > > > Audebaud&   Paulin (http://dx.doi.org/10.1016/j.scico.2007.09.002)
> > > > about the ALEA library. It starts with a gentle introduction to the
> > > > monad that is used, detailing its advantages and drawbacks.
> ...
> > > > Long story short, Pierre Corbineau's solution will work. However, I
> > > > wouldn't bother to make extraction work for such small examples:
> > > > just rewrite the code in ML yourself.
> > > >
> > > > Have fun,
> > >
> > > I also guess that if "'a" is a type which cardinality is quite
> > > small, and you can extract [0,1] to float, then you can write a
> > > short ML code:
> > >
> > > let dice_throw =
> > >     let sieve = Random.float 1. in
> > >     let rec get_rand acc l =
> > >       match l with
> > >       | v::l ->   let acc = acc + μ v in
> > >                 if acc>   sieve then v else get_rand acc l
> > >       | [] ->   failwith "Hmm, this is unexpected!"
> > >     in get_rand 0 [list of all the elements of your type]
> > >
> > > of course if you have big types for 'a such as int it can be too
> > > long and even false due to rounding errors (the failwith case can
> > > even occur).
> > > You also have to identify the "μ" function (if it is the one I think
> > > of, it is also called μ in Alea, so it should be extracted to its
> > > Unicode code, as unicode is not supported by OCaml, leading to some
> > > awfull identifier).
> > > Even in this case, I am not sure it can be very usefull, but if it
> > > is, then you won't have to rewrite all the ML code, and won't have
> > > to play with the Extraction mechanism too much (and using P.
> > > Corbineau solution seems also hazardous, as you change the
> > > extracted objects to objects which are semantically very different;
> > > furthermore, using side effects can mess a lot with extraction if
> > > some optimizations are performed ("let x = side_effect () in let y
> > > = side_effect () in ..." could be extracted to "let x = side_effect
> > > () in let y = x in ..." as well as identity due to Coq semantics,
> > > but the two functions are clearly different with OCaml)).
> >
> > Thanks for your answers!
> >
> > Actually, as I have several examples more or less simple, I would
> > like to have an automatic way which ensures that Ocaml functions
> > really corresponds to Coq functions. But as I don't know well all
> > those mechanisms, I indeed begin with writting ML code by myself and
> > having a look at extraction.
> >
> > I am not sure of my understanding of the role of the accumulator in
> > dice_throw. Is it needed to return an element of my type uniformly at
> > random thanks to the μ function? Does the failwith case not be
> > reached because the sum of the probabilities is supposed to be
> > greater than 1?
> >
> > Allyx
>
> Almost yes for the last paragraph (in fact the sum is expected to be
> exactly 1, and not "greater" than 1; in french "greater" has a large
> meaning, but in english it has a strict one).
>
> In fact I never used the Alea library, so what I wrote may be wrong,
> but if you want some commented version of the function I wrote:
> (with the notation:
> Int(a,b,P) for:
>   b
> ∫ P(x) ∂x
>   a
> and:
> f^(-1) for:
> f⁻¹
> )
>
>   let dice_throw =
>      (*First we throw a dice on [0,1] with uniform probability,
>        assuming [Random.float] is able to do so,
>        which is not the case in fact as it is a pseudo-random
>        generator and float doesn't cover all reals in [0,1],
>        but it is expected to be a good approximation of this law.
>        Let us bind this result to the "sieve" identifier.
>      *)
>      let sieve = Random.float 1. in
>      (* Now thanks to theorems which I do not remember the names,
>         we have some transfer function between uniform law and
>         any Riemann integrable law P on [a, b] such that Int(a,b,P)=1
>         (I don't recall if it can be generalized to Lebesgue,
>          but that doesn't matter here) with:
>          X |-->  (x |-->  Int(a, x, P))^(-1)(X)
>         So that now we want to compute "(x |-->  Int(a, x, P))^(-1)"
>      *)
>      let rec get_rand acc l =
>        (* [acc] is meant to be "Int(a, predecessor(List.hd l), P)"
>           where [a] is the first element of the type
>           (in fact as the type is finite, the integral is just
>            a finite sum)
>        *)
>        match l with
>        | v::l ->   let acc = acc + μ v in (*μ is our P law*)
>                  if acc>   sieve
>                  then v (*ok, we have reversed our integral*)
>                  else get_rand acc l
>                   (*we didn't reverse our integral and had to go
>                     further*)
>        | [] ->  (*if this point is reached, then by definition
>                  of probability laws, we should have that [acc]=1.
>                  and so, due to definition of sieve we should have
>                  halted earlier. Of course due to rounding,
>                  we could have [acc]=0.99988665 or so and
>                  [sieve] = 0.9999039 or so, so that this case
>                  may happen if we are very unlucky
>                *)
>                failwith "Hmm, this is unexpected!"
>     in
>     get_rand 0 [list of all the elements of your type]
>
>
> Hope these comments are not too messy!
> Once again, I guess it is what you would like, but I am not sure it is
> very convenient. And if you want to do big benchmarks to be able to
> make some observations, it may be quite slow. But I have no really
> better solution.

Thanks a lot for the explanations! It's much more clear for me now.
Allyx