Gippo
Gippo

Reputation: 65

Random permutation and random seed

I am employing the Knuth algorithm to generate a random permutation of an n-tuple. This is the code. Fixed n, it generates random permutations and collect all the different ones until it finds all the n! permutations. At the end it prints also the number of trials needed to find all the permutations. I have also inserted the initialization of the seed from the time (in a very simple and naive way, though). There are two options (A and B). A: The seed is fixed once for all in the main program. B: The seed is fixed every time a random permutation is computed (below the second option is commented).

implicit none
 

  integer :: n,ncomb
  integer :: i,h,k,x
  integer, allocatable :: list(:),collect(:,:)
  logical :: found
  integer :: trials
  !
  ! A
  !
  integer :: z,values(1:8)
  integer, dimension(:), allocatable :: seed
  call date_and_time(values=values)
  call random_seed(size=z)
  allocate(seed(1:z))
  seed(:) = values(8)
  call random_seed(put=seed)
 

  n=4
  
  
  ncomb=product((/(i,i=1,n)/))
  allocate(list(n))
  allocate(collect(n,ncomb))  
  
 trials=0 
 h=0
 do 
  trials=trials+1
  list=Shuffle(n)


  found=.false.
  do k=1,h
   x=sum(abs(list-collect(:,k))) 
   if ( x == 0 ) then
     found=.true.
     exit
   end if
  end do

  if ( .not. found ) then
   h=h+1
   collect(:,h)=list
   print*,h,')',collect(:,h)
  end if

  if ( h == ncomb ) exit
 end do 
  
 write(*,*) "Trials= ",trials
  
  
contains
 
function Shuffle(n) result(list)
  integer, allocatable :: list(:)
  integer, intent(in) :: n
  integer :: i, randpos, temp,h
  real :: r
  !
  ! B
  !
!   integer :: z,values(1:8)
!   integer, dimension(:), allocatable :: seed
!   call date_and_time(values=values)
!   call random_seed(size=z)
!   allocate(seed(1:z))
!   seed(:) = values(8)
!   call random_seed(put=seed)

 
 
  allocate(list(n))
  list = (/ (h, h=1,n) /)
  do i = n, 2, -1
    call random_number(r)
    randpos = int(r * i) + 1
    temp = list(randpos)
    list(randpos) = list(i)
    list(i) = temp
  end do
 
end function Shuffle
 
end 

You can check that the second option is not good at all. For n=4 it takes around 100 times more trials to obtain the total number of permutations and for n=5 it gets stuck.

My questions are:

  1. Why does calling random_seed multiple times give wrong results? What kind of systematic error I am introducing? Isn't it equivalent to calling random seed only once but launching the code several times (each time generating only one random permutation)?

  2. If I want to launch several times the code, computing a single permutation, I guess that if I initialize the random seed I have the same problem (regardless the position of the initialization, since now I am computing only one permutation). Correct? In this case, what I have to do in order to initialize the seed wihouth spoiling the uniform sampling? Because if I do not initialize the seed in a random way I obtain the same permutation. I guess I could print and read the seed everytime I launch the code, in order to not to start from the same pseudo-random numbers. However, this is complicated to do if I launch several instances of the code in parallel.

UPDATE

I have understood the reply. In conclusion, if I want to generate pseudorandom numbers at each call by initializing the seed, what I can do is:

A) Old gfortran

Use the subroutine init_random_seed() here

https://gcc.gnu.org/onlinedocs/gcc-5.1.0/gfortran/RANDOM_005fSEED.html

B) Most recent gfortran versions

call random_seed()

C) Fortran2018

call random_init(repeatable, image_distinct)

Questions

In the C) case, should I set repeatable=.false., image_distinct=.true. to have a different random number each time?

What could be an efficient way to write the code in a portable way, so that it works whatever the compiler is? (I mean, the code recognizes what is available and works accordingly)

Upvotes: 1

Views: 382

Answers (1)

You certainly should not ever call random_seed() repeatedly. It is supposed to be called just once. You make the problem worse by setting it in such a crude way.

Yes, one does often use the data and time to initialize it, but one must pre-process the time data through something that adds some entropy, like through some very simple random-generator. A good example can be found in the documentation of RANDOM_SEED for older versions of gfortran: https://gcc.gnu.org/onlinedocs/gcc-5.1.0/gfortran/RANDOM_005fSEED.html See how lcg() is used there to transform the data_and_time() data.

Note that more recent versions of gfortran will generate a random seed that is different every time just by calling random_seed() without any arguments. Older versions returned the same seed every time.

Also note that Fortran 2018 has random_init() where you can specify repeatable= to be true or false. With false you get a different sequence every time.

The portable thing is to use standard Fortran, that's all. But you cannot use new features and old compiler versions at the same time. This kind of portability does not exist. With old compilers you can only use old standard features. I won't even start writing about autoconf and stuff, it is not worth it.


So,

  1. you can set your random number seed to be the same every time or distinct every time (see above),

and

  1. you should always call random_seed or random_init only once.

Why does calling random_seed multiple times give wrong results?

You are re-starting the pseudorandom sequence to some unspecified state, probably with insufficient entropy. Quite easily quite close to the last starting state.

What kind of systematic error I am introducing? Isn't it equivalent to calling random seed only once but launching the code several times (each time generating only one random permutation)?

It might be similar. But your seeding using the time is way too naïve and when running in the loop the date and time is way too similar if not completely equal in most bits. Some transform as linked above might mask that problem anyway but putting the date and time itself as your seed is just not going to work.

Upvotes: 4

Related Questions