Issues with setting random seed

Question

I am attempting to write a Montecarlo algorithm to simulate interaction between agents in a population. This algorithm needs to call two random numbers at each iteration (say, 10^9 iterations).

My issue here is that everytime I change the seed (to obtain different realizations), the RNG is giving me the same output (same Montecarlo events). I have tried different ways of implementing it with to no avail.

I am completely new to Fortran and copying this code from MATLAB. Am I doing something wrong in the way I'm implementing this code?

Below is what I tried:

program Gillespie

  implicit none

  integer*8, parameter :: n_max = 10.0**8 ! max. number of iterations
  integer*8 :: t_ext, I_init, S_init, jump, S_now, I_now, i, u
  real*8 :: t, N, a0, tau, st, r1, r2

  real, parameter :: beta = 1000
  real, parameter :: gammma = 99.98
  real, parameter :: mu = 0.02
  real, parameter :: R0 = beta/(gammma+mu)
  integer :: seed = 11

  real, dimension(n_max) :: S_new ! susceptible pop. array
  real, dimension(n_max) :: I_new ! infected pop. array
  real, dimension(n_max) :: t_new ! time array
  real, dimension(5) :: events    ! events array

  open(unit=3, file='SIS_output.dat')

  t = 0                     ! initial time
  N = 40                    ! initial population size

  jump = 1                  ! time increment (save in uniform increments)
  u = 2

  t_ext = 0                 ! extiction time
  I_init = 2                ! initial infected pop.
  S_init = N-I_init         ! initial susceptible pop.

  S_now = S_init
  I_now = I_init

  S_new(1) = S_init         ! initialize susceptibles array
  I_new(1) = I_init         ! initialize infected array
  t_new(1) = t              ! initialize time array

  write(3,*) t_new(1), S_new(1), I_new(1) ! write i.c. to array 

  call random_seed(seed)

  do i=2, n_max
     call random_number(r1)
     call random_number(r2)

     events(1) = mu*N       ! Birth(S)
     events(2) = mu*S_now   ! Death(S)
     events(3) = mu*I_now   ! Death(I)
     events(4) = (beta*S_now*I_now)/N ! Infection
     events(5) = gammma*I_now ! Recovery

     a0 = events(1)+events(2)+events(3)+events(4)+events(5)

     tau = LOG(1/r1)*(1/a0) ! time increment
     t = t + tau            ! update time
     st = r2*a0             ! stochastic time???

     ! update the populations
     if (st .le. events(1)) then
        S_now = S_now + 1
     else if (st .gt. events(1) .AND. st .le. 
 #(events(1) + events(2)))  then
        S_now = S_now - 1
     else if (st .gt. (events(1) + events(2)) .AND. st .le. 
 #(events(1) + events(2) + events(3))) then
        I_now = I_now - 1
     else if (st .gt. (events(1) + events(2) + events(3)) .AND. 
 #st .le. (events(1) + events(2) + events(3) + events(4))) then
        S_now = S_now - 1
        I_now = I_now + 1
     else
        S_now = S_now + 1
        I_now = I_now - 1
     end if

     ! save time in uniform increments
     if(t .ge. jump) then
         t_new(u) = floor(t)
         S_new(u) = S_now
         I_new(u) = I_now

         write(3,*) t_new(u), S_new(u), I_new(u)

         jump = jump+1
         u = u+1
     end if

    ! write(3,*) t_new(i), S_new(i), I_new(i)

    !N = S_now + I_now   ! update population post event

     if(I_now .le. 0) then  ! if extinct, terminate
        print *, "extinct"
        goto 2
     end if 

   end do

 2     end program Gillespie

I appreciate all input. Thanks.

Hi, I think you probably need to do some more steps to set the seed for Fortran built-in RNG (e.g., this page http://stackoverflow.com/questions/21255631/gfortran-and-random-numbers ) In short, get the size of "seed array" first, allocate the seed array with that size, populate that array somehow (e.g. with system_clock()), and give the seed array to random_seed() with the "PUT" keyword (so a rather tedious procedure...) — roygvib, Aug 09 '16 at 22:58
Depending on your compiler `call random_seed()` may give you a non-repeatable sequence. Notably, gfortran won't. — francescalus, Aug 11 '16 at 18:47
Ah, thanks @VladimirF, good to know. It's one of those processor dependent things often catching people out. — francescalus, Aug 11 '16 at 18:50
@francescalus There is a very long thread on the dev site which I wasn't able t follow completely, but my impression is like that https://gcc.gnu.org/ml/gcc-patches/2015-12/msg02110.html https://gcc.gnu.org/ml/fortran/2016-08/msg00029.html — Vladimir F Героям слава, Aug 11 '16 at 18:55

score 0 · Answer 1 · answered Aug 11 '16 at 18:46

Your call

call random_seed(seed)

is strange. I thought it should not be allowed without a keyword argument, but it actually is inquiring for the size of the random seed array.

For a proper way of initializing seed see the example in

https://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html

Issues with setting random seed

1 Answers1