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()从1upto填充数组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

so..mti是“状态数组”中的位置索引，在从生成器读取每个整数后，它会加 1。

回到init_genrand——还记得吗？它一直在重置阵列mt()，但尚未将 MTI 重置回其起始位置mti = n + 1_wi。

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