*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* rhythm.f - subroutines for music analysis
*  by Andy Allinger, 2011-2017, released to the public domain
*  This program may be used by any person for any purpose.
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* OSCILL - a frequency analyzer for pulse trains
* A bank of oscillatory filters responds to a sequence of events in time.
* Each filter has the undriven form of a decaying logarithmic spiral.
* Filters respond to events by adjusting their amplitude, phase, and frequency.
* For general reference, see:
*
*   Theodosius Pavlidis
*   Biological Oscillators:  their Mathematical Analysis
*   Academic Press, 1973
*
*   Frank C. Hoppensteadt
*   "Signal Processing by Model Neural Networks"
*   SIAM Review v.34 n.3 pp.426-444
*   Society for Industrial and Applied Mathematics, September 1992
*
*   Edward W. Large and John F. Kolen
*   "Resonance and the Perception of Musical Meter"
*   Connection Science v.6 n.1 pp.177-208, 1994
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*                              I/O LIST
*___Name_________Type______In/Out___Description_________________________
*   T(M)         Real      In       Onset times of events
*   S(M)         Real      In       Saliences of events
*   M            Integer   In       Number of events
*   N            Integer   In       Number of evaluation timepoints
*   O            Integer   In       Number of oscillators
*   STEP         Real      In       Step size in time
*   FDEEP        Real      In       Center frequency of deepest oscillator
*   FSPAC        Real      In       Frequency ratio of adjacent oscillators
*   FRANG        Real      In       Ratio of min to max frequency of a single
*                                     oscillator
*   A(0:O,0:N)   Real      Out      Amplitude
*   P(0:O,0:N)   Real      Out      Phase
*   F(0:O,0:N)   Real      Out      Frequency
*   IFAULT       Integer   Out      nonzero on error
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE OSCILL (T,S,M,N,O,STEP,FDEEP,FSPAC,FRANG,A,P,F,IFAULT)
       IMPLICIT NONE
       INTEGER M, N, O, IFAULT
       REAL T, S, STEP, FDEEP, FSPAC, FRANG, A, P, F
       DIMENSION T(M), S(M), A(0:O,0:N), P(0:O,0:N), F(0:O,0:N)

*         constants
       INTEGER Q
       REAL ETA, L, PI, RINIT, RMIN, RMAX, SM, TWOPI, V, ZETA
       PARAMETER (
     $             ETA = 1.,            ! target salience per event
*     $             g = 7.,             ! sharpening coefficient
     $             L = 3.,              ! half-life, seconds
     $             PI = 3.141592654,    !
     $             Q = 24,              ! number of initial phase offsets
     $             RINIT = 0.1,         ! initial amplitude
     $             RMIN = 0.0001,       ! distance to stay away from origin
     $             RMAX = 0.9999,       ! closest approach to unit circle
     $             SM = 0.05,           ! minimal salience
     $             TWOPI  = 2 * PI,     !
     $             V = 1.0,             ! volume of the beaker
     $             ZETA = 2.618034 )    ! target events per cycle

*        local variables
       INTEGER    H, I, J, K            ! counters

       REAL
     $    BETA,                   ! decay rate
     $    B1,                     ! decay for whole step
     $    BEST,                   ! best amplitude result
     $    BIG,
     $    DS,                     ! adjustment to each salience
     $    DT, DX, DY,             ! changes in t, x and y
     $    FMIN, F0, FMAX,         ! min, central, max freq. each oscillator
     $    FM,                     ! M as Real
     $    GAMMA, MU,              ! deviation, average of salience
     $    NU,                     ! frequency, in Hertz
     $    PHI                     ! phase, measured in beats elapsed
       REAL
     $    RHO, R2,                ! amplitude, amplitude squared
     $    S1,                     ! adjusted salience
     $    SMAX,                   ! greatest entry in S
     $    THETA,                  ! phase % 1, in radians
     $    T0,                     ! time of prior event
     $    TNEXT,                  ! time of next timepoint
     $    TPREV,                  ! time last processed event or step
     $    TT,                     ! same as T(i)
     $    U,                      ! sharpening coefficient given rho, theta
     $    W, XI,                  ! weight factor of salience, crit. fraction
     $    X, Y                    ! x = rho @ cos(theta) ; y = rho @ sin(theta)

       DOUBLE PRECISION 
     $    DEPS,                                ! machine double precision
     $    C2, ST0, ST1, ST2, SP, STP, DET,     ! least-squares sums
     $    SS, SS2                              ! salience sum, sum of squares

       LOGICAL E               ! whether exponentiation is required

!!!               stats
       INTEGER processed, ignored

*            external functions       
       REAL RMACHINE,
     $  SNAP, SLOG, SSIN, SCOS, SARC  ! swift EXP, LOG, SIN, COS, ATAN2
       DOUBLE PRECISION DMACHINE
   
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*             begin
       IFAULT = 0         ! clear error codes
       BIG = RMACHINE(3)
       DEPS = DMACHINE(1)

*       validate
       IF (M < 1 .OR. N < 1 .OR. O < 1) THEN
         IFAULT = 1
         RETURN
       END IF

*       check that enough timepoints have been allocated for the given events
       IF (T(M) > N * STEP) THEN
         PRINT *, 'oscill: t(m) is ', T(M), ' n@step is ', N * STEP
         IFAULT = 2
         RETURN
       END IF

*         get f0, take root of frang
       F0 = FDEEP / FSPAC    ! OK to multiply by fspac on iteration k = 1
       FMIN = F0 / SQRT(FRANG)
       FMAX = F0 * SQRT(FRANG)
       BETA = LOG(0.5) / L
       B1 = EXP(BETA * STEP)
       
*        Total salience
       SS = 0.D0
       SS2 = 0.D0
       SMAX = 0.
       DO I = 1, M
         SS = SS + S(I)
         SS2 = SS2 + S(I) **2
         IF (S(I) > SMAX) SMAX = S(I)
       END DO
       FM = FLOAT(M)
       MU = SNGL(SS / FM)
       GAMMA = 0.6745 * SNGL(SQRT(SS2 / FM - (SS / FM) **2))
       
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*         for each oscillator
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       DO K = 1, O
         F0 = F0 * FSPAC           ! central frequency
         FMIN = FMIN * FSPAC       ! minimum frequency
         FMAX = FMAX * FSPAC       ! maximum frequency

         BEST = -BIG 
         XI = (FM - ZETA * T(M) * F0) / FM     ! fraction of events to ignore
         XI = MIN(MAX(RMIN, XI), RMAX)
         DS = MU + GAMMA * TAN(PI * (XI - 0.5))  ! Cauchy cumulative dist.
         DS = MIN(MAX(0., DS), SMAX - SM)
         SS = 0.D0
         J = 0               ! count nonzero events
         DO I = 1, M
           IF (S(I) > DS) THEN
             SS = SS + (S(I) - DS)
             J = J + 1
           END IF
         END DO
         W = ETA * J / MAX(SM, SNGL(SS))

         ignored = 0
         processed = 0

         DO H = 0, Q-1              ! for each initial phase
           E = .FALSE.               ! no stimulus event intervening
           RHO = RINIT                            ! initial amplitude
           PHI = FLOAT(H) / FLOAT(Q)              ! initial phase
           NU = F0                                ! default frequency
           A(0,0) = RHO
           P(0,0) = PHI
           F(0,0) = NU

           I = 1
           T0 = 0.
           TNEXT = 0.
           TPREV = 0.
           ST0 = 0.D0
           ST1 = 0.D0
           ST2 = 0.D0
           SP = 0.D0
           STP = 0.D0

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*          for each evaluation point
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
           DO J = 1, N
             TNEXT = TNEXT + STEP
  10         IF (I > M) GO TO 50       ! no more events ?
             TT = T(I)
             IF (TT > TNEXT) GO TO 50   ! no event before tnext ?

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*   apply the stimulus - first advance r, phi to present
*    then make an adjustment in Cartesian coordinates
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*                phi' = phi + nu @ dt
             DT = TT - TPREV
             PHI = PHI + NU * DT

*               rho' = rho @ e^(b @ dt)
             RHO = RHO * SNAP(BETA * DT)     ! approximate exponential decay
             RHO = MIN(MAX(RMIN, RHO), RMAX)

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*  Ignore stimuli below critical value.
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             S1 = S(I) - DS
             IF (S1 .LE. 0.) THEN
               ignored = ignored + 1
               I = I + 1
               TPREV = TT
               E = .TRUE.
               GO TO 10
             ELSE
               processed = processed + 1
               S1 = S1 * W
             END IF

             THETA = TWOPI * MOD(PHI, 1.0)        ! 0 .LE. theta < 2 @ pi
             X = RHO * SCOS(THETA)                ! cosine
             Y = RHO * SSIN(THETA)                ! sine

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* Use activation sharpening.  Refer to:
*   J. Devin McAuley
*   Perception of Time as Phase:  Toward an Adaptive Oscillator Model of
*   Rhythmic Pattern Processing
*   University of Indiana Ph.D. Thesis, July 1995, p. 58
*           u = ((1 + COS(theta)) / 2) ^ (g @ r)
*
* also refer to:
*   J. Devin McAuley
*   Learning to Perceive and Produce Rhythmic Patterns
*   in an Artificial Neural Network
*   Indiana University Department of Computer Science
*   Technical Report 371, 1 February 1993
* which cites:
*   Carme Torras
*   Temporal-Pattern Learning in Neural Models
*   Berlin:  Springer-Verlag, 1985
*
* Another choice might be the wrapped Cauchy distribution.  Refer to:
*   Claudio Agostinelli
*   R CRAN package circular
*   22 May 2006
*           u = (1 - r^2) / ( 2 @ pi @ (1 + r^2 - 2 @ r @ COS(theta)) )
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*           u = ((1 + COS(theta)) / 2) ** (g * rho)
             R2 = RHO **2
             U = (1. - R2) / (TWOPI * (1. + R2 - 2. * X))

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*  Consider a chemical oscillator with reactant concentrations x and y,
* subject to a stimulus of adding a substance with y = 0, and x = v,
* the maximum concentration of x.  Then the phase transition will fulfill:
*                dy/dx = -y / (v - x)
* Refer to:
*            Arthur T. Winfree
*            The Geometry of Biological Time
*            Springer-Verlag, 1980, p. 139
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             U = (U * S1) / SQRT((V-X) **2 + Y **2)
             DX = (V-X) * U
             DY = (-Y) * U
             X = MIN(X + DX, V)
             IF (ABS(DY) .GE. ABS(Y)) THEN   ! don't overcorrect
               IF (Y < 0.) PHI = PHI + 1.  ! increment phase
               Y = 0.
             ELSE
               Y = Y + DY
             END IF

             RHO = SQRT(X **2 + Y **2)
             IF (RHO < RMIN) THEN  ! don't fall on zero (the phase singularity)
               X = RMIN
               RHO = RMIN
             END IF
             IF (RHO > RMAX) RHO = RMAX

             THETA = SARC(Y, X)       ! arctangent
             PHI = AINT(PHI) + THETA / TWOPI

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*   Measure frequency as the slope of a line of phase as a function of time
*                 nu = dphi / dt
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*           decay the running sums for the least squares line
             C2 = DBLE(SNAP(BETA * (TT-T0)))     ! decay factor
             ST0 = ST0 * C2 + 1.D0               ! SUM t^0
             ST1 = ST1 * C2 + DBLE(TT)           ! SUM t^1
             ST2 = ST2 * C2 + DPROD(TT, TT)      ! SUM t^2
             SP = SP * C2 + DBLE(PHI)            ! SUM phi
             STP = STP * C2 + DPROD(TT, PHI)     ! SUM t @ phi
             DET = ST0 * ST2 - ST1 * ST1         ! determinant
             IF (DET > DEPS) NU = SNGL( (ST0 * STP - SP * ST1) / DET )
!             NU = MIN(MAX(FMIN, NU), FMAX)   ! impose limits
             IF (NU < FMIN .OR. NU > FMAX) THEN      ! reset
               NU = F0
               RHO = RMIN
               ST0 = 0.
               ST1 = 0.
               ST2 = 0.
               SP = 0.
               STP = 0.
             END IF

             I = I + 1
             T0 = TT
             TPREV = TT
             E = .TRUE.
             GO TO 10

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  50       CONTINUE   ! advance time to an evaluation point for A, P, F
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
             IF (E) THEN
               DT = TNEXT - TPREV
               RHO = RHO * SNAP(BETA * DT)  ! approximate exponential function
               PHI = PHI + NU * DT
             ELSE
               RHO = RHO * B1       ! avoid an exponentiation when possible
               PHI = PHI + NU * STEP
             END IF
             IF (RHO < RMIN) RHO = RMIN

             A(0,J) = RHO
             P(0,J) = PHI
             F(0,J) = NU
             TPREV = TNEXT
             E = .FALSE.
           END DO  ! next j

*            find sum of logarithm of amplitude
           U = 0.
           DO J = 1, N
             U = U + SLOG(A(0,J))
           END DO
           IF (U > BEST) THEN              ! new best answer
             BEST = U
             DO J = 0, N
               A(K,J) = A(0,J)
               P(K,J) = P(0,J)
               F(K,J) = F(0,J)
             END DO
           END IF
         END DO  ! next h
!         PRINT *, 'osc# ', k, ' processed ', processed, ' events ', 
!     $      'and ignored ', ignored, 
!     $ ' using ds: ', DS,
!     $    ' using w: ', W
         
       END DO  ! next k

       RETURN
      END  ! of OSCILL


*-----------------------------------------------------------------------
* orlp - logical OR of logarithms of probabilities
*-----------------------------------------------------------------------
      SUBROUTINE ORLP (X, N, RESULT)
       IMPLICIT NONE
       INTEGER N
       REAL X, RESULT
       DIMENSION X(N)
       
       INTEGER J
       REAL XMAX

*-----------------------------------------------------------------------
*       Begin.
*-----------------------------------------------------------------------
       IF (N < 1) RETURN
       XMAX = X(1)
       DO J = 2, N
         IF (X(J) > XMAX) XMAX = X(J) 
       END DO 

       RESULT = 0.
       DO J = 1, N
         RESULT = RESULT + EXP(X(J) - XMAX)
       END DO
       RESULT = LOG(RESULT) + XMAX
       RETURN
      END  ! of ORLP



*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
* strata -  The best rhythic patterns are chosen from the output of the 
* oscillators with the Viterbi algorithm.  Dropped beats are filled in.
* The phase and time estimates are approximated with a smoothing a spline 
* yielding a function for phase given time, for each of several rhythmic 
* strata.  Refer to :
*  
*   Anssi P. Klapuri, Antti J. Eronen, and Jaakko T. Astola,
*   "Analysis of the Meter of Acoustic Musical Signals"
*   IEEE Transactions on Speech and Audio Processing, 2004
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*                          TRANSFER LIST
*__Name___________Type______In/Out___Description_________________________
*  T(m)           Real      in       times of events
*  S(m)           Real      in       saliences of events
*  L              Integer   in       number of levels ("strata") to make
*  TL             Integer   in       index of the tactus level
*  m              Integer   in       number of events
*  n              Integer   in       number of evaluation timepoints 
*  NODES          Integer   in       number of spline interpolation points
*  step           Real      in       duration between oscillator evaluations
*  urpar(8)       Real      in       parameters for probability functions
*  opt_S          Logical   in       option to calculate an average tempo
*  opt_T          Logical   in       option to get tempo segments
*  full           Logical   in       calculate all the strata
*  prob           Real      out      total log probability 
*  CX(NODES+1)    Real      out      spline coefficients
*  CY(NODES,L)    Real      out          "          "
*  CA(NODES,L)    Real      out          "          "
*  CB(NODES,L)    Real      out          "          "
*  CC(NODES,L)    Real      out          "          "
*  TX(0:n)        Real      out      tempo line segment coefficients
*  TA(0:n)        Real      out             "          "
*  nt             Integer   out      number of tempo line segments
*  tempo          Real      out      average tempo
*  A(0:o,0:n,2)   Real      neither  oscillator amplitude
*  P(0:o,0:n,2)   Real      neither  oscillator phase
*  F(0:o,0:n,2)   Real      neither  oscillator output frequency
*  B(o,n)         Real      neither  back pointers 
*  C(0:n,L)       Real      neither  best paths thru the Viterbi lattice
*  W(0:n)         Real      neither  weights for least-squares spline
*  X(0:n)         Real      neither  time input for spline
*  Y(0:n)         Real      neither  phase input for spline
*  Z(7*NODES+25)  Real      neither  work array for spline
*  fault          Integer   out      error code
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE strata (T, S, L, TL, M, N, NODES, step, urpar, opt_S, 
     $ opt_T, usecfg, full, prob, CX, CY, CA, CB, CC, TX, TA, nt, 
     $ tempo, A, P, F, B, C, W, X, Y, Z, fault)
      
       IMPLICIT NONE
       INTEGER o
       PARAMETER (o = 27)          ! number of oscillators
       INTEGER L, TL, m, n, NODES, nt, B, C, fault
       REAL T, S, step, urpar, prob, CX, CY, CA, CB, CC, TX, TA, tempo,
     $  A, P, F, W, X, Y, Z
       LOGICAL opt_S, opt_T, usecfg, full
       DIMENSION T(m), S(m), urpar(8), CX(NODES+1), CY(NODES,L), 
     $  CA(NODES,L), CB(NODES,L), CC(NODES,L), TX(0:n), TA(0:n), 
     $  A(0:o,0:n,2), P(0:o,0:n,2), F(0:o,0:n,2), B(o,n), C(0:n,L), 
     $  W(0:n), X(0:n), Y(0:n), Z(7*NODES+25)

      
*         local variables
       INTEGER GN 
       PARAMETER (GN = 1024)  ! evaluation points of spline
     
       INTEGER 
     $         back,                     ! back pointer
     $         h, i, ii, j, k, k1,       ! counters
     $         jnew, jold                ! even or odd time step
     
       REAL  
     $  afact,                 ! scale down amplitude importance
     $  BIG,                   ! maximum number 
     $  CP(4),                 ! cubic polynomial coefficients
     $  dp, dt,                ! change of phase, change of time
     $  g,                     ! the Large-Klapuri g function
     $  gfact,                 ! scale down resonance importance
     $  GX(gn), GY(gn), GA(gn), GB(gn), GC(gn),  ! g function parameters
     $  mu, nu,                ! average frequency, present frequency
     $  pobs, ppre,            ! observed, predicted phase
     $  r,                     ! ratio of frequencies
     $  rpar(8),               ! rhythm parameters
     $  sigi,                  ! parameters for frequency change, initial
                               ! frequency distribution, and phase error
     $  tcum,                  ! cumulative time
     $  twovf, twovi, twovp,   ! double variance
     $  U(o,0:1),              ! cumulative log probabilities at even/odd steps
     $  ubest, utry,           ! log-prob. of likliest, other transitions
     $  V(o)                   ! log-prob. of each prior state
     
       LOGICAL first,
     $         flip(o)

       SAVE first, GX, GY, GA, GB, GC, rpar
       
*           Constants       
*         AMPIMP parameter will be scaled by STEP.  Thus, it is per
*           unit time.  Figure 4 eighth notes per second, at 5% phase error
*           gives a penalty of -0.5.  A resonator tracking at 10% amplitude
*           has a log of about -2.3
       INTEGER minstp              ! minimum timesteps for SPFIT
       REAL ampimp, gfnimp,        ! importance of amplitude, resonance
     $      fdeep, frang, fspac    ! deepest natural frequency, range, spacing
       PARAMETER (ampimp = 0.2,       ! because amplitude is not probability
     $            gfnimp = 0.3,       ! not a real probability either
     $            fdeep = 0.1,        ! center frequency of slowest oscillator
     $            frang = 4. / 3.,    ! range of each oscillator
     $            fspac = 38. / 31.,  ! chosen to avoid spurious resonances
     $            minstp = 7)         ! don't break spline interpolator

*         External functions
       INTEGER CEIL            ! round up
       REAL
     $  rmachine,              ! machine constants, single precision
     $  SLOG                   ! Swift logarithm
     
!!            timing info
       REAL TARRAY(2), RESULT

*        The default rhythm parameters are from Klapuri, Eronen, & Astola
       DATA rpar /  0.18,         ! tatum period
!     $              0.39,         ! standard deviation of log ratio
     $              0.27,         ! standard deviation of log ratio
     $              0.55,         ! tactus period
     $              0.28,         ! standard deviation of log ratio
     $              2.10,         ! measure period
     $              0.26,         ! standard deviation of log ratio
!     $              0.2,          ! std deviation of tempo changes
     $              0.1,          ! std deviation of tempo changes
     $              0.1 /         ! std deviation of phase error

       DATA first / .TRUE. / ! calculate the g function on first run
      
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*      begin
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*         define constants
       BIG = rmachine(3)
       afact = ampimp * step
       gfact = gfnimp
       IF (first) THEN
         CALL gfunc(GX, GY, GA, GB, GC, gn)
         first = .FALSE.
       END IF
       
       jnew = 0        ! avoid meaningless compiler warnings
       
*           custom config data
       IF (usecfg) THEN
         DO i = 1, 8
           rpar(i) = urpar(i)
         END DO
       END IF

*          compare number of strata against number of oscillators
       IF (L > o) THEN
         PRINT *, 'Too many rhythmic strata for oscillators allocated'
         fault = 1
         RETURN
       END IF
        
        CALL DTIME (tarray, result)
        
*           set up evaluation points
       dt = (n + 0.5) * step / NODES
       tcum = 0.
       DO j = 1, NODES+1
         CX(j) = tcum
         tcum = tcum + dt
       END DO

*          run the data thru the filterbank
       CALL OSCILL (T,S,m,n,o,step,fdeep,fspac,frang,A,P,F,fault)
       IF (fault .NE. 0) PRINT *, 'oscill returns error #', fault
       
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!          PRINT *, 'forward oscillators: '
!  41       FORMAT ('step: ', I5)
!  42       FORMAT (F9.3, $)
!  43       FORMAT ()
!           DO j = 1, n
!             print 41, j
!             DO h = 1, o
!               PRINT 42, A(h,j,1)
!             END DO
!             PRINT 43
!             DO h = 1, o
!               print 42, P(h,j,1)
!             END DO
!             PRINT 43
!             DO h = 1, o
!               PRINT 42, F(h,j,1)
!             END DO
!             PRINT 43
!           END DO
!!!!!!!!!!!!!!!!!!!!
       
*               Reverse order of the data
       tcum = n * step
!       PRINT *, 'using tcum of ', tcum
       DO i = 1, m                   ! negate time
         T(i) = -T(i) + tcum
         T(i) = MIN(MAX(0., T(i)), tcum - .001)
       END DO
       DO i = 1, m/2                 ! swap order
         r = T(i)
         T(i) = T(m-i+1)
         T(m-i+1) = r
         r = S(i)
         S(i) = S(m-i+1)
         S(m-i+1) = r
       END DO
       
*         Filter the reversed sequence.
       CALL OSCILL (T, S, m, n, o, step, fdeep, fspac, frang,
     $  A(0,0,2), P(0,0,2), F(0,0,2), fault)
       IF (fault .NE. 0) THEN
         PRINT *, 'rev. ocsill returns #', fault
         PRINT *, 'T(m): ', T(m)
         PRINT *, 'n*step: ', N * STEP
         PRINT *, 'tcum: ', tcum
       END IF
           
*            Reflect phase
       DO h = 1, o
         dp = FLOAT(CEIL(P(h,n,2)))
         DO j = 0, n
           P(h,j,2) = -P(h,j,2) + dp
         END DO
       END DO

*          Keep backwards result until forwards pass has higher amplitude
       DO h = 1, o
         flip(h) = .TRUE.
       END DO
       DO j = 0, n
         DO h = 1, o
           IF (flip(h)) THEN
             IF (A(h,j,1) < A(h,n-j,2)) THEN
               A(h,j,1) = A(h,n-j,2)
               P(h,j,1) = P(h,n-j,2)
               F(h,j,1) = F(h,n-j,2)
             ELSE
               flip(h) = .FALSE.
             END IF
           END IF
         END DO
       END DO
           
        CALL DTIME (tarray, result)
!        PRINT *, 'rhythm:  time filtering: ', RESULT

       prob = 0.
       twovf = 2 * rpar(7) **2        ! same for all strata
       twovp = 2 * rpar(8) **2

       DO h = 1, L          ! for each stratum
        CALL DTIME (tarray, result)

*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
*            Viterbi algorithm:
*    The basic equation for conditional probability is:
*                     P(A&B) = P(A|B) @ P(B)
*   Suppose each event e in a sequence of n events is conditional only on 
* the event before it.  This yields the equation for a Markov chain:
*     P(E) = P(e[1]) @ PRODUCT P(e[i+1]|e[i])              where i = 1...n-1
*  where P(E) is the probability of the entire sequence.
*  Suppose also that knowledge of each event is supplemented by an observation
*  at each step.  Applying Bayes's rule:
*               P(E&O) = P(E|O) @ P(O) = P(O|E) @ P(E)  
*  disregarding P(O) which is the same for all events,
*                 P(E|O) = P(O|E) @ P(E)
* Bayesian reasoning is less elegant when the quantities involved must
* be measured.  For the present problem, transition probability P(e[i+1]|e[i])
* is based on the assumption that frequency remains steady and therefore phase
* advances steadily.  The observation likelihood, P(O|E), is taken to be 
* proportional to the oscillator amplitude, and to be higher if a frequency
* is proportional to the tactus (first detected, (strongest?)) frequency.
* The Viterbi algorithm asks at each step in the sequence, for each possible 
* event, what event was most likely at the prior step?  The answers for each
* step are kept, and at the end of the sequence the most likely chain of
* events may be traced back.
*
*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
*          choose parameters for initial frequency distribution
         IF (h .EQ. 1) THEN                 ! tatum
           mu = rpar(1)
           sigi = rpar(2)
         ELSE IF (h .EQ. TL) THEN           ! tactus
           mu = rpar(3)
           sigi = rpar(4)
         ELSE IF (h .EQ. TL+1) THEN    ! measure
           mu = rpar(5)
           sigi = rpar(6)
         ELSE IF (h < TL) THEN         ! subtatum?
           mu = 0.09          
           sigi = 0.40
         ELSE                           ! strophe?
           mu = 4.2
           sigi = 0.40
         END IF
         twovi = 2. * sigi **2
         IF (h > 1) gfact = gfnimp / FLOAT(h - 1)
         
*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
* Log-normal initial frequency distribution:
*
*   p(T) = (1 / T sigma (2@pi)^(1/2)) @ e^((-1/2@sigma^2) @ (ln(T/m))^2)
*
*  log probability, disregarding constants is then:
*
*   ln p(T) = (-1/2@sigma^2) @ ln(T/m)^2 - ln(T)
*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
         DO k = 1, o
           nu = F(k,0,1)
           r = mu * nu                 ! frequency @ time: non-dimensional
           IF (r > 1.) r = 1. / r
           r = SLOG(r)
           r = -(r **2) / twovi 
           IF (nu > 0.) r = r + LOG(nu)
           U(k,0) = r 
         END DO
         
*          forwards phase, for each time step
         DO 30 j = 1, n
           IF (BTEST(j, 0)) THEN  ! alternate storage in array U
             jnew = 1  
             jold = 0
           ELSE       
             jnew = 0
             jold = 1
           END IF
           
*           for each oscillator at present time step
           DO 20 k = 1, o
             IF (A(k,j,1) < 0.) THEN      ! this cell has been marked 
               U(k,jnew) = -BIG
               GO TO 20
             END IF
           
             ubest = -BIG
             nu = F(k,j,1)
             ppre = P(k,j,1) - nu * step    ! predicted value for prior phase
             
*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
*                transition probabilities
*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
*              for each oscillator at previous time step
             ii = 0
             DO 10 i = 1, o
               IF (A(i,j-1,1) < 0.) GO TO 10   ! this cell has been marked
               ii = ii + 1
               utry = U(i,jold)              ! prior probability

*                 penalize change in frequency 
               r = nu / F(i,j-1,1)
               IF (r > 1.) r = 1. / r       ! restrict to 0 < r .LE. 1
               r = SLOG(r)
               utry = utry - (r **2) / twovf
               
*                penalize phase mismatch
               dp = ppre - P(i,j-1,1)
               r = dp - ANINT(dp)
               utry = utry - (r **2) / twovp               
               V(ii) = utry
               
               IF (utry > ubest) THEN 
                 ubest = utry
                 back = i
               END IF
  10         CONTINUE  ! next i
             
!!            DEBUG:
       IF (ubest .LE. -BIG) THEN 
         PRINT *, 'no way out from j=', j, ' k=', k, 'ii=', II
         PRINT *, 'unew=', U(k,jnew)
       END IF
             
*              save pointer to most likely prior state
             B(k,j) = back
             CALL ORLP (V, ii, utry)  ! add probabilities
               
*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
*               observation likelihood
*|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
*               oscillator amplitude 
             utry = utry + afact * SLOG(A(k,j,1))
             
*              is frequency a Farey relation to previous strata?
             DO i = 1, h-1
               r = F(C(j,i),j,1) / nu
               IF (r > 1.) r = 1. / r           ! evaluate g function
               CALL SCOMP(GX, GY, GA, GB, GC, gn, r, g, 1, fault)
               IF (fault .NE. 0) PRINT *, 'scomp: error#', fault
               utry = utry + gfact * g
             END DO
               
*               does frequency match prior expected frequency?
             r = mu * nu 
             IF (r > 1.) r = 1. / r
             r = SLOG(r)
             r = -(r **2) / twovi
             IF (nu > 0.) r = r + LOG(nu)
             
             utry = utry + r  
             U(k,jnew) = utry
  20       CONTINUE  ! next k
  30     CONTINUE ! next j


*         find the highest scoring path
         ubest = -BIG
         
         DO k = 1, o
           IF (U(k,jnew) > ubest) THEN
             back = k
             ubest = U(k,jnew)
           END IF
         END DO
         
*          cumulative probability 
         prob = prob + ubest
         
*          trace back pointers
         DO j = n, 1, -1
           IF (back < 1 .OR. back > o) THEN
             PRINT *, 'bad back pointer! :', back, ' j=', j, ' n=', n
             PRINT *, 'B= ', B
           END IF
           C(j,h) = back
           back = B(back,j)
         END DO
         C(0,h) = back                   ! End of Viterbi         

!!!!!!!!!!!!!! debugging dump
!         PRINT *, 'stratum ', h, ' C[]='
! 45      FORMAT (I3, $)
!         DO j = 0, n
!           PRINT 45, C(j,h)
!         END DO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        CALL DTIME (tarray, result)
!        PRINT *, 'rhythm: level ', h, ' time in DP lattice: ', result
         
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*           replace dropped beats
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         W(0) = A(back,0,1)
         X(0) = 0.
         Y(0) = P(back,0,1)
         tcum = 0. 
         k1 = back
         DO j = 1, n
           k = k1          ! index of previous
           k1 = C(j,h)    ! index of new
           ppre = P(k,j-1,1) + step * F(k1,j,1)      ! backwards integral
           pobs = P(k1,j,1)
           pobs = pobs + ANINT(ppre - pobs)   ! adjust by an integer 
           P(k1,j,1) = pobs
           
           W(j) = A(k1,j,1)
           A(k1,j,1) = -1.     ! prevent this cell from being used again
           tcum = tcum + step
           X(j) = tcum
           Y(j) = pobs
         END DO  ! next j
         
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*          fit a least-squares spline
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         IF (n < minstp) THEN            ! fit a cubic polynomial
           CALL cfit(W, X, Y, n+1, CP, fault)
           IF (fault .NE. 0) PRINT *, 'cfit: error#', fault
           CY(1,h) = CP(1)
           CA(1,h) = CP(2)
           CB(1,h) = CP(3)
           CC(1,h) = CP(4)

*                restore X
           X(0) = 0.
           DO J = 1, N
             X(J) = X(J-1) + step
           END DO
         ELSE
         
*      least squares spline fitting by deBoor and Morris (NSWC library)
           CALL SPFIT (X, Y, W, n+1, CX, NODES+1, 
     $       CY(1,h), CA(1,h), CB(1,h), CC(1,h), Z, fault)
           IF (fault .NE. 0) PRINT *, 'spfit: error #', fault           
         END IF
         
*                tactus level only
         IF (h .EQ. TL) THEN 

*              average tempo (for dance mode)  
           IF (opt_S) THEN
             CALL wlsl (W, X, Y, n+1, r, tempo, fault)  ! re-use r
             IF (fault .NE. 0) PRINT *, 'wlsl:  error#', fault
           END IF

*              line-segment approximation of tempo curve
           IF (opt_T) THEN
             CALL mullin (W, X, Y, TX, TA, n+1, nt, fault)
             IF (fault .NE. 0) PRINT *, 'mullin:  error #', fault 
           END IF ! opt_T ?
           
             CALL DTIME (tarray, result)
!            PRINT *, 'rhythm: level ', h, ' time with splines: ', result
           
           IF (.NOT. full) RETURN  ! don't need all the splines

         END IF  ! tactus level ?
       END DO ! next h
       
       RETURN
      END  ! of strata



*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*  add data to rhythm parameters
*
*__Name________Type______________In/Out___Description___________________
*  T(M)        Real              In       event times
*  S(M)        Real              In       saliences of events
*  M           Integer           In       # of events
*  N           Integer           In       # evaluation points
*  STEP        Real              In       time between evaluation points
*  A(0:o,0:n)  Real              Neither  amplitude for whole filterbank
*  P(0:o,0:n)  Real              Neither  phase for whole filterbank
*  F(0:o,0:n)  Real              Neither  frequency for whole filterbank
*  rsum(9)     Double Precision  Both     sums & sums of squares (see next f'n)
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE adrp (T, S, m, n, step, A, P, F, rsum)
       IMPLICIT NONE
       INTEGER o
       PARAMETER (o = 3)             ! how many oscillators
       INTEGER M, N
       REAL T, S, STEP, A, P, F
       DOUBLE PRECISION rsum
       DIMENSION T(m), S(m), A(0:o,0:n), P(0:o,0:n), F(0:o,0:n), rsum(9)
       
*         local variables
       INTEGER  fault, 
     $          i, j          ! counters
       
       REAL 
     $      dp,                ! change in phase
     $      e,                 ! remainder phase error
     $      fdeep, fspac, frang,              ! input to OSCILL
     $      flog,              ! logarithm of ratio of consecutive frequencies
     $      lnt,               ! logarithm of period (to match Klapuri)
     $      nu                 ! frequency
     
       SAVE fdeep, fspac, frang 
       DATA fdeep, fspac, frang / 0.5, 4., 5. /
       
*        begin
       fault = 0
       
*            one oscillator for tatum, tactus, measure
       CALL oscill (T,S,m,n,o,step,fdeep,fspac,frang,A,P,F,fault)
       IF (fault .NE. 0) PRINT *, 'adrp: error in oscillators'
       
       DO i = 1, 3                ! for each stratum
         DO j = 1, n
           nu = F(i,j)
           lnt = -LOG(nu)
           IF (i .EQ. 1) THEN          ! measure
             rsum(5) = rsum(5) + lnt
             rsum(6) = rsum(6) + DPROD(lnt, lnt)
           ELSE IF (i .EQ. 2) THEN       ! tactus
             rsum(3) = rsum(3) + lnt
             rsum(4) = rsum(4) + DPROD(lnt, lnt)
           ELSE                           ! tatum
             rsum(1) = rsum(1) + lnt
             rsum(2) = rsum(2) + DPROD(lnt, lnt)
           END IF
           
*               frequency change
           flog = LOG( nu / F(i,j-1) )
           rsum(7) = rsum(7) + DPROD(flog, flog)
           
*                  phase mismatch
           dp = P(i,j) - nu * step           ! what prior phase should be
           dp = P(i,j-1) - dp                 ! difference
           e = dp - ANINT(dp)                  ! remainder phase error
           rsum(8) = rsum(8) + DPROD(e, e)
         END DO ! next j
       END DO ! next i

*        count total events
       rsum(9) = rsum(9) + n
       RETURN
      END ! of adrp
       
       
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*  rsum   DOUBLE PRECISION  rsum(1)        sum of tatum log period
*                           rsum(2)        sum of square tatum log period
*                           rsum(3)        sum of tactus log period
*                           rsum(4)        sum of square tactus log period
*                           rsum(5)        sum of measure log period
*                           rsum(6)        sum of square measure log period
*                           rsum(7)        sum of square log frequency ratio
*                           rsum(8)        sum of square phase error
*                           rsum(9)        total evaluation points
*                    
*    rpar    REAL                   the estimated parameters:
*                    rpar(1)         mean tatum period
*                    rpar(2)         tatum standard deviation
*                    rpar(3)         tactus mean period
*                    rpar(4)         tactus standard deviation
*                    rpar(5)         measure mean period
*                    rpar(6)         measure standard deviation
*                    rpar(7)         std dev. of logarithm of frequency
*                                       ratio from one step to the next
*                                       This probably depends on the
*                                       value chosen for step
*                    rpar(8)         std. dev. of remainder phase error
*                                     ie- how far observed phase differs
*                                     from what would be predicted based
*                                     on a steady tempo 
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      SUBROUTINE MKRPAR (rsum, rpar)
      
       IMPLICIT NONE
       REAL rpar
       DOUBLE PRECISION rsum
       DIMENSION rsum(9), rpar(8)
       
       DOUBLE PRECISION ave, dev, n
       
       n = rsum(9)
       ave = rsum(1) / n
       dev = SQRT(rsum(2) / n - ave **2)
       rpar(1) = SNGL(EXP(ave))      ! units of time
       rpar(2) = SNGL(dev)
       
       ave = rsum(3) / n
       dev = SQRT(rsum(4) / n - ave **2)
       rpar(3) = SNGL(EXP(ave))
       rpar(4) = SNGL(dev)
       
       ave = rsum(5) / n
       dev = SQRT(rsum(6) / n - ave **2)
       rpar(5) = SNGL(EXP(ave))
       rpar(6) = SNGL(dev)
         
*           multiply by 3 for the three levels counted
        n = n * 3.D0
        
        dev = SQRT(rsum(7) / n)        ! assume mean of zero
        rpar(7) = SNGL(dev)
        
        dev = SQRT(rsum(8) / n)
        rpar(8) = SNGL(dev)
        RETURN
       END  ! of MKRPAR


************************************************************************
* interpolate the g function
************************************************************************
      SUBROUTINE gfunc (X, Y, A, B, C, n)
       IMPLICIT NONE
       INTEGER i, j, n, ibeg, iend, ifault
       REAL X(n), Y(n), A(n), B(n), C(n), sr2pi
       REAL farey(23), width(23), height(23), wf, hf, tot
       REAL alpha, beta

************************************************************************
* These constants represent the strength at which an oscillator latches on
* to a signal, depending on the ratio of the signal frequency to the 
* oscillator frequency.  The array farey is made up of fractions from the
* Stern-Brocott-Farey sequence.  Width and height refer to an "empirical regime
* diagram" by Large, in which width represents a tolerance around an ideal
* integer ratio, and height is a parameter that controls how strongly
* the oscillator entrains to its input.  These variables can be reinterpreted
* as intensity (height), and standard deviation (width) around a mean (ratio).
* Refer to:
*     Edward Wilson Large
*     Dynamic Representation of Musical Structure
*     Ph.D. Thesis, Ohio State University, 1994, Fig. 39
************************************************************************      
       DATA farey /  0.0, .125, .1428571429, 0.1666666667, 0.2, 0.25, 
     $   0.2857142858, 0.3333333333, 0.375, 0.400, 0.4285714286, 0.5, 
     $   0.5714285714, 0.6, 0.625, 0.6666666667, 0.7142857143, 0.75, 
     $   0.8, 0.8333333333, 0.8571428571, 0.875, 1.000 /

       DATA width / 18, 2, 3, 4, 7, 10, 3, 12, 3, 7, 3, 19, 4, 7, 2, 13,
     $   3, 10, 7, 5, 3, 2, 21 /

       DATA height / 5, 2, 2, 2.2, 2.2, 3.3, 1.9, 3.3, 1.7, 2, 1.7, 4, 
     $   1.6, 1.8, 1.6, 3.5, 1.4, 3, 1.6, 1.5, 1.2, 1.2, 5 /
     
*       scale factors
       DATA wf, hf / 2500., 3.16 /
       DATA alpha, beta / 0., 0. /       ! not used
       
************************************************************************
* evaluate the function - sum of bell curves 
* with different means and deviations
************************************************************************
       sr2pi = SQRT(2. * ACOS(-1.0))
       DO j = 1, 23
         width(j) = width(j) / wf
         height(j) = height(j) / hf
       END DO
       
       tot = 0.
       DO i = 1, n
         X(i) = FLOAT(i-1) / (n-1)
         Y(i) = 23. / FLOAT(n)
         DO j = 1, 23
           Y(i) = Y(i) + (height(j)**2 / sr2pi) * 
     $        EXP( - (X(i)-farey(j))**2 / (2 * width(j))**2 )
         END DO
         tot = tot + Y(i)
       END DO
       DO i = 1, n
         Y(i) = LOG(Y(i) / tot)
       END DO
         
*-----------------------------------------------------------------------
*     call the spline interpolator
*-----------------------------------------------------------------------
       ibeg = 0 
       iend = 0 
       ifault = 0
       CALL CBSPL(X, Y, A, B, C, N, ibeg, iend, alpha, beta, ifault)
       IF (ifault .NE. 0) PRINT *, 'cbspl: error #', ifault
       RETURN 
      END  ! of gfunc
*++++++++++++++++++++++++ End of file rhythm.f +++++++++++++++++++++++++
