*-----------------------------------------------------------------------
* puntal.f - contrapuntal analysis
* by Andy Allinger, 2012, released to the public domain
* This program may be used by any person for any purpose.
*-----------------------------------------------------------------------

*-----------------------------------------------------------------------
*  edist - distance between notes, such that they belong to different voices
*
* refer to :
*      Wai Man Szeto and Man Hon Wong
*      "Stream Segregation Algorithm for Pattern Matching
*         in Polyphonic Music Databases,"
*      draft report, 2003.
*
*___Name_______Type______In/Out___Description_________________
*   TON(N)     Real      In       Note ontimes
*   TOFF(N)    Real      In       Note offtimes
*   PITCH(N)   Integer   In       MIDI note #
*   N          Integer   In       Length of arrays
*   I          Integer   In       Index of first point
*   J          Integer   In       Index of second point
*   EDIST      Real      Out      Distance
*-----------------------------------------------------------------------
      FUNCTION EDIST (TON, TOFF, PITCH, N, I, J)
       IMPLICIT NONE
       INTEGER PITCH, N, I, J
       REAL TON, TOFF, EDIST
       DIMENSION TON(N), TOFF(N), PITCH(N)

       REAL BETA
       PARAMETER (BETA = 400.)  ! dif from A=57 to B=59 is worth 0.08 seconds

       REAL TOT, F, FMELA, FMELB
       REAL DIF1, DIF2, DIF3

*-----------------------------------------------------------------------
*          Time separation distance
*-----------------------------------------------------------------------
       DIF1 = MIN(ABS(TON(I) - TOFF(J)), ABS(TON(J) - TOFF(I)))

*          Overlap penalty
       DIF2 = DIM(MIN(TOFF(I), TOFF(J)), MAX(TON(I), TON(J)))

*          Pitch difference (use the Mel scale?)
*         Convert MIDI note number to Hertz
*                f = 440 @ 2^((n-69) / 12)
       F = 440. * 2. ** ((PITCH(I) - 69) / 12.)

*          Convert Hertz to Mel
       FMELA = 2595. * LOG10(1. + F / 700.)
       F = 440. * 2. ** ((PITCH(J) - 69) / 12.)
       FMELB = 2595. * LOG10(1. + F / 700.)
       DIF3 = ABS(FMELA - FMELB)

*          Make units comparable
       DIF3 = DIF3 / BETA

*          Combine
       TOT = DIF1 **2 + DIF2 **2 + DIF3 **2
       EDIST = SQRT(TOT)
       RETURN
      END  ! of EDIST



*-----------------------------------------------------------------------
* nextcomb - lists the next combination of n choose k in lexical order
*-----------------------------------------------------------------------
      SUBROUTINE NEXTCOMB (N, K, R, DONE)
       IMPLICIT NONE
       INTEGER N, K, R                ! I/O list
       LOGICAL DONE
       DIMENSION R(K)
       INTEGER C, I, J, P             ! local variables

*       Begin
       IF (K .EQ. 0) THEN
         DONE = .TRUE.
         RETURN
       END IF
       I = K
       C = N - K

*       Check incremented value p at position i
  10   P = R(I) + 1
       IF (P > C + I) THEN       ! carry
         I = I - 1
         IF (I < 1) THEN
           DONE = .TRUE.
           RETURN
         END IF
         GO TO 10
       END IF

*       Apply increment
       DO J = I, K            ! roll over digits
         R(J) = P
         P = P + 1
       END DO
       RETURN
      END  ! of NEXTCOMB



*-----------------------------------------------------------------------
* puntal - assign stream ID's to a set of note events
*  Refer to:
*    Elaine Chew and Xiaodan Wu
*    "Separating Voices in Polyphonic Music:  A Contig Mapping Approach"
*    Computer Music Modeling and Retrieval, pp. 1-20, 2005
*
*  see further references below
*
*___Name_______________Type__________I/O_______Description______________
*   NOTE_TIME(N)       Real          In        Note on times
*   PITCH(N)           Integer       In        Note numbers
*   OFFON(N)           Integer       In        NoteOn or NoteOff
*   N                  Integer       In        Number of note on/off events
*   V(N)               Integer       Out       Voice membership of note events
*   MSV                Integer       Out       Number of voices made
*   B(0:VLIM,N)        Integer       Neither   Backpointers, reused for voice
*   NIV(N+1)           Integer       Neither   Implicit voices each segment
*   NOV(N+1)           Integer       Neither   Observed voices each segment
*   NOTIND(PLIM,N+1)   Integer       Neither   What notes are in each sement
*   TOFF(N)            Real          Neither   Note-off time each note-on
*   PLIST(N+1)         Logical       Neither   Whether a partition is kept
*   IFAULT             Integer       Out       Error code 0 = OK
*-----------------------------------------------------------------------
      SUBROUTINE PUNTAL (NOTE_TIME, PITCH, OFFON, N, V, MSV, B,
     $  NIV, NOV, NOTIND, TOFF, PLIST, IFAULT)

       IMPLICIT NONE
       INTEGER PLIM, VLIM
       PARAMETER (PLIM = 24,   ! polyphony requirement in General MIDI Level 1
     $            VLIM = 10)   ! fingers
       INTEGER PITCH, OFFON, N, V, MSV, B, NIV, NOV, NOTIND, IFAULT
       REAL NOTE_TIME, TOFF
       LOGICAL PLIST
       DIMENSION NOTE_TIME(N), PITCH(N), OFFON(N), V(N), B(0:VLIM,N),
     $  NIV(N+1), NOV(N+1), NOTIND(PLIM,N+1), TOFF(N), PLIST(N+1)

*          Local variables
       INTEGER
     $         BACK,                   ! a back pointer to best analysis
     $         CBEST, C,               ! integer-valued connection scores
     $         DIF,
     $         H, I, J, K, L,          ! counters
     $         INEW, IOLD,             ! whether data is in even or odd column
     $         M, O,                   ! temp storage for niv(L), nov(L)
     $         N1, N2,                 ! note indices
     $         NSV,                    ! number of sounding voices
     $         NOTDEL                  ! index of a note event to delete
       INTEGER
     $         P,                      ! # partitions needing combining
     $         PALL,                   ! # of minimal partitions
     $         PSI(0:127),             ! pointer to sound inception
     $         PVECT(VLIM),            ! temporary storage
     $         R(VLIM),                ! combination vector, lexical order
     $         RTRY(VLIM),             ! trial combination
     $         REPTAB(VLIM),           ! replacement table
     $         XI, XJ                  ! indices into note data

       REAL
     $      BIG,            ! very large number
     $      COBSF, COBSM,   ! penalty for mismatch of observed vs expected
     $      CTRANS,         ! penalty for changing number of voices
     $      D(PLIM,PLIM),   ! distances between simultaneous notes
     $      EDIST,          ! external distance function
     $      RMACHINE,       ! external function
     $      PIP,            ! minimal time interval
     $      TNEW,           ! time of current event
     $      TPREV,          ! time of previous event
     $      U(0:VLIM,0:1),  ! dynamic programming table
     $      UBEST, UTRY     ! log-probabilities (arbitrary units)

       LOGICAL ALLDEC,      ! all contigs have declining voice count
     $         DONE,
     $         USEDEC      ! merge partitions with declining voice count

       SAVE COBSF, COBSM, CTRANS, PIP
       DATA COBSF, COBSM, CTRANS / -0.2, -1.0, -7.0 /
       DATA PIP / 0.035 /

*-----------------------------------------------------------------------
*        Begin
*-----------------------------------------------------------------------
       IFAULT = 0
       BIG = RMACHINE(3)
       IF (N < 3) THEN       ! trivial cases
         DO J = 1, N
           V(J) = 0
         END DO
         MSV = 0
         RETURN
       END IF

*          Limit the number of voices to make
       IF (VLIM > PLIM) THEN
         PRINT *, 'voice: Too many voices requested.  Maximum ', PLIM
         IFAULT = 2
         RETURN
       END IF

*-----------------------------------------------------------------------
*           Segment
*-----------------------------------------------------------------------
       DO K = 0, 127
         PSI(K) = 0
       END DO
       MSV = 0
       NSV = 0
       TPREV = 0
       I = 0

       DO J = 1, N
         TNEW = NOTE_TIME(J)
         IF (TNEW - TPREV > PIP) THEN
           I = I + 1
           K = 0
           DO L = 0, 127
             IF (PSI(L) > 0) THEN
               K = K + 1
               IF (K .LE. PLIM) NOTIND(K,I) = PSI(L)
             END IF
           END DO
           NOV(I) = MIN(NSV, PLIM)
           TPREV = TNEW
         END IF
         IF (OFFON(J) .EQ. 1) THEN
           IF (PSI(PITCH(J)) .EQ. 0) THEN
             PSI(PITCH(J)) = J
             NSV = NSV + 1
           ELSE
!             PRINT *, 'event ', J, ' attempt to turn on sounding note!'
             CONTINUE
           END IF
         ELSE IF (OFFON(J) .EQ. 0) THEN
           IF (PSI(PITCH(J)) > 0) THEN
             TOFF(PSI(PITCH(J))) = TNEW
             PSI(PITCH(J)) = 0
             NSV = NSV - 1
           ELSE
!             PRINT *, 'event ', j, ' attempt to hush silent note!'
             CONTINUE
           END IF
         ELSE
           PRINT *, 'event ', J, ' bad offon value: ', OFFON(J)
         END IF
       END DO
       PALL = I

*           All notes briefer than a pip.  Mark them stray.
       IF (PALL .EQ. 1) THEN
         DO I = 1, N
            V(I) = 0
         END DO
         MSV = 0
         RETURN
       END IF

*-----------------------------------------------------------------------
*            Viterbi smoothing pass
*-----------------------------------------------------------------------
       DO J = 0, VLIM                    ! uninformative distribution
         U(J,0) = 0.
       END DO

*          Forwards phase, for each time step
       DO I = 1, PALL

*           Alternate storage in array U(o,0:1) on even and odd j
         IF (BTEST(I, 0)) THEN
           INEW = 1
           IOLD = 0
         ELSE
           INEW = 0
           IOLD = 1
         END IF

*           For each # voices at present time step
         DO J = 0, VLIM
           UBEST = -BIG
           BACK = 0                 ! avoid compiler warning

*            For each # voices at previous time step
           DO 10 K = 0, VLIM
             UTRY = U(K,IOLD)              ! prior probability
             IF (UTRY < UBEST) GO TO 10          ! bail out

*              Penalize change
             IF (J .NE. K) UTRY = UTRY + CTRANS

             IF (UTRY > UBEST) THEN
               UBEST = UTRY
               BACK = K
             END IF
  10       CONTINUE  ! next k

*               Observation likelihood
           DIF = DIM(J, NOV(I))                     ! too few observed notes
           UBEST = UBEST + COBSF * DIF
           DIF = DIM(NOV(I), J)                     ! too many observed notes
           U(J,INEW) = UBEST + COBSM * DIF

*              save pointer to most likely prior state
           B(J,I) = BACK

         END DO  ! next j
       END DO ! next i

*         Find the highest scoring path
       UBEST = -BIG
       BACK = 0                          ! avoid compiler warning
       IF (BTEST(PALL, 0)) THEN          ! avoid compiler warning
         INEW = 1
       ELSE
         INEW = 0
       END IF
       DO J = 0, VLIM
         IF (U(J,INEW) > UBEST) THEN
           BACK = J
           UBEST = U(J,INEW)
         END IF
       END DO

*         Follow back pointers
       NIV(PALL) = 0
       NOV(PALL) = 0
       DO I = PALL, 2, -1
         NIV(I) = MIN(BACK, NOV(I))
         MSV = MAX(MSV, NIV(I))
         BACK = B(BACK,I)
       END DO
       NIV(1) = MIN(BACK, NOV(1))
       MSV = MAX(MSV, NIV(1))                   ! End of Viterbi

*-----------------------------------------------------------------------
*           Reduce polyphony
*-----------------------------------------------------------------------
       DO L = 1, PALL
         M = NIV(L)
         O = NOV(L)

         IF (O > M) THEN             ! remove some notes from this partition
           DO J = 1, O-1
             XJ = NOTIND(J,L)
             IF (XJ < 1) PRINT*, 'notind=0 at j=', J, ' partit:', L
             DO I = J+1, O                  ! find all similarities
               XI = NOTIND(I,L)
               IF (XI < 1) PRINT*, 'notind=0 at i=', I, ' partit:', L
               D(I,J) = EDIST(NOTE_TIME, TOFF, PITCH, N, XI, XJ)
             END DO ! next i
           END DO ! next j

*            Try all combinations to find the maximum distance
           DO K = 1, M
             RTRY(K) = K
           END DO
           UBEST = 0.
           DONE = .FALSE.
           DO WHILE (.NOT. DONE)
             UTRY = 0.
             DO J = 1, M-1
               DO I = J+1, M
                 UTRY = UTRY + D(RTRY(I),RTRY(J))
               END DO
             END DO
             IF (UTRY > UBEST) THEN
               DO K = 1, M
                 R(K) = RTRY(K)
               END DO
               UBEST = UTRY
             END IF
             CALL NEXTCOMB(O, M, RTRY, DONE)
           END DO

*             Shorten note events not marked by combination vector R
           J = 1
           DO 30 K = 1, O
             IF (J .LE. M) THEN
               IF (K .EQ. R(J)) THEN         ! skip this note
                 J = J + 1
                 GO TO 30
               END IF
             END IF

             NOTDEL = NOTIND(K,L)
             DO I = L, PALL
               DONE = .TRUE.
               DO H = 1, NOV(I)
                 IF (NOTIND(H,I) .EQ. NOTDEL) THEN   ! swap to end
                   NOTIND(H,I) = NOTIND(NOV(I),I)
                   NOV(I) = NOV(I) - 1               ! forget existence
                   NIV(I) = MIN(NIV(I), NOV(I))      ! keep consistent
                   DONE = .FALSE.
                   GO TO 20
                 END IF
               END DO  ! next h
  20           IF (DONE) GO TO 30             ! no trace of this event
             END DO  ! next i
  30       CONTINUE ! next k
         END IF ! too many events?!

*         Sort remaining events by frequency in ascending order
         DO K = 1, M
           PVECT(K) = PITCH(NOTIND(K,L))
         END DO
         CALL QSORTI(PVECT, R, M)
         CALL IORDER(NOTIND(1,L), R, M)

*         Label local voice numbers
         DO K = 1, M
           B(K,L) = M - K + 1
         END DO
         PLIST(L) = .TRUE. ! mark partition as not joined
       END DO ! next L
       PLIST(PALL) = .TRUE.

!            DEBUG:
        DO I = 1, PALL
          IF (NIV(I) .NE. NOV(I)) THEN
        PRINT *, 'at ', I, 'niv nov mismatch: ', NIV(I), NOV(I)
          END IF
        END DO

*-----------------------------------------------------------------------
*         Prioritized contig mapping
*  refer to:
*    Asako Ishigaki, Masaki Matsubara, Hiroaki Saito
*    "Prioritized Contig Combining to Segregate Voices in Polyphonic Music"
*    draft report, 2011
*-----------------------------------------------------------------------
       P = PALL - 2       ! remove all partitions except the beginning and end
       ALLDEC = .FALSE.

  41   IF (P .EQ. 0) GO TO 45

       IF (ALLDEC) THEN
         USEDEC = .TRUE.
       ELSE
         USEDEC = .FALSE.
       END IF
       ALLDEC = .TRUE.

       L = 1
  42   IF (L .EQ. PALL-1) GO TO 41             ! end of list?
       IF (PLIST(L+1)) THEN
         I = L + 1               ! find voice count at last partit of contig
         DO WHILE (.NOT. PLIST(I+1))
           I = I + 1
         END DO

         IF (NOV(L) .LE. NOV(I)) THEN     ! merge increasing
           ALLDEC = .FALSE.
           GO TO 44                         ! and merge
         ELSE IF (USEDEC) THEN
           USEDEC = .FALSE.
           ALLDEC = .FALSE.
           GO TO 44
         END IF
       END IF

  43   L = L + 1                     ! next partition
       GO TO 42

  44   IF (NOV(L) > NOV(I)) THEN    ! merge declining contigs
         H = I               ! partition with fewer notes
         I = L               ! partition with more notes
       ELSE                  ! merge increasing contigs
         H = L
       END IF

       M = NOV(H)
       O = NOV(I)
       DO K = 1, M
         RTRY(K) = K
       END DO

       CBEST = 32767
       DONE = .FALSE.
       DO WHILE (.NOT. DONE)       ! find connection costs for all combinations
         C = 0
         DO J = 1, M
           N1 = NOTIND(J,H)
           N2 = NOTIND(RTRY(J),I)
           IF (N1 .EQ. N2) THEN
             C = C - 1                      ! subtract for same note
           ELSE
             DIF = ABS(PITCH(N1) - PITCH(N2))
             C = C + DIF                    ! penalty is pitch difference
           END IF
         END DO ! next j
         IF (C < CBEST) THEN        ! lowest connection cost
           DO K = 1, M
             R(K) = RTRY(K)
           END DO
           CBEST = C
         END IF
         CALL NEXTCOMB(O, M, RTRY, DONE)
       END DO

*          Replace local note numbers with global note numbers
       DO J = 1, M
         REPTAB(B(J,H)) = B(R(J),I)     ! fill replacement table
       END DO

*          Renumber voices
       DO K = H, 1, -1
         DO J = 1, NOV(K)
           B(J,K) = REPTAB(B(J,K))
         END DO
         IF (PLIST(K)) GO TO 50
       END DO
  50   PLIST(L+1) = .FALSE.
       P = P - 1
       GO TO 43

*         Final results
  45   DO J = 1, N
         V(J) = 0
       END DO
       DO J = 1, PALL
         DO I = 1, NOV(J)
           K = B(I,J)
           L = NOTIND(I,J)
           IF (V(L) .EQ. 0) V(L) = K
         END DO
       END DO
       RETURN
      END ! of puntal
*+++++++++++++++++++++ End of file puntal.f ++++++++++++++++++++++++++++
