Line data Source code
1 : !> \file mo_boxcox.f90
2 : !> \brief \copybrief mo_boxcox
3 : !> \details \copydetails mo_boxcox
4 :
5 : !> \brief Box-Cox transformation of data.
6 : !> \details This module contains routines to calculate the Box-Cox transformation
7 : !! as well as estimating the best exponent for the Box-Cox transformation
8 : !> \changelog
9 : !! - March 2011, Matthias Cuntz
10 : !! - modified Python code of Travis Oliphant (2002): boxcox, llf_boxcox
11 : !! - Dec 2019: Robert Schweppe:
12 : !! - removed NR code (get_boxcox)
13 : !> \author Mathias Cuntz
14 : !> \date Aug 2011
15 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
16 : !! FORCES is released under the LGPLv3+ license \license_note
17 : MODULE mo_boxcox
18 :
19 : USE mo_kind, ONLY : sp, dp
20 : USE mo_utils, only : le
21 :
22 : IMPLICIT NONE
23 :
24 : PUBLIC :: boxcox
25 : PUBLIC :: invboxcox
26 :
27 : ! ------------------------------------------------------------------
28 :
29 : !> \brief Transform a positive dataset with a Box-Cox power transformation.
30 :
31 : !> \details Calculate Box-Cox power transformed values given the original values and the exponent lambda.\n
32 : !! \f[ w_\text{Box-Cox}(x) =
33 : !! \begin{cases}
34 : !! \frac{x^\lambda - 1}{\lambda}&\text{, if }\lambda \neq 0 \\
35 : !! \ln{x}&\text{, if }\lambda = 0
36 : !! \end{cases} \f]
37 : !! If an optional mask is given, then the Box-Cox transformation is only performed on
38 : !! those locations that correspond to true values in the mask.\n
39 : !! \f$x\f$ can be single or double precision. The result will have the same numerical precision.
40 : !! \f$x\f$ can be scalar or vector.\n
41 : !!
42 : !! \b Example
43 : !!
44 : !! \code{.f90}
45 : !! out = boxcox(x, lmbda, mask=mask)
46 : !! \endcode
47 : !!
48 : !! See also test folder for a detailed example, "pf_tests/test_mo_boxcox".
49 :
50 :
51 : !> \param[in] "real(sp/dp) :: x" Scalar/1D-array with input numbers (`>0.`)
52 : !> \param[in] "real(sp/dp) :: lmbda" Exponent power of Box-Cox transform (`>= 0.`)
53 : !> \param[in] "logical, optional :: mask" Scalar/1D-array of logical values with `size(x)`.
54 : !! If present, only those locations in vec corresponding
55 : !! to the true values in mask are used.
56 : !> \retval "real(sp/dp) :: boxcox" Power transformed values (at `mask=.true.`)
57 :
58 : !> \note Input values must be positive, i.e. \f$x > 0\f$.
59 :
60 : !> \author Matthias Cuntz
61 : !> \date March 2011
62 : !! - Modified Python code of Travis Oliphant (2002): boxcox, llf_boxcox, get_boxcox
63 : !! - Modified numerical recipes: brent, mnbrak, swap, shft
64 : !> \date Dec 2021
65 : !! - Updated doxygen docs
66 : INTERFACE boxcox
67 : MODULE PROCEDURE boxcox_sp, boxcox_dp
68 : END INTERFACE boxcox
69 :
70 : ! ------------------------------------------------------------------
71 :
72 : !> \brief Back-transformation of Box-Cox-transformed data.
73 :
74 : !> \details Calculate the inverse Box-Cox given the transformed values and the exponent lambda.
75 : !! \f[ w_\text{Box-Cox}^{-1}(y) =
76 : !! \begin{cases}
77 : !! (\lambda y + 1)^{\frac{1}{\lambda}}&\text{, if }\lambda \neq 0 \\
78 : !! e^{y}&\text{, if }\lambda = 0
79 : !! \end{cases} \f]
80 : !! If an optional mask is given, then the inverse Box-Cox transformation is only performed on
81 : !! those locations that correspond to true values in the mask.\n
82 : !! \f$y\f$ can be single or double precision. The result will have the same numerical precision.\n
83 : !! \f$y\f$ can be scalar or vector.
84 : !!
85 : !! \b Example
86 : !!
87 : !! \code{.f90}
88 : !! out = invboxcox(x, lmbda, mask=mask)
89 : !! \endcode
90 : !!
91 : !! See also test folder for a detailed example, "pf_tests/test_mo_boxcox".
92 :
93 : !> \param[in] "real(sp/dp) :: y" Scalar/1D-array with input numbers (`>0.`)
94 : !> \param[in] "real(sp/dp) :: lmbda" Exponent power of Box-Cox transform (`>= 0.`)
95 : !> \param[in] "optional, logical :: mask" 1D-array of logical values with `size(x)`.
96 : !! If present, only those locations in vec corresponding to the
97 : !! true values in mask are used.
98 : !! Only applicable if `x` is a 1D-array.
99 : !> \retval "real(sp/dp) :: invboxcox" Back transformed values (at `mask=.true.`)
100 :
101 : !> \author Matthias Cuntz
102 : !> \author Juliane Mai
103 : !> \author Sebastian Müller
104 : !> \date March 2011
105 : !! - Modified MC: Python code of Travis Oliphant (2002): boxcox, llf_boxcox, get_boxcox
106 : !! - Modified MC: numerical recipes: brent, mnbrak, swap, shft
107 : !! - Modified JM: scalar version of invboxcox
108 : INTERFACE invboxcox
109 : MODULE PROCEDURE invboxcox_0d_sp, invboxcox_0d_dp, invboxcox_1d_sp, invboxcox_1d_dp
110 : END INTERFACE invboxcox
111 :
112 : ! ------------------------------------------------------------------
113 :
114 : PRIVATE
115 :
116 : ! ------------------------------------------------------------------
117 :
118 : CONTAINS
119 :
120 : ! ------------------------------------------------------------------
121 :
122 9 : FUNCTION boxcox_sp(x, lmbda, mask)
123 :
124 : IMPLICIT NONE
125 :
126 : REAL(sp), DIMENSION(:), INTENT(in) :: x
127 : REAL(sp), INTENT(in) :: lmbda
128 : LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask
129 : REAL(sp), DIMENSION(size(x)) :: boxcox_sp
130 :
131 3 : LOGICAL, DIMENSION(size(x)) :: maske
132 3 : REAL(sp) :: lmbda1
133 :
134 63 : maske(:) = .true.
135 3 : if (present(mask)) then
136 2 : if (size(mask) /= size(x)) stop 'Error boxcox_sp: size(mask) /= size(x)'
137 42 : maske = mask
138 : endif
139 63 : if (any((le(x, 0.0_sp)) .and. maske)) stop 'Error boxcox_sp: x <= 0'
140 3 : if (abs(lmbda) < tiny(0.0_sp)) then
141 21 : where (maske)
142 : boxcox_sp = log(x)
143 : elsewhere
144 : boxcox_sp = x
145 : end where
146 : else
147 2 : lmbda1 = 1.0_sp / lmbda
148 42 : boxcox_sp = merge((exp(lmbda * log(x)) - 1.0_sp) * lmbda1, x, maske)
149 : endif
150 :
151 3 : END FUNCTION boxcox_sp
152 :
153 :
154 12 : FUNCTION boxcox_dp(x, lmbda, mask)
155 :
156 : IMPLICIT NONE
157 :
158 : REAL(dp), DIMENSION(:), INTENT(in) :: x
159 : REAL(dp), INTENT(in) :: lmbda
160 : LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask
161 : REAL(dp), DIMENSION(size(x)) :: boxcox_dp
162 :
163 3 : LOGICAL, DIMENSION(size(x)) :: maske
164 3 : REAL(dp) :: lmbda1
165 :
166 63 : maske(:) = .true.
167 3 : if (present(mask)) then
168 2 : if (size(mask) /= size(x)) stop 'Error boxcox_dp: size(mask) /= size(x)'
169 42 : maske = mask
170 : endif
171 63 : if (any((le(x, 0.0_dp)) .and. maske)) then
172 0 : stop 'Error boxcox_dp: x <= 0'
173 : end if
174 3 : if (abs(lmbda) < tiny(0.0_dp)) then
175 21 : where (maske)
176 : boxcox_dp = log(x)
177 : elsewhere
178 : boxcox_dp = x
179 : end where
180 : else
181 2 : lmbda1 = 1.0_dp / lmbda
182 42 : boxcox_dp = merge((exp(lmbda * log(x)) - 1.0_dp) * lmbda1, x, maske)
183 : endif
184 :
185 3 : END FUNCTION boxcox_dp
186 :
187 : ! ------------------------------------------------------------------
188 :
189 12 : FUNCTION invboxcox_1d_sp(x, lmbda, mask)
190 :
191 : IMPLICIT NONE
192 :
193 : REAL(sp), DIMENSION(:), INTENT(in) :: x
194 : REAL(sp), INTENT(in) :: lmbda
195 : LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask
196 : REAL(sp), DIMENSION(size(x)) :: invboxcox_1d_sp
197 :
198 3 : LOGICAL, DIMENSION(size(x)) :: maske
199 63 : REAL(sp) :: lmbda1
200 66 : REAL(sp), DIMENSION(size(x)) :: temp
201 :
202 63 : maske(:) = .true.
203 3 : if (present(mask)) then
204 2 : if (size(mask) /= size(x)) stop 'Error invboxcox_1d_sp: size(mask) /= size(x)'
205 42 : maske = mask
206 : endif
207 3 : if (abs(lmbda) < tiny(0.0_sp)) then
208 21 : where (maske)
209 : invboxcox_1d_sp = exp(x)
210 : elsewhere
211 1 : invboxcox_1d_sp = x
212 : end where
213 : else
214 2 : lmbda1 = 1.0_sp / lmbda
215 42 : temp = lmbda * x + 1.0_sp
216 122 : where (temp > 0.0_sp)
217 : temp = exp(lmbda1 * log(temp))
218 : elsewhere
219 : temp = 0.0_sp
220 : end where
221 42 : invboxcox_1d_sp = merge(temp, x, maske)
222 : endif
223 :
224 3 : END FUNCTION invboxcox_1d_sp
225 :
226 12 : FUNCTION invboxcox_1d_dp(x, lmbda, mask)
227 :
228 : IMPLICIT NONE
229 :
230 : REAL(dp), DIMENSION(:), INTENT(in) :: x
231 : REAL(dp), INTENT(in) :: lmbda
232 : LOGICAL, DIMENSION(:), INTENT(in), OPTIONAL :: mask
233 : REAL(dp), DIMENSION(size(x)) :: invboxcox_1d_dp
234 :
235 3 : LOGICAL, DIMENSION(size(x)) :: maske
236 63 : REAL(dp) :: lmbda1
237 66 : REAL(dp), DIMENSION(size(x)) :: temp
238 :
239 63 : maske(:) = .true.
240 3 : if (present(mask)) then
241 2 : if (size(mask) /= size(x)) stop 'Error invboxcox_1d_dp: size(mask) /= size(x)'
242 42 : maske = mask
243 : endif
244 3 : if (abs(lmbda) < tiny(0.0_dp)) then
245 21 : where (maske)
246 : invboxcox_1d_dp = exp(x)
247 : elsewhere
248 1 : invboxcox_1d_dp = x
249 : end where
250 : else
251 2 : lmbda1 = 1.0_dp / lmbda
252 42 : temp = lmbda * x + 1.0_dp
253 122 : where (temp > 0.0_dp)
254 : temp = exp(lmbda1 * log(temp))
255 : elsewhere
256 : temp = 0.0_dp
257 : end where
258 42 : invboxcox_1d_dp = merge(temp, x, maske)
259 : endif
260 :
261 3 : END FUNCTION invboxcox_1d_dp
262 :
263 2 : FUNCTION invboxcox_0d_sp(x, lmbda)
264 :
265 : IMPLICIT NONE
266 :
267 : REAL(sp), INTENT(in) :: x
268 : REAL(sp), INTENT(in) :: lmbda
269 : REAL(sp) :: invboxcox_0d_sp
270 :
271 2 : REAL(sp) :: lmbda1
272 2 : REAL(sp) :: temp
273 :
274 2 : if (abs(lmbda) < tiny(0.0_sp)) then
275 1 : invboxcox_0d_sp = exp(x)
276 : else
277 1 : lmbda1 = 1.0_sp / lmbda
278 1 : temp = lmbda * x + 1.0_sp
279 1 : if (temp > 0.0_sp) then
280 1 : temp = exp(lmbda1 * log(temp))
281 : else
282 : temp = 0.0_sp
283 : end if
284 : invboxcox_0d_sp = temp
285 : endif
286 :
287 3 : END FUNCTION invboxcox_0d_sp
288 :
289 2 : FUNCTION invboxcox_0d_dp(x, lmbda)
290 :
291 : IMPLICIT NONE
292 :
293 : REAL(dp), INTENT(in) :: x
294 : REAL(dp), INTENT(in) :: lmbda
295 : REAL(dp) :: invboxcox_0d_dp
296 :
297 2 : REAL(dp) :: lmbda1
298 2 : REAL(dp) :: temp
299 :
300 2 : if (abs(lmbda) < tiny(0.0_dp)) then
301 1 : invboxcox_0d_dp = exp(x)
302 : else
303 1 : lmbda1 = 1.0_dp / lmbda
304 1 : temp = lmbda * x + 1.0_dp
305 1 : if (temp > 0.0_dp) then
306 1 : temp = exp(lmbda1 * log(temp))
307 : else
308 : temp = 0.0_dp
309 : end if
310 : invboxcox_0d_dp = temp
311 : endif
312 :
313 2 : END FUNCTION invboxcox_0d_dp
314 :
315 : END MODULE mo_boxcox
|
