simplified halto procedure (still needs testing)
This commit is contained in:
Martin Diehl
parent
98df2d1427
commit
ae27660e86
252
src/math.f90
252
src/math.f90
|
|
@ -161,10 +161,8 @@ module math
|
|||
math_rotate_forward3333, &
|
||||
math_limit
|
||||
private :: &
|
||||
halton, &
|
||||
halton_memory, &
|
||||
halton_ndim_set, &
|
||||
halton_seed_set
|
||||
math_check, &
|
||||
halton
|
||||
|
||||
contains
|
||||
|
||||
|
|
@ -217,8 +215,7 @@ subroutine math_init
|
|||
|
||||
call random_seed(put = randInit)
|
||||
|
||||
call halton_seed_set(int(randInit(1), pInt))
|
||||
call halton_ndim_set(3_pInt)
|
||||
randTest = halton(4,seed = randTest(1))
|
||||
|
||||
call math_check()
|
||||
|
||||
|
|
@ -283,7 +280,7 @@ subroutine math_check
|
|||
endif
|
||||
|
||||
! +++ check rotation sense of q and R +++
|
||||
call halton(3_pInt,v) ! random vector
|
||||
v = halton(3_pInt) ! random vector
|
||||
R = math_qToR(q)
|
||||
if (any(abs(math_mul33x3(R,v) - math_qRot(q,v)) > tol_math_check)) then
|
||||
write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v'
|
||||
|
|
@ -1237,7 +1234,7 @@ function math_qRand()
|
|||
real(pReal), dimension(4) :: math_qRand
|
||||
real(pReal), dimension(3) :: rnd
|
||||
|
||||
call halton(3_pInt,rnd)
|
||||
rnd = halton(3_pInt)
|
||||
math_qRand = [cos(2.0_pReal*PI*rnd(1))*sqrt(rnd(3)), &
|
||||
sin(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), &
|
||||
cos(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), &
|
||||
|
|
@ -1760,7 +1757,7 @@ function math_sampleRandomOri()
|
|||
implicit none
|
||||
real(pReal), dimension(3) :: math_sampleRandomOri, rnd
|
||||
|
||||
call halton(3_pInt,rnd)
|
||||
rnd = halton(3_pInt)
|
||||
math_sampleRandomOri = [rnd(1)*2.0_pReal*PI, &
|
||||
acos(2.0_pReal*rnd(2)-1.0_pReal), &
|
||||
rnd(3)*2.0_pReal*PI]
|
||||
|
|
@ -1793,7 +1790,7 @@ function math_sampleGaussOri(center,FWHM)
|
|||
cosScatter = cos(scatter)
|
||||
|
||||
do
|
||||
call halton(5_pInt,rnd)
|
||||
rnd = halton(5_pInt)
|
||||
rnd(1:3) = 2.0_pReal*rnd(1:3)-1.0_pReal ! expand 1:3 to range [-1,+1]
|
||||
disturb = [ scatter * rnd(1), & ! phi1
|
||||
sign(1.0_pReal,rnd(2))*acos(cosScatter+(1.0_pReal-cosScatter)*rnd(4)), & ! Phi
|
||||
|
|
@ -1853,7 +1850,7 @@ function math_sampleFiberOri(alpha,beta,FWHM)
|
|||
|
||||
! ---# rotation matrix about fiber axis (random angle) #---
|
||||
do
|
||||
call halton(6_pInt,rnd)
|
||||
rnd = halton(6_pInt)
|
||||
fRot = math_EulerAxisAngleToR(fiberInS,rnd(1)*2.0_pReal*pi)
|
||||
|
||||
! ---# rotation about random axis perpend to fiber #---
|
||||
|
|
@ -1909,7 +1906,7 @@ real(pReal) function math_sampleGaussVar(meanvalue, stddev, width)
|
|||
myWidth = merge(width,3.0_pReal,present(width)) ! use +-3*sigma as default value for scatter if not given
|
||||
|
||||
do
|
||||
call halton(2_pInt, rnd)
|
||||
rnd = halton(2_pInt)
|
||||
scatter = myWidth * (2.0_pReal * rnd(1) - 1.0_pReal)
|
||||
if (rnd(2) <= exp(-0.5_pReal * scatter ** 2.0_pReal)) exit ! test if scattered value is drawn
|
||||
enddo
|
||||
|
|
@ -2286,32 +2283,6 @@ pure function math_invariantsSym33(m)
|
|||
end function math_invariantsSym33
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief computes the next element in the Halton sequence.
|
||||
!> @author John Burkardt
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine halton(ndim, r)
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: ndim !< dimension of the element
|
||||
real(pReal), intent(out), dimension(ndim) :: r !< next element of the current Halton sequence
|
||||
integer(pInt), dimension(ndim) :: base
|
||||
integer(pInt) :: seed
|
||||
integer(pInt), dimension(1) :: value_halton
|
||||
|
||||
call halton_memory ('GET', 'SEED', 1_pInt, value_halton)
|
||||
seed = value_halton(1)
|
||||
|
||||
call halton_memory ('GET', 'BASE', ndim, base)
|
||||
|
||||
call i_to_halton (seed, base, ndim, r)
|
||||
|
||||
value_halton(1) = 1_pInt
|
||||
call halton_memory ('INC', 'SEED', 1_pInt, value_halton)
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
contains
|
||||
|
||||
!-------------------------------------------------------------------------------------------------
|
||||
!> @brief computes an element of a Halton sequence.
|
||||
!> @details Only the absolute value of SEED is considered. SEED = 0 is allowed, and returns R = 0.
|
||||
|
|
@ -2322,136 +2293,27 @@ subroutine halton(ndim, r)
|
|||
!> @details multi-dimensional integrals, Numerische Mathematik, Volume 2, pages 84-90, 1960.
|
||||
!> @author John Burkardt
|
||||
!-------------------------------------------------------------------------------------------------
|
||||
subroutine i_to_halton (seed, base, ndim, r)
|
||||
use IO, only: &
|
||||
IO_error
|
||||
function halton(ndim, seed)
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: &
|
||||
ndim, & !< dimension of the sequence
|
||||
seed !< index of the desired element
|
||||
integer(pInt), intent(in), dimension(ndim) :: base !< Halton bases
|
||||
real(pReal), intent(out), dimension(ndim) :: r !< the SEED-th element of the Halton sequence for the given bases
|
||||
|
||||
real(pReal), dimension(ndim) :: base_inv
|
||||
ndim !< dimension of the sequence
|
||||
real(pReal), intent(in), optional :: &
|
||||
seed !< seed value [0.0,1.0]
|
||||
real(pReal), dimension(ndim) :: &
|
||||
halton
|
||||
integer(pInt), save :: &
|
||||
start, &
|
||||
current
|
||||
real(pReal), dimension(ndim) :: &
|
||||
base_inv
|
||||
integer(pInt), dimension(ndim) :: &
|
||||
digit, &
|
||||
seed2
|
||||
|
||||
seed2 = abs(seed)
|
||||
r = 0.0_pReal
|
||||
|
||||
if (any (base(1:ndim) <= 1_pInt)) call IO_error(error_ID=405_pInt)
|
||||
|
||||
base_inv(1:ndim) = 1.0_pReal / real (base(1:ndim), pReal)
|
||||
|
||||
do while ( any ( seed2(1:ndim) /= 0_pInt) )
|
||||
digit(1:ndim) = mod ( seed2(1:ndim), base(1:ndim))
|
||||
r(1:ndim) = r(1:ndim) + real ( digit(1:ndim), pReal) * base_inv(1:ndim)
|
||||
base_inv(1:ndim) = base_inv(1:ndim) / real ( base(1:ndim), pReal)
|
||||
seed2(1:ndim) = seed2(1:ndim) / base(1:ndim)
|
||||
enddo
|
||||
|
||||
end subroutine i_to_halton
|
||||
|
||||
|
||||
end subroutine halton
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief sets or returns quantities associated with the Halton sequence.
|
||||
!> @details If action_halton is 'SET' and action_halton is 'BASE', then NDIM is input, and
|
||||
!> @details is the number of entries in value_halton to be put into BASE.
|
||||
!> @details If action_halton is 'SET', then on input, value_halton contains values to be assigned
|
||||
!> @details to the internal variable.
|
||||
!> @details If action_halton is 'GET', then on output, value_halton contains the values of
|
||||
!> @details the specified internal variable.
|
||||
!> @details If action_halton is 'INC', then on input, value_halton contains the increment to
|
||||
!> @details be added to the specified internal variable.
|
||||
!> @author John Burkardt
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine halton_memory (action_halton, name_halton, ndim, value_halton)
|
||||
use IO, only: &
|
||||
IO_lc
|
||||
|
||||
implicit none
|
||||
character(len = *), intent(in) :: &
|
||||
action_halton, & !< desired action: GET the value of a particular quantity, SET the value of a particular quantity, INC the value of a particular quantity (only for SEED)
|
||||
name_halton !< name of the quantity: BASE: Halton base(s), NDIM: spatial dimension, SEED: current Halton seed
|
||||
integer(pInt), dimension(*), intent(inout) :: value_halton
|
||||
integer(pInt), allocatable, save, dimension(:) :: base
|
||||
logical, save :: first_call = .true.
|
||||
integer(pInt), intent(in) :: ndim !< dimension of the quantity
|
||||
integer(pInt), save :: ndim_save = 0_pInt, seed = 1_pInt
|
||||
integer(pInt) :: i
|
||||
|
||||
if (first_call) then
|
||||
ndim_save = 1_pInt
|
||||
allocate(base(ndim_save))
|
||||
base(1) = 2_pInt
|
||||
first_call = .false.
|
||||
endif
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! Set
|
||||
actionHalton: if(IO_lc(action_halton(1:1)) == 's') then
|
||||
|
||||
nameSet: if(IO_lc(name_halton(1:1)) == 'b') then
|
||||
if(ndim_save /= ndim) ndim_save = ndim
|
||||
base = value_halton(1:ndim)
|
||||
elseif(IO_lc(name_halton(1:1)) == 'n') then nameSet
|
||||
if(ndim_save /= value_halton(1)) then
|
||||
ndim_save = value_halton(1)
|
||||
base = [(prime(i),i=1_pInt,ndim_save)]
|
||||
else
|
||||
ndim_save = value_halton(1)
|
||||
endif
|
||||
elseif(IO_lc(name_halton(1:1)) == 's') then nameSet
|
||||
seed = value_halton(1)
|
||||
endif nameSet
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! Get
|
||||
elseif(IO_lc(action_halton(1:1)) == 'g') then actionHalton
|
||||
nameGet: if(IO_lc(name_halton(1:1)) == 'b') then
|
||||
if(ndim /= ndim_save) then
|
||||
ndim_save = ndim
|
||||
base = [(prime(i),i=1_pInt,ndim_save)]
|
||||
endif
|
||||
value_halton(1:ndim_save) = base(1:ndim_save)
|
||||
elseif(IO_lc(name_halton(1:1)) == 'n') then nameGet
|
||||
value_halton(1) = ndim_save
|
||||
elseif(IO_lc(name_halton(1:1)) == 's') then nameGet
|
||||
value_halton(1) = seed
|
||||
endif nameGet
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! Increment
|
||||
elseif(IO_lc(action_halton(1:1)) == 'i') then actionHalton
|
||||
if(IO_lc(name_halton(1:1)) == 's') seed = seed + value_halton(1)
|
||||
endif actionHalton
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
contains
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief returns any of the first 1500 prime numbers.
|
||||
!> @details n = 0 is legal, returning PRIME = 1.
|
||||
!> @details Reference:
|
||||
!> @details Milton Abramowitz and Irene Stegun: Handbook of Mathematical Functions,
|
||||
!> @details US Department of Commerce, 1964, pages 870-873.
|
||||
!> @details Daniel Zwillinger: CRC Standard Mathematical Tables and Formulae,
|
||||
!> @details 30th Edition, CRC Press, 1996, pages 95-98.
|
||||
!> @author John Burkardt
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
integer(pInt) function prime(n)
|
||||
use IO, only: &
|
||||
IO_error
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: n !< index of the desired prime number
|
||||
base, &
|
||||
t
|
||||
integer(pInt) :: &
|
||||
i
|
||||
integer(pInt), dimension(0:1600), parameter :: &
|
||||
npvec = int([&
|
||||
prime = int([&
|
||||
1, &
|
||||
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, &
|
||||
31, 37, 41, 43, 47, 53, 59, 61, 67, 71, &
|
||||
|
|
@ -2629,56 +2491,38 @@ subroutine halton_memory (action_halton, name_halton, ndim, value_halton)
|
|||
13313, 13327, 13331, 13337, 13339, 13367, 13381, 13397, 13399, 13411, &
|
||||
13417, 13421, 13441, 13451, 13457, 13463, 13469, 13477, 13487, 13499],pInt)
|
||||
|
||||
if (n < size(npvec)) then
|
||||
prime = npvec(n)
|
||||
|
||||
if (present(seed)) then
|
||||
start = int(seed * 1599.0_pReal,pInt) + 1_pInt ! select one prime number to start with
|
||||
current = 1_pInt
|
||||
else
|
||||
call IO_error(error_ID=406_pInt)
|
||||
current = current + 1_pInt
|
||||
endif
|
||||
|
||||
end function prime
|
||||
base = [(prime(mod(start,1600_pInt)+i), i=1_pInt, ndim)]
|
||||
base_inv = 1.0_pReal/real(base,pReal)
|
||||
|
||||
end subroutine halton_memory
|
||||
halton = 0.0_pReal
|
||||
t = current
|
||||
|
||||
do while (any( t /= 0_pInt) )
|
||||
halton = halton + real(mod(t,base), pReal) * base_inv
|
||||
base_inv = base_inv / real(base, pReal)
|
||||
t = t / base
|
||||
enddo
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief sets the dimension for a Halton sequence
|
||||
!> @brief returns any of the first 1600 prime numbers.
|
||||
!> @details n = 0 is legal, returning PRIME = 1.
|
||||
!> @details 0 > n > 1600 returns -1
|
||||
!> @details Reference:
|
||||
!> @details Milton Abramowitz and Irene Stegun: Handbook of Mathematical Functions,
|
||||
!> @details US Department of Commerce, 1964, pages 870-873.
|
||||
!> @details Daniel Zwillinger: CRC Standard Mathematical Tables and Formulae,
|
||||
!> @details 30th Edition, CRC Press, 1996, pages 95-98.
|
||||
!> @author John Burkardt
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine halton_ndim_set(ndim)
|
||||
|
||||
implicit none
|
||||
integer(pInt), intent(in) :: ndim !< dimension of the Halton vectors
|
||||
integer(pInt) :: value_halton(1)
|
||||
|
||||
value_halton(1) = ndim
|
||||
call halton_memory ('SET', 'NDIM', 1_pInt, value_halton)
|
||||
|
||||
end subroutine halton_ndim_set
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief sets the seed for the Halton sequence.
|
||||
!> @details Calling HALTON repeatedly returns the elements of the Halton sequence in order,
|
||||
!> @details starting with element number 1.
|
||||
!> @details An internal counter, called SEED, keeps track of the next element to return. Each time
|
||||
!> @details is computed, and then SEED is incremented by 1.
|
||||
!> @details To restart the Halton sequence, it is only necessary to reset SEED to 1. It might also
|
||||
!> @details be desirable to reset SEED to some other value. This routine allows the user to specify
|
||||
!> @details any value of SEED.
|
||||
!> @details The default value of SEED is 1, which restarts the Halton sequence.
|
||||
!> @author John Burkardt
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine halton_seed_set(seed)
|
||||
implicit none
|
||||
|
||||
integer(pInt), parameter :: NDIM = 1_pInt
|
||||
integer(pInt), intent(in) :: seed !< seed for the Halton sequence.
|
||||
integer(pInt) :: value_halton(ndim)
|
||||
|
||||
value_halton(1) = seed
|
||||
call halton_memory ('SET', 'SEED', NDIM, value_halton)
|
||||
|
||||
end subroutine halton_seed_set
|
||||
end function halton
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
|
|||
Loading…
Reference in New Issue