************************************************************************
* melody.f - routines for analyzing melodies
*              by Andy Allinger, 2012-2017, public domain
*
*    refer to:
*      David Temperley
*      Music and Probability, ch. 3-6
*      MIT Press, 2007
*
************************************************************************
*                             I/O LIST
*___Name___________Type______In/Out___Description_______________________
*   timnot(n)      Real      In       Times of NoteOn events
*   pitch(n)       Integer   In       Pitches of notes as MIDI #
*   n              Integer   In       Number of note events
*   key_time(k)    Real      In       Times of key changes
*   sharps(k)      Integer   In       # sharps in the key signature
*   majmin(k)      Integer   In       Major or minor key
*   k              Integer   In       # key changes
*   beat(blim,L)   Real      In       Times of beats in various strata
*   blim           Integer   In       Row bound of beat
*   pkarr(0:127, 0:127, 0:11, 0:1)   
*                  Real      In       Log probability of a pitch, given
*                                      previous pitch, its scalar degree,
*                                      and a major or minor scalar mode
*   p              Real      Out      Log probability
*   fault          Integer   Out      0 = no errors
************************************************************************
      SUBROUTINE melody (timnot, pitch, n, key_time, sharps, majmin,
     &  k, beat, blim, pkarr, p, fault)
       
       IMPLICIT NONE
       INTEGER L                   ! # of rhythmic strata
       PARAMETER (L = 5)
       INTEGER pitch, n, sharps, majmin, k, blim, fault
       REAL timnot, key_time, beat, pkarr, p
       DIMENSION timnot(n), pitch(n), key_time(k), sharps(k), 
     &  majmin(k), beat(blim,L), pkarr(0:127, 0:127, 0:11, 0:1)

*                constants
       INTEGER TL
       PARAMETER (TL = 3)                ! which stratum is the tactus level

*        local variables
       INTEGER
     &    d,                     ! scale degree
     &    h, i, j,               ! counters
     &    missed,                ! beats without notes
     &    m, mold,               ! MIDI note number
     &    phiold(L), phinew(L),  ! phase of prev, current beat that had a note
     &    scale,                 ! 0=major, 1=minor
     &    stot(-7:7,0:1),        ! translate key signature to tonic-12
     &    tonic                  ! root of current key
     
       REAL 
     &      e,                    ! remainder phase error
     &      phase,                ! phase of a note at a level
     &      lnprof(0:11,0:1),     ! Logarithms of Essen key profiles
     &      thresh                ! threshhold for a note coinciding with beat
     
       LOGICAL found,             ! note occurs off any beat
     &         ISINF              ! is infinity external function
      
       SAVE lnprof, stot, thresh

       DATA lnprof / -1.693, -6.908, -1.864, -5.809, -1.655, -2.216,  ! major
     &               -5.298, -1.542, -6.908, -2.551, -5.521, -2.900,
     &               -1.650, -5.298, -1.904, -1.720, -6.215, -1.938,  ! minor
     &               -6.215, -1.604, -3.270, -4.423, -2.937, -3.817 /

       DATA stot / 11, 6, 1, 8, 3, 10, 5, 0, 7, 2, 9, 4, 11, 6, 1,    ! major
     &              8, 3, 10, 5, 0, 7, 2, 9, 4, 11, 6, 1, 8, 3, 10 /  ! minor

       DATA thresh / 0.05 /

************************************************************************
*           validate
************************************************************************
       fault = 0
       p = 0.
       IF (n < 1) THEN
         fault = 1
         RETURN
       ELSE IF (n < 2) THEN
         RETURN
       END IF
       
************************************************************************
*           beats
************************************************************************
*         phase of first note
       DO i = 1, L 
         phinew(i) = NINT(beat(1,i))
       END DO
       DO j = 1, n                        ! for each event
         found = .FALSE.

*    log probability of a note on a beat, LN(P(N|B)), combined with
*   log probability of the beats in which no note occurred, LN(1-P(N|B))
         DO i = L, 1, -1                    ! for each stratum
           phase = beat(j,i)
!!!           PRINT *, 'melody: beat(', j, ',', i, ') is ', beat(j,i)
           e = ABS(phase - ANINT(phase))      ! remainder phase error
           IF (e < thresh) THEN
             phiold(i) = phinew(i)
             phinew(i) = NINT(phase)
             missed = DIM(phinew(i) - phiold(i), 1)  ! beats w/o a note
             IF (ISINF(REAL(missed))) PRINT *, 'melody: missed is INF!'
             IF (i > TL) THEN                ! strophe level
               IF (.NOT. found) p = p - 0.0513 
               p = p - missed * 2.996
!!!!               PRINT *, 'dec2996rementing p by: ', missed * 2.996
             ELSE IF (i .EQ. TL) THEN        ! tactus level
               IF (.NOT. found) p = p - 0.301
               p = p - missed * 1.347
!!!               PRINT *, 'dec1347rementing p by: ', missed * 1.347
             ELSE                            ! tatum level
               IF (.NOT. found) p = p - 0.968
               p = p - missed * 0.478
!!!               PRINT *, 'dec478rementing p by: ', missed * .478
             END IF
             found = .TRUE.
           END IF
         END DO ! loop i

*              not a beat at any level
         IF (.NOT. found) p = p - 4.605
!!!         PRINT *, 'dec4605rementing p by ', 4.605
       END DO ! loop j

************************************************************************
*          pitches
************************************************************************
       h = 1
       scale = majmin(1)
       tonic = stot(sharps(1), scale)
       m = pitch(1)
       
*         P(N[1])
       d = MOD( (m - tonic + 12), 12)
       p = p + lnprof(d, scale)
       IF (ISINF(p)) PRINT *, 'Melody: INF adding in lnprof!!'
       
*         PRODUCT( P(N[j+1] | N[j])
       DO j = 2, n
         IF (h < k) THEN              ! key change
           IF (timnot(j) .GE. key_time(h+1)) THEN
             h = h + 1
             scale = majmin(h)
             tonic = stot(sharps(h), scale)
           END IF
         END IF
         mold = m
         m = pitch(j)
         p = p + pkarr(m, mold, tonic, scale)
!!!         PRINT *, 'decrementing p by ', pkarr(m,mold,tonic,scale)
         IF (ISINF(P)) PRINT *, 'melody:  INF adding in pkarr!!!'
       END DO          ! it only looks easy 
       RETURN
      END ! of melody



************************************************************************
*                initialize an array of constants
*    uprofile[0:11, 0:1]              user's custom key profile
*    uvar                             user's custom variance of pitch proximity
*    usecfg                           override defaults
*    pkarr                            log probabilities of a pitch transition
************************************************************************
      SUBROUTINE MKPKARR (uprofile, uvar, usecfg, pkarr)
       IMPLICIT NONE
       REAL uprofile, uvar, pkarr              ! the proximity-key array
       LOGICAL usecfg
       DIMENSION uprofile(0:11, 0:1), pkarr(0:127, 0:127, 0:11, 0:1)

       INTEGER i, j, k, L
       REAL default, dif, p, profile(0:11, 0:1), twovar
       DOUBLE PRECISION tot
       LOGICAL ISINF               ! external function

*        The key-profiles derived from the Essen folk music collection
       SAVE profile, twovar  ! the routine shouldn't have to be called twice
       DATA profile / 0.184, 0.001, 0.155, 0.003, 0.191, 0.109, ! major keys
     &                0.005, 0.214, 0.001, 0.078, 0.004, 0.055, 
     &                0.192, 0.005, 0.149, 0.179, 0.002, 0.144, ! minor keys
     &                0.002, 0.201, 0.038, 0.012, 0.053, 0.022 /
       DATA twovar / 14.4 /      ! twice the variance
       
************************************************************************
       default = 0.000 001
!       default = LOG(0.000 001)
!       PRINT *, 'using default=', default
       
*              don't use the Essen profiles?
       IF (usecfg) THEN
          DO L = 0, 1
            DO k = 0, 11
              profile(k,L) = uprofile(k,L)
            END DO
          END DO
          twovar = 2 * uvar
        END IF
        
*              multiply
       DO L = 0, 1               ! major, minor
         DO k = 0, 11             ! each key
           DO j = 0, 127           ! prior pitch
             tot = 0.
             DO i = 0, 127
               dif = ABS(i - j)
               p = EXP(-(dif **2) / twovar)          ! proximity
               p = p * profile(MOD(i-k+12, 12), L)              ! key
               pkarr(i,j,k,L) = p
               tot = tot + p
             END DO  ! next i
             IF (tot .LE. 0.D0) tot = default
             
*                 make probabilities sum to 1
             DO i = 0, 127
               pkarr(i,j,k,L) = SNGL(DBLE(pkarr(i,j,k,L)) / tot)
               IF (pkarr(i,j,k,L) > default) THEN
                 pkarr(i,j,k,L) = LOG(pkarr(i,j,k,L))
               ELSE
                 pkarr(i,j,k,L) = LOG(default)
               END IF
               
!!!!!!!!!!!!!!!!!!!!!!!! debugging
               IF (pkarr(i,j,k,L) .NE. pkarr(i,j,k,L)) 
     &           PRINT *, 'mkpkarr: NaN! i,j,k,L: ', i, j, k, L
               IF (ISINF(pkarr(i,j,k,L)))
     &           PRINT *, 'mkpkarr: INF! i,j,k,L: ', i, j, k, L
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     
     
             END DO  ! loop i
!             PRINT  *, 'k=', k, ' total in column ', j, ' is ', tot
           END DO  ! next j             

         END DO ! next k
       END DO  ! next L
       
!!!!            dump result
!       DO L = 0, 1
!         DO k = 0, 11
!           DO j = 0, 127
!             DO i = 0, 127
!               IF (k .EQ. 2) THEN
!        PRINT *, 'at:', i, j, k, L, ' %12:', MOD(i,12), pkarr(i,j,k,L)
!               END IF
!             END DO
!           END DO
!         END DO
!       END DO
      
       RETURN       
      END  ! of MKPKARR 
             

************************************************************************
*         Add Data To Profile
************************************************************************
      SUBROUTINE ADTP (timnot, pitch, n, key_time, sharps,
     & majmin, k, profile, ms, ms2)
     
       IMPLICIT NONE
       INTEGER pitch, n, sharps, majmin, k
       REAL timnot, key_time, profile
       DOUBLE PRECISION ms, ms2
       DIMENSION timnot(n), pitch(n), key_time(k), sharps(k), 
     &  majmin(k), profile(0:11, 0:1)
     
*          local
       INTEGER d, dif, h, hsteps(0:23), j, m, mold,
     &    scale, skip, stot(-7:7, 0:1), tm12
       
       SAVE hsteps, stot
       DATA hsteps /0,1,2,3,4,5,6,7,8,9,10,11,0,1,2,3,4,5,6,7,8,9,10,11/
       DATA stot / -1,-6,-11,-4,-9,-2,-7,-12,-5,-10,-3,-8,-1,-6,-11, ! major
     &             -4,-9,-2,-7,-12,-5,-10,-3,-8,-1,-6,-11,-4,-9,-2 / ! minor
       
************************************************************************
       h = 1
       scale = majmin(1)
       tm12 = stot(sharps(1),scale)
       m = pitch(1)

       d = hsteps(MOD(m, 12) - tm12)
       profile(d,scale) = profile(d,scale) + 1.
       
       skip = 0
       DO j = 2, n
         IF (h < k) THEN
           IF (timnot(j) .GE. key_time(h+1)) THEN
             h = h + 1
             scale = majmin(h)
             tm12 = stot(sharps(h),scale)
             skip = 1
           END IF
         END IF
         
*           scale degree frequency
         mold = m
         m = pitch(j)
         IF (skip < 1) THEN           ! omit pitches straddling a key change
           d = hsteps(MOD(m, 12) - tm12)
           profile(d,scale) = profile(d,scale) + 1.
         ELSE
           skip = 0
         END IF

*            proximity variance
         dif = m - mold
         ms = ms + dif
         ms2 = ms2 + dif **2
       END DO ! loop j
       RETURN
      END  ! of ADTP



************************************************************************
*       finish making the profile
************************************************************************
      SUBROUTINE MKPROF (profile, var, ms, ms2)
       IMPLICIT NONE
       REAL profile, var
       DOUBLE PRECISION ms, ms2
       DIMENSION profile(0:11, 0:1)

       INTEGER i, L
       REAL default, neach, ntotal
       DOUBLE PRECISION nn
       DATA default / 0.083333333 /
       
       ntotal = 0.
       DO L = 0, 1
         neach = 0.
         DO i = 0, 11
           neach = neach + profile(i,L)
         END DO
         DO i = 0, 11
           IF (neach > 0.) THEN
             profile(i,L) = profile(i,L) / neach
           ELSE
             profile(i,L) = default
           END IF
         END DO
         ntotal = ntotal + neach
       END DO
       
       nn = DBLE(ntotal)
       var = SNGL(ms2 / nn - (ms / nn) **2)
       RETURN
      END  ! of MKPROF

*++++++++++++++++++++++++ End of file melody.f +++++++++++++++++++++++++
