From e84f14bea4220b71ea8c00cc15237006a50d7224 Mon Sep 17 00:00:00 2001 From: Yen-Chen Chen <yen-chen.chen@kit.edu> Date: Fri, 19 Jul 2024 09:24:50 +0000 Subject: [PATCH] Move random number generator to math-support (icon-libraries/libfortran-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: Daniel Hupp <daniel.hupp@meteoswiss.ch> Merged-by: Yen-Chen Chen <yen-chen.chen@kit.edu> Changelog: feature --- README.md | 1 - src/CMakeLists.txt | 1 - src/mo_random_number_generators.F90 | 152 ------------------ test/fortran/helpers.f90 | 52 ------ .../fortran/test_random_number_generators.f90 | 101 ------------ 5 files changed, 307 deletions(-) delete mode 100644 src/mo_random_number_generators.F90 delete mode 100644 test/fortran/test_random_number_generators.f90 diff --git a/README.md b/README.md index d341e99..52e0b10 100644 --- a/README.md +++ b/README.md @@ -51,7 +51,6 @@ The `libfortran-support` library includes some general Fortran supporting module - `mo_io_units`: io unit definitions - `mo_namelist`: open/close namelist files - `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_util_backtrace`: function backtrace - `mo_util_file`: file operations diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 18d8054..fc3c7d6 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -34,7 +34,6 @@ add_library( mo_io_units.F90 mo_namelist.F90 mo_octree.F90 - mo_random_number_generators.F90 mo_simple_dump.F90 mo_util_backtrace.F90 mo_util_file.F90 diff --git a/src/mo_random_number_generators.F90 b/src/mo_random_number_generators.F90 deleted file mode 100644 index e77e04a..0000000 --- a/src/mo_random_number_generators.F90 +++ /dev/null @@ -1,152 +0,0 @@ -! 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 diff --git a/test/fortran/helpers.f90 b/test/fortran/helpers.f90 index 602f262..6c260ba 100644 --- a/test/fortran/helpers.f90 +++ b/test/fortran/helpers.f90 @@ -53,56 +53,4 @@ CONTAINS 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 diff --git a/test/fortran/test_random_number_generators.f90 b/test/fortran/test_random_number_generators.f90 deleted file mode 100644 index 78f4932..0000000 --- a/test/fortran/test_random_number_generators.f90 +++ /dev/null @@ -1,101 +0,0 @@ -! 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 -- GitLab