Line data Source code
1 : !> \file mo_mad.f90
2 : !> \brief \copybrief mo_mad
3 : !> \details \copydetails mo_mad
4 :
5 : !> \brief Median absolute deviation test.
6 : !> \details This module provides a median absolute deviation (MAD) test.
7 : !> \authors Matthias Cuntz, Mathias Zink
8 : !> \date 2011 - 2012
9 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
10 : !! FORCES is released under the LGPLv3+ license \license_note
11 : MODULE mo_mad
12 :
13 : USE mo_kind, ONLY: i4, sp, dp
14 : USE mo_percentile, ONLY: median
15 :
16 : Implicit NONE
17 :
18 : PUBLIC :: mad ! Mean absolute deviation test
19 :
20 : ! ------------------------------------------------------------------
21 :
22 : !> \brief Mean absolute deviation test.
23 :
24 : !> \details Mean absolute deviation test with optional z-value (default: 7) and mask,
25 : !! and 3 optional variants:
26 : !!
27 : !! - 0: raw values (default)
28 : !! - 1: 1st derivative
29 : !! - 2: 2nd derivative
30 : !!
31 : !! Returns mask with true everywhere except where `0<(median-MAD*z/0.6745)` or `>(md+MAD*z/0.6745)` \n
32 : !!
33 : !! If an optional mask is given, the mad test is only performed on those locations that correspond
34 : !! to true values in the mask.\n
35 : !!
36 : !! If tout is given mad returns the array with the enteries exceeding the treshold
37 : !! being set to the threshold. With this setting arrays are accepted. \n
38 : !! tout accepts:
39 : !!
40 : !! - u: upper values are cut at the threshold
41 : !! - l: lower values are cut at the threshold
42 : !! - b: upper and lower values are cut at the threshold
43 : !!
44 : !! With this setting only the variant 0 is available (no argument implemented).
45 : !!
46 : !! \b Example
47 : !!
48 : !! \code{.f90}
49 : !! vec = (/ -0.25,0.68,0.94,1.15,2.26,2.35,2.37,2.40,2.47,2.54,2.62, &
50 : !! 2.64,2.90,2.92,2.92,2.93,3.21,3.26,3.30,3.59,3.68,4.30, &
51 : !! 4.64,5.34,5.42,8.01 /)
52 : !! mask(:) = true
53 : !! ! Sets last entry false
54 : !! mask = mask .and. mad(vec, z=4., deriv=0, mask=mask)
55 : !! \endcode
56 : !!
57 : !! See also example in test directory.
58 :
59 :
60 : !> \param[in] "real(sp/dp) :: vec(:)" 1D-array with input numbers
61 : !> \param[in] "real(sp/dp) :: arr" nD-array with input numbers
62 : !> \param[in] "real(sp/dp), optional :: z" Input is allowed to deviate maximum z standard deviations
63 : !! from the median (default: 7).
64 : !> \param[in] "integer(i4), optional :: deriv" 0: Act on raw input (default: 0)\n
65 : !! 1: Use first derivatives\n
66 : !! 2: Use 2nd derivatives
67 : !> \param[in] "logical, optional :: mask(:)" 1D-array of logical values with size(vec).
68 : !! If present, only those locations in vec corresponding to the true
69 : !! values in mask are used. nD-array if tout is used.
70 : !> \param[in] "character(1):: tout" u: Trim only values above mad\n
71 : !! l: Trim only values below mad\n
72 : !! b: Trim values below and above mad
73 :
74 : !> \retval "logical or real(sp/dp), dimension(size(arr)) :: out" mask with true everywhere except where input
75 : !! deviates more than z standard deviations from
76 : !! median. If tout five, threshold value is returned.
77 :
78 : !> \note
79 : !! 1st derivative is
80 : !!
81 : !! d = datin[1:n]-datin[0:n-1]
82 : !!
83 : !! because mean of left and right would give 0 for spikes.
84 : !!
85 : !> \author Matthias Cuntz
86 : !> \date Mar 2011
87 :
88 : !> \author Mathhias Kelbling
89 : !> \date May 2018
90 : !! - mad_val added by Matthias Kelbling, May 2018
91 :
92 : ! ------------------------------------------------------------------
93 :
94 : INTERFACE mad
95 : MODULE PROCEDURE mad_sp, mad_dp, mad_val_dp, mad_val_sp
96 : END INTERFACE mad
97 :
98 : PRIVATE
99 :
100 : ! ------------------------------------------------------------------
101 :
102 : CONTAINS
103 :
104 : ! ------------------------------------------------------------------
105 :
106 0 : FUNCTION mad_dp(arr, z, mask, deriv)
107 :
108 : IMPLICIT NONE
109 :
110 : REAL(dp), DIMENSION(:), INTENT(IN) :: arr
111 : REAL(dp), OPTIONAL, INTENT(IN) :: z
112 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
113 : INTEGER(i4), OPTIONAL, INTENT(IN) :: deriv
114 : LOGICAL, DIMENSION(size(arr)) :: mad_dp
115 :
116 0 : LOGICAL, DIMENSION(size(arr)) :: maske
117 0 : REAL(dp), DIMENSION(size(arr)) :: d
118 0 : LOGICAL, DIMENSION(size(arr)) :: dmask
119 : INTEGER(i4) :: n, ideriv
120 : !INTEGER(i4) :: m
121 0 : REAL(dp) :: iz, med, mabsdev, thresh
122 :
123 0 : n = size(arr)
124 0 : maske(:) = .true.
125 0 : if (present(mask)) then
126 0 : if (size(mask) /= n) stop 'Error mad_dp: size(mask) /= size(arr)'
127 0 : maske = mask
128 : endif
129 0 : if (present(z)) then
130 0 : iz = z
131 : else
132 : iz = 7.0_dp
133 : endif
134 0 : if (present(deriv)) then
135 0 : ideriv = deriv
136 : else
137 : ideriv = 0
138 : endif
139 :
140 0 : select case(ideriv)
141 : case(0) ! pure values
142 : !m = count(maske)
143 0 : med = median(arr,mask=maske)
144 0 : mabsdev = median(abs(arr-med),mask=maske)
145 0 : thresh = mabsdev * iz/0.6745_dp
146 0 : mad_dp = (arr .ge. (med-thresh)) .and. (arr .le. (med+thresh)) .and. maske
147 : case(1) ! 1st derivative
148 : ! d(1:n-2) = 0.5_dp* (arr(3:n) - arr(1:n-2)) ! does not work -> ask Clemens & Matthias M
149 0 : d(1:n-1) = arr(2:n) - arr(1:n-1)
150 0 : dmask(1:n-1) = maske(2:n) .and. maske(1:n-1)
151 : !m = count(dmask(1:n-1))
152 0 : med = median(d(1:n-1),mask=dmask(1:n-1))
153 0 : mabsdev = median(abs(d(1:n-1)-med),mask=dmask(1:n-1))
154 0 : thresh = mabsdev * iz/0.6745_dp
155 0 : mad_dp(n) = .true.
156 0 : mad_dp(1:n-1) = (d(1:n-1) .ge. (med-thresh)) .and. (d(1:n-1) .le. (med+thresh)) .and. dmask(1:n-1)
157 : case(2) ! -2nd derivative
158 0 : d(1:n-2) = arr(2:n-1) + arr(2:n-1) - arr(1:n-2) - arr(3:n)
159 0 : dmask(1:n-2) = maske(2:n-1) .and. maske(1:n-2) .and. maske(3:n)
160 : !m = count(dmask(1:n-2))
161 0 : med = median(d(1:n-2),mask=dmask(1:n-2))
162 0 : mabsdev = median(abs(d(1:n-2)-med),mask=dmask(1:n-2))
163 0 : thresh = mabsdev * iz/0.6745_dp
164 0 : mad_dp(1) = .true.
165 0 : mad_dp(n) = .true.
166 0 : mad_dp(2:n-1) = (d(1:n-2) .ge. (med-thresh)) .and. (d(1:n-2) .le. (med+thresh)) .and. dmask(1:n-2)
167 : case default
168 0 : stop 'Unimplemented option in mad_dp'
169 : end select
170 :
171 0 : END FUNCTION mad_dp
172 :
173 :
174 0 : FUNCTION mad_sp(arr, z, mask, deriv)
175 :
176 : IMPLICIT NONE
177 :
178 : REAL(sp), DIMENSION(:), INTENT(IN) :: arr
179 : REAL(sp), OPTIONAL, INTENT(IN) :: z
180 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
181 : INTEGER(i4), OPTIONAL, INTENT(IN) :: deriv
182 : LOGICAL, DIMENSION(size(arr)) :: mad_sp
183 :
184 0 : LOGICAL, DIMENSION(size(arr)) :: maske
185 0 : REAL(sp), DIMENSION(size(arr)) :: d
186 0 : LOGICAL, DIMENSION(size(arr)) :: dmask
187 : INTEGER(i4) :: n, ideriv
188 : !INTEGER(i4) :: m
189 0 : REAL(sp) :: iz, med, mabsdev, thresh
190 :
191 0 : n = size(arr)
192 0 : maske(:) = .true.
193 0 : if (present(mask)) then
194 0 : if (size(mask) /= n) stop 'Error mad_sp: size(mask) /= size(arr)'
195 0 : maske = mask
196 : endif
197 0 : if (present(z)) then
198 0 : iz = z
199 : else
200 : iz = 7.0_sp
201 : endif
202 0 : if (present(deriv)) then
203 0 : ideriv = deriv
204 : else
205 : ideriv = 0
206 : endif
207 :
208 0 : select case(ideriv)
209 : case(0) ! pure values
210 : !m = count(maske)
211 0 : med = median(arr,mask=maske)
212 0 : mabsdev = median(abs(arr-med),mask=maske)
213 0 : thresh = mabsdev * iz/0.6745_sp
214 0 : mad_sp = (arr .ge. (med-thresh)) .and. (arr .le. (med+thresh)) .and. maske
215 : case(1) ! 1st derivative
216 : ! d(1:n-2) = 0.5_sp* (arr(3:n) - arr(1:n-2)) ! does not work -> ask Clemens & Matthias M
217 0 : d(1:n-1) = arr(2:n) - arr(1:n-1)
218 0 : dmask(1:n-1) = maske(2:n) .and. maske(1:n-1)
219 : !m = count(dmask(1:n-1))
220 0 : med = median(d(1:n-1),mask=dmask(1:n-1))
221 0 : mabsdev = median(abs(d(1:n-1)-med),mask=dmask(1:n-1))
222 0 : thresh = mabsdev * iz/0.6745_sp
223 0 : mad_sp(n) = .true.
224 0 : mad_sp(1:n-1) = (d(1:n-1) .ge. (med-thresh)) .and. (d(1:n-1) .le. (med+thresh)) .and. dmask(1:n-1)
225 : case(2) ! -2nd derivative
226 0 : d(1:n-2) = arr(2:n-1) + arr(2:n-1) - arr(1:n-2) - arr(3:n)
227 0 : dmask(1:n-2) = maske(2:n-1) .and. maske(1:n-2) .and. maske(3:n)
228 : !m = count(dmask(1:n-2))
229 0 : med = median(d(1:n-2),mask=dmask(1:n-2))
230 0 : mabsdev = median(abs(d(1:n-2)-med),mask=dmask(1:n-2))
231 0 : thresh = mabsdev * iz/0.6745_sp
232 0 : mad_sp(1) = .true.
233 0 : mad_sp(n) = .true.
234 0 : mad_sp(2:n-1) = (d(1:n-2) .ge. (med-thresh)) .and. (d(1:n-2) .le. (med+thresh)) .and. dmask(1:n-2)
235 : case default
236 0 : stop 'Unimplemented option in mad_sp'
237 : end select
238 :
239 0 : END FUNCTION mad_sp
240 :
241 : ! ------------------------------------------------------------------
242 :
243 0 : FUNCTION mad_val_dp(arr, z, mask, tout, mval)
244 :
245 : IMPLICIT NONE
246 :
247 : REAL(dp), DIMENSION(:), INTENT(IN) :: arr
248 : REAL(dp), OPTIONAL, INTENT(IN) :: z, mval
249 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
250 : REAL(dp), DIMENSION(size(arr)) :: mad_val_dp
251 : CHARACTER(1) :: tout ! type out
252 : ! u : cut upper; l : cut lower; b : cut upper and lower
253 :
254 0 : LOGICAL, DIMENSION(size(arr)) :: maske
255 : INTEGER(i4) :: n
256 0 : REAL(dp) :: iz, med, mabsdev, thresh
257 :
258 0 : n = size(arr)
259 0 : maske(:) = .true.
260 0 : if (present(mask)) then
261 0 : if (size(mask) /= n) stop 'Error mad_val_dp: size(mask) /= size(arr)'
262 0 : maske = mask
263 : endif
264 0 : if (present(z)) then
265 0 : iz = z
266 : else
267 : iz = 7.0_dp
268 : endif
269 :
270 0 : if (present(mval)) then
271 0 : where (abs(arr - mval) .lt. tiny(1._dp) ) maske = .false.
272 : ! reset if no values remain
273 0 : if (.not. any(maske)) then
274 0 : where ( abs(arr - mval) .lt. tiny(1._dp) ) maske = .true.
275 : end if
276 : endif
277 :
278 0 : med = median(arr,mask=maske)
279 0 : mabsdev = median(abs(arr-med),mask=maske)
280 0 : thresh = mabsdev * iz/0.6745_dp
281 0 : mad_val_dp = arr
282 :
283 : select case(tout)
284 : case("u")
285 : ! print *, "The threshold is set to", med, "+", thresh
286 : where ((mad_val_dp .gt. (med+thresh)) &
287 0 : .and. maske) mad_val_dp = med+thresh
288 : case("l")
289 : ! print *, "The threshold is set to", med, "-", thresh
290 : where ((mad_val_dp .lt. (med-thresh)) &
291 0 : .and. maske) mad_val_dp = med-thresh
292 : case("b")
293 : ! print *, "The threshold is set to", med, "+/-", thresh
294 : where ((mad_val_dp .gt. (med+thresh)) &
295 0 : .and. maske) mad_val_dp = med+thresh
296 : where ((mad_val_dp .lt. (med-thresh)) &
297 0 : .and. maske) mad_val_dp = med-thresh
298 : case default
299 0 : stop 'Unimplemented option in mad_val_dp'
300 : end select
301 :
302 0 : END FUNCTION mad_val_dp
303 :
304 : ! ------------------------------------------------------------------
305 :
306 0 : FUNCTION mad_val_sp(arr, z, mask, tout, mval)
307 :
308 : IMPLICIT NONE
309 :
310 : REAL(sp), DIMENSION(:), INTENT(IN) :: arr
311 : REAL(sp), OPTIONAL, INTENT(IN) :: z, mval
312 : LOGICAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: mask
313 : REAL(sp), DIMENSION(size(arr)) :: mad_val_sp
314 : CHARACTER(1) :: tout ! type out
315 : ! u : cut upper; l : cut lower; b : cut upper and lower
316 :
317 0 : LOGICAL, DIMENSION(size(arr)) :: maske
318 : INTEGER(i4) :: n
319 0 : REAL(sp) :: iz, med, mabsdev, thresh
320 :
321 0 : n = size(arr)
322 0 : maske(:) = .true.
323 0 : if (present(mask)) then
324 0 : if (size(mask) /= n) stop 'Error mad_val_sp: size(mask) /= size(arr)'
325 0 : maske = mask
326 : endif
327 0 : if (present(z)) then
328 0 : iz = z
329 : else
330 : iz = 7.0_sp
331 : endif
332 :
333 0 : if (present(mval)) then
334 0 : where (abs(arr - mval) .lt. tiny(1._sp)) maske = .false.
335 : ! reset if no values remain
336 0 : if (.not. any(maske)) then
337 0 : where ( abs(arr - mval) .lt. tiny(1._dp) ) maske = .true.
338 : end if
339 : endif
340 :
341 0 : med = median(arr,mask=maske)
342 0 : mabsdev = median(abs(arr-med),mask=maske)
343 0 : thresh = mabsdev * iz/0.6745_sp
344 0 : mad_val_sp = arr
345 0 : select case(tout)
346 : case("u")
347 0 : print *, "The threshold is set to", med, "+", thresh
348 : where ((mad_val_sp .gt. (med+thresh)) &
349 0 : .and. maske) mad_val_sp = med+thresh
350 : case("l")
351 0 : print *, "The threshold is set to", med, "-", thresh
352 : where ((mad_val_sp .lt. (med-thresh)) &
353 0 : .and. maske) mad_val_sp = med-thresh
354 : case("b")
355 0 : print *, "The threshold is set to", med, "+/-", thresh
356 : where ((mad_val_sp .gt. (med+thresh)) &
357 0 : .and. maske) mad_val_sp = med+thresh
358 : where ((mad_val_sp .lt. (med-thresh)) &
359 0 : .and. maske) mad_val_sp = med-thresh
360 : case default
361 0 : stop 'Unimplemented option in mad_val_sp'
362 : end select
363 :
364 0 : END FUNCTION mad_val_sp
365 :
366 : ! ------------------------------------------------------------------
367 :
368 :
369 : END MODULE mo_mad
|
