!! compute statistics
!|author:  <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a>
! license: <a href="http://www.gnu.org/licenses/">GPL</a>
!    
!### History
!
! current version 1.1 - 19th May 2021 
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 21/Jun/2017 | Original code |
! | 1.1      | 19/May/2021 |linearRegressionSlope compute R2 as optional argument |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
! This file is part of 
!
!   MOSAICO -- MOdular library for raSter bAsed hydrologIcal appliCatiOn.
! 
!   Copyright (C) 2011 Giovanni Ravazzani
!
!### Module Description 
!  library to compute statistics
!
MODULE Statistics         
			
! Modules used: 

USE DataTypeSizes, ONLY: &
  ! Imported type definitions:
  short, long, float 

USE LogLib, ONLY: &
  ! Imported routines:
  Catch
 

IMPLICIT NONE 

! Global (Public) 
INTERFACE GetMeanVector
   MODULE PROCEDURE GetMeanFloat
END INTERFACE

PUBLIC :: GetMeanVector
PUBLIC :: LinearRegressionSlope

! Local (i.e. private) 

PRIVATE :: GetMeanFloat


!=======         
CONTAINS
!======= 
! Define procedures contained in this module. 


!==============================================================================
!| Description:
!  compute mean of a vector of real numbers. If nodata is passed
!   numbers are filetered before computing mean
FUNCTION GetMeanFloat &
!
(vector, nodata) &
!
RESULT (mean)

IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT(IN) :: vector (:)
REAL (KIND = float), OPTIONAL, INTENT(IN) :: nodata


!Local declarations:
REAL (KIND = float) :: mean
REAL (KIND = float) :: count
INTEGER (KIND = short) :: i

!---------------------------end of declarations--------------------------------
mean = 0.
count = 0.

!cumulate
DO i = 1, SIZE(vector)
  IF ( PRESENT(nodata) ) THEN
    IF (vector(i) /= nodata) THEN
       mean = mean + vector (i)
       count = count + 1.
    END IF
  ELSE
    mean = mean + vector (i)
    count = count + 1.
  END IF
END DO

!compute mean
IF ( PRESENT(nodata) ) THEN 
  IF (count == 0.) THEN
     mean = 0.
     CALL Catch ('warning', 'Statistics',  &
        'calculate mean: no valid numbers in vector: ', argument = &
        'mean set to zero' )
  ELSE
     mean = mean / count
  END IF
ELSE
   mean = mean / count
END IF


RETURN
END FUNCTION GetMeanFloat


!==============================================================================
!| Description:
!   compute slope of linear regression between vector x and y
FUNCTION LinearRegressionSlope &
!
(x, y, nodata, r2) &
!
RESULT (lrs)

IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT(IN) :: x (:)
REAL (KIND = float), INTENT(IN) :: y (:)
REAL (KIND = float), OPTIONAL, INTENT(IN) :: nodata

!Arguments with intent(out)
REAL (KIND = float), OPTIONAL, INTENT(OUT) :: r2


!Local declarations:
REAL (KIND = float) :: lrs !!slope of linear regression
REAL (KIND = float) :: meanx, meany !!mean of x and y vectors
REAL (KIND = float), ALLOCATABLE :: vx(:), vy(:) !!vectors containing  pairs of valid data
REAL (KIND = float) :: num, den !!numerator and denominator for computing slope
REAL (KIND = float) :: sumx, sumy, sumxy, sumx2, sumy2 !!used to compute r2
INTEGER (KIND = short) :: i, j, count

!---------------------------end of declarations--------------------------------

!check x and y have the same size
IF ( SIZE(x) /= SIZE(y) ) THEN
  CALL Catch ('error', 'Statistics',  &
        'calculate linear regression slope: ', argument = &
        'vectors have different size' ) 
END IF

!remove nodata, keep only pairs of valid data
!count pairs of valid data
IF ( PRESENT (nodata) )  THEN
    count = 0
    DO i = 1, SIZE (x)
      IF ( x (i) /= nodata .AND. &
           y (i) /= nodata ) THEN !found one pair of valid data
           count = count + 1
      END IF
    END DO

    IF (count == 0.) THEN
         lrs = 0.
         CALL Catch ('warning', 'Statistics',  &
            'calculate slope of linear regression: no valid numbers in vectors, ', &
              argument = 'slope set to zero' )
         RETURN
    END IF
END IF

!allocate memory of temporary files and store valid pairs
IF ( PRESENT (nodata) .AND. count > 0 ) THEN
    ALLOCATE ( vx(count), vy(count) )
    j = 1
    DO i = 1, SIZE (x)
      IF ( x (i) /= nodata .AND. &
           y (i) /= nodata ) THEN !found one pair of valid data
           vx(j) = x(i)
           vy(j) = y(i)
           j = j + 1
      END IF
    END DO
END IF

!compute mean of x and y vectors
IF ( PRESENT (nodata) .AND. count > 0 ) THEN
   meanx = GetMeanVector (vx)
   meany = GetMeanVector (vy)
ELSE
   meanx = GetMeanVector (x)
   meany = GetMeanVector (y)
END IF

!cumulate
num = 0.
den = 0.

IF ( PRESENT (nodata) .AND. count > 0 ) THEN
    DO i = 1, count
      num = num + ( vx(i) - meanx ) * ( vy(i) - meany )
      den = den + ( vx(i) - meanx ) ** 2.
    END DO

ELSE
    DO i = 1, SIZE(x)
      num = num + ( x(i) - meanx ) * ( y(i) - meany )
      den = den + ( x(i) - meanx ) ** 2.
    END DO

END IF

!compute slope
IF ( den == 0. ) THEN
  CALL Catch ('warning', 'Statistics',  &
  'calculate slope of linear regression: slope is infinite, ', argument = &
  'slope set to huge number' )
  lrs = HUGE (lrs)
ELSE
  lrs = num / den
END IF

!optional: compute R2 (the square of the correlation coefficient)

IF (PRESENT (r2) ) THEN
    count = SIZE (vx)
    sumx = 0.
    sumy = 0.
    sumxy = 0.
    sumx2 = 0.
    sumy2 = 0.
    DO i = 1, count
        sumx = sumx + vx (i)
        sumy = sumy + vy (i)
        sumxy = sumxy + vx (i) * vy (i)
        sumx2 = sumx2 + ( vx (i) )**2.
        sumy2 = sumy2 + ( vy (i) )**2.
    END DO
   
    r2 = (  count * sumxy - sumx * sumy ) / &
         (  (count * sumx2 - sumx**2. )  * (count * sumy2 - sumy**2. )    )**0.5
    r2 = r2 * r2
    
END IF

IF ( PRESENT (nodata) .AND. count > 0 ) THEN
  DEALLOCATE (vx, vy)
END IF


RETURN
END FUNCTION LinearRegressionSlope



END MODULE Statistics