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

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

0 голосов
/ 13 июня 2018

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

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

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

1 Ответ

0 голосов
/ 13 июня 2018

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

...