*-----------------------------------------------------------------------
* mullin.f - multiple line-segments approximation of a data series
* by Andy Allinger, 2012-2016, released to the public domain
* This program may be used by any person for any purpose.
*-----------------------------------------------------------------------
      
*-----------------------------------------------------------------------
*   find slope from least-square accumulations
*-----------------------------------------------------------------------
      SUBROUTINE CALC_SLOPE (i, k, n, CW, CWX, CWY, CWX2, CWXY, slope)     
       IMPLICIT NONE
       INTEGER i, k, n
       REAL slope
       DOUBLE PRECISION CW, CWX, CWY, CWX2, CWXY
       DIMENSION CW(0:n), CWX(0:n), CWY(0:n), CWX2(0:n),
     &  CWXY(0:n), slope(n)

*          local variables
       INTEGER h
       DOUBLE PRECISION det, sw, swx, swy, swx2, swxy
       LOGICAL DALMEQ
       
*          begin 
       h = i - 1
       sw = CW(k) - CW(h)              ! differences of accumulated sums
       swx = CWX(k) - CWX(h)
       swy = CWY(k) - CWY(h)
       swx2 = CWX2(k) - CWX2(h)
       swxy = CWXY(k) - CWXY(h)
       det = swx2 * sw - swx * swx
       IF (DALMEQ(det, 0D0)) THEN
         PRINT *, 'at hik=', h,i,k, 'Warning: zero determinant'
       ELSE
         slope(i) = SNGL((sw * swxy - swx * swy) / det)
         IF (ISNAN(slope(i)))
     &     PRINT *, 'slope: nan encountered!! sw: ', sw, ' swx: ', swx,
     &      'swy: ', swy, ' swxy: ', swxy, ' det: ', det
       END IF
       RETURN
      END  ! of CALC_SLOPE
      

*-----------------------------------------------------------------------
*     probability slopes are distinct, by a 1-tailed t-test
*-----------------------------------------------------------------------
      SUBROUTINE CALC_P (i, j, k, n, CM, CM2, slope, p)
       IMPLICIT NONE
       INTEGER i, j, k, n
       REAL slope, p
       DOUBLE PRECISION CM, CM2
       DIMENSION CM(n), CM2(n), slope(n)
       
*                local variables
       INTEGER n1, n2, nu       
       REAL CDFT,                ! function from STARPAC
     &      EPS, rmachine,       ! macine precision
     &      t, var1, var2        ! t-statistic, variances
       DOUBLE PRECISION sm, sm2
       LOGICAL first 
       SAVE first, EPS
       DATA first / .TRUE. /
       
*              begin
       IF (first) THEN
         first = .FALSE.
         EPS = rmachine(1)
       END IF       
       
*                  n1, n2 @ variances of proposed segments
       n1 = j - i
       sm = CM(j) - CM(i)
       sm2 = CM2(j) - CM2(i)
       var1 = SNGL(sm2 - sm **2 / DBLE(n1))   ! n[1] @ s[1]^2
       n2 = k - j
       sm = CM(k) - CM(j)
       sm2 = CM2(k) - CM2(j)
       var2 = SNGL(sm2 - sm **2 / DBLE(n2))   ! n[2] @ s[2]^2
           
*             probability segments are distinct
       IF (var1 + var2 > 0.) THEN
         nu = n1 + n2 - 2
         IF (nu < 1) PRINT *, 'calc_p: nu ', nu, ' with ijk ', i, j, k
         t = ABS(slope(j) - slope(i))  
     &        / SQRT( ((var1 + var2) / FLOAT(nu)) 
     &                      * (1. / FLOAT(n1) + 1. / FLOAT(n2)) )
         p = CDFT(t, nu)       ! (im)probability of segments being same
       ELSE
         IF (ABS(slope(j) - slope(i)) > EPS) THEN
           p = 1.
         ELSE
           p = 0.
         END IF
       END IF
       RETURN
      END  ! of CALC_P


*-----------------------------------------------------------------------
* Departition intervals of the data with similar slopes.  
* Inspired by the HarmAn algorithm of Pardo & Birmingham
*                               I/O LIST
*_Name_______Type__________In/Out_____Description_______________________
*  W[n]       REAL          in         importance weights of data points
*  X[n]       REAL          in         data independent variable
*  Y[n]       REAL          in         data dependent variable
*  scis[n]    REAL          out        X coordinate of segment beginnings
*  slope[n]   REAL          out        slope of line segments
*  n          INTEGER       in         number of data points
*  v          INTEGER       out        number of approximating segments
*  ierr       INTEGER       out        error code 0=no errors
*-----------------------------------------------------------------------
      SUBROUTINE mullin (W, X, Y, scis, slope, n, v, ierr)
      
       IMPLICIT NONE
       INTEGER n, v, ierr 
       REAL W, X, Y                     ! input
       REAL scis, slope                 ! output
       DIMENSION W(n), X(n), Y(n), scis(n), slope(n)
       
*                 constants
       REAL pcrit                       ! cutoff 1-tail t-test
       PARAMETER (pcrit = 0.999999) 
       
*                 local variables
       INTEGER g, h, i, j, k, vv              ! counters

       REAL
     &  EPS, rmachine,                 ! machine precision
     &  p_h, p_i,                      ! p-values merging prev or next
     &  pmin,                          ! least p-value
     &  r,                             ! a ratio
     &  tau,                           ! error tolerance in HFTI
     &  wnew, xold, xnew, yold, ynew   ! working values of w, x, y

*           least-squares sums and accumulations
       DOUBLE PRECISION
     &  m,                                ! slope
     &  wx                                ! temp for least squares sums

       LOGICAL done, neg
     
*              allocatable arrays
       INTEGER, DIMENSION (:), ALLOCATABLE :: prev, next, IP
       REAL, DIMENSION (:), ALLOCATABLE :: cept, peach
       DOUBLE PRECISION, DIMENSION (:), ALLOCATABLE ::
     &      CM, CM2, CW, CWX, CWY, CWX2, CWXY
       REAL, DIMENSION (:,:), ALLOCATABLE :: A
     
*-----------------------------------------------------------------------
*           begin
*-----------------------------------------------------------------------
       ierr = 0
       EPS = rmachine(1)
       
       IF (n < 2) THEN                    ! not enough data
         ierr = 3
         RETURN
       END IF
 
*-----------------------------------------------------------------------
*          initial values
*-----------------------------------------------------------------------
       ALLOCATE (prev(n), next(n), peach(n), cept(n), CM(n), CM2(n), 
     &  CW(0:n), CWX(0:n), CWY(0:n), CWX2(0:n), CWXY(0:n), STAT=ierr)
       IF (ierr .NE. 0) RETURN
         
       CW(0) = 0.; CWX(0) = 0.; CWY(0) = 0.; CWX2(0) = 0.; CWXY(0) = 0.
       xnew = -1E9

       DO j = 1, n
         xold = xnew  
         wnew = W(j) 
         xnew = X(j)        
         ynew = Y(j)

*           check that input X is increasing
         IF (xnew - xold < EPS) THEN
           ierr = 1
           PRINT *, 'x not distinct or not sorted'
           RETURN
         END IF
         
*           accumulate weighted least squares sums         
         CW(j) = CW(j-1) + wnew
         wx = DPROD(wnew, xnew)
         CWX(j) = CWX(j-1) + wx
         CWY(j) = CWY(j-1) + DPROD(wnew, ynew)
         CWX2(j) = CWX2(j-1) + wx * xnew
         CWXY(j) = CWXY(j-1) + wx * ynew
         
          prev(j) = j - 1
          next(j) = j + 1
       END DO

*           accumulate standard deviation sums 
       CM(1) = 0. ; CM2(1) = 0.
       xnew = X(1) ; ynew = Y(1)
       DO j = 2, n
         xold = xnew ; xnew = X(j)
         yold = ynew ; ynew = Y(j)
         m = (ynew - yold) / (xnew - xold)
         CM(j) = CM(j-1) + m
         CM2(j) = CM2(j-1) + m **2
       END DO

*-----------------------------------------------------------------------
*        departition the line segments into longer, least-squares lines
*-----------------------------------------------------------------------
       v = n - 1    ! all partitions

*          merge pairs to reach nonzero variance
       DO j = 2, n-1, 2
         i = prev(j)
         k = next(j)
         CALL CALC_SLOPE (i, k, n, CW, CWX, CWY, CWX2, CWXY, slope)
         IF (ISNAN(slope(I))) PRINT *, 'mullin: Nan at begin phase!!'
         prev(k) = i
         next(i) = k
         v = v - 1       ! one fewer segment
       END DO  ! next j
       
       IF (BTEST(n-1,0)) THEN           ! required for even n
         i = n-1
         j = n
         CALL CALC_SLOPE (i, j, n, CW, CWX, CWY, CWX2, CWXY, slope)
         IF (ISNAN(slope(i))) PRINT *, 'mullin: NaN on even N!!!'
       END IF
       
*----------------------------------------------------------------------- 
*           evaluate other candidates to merge
*----------------------------------------------------------------------- 
      DO j = 3, n-1, 2
         i = prev(j)
         k = next(j)
         CALL CALC_P (i, j, k, n, CM, CM2, slope, peach(j))
       END DO
       
*----------------------------------------------------------------------- 
*         merge segments one at a time
*-----------------------------------------------------------------------
       done = .FALSE.
       DO WHILE (.NOT. (done .OR. v .EQ. 1))
         done = .TRUE.
         pmin = pcrit
         g = 0
         DO j = 2, n-1             ! find least p
           i = prev(j)
           IF (next(i) .NE. j) CYCLE       ! node already deleted
           
           IF (peach(j) < pmin) THEN
             pmin = peach(j)
             g = j
           END IF           
         END DO  ! next j
         
         IF (pmin < pcrit) THEN     ! merge slope
           done = .FALSE.
           i = prev(g)
           k = next(g)
           CALL CALC_SLOPE (i, k, n, CW, CWX, CWY, CWX2, CWXY, slope)
           IF (ISNAN(slope(i))) PRINT *, 'mullin:  NaN at merge!!!'
           prev(k) = i
           next(i) = k
           v = v - 1       ! one fewer segment           

*             consider merging new segment with previous
           IF (i > 1)
     &         CALL CALC_P (prev(i), i, k, n, CM, CM2, slope, peach(i))
           
*             consider merging new segment with next
           IF (k < n)
     &         CALL CALC_P (i, k, next(k), n, CM, CM2, slope, peach(k))
           
         END IF  ! p value less than threshhold?
       END DO  ! while merges still to make
       
*           arrays with bound depending on v
       ALLOCATE (IP(v+1), A(n,v+1), STAT=ierr)
       IF (ierr .NE. 0) RETURN
       
*-----------------------------------------------------------------------
  50   CONTINUE  ! least squares problem for piecewise line
*-----------------------------------------------------------------------
       DO j = 1, v+1
         DO i = 1, n
           A(i,j) = 0.
         END DO
       END DO
       
*                 the coefficient matrix
       j = 1
       h = 1
       k = next(h)
       DO i = 1, n
         IF (i .EQ. k .AND. i < n) THEN
           j = j + 1
           h = k
           k = next(k)
         END IF
         r = (X(i) - X(h)) / (X(k) - X(h))
         A(i,j) = 1. - r
         A(i,j+1) = r
       END DO
       
*                right-hand side
       DO i = 1, n
         cept(i) = Y(i)
       END DO
              
*    norm of A:  maximum of all column sums of magnitude
       tau = 0.
       DO j = 1, v+1
         r = 0.
         DO i = 1, n
           r = r + ABS(A(i,j))
         END DO
         tau = MAX(tau, r)
       END DO
       tau = tau * EPS
       
*           solve least-squares problem
       CALL HFTI (A, n, n, v+1, cept, n, 1, tau, k, 
     &               peach, scis, slope, IP)
       IF (k .EQ. 0) THEN
         PRINT *, 'hfti:  zero rank!'
         ierr = 4
         RETURN
       END IF
     
*               endpoint interpolation for slope
       i = 0
       j = 1
       DO WHILE (j .LE. n)
         i = i + 1
         scis(i) = X(j)
         j = next(j)
       END DO
       
       xnew = scis(1) ; ynew = cept(1)
       DO j = 1, v
         xold = xnew ; xnew = scis(j+1)
         yold = ynew ; ynew = cept(j+1)
         slope(j) = (ynew - yold) / (xnew - xold)
         IF (ISNAN(slope(j))) PRINT *, 'mullin: NaN at interp!!!'
       END DO
       
*-----------------------------------------------------------------------
*        require increasing slope in each segment
*-----------------------------------------------------------------------
       neg = .FALSE.
       vv = v
       h = 1
       DO k = 1, v
         IF (slope(k) .LE. 0.) THEN
           neg = .TRUE.
           PRINT *, 'negative slope at node ', k, ' of ', v
           IF (vv .EQ. 1) THEN 
             PRINT *, 'mullin:  Negative slope, unrecoverable error'
             ierr = 2
             RETURN
           ELSE IF (k .EQ. 1) THEN    ! merge with next segment
             j = next(next(1))
             prev(j) = 1
             next(1) = j
           ELSE IF (next(h) .EQ. n) THEN   ! merge with previous segment
             g = prev(h)             
             prev(n) = g
             next(g) = n
           ELSE   ! decide which segment to merge
             g = prev(h)
             i = next(h)
             j = next(i)
             CALL CALC_P (g, h, j, n, CM, CM2, slope, p_h)
             CALL CALC_P (g, i, j, n, CM, CM2, slope, p_i)
             
             IF (p_h < p_i) THEN ! merge with previous
               prev(i) = g
               next(g) = i
             ELSE                               ! merge with next
               prev(j) = h
               next(h) = j
             END IF
           END IF       ! j = {1, v} ?
           vv = vv - 1
         END IF    ! negative slope?
         IF (next(h) .GE. n) EXIT
         h = next(h)
       END DO    ! next node 
       v = vv
         
       IF (neg) GO TO 50   ! and solve the least-squares problem
       DEALLOCATE (prev, next, peach, cept, CM, CM2, CW, CWX, CWY, 
     &  CWX2, CWXY, IP, A, STAT=ierr)
       RETURN
         
      END ! of mullin
*++++++++++++++++++++++++ End of file mullin.f +++++++++++++++++++++++++
