MT19937 НЕ воспроизводит одну и ту же псевдослучайную последовательность, сохраняя начальное значение постоянным

Question

Я пишу функцию контрольной точки в моем симуляции Монте-Карло в Fortran 90/95, компилятор, который я использую, это 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 будет считывать данные из другой начальной точки. Я сомневаюсь, что это было задумано, так что это, вероятно, крошечная ошибка, которую легко не заметить.