Modern fortran implementation of a Sobol sequence

This is a fork of the module sobseq by Dann van Vugt. I've added more intuitive functionality and packaging.

The Sobol sequence is a low-discrepancy sequence that can be used for more efficient (Quasi-)Monte-Carlo integration.

This repository contains a module, sobseq.f90, which can be used to generate points in Sobol series, in a continuous or strided (2**n) manner.

Setup using the Fortran Package Manager (fpm)

To use sobseq for testing in your fpm project without making it a build dependency, just add the following to your package manifest file (fpm.toml):

[dependencies]
sobseq = { git = "https://github.com/renatomatz/sobseq.git" }

You can then use the package in your testing programs with use sobseq.

Examples

To generate sobol sequences in many dimensions, direction numbers are needed. A good set of these was made by Joe and Kuo. The direction numbers for the first eight are reproduced here, and more can be found at the link above.

integer, parameter, dimension(1:12)   :: s = (/1,2,3,3,4,4,5,5,5,5,5,5/)
integer, parameter, dimension(1:12)   :: a = (/0,1,1,2,1,4,2,4,7,11,13,14/)
integer, parameter, dimension(12,5) :: m = &
    transpose(reshape( &
        [1, 0, 0, 0,  0,  &
         1, 3, 0, 0,  0,  &
         1, 3, 1, 0,  0,  &
         1, 1, 1, 0,  0,  &
         1, 1, 3, 3,  0,  &
         1, 3, 5, 13, 0,  &
         1, 1, 5, 5,  17, &
         1, 1, 5, 5,  5,  &
         1, 1, 7, 11, 19, &
         1, 1, 5, 1,  1,  &
         1, 1, 1, 3,  11, &
         1, 3, 5, 5,  31], [12,5])

The following code snippet contains a small example.

integer, parameter :: n_dim=9
integer, parameter :: n_samples=200
type(sobol_state), dimension(n_dim) :: rng
integer :: i, j
real :: u(n_dim)

! Initialization
do i=1,n_dim
  call rng%initialize(s(i), a(i), m(:,i))
end do

! Generation
do j=1,n_samples
  do i=1,n_dim
    u(i) = rng(i)%next()
  end do

  ! do something with u
end do

Strided operation

The module also contains functions for generating numbers in a strided manner, which is useful for parallel operation. To use this, set n_streams to the 2-logarithm of the number of streams you want (if n_streams is not a power of 2 you might get distribution problems) and use the rng%next_strided() generator.

Original Authors

Daan van Vugt daanvanvugt@gmail.com
Koen Beljaars k.p.beljaars@tue.nl

License

This work is released under the MIT license. A copy can be found in the file LICENSE in the repository.