[R-pkg-devel] For reproducibility issue

Martin Maechler m@ech|er @end|ng |rom @t@t@m@th@ethz@ch
Mon Jan 20 09:46:37 CET 2020


>>>>> وليد خلف معوض المطيرى 
>>>>>     on Fri, 17 Jan 2020 19:10:53 +0000 writes:

    > Hi all,

    > I think what I�ve done is something different.

yes, but ...

    > Inside the Fortran subroutine, I have a subroutine for setting the seed value of the RNG of GNU Fortran which I think it is not related to the R RNG like the one below:

    > subroutine initrandomseedsinr(temp)
    > implicit none
    > integer :: n
    > integer, intent(in):: temp
    > integer, dimension(:), allocatable :: seed

    > call random_seed(size = n)
    > allocate(seed(n))
    > seed = temp
    > call random_seed(PUT = seed)
    > deallocate(seed)

    > end subroutine initrandomseedsinr
    > , 
    > where temp is an argument of the Fortran subroutine as well as in the wrapper R function. This will be related to the RNG method used in the GNU Fortran that build on GCC not to the R. I am not sure if I am right on this, but tried with using RNGkind(sample.kind = "Rounding") and it doesn�t help.

indeed, it has no effect on your Fortran / run time random
generator.

But really really truly :  a "decent" R package should use R's
random number generators rather than (C library / Fortran
library / ...).
When using R's RNGs, your package behaves "R like" e.g.,
set.seed() has the desired effect, etc etc.

What Avi Adler describes below is what I strongly believe you
*really* should change your package to do as well.

Best regards,

Martin Maechler
ETH Zurich and R Core Team


    > The difference in the results were not major. The output at the end of running the functions kept having very similar results, but still have the issue of reproducing exact results which I need it for relating work that is based on the package.

    > Many thanks,

    > Waleed




    > Sent from Mail<https://go.microsoft.com/fwlink/?LinkId=550986> for Windows 10

    > From: Avraham Adler<mailto:avraham.adler using gmail.com>
    > Sent: �Friday,� �January �17,� �2020 �06:30 ��
    > To: Ivan Krylov<mailto:krylov.r00t using gmail.com>
    > Cc: ���� ��� ���� �������<mailto:wkmtierie using qu.edu.sa>; R Package Development<mailto:r-package-devel using r-project.org>
    > Subject: Re: [R-pkg-devel] For reproducibility issue

    > Hi.

    > If it helps, I call the R RNG from Fortran in my Delaporte package
    > [1], also using iso_c_bindings. Specifically, I have the following C
    > code [2]:

    > void F77_SUB(unifrnd) (int *n, double *x){
    > GetRNGstate();
    > for (int i = 0; i < *n; ++i){
    > *(x + i) = unif_rand();
    > }
    > PutRNGstate();
    > }
    > and call it in Fortran like so [3]:

    > subroutine rdelap_f(n, a, na, b, nb, l, nl, vars) bind(C, name="rdelap_f_")

    > external unifrnd

    > integer(kind = c_int), intent(in), value :: n, na, nb, nl
    > real(kind = c_double), intent(out), dimension(n) :: vars
    > real(kind = c_double), intent(in) :: a(na), b(nb), l(nl)
    > real(kind = c_double), dimension(n) :: p
    > integer(kind = c_int) :: lg, lt

    > call unifrnd(n, p)
    > lt = 1_c_int
    > lg = 0_c_int
    > call qdelap_f(p, n, a, na, b, nb, l, nl, lt, lg, vars)

    > end subroutine rdelap_f

    > The package passes CRAN tests (at least as of now) on anything between
    > GCC 4 and GCC9 [4].

    > Hope that helps,

    > Avi

    > [1] https://bitbucket.org/aadler/delaporte/src/master/
    > [2] https://bitbucket.org/aadler/delaporte/src/master/src/utils_and_wrappers.c
    > [3] https://bitbucket.org/aadler/delaporte/src/master/src/delaporte.f95
    > [4] https://cran.r-project.org/web/checks/check_results_Delaporte.html


    > On Fri, Jan 17, 2020 at 2:39 PM Ivan Krylov <krylov.r00t using gmail.com> wrote:
    >> 
    >> On Fri, 17 Jan 2020 13:55:39 +0000
    >> ���� ��� ���� ������� <wkmtierie using qu.edu.sa> wrote:
    >> 
    >> > So, does anyone have an idea of how to solve this issue.
    >> 
    >> "Writing R Extensions", 1.6. Writing portable packages:
    >> 
    >> >> Compiled code should not call the system random number generators
    >> >> such as rand, drand48 and random, but rather use the interfaces to
    >> >> R�s RNGs described in Random numbers. In particular, if more than
    >> >> one package initializes the system RNG (e.g. via srand), they will
    >> >> interfere with each other.
    >> 
    >> >> Nor should the C++11 random number library be used, nor any other
    >> >> third-party random number generators such as those in GSL.
    >> 
    >> It somewhat less convenient to call the R random number generator from
    >> Fortran than it would be from C or C++, but still possible. There is a
    >> F77-style example of such use [1], but since you are already using
    >> iso_c_binding, you should be able to declare the C API [2] right in the
    >> Fortran source:
    >> 
    >> subroutine GetRNGState() bind(c)
    >> end subroutine
    >> 
    >> subroutine PutRNGstate() bind(c)
    >> end subroutine
    >> 
    >> As a bonus, you get to use the R distribution functions [3], without
    >> the need to implement them yourself from uniformly distributed samples:
    >> 
    >> function rnorm(mu, sigma) bind(c) result(ret)
    >> use intrinsic, iso_c_binding, only: c_double
    >> real(c_double), value :: mu, sigma
    >> real(c_double) ret
    >> end function
    >> 
    >> function rgamma(shape, scale) bind(c) result(ret)
    >> use intrinsic, iso_c_binding, only: c_double
    >> real(c_double), value :: shape, scale
    >> real(c_double) ret
    >> end function
    >> 
    >> (The prototypes above are unchecked; I haven't written any Fortran 2003
    >> in more than a year.)
    >> 
    >> --
    >> Best regards,
    >> Ivan
    >> 
    >> [1]:
    >> https://cran.r-project.org/doc/manuals/R-exts.html#index-Random-numbers-in-Fortran
    >> 
    >> [2]: https://cran.r-project.org/doc/manuals/R-exts.html#Random-numbers
    >> 
    >> [3]:
    >> https://cran.r-project.org/doc/manuals/R-exts.html#Distribution-functions
    >> 
    >> ______________________________________________
    >> R-package-devel using r-project.org mailing list
    >> https://stat.ethz.ch/mailman/listinfo/r-package-devel


    > [[alternative HTML version deleted]]

    > ______________________________________________
    > R-package-devel using r-project.org mailing list
    > https://stat.ethz.ch/mailman/listinfo/r-package-devel



More information about the R-package-devel mailing list