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

Я пишу функцию контрольной точки в своей симуляции Монте-Карло на Фортране 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() выходы не идентичны.

Из-за этого неожиданного результата я изо всех сил пытался возобновить симуляцию в произвольном состоянии системы, поскольку симуляция не воспроизводит последнее состояние конфигурации моделируемой системы.

Я не уверен, что предоставил достаточно информации. Пожалуйста, дайте мне знать, если какая-либо часть этого вопроса должна быть более конкретной?


random fortran mersenne-twister
person Gvxfjørt    schedule 13.06.2018    source источник
comment
Пожалуйста, используйте тег [tagnfortran] для всех вопросов по Fortran. Нет смысла маркировать fortran90 и fortran95, ваш вопрос в любом случае не зависит от версии.   -  person Vladimir F    schedule 13.06.2018
comment
Как init_genrand заполняет полное состояние PRNG?   -  person francescalus    schedule 13.06.2018
comment
@francescalus, пожалуйста, см. [mt19937ar]{web.mst.edu/~vojtat/ class_5403/mt19937/mt19937ar.f90}, мой первый тест показал, что он может воспроизвести последовательность, задав ей постоянное начальное значение, но во втором тесте, когда я повторно инициализировал последовательность с тем же начальным значением , он каким-то образом не сбросил состояние последовательности до того места, где оно началось. функция init_genrand, когда я читаю подпрограмму, мне кажется, что init_genrand() инициализируется там, где последовательность начинается с заданного начального значения, а genrand_real3() предоставляет значение в (0,1). Есть шанс, что вы видите какие-либо разногласия?   -  person Gvxfjørt    schedule 13.06.2018

Ответы (1)


arrow_upward
1
arrow_downward

Из предоставленного вами исходного кода (см. [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 будет читать из другой начальной точки. Я сомневаюсь, что это было задумано, так что, вероятно, это крошечная ошибка, которую легко не заметить.

person quetzalcoatl    schedule 13.06.2018
comment
видимо, я не единственный, кто столкнулся с этой ошибкой, это известная проблема реализации перезапуска контрольной точки на основе mt19937, сброс позиционного индекса требует знания точного общего количества вызовов genrand_real() и genrand_int(), и нет не является общим решением для этого, так как это зависит от приложения. Но спасибо за указание на проблему! - person Gvxfjørt; 27.06.2018