Commit 2e7e994a authored by Sebastian Müller's avatar Sebastian Müller 🐈
Browse files

Merge branch 'mo_kernel_to_forces' into 'develop'

include SMI dependencies

See merge request !42
parents 8951ca11 0aca1ba3
Pipeline #61583 passed with stages
in 9 minutes and 37 seconds
!> \file mo_integrate.f90
!> \copydoc mo_integrate
!> \brief Provides integration routines
!> \details This module provides routine for numerical integration such a Newton-Cotes formulas, etc.
!> \authors Matthias Cuntz
!> \date Mar 2013
MODULE mo_integrate
USE mo_kind, ONLY: i4, sp, dp
IMPLICIT NONE
PUBLIC :: int_regular ! Integrate regularily spaced data
! ------------------------------------------------------------------
!> \brief Integrate regularily spaced data.
!> \details Integrates regularily spaced data with a 5-point Newton-Cotes formula:
!! \f[ \int y dx = \frac{2}{45} \sum_{i=5,n-4,4} 7 y_{i-4} + 32 y_{i-3} + 12 y_{i-2} + 32 y_{i-1} + 7 y_{i} \f]
!!
!! dx=1 if not given.
!!
!! \b Example
!! \code{.f90}
!! vec = (/ 1., 2, 3., 4., 5., 6., 7., 8., 9. /)
!! m = int_regular(vec)
!! -> see also example in test directory
!! \endcode
!!
!> \param[in] "real(sp/dp) :: dat(:)" \f$ y_i \f$ 1D-array with y-values.
!> \param[in] "real(sp/dp) :: dx" x-spacing (default=1.)
!> \return real(sp/dp) :: integral — \f$ \int y \f$ integral of y values
!> \author Matthias Cuntz
!> \date Mar 2013
INTERFACE int_regular
MODULE PROCEDURE int_regular_sp, int_regular_dp
END INTERFACE int_regular
! ------------------------------------------------------------------
PRIVATE
! ------------------------------------------------------------------
CONTAINS
! ------------------------------------------------------------------
! integrates tabulated function with equal distant points usign 5-point Newton-Cotes formula
FUNCTION int_regular_dp(dat, dx)
IMPLICIT NONE
REAL(dp), DIMENSION(:), INTENT(IN) :: dat
REAL(dp), OPTIONAL, INTENT(IN) :: dx
REAL(dp) :: int_regular_dp
INTEGER(i4) :: n, n0
REAL(dp) :: ddx
if (size(dat,1) < 5) stop 'Error int_regular_dp: size(dat) < 5'
if (present(dx)) then
ddx = dx*2.0_dp/45.0_dp
else
ddx = 2.0_dp/45.0_dp
endif
n0 = 5
n = size(dat,1)
if (ddx .gt. 0.0_dp) then
int_regular_dp = sum( &
(7.0_dp*(dat(n0-4:n-4:4) + dat(n0:n:4)) + &
32.0_dp*(dat(n0-3:n-3:4) + dat(n0-1:n-1:4)) + &
12.0_dp*dat(n0-2:n-2:4)) )
! to avoid underflow issues
if ( ddx .lt. 1.0_dp ) then
if ( int_regular_dp .gt. tiny(1.0_dp)/ddx ) then
int_regular_dp = ddx * int_regular_dp
else
int_regular_dp = tiny(1.0_dp)
end if
else
int_regular_dp = ddx * int_regular_dp
end if
else
int_regular_dp = 0.0_dp
end if
END FUNCTION int_regular_dp
FUNCTION int_regular_sp(dat, dx)
IMPLICIT NONE
REAL(sp), DIMENSION(:), INTENT(IN) :: dat
REAL(sp), OPTIONAL, INTENT(IN) :: dx
REAL(sp) :: int_regular_sp
INTEGER(i4) :: n, n0
REAL(sp) :: ddx
if (size(dat,1) < 5) stop 'Error int_regular_sp: size(dat) < 5'
if (present(dx)) then
ddx = dx*2.0_sp/45.0_sp
else
ddx = 2.0_sp/45.0_sp
endif
n0 = 5
n = size(dat,1)
if (ddx .gt. 0.0_sp) then
int_regular_sp = sum( &
(7.0_sp*(dat(n0-4:n-4:4) + dat(n0:n:4)) + &
32.0_sp*(dat(n0-3:n-3:4) + dat(n0-1:n-1:4)) + &
12.0_sp*dat(n0-2:n-2:4)) )
! to avoid underflow issues
if ( ddx .lt. 1.0_sp ) then
if ( int_regular_sp .gt. tiny(1.0_sp)/ddx ) then
int_regular_sp = ddx * int_regular_sp
else
int_regular_sp = tiny(1.0_sp)
end if
else
int_regular_sp = ddx * int_regular_sp
end if
else
int_regular_sp = 0.0_sp
end if
END FUNCTION int_regular_sp
END MODULE mo_integrate
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
module test_mo_integrate
use funit
use mo_kind, only: i4, sp, dp
use mo_integrate, only: int_regular
implicit none
integer(i4), parameter :: n=9
contains
@test
subroutine test_integrate_dp()
real(dp), dimension(n) :: dat
real(dp) :: t=1.0E-5_dp
integer(i4) :: i
forall(i=1:n) dat(i) = real(i, dp)
@assertEqual(int_regular(dat), 40.0_dp, tolerance=t, message='int_regular double precision')
end subroutine test_integrate_dp
@test
subroutine test_integrate_sp()
real(sp), dimension(n) :: dat
real(sp) :: t=1.0E-5_sp
integer(i4) :: i
forall(i=1:n) dat(i) = real(i, sp)
@assertEqual(int_regular(dat), 40.0_sp, tolerance=t, message='int_regular single precision')
end subroutine test_integrate_sp
end module test_mo_integrate
module test_mo_kernel
use funit
use mo_kind, only: i4, dp
use mo_kernel, only: kernel_regression_h, kernel_regression
use mo_kernel, only: kernel_density_h, kernel_density
use mo_kernel, only: kernel_cumdensity
implicit none
integer(i4) :: ii, nn
real(dp), dimension(20) :: yin_1d
real(dp), dimension(10) :: test1, y
real(dp) :: hh
real(dp), dimension(2) :: h
real(dp), dimension(10,2) :: x
!> Test data: 20 temperature data of Guayaquil
real(dp), dimension(20), parameter :: xin_1d = (/ &
26.1_dp, 24.5_dp, 24.8_dp, 24.5_dp, 24.1_dp, &
24.3_dp, 26.4_dp, 24.9_dp, 23.7_dp, 23.5_dp, &
24.0_dp, 24.1_dp, 23.7_dp, 24.3_dp, 26.6_dp, &
24.6_dp, 24.8_dp, 24.4_dp, 26.8_dp, 25.2_dp /)
!> tolerance
real(dp), parameter :: tt=1.0E-4_dp
contains
@test
subroutine test_approx_bandwidth()
nn = size(x,1)
forall(ii=1:nn) x(ii,1) = real(ii-1,dp)/9.0_dp
forall(ii=1:nn) x(ii,2) = 1.0_dp / (real(ii-1,dp)/9.0_dp + 0.1_dp)
y(:) = 1.0_dp + x(:,1)**2 - sin(x(:,2))** 2
h = kernel_regression_h(x,y,silverman=.true.)
@assertEqual(h(1), 0.2291_dp, tolerance=tt, message='Bandwith h (Silvermans rule) Ha01')
@assertEqual(h(2), 1.9033_dp, tolerance=tt, message='Bandwith h (Silvermans rule) Ha02')
h = kernel_regression_h(x,y,silverman=.false.)
@assertEqual(h(1), 0.1726_dp, tolerance=tt, message='Bandwith h (Cross-validation) Ha03')
@assertEqual(h(2), 9.5169_dp, tolerance=tt, message='Bandwith h (Cross-validation) Ha04')
end subroutine test_approx_bandwidth
@test
subroutine test_regression()
nn = size(x,1)
forall(ii=1:nn) x(ii,1) = real(ii-1,dp)/9.0_dp
forall(ii=1:nn) x(ii,2) = 1.0_dp / (real(ii-1,dp)/9.0_dp + 0.1_dp)
y(:) = 1.0_dp + x(:,1)**2 - sin(x(:,2))** 2
y = kernel_regression(x,y,silverman=.false.)
test1 = (/ 5224.0_dp, 5256.0_dp, 5417.0_dp, 5178.0_dp, 4764.0_dp, 4923.0_dp, 6034.0_dp, 7774.0_dp, 9545.0_dp, 10960.0_dp /)
@assertEqual(y, test1 * tt, tolerance=tt, message='kernel regression using bandwith h')
end subroutine test_regression
@test
subroutine test_approx_bandwidth_density()
hh = kernel_density_h(xin_1d(:),silverman=.false.)
@assertEqual(hh, 0.4582_dp, tolerance=tt, message='bandwith h for density estimation')
end subroutine test_approx_bandwidth_density
@test
subroutine test_kernel_density()
hh = kernel_density_h(xin_1d(:),silverman=.false.)
yin_1d = kernel_density(xin_1d, h=hh, silverman=.false.)
test1 = (/ 1260.0_dp, 4713.0_dp, 3900.0_dp, 4713.0_dp, 4443.0_dp, 4773.0_dp, 1498.0_dp, 3494.0_dp, 3035.0_dp, 2214.0_dp /)
@assertEqual(yin_1d(1:10), test1 * tt, tolerance=tt, message='kernel prob. density (PDF) estimation')
end subroutine test_kernel_density
@test
subroutine test_kernel_cumdensity()
hh = kernel_density_h(xin_1d(:),silverman=.false.)
yin_1d = kernel_cumdensity(xin_1d, h=hh, silverman=.false., romberg=.false.)
test1 = (/ 8461.0_dp, 4744.0_dp, 6055.0_dp, 4744.0_dp, 2861.0_dp, 3789.0_dp, 8880.0_dp, 6425.0_dp, 1344.0_dp, 819.0_dp /)
@assertEqual(yin_1d(1:10), test1 * tt, tolerance=tt, message='kernel cum. density (CDF) estimation')
end subroutine test_kernel_cumdensity
end module test_mo_kernel
module test_mo_nelmin
use funit
use mo_kind, only: i4, dp
use mo_nelmin, only: nelmin, nelminxy
use mo_opt_functions, only: rosenbrock
implicit none
integer(i4), parameter :: n = 2
integer(i4) :: icount
integer(i4) :: ifault
integer(i4) :: kcount
integer(i4) :: konvge
integer(i4) :: numres
real(dp) :: reqmin
real(dp) :: start(n)
real(dp) :: step(n)
real(dp) :: xmin(n)
real(dp) :: ynewlo
real(dp) :: xx(1)
real(dp) :: yy(1)
real(dp), allocatable :: history(:)
integer(i4) :: i
!> tolerance
real(dp), parameter :: tt=1.0E-2_dp
contains
@test
subroutine test_nelmin()
start(1:n) = (/ -1.2_dp, 1.0_dp /)
reqmin = 1.0E-08_dp
step(1:n) = (/ 1.0_dp, 1.0_dp /)
konvge = 10
kcount = 500
xmin = nelmin( rosenbrock, start, funcmin=ynewlo, varmin=reqmin, step=step, &
konvge=konvge, maxeval=kcount, neval=icount, numrestart=numres, ierror=ifault, history=history)
@assertEqual(xmin, [1._dp, 1._dp], tolerance=tt, message='rosenbrock minimum 1')
start(1:n) = (/ -1.2_dp, 1.0_dp /)
reqmin = 1.0E-08_dp
step(1:n) = (/ 1.0_dp, 1.0_dp /)
konvge = 10
kcount = 500
xmin = nelmin(rosenbrock, start, varmin=reqmin, step=step, konvge=konvge, maxeval=kcount)
@assertEqual(xmin, [1._dp, 1._dp], tolerance=tt, message='rosenbrock minimum 2')
start(1:n) = (/ -1.2_dp, 1.0_dp /)
reqmin = 1.0E-08_dp
konvge = 10
kcount = 500
xmin = nelmin(rosenbrock, start, varmin=reqmin, konvge=konvge, maxeval=kcount)
@assertEqual(xmin, [1._dp, 1._dp], tolerance=tt, message='rosenbrock minimum 3')
start(1:n) = (/ -1.2_dp, 1.0_dp /)
konvge = 10
kcount = 500
xmin = nelmin(rosenbrock, start, konvge=konvge, maxeval=kcount)
@assertEqual(xmin, [1._dp, 1._dp], tolerance=tt, message='rosenbrock minimum 4')
start(1:n) = (/ -1.2_dp, 1.0_dp /)
xmin = nelmin(rosenbrock, start)
@assertEqual(xmin, [1._dp, 1._dp], tolerance=tt, message='rosenbrock minimum 5')
start(1:n) = (/ -12.0_dp, 10.0_dp /)
xmin = nelmin(rosenbrock, start)
@assertEqual(xmin, [1._dp, 1._dp], tolerance=tt, message='rosenbrock minimum 6')
start(1:n) = (/ -12.0_dp, 10.0_dp /)
xx(1) = 1.0_dp
yy(1) = 100.0_dp
xmin = nelminxy(rosenbrockxy, start, xx, yy, neval=icount, history=history)
@assertEqual(xmin, [1._dp, 1._dp], tolerance=tt, message='rosenbrock minimum 7')
end subroutine test_nelmin
function rosenbrockxy(p, xx, yy)
real(dp), dimension(:), intent(in) :: p
real(dp), dimension(:), intent(in) :: xx
real(dp), dimension(:), intent(in) :: yy
real(dp) :: rosenbrockxy
rosenbrockxy = rosenbrock(p)
end function rosenbrockxy
end module test_mo_nelmin
module test_mo_utils
use funit
use mo_kind, only: sp, dp, i4, i8, i1
use mo_utils, only: eq, ge, le, ne
use mo_utils, only: swap, locate !, cumsum, arange, linspace
use mo_utils, only: swap, locate, cumsum, arange, linspace
use mo_utils, only: is_finite, is_nan, is_normal, special_value, unpack_chunkwise
implicit none
......@@ -34,16 +34,19 @@ module test_mo_utils
integer(i4) :: i
logical :: compare
!> tolerance
real(dp), parameter :: tt=1.0E-5_dp
contains
@test
subroutine test_utils()
! -----------------------------------------------------
! DOUBLE PRECISON
! -----------------------------------------------------
! -----------------------------------------------------
! DOUBLE PRECISON
! -----------------------------------------------------
write(*,*) ''
write(*,*) 'Test: eq/ equal: dp'
......@@ -139,9 +142,9 @@ contains
write(*,'(E24.17,A4,E24.17,A5,L2)') a_dp,' >= ',b_dp,' --> ',compare
@assertFalse(compare)
! -----------------------------------------------------
! SINGLE PRECISON
! -----------------------------------------------------
! -----------------------------------------------------
! SINGLE PRECISON
! -----------------------------------------------------
write(*,*) ''
write(*,*) 'Test: eq/ equal: sp'
......@@ -237,8 +240,8 @@ contains
write(*,'(E15.8,A4,E15.8,A5,L2)') a_sp,' >= ',b_sp,' --> ',compare
@assertFalse(compare)
! -----------------------------------------------------
! Swap
! -----------------------------------------------------
! Swap
write(*,*) ''
write(*,*) 'Test: swap'
......@@ -274,9 +277,9 @@ contains
call swap(iat2, nn, 1)
@assertEqual(iat2, iat3)
! -----------------------------------------------------
! Locate
! -----------------------------------------------------
! Locate
write(*,*) ''
write(*,*) 'Test: locate'
......@@ -301,10 +304,10 @@ contains
s5 = (/ 0.1_sp, 5.5_sp, 10.1_sp, 50.5_sp, 200.1_sp /)
ii5 = locate(sat1, s5)
@assertEqual(ii5, (/0, 5, 10, 50, 100/))
! -----------------------------------------------------
! is_finite, is_nan, is_normal
! -----------------------------------------------------
! is_finite, is_nan, is_normal
write(*,*) ''
write(*,*) 'Test: is_finite, is_nan, is_normal'
......@@ -346,9 +349,9 @@ contains
!@assertFalse(is_normal(sat1(2)))
!@assertFalse(any(is_normal(sat1(1:2))))
! -----------------------------------------------------
! special_value
! -----------------------------------------------------
! special_value
write(*,*) ''
write(*,*) 'Test: special_value dp'
! TODO: fix this, it fails in gfortran debug (option "-ffpe-trap=zero,overflow,underflow")
......@@ -396,79 +399,75 @@ contains
@assertEqual(abs(special_value(1.0_sp, 'IEEE_POSITIVE_ZERO')), 0.0_sp)
@assertEqual(abs(special_value(1.0_sp, 'IEEE_NEGATIVE_ZERO')), 0.0_sp)
! ! -----------------------------------------------------
! ! COMMENTED since cumsum, arange and linspace are not existing anymore
! ! Cumsum
! ! double precision
! dat1(:) = 1.0_dp
! dat2 = cumsum(dat1)
! isgood = isgood .and. (eq(dat2(1),1.0_dp))
! isgood = isgood .and. (eq(dat2(nn),real(nn,dp)))
! ! single precision
! sat1(:) = 2.0_sp
! sat2 = cumsum(sat1)
! isgood = isgood .and. (eq(sat2(1),2.0_sp))
! isgood = isgood .and. (eq(sat2(nn),real(2*nn,sp)))
! ! integer
! iat1(:) = 3
! iat2 = cumsum(iat1)
! isgood = isgood .and. (iat2(1) == 3)
! isgood = isgood .and. (iat2(nn) == 3*nn)
! ! complex
! cat1(:) = (1.0_dp,1.0_dp)
! cat2 = cumsum(cat1)
! isgood = isgood .and. (eq(real(cat2(1)),1.0_dp))
! isgood = isgood .and. (eq(aimag(cat2(nn)),real(nn,dp)))
!
! ! -----------------------------------------------------
! ! Range
! ! double precision
! dat1(:) = 1.0_dp
! dat2 = cumsum(dat1)
! dat1 = arange(real(nn,dp))
! isgood = isgood .and. all(eq(dat1,dat2))
! ! single precision
! sat1(:) = 2.0_dp
! sat2 = cumsum(sat1)
! sat1 = arange(real(nn,sp)) * 2.0_sp
! isgood = isgood .and. all(eq(sat1,sat2))
! ! integer
! iat1(:) = 1
! iat2 = cumsum(iat1) - 2
! iat1 = arange(-1,nn-2)
! isgood = isgood .and. all(iat1==iat2)
! ! allocatable out
! adat1 = arange(real(nn,dp))
!
! ! -----------------------------------------------------
! ! Linspace
! ! double precision
! dat1 = arange(real(nn,dp))
! dat2 = linspace(1.0_dp,real(nn,dp),nn)
! isgood = isgood .and. all(eq(dat1,dat2))
! ! single precision
! sat1 = arange(real(nn,sp))/real(nn,sp)
! sat2 = linspace(0.01_sp,1.0_sp,nn)
! isgood = isgood .and. all(eq(sat1,sat2))
! ! integer
! iat1(:) = 3
! iat2 = cumsum(iat1)
! iat1 = linspace(3,3*nn,nn)
! isgood = isgood .and. all(iat1==iat2)
! -----------------------------------------------------
! unpack_chunkwise
! -----------------------------------------------------
write(*,*) ''
write(*,*) 'Test: unpack_chunkwise'
! -----------------------------------------------------
! Cumsum
! double precision
dat1(:) = 1.0_dp
dat2 = cumsum(dat1)
@assertEqual(dat2(1),1.0_dp, tolerance=tt)
@assertEqual(dat2(nn),real(nn,dp), tolerance=tt)
! single precision
sat1(:) = 2.0_sp
sat2 = cumsum(sat1)
@assertEqual(sat2(1),2.0_sp, tolerance=real(tt, sp))
@assertEqual(sat2(nn),real(2*nn,sp), tolerance=real(tt, sp))
! integer
iat1(:) = 3
iat2 = cumsum(iat1)
@assertEqual(iat2(1),3)
@assertEqual(iat2(nn),3*nn)
! complex
cat1(:) = (1.0_dp,1.0_dp)
cat2 = cumsum(cat1)
@assertEqual(real(cat2(1)),1.0_dp, tolerance=tt)
@assertEqual(aimag(cat2(nn)),real(nn,dp), tolerance=tt)
! -----------------------------------------------------
! Range
! double precision
dat1(:) = 1.0_dp
dat2 = cumsum(dat1)
dat1 = arange(real(nn,dp))
@assertEqual(dat1,dat2, tolerance=tt)
! single precision
sat1(:) = 2.0_dp
sat2 = cumsum(sat1)
sat1 = arange(real(nn,sp)) * 2.0_sp
@assertEqual(sat1,sat2, tolerance=real(tt, sp))
! integer
iat1(:) = 1
iat2 = cumsum(iat1) - 2
iat1 = arange(-1,nn-2)
@assertEqual(iat1,iat2)
! allocatable out
adat1 = arange(real(nn,dp))
! -----------------------------------------------------
! Linspace
! double precision
dat1 = arange(real(nn,dp))
dat2 = linspace(1.0_dp,real(nn,dp),nn)
@assertEqual(dat1,dat2, tolerance=tt)
! single precision
sat1 = arange(real(nn,sp))/real(nn,sp)
sat2 = linspace(0.01_sp,1.0_sp,nn)
@assertEqual(sat1,sat2, tolerance=real(tt, sp))
! integer
iat1(:) = 3
iat2 = cumsum(iat1)
iat1 = linspace(3,3*nn,nn)
@assertEqual(iat1,iat2)
! -----------------------------------------------------
! unpack_chunkwise
! -----------------------------------------------------
a_dp = -9999.0_dp
allocate(dalloc(5), lalloc(10))
dalloc = [ 0.1_dp, 5.5_dp, 10.1_dp, 50.5_dp, 200.1_dp ]
lalloc = [ .false., .true., .false., .false., .false., .false., .true., .true., .true., .true. ]
@assertEqual(unpack_chunkwise(dalloc, lalloc, a_dp, 2_i8), unpack(dalloc, lalloc, a_dp))
deallocate(lalloc, dalloc)
end subroutine test_utils
! TODO: this is commented as it takes really long to test, any ideas?
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment