[R-pkg-devel] For reproducibility issue

Hi all,

I think what I?ve done is something different. 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. 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 at gmail.com>
Sent: ?Friday,? ?January ?17,? ?2020 ?06:30 ??
To: Ivan Krylov<mailto:krylov.r00t at gmail.com>
Cc: ???? ??? ???? ???????<mailto:wkmtierie at qu.edu.sa>; R Package Development<mailto:r-package-devel at 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, 17 Jan 2020 13:55:39 +0000
???? ??? ???? ??????? <wkmtierie at 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 at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-package-devel