Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

MT19937 does NOT reproduce the same pseudo-random sequence by holding the seed value a constant

I'm writing a checkpoint function in my Monte Carlo simulation in Fortran 90/95, the compiler I'm using is ifort 18.0.2, before going through detail just to clarify the version of pseudo-random generator I'm using:

A C-program for MT19937, with initialization, improved 2002/1/26.
Coded by Takuji Nishimura and Makoto Matsumoto.

Code converted to Fortran 95 by Josi Rui Faustino de Sousa
Date: 2002-02-01

See mt19937 for the source code.

The general structure of my Monte Carlo simulation code is given below:

program montecarlo
 call read_iseed(...)
 call mc_subroutine(...)
end 

Within the read_iseed

subroutine read_iseed(...)
  use mt19937

    if (Restart == 'n') then

    call system('od -vAn -N4 -td4 < /dev/urandom > '//trim(IN_ISEED)
    open(unit=7,file=trim(IN_ISEED),status='old')
    read(7,*) i
    close(7)

    !This is only used to initialise the PRNG sequence
    iseed = abs(i)
    else if (Restart == 'y') then

    !Taking seed value from the latest iteration of previous simulation
    iseed = RestartSeed

    endif

    call init_genrand(iseed)
    print *, 'first pseudo-random value ',genrand_real3(), 'iseed ',iseed

    return
end subroutine

Based on my understanding, if the seed value holds a constant, the PRNG should be able to reproduce the pseudo-random sequence every time?

In order to prove this is the case, I ran two individual simulations by using the same seed value, they are able to reproduce the exact sequence. So far so good!

Based on the previous test, I'd further assume that regardless the number of times init_genrand() being called within one individual simulation, the PRNG should also be able to reproduce the pseudo-random value sequence? So I did a little modification to my read_iseed() subroutine

subroutine read_iseed(...)
  use mt19937

    if (Restart == 'n') then

    call system('od -vAn -N4 -td4 < /dev/urandom > '//trim(IN_ISEED)
    open(unit=7,file=trim(IN_ISEED),status='old')
    read(7,*) i
    close(7)

    !This is only used to initialise the PRNG sequence
    iseed = abs(i)
    else if (Restart == 'y') then

    !Taking seed value from the latest iteration of the previous simulation
    iseed = RestartSeed

    endif

    call init_genrand(iseed)
    print *, 'first time initialisation ',genrand_real3(), 'iseed ',iseed

    call init_genrand(iseed)
    print *, 'second time initialisation ',genrand_real3(), 'iseed ',iseed

    return
end subroutine

The output is surprisingly not the case I thought would be, by all means iseed outputs are identical in between two initializations, however, genrand_real3() outputs are not identical.

Because of this unexpected result, I struggled with resuming the simulation at an arbitrary state of the system since the simulation is not reproducing the latest configuration state of the system I'm simulating.

I'm not sure if I've provided enough information, please let me know if any part of this question needs to be more specific?

like image 381
Gvxfjørt Avatar asked Nov 08 '22 05:11

Gvxfjørt


1 Answers

From the source code you've provided (See [mt19937]{http://web.mst.edu/~vojtat/class_5403/mt19937/mt19937ar.f90} for the source code.), the init_genrand does not clear the whole state.

There are 3 critical state variables:

integer( kind = wi )  :: mt(n)            ! the array for the state vector
logical( kind = wi )  :: mtinit = .false._wi   ! means mt[N] is not initialized
integer( kind = wi )  :: mti = n + 1_wi   ! mti==N+1 means mt[N] is not initialized

The first one is the "array for the state vector", second one is a flag that ensures we don't start with uninitialized array, and the third one is some position marker, as I guess from the condition stated in the comment.

Looking at subroutine init_genrand( s ), it sets mtinit flag, and fills the mt() array from 1 upto n. Alright.

Looking at genrand_real3 it's based on genrand_int32.

Looking at genrand_int32, it starts up with

if ( mti > n ) then ! generate N words at one time
  ! if init_genrand() has not been called, a default initial seed is used
  if ( .not. mtinit ) call init_genrand( seed_d )

and does its arithmetic magic and then starts getting the result:

y = mt(mti)
mti = mti + 1_wi

so.. mti is a positional index in the 'state array', and it is incremented by 1 after each integer read from the generator.

Back to init_genrand - remember? it have been resetting the array mt() but it has not resetted the MTI back to its starting mti = n + 1_wi.

I bet this is the cause of the phenomenon you've observed, since after re-initializing with the same seed, the array would be filled with the same set of values, but later the int32 generator would read from a different starting point. I doubt it was intended, so it's probably a tiny bug easy to overlook.

like image 57
quetzalcoatl Avatar answered Jan 01 '23 21:01

quetzalcoatl