0

I'm running a Fortran code which performs a stochastic simulation of a marked Poisson cluster process. In practice, event properties (eg. time of occurrences) are generated by inversion method, i.e. by random sampling of the cumulative distribution function. Because of the Poissonian randomness, I expect each generated sequence to be different, but this is not the case. I guess the reason is that the seed for the pseudorandom number generator is the same at each simulation. I do not know Fortran, so I have no idea how to solve this issue. Here is the part of the code involved with the pseudorandom number generator, any idea?

subroutine pseud0(r)
c     generation of pseudo-random numbers
c     data ir/584287/
      data ir/574289/
      ir=ir*48828125
      if(ir) 10,20,20
   10 ir=(ir+2147483647)+1
   20 r=float(ir)*0.4656613e-9
      return
      end
      subroutine pseudo(random)
c     wichmann+hill (1982) Appl. Statist 31
      data ix,iy,iz /1992,1111,1151/
      ix=171*mod(ix,177)-2*(ix/177)
      iy=172*mod(iy,176)-35*(iy/176)
      iz=170*mod(iz,178)-63*(iz/178)
      if (ix.lt.0) ix=ix+30269
      if (iy.lt.0) iy=iy+30307
      if (iz.lt.0) iz=iz+30323
      random=mod(float(ix)/30269.0+float(iy)/30307.0+
     &            float(iz)/30323.0,1.0)
      return
      end
Roland
  • 427
  • 1
  • 4
  • 15
  • 1
    The seed is hardcoded in to the source. You can change the values there simply (in the `data` statements), but if you want to do something more exotic then this will require quite a redesign. – francescalus Oct 09 '19 at 15:37
  • The problem is I need to run this code at least 100 times, which would make it a hassle to manually edit the seed each time. – Roland Oct 09 '19 at 15:54
  • One method is to pass the seed from the command-line, so that your code does not require changes and recompilation. Another method is to read "randomness" from `/dev/urandom` (posix systems). – Pierre de Buyl Oct 09 '19 at 21:00
  • Thank you @PierredeBuyl, shouldn't I edit the code for it to accept input from the command line? – Roland Oct 11 '19 at 14:48

1 Answers1

2

First, I would review the modern literature for PRNG and pick a modern implementation. Second, I would rewrite the code in modern Fortran.

You need to follow @francescalus advice and have a method for updating the seed. Without attempting to modernizing your code, here is one method for the pseud0 prng

 subroutine init0(i)
    integer, intent(in) :: i
    common /myseed0/iseed
    iseed = i
 end subroutine init0

subroutine pseud0(r)
   common /myseed0/ir
   ir = ir * 48828125
   if (ir) 10,20,20
10 ir = (ir+2147483647)+1
20 r = ir*0.4656613e-9
end subroutine pseud0

program foo
   integer i
   real r1
   call init0(574289)  ! Original seed
   do i = 1, 10
      call pseud0(r1)
      print *, r1
   end do
   print *
   call init0(289574)  ! New seed
   do i = 1, 10
      call pseud0(r1)
      print *, r1
   end do
   print *
end program foo
steve
  • 657
  • 1
  • 4
  • 7
  • Thank you @steve. Indeed, this is a very old code, that is why it's even more difficult for those who do not know Fortran. Do you have any idea why there is the `pseud0` prng and the `pseudo(random)` prng? What is the difference? As far as I can see, in the whole code, only the `pseud0` prng is called. – Roland Oct 09 '19 at 19:25
  • @francescalus, of course, "modern" for some means to avoid `common`. I wasn't attempting to be modern, only showing the OP one method to meet OP's need. Hopefully, OP follows my first recommendation, and looks for a modern PRNG. – steve Oct 09 '19 at 19:44
  • @Angela You could do worse than Wichmann-Hill for the generator, but you might look into the FORTRAN implementation of [xoroshiro128+](https://github.com/jannisteunissen/rng_fortran). It does well on standard test suites and is very fast. – pjs Oct 09 '19 at 20:17
  • @francescalus, I updated my answer to note that I have made no attempt at modernizing the code. Just reading the old code makes my eyes hurt. – steve Oct 09 '19 at 20:52
  • @steve , the problem is that I do not see why the Wichmann-Hill generator is there, the pseudo(random) subroutine is never called by the program. – Roland Oct 11 '19 at 14:52
  • If the routine is never used, a good compiler will simple reap the code. My guess is that whoever wrote the original code would switch between the two PRNG routine for testing purposes. – steve Oct 11 '19 at 17:02