random - MT19937 不会通过将种子值保持为常数来重现相同的伪随机序列

Question

我正在 Fortran 90/95 的 Monte Carlo 模拟中编写一个检查点函数，我使用的编译器是 ifort 18.0.2，在详细说明我正在使用的伪随机生成器的版本之前:

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

见 mt19937对于源代码。

我的蒙特卡罗模拟代码的一般结构如下:

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

内 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

根据我的理解，如果种子值保持不变，PRNG 应该每次都能重现伪随机序列？

为了证明情况确实如此，我使用相同的种子值运行了两个单独的模拟，它们能够重现准确的序列。到现在为止还挺好!

基于之前的测试，我会进一步假设无论次数如何 init_genrand()被叫 内一个单独的模拟，PRNG 也应该能够重现伪随机值序列？所以我对我的 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 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

输出出乎意料地不是我想的那样，无论如何 iseed两次初始化之间的输出是相同的，但是， genrand_real3()输出不相同。

由于这个意外的结果，我努力在系统的任意状态下恢复模拟，因为模拟没有重现我正在模拟的系统的最新配置状态。

我不确定我是否提供了足够的信息，请让我知道这个问题的任何部分是否需要更具体？

Answer 1

从您提供的源代码(请参阅 [mt19937]{ http://web.mst.edu/~vojtat/class_5403/mt19937/mt19937ar.f90 } 获取源代码。)，init_genrand 不清楚 整个州。

有 3 个临界状态变量:

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

第一个是“状态向量的数组”，第二个是确保我们不会以未初始化的数组开始的标志，第三个是一些位置标记，正如我从评论中陈述的条件猜测的那样。

看着 subroutine init_genrand( s ) ，它设置 mtinit标志，并填充 mt()来自 1 的数组高达 n .好吧。

看着 genrand_real3它基于 genrand_int32 .

看着 genrand_int32 ，它开始于

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 )

并执行其算术魔术，然后开始获得结果:

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

所以.. mti是“状态数组”中的位置索引，每次从生成器读取整数后，它都会增加 1。

返回 init_genrand - 记住？它一直在重置阵列 mt()但它 尚未重置 MTI 回到起点 mti = n + 1_wi .

我敢打赌，这就是您观察到的现象的原因，因为在使用相同的种子重新初始化后，数组将填充相同的一组值，但稍后 int32 生成器将从不同的起点读取。我怀疑它是故意的，所以它可能是一个容易被忽视的小错误。