libsim  Versione7.1.6
simple_stat.f90
1 ! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2 ! authors:
3 ! Davide Cesari <dcesari@arpa.emr.it>
4 ! Paolo Patruno <ppatruno@arpa.emr.it>
5 
6 ! This program is free software; you can redistribute it and/or
7 ! modify it under the terms of the GNU General Public License as
8 ! published by the Free Software Foundation; either version 2 of
9 ! the License, or (at your option) any later version.
10 
11 ! This program is distributed in the hope that it will be useful,
12 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
13 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 ! GNU General Public License for more details.
15 
16 ! You should have received a copy of the GNU General Public License
17 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
29 IMPLICIT NONE
30 
39 INTERFACE stat_average
40  MODULE PROCEDURE stat_averager, stat_averaged
41 END INTERFACE
42 
54 INTERFACE stat_variance
55  MODULE PROCEDURE stat_variancer, stat_varianced
56 END INTERFACE
57 
69 INTERFACE stat_stddev
70  MODULE PROCEDURE stat_stddevr, stat_stddevd
71 END INTERFACE
72 
90  MODULE PROCEDURE stat_linear_corrr, stat_linear_corrd
91 END INTERFACE
92 
108  MODULE PROCEDURE stat_linear_regressionr, stat_linear_regressiond
109 END INTERFACE
110 
122  MODULE PROCEDURE stat_percentiler, stat_percentiled
123 END INTERFACE
124 
141 INTERFACE stat_bin
142  MODULE PROCEDURE stat_binr, stat_bind
143 END INTERFACE stat_bin
144 
161  MODULE PROCEDURE stat_mode_histogramr, stat_mode_histogramd
162 END INTERFACE stat_mode_histogram
163 
164 PRIVATE
167  normalizeddensityindex
168 
169 CONTAINS
170 
171 
172 FUNCTION stat_averager(sample, mask, nomiss) RESULT(average)
173 REAL,INTENT(in) :: sample(:)
174 LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
175 LOGICAL,OPTIONAL,INTENT(in) :: nomiss ! informs that the sample does not contain missing data and is not of zero length (for speedup)
176 
177 REAL :: average
178 
179 INTEGER :: sample_count
180 LOGICAL :: sample_mask(size(sample))
181 
182 IF (optio_log(nomiss)) THEN
183  average = sum(sample)/SIZE(sample)
184 ELSE
185  sample_mask = (sample /= rmiss)
186  IF (present(mask)) sample_mask = sample_mask .AND. mask
187  sample_count = count(sample_mask)
188 
189  IF (sample_count > 0) THEN
190 ! compute average
191  average = sum(sample, mask=sample_mask)/sample_count
192  ELSE
193  average = rmiss
194  ENDIF
195 ENDIF
196 
197 END FUNCTION stat_averager
198 
199 
200 FUNCTION stat_averaged(sample, mask, nomiss) RESULT(average)
201 DOUBLE PRECISION,INTENT(in) :: sample(:)
202 LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
203 LOGICAL,OPTIONAL,INTENT(in) :: nomiss ! informs that the sample does not contain missing data and is not of zero length (for speedup)
204 
205 DOUBLE PRECISION :: average
206 
207 INTEGER :: sample_count
208 LOGICAL :: sample_mask(size(sample))
209 
210 IF (optio_log(nomiss)) THEN
211  average = sum(sample)/SIZE(sample)
212 ELSE
213  sample_mask = (sample /= dmiss)
214  IF (present(mask)) sample_mask = sample_mask .AND. mask
215  sample_count = count(sample_mask)
216 
217  IF (sample_count > 0) THEN
218 ! compute average
219  average = sum(sample, mask=sample_mask)/sample_count
220  ELSE
221  average = dmiss
222  ENDIF
223 ENDIF
224 
225 END FUNCTION stat_averaged
226 
227 
228 FUNCTION stat_variancer(sample, average, mask, nomiss, nm1) RESULT(variance)
229 REAL,INTENT(in) :: sample(:)
230 REAL,OPTIONAL,INTENT(out) :: average
231 LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
232 LOGICAL,OPTIONAL,INTENT(in) :: nomiss
233 LOGICAL,OPTIONAL,INTENT(in) :: nm1
234 
235 REAL :: variance
236 
237 REAL :: laverage
238 INTEGER :: sample_count, i
239 LOGICAL :: sample_mask(size(sample))
240 
241 IF (optio_log(nomiss)) THEN
242 ! compute average
243  laverage = sum(sample)/SIZE(sample)
244  IF (present(average)) average = laverage
245  IF (optio_log(nm1)) THEN
246  variance = sum((sample-laverage)**2)/max(SIZE(sample)-1,1)
247  ELSE
248  variance = sum((sample-laverage)**2)/SIZE(sample)
249  ENDIF
250 
251 ELSE
252  sample_mask = (sample /= rmiss)
253  IF (present(mask)) sample_mask = sample_mask .AND. mask
254  sample_count = count(sample_mask)
255 
256  IF (sample_count > 0) THEN
257 ! compute average
258  laverage = sum(sample, mask=sample_mask)/sample_count
259  IF (present(average)) average = laverage
260 ! compute sum of squares and variance
261  variance = 0.
262  DO i = 1, SIZE(sample)
263  IF (sample_mask(i)) variance = variance + (sample(i)-laverage)**2
264  ENDDO
265  IF (optio_log(nm1)) THEN
266  variance = variance/max(sample_count-1,1)
267  ELSE
268  variance = variance/sample_count
269  ENDIF
270  ELSE
271  IF (present(average)) average = rmiss
272  variance = rmiss
273  ENDIF
274 
275 ENDIF
276 
277 END FUNCTION stat_variancer
278 
279 
280 FUNCTION stat_varianced(sample, average, mask, nomiss, nm1) RESULT(variance)
281 DOUBLE PRECISION,INTENT(in) :: sample(:)
282 DOUBLE PRECISION,OPTIONAL,INTENT(out) :: average
283 LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
284 LOGICAL,OPTIONAL,INTENT(in) :: nomiss
285 LOGICAL,OPTIONAL,INTENT(in) :: nm1
286 
287 DOUBLE PRECISION :: variance
288 
289 DOUBLE PRECISION :: laverage
290 INTEGER :: sample_count, i
291 LOGICAL :: sample_mask(size(sample))
292 
293 IF (optio_log(nomiss)) THEN
294 ! compute average
295  laverage = sum(sample)/SIZE(sample)
296  IF (present(average)) average = laverage
297  IF (optio_log(nm1)) THEN
298  variance = sum((sample-laverage)**2)/max(SIZE(sample)-1,1)
299  ELSE
300  variance = sum((sample-laverage)**2)/SIZE(sample)
301  ENDIF
302 
303 ELSE
304  sample_mask = (sample /= dmiss)
305  IF (present(mask)) sample_mask = sample_mask .AND. mask
306  sample_count = count(sample_mask)
307 
308  IF (sample_count > 0) THEN
309 ! compute average
310  laverage = sum(sample, mask=sample_mask)/sample_count
311  IF (present(average)) average = laverage
312 ! compute sum of squares and variance
313  variance = 0.
314  DO i = 1, SIZE(sample)
315  IF (sample_mask(i)) variance = variance + (sample(i)-laverage)**2
316  ENDDO
317  IF (optio_log(nm1)) THEN
318  variance = variance/max(sample_count-1,1)
319  ELSE
320  variance = variance/sample_count
321  ENDIF
322  ELSE
323  IF (present(average)) average = dmiss
324  variance = dmiss
325  ENDIF
326 
327 ENDIF
328 
329 END FUNCTION stat_varianced
330 
331 
332 FUNCTION stat_stddevr(sample, average, mask, nomiss, nm1) RESULT(stddev)
333 REAL,INTENT(in) :: sample(:)
334 REAL,OPTIONAL,INTENT(out) :: average
335 LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
336 LOGICAL,OPTIONAL,INTENT(in) :: nomiss
337 LOGICAL,OPTIONAL,INTENT(in) :: nm1
338 
339 REAL :: stddev
340 
341 stddev = stat_variance(sample, average, mask, nomiss, nm1)
342 IF (c_e(stddev)) stddev = sqrt(stddev)
343 
344 END FUNCTION stat_stddevr
345 
346 
347 FUNCTION stat_stddevd(sample, average, mask, nomiss, nm1) RESULT(stddev)
348 DOUBLE PRECISION,INTENT(in) :: sample(:)
349 DOUBLE PRECISION,OPTIONAL,INTENT(out) :: average
350 LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
351 LOGICAL,OPTIONAL,INTENT(in) :: nomiss
352 LOGICAL,OPTIONAL,INTENT(in) :: nm1
353 
354 DOUBLE PRECISION :: stddev
355 
356 stddev = stat_variance(sample, average, mask, nomiss, nm1)
357 IF (c_e(stddev)) stddev = sqrt(stddev)
358 
359 END FUNCTION stat_stddevd
360 
361 
362 FUNCTION stat_linear_corrr(sample1, sample2, average1, average2, &
363  variance1, variance2, mask, nomiss) result(linear_corr)
364 REAL,INTENT(in) :: sample1(:)
365 REAL,INTENT(in) :: sample2(:)
366 REAL,OPTIONAL,INTENT(out) :: average1
367 REAL,OPTIONAL,INTENT(out) :: average2
368 REAL,OPTIONAL,INTENT(out) :: variance1
369 REAL,OPTIONAL,INTENT(out) :: variance2
370 LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
371 LOGICAL,OPTIONAL,INTENT(in) :: nomiss
372 
373 REAL :: linear_corr
374 
375 REAL :: laverage1, laverage2, lvariance1, lvariance2
376 INTEGER :: sample_count, i
377 LOGICAL :: sample_mask(size(sample1))
378 
379 IF (SIZE(sample1) /= SIZE(sample2)) THEN
380  IF (present(average1)) average1 = rmiss
381  IF (present(average2)) average2 = rmiss
382  IF (present(variance1)) variance1 = rmiss
383  IF (present(variance2)) variance2 = rmiss
384  linear_corr = rmiss
385  RETURN
386 ENDIF
387 
388 sample_mask = (sample1 /= rmiss .AND. sample2 /= rmiss)
389 IF (present(mask)) sample_mask = sample_mask .AND. mask
390 sample_count = count(sample_mask)
391 IF (sample_count > 0) THEN
392 ! compute averages
393  laverage1 = sum(sample1, mask=sample_mask)/sample_count
394  laverage2 = sum(sample2, mask=sample_mask)/sample_count
395  IF (present(average1)) average1 = laverage1
396  IF (present(average2)) average2 = laverage2
397 ! compute sum of squares and variances
398  lvariance1 = 0.
399  lvariance2 = 0.
400  DO i = 1, SIZE(sample1)
401  IF (sample_mask(i))THEN
402  lvariance1 = lvariance1 + (sample1(i)-laverage1)**2
403  lvariance2 = lvariance2 + (sample2(i)-laverage2)**2
404  ENDIF
405  ENDDO
406  lvariance1 = lvariance1/sample_count
407  lvariance2 = lvariance2/sample_count
408  IF (present(variance1)) variance1 = lvariance1
409  IF (present(variance2)) variance2 = lvariance2
410 ! compute correlation
411  linear_corr = 0.
412  DO i = 1, SIZE(sample1)
413  IF (sample_mask(i)) linear_corr = linear_corr + &
414  (sample1(i)-laverage1)*(sample2(i)-laverage2)
415  ENDDO
416  linear_corr = linear_corr/sample_count / sqrt(lvariance1*lvariance2)
417 ELSE
418  IF (present(average1)) average1 = rmiss
419  IF (present(average2)) average2 = rmiss
420  IF (present(variance1)) variance1 = rmiss
421  IF (present(variance2)) variance2 = rmiss
422  linear_corr = rmiss
423 ENDIF
424 
425 END FUNCTION stat_linear_corrr
426 
427 
428 FUNCTION stat_linear_corrd(sample1, sample2, average1, average2, &
429  variance1, variance2, mask, nomiss) result(linear_corr)
430 DOUBLE PRECISION, INTENT(in) :: sample1(:)
431 DOUBLE PRECISION, INTENT(in) :: sample2(:)
432 DOUBLE PRECISION, OPTIONAL, INTENT(out) :: average1
433 DOUBLE PRECISION, OPTIONAL, INTENT(out) :: average2
434 DOUBLE PRECISION, OPTIONAL, INTENT(out) :: variance1
435 DOUBLE PRECISION, OPTIONAL, INTENT(out) :: variance2
436 LOGICAL, OPTIONAL, INTENT(in) :: mask(:)
437 LOGICAL, OPTIONAL, INTENT(in) :: nomiss
438 
439 DOUBLE PRECISION :: linear_corr
440 
441 DOUBLE PRECISION :: laverage1, laverage2, lvariance1, lvariance2
442 INTEGER :: sample_count, i
443 LOGICAL :: sample_mask(size(sample1))
444 
445 IF (SIZE(sample1) /= SIZE(sample2)) THEN
446  IF (present(average1)) average1 = dmiss
447  IF (present(average2)) average2 = dmiss
448  IF (present(variance1)) variance1 = dmiss
449  IF (present(variance2)) variance2 = dmiss
450  linear_corr = dmiss
451  RETURN
452 ENDIF
453 
454 sample_mask = (sample1 /= dmiss .AND. sample2 /= dmiss)
455 IF (present(mask)) sample_mask = sample_mask .AND. mask
456 sample_count = count(sample_mask)
457 IF (sample_count > 0) THEN
458 ! compute averages
459  laverage1 = sum(sample1, mask=sample_mask)/sample_count
460  laverage2 = sum(sample2, mask=sample_mask)/sample_count
461  IF (present(average1)) average1 = laverage1
462  IF (present(average2)) average2 = laverage2
463 ! compute sum of squares and variances
464  lvariance1 = 0.
465  lvariance2 = 0.
466  DO i = 1, SIZE(sample1)
467  IF (sample_mask(i))THEN
468  lvariance1 = lvariance1 + (sample1(i)-laverage1)**2
469  lvariance2 = lvariance2 + (sample2(i)-laverage2)**2
470  ENDIF
471  ENDDO
472  lvariance1 = lvariance1/sample_count
473  lvariance2 = lvariance2/sample_count
474  IF (present(variance1)) variance1 = lvariance1
475  IF (present(variance2)) variance2 = lvariance2
476 ! compute correlation
477  linear_corr = 0.0d0
478  DO i = 1, SIZE(sample1)
479  IF (sample_mask(i)) linear_corr = linear_corr + &
480  (sample1(i)-laverage1)*(sample2(i)-laverage2)
481  ENDDO
482  linear_corr = linear_corr/sample_count / sqrt(lvariance1*lvariance2)
483 ELSE
484  IF (present(average1)) average1 = dmiss
485  IF (present(average2)) average2 = dmiss
486  IF (present(variance1)) variance1 = dmiss
487  IF (present(variance2)) variance2 = dmiss
488  linear_corr = dmiss
489 ENDIF
490 
491 END FUNCTION stat_linear_corrd
492 
493 
494 SUBROUTINE stat_linear_regressionr(sample1, sample2, alpha0, alpha1, mask)
495 REAL,INTENT(in) :: sample1(:)
496 REAL,INTENT(in) :: sample2(:)
497 REAL,INTENT(out) :: alpha0
498 REAL,INTENT(out) :: alpha1
499 LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
500 
501 REAL :: laverage1, laverage2
502 INTEGER :: sample_count
503 LOGICAL :: sample_mask(size(sample1))
504 
505 IF (SIZE(sample1) /= SIZE(sample2)) THEN
506  alpha0 = rmiss
507  alpha1 = rmiss
508  RETURN
509 ENDIF
510 
511 sample_mask = (sample1 /= rmiss .AND. sample2 /= rmiss)
512 IF (present(mask)) sample_mask = sample_mask .AND. mask
513 sample_count = count(sample_mask)
514 
515 IF (sample_count > 0) THEN
516  laverage1 = sum(sample1, mask=sample_mask)/sample_count
517  laverage2 = sum(sample2, mask=sample_mask)/sample_count
518  alpha1 = sum((sample1-laverage1)*(sample2-laverage2), mask=sample_mask)/ &
519  sum((sample1-laverage1)**2, mask=sample_mask)
520  alpha0 = laverage1 - alpha1*laverage2
521 
522 ELSE
523  alpha0 = rmiss
524  alpha1 = rmiss
525 ENDIF
526 
527 END SUBROUTINE stat_linear_regressionr
528 
529 
530 SUBROUTINE stat_linear_regressiond(sample1, sample2, alpha0, alpha1, mask)
531 DOUBLE PRECISION,INTENT(in) :: sample1(:)
532 DOUBLE PRECISION,INTENT(in) :: sample2(:)
533 DOUBLE PRECISION,INTENT(out) :: alpha0
534 DOUBLE PRECISION,INTENT(out) :: alpha1
535 LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
536 
537 DOUBLE PRECISION :: laverage1, laverage2
538 INTEGER :: sample_count
539 LOGICAL :: sample_mask(size(sample1))
540 
541 IF (SIZE(sample1) /= SIZE(sample2)) THEN
542  alpha0 = dmiss
543  alpha1 = dmiss
544  RETURN
545 ENDIF
546 
547 sample_mask = (sample1 /= dmiss .AND. sample2 /= dmiss)
548 IF (present(mask)) sample_mask = sample_mask .AND. mask
549 sample_count = count(sample_mask)
550 
551 IF (sample_count > 0) THEN
552  laverage1 = sum(sample1, mask=sample_mask)/sample_count
553  laverage2 = sum(sample2, mask=sample_mask)/sample_count
554  alpha1 = sum((sample1-laverage1)*(sample2-laverage2), mask=sample_mask)/ &
555  sum((sample1-laverage1)**2, mask=sample_mask)
556  alpha0 = laverage1 - alpha1*laverage2
557 
558 ELSE
559  alpha0 = dmiss
560  alpha1 = dmiss
561 ENDIF
562 
563 END SUBROUTINE stat_linear_regressiond
564 
565 
566 FUNCTION stat_percentiler(sample, perc_vals, mask, nomiss) RESULT(percentile)
567 REAL,INTENT(in) :: sample(:)
568 REAL,INTENT(in) :: perc_vals(:)
569 LOGICAL,OPTIONAL,INTENT(in) :: mask(:)
570 LOGICAL,OPTIONAL,INTENT(in) :: nomiss
571 REAL :: percentile(size(perc_vals))
572 
573 REAL :: lsample(size(sample)), rindex
574 INTEGER :: sample_count, j, iindex
575 LOGICAL :: sample_mask(size(sample))
576 
577 percentile(:) = rmiss
578 IF (.NOT.optio_log(nomiss)) THEN
579  sample_mask = c_e(sample)
580  IF (present(mask)) sample_mask = sample_mask .AND. mask
581  sample_count = count(sample_mask)
582  IF (sample_count == 0) RETURN ! special case
583 ELSE
584  sample_count = SIZE(sample)
585 ENDIF
586 
587 IF (sample_count == SIZE(sample)) THEN
588  lsample = sample
589 ELSE
590  lsample(1:sample_count) = pack(sample, mask=sample_mask)
591 ENDIF
592 
593 IF (sample_count == 1) THEN ! other special case
594  percentile(:) = lsample(1)
595  RETURN
596 ENDIF
597 
598 
599 ! this sort is very fast but with a lot of equal values it is very slow and fails
600 CALL sort(lsample(1:sample_count))
601 
602 ! this sort is very very slow
603 !!$DO j = 2, sample_count
604 !!$ v = lsample(j)
605 !!$ DO i = j-1, 1, -1
606 !!$ IF (v >= lsample(i)) EXIT
607 !!$ lsample(i+1) = lsample(i)
608 !!$ ENDDO
609 !!$ lsample(i+1) = v
610 !!$ENDDO
611 
612 DO j = 1, SIZE(perc_vals)
613  IF (perc_vals(j) >= 0. .AND. perc_vals(j) <= 100.) THEN
614 ! compute real index of requested percentile in sample
615  rindex = REAL(sample_count-1, kind=kind(rindex))*perc_vals(j)/100.+1.
616 ! compute integer index of previous element in sample, beware of corner cases
617  iindex = min(max(int(rindex), 1), sample_count-1)
618 ! compute linearly interpolated percentile
619 
620  percentile(j) = lsample(iindex)*(REAL(iindex+1, kind=kind(rindex))-rindex) &
621  + lsample(iindex+1)*(rindex-REAL(iindex, kind=kind(rindex)))
622 
623  ENDIF
624 ENDDO
625 
626 END FUNCTION stat_percentiler
627 
628 
629 FUNCTION stat_percentiled(sample, perc_vals, mask, nomiss) RESULT(percentile)
630 DOUBLE PRECISION, INTENT(in) :: sample(:)
631 DOUBLE PRECISION, INTENT(in) :: perc_vals(:)
632 LOGICAL, OPTIONAL, INTENT(in) :: mask(:)
633 LOGICAL, OPTIONAL, INTENT(in) :: nomiss
634 DOUBLE PRECISION :: percentile(size(perc_vals))
635 
636 DOUBLE PRECISION :: lsample(size(sample)), rindex
637 INTEGER :: sample_count, j, iindex
638 LOGICAL :: sample_mask(size(sample))
639 
640 
641 percentile(:) = dmiss
642 IF (.NOT.optio_log(nomiss)) THEN
643  sample_mask = (sample /= dmiss)
644  IF (present(mask)) sample_mask = sample_mask .AND. mask
645  sample_count = count(sample_mask)
646  IF (sample_count == 0) RETURN ! particular case
647 ELSE
648  sample_count = SIZE(sample)
649 ENDIF
650 
651 IF (sample_count == SIZE(sample)) THEN
652  lsample = sample
653 ELSE
654  lsample(1:sample_count) = pack(sample, mask=sample_mask)
655 ENDIF
656 
657 IF (sample_count == 1) THEN ! other particular case
658  percentile(:) = lsample(1)
659  RETURN
660 ENDIF
661 
662 ! this sort is very fast but with a lot of equal values it is very slow and fails
663 CALL sort(lsample(1:sample_count))
664 
665 ! this sort is very very slow
666 !DO j = 2, sample_count
667 ! v = lsample(j)
668 ! DO i = j-1, 1, -1
669 ! IF (v >= lsample(i)) EXIT
670 ! lsample(i+1) = lsample(i)
671 ! ENDDO
672 ! lsample(i+1) = v
673 !ENDDO
674 
675 DO j = 1, SIZE(perc_vals)
676  IF (perc_vals(j) >= 0.d0 .AND. perc_vals(j) <= 100.d0) THEN
677 ! compute real index of requested percentile in sample
678  rindex = REAL(sample_count-1, kind=kind(rindex))*perc_vals(j)/100.d0+1.d0
679 ! compute integer index of previous element in sample, beware of corner cases
680  iindex = min(max(int(rindex), 1), sample_count-1)
681 ! compute linearly interpolated percentile
682  percentile(j) = lsample(iindex)*(REAL(iindex+1, kind=kind(rindex))-rindex) &
683  + lsample(iindex+1)*(rindex-REAL(iindex, kind=kind(rindex)))
684  ENDIF
685 ENDDO
686 
687 END FUNCTION stat_percentiled
688 
689 
690 SUBROUTINE stat_binr(sample, bin, nbin, start, finish, mask, binbounds)
691 REAL,INTENT(in) :: sample(:)
692 INTEGER,INTENT(out),ALLOCATABLE :: bin(:)
693 INTEGER,INTENT(in) :: nbin
694 REAL,INTENT(in),OPTIONAL :: start
695 REAL,INTENT(in),OPTIONAL :: finish
696 LOGICAL,INTENT(in),OPTIONAL :: mask(:)
697 REAL,INTENT(out),ALLOCATABLE,OPTIONAL :: binbounds(:)
698 
699 INTEGER :: i, ind
700 REAL :: lstart, lfinish, incr
701 REAL,ALLOCATABLE :: lbinbounds(:)
702 LOGICAL :: lmask(size(sample))
703 
704 ! safety checks
705 IF (nbin < 1) RETURN
706 IF (present(mask)) THEN
707  lmask = mask
708 ELSE
709  lmask =.true.
710 ENDIF
711 lmask = lmask .AND. c_e(sample)
712 IF (count(lmask) < 1) RETURN
713 
714 lstart = optio_r(start)
715 IF (.NOT.c_e(lstart)) lstart = minval(sample, mask=lmask)
716 lfinish = optio_r(finish)
717 IF (.NOT.c_e(lfinish)) lfinish = maxval(sample, mask=lmask)
718 IF (lfinish <= lstart) RETURN
719 
720 incr = (lfinish-lstart)/nbin
721 ALLOCATE(bin(nbin))
722 
723 ALLOCATE(lbinbounds(nbin+1))
724 
725 DO i = 1, nbin
726  lbinbounds(i) = lstart + (i-1)*incr
727 ENDDO
728 lbinbounds(nbin+1) = lfinish ! set exact value to avoid truncation error
729 
730 DO i = 1, nbin-1
731  bin(i) = count(sample >= lbinbounds(i) .AND. sample < lbinbounds(i+1) .AND. lmask)
732 !, .AND. c_e(sample))
733 ENDDO
734 bin(nbin) = count(sample >= lbinbounds(nbin) .AND. sample <= lbinbounds(nbin+1) .AND. lmask)
735 !, .AND. c_e(sample)) ! special case, include upper limit
736 
737 IF (present(binbounds)) binbounds = lbinbounds
738 
739 END SUBROUTINE stat_binr
740 
741 
742 ! alternative algorithm, probably faster with big samples, to be tested
743 SUBROUTINE stat_binr2(sample, bin, nbin, start, finish, mask, binbounds)
744 REAL,INTENT(in) :: sample(:)
745 INTEGER,INTENT(out),ALLOCATABLE :: bin(:)
746 INTEGER,INTENT(in) :: nbin
747 REAL,INTENT(in),OPTIONAL :: start
748 REAL,INTENT(in),OPTIONAL :: finish
749 LOGICAL,INTENT(in),OPTIONAL :: mask(:)
750 REAL,INTENT(out),ALLOCATABLE,OPTIONAL :: binbounds(:)
751 
752 INTEGER :: i, ind
753 REAL :: lstart, lfinish, incr
754 LOGICAL :: lmask(size(sample))
755 
756 ! safety checks
757 IF (nbin < 1) RETURN
758 IF (present(mask)) THEN
759  lmask = mask
760 ELSE
761  lmask =.true.
762 ENDIF
763 lmask = lmask .AND. c_e(sample)
764 IF (count(lmask) < 1) RETURN
765 
766 lstart = optio_r(start)
767 IF (.NOT.c_e(lstart)) lstart = minval(sample, mask=c_e(sample))
768 lfinish = optio_r(finish)
769 IF (.NOT.c_e(lfinish)) lfinish = maxval(sample, mask=c_e(sample))
770 IF (lfinish <= lstart) RETURN
771 
772 incr = (lfinish-lstart)/nbin
773 ALLOCATE(bin(nbin))
774 
775 bin(:) = 0
776 
777 DO i = 1, SIZE(sample)
778  IF (lmask(i)) THEN
779  ind = int((sample(i)-lstart)/incr) + 1
780  IF (ind > 0 .AND. ind <= nbin) THEN
781  bin(ind) = bin(ind) + 1
782  ELSE
783  IF (sample(i) == finish) bin(nbin) = bin(nbin) + 1 ! special case, include upper limit
784  ENDIF
785  ENDIF
786 ENDDO
787 
788 IF (present(binbounds)) THEN
789  ALLOCATE(binbounds(nbin+1))
790  DO i = 1, nbin
791  binbounds(i) = lstart + (i-1)*incr
792  ENDDO
793  binbounds(nbin+1) = lfinish ! set exact value to avoid truncation error
794 ENDIF
795 
796 END SUBROUTINE stat_binr2
797 
798 
799 SUBROUTINE stat_bind(sample, bin, nbin, start, finish, mask, binbounds)
800 DOUBLE PRECISION,INTENT(in) :: sample(:)
801 INTEGER,INTENT(out),ALLOCATABLE :: bin(:)
802 INTEGER,INTENT(in) :: nbin
803 DOUBLE PRECISION,INTENT(in),OPTIONAL :: start
804 DOUBLE PRECISION,INTENT(in),OPTIONAL :: finish
805 LOGICAL,INTENT(in),OPTIONAL :: mask(:)
806 DOUBLE PRECISION,INTENT(out),ALLOCATABLE,OPTIONAL :: binbounds(:)
807 
808 INTEGER :: i, ind
809 DOUBLE PRECISION :: lstart, lfinish, incr
810 DOUBLE PRECISION,ALLOCATABLE :: lbinbounds(:)
811 LOGICAL :: lmask(size(sample))
812 
813 ! safety checks
814 IF (nbin < 1) RETURN
815 IF (present(mask)) THEN
816  lmask = mask
817 ELSE
818  lmask =.true.
819 ENDIF
820 lmask = lmask .AND. c_e(sample)
821 IF (count(lmask) < 1) RETURN
822 
823 lstart = optio_d(start)
824 IF (.NOT.c_e(lstart)) lstart = minval(sample, mask=lmask)
825 lfinish = optio_d(finish)
826 IF (.NOT.c_e(lfinish)) lfinish = maxval(sample, mask=lmask)
827 IF (lfinish <= lstart) RETURN
828 
829 incr = (lfinish-lstart)/nbin
830 ALLOCATE(bin(nbin))
831 
832 ALLOCATE(lbinbounds(nbin+1))
833 
834 DO i = 1, nbin
835  lbinbounds(i) = lstart + (i-1)*incr
836 ENDDO
837 lbinbounds(nbin+1) = lfinish ! set exact value to avoid truncation error
838 
839 DO i = 1, nbin-1
840  bin(i) = count(sample >= lbinbounds(i) .AND. sample < lbinbounds(i+1) .AND. lmask)
841 !, .AND. c_e(sample))
842 ENDDO
843 bin(nbin) = count(sample >= lbinbounds(nbin) .AND. sample <= lbinbounds(nbin+1) .AND. lmask)
844 !, .AND. c_e(sample)) ! special case, include upper limit
845 
846 IF (present(binbounds)) binbounds = lbinbounds
847 
848 END SUBROUTINE stat_bind
849 
850 
851 FUNCTION stat_mode_histogramr(sample, nbin, start, finish, mask) RESULT(mode)
852 REAL,INTENT(in) :: sample(:)
853 INTEGER,INTENT(in) :: nbin
854 REAL,INTENT(in),OPTIONAL :: start
855 REAL,INTENT(in),OPTIONAL :: finish
856 LOGICAL,INTENT(in),OPTIONAL :: mask(:)
857 
858 REAL :: mode
859 
860 INTEGER :: loc(1)
861 INTEGER,ALLOCATABLE :: bin(:)
862 REAL,ALLOCATABLE :: binbounds(:)
863 
864 CALL stat_bin(sample, bin, nbin, start, finish, mask, binbounds)
865 mode = rmiss
866 IF (ALLOCATED(bin)) THEN
867  loc = maxloc(bin)
868  mode = (binbounds(loc(1)) + binbounds(loc(1)+1))*0.5
869 ENDIF
870 
871 END FUNCTION stat_mode_histogramr
872 
873 
874 FUNCTION stat_mode_histogramd(sample, nbin, start, finish, mask) RESULT(mode)
875 DOUBLE PRECISION,INTENT(in) :: sample(:)
876 INTEGER,INTENT(in) :: nbin
877 DOUBLE PRECISION,INTENT(in),OPTIONAL :: start
878 DOUBLE PRECISION,INTENT(in),OPTIONAL :: finish
879 LOGICAL,INTENT(in),OPTIONAL :: mask(:)
880 
881 DOUBLE PRECISION :: mode
882 
883 INTEGER :: loc(1)
884 INTEGER,ALLOCATABLE :: bin(:)
885 DOUBLE PRECISION,ALLOCATABLE :: binbounds(:)
886 
887 CALL stat_bin(sample, bin, nbin, start, finish, mask, binbounds)
888 mode = rmiss
889 IF (ALLOCATED(bin)) THEN
890  loc = maxloc(bin)
891  mode = (binbounds(loc(1)) + binbounds(loc(1)+1))*0.5d0
892 ENDIF
893 
894 END FUNCTION stat_mode_histogramd
895 
896 !!$ Calcolo degli NDI calcolati
897 !!$
898 !!$ Gli NDI vengono calcolati all’interno degli intervalli dei percentili.
899 !!$ Questo significa che il numero di NDI è pari al numero di percentili
900 !!$ non degeneri meno 1; il metodo di calcolo è il seguente:
901 !!$
902 !!$ si calcola il Density Index (DI) per ogni intervallo di percentili
903 !!$ come rapporto fra quanti intervalli di percentile vengono utilizzati
904 !!$ per definire la larghezza dell’intervallo, più di uno se si hanno
905 !!$ percentili degeneri, e la larghezza dell’intervallo stesso.
906 !!$
907 !!$ Si calcola il valore dell’ADI (Actual Density Index) dividendo i DI
908 !!$ ottenuti in intervalli aventi percentili degeneri per il numero di
909 !!$ intervalli che hanno contribuito a definire quell’intervallo.
910 !!$
911 !!$ si calcola il valore del 50-esimo percentile dei valori degli ADI
912 !!$ precedentemente ottenuti;
913 !!$
914 !!$ si divide il valore degli ADI per il 50-percentile degli ADI in
915 !!$ maniera tale da normalizzarli ottenendo così gli NADI (Normalized
916 !!$ Actual Density Index) calcolati
917 !!$
918 !!$ Si sommano i NADI relativi ad intervalli di percentili degeneri
919 !!$ ottenendo gli NDI
920 
921 
922 
923 !!$! Sobroutine calculating the number of NDI values
924 !!$SUBROUTINE NOFNDI_old(npclint, pcl0, pclint, pcl100, nndi)
925 !!$
926 !!$INTEGER, INTENT(IN) :: npclint !< number of pclint
927 !!$REAL, INTENT(IN) :: pcl0 !< the 0-th percentile = min(sample)
928 !!$REAL, DIMENSION(npclint), INTENT(IN) :: pclint !< the percentiles between the 0-th percentile and the 100-th percentile
929 !!$REAL, INTENT(IN) :: pcl100 !< the 100-th percentile = max(sample)
930 !!$INTEGER, INTENT(OUT) :: nndi !< number of NDI
931 !!$
932 !!$REAL, DIMENSION(:), allocatable :: pcl !< array inglobing in itself pcl0, pclint and pcl100
933 !!$INTEGER :: &
934 !!$ npcl,&
935 !!$ ipclint,&! integer loop variable on npclint
936 !!$ ipcl, & ! integer loop variable on npcl=2+npclint
937 !!$ npcl_act ! number of non redundant values in pcl equal to the number of
938 !!$ ! pcl_act values
939 !!$
940 !!$! Calculation of npcl
941 !!$npcl=SIZE(pclint)+2
942 !!$
943 !!$! Allocation of pcl and its initialisation
944 !!$ALLOCATE(pcl(npcl))
945 !!$pcl(1)=pcl0
946 !!$DO ipclint=1,npclint
947 !!$ pcl(ipclint+1)=pclint(ipclint)
948 !!$ENDDO
949 !!$pcl(SIZE(pclint)+2)=pcl100
950 !!$
951 !!$IF ( ANY(pcl == rmiss) ) THEN
952 !!$
953 !!$ nndi=0
954 !!$
955 !!$ELSE
956 !!$
957 !!$ ! Calculation of npcl_act
958 !!$ npcl_act=1
959 !!$ DO ipcl=2,npcl
960 !!$ IF (pcl(ipcl) /= pcl(ipcl-1)) THEN
961 !!$ npcl_act=npcl_act+1
962 !!$ ENDIF
963 !!$ ENDDO
964 !!$
965 !!$ ! Definition of number of ndi that is the wanted output value
966 !!$ nndi=npcl_act-1
967 !!$
968 !!$ENDIF
969 !!$
970 !!$
971 !!$END SUBROUTINE NOFNDI_old
972 !!$
973 !!$
974 !!$
975 !!$! Sobroutine calculating the ndi values
976 !!$SUBROUTINE NDIC_old(npclint, pcl0, pclint, pcl100, nndi, ndi, limbins)
977 !!$
978 !!$INTEGER, INTENT(IN) :: npclint !< number of pclint
979 !!$REAL, INTENT(IN) :: pcl0 !< the 0-th percentile = min(sample)
980 !!$REAL, DIMENSION(npclint), INTENT(IN) :: pclint !< the percentiles between the 0-th percentile and the 100-th percentile
981 !!$REAL, INTENT(IN) :: pcl100 !< the 100-th percentile = max(sample)
982 !!$INTEGER, INTENT(IN) :: nndi !< number of ndi
983 !!$REAL, DIMENSION(nndi), INTENT(OUT) :: ndi !< ndi that I want to calculate
984 !!$REAL, DIMENSION(nndi+1), INTENT(OUT) :: limbins !< ndi that I want to calculate
985 !!$
986 !!$REAL, DIMENSION(:), allocatable :: &
987 !!$ pcl, & ! array inglobing in itself pcl0, pclint and pcl100
988 !!$ pcl_act, & ! actual values of pcls removing redundat information
989 !!$ ranges, & ! distances between the different pcl_act(s)
990 !!$ di, & ! density indexes
991 !!$ adi, & ! actual density indexes
992 !!$ nadi ! normalized actual density indexes
993 !!$
994 !!$INTEGER, DIMENSION(:), allocatable :: weights ! weights = number of redundance times of a certain values
995 !!$ ! in pcl values
996 !!$
997 !!$REAL, DIMENSION(1) :: &
998 !!$ perc_vals, & ! perc_vals contains the percentile position (0-100)
999 !!$ med ! mediana value
1000 !!$
1001 !!$REAL, DIMENSION(:), allocatable :: infopclval
1002 !!$INTEGER, DIMENSION(:), allocatable :: infopclnum
1003 !!$
1004 !!$INTEGER :: &
1005 !!$ nbins,& ! number of intervals
1006 !!$ ibins,& ! integer loop variable on nbins_act
1007 !!$ nbins_act,& ! number of non redundant intervals
1008 !!$ ibins_act,& ! integer loop variable on nbins_act
1009 !!$ ipclint,& ! integer loop variable on npclint
1010 !!$ npcl, & ! number of percentiles
1011 !!$ ipcl, & ! integer loop variable on npcl
1012 !!$ npcl_act,& ! number of non redundant percentiles
1013 !!$ ipcl_act,& ! integer loop variable on npcl_act
1014 !!$ npclplus, & ! plus number information
1015 !!$ ncount ! number of redundant percentiles values
1016 !!$
1017 !!$
1018 !!$! Actual number of percentiles
1019 !!$npcl=SIZE(pclint)+2
1020 !!$
1021 !!$! Allocation of infopclval and infopclnum
1022 !!$ALLOCATE(infopclval(npcl))
1023 !!$ALLOCATE(infopclnum(npcl))
1024 !!$infopclval(:)=0
1025 !!$infopclnum(:)=0
1026 !!$
1027 !!$! Allocation of pcl
1028 !!$ALLOCATE(pcl(npcl))
1029 !!$! and storing of pcl values
1030 !!$pcl(1)=pcl0
1031 !!$DO ipclint=1,npclint
1032 !!$ pcl(ipclint+1)=pclint(ipclint)
1033 !!$ENDDO
1034 !!$pcl(SIZE(pclint)+2)=pcl100
1035 !!$
1036 !!$ ! Calculation of non redundant values of percentiles
1037 !!$
1038 !!$npcl_act=1
1039 !!$infopclval(1) = pcl(1)
1040 !!$infopclnum(1) = 1
1041 !!$
1042 !!$DO ipcl=2,npcl
1043 !!$ infopclval(ipcl) = pcl(ipcl)
1044 !!$ IF ( pcl(ipcl) /= pcl(ipcl-1) ) THEN
1045 !!$ npcl_act = npcl_act + 1
1046 !!$ ENDIF
1047 !!$ infopclnum(ipcl) = npcl_act
1048 !!$ENDDO
1049 !!$
1050 !!$ ! Allocation of pcl_act
1051 !!$ALLOCATE(pcl_act(npcl_act))
1052 !!$ ! and storing in pcl_act of percentiles values
1053 !!$DO ipcl_act=1,npcl_act
1054 !!$ DO ipcl=1,npcl
1055 !!$ IF (ipcl_act == infopclnum(ipcl)) THEN
1056 !!$ pcl_act(ipcl_act) = pcl(ipcl)
1057 !!$ CYCLE
1058 !!$ ENDIF
1059 !!$ ENDDO
1060 !!$ENDDO
1061 !!$
1062 !!$
1063 !!$ ! Allocation of ranges, weights and di and their initialisation
1064 !!$ALLOCATE(ranges(npcl_act-1))
1065 !!$ALLOCATE(weights(npcl_act-1))
1066 !!$ALLOCATE(di(npcl_act-1))
1067 !!$ranges(:)=0
1068 !!$di(:)=0
1069 !!$weights(:)=0
1070 !!$
1071 !!$ ! Definition of nbins_act
1072 !!$nbins_act=npcl_act-1
1073 !!$ ! Cycle on ibins_act for calculating ranges and weights
1074 !!$ ! and consequently the di values
1075 !!$DO ibins_act=1,nbins_act
1076 !!$ ranges(ibins_act)=pcl_act(ibins_act+1) - pcl_act(ibins_act)
1077 !!$
1078 !!$ IF ( pcl_act(ibins_act+1) == pcl_act(npcl_act) ) THEN
1079 !!$ weights(ibins_act) = COUNT( ibins_act == infopclnum ) + &
1080 !!$ COUNT( ibins_act+1 == infopclnum ) - 1
1081 !!$ ELSE
1082 !!$ weights(ibins_act) = COUNT ( ibins_act == infopclnum )
1083 !!$ ENDIF
1084 !!$ di(ibins_act) = weights(ibins_act)/ranges(ibins_act)
1085 !!$ENDDO
1086 !!$
1087 !!$ ! Allocation of adi and its initialisation
1088 !!$ALLOCATE(adi(npcl-1))
1089 !!$adi(:)=0
1090 !!$npclplus=0
1091 !!$DO ibins_act=1,nbins_act
1092 !!$ ncount=weights(ibins_act)
1093 !!$ ! Calculation of adi
1094 !!$ DO ibins=npclplus + 1,npclplus + ncount
1095 !!$ adi(ibins)=di(ibins_act)/ncount
1096 !!$ ENDDO
1097 !!$ npclplus=npclplus+ncount
1098 !!$ENDDO
1099 !!$
1100 !!$ ! Mediana calculation for perc_vals
1101 !!$perc_vals(1)=50
1102 !!$med = stat_percentile(sample=adi, perc_vals=perc_vals)
1103 !!$ ! Allocation of nadi and its initialisation
1104 !!$ALLOCATE(nadi(npcl-1))
1105 !!$nadi(:)=0
1106 !!$ ! Definition of values of nadi
1107 !!$nadi(:) = adi(:)/med(1)
1108 !!$
1109 !!$ ! Initialisation of ndi
1110 !!$ndi(:)=0
1111 !!$ ! Calculation of the ndi values
1112 !!$ipcl_act=1
1113 !!$nbins=npcl-1
1114 !!$DO ibins=1,nbins
1115 !!$ ndi(ipcl_act)=nadi(ibins)+ndi(ipcl_act)
1116 !!$ IF ( ( pcl(ibins+1) /= pcl(ibins) ) .AND. &
1117 !!$ ( pcl(ibins+1) /= pcl(npcl) ) ) THEN
1118 !!$ ipcl_act=ipcl_act+1
1119 !!$ ENDIF
1120 !!$ENDDO
1121 !!$
1122 !!$DO ipcl_act=1,npcl_act
1123 !!$ limbins(ipcl_act)=pcl_act(ipcl_act)
1124 !!$ENDDO
1125 !!$
1126 !!$ ! Deallocation part
1127 !!$DEALLOCATE(infopclval)
1128 !!$DEALLOCATE(infopclnum)
1129 !!$DEALLOCATE(pcl)
1130 !!$DEALLOCATE(pcl_act)
1131 !!$DEALLOCATE(ranges)
1132 !!$DEALLOCATE(weights)
1133 !!$DEALLOCATE(di)
1134 !!$DEALLOCATE(adi)
1135 !!$DEALLOCATE(nadi)
1136 !!$
1137 !!$END SUBROUTINE NDIC_old
1138 !!$
1139 !!$
1140 !!$!> Calculate the number of NDI values
1141 !!$SUBROUTINE NOFNDI(pcl, nndi)
1142 !!$
1143 !!$REAL, DIMENSION(:), INTENT(IN) :: pcl !< the percentiles between the 0-th percentile and the 100-th percentile
1144 !!$INTEGER, INTENT(OUT) :: nndi !< number of NDI
1145 !!$
1146 !!$IF ( ANY(pcl == rmiss) ) then
1147 !!$ nndi=0
1148 !!$else
1149 !!$ nndi=count_distinct(pcl)-1
1150 !!$end IF
1151 !!$
1152 !!$END SUBROUTINE NOFNDI
1153 
1154 
1155 
1156 !!$!example to manage exceptions
1157 !!$
1158 !!$use,intrinsic :: IEEE_EXCEPTIONS
1159 !!$
1160 !!$logical fail_o,fail_zero
1161 !!$
1162 !!$a=10.
1163 !!$b=tiny(0.)
1164 !!$
1165 !!$call safe_divide(a,b,c,fail_o,fail_zero)
1166 !!$
1167 !!$print*,fail_o,fail_zero
1168 !!$print *,a,b,c
1169 !!$
1170 !!$b=0.
1171 !!$call safe_divide(a,b,c,fail_o,fail_zero)
1172 !!$
1173 !!$print*,fail_o,fail_zero
1174 !!$print *,a,b,c
1175 !!$
1176 !!$contains
1177 !!$
1178 !!$subroutine safe_divide(a, b, c, fail_o,fail_zero)
1179 !!$
1180 !!$real a, b, c
1181 !!$logical fail_o,fail_zero
1182 !!$type(IEEE_STATUS_TYPE) status
1183 !!$! save the current floating-point environment, turn halting for
1184 !!$! divide-by-zero off, and clear any previous divide-by-zero flag
1185 !!$call IEEE_GET_STATUS(status)
1186 !!$call IEEE_SET_HALTING_MODE(IEEE_DIVIDE_BY_ZERO, .false.)
1187 !!$call IEEE_SET_HALTING_MODE(ieee_overflow, .false.)
1188 !!$
1189 !!$call IEEE_SET_FLAG(ieee_overflow, .false.)
1190 !!$call IEEE_SET_FLAG(IEEE_DIVIDE_BY_ZERO, .false.)
1191 !!$! perform the operation
1192 !!$c = a/b
1193 !!$! determine if a failure occurred and restore the floating-point environment
1194 !!$call IEEE_GET_FLAG(ieee_overflow, fail_o)
1195 !!$call IEEE_GET_FLAG(IEEE_DIVIDE_BY_ZERO, fail_zero)
1196 !!$call IEEE_SET_STATUS(status)
1197 !!$end subroutine safe_divide
1198 !!$end program
1199 !!$
1200 
1201 
1202 ! here you have to adopt the example above in the code below the use the wrong logic isnan()
1203 
1204 !!$!check NaN
1205 !!$if (any(isnan(di)) .and. .not. all(isnan(di))) then
1206 !!$
1207 !!$! from left
1208 !!$ do i=2,size(delta)
1209 !!$ if (isnan(di(i)) .and. .not. isnan(di(i-1))) then
1210 !!$ delta(i) = delta(i-1)
1211 !!$ w(i) = w(i) + w(i-1)
1212 !!$ end if
1213 !!$ end do
1214 !!$
1215 !!$! recompute
1216 !!$ print *,"WW=",w
1217 !!$ di = w/delta
1218 !!$
1219 !!$ if (any(isnan(di)) .and. .not. all(isnan(di))) then
1220 !!$
1221 !!$! from right
1222 !!$ do i=size(delta)-1,1,-1
1223 !!$ if (isnan(di(i)) .and. .not. isnan(di(i+1))) then
1224 !!$ delta(i) = delta(i+1)
1225 !!$ w(i) = w(i) + w(i+1)
1226 !!$ end if
1227 !!$ end do
1228 !!$
1229 !!$! one more step
1230 !!$ call DensityIndex(di,perc_vals,limbins)
1231 !!$ end if
1232 !!$end if
1233 
1234 
1235 !!$subroutine DensityIndex_old(di,perc_vals,limbins)
1236 !!$real,intent(inout) :: di(:)
1237 !!$real,intent(in) :: perc_vals(:)
1238 !!$real,intent(in) :: limbins(:)
1239 !!$
1240 !!$real :: delta(size(di)),w(size(di))
1241 !!$integer :: i
1242 !!$
1243 !!$do i=1,size(delta)
1244 !!$ delta(i) = limbins(i+1) - limbins(i)
1245 !!$ w(i) = perc_vals(i+1) - perc_vals(i)
1246 !!$end do
1247 !!$
1248 !!$di=rmiss
1249 !!$
1250 !!$if ( .not. all(delta == 0.)) then
1251 !!$ call DensityIndex_recurse(delta,w)
1252 !!$ di = w/delta
1253 !!$end if
1254 !!$
1255 !!$end subroutine DensityIndex_old
1256 !!$
1257 !!$recursive subroutine DensityIndex_recurse(delta,w)
1258 !!$real :: delta(:),w(:)
1259 !!$integer :: i
1260 !!$
1261 !!$!check divide by 0
1262 !!$if (any(delta == 0.0)) then
1263 !!$
1264 !!$! from left
1265 !!$ do i=2,size(delta)
1266 !!$ if (delta(i) == 0.0 .and. delta(i-1) > 0.0 ) then
1267 !!$ delta(i) = delta(i-1)
1268 !!$ w(i) = w(i) + w(i-1)
1269 !!$ end if
1270 !!$ end do
1271 !!$
1272 !!$! from right
1273 !!$ do i=size(delta)-1,1,-1
1274 !!$ if (delta(i) == 0.0 .and. delta(i+1) > 0.0 )then
1275 !!$ delta(i) = delta(i+1)
1276 !!$ w(i) = w(i) + w(i+1)
1277 !!$ end if
1278 !!$ end do
1279 !!$
1280 !!$end if
1281 !!$
1282 !!$! one more step
1283 !!$if (any(delta == 0.0)) then
1284 !!$ call DensityIndex_recurse(delta,w)
1285 !!$end if
1286 !!$
1287 !!$end subroutine DensityIndex_recurse
1288 !!$
1289 
1290 subroutine densityindex(di,nlimbins,occu,rnum,limbins)
1291 real,intent(out) :: di(:)
1292 real,intent(out) :: nlimbins(:)
1293 integer,intent(out) :: occu(:)
1294 REAL, DIMENSION(:), INTENT(IN) :: rnum
1295 real,intent(in) :: limbins(:)
1296 
1297 real :: nnum(size(rnum))
1298 integer :: i,k,sample_count
1299 logical :: sample_mask(size(rnum))
1300 
1301 di=rmiss
1302 occu=imiss
1303 nlimbins=rmiss
1304 
1305 nlimbins(1)=limbins(1) ! compute unique limits
1306 k=1
1307 do i=2,size(limbins)
1308  if (limbins(i) /= limbins(k)) then
1309  k=k+1
1310  nlimbins(k)= limbins(i)
1311  end if
1312 end do
1313 
1314 di=rmiss
1315 if (k == 1) return
1316 
1317 sample_mask = (rnum /= rmiss) ! remove missing values
1318 sample_count = count(sample_mask)
1319 IF (sample_count == 0) RETURN
1320 nnum(1:sample_count) = pack(rnum, mask=sample_mask)
1321 
1322 do i=1,k-2 ! compute occorrence and density index
1323  occu(i)=count(nnum>=nlimbins(i) .and. nnum<nlimbins(i+1))
1324  di(i) = float(occu(i)) / (nlimbins(i+1) - nlimbins(i))
1325 end do
1326 
1327 i=k-1 ! the last if is <=
1328 occu(i)=count(nnum>=nlimbins(i) .and. nnum<=nlimbins(i+1))
1329 di(i) = float(occu(i)) / (nlimbins(i+1) - nlimbins(i))
1330 
1331 end subroutine densityindex
1332 
1333 
1335 SUBROUTINE normalizeddensityindex(rnum, perc_vals, ndi, nlimbins)
1336 
1337 REAL, DIMENSION(:), INTENT(IN) :: rnum
1338 REAL, DIMENSION(:), INTENT(IN) :: perc_vals
1339 REAL, DIMENSION(:), INTENT(OUT) :: ndi
1340 REAL, DIMENSION(:), INTENT(OUT) :: nlimbins
1341 
1342 REAL, DIMENSION(size(ndi)) :: di
1343 INTEGER, DIMENSION(size(ndi)) :: occu
1344 REAL, DIMENSION(size(nlimbins)) :: limbins
1345 real :: med
1346 integer :: i,k,middle
1347 
1348 ndi=rmiss
1349 limbins = stat_percentile(rnum,perc_vals) ! compute percentile
1350 
1351 call densityindex(di,nlimbins,occu,rnum,limbins)
1352 
1353 ! Mediana calculation for density index
1354 k=0
1355 middle=count(c_e(rnum))/2
1356 
1357 do i=1,count(c_e(occu))
1358  k=k+occu(i)
1359  if (k > middle) then
1360  if (k > 1 .and. (k - occu(i)) == middle) then
1361  med = (di(i-1) + di(i)) / 2.
1362  else
1363  med = di(i)
1364  end if
1365  exit
1366  end if
1367 end do
1368 
1369 !weighted density index
1370 ndi(:count(c_e(di))) = min(pack(di,mask=c_e(di))/med,1.0)
1371 
1372 END SUBROUTINE normalizeddensityindex
1373 
1374 
1375 !! TODO translate from python
1376 
1377 !!$def gauss(x, A=1, mu=1, sigma=1):
1378 !!$ """
1379 !!$ Evaluate Gaussian.
1380 !!$
1381 !!$ Parameters
1382 !!$ ----------
1383 !!$ A : float
1384 !!$ Amplitude.
1385 !!$ mu : float
1386 !!$ Mean.
1387 !!$ std : float
1388 !!$ Standard deviation.
1389 !!$
1390 !!$ """
1391 !!$ return np.real(A * np.exp(-(x - mu)**2 / (2 * sigma**2)))
1392 !!$
1393 !!$def fit_direct(x, y, F=0, weighted=True, _weights=None):
1394 !!$ """Fit a Gaussian to the given data.
1395 !!$
1396 !!$ Returns a fit so that y ~ gauss(x, A, mu, sigma)
1397 !!$
1398 !!$ Parameters
1399 !!$ ----------
1400 !!$ x : ndarray
1401 !!$ Sampling positions.
1402 !!$ y : ndarray
1403 !!$ Sampled values.
1404 !!$ F : float
1405 !!$ Ignore values of y <= F.
1406 !!$ weighted : bool
1407 !!$ Whether to use weighted least squares. If True, weigh
1408 !!$ the error function by y, ensuring that small values
1409 !!$ has less influence on the outcome.
1410 !!$
1411 !!$ Additional Parameters
1412 !!$ ---------------------
1413 !!$ _weights : ndarray
1414 !!$ Weights used in weighted least squares. For internal use
1415 !!$ by fit_iterative.
1416 !!$
1417 !!$ Returns
1418 !!$ -------
1419 !!$ A : float
1420 !!$ Amplitude.
1421 !!$ mu : float
1422 !!$ Mean.
1423 !!$ std : float
1424 !!$ Standard deviation.
1425 !!$
1426 !!$ """
1427 !!$ mask = (y > F)
1428 !!$ x = x[mask]
1429 !!$ y = y[mask]
1430 !!$
1431 !!$ if _weights is None:
1432 !!$ _weights = y
1433 !!$ else:
1434 !!$ _weights = _weights[mask]
1435 !!$
1436 !!$ # We do not want to risk working with negative values
1437 !!$ np.clip(y, 1e-10, np.inf, y)
1438 !!$
1439 !!$ e = np.ones(len(x))
1440 !!$ if weighted:
1441 !!$ e = e * (_weights**2)
1442 !!$
1443 !!$ v = (np.sum(np.vander(x, 5) * e[:, None], axis=0))[::-1]
1444 !!$ A = v[sl.hankel([0, 1, 2], [2, 3, 4])]
1445 !!$
1446 !!$ ly = e * np.log(y)
1447 !!$ ls = np.sum(ly)
1448 !!$ x_ls = np.sum(ly * x)
1449 !!$ xx_ls = np.sum(ly * x**2)
1450 !!$ B = np.array([ls, x_ls, xx_ls])
1451 !!$
1452 !!$ (a, b, c), res, rank, s = np.linalg.lstsq(A, B)
1453 !!$
1454 !!$ A = np.exp(a - (b**2 / (4 * c)))
1455 !!$ mu = -b / (2 * c)
1456 !!$ sigma = sp.sqrt(-1 / (2 * c))
1457 !!$
1458 !!$ return A, mu, sigma
1459 !!$
1460 !!$def fit_iterative(x, y, F=0, weighted=True, N=10):
1461 !!$ """Fit a Gaussian to the given data.
1462 !!$
1463 !!$ Returns a fit so that y ~ gauss(x, A, mu, sigma)
1464 !!$
1465 !!$ This function iteratively fits using fit_direct.
1466 !!$
1467 !!$ Parameters
1468 !!$ ----------
1469 !!$ x : ndarray
1470 !!$ Sampling positions.
1471 !!$ y : ndarray
1472 !!$ Sampled values.
1473 !!$ F : float
1474 !!$ Ignore values of y <= F.
1475 !!$ weighted : bool
1476 !!$ Whether to use weighted least squares. If True, weigh
1477 !!$ the error function by y, ensuring that small values
1478 !!$ has less influence on the outcome.
1479 !!$ N : int
1480 !!$ Number of iterations.
1481 !!$
1482 !!$ Returns
1483 !!$ -------
1484 !!$ A : float
1485 !!$ Amplitude.
1486 !!$ mu : float
1487 !!$ Mean.
1488 !!$ std : float
1489 !!$ Standard deviation.
1490 !!$
1491 !!$ """
1492 !!$ y_ = y
1493 !!$ for i in range(N):
1494 !!$ p = fit_direct(x, y, weighted=True, _weights=y_)
1495 !!$ A, mu, sigma = p
1496 !!$ y_ = gauss(x, A, mu, sigma)
1497 !!$
1498 !!$ return np.real(A), np.real(mu), np.real(sigma)
1499 !!$
1500 
1501 
1502 
1503 END MODULE simple_stat
Compute the average of the random variable provided, taking into account missing data.
Definition: simple_stat.f90:39
Definitions of constants and functions for working with missing values.
Bin a sample into equally spaced intervals to form a histogram.
Compute the mode of the random variable provided taking into account missing data.
Compute a set of percentiles for a random variable.
Module for basic statistical computations taking into account missing data.
Definition: simple_stat.f90:25
Compute the standard deviation of the random variable provided, taking into account missing data...
Definition: simple_stat.f90:69
Compute the linear regression coefficients between the two random variables provided, taking into account missing data.
Compute the variance of the random variable provided, taking into account missing data...
Definition: simple_stat.f90:54
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
This module defines usefull general purpose function and subroutine.
Compute the linear correlation coefficient between the two random variables provided, taking into account missing data.
Definition: simple_stat.f90:89

Generated with Doxygen.