*-----------------------------------------------------------------------
* average.f - routines for ordinary, cyclic, and weighted averages
* by Andy Allinger, 2013-2017, released to the public domain
* This file may be used by any person for any purpose.
*-----------------------------------------------------------------------

*-----------------------------------------------------------------------
* MEAN - arithmetic mean
*   X[N]       Real      In        data array
*   N          Integer   In        length of X
*   AVE        Real      Out       average
*-----------------------------------------------------------------------
      SUBROUTINE MEAN (X, N, AVE)
       IMPLICIT NONE
       INTEGER N
       REAL X, AVE
       DIMENSION X(N)
       
       IF (N < 1) RETURN
       CALL ADD (X, N, AVE)
       AVE = AVE / FLOAT(N)
       RETURN
      END  ! of MEAN
      
 
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* CYMEAN - cyclic mean
*   by Andy Allinger, 2013, released to the public domain
* origin:  David Radcliffe, 2005 (or earlier)
*
*      He proposed the algorithm:
* 1. Find the average and standard deviation of the given numbers.
* 2. Increase the smallest number by 360.
* 3. Repeat steps 1 and 2 until all numbers are greater than 360.
* 4. Choose the average that yields the smallest standard deviation.
* 5. If the average is greater than 360 then subtract 360.
*
* Notice that unlike the linear mean, the cyclic mean does not
* always have a distinct value.  Instead, the array ave[] contains in
* its beginning entries the various averages, and the integer ties denotes
* how many averages were found.  
* WARNING:  X is returned shifted and sorted
*
*     X [n]       Real     In     data array
*     n           Integer  In     length of X
*     t           Real     In     period 
*     AVE [n]     Real     Out    average
*     ties        Integer  Out    how many values in ave[]
*     dev         Real     Out    standard deviation
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE CYMEAN (X, n, t, AVE, ties, dev)
       IMPLICIT NONE
       INTEGER n, ties
       REAL X, t, AVE, dev
       DIMENSION X(n), AVE(n)
       
       INTEGER i
       REAL
     &      rmachine,       ! function for EPSILON()
     &      r,              ! temporary scalar X[i]
     &      tol             ! a small number
       DOUBLE PRECISION 
     &                  nn,        ! double precision of n
     &                  s, s2,     ! sum, sum of squares
     &                  v, vmin,   ! n @ variance, minimal n @ variance
     &                  tt, t2,    ! period, period squared
     &                  rr         ! double precision of r
       
*        begin
       IF (n .EQ. 1) THEN 
         ties = 1
         AVE(1) = X(1)
         dev = 0.
         RETURN
       END IF
 
       nn = DBLE(n) 
       tt = DBLE(t) 
!       tol = SQRT(t * rmachine(1))
       tol = t * rmachine(1)
       t2 = tt **2
       s = 0. 
       s2 = 0. 
              
*        move all data to the range 0...t, sort
       DO i = 1, n
         r = X(i)
         r = MOD(r, t)
         IF (r < 0.) r = r + t
         rr = DBLE(r)
         s = s + rr 
         s2 = s2 + rr **2        ! accumulate sums
         X(i) = r
       END DO
       CALL SSORT (X, n)
       
*           find phase of seam of lowest deviation
       vmin = rmachine(3) 
       ties = 0
       DO i = 1, n
         v = s2 - s **2 / nn              ! n @ variance
         IF (vmin - v > tol) THEN          ! new low deviation
           ties = 0
           vmin = v
         END IF
         IF (ABS(v - vmin) .LE. tol) THEN  ! tie for new low deviation
           ties = ties + 1
           AVE(ties) = SNGL(s / nn)
         END IF
         
*           apply update formula to sum and sum of squares
         s = s + tt
         s2 = s2 + 2D0 * DBLE(X(i)) * tt + t2
       END DO
       
*           subtract if out of 0...t
       DO i = 1, ties
         IF (AVE(i) .GE. t) AVE(i) = AVE(i) - t      
       END DO
       dev = SNGL(SQRT(vmin / nn))                ! (n@variance / n)^(1/2)

       RETURN
      END ! of CYMEAN


*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* CYMED - cyclic median, a variation of Radcliffe's algorithm
*   by Andy Allinger, 2013, released to the public domain
* The cyclic median does not always have a distinct value.
*   WARNING:  X is returned shifted and sorted
*
*     X [n]       Real     In     data array
*     n           Integer  In     length of X
*     t           Real     In     period
*     AVE [n]     Real     Out    average
*     ties        Integer  Out    how many values in AVE[]
*     dev         Real     Out    deviation
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE CYMED (X, n, t, AVE, ties, dev)
       IMPLICIT NONE
       INTEGER n, ties
       REAL X, t, AVE, dev
       DIMENSION X(n), AVE(n)
       
       INTEGER i,         ! count
     &         j0, j1     ! positions where to find current estimate of median
       REAL
     &      rmachine,       ! function when EPSILON() not available
     &      r,              ! temporary copy of X
     &      tol             ! a small number
       DOUBLE PRECISION 
     &                  s1, smin,  ! sum of deviations, minimum sum deviation
     &                  tt         ! period
       LOGICAL odd              ! even or odd n
       
*        begin
       IF (n .EQ. 1) THEN
         ties = 1
         AVE(1) = X(1)
         dev = 0.         
         RETURN
       END IF
       j1 = n / 2 + 1
       odd = MOD(n, 2) .EQ. 1              ! is n an odd number?
       tt = DBLE(t)
       tol = t * rmachine(1)
              
*        move all data to the range 0...t, sort
       DO i = 1, n
         r = MOD(X(i), t)
         IF (r < 0.) r = r + t
         X(i) = r
       END DO
       CALL SSORT (X, n)

*            initial value of cyclic median
       IF (odd) THEN
         r = X(j1)
       ELSE
         j0 = j1 - 1
         r = (X(j0) + X(j1)) / 2.     ! split the difference
       END IF

*         sum of absolute deviations
       s1 = 0. 
       DO i = 1, n
         s1 = s1 + DBLE(ABS(X(i) - r))
       END DO
       
*           find phase of seam of lowest deviation
       smin = rmachine(3)
       ties = 0
       DO i = 1, n
         IF (smin - s1 > tol) THEN          ! new low deviation
           ties = 0
           smin = s1
         END IF
         IF (ABS(s1 - smin) .LE. tol) THEN  ! tie for new low deviation
           ties = ties + 1
           IF (odd) THEN
             AVE(ties) = X(j1)
           ELSE              ! two points required for even case 
             j0 = j1 - 1 
             IF (j0 .EQ. 0) j0 = n   ! j0 is the prior point 
             AVE(ties) = (X(j0) + X(j1)) / 2. 
           END IF 
         END IF
         
*           apply update formula        
         r = X(i)
         IF (odd) THEN
           j0 = j1 + 1 
           IF (j0 > n) j0 = 1        ! new and current medians
           s1 = s1 + 2. * r + tt - DBLE(X(j0) + X(j1))   ! j0 is the new median
         ELSE
           s1 = s1 + 2. * r + tt - 2. * X(j1)
         END IF
         X(i) = r + t               ! advance seam
         j1 = j1 + 1 
         IF (j1 > n) j1 = 1         
       END DO

*         bring into 0...t
       DO i = 1, ties       
         IF (AVE(i) .GE. t) AVE(i) = AVE(i) - t
       END DO       
       dev = SNGL(smin) / FLOAT(n)                     ! deviation

       RETURN
      END ! of CYMED



*----------------------------------------------------------------------- 
* imode - the mode of an array of integers 
* 
*___Name_______Type______In/Out____Description__________________________
*   JX(N)      Integer   In        An array 
*   JAVE(N)    Integer   Out       Most frequent value in JX
*   JPERM(N)   Integer   Neither   Permutation vector
*   N          Integer   In        Length of array 
*   JTIES      Integer   Out       How many ties for mode
*----------------------------------------------------------------------- 
      SUBROUTINE IMODE (JX, JAVE, JPERM, N, JTIES)      
       IMPLICIT NONE
       INTEGER JX, JAVE, JPERM, N, JTIES
       DIMENSION JX(N), JAVE(N), JPERM(N)
      
       INTEGER I, J, JI, JJ, L, M
      
*        begin 
       IF (N < 1) RETURN
       IF (N .EQ. 1) THEN                ! special case
         JAVE(1) = JX(1)
         JTIES = 1
         RETURN
       END IF            
      
*        find duplicates, mark positions and populations
       CALL QSORTI (JX, JPERM, N)
       M = 0
       J = 1
       JJ = JPERM(1)
       JI = JPERM(2)
       DO I = 2, N
         JI = JPERM(I)
         IF (JX(JI) .NE. JX(JJ)) THEN
           L = I - J
           IF (L > M) THEN
             M = L
             JTIES = 0
           END IF
           IF (L .EQ. M) THEN
             JTIES = JTIES + 1
             JAVE(JTIES) = JX(JJ)
           END IF
           J = I
           JJ = JI
         END IF
       END DO
       L = N - J + 1
       IF (L > M) THEN
         M = L
         JTIES = 0
       END IF
       IF (L .EQ. M) THEN
         JTIES = JTIES + 1
         JAVE(JTIES) = JX(JI)
       END IF

       RETURN
      END  ! of IMODE


************************************************************************
*  discrete mode of Real approximate integers
*
*                              I/O list
*_Name__________Type_________I/O________Description______________________
* X[N]          Real         In         data array
* N             Integer      In         length of array
* R             Integer      In         maximum value in X
* AVE[R]        Real         Out        modes
* TIES          Integer      Out        how many entries in AVE[]
************************************************************************
      SUBROUTINE RMODE (X, N, R, AVE, TIES)
       IMPLICIT NONE
       INTEGER N, R, TIES
       REAL X, AVE
       DIMENSION X(N), AVE(R)
       
*         local variables
       INTEGER I, J   
       
*           begin
       IF (N .EQ. 1) THEN
         TIES = 1
         AVE(1) = X(1)
         RETURN
       END IF
       
       DO J = 1, R
         AVE(J) = 0.
       END DO
       DO I = 1, N                       ! count
         J = NINT(X(I))
         AVE(J) = AVE(J) + 1.
       END DO
       
       I = 0
       TIES = 0
       DO J = 1, R                         ! find max
         IF (NINT(AVE(J)) > I) THEN
           TIES = 0
           I = NINT(AVE(J))
         END IF
         IF (NINT(AVE(J)) .EQ. I) THEN
           TIES = TIES + 1
           AVE(TIES) = FLOAT(J)
         END IF
       END DO
       
       RETURN 
      END ! of RMODE

      

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Weighted averages and weighted circular averages
!  by Andy Allinger, 2016, released to the public domain
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!            weighted mean
!__Name____Type_____In/Out________Description___________________________
!  X[N]    Real     In            data array
!  W[N]    Real     In            importance weights
!  N       Integer  In            length of arrays
!  XWM     Real     Out           weighted mean
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE WMEAN (X, W, N, XWM)
       IMPLICIT NONE
       INTEGER N
       REAL X, W, XWM
       DIMENSION X(N), W(N)
       
!          local variables
       INTEGER I
       DOUBLE PRECISION SW, SWX, C, Y, T
       
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!              begin
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       IF (N < 1) RETURN

       SW = 0.
       C = 0.
       DO I = 1, N              ! compensated addition for SW
         Y = W(I) - C
         T = SW + Y
         C = (T - SW) - Y
         SW = T
       END DO
       
       SWX = 0.
       C = 0.
       DO I = 1, N              ! compensated addition for SWX
         Y = W(I) * X(I) - C
         T = SWX + Y
         C = (T - SWX) - Y
         SWX = T
       END DO
       
       IF (SW > 0.) XWM = SNGL(SWX / SW)
       RETURN
       
      END ! of WMEAN
       

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!              weighted median
!__Name____Type_____In/Out________Description___________________________
!  X[N]    Real     In            data array
!  W[N]    Real     In            importance weights
!  IND[N]  Integer  Work          index array 
!  N       Integer  In            length of arrays
!  XWM     Real     Out           weighted median
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE WMED (X, W, IND, N, XWM)
       IMPLICIT NONE
       INTEGER IND, N
       REAL X, W, XWM
       DIMENSION X(N), W(N), IND(N)
       
!         local variables
       INTEGER I
       REAL SW, HSW
       
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!            begin
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       IF (N < 1) RETURN
       IF (N .EQ. 1) THEN
         XWM = X(1)
         RETURN
       END IF
       
       SW = 0.
       DO I = 1, N
         SW = SW + W(I)
       END DO       
       IF (SW .LE. 0.) RETURN
       HSW = SW / 2.       
       CALL QSORTR (X, IND, N)
       CALL RORDER (W, IND, N)
       CALL RORDER (X, IND, N)
       
       SW = 0.
       DO I = 1, N
         SW = SW + W(I)
         IF (SW > HSW) THEN
           XWM = X(I)
           RETURN
         ELSE IF (SW .EQ. HSW) THEN
           XWM = 0.5 * (X(I) + X(I+1))
           RETURN
         END IF
       END DO
       
      END ! of WMED
      

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!         weighted cyclic mean for phase only
!
!__Name____Type_____In/Out__Description_________________________________
!  X[N]    Real     In      data array
!  W[N]    Real     In      importance weights    
!  IND[N]  Integer  Work    permutation vector
!  N       Integer  In      length of arrays
!  T       Real     In      period 
!  AVE[N]  Real     Out     average
!  TIES    Integer  Out     how many values in AVE[]
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE WCYMEN (X, W, IND, N, T, AVE, TIES)
       IMPLICIT NONE
       INTEGER IND, N, TIES
       REAL X, W, T, AVE
       DIMENSION X(N), W(N), IND(N), AVE(N)
       
       INTEGER I
       REAL 
     &      RMACHINE,       ! function for EPSILON and HUGE
     &      R,              ! temporary scalar X[i]
     &      TOL             ! a small number
       DOUBLE PRECISION 
     &                  SW, SWX, SWX2,     ! sum of weights, sum of squares
     &                  V, VMIN,           ! n @ variance, minimal n @ variance
     &                  TT, T2,            ! period, period squared
     &                  RR,                ! double precision of R
     &                  WW                 ! double precision of W(I)
       
!        begin
       IF (N .EQ. 1) THEN 
         TIES = 1
         AVE(1) = X(1)
         RETURN
       END IF
 
       TOL = T * RMACHINE(1)
       TT = DBLE(T)
       T2 = TT **2
       SW = 0.
       SWX = 0.
       SWX2 = 0.
              
!        move all data to the range 0...T, sort
       DO I = 1, N
         R = MOD(X(I), T)
         IF (R < 0.) R = R + T
         WW = DBLE(W(I))
         RR = DBLE(R)
         SW = SW + WW                             ! accumulate sums
         SWX = SWX + WW * RR
         SWX2 = SWX2 + WW * RR **2
         X(I) = R
       END DO
       CALL QSORTR (X, IND, N)
       CALL RORDER (W, IND, N)
       CALL RORDER (X, IND, N)
       
!           find phase of seam of lowest weighted deviation
       VMIN = RMACHINE(3)
       TIES = 0
       DO I = 1, N
         V = SWX2 - SWX **2 / SW           ! SW @ weighted variance
         IF (VMIN - V > TOL) THEN          ! new low deviation
           TIES = 0
           VMIN = V
         END IF
         IF (ABS(V - VMIN) .LE. TOL) THEN  ! tie for new low deviation
           TIES = TIES + 1
           AVE(TIES) = SNGL(SWX / SW)    ! weighted mean
         END IF
         
!           apply update formula to sum and sum of squares
         WW = DBLE(W(I))
         RR = DBLE(X(I))
         SWX = SWX + WW * TT
         SWX2 = SWX2 + WW * (2D0 * RR * TT + T2)
       END DO
       
       DO I = 1, TIES
         R = AVE(I)
         IF (R .GE. T) R = R - T      ! subtract if out of 0...T
         AVE(I) = R
       END DO
       
       RETURN
      END ! of WCYMEN
*!!!!!!!!!!!!!!!!!!!!!! End of file average.f !!!!!!!!!!!!!!!!!!!!!!!!!!
