Skip to content
GitLab
Select Git revision
  • master default
  • fortran_module
  • dev
  • 1.1.0
  • 1.0.11
  • 1.0.7
  • 1.0.6
  • 1.0.2
  • 1.0.1
  • 1.0.0
  • 0.1.2
  • 0.1.1
  • 0.1.0
Find file BlameHistoryPermalink
similarity_module.f90 5.89 KiB
1 
module similarity_module
2 
  use iso_c_binding, only: c_int, c_double
3 
  use, intrinsic :: ieee_arithmetic
4 
  implicit none
5
6 
  abstract interface
7 
    function objective_fun(sim, obs, n) result(val)
8 
      use iso_c_binding, only: c_int, c_double
9 
      implicit none
10 
      integer(c_int), intent(in) :: n
11 
      real(c_double), intent(in) :: sim(n), obs(n)
12 
      real(c_double) :: val
13 
    end function objective_fun
14 
  end interface
15
16 
  procedure(objective_fun), pointer :: FUN => null()
17
18 
contains
19
20 
  subroutine set_objective_function(crit)
21 
    implicit none
22 
    integer(c_int), intent(in) :: crit
23
24 
    select case (crit)
25 
      case(1)
26 
        FUN => kge
27 
      case(2)
28 
        FUN => rmse
29 
      case(3)
30 
        FUN => invkge
31 
      case(4)
32 
        FUN => invrmse
33 
      case default
34 
        FUN => invkge
35 
    end select
36
37 
  end subroutine set_objective_function
38
39 
  subroutine similarity(Rn, nt, nb, crit, nthreads, sim_matrix) bind(C, name="similarity")
40 
    use, intrinsic :: ieee_arithmetic
41 
    use omp_lib
42 
    implicit none
43 
    integer(c_int), value :: nt, nb, crit, nthreads
44 
    real(c_double), intent(in) :: Rn(nt,nb)
45 
    real(c_double), intent(out) :: sim_matrix(nb,nb)
46
47 
    integer(c_int) :: i, j
48 
    real(c_double) :: val, nan_val
49 
    nan_val = ieee_value(0.0_c_double, ieee_quiet_nan)
50
51 
    call set_objective_function(crit)
52 
    if (.not. associated(FUN)) then
53 
      sim_matrix = nan_val
54 
      return
55 
    end if
56
57 
    select case (crit)
58 
      case(2, 4) ! if symetric, like RMSE
59 
        !$omp parallel do num_threads(nthreads) default(shared) private(i, j) schedule(static)
60 
        do i = 1, nb
61 
          do j = i, nb
62 
            sim_matrix(i,j) = FUN(Rn(:,i), Rn(:,j), nt)
63 
            sim_matrix(j,i) = sim_matrix(i,j)
64 
          end do
65 
        end do
66 
        !$omp end parallel do
67 
      case default ! if not symetric, like KGE
68 
        !$omp parallel do num_threads(nthreads) default(shared) private(i, j) schedule(static)
69 
        do i = 1, nb
70 
          do j = 1, nb
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
sim_matrix(i,j) = FUN(Rn(:,i), Rn(:,j), nt) end do end do !$omp end parallel do end select end subroutine similarity function kge(sim, obs, n) result(kge_val) use, intrinsic :: ieee_arithmetic implicit none integer, intent(in) :: n real(c_double), dimension(n), intent(in) :: sim, obs real(c_double) :: kge_val integer :: i, count real(c_double) :: sum_obs, sum_sim, sum_obs2, sum_sim2, sum_obs_sim real(c_double) :: mean_obs, mean_sim, std_obs, std_sim, r, beta, gamma, diff real(c_double) :: nan_val, eps eps = 1.0e-12_c_double count = 0 sum_obs = 0.0_c_double sum_sim = 0.0_c_double sum_obs2 = 0.0_c_double sum_sim2 = 0.0_c_double sum_obs_sim = 0.0_c_double nan_val = ieee_value(0.0_c_double, ieee_quiet_nan) ! Loop to filter NaN values and accumulate sums do i = 1, n if (.not. ieee_is_nan(obs(i)) .and. .not. ieee_is_nan(sim(i))) then count = count + 1 sum_obs = sum_obs + obs(i) sum_sim = sum_sim + sim(i) sum_obs2 = sum_obs2 + obs(i) * obs(i) sum_sim2 = sum_sim2 + sim(i) * sim(i) sum_obs_sim = sum_obs_sim + obs(i) * sim(i) end if end do ! Return NA if the common period is not long enough if (count <= 1) then kge_val = nan_val return end if ! Compute means mean_obs = sum_obs / count mean_sim = sum_sim / count ! Check to avoid division by zero if (abs(mean_obs) < eps) then kge_val = nan_val return end if ! Compute standard deviations (using Welford's correction) std_obs = sqrt((sum_obs2 - sum_obs * mean_obs) / (count - 1)) std_sim = sqrt((sum_sim2 - sum_sim * mean_sim) / (count - 1)) ! Check standard deviations if (std_obs < eps .or. std_sim < eps) then kge_val = nan_val return end if ! Compute correlation r = ((sum_obs_sim - count * mean_obs * mean_sim) / (count - 1)) / (std_obs * std_sim) ! Compute KGE components
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
beta = mean_sim / mean_obs gamma = (std_sim / mean_sim) / (std_obs / mean_obs) ! Final KGE computation diff = (r - 1.0_c_double)**2 + (beta - 1.0_c_double)**2 + (gamma - 1.0_c_double)**2 kge_val = 1.0_c_double - sqrt(diff) if (abs(kge_val - 1.0_c_double) < eps) kge_val = 1.0_c_double end function kge function rmse(sim, obs, n) result(res) integer, intent(in) :: n real(c_double), dimension(n), intent(in) :: sim, obs real(c_double) :: res integer :: i, count real(c_double) :: sum_sq, diff real(c_double) :: nan_val, eps nan_val = ieee_value(0.0_c_double, ieee_quiet_nan) eps = 1.0e-12_c_double count = 0 sum_sq = 0.0_c_double do i = 1, n if (.not. ieee_is_nan(obs(i)) .and. .not. ieee_is_nan(sim(i))) then count = count + 1 diff = sim(i) - obs(i) sum_sq = sum_sq + diff**2 end if end do if (count == 0) then res = nan_val return end if res = sqrt(sum_sq / count) if (res < eps) res = 0.0_c_double end function rmse function invrmse(sim, obs, n) result(res) integer, intent(in) :: n real(c_double), dimension(n), intent(in) :: sim, obs real(c_double) :: rmse_val, res real(c_double) :: inf_val, eps inf_val = ieee_value(0.0_c_double, ieee_positive_inf) eps = 1.0e-12_c_double rmse_val = rmse(obs, sim, n) if (rmse_val < eps) then res = inf_val else res = 1.0_c_double / rmse_val end if end function invrmse function invkge(sim, obs, n) result(res) integer, intent(in) :: n real(c_double), dimension(n), intent(in) :: sim, obs real(c_double) :: kge_val, res real(c_double) :: inf_val, eps inf_val = ieee_value(0.0_c_double, ieee_positive_inf) eps = 1.0e-12_c_double kge_val = kge(obs, sim, n) if (abs(kge_val - 1.0_c_double) < eps) then res = inf_val else res = 1.0_c_double / (1-kge_val) end if end function invkge
211212213214
end module similarity_module