*======================================================================= 
* affine.f - search for self-affinity in a musical sequence 
* by Andy Allinger, 2012, released to the public domain 
* This program may be used by any person for any purpose. 
*
* refer to:
*     Dimitrios Rafailidis and Yannis Manolopoulos
*     "The Power of Music:  Seaching for Power-Laws in Symbolic Musical Data"
*     draft report, 2006.
*
*======================================================================= 
*                               I/O LIST 
* Name          Type       In/Out    Description 
* note_time[m]  REAL       in        times of note onsets 
* pitch[m]      INTEGER    in        MIDI note numbers 
* m             INTEGER    in        how many notes in the sequence 
* key_time[k]   REAL       in        times of key change events 
* sharps[k]     INTEGER    in        MIDI key signature 
* majmin[k]     INTEGER    in        major or minor keys 
* k             INTEGER    in        # of key changes 
* beat[m]       REAL       in        phase of each note event 
* n             INTEGER    in        number of beats
* aff           REAL       out       estimated self-affinity
* fault         INTEGER    out       error code 0 = no errors
*
*======================================================================= 
      SUBROUTINE affine(note_time, pitch, m, key_time, sharps, majmin,
     &     k, beat, n, aff, fault)
     
       IMPLICIT NONE
       INTEGER pitch, m, sharps, majmin, k, n, fault
       REAL note_time, key_time, beat, aff
       DIMENSION note_time(m), pitch(m), key_time(k), sharps(k), 
     &   majmin(k), beat(m)
     
*              constants
        INTEGER L
        PARAMETER (L = 5)      ! length of a segment to compare
      
*            local variables       
       INTEGER 
     &    D(m),                            ! pitch, in scale degrees
     &    E(n),                            ! pitch, at each beat
     &    G(0:11, 0:1),                    ! translate note # to Do Re Mi
     &    h, i, j, o,                      ! counters
     &    ind(n-L+1),                      ! permutation vector for IORDER
     &    JSW(2*(n-L+1)),                  ! workspace for radix sort
     &    p,                               ! ocurrences of an identical segment
     &    q,                               ! number of distinct segments
     &    scale,                           ! a value of majmin[]
     &    stot(-7:7, 0:1),                 ! convert keysig to tonic pitch
     &    tm12,                            ! a value of stot[], tonic - 12
     &    X(n-L+1, L),                     ! matrix of segments
     &    xlim                             ! # bits, largest entry in X[]
     
       REAL
     &    ave,                          ! average of log population
     &    b,                            ! time of a beat
     &    BIG,                          ! very large number
     &    dif,                          ! time from an event to a beat
     &    idzdz,                        ! external function inverse(Z'/Z)
     &    rmachine,                     ! external function
     &    near                          ! shortest value of dif      
     
       DOUBLE PRECISION sum             ! sum of log populations
      
       SAVE G, stot
       DATA G / 0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12,          ! major
     &          0, 1, 2, 4, 5, 6, 7, 8, 10, 11, 12, 13 /        ! minor
       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
 
*======================================================================= 
*        begin
*======================================================================= 
       fault = 0
       BIG = rmachine(3)
       
       IF (m < 10) THEN      ! never mind the short stuff
         aff = 6.     ! no affinity on scale 1...6
!!!         fault = -1
         RETURN
       END IF
 
*======================================================================= 
*         transpose
*======================================================================= 
       h = 1
       scale = majmin(1)
       tm12 = stot(sharps(1),scale)
       
*        make invariant to tranpositions within the scale and between keys
       DO i = 1, m
         IF (h < k) THEN               ! change key
           IF (key_time(h+1) .LE. note_time(i)) THEN
             h = h + 1
             scale = majmin(h)
             tm12 = stot(sharps(h), scale)
           END IF
         END IF
         j = pitch(i) - tm12        ! Roman numeral form
         D(i) = 14 * (j / 12) + G( MOD(j, 12), scale ) 
       END DO

*======================================================================= 
*       difference
*======================================================================= 
       xlim = 1
       h = D(m)
       DO i = m, 2, -1
         j = h
         h = D(i-1)
         o = j - h
         D(i) = o
         xlim = MAX(xlim, ABS(o))
       END DO
       D(1) = 0
         
*======================================================================= 
*       quantize, so inserting or deleting brief notes should have no effect
*======================================================================= 
       h = 1
       i = 1
       b = beat(1)
       sum = 0.
       DO j = 1, n                         ! for each beat
         near = BIG                   
  10     IF (NINT(b) .LE. j) THEN          ! is beat j the nearest to event i ?
           dif = ABS(b - j)
           IF (dif < near) THEN            ! is event i the nearest to beat j ?
             near = dif
             h = i
           END IF
           IF (i < m) THEN                  ! more events ?
             i = i + 1
             b = beat(i)
             GO TO 10
           END IF
         END IF
         E(j) = D(h)
       END DO ! next j
 
*======================================================================= 
*             look for repeated patterns
*======================================================================= 
*           put segments into matrix
       DO i = 1, L                  ! for each pitch in segment
         DO h = 1, n-L+1               ! for each starting offset
           X(h,i) = E(i+h-1)
         END DO
       END DO

*             bits to represent largest entry in array
       i = BIT_SIZE(xlim)
       DO WHILE (.NOT. BTEST(xlim, i-1))
         i = i - 1
       END DO
       xlim = i

*          sort each column to bring matches together
       DO i = 1, L
         DO h = 1, n-L+1
           ind(h) = h
         END DO
         CALL RSORTI(X(1,i), ind, n-L+1, -xlim, JSW)
         DO h = 1, L
           IF (h .EQ. i) CYCLE     ! already sorted
           CALL IORDER(X(1,h), ind, n-L+1)
         END DO
       END DO
         
*            compare ajacent rows
       sum = 0.
       p = 1
       q = 1
       DO i = 2, n-L+1
         DO h = 1, L
           IF (X(i,h) .NE. X(i-1,h)) GO TO 20
         END DO
         p = p + 1           ! same
         CYCLE

  20     CONTINUE              ! different
         q = q + 1    
         IF (p .NE. 1) sum = sum + LOG(FLOAT(p))
         p = 1
       END DO
       IF (p .NE. 1) sum = sum + LOG(FLOAT(p))

*           the Zipf exponent
       ave = SNGL(sum / DBLE(q))
       IF (ave < 0.0223) THEN
         aff = 5.
       ELSE IF (ave > 10 000.) THEN
         aff = 1.
       ELSE
         aff = idzdz(-ave)
       END IF
         
      END ! of affine


         
*ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
*   inverse of the derivative of the Riemann zeta function divided 
* by the Riemann zeta function, for use in calculating the Zipf distribution 
*
* an approximation of an approximation of an approximation!, by fitting a 
* spline to the CERN MATHLIB zeta function, differentiating the spline, 
* dividing, swapping X and Y, and fitting a quartic rational algebraic 
* function.
* If you know a better way to invert Z'/Z, use it.  
*
*     Recommended domain:   -2500 < X < -0.025           error < 5%  ?
*
*ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
      FUNCTION IDZDZ(X)
        
       IMPLICIT NONE
       REAL IDZDZ, X
       
       REAL C(10), numer, denom
       SAVE C
       DATA C / 0.2078623159E-03, 0.2160934794, 16.62735147, 
     &   -31.33392134, 5.937427891, 0.2078622858E-03, 0.2163012693, 
     &    16.84372205, -14.37591128, 1.000000000  /
      
*        synthetic division
       numer = (((C(1) * X + C(2)) * X + C(3)) * X + C(4)) * X + C(5)
       denom = (((C(6) * X + C(7)) * X + C(8)) * X + C(9)) * X + C(10)
       idzdz = numer / denom
       
      END ! of IDZDZ
*++++++++++++++++++++++++ End of file affine.f +++++++++++++++++++++++++
