*-----------------------------------------------------------------------
* harmony.f - key and chord music analysis
*      by Andy Allinger, 2011-2015
*      andy_a@users.sourceforge.net
*
*   This program is free software, released to the public domain.  
*   It may be used by any person for any purpose.
*-----------------------------------------------------------------------

*-----------------------------------------------------------------------
* segment - split a note list into minimal partitions
*-----------------------------------------------------------------------
*
*                                 I/O LIST
*_Name______________Type________In/Out____Description___________________
* note_time[n]      REAL        In        times of note on/off events
*                                           (assumed already sorted)
* pitch[n]          INTEGER     In        MIDI note #
* offon[n]          INTEGER     In        0=OFF, 1=ON
* n                 INTEGER     In        number of note events
* pall              INTEGER     Out       # of minimal partitions
* partit[n]         REAL        Out       times of partitions
* seg[0:11,0:n]     INTEGER     Out       notes each pitch class, each segment
* cumseg[0:11,0:n]  INTEGER     Out       cumulative seg[]
* fault             INTEGER     Out       error code, 0=success
*
*-----------------------------------------------------------------------
      SUBROUTINE SEGMENT (note_time, pitch, offon, n, 
     &  pall, partit, seg, cumseg, fault)
       
       IMPLICIT NONE
       INTEGER pitch, offon, n, pall, seg, cumseg, fault
       REAL note_time, partit
       DIMENSION note_time(n), pitch(n), offon(n), partit(n), 
     &  seg(0:11,n), cumseg(0:11,n)
        
*          local variables
       INTEGER 
     &        i, j, k,                         ! counters
     &        npc,                             ! pitch % 12
     &        state(0:11)                      ! pitch vector
     
       REAL  
     &       pip,                               ! very brief duration
     &       t, tprev                           ! event times
       
      PARAMETER (pip = 0.035)                   ! negligible time
                               
*-----------------------------------------------------------------------
*         initialize, validate
*-----------------------------------------------------------------------
       fault = 0                 ! clear error codes
       
       IF (n < 2) THEN
         fault = 1               ! at least one NoteOff and NoteOn
         RETURN
       END IF
       
*-----------------------------------------------------------------------
* Find minimal segments -initialize states, times, and partition times
*       seg[i]        the pitch vector for the segment beginning at partit[i]
*       cumseg[i]     the accumulation of seg up to BUT NOT INCLUDING seg[i]
*       partit[1]     time of the first note on
*       partit[pall]  time of the last note off
* In the case of non-overlapping monophony, n = pall
*-----------------------------------------------------------------------
       DO j = 0, 11 
         state(j) = 0 
       END DO
       
*-----------------------------------------------------------------------
*          consider each NoteOff or NoteOn as a possible partition point
*-----------------------------------------------------------------------
       tprev = note_time(1)        ! indifferent to beginning silence
       t = tprev                     ! avoid compiler warning
       i = 0
       DO j = 1, n
         npc = MOD(pitch(j), 12)
         t = note_time(j)
         IF (t - tprev > pip) THEN                        ! different time?
           i = i + 1
           DO k = 0, 11 
             seg(k, i) = state(k) 
           END DO
           partit(i) = tprev
           tprev = t
         END IF
         IF (offon(j) .EQ. 1) THEN
           state(npc) = state(npc) + 1               ! increment pitch vector
         ELSE IF (offon(j) .EQ. 0) THEN
           state(npc) = state(npc) - 1               ! decrement pitch vector
         ELSE 
           fault = 2                  ! bad input
           RETURN
         END IF
       END DO

*          last note off (bogus silent segment)
       i = i + 1
       pall = i
       partit(i) = t
       DO j = 0, 11 
         seg(j,i) = state(j)
         IF (state(j) .NE. 0) PRINT *, 'segment: note sounding at end!!'
       END DO

*-----------------------------------------------------------------------
*       Accumulate the pitch vectors.  
*-----------------------------------------------------------------------
       DO j = 0, 11
         cumseg(j,1) = 0
       END DO
       
       DO i = 2, pall              ! over all minimal segments
         DO j = 0, 11
           cumseg(j,i) = cumseg(j,i-1) + seg(j,i-1)
         END DO
       END DO
       RETURN
      END ! of SEGMENT


*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
* keyal - find key changes using the chi-squared test and cross-entropy
*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
*
*                         TRANSFER LIST
*__Name________________Type_____I/O___Description_______________________
*  partit[n]           REAL     In    times of minimal partitions
*  cumseg[0:11,n]      INTEGER  In    cumulative pitch vector
*  n                   INTEGER  In    length of partit and cumseg
*  key_t[n]            REAL     Out   times of key changes
*  sharps[n]           INTEGER  Out   # sharps in key signature
*  majmin[n]           INTEGER  Out   major or minor key
*  k                   INTEGER  Out   # key changes detected
*  plist[n]            INTEGER  Work  0=Not a partition, -1=is a partition
*  cumw[n]             INTEGER  Work  sum of cumseg each partition
*  uprofile[0:11,0:1]  REAL     In    custom pitch profile
*  usecfg              LOGICAL  In    user has set a custom profile
*  fault               INTEGER  Out   0=no errors
*
*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
      SUBROUTINE KEYAL (partit, cumseg, n, key_t, sharps, majmin, k,
     & plist, cumw, uprofile, usecfg, fault)
       IMPLICIT NONE
       INTEGER cumseg, n, sharps, majmin, k, plist, cumw, fault
       REAL partit, key_t, uprofile
       LOGICAL usecfg
       DIMENSION partit(n), cumseg(0:11,n), key_t(n), sharps(n),
     &  majmin(n), plist(n), cumw(n), uprofile(0:11,0:1)
     
*           local variables 
       INTEGER 
     &         ba(0:11,2),       ! before-after contingency table
     &         i0, i, j,                   ! counters
     &         iearly(32), ilate(32),     ! indices of partition points
     &         L, m,            ! count NPC's (Neutral Pitch Classes), modes
     &         Ldir,            ! -1 descending tree, +1 ascending tree
     &         root, scale,     ! detected root and scale
     &         wmin,            ! limit for performing chi-square test
     &         wt(2)            ! cumw before and after proposed split
     
       REAL  
     &       default,       ! small probability
     &       e,             ! expected contingency cell value
     &       rmachine,      ! external function
     &       obs(0:11),     ! real-valued pitch vector
     &       profile(0:11,0:11,0:1),     ! key-finding template
     &       s, sbest,      ! entropy
     &       tprev, 
     &       wtot,          ! total of contingency table entries
     &       xcrit,         ! cutoff value for chi-square test
     &       x2, x2big      ! chi-square statistic
     
       LOGICAL first

       SAVE first, profile
       PARAMETER (xcrit = 33.0,       ! p > .9995 with 11 degrees of freedom
     &            wmin = 60,          ! 5 objects x 12 cells
     &            default = 0.000 001) ! 1 millionth

       DATA first / .TRUE. /

*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
*        The Kostka-Payne key profiles
*
*     Refer to:
*         David Temperley
*         Music and Probability, ch. 6:  "A Polyphonic Key-Finding Model"
*         MIT Press, 2007        
*
*          profile[npc, root, mode]
*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
       DATA profile( 0,0,0)  / 0.748 /
       DATA profile( 1,0,0)  / 0.060 /
       DATA profile( 2,0,0)  / 0.488 /
       DATA profile( 3,0,0)  / 0.082 /
       DATA profile( 4,0,0)  / 0.670 / 
       DATA profile( 5,0,0)  / 0.460 /
       DATA profile( 6,0,0)  / 0.096 /
       DATA profile( 7,0,0)  / 0.715 /
       DATA profile( 8,0,0)  / 0.104 /
       DATA profile( 9,0,0)  / 0.366 /
       DATA profile(10,0,0)  / 0.057 /
       DATA profile(11,0,0)  / 0.400 /

       DATA profile( 0,0,1)  / 0.712 /
       DATA profile( 1,0,1)  / 0.084 /
       DATA profile( 2,0,1)  / 0.474 /
       DATA profile( 3,0,1)  / 0.618 /
       DATA profile( 4,0,1)  / 0.049 /
       DATA profile( 5,0,1)  / 0.460 /
       DATA profile( 6,0,1)  / 0.105 /
       DATA profile( 7,0,1)  / 0.747 /
       DATA profile( 8,0,1)  / 0.404 /
       DATA profile( 9,0,1)  / 0.067 /
       DATA profile(10,0,1)  / 0.133 /
       DATA profile(11,0,1)  / 0.330 /
        
*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
*         begin
*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
       IF (n < 2) THEN
         fault = -1 
         RETURN
       ELSE
         fault = 0
       END IF
       
*        prepare key profiles
       IF (first) THEN
         first = .FALSE.
         IF (usecfg) THEN     ! copy user's data
           DO m = 0, 1
             DO i = 0, 11
               IF (uprofile(i,m) > 0.) THEN
                 profile(i,0,m) = uprofile(i,m)
               ELSE
                 profile(i,0,m) = default
               END IF
             END DO
           END DO
         END IF
         
*               log probability
         DO m = 0, 1
           DO i = 0, 11
             profile(i,0,m) = LOG(profile(i,0,m))
           END DO
         END DO        

*             transpose profiles
         DO m = 0, 1             ! major, minor
           DO j = 1, 11          ! roots 1...11
             DO i = 0, 11        ! each pitch class
               profile(i,j,m) = profile(MOD(i+11,12), j-1, m)
             END DO
           END DO
         END DO
       END IF  ! first ?
       
*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
*           initial values
*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
       DO i = 1, n
         L = 0
         DO j = 0, 11
           L = L + cumseg(j,i)
         END DO
         cumw(i) = L
         plist(i) = 0
       END DO
           
*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
*               recursively partition
*XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
       L = 1                 ! at trunk of tree       
       Ldir = +1
       iearly(L) = 1
       ilate(L) = n
       plist(1) = -1
       plist(n) = -1
       
       DO WHILE (L > 0)             ! until back down from tree
         IF (Ldir < 0) THEN          ! already considered this area ?
           IF (iearly(L) .EQ. iearly(L+1)) THEN    ! came from the early side?
             iearly(L+1) = ilate(L+1)              ! do the late sub-region
             ilate(L+1) = ilate(L)
             L = L + 1                    ! ascend
             Ldir = +1
           ELSE                         ! came from late sub-region
             L = L - 1              ! continue slide
           END IF
         ELSE               ! find optimal split
           x2big = 0.
           i0 = 1     ! avoid compiler warning
           DO 40 i = iearly(L)+1, ilate(L)-1
             wt(1) = cumw(i) - cumw(iearly(L))
             wt(2) = cumw(ilate(L)) - cumw(i)
             wtot = FLOAT(wt(1) + wt(2))
             IF (wt(1) < wmin .OR. wt(2) < wmin) GO TO 40    ! can't apply test
             DO m = 1, 2                 ! clear contingency table
               DO j = 0, 11
                 ba(j,m) = 0
               END DO
             END DO
             DO j = 0, 11              ! count
               ba(j,1) = cumseg(j,i) - cumseg(j,iearly(L))
             END DO
             DO j = 0, 11
               ba(j,2) = cumseg(j,ilate(L)) - cumseg(j,i)
             END DO
           
*              chi-square
             x2 = 0.
             DO m = 1, 2
               DO j = 0, 11
                 e = FLOAT(wt(m) * (ba(j,1) + ba(j,2))) / wtot ! expected count
                 x2 = x2 + (ba(j,m) - e) **2 / e
               END DO
             END DO
             IF (x2 > x2big) THEN 
               x2big = x2          ! maximum chi-square
               i0 = i                ! best split point
             END IF
  40       END DO ! next i
          
           IF (x2big > xcrit) THEN      ! split?
             plist(i0) = -1               ! true
             iearly(L+1) = iearly(L)      ! consider subdividing this part also
             ilate(L+1) = i0
             L = L + 1                      ! ascend
           ELSE                         ! leave intact
             Ldir = -1
             L = L - 1                   ! descend
           END IF 
         END IF ! ascend or descend ?
       END DO ! while loop
         
*-----------------------------------------------------------------------
*  Simple key-finding based on the Krumhansl-Schmuckler algorithm.
*  Cross-entropy is used instead of the correlation coefficient.       
*-----------------------------------------------------------------------
       i0 = 1
       k = 0
       tprev = 0.
       DO i = 2, n
         IF (plist(i) .EQ. 0) GO TO 60 ! not a partition
           
*             get an observation vector
         DO j = 0, 11 
           obs(j) = FLOAT(cumseg(j,i) - cumseg(j,i0))
         END DO
      
*-----------------------------------------------------------------------
*     measure cross-entropy      S = -SUM (No[i] @ ln(Pm[i])
*          where     No   observed frequency
*                    Pm   model probability
*-----------------------------------------------------------------------
         sbest = rmachine(3)        ! very large number
         root = 0                   ! avoid uninitialized warnings
         scale = 0
         DO m = 0, 1                   ! major, minor
           DO L = 0, 11                  ! each root
             s = 0.             
             DO j = 0, 11                   ! each pitch class
               s = s - obs(j) * profile(j,L,m)
             END DO
             IF (s < sbest) THEN
               sbest = s
               root = L
               scale = m
             END IF
           END DO
        END DO
      
*      convert root to number of sharps in key signature
         j = MOD(7 * root + 5, 12) - 5
         IF (k > 0) THEN      ! check for same key
           IF (j .EQ. sharps(k) .AND. scale .EQ. majmin(k)) GO TO 60
         END IF
         
         k = k + 1
         key_t(k) = tprev
         sharps(k) = j
         majmin(k) = scale

         i0 = i
         tprev = partit(i)
  60   END DO    ! next i

       RETURN
      END   ! of KEYAL
 

*-----------------------------------------------------------------------
* chordal - detect chord changes
*
* based on the algorithm published in:
*     Bryan Pardo & William P. Birmingham, 
*     "The Chordal Analysis of Tonal Music" 
*     University of Michigan technical report, 28 March 2001
*
* and:
*     Bryan Pardo & William P. Birmingham
*     "Algorithms for Chordal Analysis"
*     Computer Music Journal, v. 26 n. 2 pp. 27-49, Summer 2002
*
*-----------------------------------------------------------------------
*
*                                 I/O LIST
*_Name______________Type________In/Out____Description___________________
* partit[n]         Real        In        times of minimal partitions
* seg[0:11,n]       Integer     In        pitch vector each partition
* cumseg[0:11,n]    Integer     In        cumulative seg[]
* n                 Integer     In        length of partit, seg, and cumseg
* key_time[k]       Real        In        times of key change events
*                                           (assumed already sorted)
* sharps[k]         Integer     In        # sharps in key signature {-7...7}
*                                           (matches MIDI representation)
* majmin[k]         Integer     In        0=major key, 1=minor key
* k                 Integer     In        number of key change events
* job               Character   In        {L, T, E}
*                                          L => label chords
*                                          T => label and train
*                                          E => label and evaluate
* c                 Integer     Out       # of chord segments detected
* ch_t[n]           Real        Out       times of chord changes [seconds]
* q[n]              Integer     Out       "qualities" of chords
*                                           (as in the templates, see below)
* r[n]              Integer     Out       roots of chords {0...11}
*                                           as pitch % 12
* p                 Real        Out       log probability of chord sequence
* ccarr[num_q,0:11,0:11,0:1]
*                   Real        In        log probability of chord given chord
*                                          in Roman numeral notation.  The 
*                                          table is given twice, once
*                                          for major keys, and again for 
*                                          minor keys.  num_q is the number 
*                                          of chord qualities.  A separate 
*                                          routine is needed to define ccarr[].
* prox[0:6]         Real        Out       circle-of-fifth proximity of roots
* qual[num_q, 0:1]  Real        Out       count of each quality, each mode
* root[0:11, 0:1]   Real        Out       count of each root, each mode
* fault             Integer     Out       error code, 0=no errors
*
*-----------------------------------------------------------------------
      SUBROUTINE CHORDAL (partit, seg, cumseg, n, key_time, sharps,
     & majmin, k, job, c, ch_t, q, r, p, ccarr, prox, qual, root, fault)
       
       IMPLICIT NONE
       INTEGER num_q 
       PARAMETER (num_q = 15)        ! how many chord qualities
       
       INTEGER seg, cumseg, n, sharps, majmin, k, c, q, r, fault
       REAL partit, key_time, ch_t, p, ccarr, prox, qual, root
       CHARACTER job
       DIMENSION partit(n), seg(0:11,n), cumseg(0:11,n), 
     &  key_time(k), sharps(k), majmin(k), ch_t(n), q(n), r(n), 
     &  ccarr(num_q, 0:11, 0:11, 0:1), 
     &  prox(0:6), qual(num_q, 0:1), root(0:11, 0:1)
        
*          local variables
       INTEGER 
     &        C5(0:23),                        ! half-steps to circle of fifths
     &        hsteps(0:23),                    ! shortcut for % 12
     &        h, i, j,                         ! counters
     &        prev, pprev,                     ! index of prev seg, past prev
     &        qdel, qkeep, qold, qnew,         ! qualities
     &        rdel, rkeep, rold, rnew,         ! roots of chords
     &        romold, romnew,                  ! Roman Numeral version of roots
     &        rpp,                             ! an earlier root value
     &        scale,                           ! 0=major or 1=minor
     &        skip, 
     &        stot(-7:7,0:1),           ! convert sharps to tonic -12...-1
     &        template(0:11, num_q),           ! chord templates
     &        tm12,                            ! tonic - 12
     &        weight(0:11)                     ! pitch vector to analyze
     
       REAL  
     &       sdel, skeep, sprev                 ! possible segment scores
       
       LOGICAL 
     &         ALMEQ,                  ! an almost-equal function
     &         shifted                 ! whether templates have been transposed
       
       SAVE hsteps, C5, shifted, stot, template
       
*-----------------------------------------------------------------------
*       chord templates:  bits 0...11 of an integer represent pitches as
*         MIDI note number % 12
*      If a note is present in the chord, the corresponding bit is set
*      in the template.  Template data is included for a root note of C
*      and transposed to the other keys.  For example, C Major has the notes
*      C, E, G and the template {0, 4, 7}, which as an integer is:
*        2^0 + 2^4 + 2^7 = 145
* 
*      The order of the array is significant, lower numbered qualities
*      have precedence in case of ties during template-matching.
*    
* It is possible to redefine which templates are included.  There is no 
* authority as to what is or is not a real chord.  This list should be 
* adequate for most purposes.
*-----------------------------------------------------------------------
*                           integer  name                     template   symbol

       DATA template(0, 1) / 145/  ! major                    {0,4,7}    C
       DATA template(0, 2) / 137/  ! minor                    {0,3,7}    Cm
       DATA template(0, 3) /1169/  ! seventh                  {0,4,7,10} C7 
       DATA template(0, 4) /2193/  ! major seventh            {0,4,7,11} Cmaj7
       DATA template(0, 5) /1161/  ! minor seventh            {0,3,7,10} Cm7
       DATA template(0, 6) / 657/  ! sixth                    {0,4,7,9}  C6
       DATA template(0, 7) / 649/  ! minor sixth              {0,3,7,9}  Cm6
       DATA template(0, 8) /  73/  ! diminished               {0,3,6}    Cdim
       DATA template(0, 9) / 273/  ! augmented                {0,4,8}    Caug
       DATA template(0,10) / 585/  ! diminished seventh       {0,3,6,9}  Cdim7
       DATA template(0,11) /1105/  ! seventh flat fifth       {0,4,6,10} C7-5 
       DATA template(0,12) /1185/  ! seventh suspended fourth {0,5,7,10} C7sus4
       DATA template(0,13) /1173/  ! ninth                {0,2,4,7,10}   C9
       DATA template(0,14) /1205/  ! eleventh             {0,2,4,5,7,10} C11
       DATA template(0,15) /1685/  ! thirteenth           {0,2,4,7,9,10} C13
 
       DATA shifted /.FALSE./
       DATA C5 / 0,5,2,3,4,1,6,1,4,3,2,5,0,5,2,3,4,1,6,1,4,3,2,5 /
       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
                               
*-----------------------------------------------------------------------
*         initialize, validate
*-----------------------------------------------------------------------
       fault = 0                 ! clear error codes
       
*         transpose chord templates into all keys
       IF (.NOT. shifted) THEN
         DO i = 1, num_q
           DO j = 1, 11
             template(j,i) = ISHFTC(template(j-1,i), 1, 12)
           END DO
         END DO
         shifted = .TRUE.
       END IF
       
*            check for valid job
       IF (job .NE. 'L' .AND. job .NE. 'E' .AND. job .NE. 'T') THEN
         fault = 3
         RETURN
       END IF       

*-----------------------------------------------------------------------
*          HarmAn search.  Consider each partition point in order, as 
* a possible point for a chord change, based on whether keeping the 
* partition or deleting it yields a better match to a chord template. 
*-----------------------------------------------------------------------
       h = 1
       tm12 = stot(sharps(h), majmin(h))
       prev = 0
       DO j = 0, 11 
         weight(j) = seg(j,1)
       END DO
       CALL score_seg (weight, template, tm12, prev, 0, 
     &        q(1), r(1), sprev)
       
*-----------------------------------------------------------------------
*         main loop               
*-----------------------------------------------------------------------
       pprev = 0
       prev = 1
       rpp = 0
       DO i = 2, n-1

         IF (h < k) THEN                       ! change key
           IF (partit(i) .GE. key_time(h+1)) THEN
             h = h + 1
             tm12 = stot(sharps(h), majmin(h))
           END IF
         END IF
       
*-----------------------------------------------------------------------
*             score segment i
*-----------------------------------------------------------------------
         DO j = 0, 11  
           weight(j) = seg(j,i) 
         END DO
         CALL score_seg (weight, template, tm12, prev, r(prev), 
     &        qkeep, rkeep, skeep)
         
*-----------------------------------------------------------------------
*            score segment formed by deleting partition i
*-----------------------------------------------------------------------
         DO j = 0, 11  
           weight(j) = cumseg(j,i+1) - cumseg(j,prev) 
         END DO
         IF (pprev > 0) rpp = r(pprev)
         CALL score_seg (weight, template, tm12, pprev, rpp, 
     &        qdel, rdel, sdel)
         
*-----------------------------------------------------------------------
*             keep or delete partition point ?
*-----------------------------------------------------------------------
         IF (skeep + sprev > sdel) THEN       ! keep partition
           q(i) = qkeep
           r(i) = rkeep
           sprev = skeep
           pprev = prev
           prev = i
         ELSE                                 ! delete partition
           q(i) = -1                           ! delete labeling
           r(i) = -1                           !  "      "
           q(prev) = qdel
           r(prev) = rdel
           sprev = sdel
         END IF
       END DO

*-----------------------------------------------------------------------
*        remove values of q[] and r[] at deleted partition points
*-----------------------------------------------------------------------
       c = 1
       qold = q(1) 
       rold = r(1)
       DO 30 i = 2, n-1
         IF (q(i) < 0) GO TO 30           ! not a partition point
         qnew = q(i)
         rnew = r(i)
         IF ( qold .EQ. qnew .AND. rold .EQ. rnew ) GO TO 30 ! skip duplicates
         
         c = c + 1                         ! copy data to lower position
         ch_t(c) = partit(i)
         q(c) = qnew
         r(c) = rnew
         
         qold = qnew
         rold = rnew
  30   END DO ! next i
       
       IF (job .EQ. 'L') RETURN
       IF (job .EQ. 'T') GO TO 70    ! count chord change probabilities

*-----------------------------------------------------------------------
*   If JOB is set to evaluate, find probability of the sequence of chords, 
*   figured as a Markov chain of events E[1]...E[n] :        
*         P(E) = P(E[n]) * PROD ( P(E[i-1] | E[i]) )           2...i...n
*-----------------------------------------------------------------------
       p = 0.
       h = 1
       rold = r(1)
       scale = majmin(1)
       tm12 = stot(sharps(1), scale)
       romold = hsteps(r(1) - tm12)
       
       DO i = 2, c
         IF (h < k) THEN                       ! change key
           IF (ch_t(i) .GE. key_time(h+1)) THEN
             h = h + 1
             scale = majmin(h)
             tm12 = stot(sharps(h), scale)
           END IF
         END IF
       
         qnew = q(i)
         rnew = r(i)
         romnew = hsteps(rnew - tm12)
         
*               Markov chain of log-probability
         p = p + ccarr(qnew, romnew, romold, scale)  
         rold = rnew
         romold = romnew
       END DO
       RETURN       ! Done evaluating probability
             
*-----------------------------------------------------------------------       
  70   CONTINUE  ! count chord proximity, quality and root
*-----------------------------------------------------------------------
       h = 1
       rold = r(1)
       scale = majmin(1)
       tm12 = stot(sharps(1), scale)
       romold = hsteps(r(1) - tm12)
       skip = 0
       DO i = 2, c
         IF (h < k) THEN                       ! change key
           IF (ch_t(i) .GE. key_time(h+1)) THEN
             h = h + 1
             scale = majmin(h)
             tm12 = stot(sharps(h),scale)
             IF (ALMEQ(ch_t(i), key_time(h))) THEN
               skip = 1          ! skip one transition on a neat key change
             ELSE
               skip = 2        ! skip two transitions on a mid-chord key change
             END IF
           END IF
         END IF
         qnew = q(i) 
         rnew = r(i)
         romnew = hsteps(rnew - tm12)
         
         IF (skip .LE. 0) THEN
           j = C5(romnew - romold + 12)       ! proximity
           prox(j) = prox(j) + 1.
           qual(qnew,scale) = qual(qnew,scale) + 1.       ! quality
           root(romnew,scale) = root(romnew,scale) + 1.   ! root
         END IF
        
         rold = rnew
         romold = romnew 
         skip = skip - 1
       END DO
       
       RETURN 
      END ! of CHORDAL
 
       
*-----------------------------------------------------------------------
* score_seg - score a segment
*
*__Name__________________Type______In/Out___Description_________________
*  weight[0:11]          Integer   In       observed pitch vector
*  template[0:11,num_q]  Integer   In       chord template data
*  tm12                  Integer   In       tonic key as NPC, minus 12 (!)
*  prev                  Integer   In       index of previous segment
*  rprev                 Integer   In       root chosen at previous segment
*  qbest                 Integer   Out      quality this segment
*  rbest                 Integer   Out      root this segment
*  sbest                 Real      Out      score this segment labeling
*                                            bigger numbers are better
*-----------------------------------------------------------------------
      SUBROUTINE score_seg (weight, template, tm12, prev, rprev,
     &                           qbest, rbest, sbest)
       IMPLICIT NONE
       INTEGER num_q 
       PARAMETER (num_q = 15)
       INTEGER weight, template, tm12, prev, rprev, qbest, rbest
       REAL sbest
       DIMENSION weight(0:11), template(0:11,num_q)
       
*        local variables       
       INTEGER C5(0:23),                ! half steps to circle of fifths dist.
     &         ctPP, ctPA, ctAP, ctAA,  ! contingency table entries
     &         h, i, j,                 ! count qualities, roots, pitches
     &         o, n,                    ! temporary integers
     &         wsum                     ! total of weight[]
       REAL rmachine, score
       
       SAVE C5
       DATA C5 / 0,5,2,3,4,1,6,1,4,3,2,5,0,5,2,3,4,1,6,1,4,3,2,5 /

*-----------------------------------------------------------------------
       wsum = 0
       DO j = 0, 11
         wsum = wsum + weight(j)
       END DO
       IF (wsum .EQ. 0) THEN
         qbest = 1
         rbest = 0
         sbest = 0.
         RETURN
       END IF
         
       sbest = -rmachine(3)                ! real-valued scores
       DO h = 1, num_q                     ! for each quality
         DO i = 0, 11                         ! for each root
           score = 0.
           ctPP = 0
           ctPA = 0
           ctAP = 0
           ctAA = 0
           
           DO j = 0, 11                            ! for each pitch class
             IF ( BTEST(template(i,h), j) ) THEN
               IF (weight(j) > 0) THEN                 
                 ctPP = ctPP + weight(j)     ! in template
               ELSE  
!!                 ctPA = ctPA + wsum    ! don't know why, but it helps...
                 ctPA = ctPA + 1
               END IF
             ELSE                            ! not in template
               IF (weight(j) > 0) THEN
                 ctAP = ctAP + weight(j)
               ELSE
!!                 ctAA = ctAA + wsum    ! ...to add wsum rather than 1 !
                 ctAA = ctAA + 1
               END IF
             END IF
           END DO
           
*-----------------------------------------------------------------------
*            contingency table formula 
*  Don't know why it's the fourth root of the denominator, but it helps
*-----------------------------------------------------------------------
           score = FLOAT(ctPP * ctAA - ctPA * ctAP) * SQRT(FLOAT(wsum))
     &            / SQRT(SQRT(FLOAT(
     & (ctPP + ctPA) * (ctPA + ctAA) * (ctAA + ctAP) * (ctAP + ctPP) )))
                      
           IF (score - sbest) 50, 30, 40

  30       CONTINUE     ! tiebreakers

*                    prefer more root pitch
           IF (weight(i) > weight(rbest)) THEN    
             qbest = h
             rbest = i
             sbest = score
             GO TO 50
           ELSE IF (weight(i) < weight(rbest)) THEN
             GO TO 50
           END IF

*                prefer clearer harmonies
           IF  (h < qbest) THEN                 
             qbest = h
             rbest = i
             sbest = score
             GO TO 50
           ELSE IF (h > qbest) THEN
             GO TO 50
           END IF

*                     nearest root on the circle of fifths
           IF (prev > 0) THEN
             o = C5(rbest - rprev + 12)     ! old
             n = C5(i - rprev + 12)         ! new
             IF (n < o) THEN
               qbest = h 
               rbest = i
               sbest = score
               GO TO 50
             ELSE IF (n > o) THEN
               GO TO 50
             END IF
           END IF
                        
*                    prefer roots near to the tonic on the circle of fifths
           o = C5(rbest - tm12)     ! old
           n = C5(i - tm12)         ! new
           IF (n < o) THEN
             qbest = h
             rbest = i
             sbest = score
             GO TO 50
           ELSE IF (n > o) THEN
             GO TO 50
           END IF

*                resolve aug, dim7, 7b5 based on prior chord
*                 This rule is intended to make Bdim7 follow C
           IF ( (h .EQ. 9) .OR. (h .EQ. 10) .OR. (h .EQ. 11)
     &                    .AND. (prev > 0) ) THEN
             o = MOD(rprev - rbest + 12, 12)         ! half steps DESCENDED
             n = MOD(rprev - i + 12, 12)
             IF (n < o) THEN
               qbest = h
               rbest = i
               sbest = score
               GO TO 50
             ELSE IF (n > o) THEN
               GO TO 50
             END IF
           END IF  ! dim7 resolution

!           PRINT *, 'chordal: warning - unbreakable tie: ', 
!     &         h, i, ' vs ', qbest, rbest
           GO TO 50
 
  40       CONTINUE        ! new preferred chord
           qbest = h
           rbest = i
           sbest = score

  50     END DO ! next i
       END DO  ! next h

       RETURN
      END ! of score_seg
 

************************************************************************
*        finish making the user custom chord profile
************************************************************************
      SUBROUTINE MKPQR (prox, qual, root)
       IMPLICIT NONE
       INTEGER num_q
       PARAMETER (num_q = 15)
       REAL prox, qual, root
       DIMENSION prox(0:6), qual(num_q,0:1), root(0:11,0:1)

       INTEGER i, L
       REAL seventh, twelfth, tot
       LOGICAL safdiv           ! safe-to-divide function
       SAVE seventh, twelfth
       DATA seventh, twelfth / 0.14285714, 0.083333333 /
       
*            proximity
       tot = 0.
       DO i = 0, 6
         tot = tot + prox(i)
       END DO
       DO i = 0, 6
         IF (SAFDIV(prox(i), tot)) THEN
           prox(i) = prox(i) / tot
         ELSE
           prox(i) = seventh
         END IF
       END DO

*             quality
       DO L = 0, 1
         tot = 0.
         DO i = 1, num_q
           tot = tot + qual(i,L)
         END DO
         DO i = 1, num_q
           IF (SAFDIV(qual(i,L), tot)) THEN
             qual(i,L) = qual(i,L) / tot
           ELSE
             qual(i,L) = twelfth
           END IF
         END DO
       END DO

*              root
       DO L = 0, 1
         tot = 0.
         DO i = 0, 11
           tot = tot + root(i,L)
         END DO
         DO i = 0, 11
           IF (SAFDIV(root(i,L), tot)) THEN
             root(i,L) = root(i,L) / tot
           ELSE
             root(i,L) = twelfth
           END IF
         END DO
       END DO
       
       RETURN
      END ! of MKPQR

        
*-----------------------------------------------------------------------
* Make a chord change probability matrix
* The chord data here has been collected from 10 tunes, 5 in major keys 
* and 5 in minor keys, selected at whim from:
*    Songs of the Dragon
*    WNGGA, 1993.
*
*                          TRANSFER LIST
*_Name______________Type_____I/O__Description______________________
* uprox[0:6]        Real     In   probability circle-of-fifths root proximity
* uqual[num_q, 0:1] Real     In   probability each chord quality, each mode
* uroot[0:11, 0:1]  Real     In   probability of Roman numeral root, each mode
* usecfg            Logical  In   user supplies values for prox, qual, root
*
* ccarr[num_q, 0:11, 0:11, 0:1]
*                   Real     Out  log probability each successive quality,
*                                     each succesive (Roman numeral) root;
*                                     given each preceding root,
*                                     for major and minor modes
*-----------------------------------------------------------------------
      SUBROUTINE MKCCARR (uprox, uqual, uroot, usecfg, ccarr)
       IMPLICIT NONE
       INTEGER num_q 
       PARAMETER (num_q = 15)
       REAL ccarr, uqual, uprox, uroot
       LOGICAL usecfg
       DIMENSION ccarr(num_q, 0:11, 0:11, 0:1),
     &  uprox(0:6), uqual(num_q, 0:1), uroot(0:11, 0:1)
       
*        local variables
       INTEGER i, j, k, m,                   ! counters
     &         C5(0:23),                     ! circle of fifths distance
     &         dif

       REAL 
     &      default,                        ! arbitrary small probability
     &      prox(0:6),                      ! prob. of root proximity 
     &      qual(num_q, 0:1),               ! prob. of each quality, each mode
     &      root(0:11, 0:1)                 ! prob. of each root, each mode
     
       DOUBLE PRECISION tot
       
       SAVE C5, prox, qual, root
       DATA C5 / 0,5,2,3,4,1,6,1,4,3,2,5,0,5,2,3,4,1,6,1,4,3,2,5 /
       
*  Casual inspection showed no difference between major and minor modes,
*   so the data was pooled.  A Cauchy distribution with average 1.161 
*   and dispersion .471 extrapolated the rare/unobserved values.  
       DATA prox / .102, .644, .173, .044, .019, .011, .007 /
       
*    The quality probabilities also take into account chords noted in 
*  the Kostka-Payne data set, as reported by Pardo & Birmingham.
*  This array must be arranged to match the chord templates!
       DATA qual / .600, .056, .222, .003, .003, .003, .035, .019, ! Major 
     &             .003, .041, .003, .003, .003, .003, .003, 
     &             .186, .526, .134, .003, .003, .003, .068, .019, 
     &             .003, .041, .003, .003, .003, .003, .003 /

*         probability of each root, nominal value for unobserved roots
       DATA root / 0.387, 0.001, 0.114, 0.001, 0.001, 0.114,       ! major
     &              0.001, 0.334, 0.001, 0.030, 0.015, 0.001,
     &              0.441, 0.001, 0.015, 0.052, 0.001, 0.150,       ! minor
     &              0.001, 0.269, 0.007, 0.001, 0.060, 0.001 /

       DATA default / 0.000 001 /             ! one in a million
            
*-----------------------------------------------------------------------
*         multiply
*-----------------------------------------------------------------------
       IF (usecfg)  THEN           ! use defaults
         DO i = 0, 6 
           prox(i) = uprox(i) 
         END DO
         DO j = 0, 1
           DO i = 1, num_q
             qual(i,j) = uqual(i,j)
           END DO
         END DO
         DO j = 0, 1
           DO i = 0, 11
             root(i,j) = uroot(i,j)
           END DO
         END DO
       END IF ! default chord data?
           
       DO m = 0, 1             ! major and minor modes
         DO k = 0, 11            ! each prior root
           DO j = 0, 11            ! each present root
             dif = C5(ABS(j - k))
             DO i = 1, num_q         ! each present quality
               ccarr(i,j,k,m) = prox(dif) * qual(i,m) * root(j,m)
             END DO
           END DO
         END DO
       END DO
       
*-----------------------------------------------------------------------
*        normalize
*-----------------------------------------------------------------------
       DO m = 0, 1                     ! each mode
         DO k = 0, 11                    ! each prior root

*                add           
           tot = 0.                        ! sum over priors
           DO j = 0, 11                      ! each present root
             DO i = 1, num_q                   ! each present quality
               tot = tot + ccarr(i,j,k,m) 
             END DO
           END DO
           IF (tot .LE. 0.) tot = default 
             
*              divide and take logarithm
           DO j = 0, 11
             DO i = 1, num_q
               ccarr(i,j,k,m) = SNGL(DBLE(ccarr(i,j,k,m)) / tot)
               IF (ccarr(i,j,k,m) > default) THEN
                 ccarr(i,j,k,m) = LOG(ccarr(i,j,k,m))
               ELSE
                 ccarr(i,j,k,m) = LOG(default)
               END IF
             END DO
           END DO
             
         END DO
       END DO
       
       RETURN
      END ! of mkccarr           
*+++++++++++++++++++++++ End of file harmony.f +++++++++++++++++++++++++ 
