Skip to content
Snippets Groups Projects
Commit e84f14be authored by Yen-Chen Chen's avatar Yen-Chen Chen
Browse files

Move random number generator to math-support (!92)

## What is the new feature
Move random number generator to `libmath-support`. (Adds in https://gitlab.dkrz.de/icon-libraries/libmath-support/-/merge_requests/33

)
## How is it implemented
The random number generator module is more suitable to be put into `libmath-support`. Thus we are moving it from `libfortran-support`.

**The next versions `libfortran-support` 2.0.0 and `libmath-support` 1.1.0 should be merged at the same time.**

Approved-by: default avatarDaniel Hupp <daniel.hupp@meteoswiss.ch>
Merged-by: default avatarYen-Chen Chen <yen-chen.chen@kit.edu>
Changelog: feature
parent a66da626
No related branches found
No related tags found
1 merge request!92Move random number generator to math-support
Pipeline #74119 passed
...@@ -51,7 +51,6 @@ The `libfortran-support` library includes some general Fortran supporting module ...@@ -51,7 +51,6 @@ The `libfortran-support` library includes some general Fortran supporting module
- `mo_io_units`: io unit definitions - `mo_io_units`: io unit definitions
- `mo_namelist`: open/close namelist files - `mo_namelist`: open/close namelist files
- `mo_octree`: octree data structure and operations - `mo_octree`: octree data structure and operations
- `mo_random_number_generators`: generators for pseudo random numbers (should be moved to math-support once available)
- `mo_simple_dump`: array value dumping - `mo_simple_dump`: array value dumping
- `mo_util_backtrace`: function backtrace - `mo_util_backtrace`: function backtrace
- `mo_util_file`: file operations - `mo_util_file`: file operations
......
...@@ -34,7 +34,6 @@ add_library( ...@@ -34,7 +34,6 @@ add_library(
mo_io_units.F90 mo_io_units.F90
mo_namelist.F90 mo_namelist.F90
mo_octree.F90 mo_octree.F90
mo_random_number_generators.F90
mo_simple_dump.F90 mo_simple_dump.F90
mo_util_backtrace.F90 mo_util_backtrace.F90
mo_util_file.F90 mo_util_file.F90
......
! ICON
!
! ---------------------------------------------------------------
! Copyright (C) 2004-2024, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
! Contact information: icon-model.org
!
! See AUTHORS.TXT for a list of authors
! See LICENSES/ for license information
! SPDX-License-Identifier: BSD-3-Clause
! ---------------------------------------------------------------
MODULE mo_random_number_generators
USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64
USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: int32, int64, real64
IMPLICIT NONE
PRIVATE
PUBLIC :: random_normal_values, clip, initialize_random_number_generator, generate_uniform_random_number
INTERFACE generate_uniform_random_number
MODULE PROCEDURE generate_uniform_random_number_wp
END INTERFACE generate_uniform_random_number
INTERFACE random_normal_values
MODULE PROCEDURE random_normal_values_wp
END INTERFACE random_normal_values
INTERFACE clip
MODULE PROCEDURE clip_wp
END INTERFACE clip
! integer parameters for using revised minstd parameters as in
! https://en.wikipedia.org/wiki/Lehmer_random_number_generator#Parameters_in_common_use
INTEGER(KIND=int64), PARAMETER :: minstd_A0 = 16807
INTEGER(KIND=int64), PARAMETER :: minstd_A = 48271
INTEGER(KIND=int64), PARAMETER :: minstd_M = 2147483647
REAL(KIND=real64), PARAMETER :: minstd_Scale = 1.0_real64/REAL(minstd_M, real64)
CONTAINS
! initializes Lehmer random number generator
PURE FUNCTION initialize_random_number_generator(iseed, jseed) RESULT(minstd_state)
INTEGER(KIND=int32), INTENT(IN) :: iseed
INTEGER(KIND=int32), INTENT(IN) :: jseed
INTEGER(KIND=int32) :: minstd_state
! doing kind of an xorshift32 see https://en.wikipedia.org/wiki/Xorshift
minstd_state = iseed
minstd_state = IEOR(minstd_state, jseed*2**13)
minstd_state = IEOR(minstd_state, jseed/2**17)
minstd_state = IEOR(minstd_state, jseed*2**5)
minstd_state = MOD(minstd_A0*minstd_state, minstd_M)
minstd_state = MOD(minstd_A*minstd_state, minstd_M)
END FUNCTION initialize_random_number_generator
! generated Lehmer random number generator https://en.wikipedia.org/wiki/Lehmer_random_number_generator
FUNCTION generate_uniform_random_number_wp(minstd_state) RESULT(random_number)
INTEGER(KIND=int32), INTENT(INOUT) :: minstd_state
REAL(KIND=real64) :: random_number
minstd_state = MOD(minstd_A*minstd_state, minstd_M)
random_number = ABS(minstd_Scale*minstd_state)
END FUNCTION generate_uniform_random_number_wp
SUBROUTINE clip_wp(values_range, values)
! Subroutine arguments (in/out/inout)
REAL(wp), INTENT(IN) :: values_range ! allowed range of random numbers
REAL(wp), DIMENSION(:), INTENT(INOUT) :: values ! array of value to be filled with random numbers
! Limit random number range
WHERE (values(:) > values_range)
values(:) = values_range
ELSEWHERE(values(:) < -values_range)
values(:) = -values_range
END WHERE
END SUBROUTINE clip_wp
SUBROUTINE random_normal_values_wp(seed, values_range, values, mean, stdv)
! Subroutine arguments (in/out/inout)
INTEGER, INTENT(IN) :: seed ! seed to be used
REAL(wp), INTENT(IN) :: values_range ! allowed range of random numbers
REAL(wp), DIMENSION(:), INTENT(INOUT) :: values ! array of value to be filled with random numbers
REAL(wp), INTENT(IN), OPTIONAL :: mean
REAL(wp), INTENT(IN), OPTIONAL :: stdv
REAL(wp), PARAMETER :: pi = 3.141592653589793238_wp ! replace with math constants from math-support once available
! Local variables
INTEGER(KIND=int32) :: i, rng_state, n
REAL(KIND=wp) :: u1, u2, z0, z1, magnitude, theta, mu, sigma
IF (PRESENT(mean)) THEN
mu = mean
ELSE
mu = 0.0_wp
END IF
IF (PRESENT(stdv)) THEN
sigma = stdv
ELSE
sigma = 1.0_wp
END IF
n = SIZE(values)
! -----------------------------------------------------------------------
! Initialize a random number generator, by using the MINSTD linear
! congruential generator (LCG). "nmaxstreams" indicates that random
! numbers will be requested in blocks of this length. The generator
! is seeded with "seed".
!-------------------------------------------------------------------------
DO i = 1, n/2 + 1
rng_state = initialize_random_number_generator(seed, i)
u1 = generate_uniform_random_number(rng_state)
DO WHILE (u1 <= EPSILON(0.0_wp))
u1 = generate_uniform_random_number(rng_state)
END DO
u2 = generate_uniform_random_number(rng_state)
! using the Box-Muller transform https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
magnitude = sigma*SQRT(-2.0_wp*LOG(u1));
theta = 2.0_wp*pi*u2
z0 = magnitude*COS(theta) + mu;
z1 = magnitude*SIN(theta) + mu;
values(2*i - 1) = z0
IF (2*i <= n) values(2*i) = z1
END DO
CALL clip(values_range, values)
END SUBROUTINE random_normal_values_wp
END MODULE mo_random_number_generators
...@@ -53,56 +53,4 @@ CONTAINS ...@@ -53,56 +53,4 @@ CONTAINS
END SUBROUTINE open_logfile END SUBROUTINE open_logfile
FUNCTION calculate_mean_wp(array) RESULT(mean)
REAL(wp), INTENT(IN) :: array(:)
REAL(wp) :: mean
INTEGER :: i, n
mean = 0.0_wp
variance = 0.0_wp
n = SIZE(array)
DO i = 1, SIZE(array)
mean = mean + array(i)
END DO
mean = mean/n
END FUNCTION
FUNCTION calculate_variance_wp(array) RESULT(variance)
REAL(wp), INTENT(IN) :: array(:)
REAL(wp) :: mean, variance
INTEGER :: i, n
mean = 0.0_wp
variance = 0.0_wp
n = SIZE(array)
mean = calculate_mean_wp(array)
DO i = 1, SIZE(array)
variance = variance + (array(i) - mean)**2
END DO
variance = variance/n
END FUNCTION
SUBROUTINE assert_statistics(test_name, array, mean, variance, tol, max, min)
USE FORTUTF
CHARACTER(LEN=*), INTENT(IN) :: test_name
REAL(wp), INTENT(IN) :: array(:)
REAL(wp), INTENT(IN) :: mean, variance, tol, max, min
CALL TAG_TEST(test_name//"__mean")
CALL ASSERT_ALMOST_EQUAL(calculate_mean_wp(array), mean, tol)
CALL TAG_TEST(test_name//"__variance")
CALL ASSERT_ALMOST_EQUAL(calculate_variance_wp(array), variance, tol)
CALL TAG_TEST(test_name//"__max")
CALL ASSERT_LESS_THAN_EQUAL(MAXVAL(array), max)
CALL TAG_TEST(test_name//"__min")
CALL ASSERT_GREATER_THAN_EQUAL(MINVAL(array), min)
END SUBROUTINE
END MODULE helpers END MODULE helpers
! ICON
!
! ---------------------------------------------------------------
! Copyright (C) 2004-2024, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
! Contact information: icon-model.org
!
! See AUTHORS.TXT for a list of authors
! See LICENSES/ for license information
! SPDX-License-Identifier: BSD-3-Clause
! ---------------------------------------------------------------
MODULE TEST_mo_random_number_generators
USE FORTUTF
USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: wp => real64
USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: int32
USE helpers, ONLY: assert_statistics, calculate_mean_wp, calculate_variance_wp
CONTAINS
SUBROUTINE TEST_CALCULATE_STATISTICS_WP
REAL(wp) :: array(8) = (/2.0_wp, 4.0_wp, 4.0_wp, 4.0_wp, 5.0_wp, 5.0_wp, 7.0_wp, 9.0_wp/)
CALL TAG_TEST("TEST_calculate_mean_wp")
CALL ASSERT_ALMOST_EQUAL(calculate_mean_wp(array), 5.0_wp)
CALL TAG_TEST("TEST_calculate_variance_wp")
CALL ASSERT_ALMOST_EQUAL(calculate_variance_wp(array), 4.0_wp)
END SUBROUTINE
SUBROUTINE TEST_CLIP_WP
USE mo_random_number_generators
REAL(wp) :: array(6) = (/4.1_wp, -4.1_wp, 2.0_wp, -2.0_wp, 1.0_wp, -1.0_wp/)
CALL TAG_TEST("TEST_clip_wp")
CALL clip(2.0_wp, array)
CALL ASSERT_EQUAL(array, (/2.0_wp, -2.0_wp, 2.0_wp, -2.0_wp, 1.0_wp, -1.0_wp/))
END SUBROUTINE
SUBROUTINE TEST_UNIFORM_RANDOM_NUMBER_WP_BASIC
USE mo_random_number_generators
INTEGER, PARAMETER :: n = 1e5
INTEGER :: i
INTEGER(int32) :: rng_state
REAL(wp) :: array(n)
rng_state = initialize_random_number_generator(11, 17)
DO i = 1, n
array(i) = generate_uniform_random_number(rng_state)
END DO
CALL assert_statistics("TEST_UNIFORM_RANDOM_NUMBER_WP_BASIC", &
array, mean=0.5_wp, variance=(1.0_wp/12.0_wp), tol=1.0e-2_wp, max=1.0_wp, min=0.0_wp)
END SUBROUTINE
SUBROUTINE TEST_UNIFORM_RANDOM_NUMBER_WP_ISEED
USE mo_random_number_generators
INTEGER, PARAMETER :: n = 1e5
INTEGER :: i
INTEGER(int32) :: rng_state
REAL(wp) :: array(n)
DO i = 1, n
rng_state = initialize_random_number_generator(i, 11)
array(i) = generate_uniform_random_number(rng_state)
END DO
CALL assert_statistics("TEST_UNIFORM_RANDOM_NUMBER_WP_ISEED", &
array, mean=0.5_wp, variance=(1.0_wp/12.0_wp), tol=1.0e-2_wp, max=1.0_wp, min=0.0_wp)
END SUBROUTINE
SUBROUTINE TEST_UNIFORM_RANDOM_NUMBER_WP_JSEED
USE mo_random_number_generators
INTEGER, PARAMETER :: n = 1e5
INTEGER :: i
INTEGER(int32) :: rng_state
REAL(wp) :: array(n)
DO i = 1, n
rng_state = initialize_random_number_generator(11, i)
array(i) = generate_uniform_random_number(rng_state)
END DO
CALL assert_statistics("TEST_UNIFORM_RANDOM_NUMBER_WP_JSEED", &
array, mean=0.5_wp, variance=(1.0_wp/12.0_wp), tol=1.0e-2_wp, max=1.0_wp, min=0.0_wp)
END SUBROUTINE
SUBROUTINE TEST_RANDOM_NORMAL_VALUES_WP
USE mo_random_number_generators
INTEGER, PARAMETER :: n = 1e5 + 1
REAL(wp) :: array(n)
REAL(wp) :: mean = -5.0_wp
REAL(wp) :: variance = 7.1_wp
REAL(wp) :: value_range = 1.0e32_wp
CALL random_normal_values(seed=11, values_range=value_range, values=array, mean=mean, stdv=SQRT(variance))
CALL assert_statistics("TEST_RANDOM_NORMAL_VALUES_WP", &
array, mean=mean, variance=variance, tol=1.0e-2_wp, max=value_range, min=-value_range)
END SUBROUTINE
END MODULE TEST_mo_random_number_generators
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment