Line data Source code
1 : !> \file mo_likelihood.f90
2 : !> \brief \copybrief mo_likelihood
3 : !> \details \copydetails mo_likelihood
4 :
5 : !> \brief Added for testing purposes of test_mo_mcmc
6 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
7 : !! FORCES is released under the LGPLv3+ license \license_note
8 : module mo_likelihood
9 :
10 : USE mo_kind, only: i4, dp
11 : USE mo_moment, only: stddev
12 : use mo_optimization_utils, only: eval_interface
13 :
14 : Implicit NONE
15 :
16 : PRIVATE
17 :
18 : PUBLIC :: data
19 : PUBLIC :: model_dp
20 : PUBLIC :: likelihood_dp, loglikelihood_dp ! "real" likelihood (sigma is an error model or given)
21 : PUBLIC :: likelihood_stddev_dp, loglikelihood_stddev_dp ! "faked" likelihood (sigma is computed by obs vs model)
22 : PUBLIC :: setmeas
23 :
24 : ! -------------------------------
25 : !> \brief synthetic data
26 : !> \details Data generated with
27 : !! paraset(1) = 1.0
28 : !! paraset(2) = 2.0
29 : !! paraset(3) = 3.0
30 : !! plus additive, Gaussian distributed error with mu=0.0 and sigma=0.5
31 : INTERFACE data
32 : MODULE PROCEDURE data_dp
33 : END INTERFACE data
34 :
35 : !> \brief set meas
36 : INTERFACE setmeas
37 : MODULE PROCEDURE setmeas_dp
38 : END INTERFACE setmeas
39 :
40 : REAL(DP),DIMENSION(100,2) :: meas ! measurements
41 : REAL(DP),PARAMETER :: stddev_global=0.5_dp ! standard deviation of measurement error
42 :
43 : ! ------------------------------------------------------------------
44 :
45 : CONTAINS
46 :
47 : ! -------------------------------
48 : !> \brief A Likelihood function: "real" likelihood (sigma is an error model or given)
49 0 : function likelihood_dp(paraset, eval, stddev_in, stddev_new, likeli_new)
50 : REAL(DP), DIMENSION(:), INTENT(IN) :: paraset ! parameter set
51 : procedure(eval_interface), INTENT(IN), pointer :: eval
52 : REAL(DP), INTENT(IN), optional :: stddev_in ! standard deviation of data
53 : REAL(DP), INTENT(OUT), OPTIONAL :: stddev_new ! standard deviation of errors using paraset
54 : REAL(DP), INTENT(OUT), OPTIONAL :: likeli_new ! likelihood using stddev_new,
55 : ! ! i.e. using new parameter set
56 : REAL(DP) :: likelihood_dp
57 :
58 : ! local
59 0 : REAL(DP), DIMENSION(size(meas,1)) :: errors
60 0 : real(dp), dimension(:,:), allocatable :: runoff
61 :
62 0 : call eval(paraset, runoff=runoff)
63 0 : errors(:) = runoff(:,1)-data()
64 0 : likelihood_dp = exp(-0.5_dp * sum( errors(:) * errors(:) / stddev_global**2 ))
65 :
66 0 : end function likelihood_dp
67 :
68 : ! -------------------------------
69 : !> \brief A Log-Likelihood function: "real" likelihood (sigma is an error model or given)
70 86763 : function loglikelihood_dp(paraset, eval, stddev_in, stddev_new, likeli_new)
71 : REAL(DP), DIMENSION(:), INTENT(IN) :: paraset ! parameter set
72 : procedure(eval_interface), INTENT(IN), pointer :: eval
73 : REAL(DP), INTENT(IN), optional :: stddev_in ! standard deviation of data
74 : REAL(DP), INTENT(OUT), OPTIONAL :: stddev_new ! standard deviation of errors using paraset
75 : REAL(DP), INTENT(OUT), OPTIONAL :: likeli_new ! likelihood using stddev_new,
76 : ! ! i.e. using new parameter set
77 : REAL(DP) :: loglikelihood_dp
78 :
79 : ! local
80 8763063 : REAL(DP), DIMENSION(size(meas,1)) :: errors
81 86763 : real(dp), dimension(:,:), allocatable :: runoff
82 :
83 86763 : call eval(paraset, runoff=runoff)
84 8763063 : errors(:) = runoff(:,1)-data()
85 8763063 : loglikelihood_dp = -0.5_dp * sum( errors(:) * errors(:) / stddev_global**2 )
86 :
87 86763 : end function loglikelihood_dp
88 :
89 : ! -------------------------------
90 : !> \brief A Likelihood function: "faked" likelihood (sigma is computed by obs vs model)
91 0 : function likelihood_stddev_dp(paraset, eval, stddev_in, stddev_new, likeli_new)
92 : REAL(DP), DIMENSION(:), INTENT(IN) :: paraset ! parameter set
93 : procedure(eval_interface), INTENT(IN), pointer :: eval
94 : REAL(DP), INTENT(IN), optional :: stddev_in ! standard deviation of data
95 : REAL(DP), INTENT(OUT), OPTIONAL :: stddev_new ! standard deviation of errors using paraset
96 : REAL(DP), INTENT(OUT), OPTIONAL :: likeli_new ! likelihood using stddev_new,
97 : ! ! i.e. using new parameter set
98 : REAL(DP) :: likelihood_stddev_dp
99 :
100 : ! local
101 0 : REAL(DP), DIMENSION(size(meas,1)) :: errors
102 0 : REAL(DP) :: stddev_err
103 0 : real(dp), dimension(:,:), allocatable :: runoff
104 :
105 0 : call eval(paraset, runoff=runoff)
106 0 : errors(:) = runoff(:,1)-data()
107 0 : likelihood_stddev_dp = exp(-0.5_dp * sum( errors(:) * errors(:) / stddev_in**2 ))
108 :
109 : ! optional out
110 0 : stddev_err = stddev(errors)
111 0 : if (present( stddev_new )) then
112 0 : stddev_new = stddev_err
113 : end if
114 0 : if (present( likeli_new )) then
115 0 : likeli_new = exp(-0.5_dp * sum( errors(:) * errors(:) / stddev_err**2 ))
116 : end if
117 :
118 86763 : end function likelihood_stddev_dp
119 :
120 : ! -------------------------------
121 : !> \brief A Log-Likelihood_stddev function: "faked" likelihood (sigma is computed by obs vs model)
122 182156 : function loglikelihood_stddev_dp(paraset, eval, stddev_in, stddev_new, likeli_new)
123 : REAL(DP), DIMENSION(:), INTENT(IN) :: paraset ! parameter set
124 : procedure(eval_interface), INTENT(IN), pointer :: eval
125 : REAL(DP), INTENT(IN), optional :: stddev_in ! standard deviation of data
126 : REAL(DP), INTENT(OUT), OPTIONAL :: stddev_new ! standard deviation of errors using paraset
127 : REAL(DP), INTENT(OUT), OPTIONAL :: likeli_new ! likelihood using stddev_new,
128 : ! ! i.e. using new parameter set
129 : REAL(DP) :: loglikelihood_stddev_dp
130 :
131 : ! local
132 18397756 : REAL(DP), DIMENSION(size(meas,1)) :: errors
133 182156 : REAL(DP) :: stddev_err
134 182156 : real(dp), dimension(:,:), allocatable :: runoff
135 :
136 182156 : call eval(paraset, runoff=runoff)
137 18397756 : errors(:) = runoff(:,1)-data()
138 18397756 : loglikelihood_stddev_dp = -0.5_dp * sum( errors(:) * errors(:) / stddev_in**2 )
139 :
140 : ! optional out
141 182156 : stddev_err = stddev(errors)
142 182156 : if (present( stddev_new )) then
143 39002 : stddev_new = stddev_err
144 : end if
145 182156 : if (present( likeli_new )) then
146 3939202 : likeli_new = -0.5_dp * sum( errors(:) * errors(:) / stddev_err**2 )
147 : end if
148 :
149 182156 : end function loglikelihood_stddev_dp
150 :
151 : ! -------------------------------
152 : !> \brief A Model: p1*x^2 + p2*x + p3
153 268919 : subroutine model_dp(parameterset, opti_domain_indices, runoff, smOptiSim, neutronsOptiSim, etOptiSim, twsOptiSim, &
154 : lake_level, lake_volume, lake_area, lake_spill, lake_outflow, BFI)
155 :
156 182156 : use mo_kind, only: dp
157 : use mo_optimization_types, only: optidata_sim
158 : !! !$ USE omp_lib, only: OMP_GET_THREAD_NUM
159 :
160 : real(dp), dimension(:), intent(in) :: parameterset
161 : integer(i4), dimension(:), optional, intent(in) :: opti_domain_indices
162 : real(dp), dimension(:, :), allocatable, optional, intent(out) :: runoff ! dim1=time dim2=gauge
163 : type(optidata_sim), dimension(:), optional, intent(inout) :: smOptiSim ! dim1=ncells, dim2=time
164 : type(optidata_sim), dimension(:), optional, intent(inout) :: neutronsOptiSim ! dim1=ncells, dim2=time
165 : type(optidata_sim), dimension(:), optional, intent(inout) :: etOptiSim ! dim1=ncells, dim2=time
166 : type(optidata_sim), dimension(:), optional, intent(inout) :: twsOptiSim ! dim1=ncells, dim2=time
167 : real(dp), dimension(:, :), allocatable, optional, intent(out) :: lake_level !< dim1=time dim2=lake
168 : real(dp), dimension(:, :), allocatable, optional, intent(out) :: lake_volume !< dim1=time dim2=lake
169 : real(dp), dimension(:, :), allocatable, optional, intent(out) :: lake_area !< dim1=time dim2=lake
170 : real(dp), dimension(:, :), allocatable, optional, intent(out) :: lake_spill !< dim1=time dim2=lake
171 : real(dp), dimension(:, :), allocatable, optional, intent(out) :: lake_outflow !< dim1=time dim2=lake
172 : real(dp), dimension(:), allocatable, optional, intent(out) :: BFI !< baseflow index, dim1=domainID
173 : integer(i4) :: i, n
174 : ! for OMP
175 : !! !$ integer(i4) :: n_threads, is_thread
176 :
177 268919 : n = size(meas,1)
178 268919 : allocate(runoff(n, 1))
179 :
180 : !! !$ is_thread = OMP_GET_THREAD_NUM()
181 : !! !$ write(*,*) 'OMP_thread: ', is_thread
182 :
183 : !$OMP parallel default(shared) &
184 : !$OMP private(i)
185 : !$OMP do
186 27160819 : do i=1, n
187 : !! !$ if (is_thread /= 0) write(*,*) ' OMP_thread-1: ', is_thread
188 27160819 : runoff(i,1) = parameterset(1) * meas(i,1) * meas(i,1) + parameterset(2) * meas(i,1) + parameterset(3)
189 : end do
190 : !$OMP end do
191 : !$OMP end parallel
192 :
193 268919 : end subroutine model_dp
194 :
195 537838 : function data_dp()
196 268919 : use mo_kind
197 : REAL(DP), DIMENSION(size(meas,1)) :: data_dp
198 :
199 27160819 : data_dp = meas(:,2)
200 :
201 268919 : end function data_dp
202 :
203 1 : subroutine setmeas_dp()
204 :
205 : integer(i4) :: i
206 :
207 101 : do i=1,100
208 101 : meas(i,1) = real(i,dp)
209 : end do
210 :
211 : meas(:,2) = (/ 5.49537_dp, 10.7835_dp, 17.6394_dp, 26.8661_dp, 36.9247_dp, 50.9517_dp, 66.2058_dp, &
212 : 82.9703_dp, 101.26_dp, 123.076_dp, 145.457_dp, 171.078_dp, 198.349_dp, 227.23_dp, 257.922_dp, &
213 : 290.098_dp, 325.775_dp, 362.724_dp, 402.669_dp, 442.461_dp, 486.122_dp, 531.193_dp, &
214 : 577.931_dp, 627.091_dp, 678.096_dp, 731.364_dp, 786.039_dp, 843.531_dp, 903.126_dp, &
215 : 963.037_dp, 1025.85_dp, 1091.85_dp, 1158.47_dp, 1226.65_dp, 1298.25_dp, 1370.87_dp, &
216 : 1444.96_dp, 1522.6_dp, 1602.6_dp, 1684.38_dp, 1765.15_dp, 1850.74_dp, 1937.85_dp, 2027.41_dp, &
217 : 2118.44_dp, 2210.62_dp, 2306.9_dp, 2403.27_dp, 2501.83_dp, 2602.96_dp, 2705.29_dp, &
218 : 2811.44_dp, 2917.82_dp, 3027.09_dp, 3137.64_dp, 3250.86_dp, 3366.67_dp, 3482.56_dp, &
219 : 3602.37_dp, 3722.55_dp, 3845.8_dp, 3970.15_dp, 4098.15_dp, 4227.27_dp, 4357.77_dp, &
220 : 4491.49_dp, 4626.23_dp, 4762.95_dp, 4901.15_dp, 5042.95_dp, 5184.86_dp, 5330.4_dp, &
221 : 5478.11_dp, 5626.97_dp, 5776.94_dp, 5931.15_dp, 6085.9_dp, 6242.7_dp, 6402.17_dp, 6563.46_dp, &
222 : 6726.37_dp, 6890.34_dp, 7058.4_dp, 7227.17_dp, 7397.09_dp, 7570.77_dp, 7745.95_dp, &
223 : 7923.72_dp, 8102.21_dp, 8282.98_dp, 8465.42_dp, 8651.37_dp, 8838.43_dp, 9027.57_dp, &
224 101 : 9219.21_dp, 9411.24_dp, 9606.31_dp, 9802.88_dp, 10002.3_dp, 10202.6_dp /)
225 268919 : end subroutine setmeas_dp
226 :
227 : end module mo_likelihood
|
