*-----------------------------------------------------------------------
* lils.f - Linearizing Integral Least Squares
* by Andy Allinger, 2025, released to the public domain
*
*    Permission  to  use, copy, modify, and distribute this software and
*    its documentation  for  any  purpose  and  without  fee  is  hereby
*    granted,  without any conditions or restrictions.  This software is
*    provided "as is" without express or implied warranty.
*
* Parameter estimates for models of growth and vibration.
* Specifically:
*      Exponential growth
*      Diminishing returns / monomolecular model
*      Gompertz growth
*      Logistic growth
*      Korf growth model
*      Sinusoid oscillation
*      Damped vibration
*      Logarithmic spiral
*
*-----------------------------------------------------------------------

*-----------------------------------------------------------------------
* ratest - Estimate the growth rate.
* Input must be sorted in order of increasing X.
*
*___Name_____Type_______________In/Out____Description___________________
*   X(N)     Double precision   In        Independent variable (time)
*   Y(N)     Double precision   In        Quantity that grows
*   N        Integer            In        Number of data points
*   R        Double precision   Out       Growth rate
*   S(N)     Double precision   Neither   Work array
*   IERR     Integer            Neither   Error code
*-----------------------------------------------------------------------
      SUBROUTINE RATEST (X, Y, N, R, S, IERR)
       IMPLICIT NONE
       INTEGER N, IERR
       DOUBLE PRECISION X(N), Y(N), R, S(N)

       DOUBLE PRECISION ZERO, HALF
       PARAMETER (ZERO = 0.0D0, HALF = 0.5D0)
       INTEGER J
       DOUBLE PRECISION TS2, TSY
       LOGICAL DAFDIV

*-----------------------------------------------------------------------
*             Begin.
*-----------------------------------------------------------------------
       S(1) = ZERO
       DO J = 2, N
         IF (X(J) .LE. X(J-1)) THEN
           IERR = 2
           RETURN
         END IF
         S(J) = S(J-1) + HALF * (Y(J) + Y(J-1)) * (X(J) - X(J-1))
       END DO
       CALL DTRND1 (X, S, N)

       TS2 = ZERO
       TSY = ZERO
       DO J = 1, N
         TS2 = TS2 + S(J)**2
         TSY = TSY + S(J) * Y(J)
       END DO

       IF (DAFDIV(TSY, TS2)) THEN
         R = TSY / TS2
       ELSE
         IERR = 1
         RETURN
       END IF

       IERR = 0
       RETURN
      END  ! of ratest


*-----------------------------------------------------------------------
* freest - estimate the frequency
* Input must be sorted in order of increasing X.
*
*___Name_____Type_______________In/Out____Description___________________
*   X(N)     Double precision   In        Independent variable (time)
*   Y(N)     Double precision   In        Quantity that oscillates
*   N        Integer            In        Number of data points
*   F        Double precision   Out       Frequency
*   S(N)     Double precision   Neither   Work array
*   IERR     Integer            Neither   Error code
*-----------------------------------------------------------------------
      SUBROUTINE FREEST (X, Y, N, F, S, IERR)
       IMPLICIT NONE
       INTEGER N, IERR
       DOUBLE PRECISION X(N), Y(N), F, S(N)

       DOUBLE PRECISION ZERO, HALF, TWOPI
       PARAMETER (ZERO = 0.0D0,
     $            HALF = 0.5D0,
     $            TWOPI = 6.283185307179586477D0)

       INTEGER J
       DOUBLE PRECISION SP, SPP, TS22, TS2Y, W2
       LOGICAL DAFDIV

*-----------------------------------------------------------------------
*             Begin.
*-----------------------------------------------------------------------
       IF (N < 4) THEN
         IERR = 2
         RETURN
       END IF

       S(1) = ZERO
       DO J = 2, N
         IF (X(J) .LE. X(J-1)) THEN
           IERR = 4
           RETURN
         END IF
         S(J) = S(J-1) + HALF * (Y(J) + Y(J-1)) * (X(J) - X(J-1))
       END DO
       SP = ZERO
       DO J = 2, N
         SPP = SP
         SP = S(J)
         S(J) = S(J-1) + HALF * (S(J) + SPP) * (X(J) - X(J-1))
       END DO
       CALL DTRND2 (X, S, N)

       TS22 = ZERO
       TS2Y = ZERO
       DO J = 1, N
         TS22 = TS22 + S(J)**2
         TS2Y = TS2Y + S(J) * Y(J)
       END DO

       IF (DAFDIV(TS2Y, TS22)) THEN
         W2 = -TS2Y / TS22
       ELSE
         IERR = 1
         RETURN
       END IF

       IF (W2 < ZERO) THEN
         IERR = 3
         RETURN
       END IF

       F = SQRT(W2) / TWOPI
       IERR = 0
       RETURN
      END  ! of freest


*-----------------------------------------------------------------------
* korest - Estimate the growth rate according to the Korf model.
* Input must be sorted in order of increasing X.
*
*___Name_____Type_______________In/Out____Description___________________
*   X(N)     Double precision   In        Independent variable (time)
*   Y(N)     Double precision   In        Quantity that grows
*   N        Integer            In        Number of data points
*   R        Double precision   Out       Growth rate
*   LNY(N)   Double precision   Neither   Work array
*   S(N)     Double precision   Neither   Work array
*   IERR     Integer            Neither   Error code
*-----------------------------------------------------------------------
      SUBROUTINE KOREST (X, Y, N, R, LNY, S, IERR)
       IMPLICIT NONE
       INTEGER N, IERR
       DOUBLE PRECISION X(N), Y(N), R, LNY(N), S(N)

       DOUBLE PRECISION ZERO, HALF, ONE
       PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0)
       INTEGER J
       DOUBLE PRECISION DENOM, DN, NUMER, ODX, Q, SOX,
     $                   TS2OX2, TSOX2, TSOX, TOX, TOX2, TSOXY, TOXY, TY
       LOGICAL DAFDIV

*-----------------------------------------------------------------------
*             Begin.
*-----------------------------------------------------------------------
       DO J = 1, N
         IF (Y(J) .LE. ZERO) THEN
           IERR = 2
           RETURN
         END IF
         LNY(J) = LOG(Y(J))
       END DO
       S(1) = ZERO
       DO J = 2, N
         IF (X(J) - X(J-1) .LE. ZERO) THEN
           IERR = 3
           RETURN
         END IF
         S(J) = S(J-1) + HALF * (LNY(J) + LNY(J-1)) * (X(J) - X(J-1))
       END DO

       TS2OX2 = ZERO                             ! Collect sums.
       TSOX2 = ZERO
       TSOX = ZERO
       TOX = ZERO
       TOX2 = ZERO
       TSOXY = ZERO
       TOXY = ZERO
       TY = ZERO
       DN = DBLE(N)
       DO J = 1, N
         ODX = ONE / X(J)
         SOX = S(J) * ODX
         TOX = TOX + ODX
         TOX2 = TOX2 + ODX * ODX
         TSOX = TSOX + SOX
         TSOX2 = TSOX2 + SOX * ODX
         TS2OX2 = TS2OX2 + SOX * SOX
         TSOXY = TSOXY + SOX * LNY(J)
         TOXY = TOXY + ODX * LNY(J)
         TY = TY + LNY(J)
       END DO

*-----------------------------------------------------------------------
*             Solve system by Cramer's rule.
*-----------------------------------------------------------------------
       NUMER =   TSOXY * (DN*TOX2 - TOX**2)
     $         - TSOX2 * (DN*TOXY - TOX*TY)
     $         + TSOX  * (TOX*TOXY - TOX2*TY)

       DENOM =   TS2OX2 * (DN*TOX2 - TOX**2)
     $         - TSOX2  * (DN*TSOX2 - TOX*TSOX)
     $         + TSOX   * (TOX*TSOX2 - TOX2*TSOX)

       IF (DAFDIV(NUMER,DENOM)) THEN
         Q = NUMER / DENOM
       ELSE
         IERR = 1
         RETURN
       END IF
       R = ONE - Q

       IERR = 0
       RETURN
      END  ! of korest


*-----------------------------------------------------------------------
* damest - growth rate and frequency of a damped vibration
* Input must be sorted in order of increasing X.
*
* Given a dataset (x,y) and a model
*
*     y = c e^(r x) cos(2 pi f x - p) + k
*
* Determine the parameters r and f.  Amplitude, phase shift, and offset
* may then be found be linear regression.
*
*___Name_____Type_______________In/Out____Description___________________
*   X(N)     Double precision   In        Independent variable (time)
*   Y(N)     Double precision   In        Quantity that oscillates
*   N        Integer            In        Number of data points
*   R        Double precision   Out       Rate, negative for decay
*   F        Double precision   Out       Frequency
*   S1(N)    Double precision   Neither   Work array
*   S2(N)    Double precision   Neither   Work array
*   IERR     Integer            Neither   Error code
*-----------------------------------------------------------------------
      SUBROUTINE DAMEST (X, Y, N, R, F, S1, S2, IERR)
       IMPLICIT NONE
       INTEGER N, IERR
       DOUBLE PRECISION X(N), Y(N), R, F, S1(N), S2(N)

       INTEGER P
       DOUBLE PRECISION ZERO, HALF, TWOPI
       PARAMETER (P = 5,
     $            ZERO = 0.D0,
     $            HALF = 0.5D0,
     $            TWOPI = 6.283185307179586477D0)

       INTEGER I, J, K, IND(P)
       DOUBLE PRECISION DIF, TOL, X2, W2
       DOUBLE PRECISION A(P,P), B(P), C(P), SOLN(P), W(P)

       DOUBLE PRECISION DPMPAR

*-----------------------------------------------------------------------
*             Integrate.
*-----------------------------------------------------------------------
       IF (N < 5) THEN
         IERR = 2
         RETURN
       END IF
       TOL = DPMPAR(1)              ! Tolerance is machine epsilon.

       S1(1) = ZERO
       S2(1) = ZERO
       DO J = 2, N
         DIF = X(J) - X(J-1)
         IF (DIF .LE. ZERO) THEN
           IERR = 4
           RETURN
         END IF
         S1(J) = S1(J-1) + HALF * (Y(J) + Y(J-1)) * DIF
         S2(J) = S2(J-1) + HALF * (S1(J) + S1(J-1)) * DIF
       END DO

*-----------------------------------------------------------------------
*             Solve the system of linear equations.
*-----------------------------------------------------------------------
       DO J = 1, P
         DO I = 1, P
           A(I,J) = ZERO
         END DO
         B(J) = ZERO
       END DO

       DO J = 1, N                          ! Collect sums.
         X2 = X(J)**2
         A(1,1) = A(1,1) + S2(J)**2
         A(2,1) = A(2,1) + S1(J) * S2(J)
         A(3,1) = A(3,1) + X2 * S2(J)
         A(4,1) = A(4,1) + X(J) * S2(J)
         A(5,1) = A(5,1) + S2(J)
         A(2,2) = A(2,2) + S1(J)**2
         A(3,2) = A(3,2) + X2 * S1(J)
         A(4,2) = A(4,2) + X(J) * S1(J)
         A(5,2) = A(5,2) + S1(J)
         A(3,3) = A(3,3) + X2**2
         A(4,3) = A(4,3) + X2 * X(J)
         A(5,3) = A(5,3) + X2
         A(4,4) = A(4,4) + X2
         A(5,4) = A(5,4) + X(J)
         B(1) = B(1) + S2(J) * Y(J)
         B(2) = B(2) + S1(J) * Y(J)
         B(3) = B(3) + X2 * Y(J)
         B(4) = B(4) + X(J) * Y(J)
         B(5) = B(5) + Y(J)
       END DO

       DO J = 1, P-1                        ! Copy to upper triangle.
         DO I = J+1, P
           A(J,I) = A(I,J)
         END DO
       END DO
       A(5,5) = DBLE(N)

       CALL OSOLVD (P, P, A, P, B, TOL, K, SOLN, W, C, IND, IERR)
       IF (IERR .NE. 0) RETURN
       IF (K .NE. P) THEN
         IERR = 1
         RETURN
       END IF

       R = HALF * SOLN(2)
       W2 = -(SOLN(1) + R**2)

       IF (W2 < ZERO) THEN
         IERR = 3
         RETURN
       END IF

       F = SQRT(W2) / TWOPI
       IERR = 0
       RETURN
      END  ! of damest


*-----------------------------------------------------------------------
* spiest - estimate parameters for a logarithmic spiral
* Input must be sorted in order of increasing X.
*
* Given a dataset (x,y) and a model
*
*     y = c e^((r + iw)x - ip) + y_0
*
* where x is a real parameter and y is complex,
* determine the parameters r, f=w/2pi, c, p, and y0.
*
*___Name_______Type______In/Out____Description__________________________
*   X(N)       Double    In        Independent variable (time)
*   YRI(2*N)   Double    In        Coordinates in the complex plane
*                                     YRI(1:N) the real parts
*                                     YRI(N+1:2*N) the imaginary parts
*   N          Integer   In        Number of data points
*   RATE       Double    Out       Rate, negative for decay
*   FREQ       Double    Out       Frequency, cycles per unit X
*   AMP        Double    Out       Amplitude at initial time
*   PHASE      Double    Out       Phase offset, radians
*   Y0R        Double    Out       Real part of spiral center
*   Y0I        Double    Out       Imaginary part of spiral center
*   W(2*N,4)   Double    Neither   Work array
*   IERR       Integer   Neither   Error code
*-----------------------------------------------------------------------
      SUBROUTINE SPIEST (X,YRI,N,RATE,FREQ,AMP,PHASE,Y0R,Y0I,W,IERR)
       IMPLICIT NONE
       INTEGER N, IERR
       DOUBLE PRECISION X(N), YRI(2*N), RATE, FREQ, AMP, PHASE, Y0R,
     $                   Y0I, W(2*N,4)

       INTEGER Q
       DOUBLE PRECISION ZERO, ONE, TWOPI
       PARAMETER (Q = 4,
     $            ZERO = 0.D0,
     $            ONE = 1.D0,
     $            TWOPI = 6.283185307179586477D0)

       INTEGER J, K, IND(Q)
       DOUBLE PRECISION E, IMSIMY, IMSREY, RESIMY, RESREY, T, TSQ
       DOUBLE PRECISION COL(Q), SOLN(Q), TOL, WORK(Q)

       DOUBLE PRECISION DPMPAR
       LOGICAL DAFDIV

*-----------------------------------------------------------------------
*             Integrate real and imaginary parts.
*-----------------------------------------------------------------------
       IF (N < 5) THEN
         IERR = 2
         RETURN
       END IF
       TOL = DPMPAR(1)              ! Tolerance is machine epsilon.

       W(1,1) = ZERO
       W(1,2) = ZERO
       DO J = 2, N
         T = X(J) - X(J-1)
         IF (T .LE. ZERO) THEN
           IERR = 5
           RETURN
         END IF
         W(J,1) = W(J-1,1) + 0.5D0 * (YRI(J) + YRI(J-1)) * T
         W(J,2) = W(J-1,2) + 0.5D0 * (YRI(N+J) + YRI(N+J-1)) * T
       END DO
       CALL DTRND1 (X, W(1,1), N)                ! Detrend.
       CALL DTRND1 (X, W(1,2), N)

*-----------------------------------------------------------------------
*             Collect sums.
*-----------------------------------------------------------------------
       RESREY = ZERO
       IMSIMY = ZERO
       RESIMY = ZERO
       IMSREY = ZERO
       TSQ = ZERO
       DO J = 1, N
         RESREY = RESREY + W(J,1) * YRI(J)
         IMSIMY = IMSIMY + W(J,2) * YRI(N+J)
         RESIMY = RESIMY + W(J,1) * YRI(N+J)
         IMSREY = IMSREY + W(J,2) * YRI(J)
         TSQ = TSQ + W(J,1)**2 + W(J,2)**2
       END DO

*-----------------------------------------------------------------------
*             Growth rate and frequency
*-----------------------------------------------------------------------
       T = RESREY + IMSIMY
       IF (DAFDIV(T, TSQ)) THEN
         RATE = T / TSQ
       ELSE
         IERR = 3
         RETURN
       END IF
       T = RESIMY - IMSREY
       IF (DAFDIV(T, TSQ)) THEN
         FREQ = T / TSQ                          ! angular frequency
       ELSE
         IERR = 4
         RETURN
       END IF

*-----------------------------------------------------------------------
*             Solve the system of linear equations.
*-----------------------------------------------------------------------
       DO J = 1, N
         E = EXP(RATE * X(J))
         W(J,1) = E * COS(FREQ * X(J))
         W(N+J,1) = E * SIN (FREQ * X(J))
         W(J,2) = E * SIN(FREQ * X(J))
         W(N+J,2) = -E * COS(FREQ * X(J))
         W(J,3) = ONE
         W(N+J,3) = ZERO
         W(J,4) = ZERO
         W(N+J,4) = ONE
       END DO

       CALL OSOLVD (2*N,Q,W,2*N,YRI,TOL,K,SOLN,WORK,COL,IND,IERR)
       IF (IERR .NE. 0) RETURN
       IF (K .NE. Q) THEN
         IERR = 1
         RETURN
       END IF

       FREQ = FREQ / TWOPI
       AMP = SQRT(SOLN(1)**2 + SOLN(2)**2)
       PHASE = ATAN2(SOLN(2), SOLN(1))
       Y0R = SOLN(3)
       Y0I = SOLN(4)

       IERR = 0
       RETURN
      END  ! of spiest


*-----------------------------------------------------------------------
* dtrnd1 - remove linear trends
*
*___Name_____Type_______________In/Out___Description____________________
*   X(N)     Double precision   In       Independent variable
*   Y(N)     Double precision   Both     Dependent variable
*   N        Integer            In       # points in approximation
*-----------------------------------------------------------------------
       SUBROUTINE DTRND1 (X, Y, N)
       IMPLICIT NONE
       INTEGER N
       DOUBLE PRECISION X(N), Y(N)

       DOUBLE PRECISION ZERO
       PARAMETER (ZERO = 0.0D0)
       INTEGER I
       DOUBLE PRECISION A0, A1, DET, SX, SX2, SXY, SY, X1, X2, Y1

       LOGICAL DAFDIV             ! external function

*-----------------------------------------------------------------------
*             Begin
*-----------------------------------------------------------------------
       IF (N < 2) RETURN

*-----------------------------------------------------------------------
*             Least-squares sums
*-----------------------------------------------------------------------
       SX = ZERO
       SX2 = ZERO
       SY = ZERO
       SXY = ZERO
       DO I = 1, N
         X1 = X(I)
         Y1 = Y(I)
         X2 = X1 * X1
         SX = SX + X1             ! SUM x
         SX2 = SX2 + X2           ! SUM x^2
         SY = SY + Y1             ! SUM y
         SXY = SXY + X1 * Y1      ! SUM x @ y
       END DO

*-----------------------------------------------------------------------
*             Solve system of equations.
*-----------------------------------------------------------------------
       DET = N * SX2 - SX * SX
       A0 = SX2 * SY - SX * SXY
       A1 = N * SXY - SX * SY
       IF (.NOT. DAFDIV(A0, DET) .OR. .NOT. DAFDIV(A1, DET)) RETURN
       A0 = A0 / DET
       A1 = A1 / DET

*-----------------------------------------------------------------------
*             Detrend.
*-----------------------------------------------------------------------
       DO I = 1, N
         Y(I) = Y(I) - A0 - X(I) * A1
       END DO
       RETURN
      END  ! of dtrnd1


*-----------------------------------------------------------------------
* dtrnd2 - remove second order trend
*
*___Name_____Type_______________In/Out___Description____________________
*   X(N)     Double precision   In       Independent variable
*   Y(N)     Double precision   Both     Dependent variable
*   N        Integer            In       # points in approximation
*-----------------------------------------------------------------------
       SUBROUTINE DTRND2 (X, Y, N)
       IMPLICIT NONE
       INTEGER N
       DOUBLE PRECISION X(N), Y(N)

       DOUBLE PRECISION ZERO
       PARAMETER (ZERO = 0.0D0)
       INTEGER I
       DOUBLE PRECISION A0, A1, A2,
     $  DET, SUB1, SUB2, SUB3, SUB4, SUB5, SUB6,
     $  SX, SX2, SX3, SX4, SY, SXY, SX2Y, X1, X2, Y1

       LOGICAL DAFDIV             ! external function

*-----------------------------------------------------------------------
*           Begin
*-----------------------------------------------------------------------
       IF (N < 3) RETURN

*-----------------------------------------------------------------------
*             Least-squares sums.
*-----------------------------------------------------------------------
       SX = ZERO
       SX2 = ZERO
       SX3 = ZERO
       SX4 = ZERO
       SY = ZERO
       SXY = ZERO
       SX2Y = ZERO
       DO I = 1, N
         X1 = X(I)
         Y1 = Y(I)
         X2 = X1 * X1
         SX = SX + X1             ! SUM x
         SX2 = SX2 + X2           ! SUM x^2
         SX3 = SX3 + X2 * X1      ! SUM x^3
         SX4 = SX4 + X2 * X2      ! SUM x^4
         SY = SY + Y1             ! SUM y
         SXY = SXY + X1 * Y1      ! SUM x @ y
         SX2Y = SX2Y + X2 * Y1    ! SUM x^2 @ y
       END DO

*-----------------------------------------------------------------------
*             Parabolic approximation
*-----------------------------------------------------------------------
       SUB1 = SX2 * SX4 - SX3 * SX3
       SUB2 = SX * SX4 - SX2 * SX3
       SUB3 = SX * SX3 - SX2 * SX2
       SUB4 = SXY * SX4 - SX2Y * SX3
       SUB5 = SXY * SX3 - SX2Y * SX2
       SUB6 = SX2Y * SX - SXY * SX2
       DET = N*SUB1 - SX*SUB2 + SX2*SUB3
       A0 = SY*SUB1 - SX*SUB4 + SX2*SUB5
       A1 = N*SUB4 - SY*SUB2 + SX2*SUB6
       A2 = N*(-SUB5) - SX*SUB6 + SY*SUB3
       IF (     .NOT. DAFDIV(A0, DET)
     $     .OR. .NOT. DAFDIV(A1, DET)
     $     .OR. .NOT. DAFDIV(A2, DET) ) RETURN
       A0 = A0 / DET
       A1 = A1 / DET
       A2 = A2 / DET

*-----------------------------------------------------------------------
*             Detrend
*-----------------------------------------------------------------------
       DO I = 1, N
         Y(I) = Y(I) - A0 - X(I) * (A1 + X(I) * A2)
       END DO
       RETURN
      END  ! of dtrnd2