//
// Copyright:     Copyright 2003 Craig Stuart Sapp
// Programmer:    Craig Stuart Sapp <craig@ccrma.stanford.edu>
// Creation Date: Fri Feb 28 06:53:11 PST 2003
// Last Modified: Fri Feb 28 10:23:09 PST 2003
// Filename:      settheory.c
// Web Address:   http://peabody.sapp.org/class/dmp2/lab/settheory/settheory.c
// Syntax:        C; Max4/MSP2 External Object; CodeWarrior 6.0
// OS:            Mac OS 9; PPC
//
// Description:   Real-time Set Theory calculator inspired by Matthew McCabe's
//                SetCalc external object for Max:
//                   http://www.euph0r1a.net/projects/?handler=mumble
//
// Bugs to fix:   * prime form does not always reduce properly: 
//                  try D-A which should be 0 5.
//                * normal form and pitch class set display the last 
//                  pitch class played.  Would be nice to set them to zero 
//                  or clear them.
//                * pitch class count gets stuck notes for some reason, 
//                  even with overdrive.
//                * perhaps add Forte number analysis outlet
//

#include "ext.h"
#include <cstring>        /* included for memcpy() function */
#include <cstdlib>        /* included for qsort() function */

typedef struct {
   t_object max_data;     // Max/MSP data, MUST come first in struct

   long     invel;        // temporary storage for incoming note velocity
   long     pc[12];       // count of the pitch classes being played

   long     pcset[12];    // storage for pitch class set
   long     normal[12];   // storage for normal form of PC set
   long     prime[12];    // storage for prime form of PC set
   long     iclass[12];   // storage for interval class vector

   long     pcsetSize;    // number of elements in pcset array
   long     normalSize;   // number of elements in normal array
   long     primeSize;    // number of elements in prime array
   long     iclassSize;   // number of elements in interval class vector array

   void*    outputPccount; // output for Pitch Class note counts
   void*    outputPcset;   // output for Pitch Class sets
   void*    outputNormal;  // output for Normal form for PC set
   void*    outputPrime;   // output for Prime form of PC set
   void*    outputIclass;  // output for Interval Class Vector 

} MyObject;

void* object_data = NULL;

// function declarations:
void   main                (void);
void*  create_object       (void);
void   InputBang           (MyObject* mo);
void   InputKeyNumber      (MyObject* mo, long value);
void   InputVelocity       (MyObject* mo, long value);
void   MessageClear        (MyObject* mo);
void   clear               (MyObject* mo);
long   midilimit           (long value);
long   getIntervalClass    (long input);
int    longcompare         (const void* a, const void* b);
void   createOutputList    (t_atom *outputList, long* input, long size);
long   returnSmallest      (long *list, long count);
long   calculatePCSet      (long *output, long* input);
long   calculateNormalForm (long *output, long *input, long size);
int    calculatePrimeForm  (long *output, long *input, long size);
long   calculateICV        (long *output, long *input, long size);
void   doAnalysis          (MyObject* mo);


/////////////////////////////////////////////////////////////////////////
//
// Initialization functions
//

//////////////////////////////
//
// main -- called once when the object is created in a patcher window.
//

void main(void) {
   setup((t_messlist**)&object_data, (method)create_object,
         NULL, sizeof(MyObject), NULL, A_NOTHING);
   addbang((method)InputBang);
   addint ((method)InputKeyNumber);       // inlet 1
   addinx ((method)InputVelocity, 1);     // inlet 2
   addmess((method)MessageClear, "clear", A_NOTHING);  
}



//////////////////////////////
//
// create_object -- create the data storage for the mydiff object and
//     and setup input 1.
//

void* create_object(void) {
   MyObject* mo = (MyObject*)newobject(object_data);
   mo->invel    = 0;
   MessageClear(mo);
   mo->outputIclass  = listout(mo);   // outlet 5
   mo->outputPrime   = listout(mo);   // outlet 4
   mo->outputNormal  = listout(mo);   // outlet 3
   mo->outputPcset   = listout(mo);   // outlet 2
   mo->outputPccount = listout(mo);   // outlet 1
   intin(mo, 1);                      // inlet 2
   return mo;
}


/////////////////////////////////////////////////////////////////////////
//
// Behavior functions
//


//////////////////////////////
//
// InputKeyNumber -- what to do when a new note number comes into
//     the object.
//

void InputKeyNumber(MyObject* mo, long value) {
   long pitchclass;

   value = midilimit(value);
   pitchclass = value % 12;
   
   if (mo->invel == 0) {
      mo->pc[pitchclass]--;
      if (mo->pc[pitchclass] < 0) {
         mo->pc[pitchclass] = 0;
      }
   } else {
      mo->pc[pitchclass]++;
   }

   InputBang(mo); 
} 
   
   
   
//////////////////////////////
//
// InputBang -- do the analyses
//

void InputBang(MyObject *mo) {
   t_atom outputList[12];
   
   doAnalysis(mo);
   
   createOutputList(outputList, mo->iclass, mo->iclassSize);
   outlet_list(mo->outputIclass,  NULL, mo->iclassSize,  outputList);
     
   createOutputList(outputList, mo->prime, mo->primeSize == 0 ? 1 : mo->primeSize);
   outlet_list(mo->outputPrime,  NULL, mo->primeSize,  outputList);

   createOutputList(outputList, mo->normal, mo->normalSize == 0 ? 1 : mo->normalSize);
   outlet_list(mo->outputNormal, NULL, mo->normalSize, outputList);

   createOutputList(outputList, mo->pcset, mo->pcsetSize == 0 ? 1 : mo->pcsetSize);
   outlet_list(mo->outputPcset, NULL, mo->pcsetSize, outputList);
 
   createOutputList(outputList, mo->pc, 12);
   outlet_list(mo->outputPccount, NULL, 12, outputList);
}



//////////////////////////////
//
// doAnalysis -- do the set theory calculations to find the
//    normal form, prime form, and interval class.
//

void doAnalysis(MyObject* mo) {
   int count = 0;
   int i;
   for (i=0; i<12; i++) {
      count += mo->pc[i];
   }
   if (count > 0) {
      mo->pcsetSize  = calculatePCSet     (mo->pcset,  mo->pc);
      mo->normalSize = calculateNormalForm(mo->normal, mo->pcset,  mo->pcsetSize);
      mo->primeSize  = calculatePrimeForm (mo->prime,  mo->normal, mo->normalSize);
      mo->iclassSize = calculateICV       (mo->iclass, mo->prime,  mo->primeSize);
   } else {
      mo->pcsetSize  = 0;
      mo->normalSize = 0;
      mo->primeSize  = 0;
      mo->iclassSize = 0;
      
      mo->pcset[0]  = 0;
      mo->normal[0] = 0;
      mo->prime[0]  = 0;
      mo->iclass[0] = mo->iclass[1] = mo->iclass[2] = mo->iclass[3] = 
         mo->iclass[4] = mo->iclass[5] = 0;
   }
}



//////////////////////////////
//
// createOutputList -- store the integer arrays for output into the
//   outputList.
//

void createOutputList(t_atom *outputList, long* input, long size) {
   int i;
   for (i=0; i<size; i++) {
      outputList[i].a_type     = A_LONG; 
      outputList[i].a_w.w_long = input[i];
   }
}



//////////////////////////////
//
// InputVelocity -- what to do when a new velocity value comes into
//     the object.
//

void InputVelocity(MyObject* mo, long value) {
   mo->invel = value;
}



//////////////////////////////
//
// MessageClear -- clear the contents of the object.
// 

void MessageClear(MyObject* mo) {
     clear(mo);
     InputBang(mo);
}


/////////////////////////////////////////////////////////////////////////
//
// Non-interface functions:
//

void clear(MyObject* mo) {
   int i;
   for (i=0; i<12; i++) {
      mo->pc[i]     = 0;
      mo->pcset[i]  = 0;
      mo->normal[i] = 0;
      mo->prime[i]  = 0;
      mo->iclass[i] = 0;
   }
   mo->pcsetSize  = 0;
   mo->normalSize = 0;
   mo->primeSize  = 0;
   mo->iclassSize = 0;
}


//////////////////////////////
//
// midilimit -- limit a number to the range from 0 to 127.
//    if the input is less than 0, return 0.  
//    if the input is greater than 127, return 127.  
//

long midilimit(long value) {
   if (value < 0)     return   0;
   if (value > 127)   return 127;
   return value;
}



//////////////////////////////
//
// calculatePCSet -- given a list of twelve-tone states, generate
//   a list of the pitch classes contained in the array.  Returns
//   the number of elements in the calculated pitch set vector
//

long calculatePCSet(long *output, long* input) {
   long i;
   long j = 0;
   for (i=0; i<12; i++) {
      if (input[i] > 0) {
         output[j++] = i;
      }
   }
   return j;
}



//////////////////////////////
//
// calculateNormalForm -- given a input pitch class set, convert it to
//   normal form.  Returns the size of the output vector, which should
//   match the size of the input vector.
//

long calculateNormalForm(long *output, long *input, long size) {
   long matrix[12][12];
   long i,j,k,smallest;
   long min, minIndex = 0;
   long dist, first[12];
   long duplicates[12][12];
   long duplicateCount = 0;
   long rowCount = size;
   long colCount = size;

   //  create a matrix of all possible orderings of this set
   for (i=0; i<size; i++)  {
      for (j=0; j<size; j++)  {
         matrix[i][j] = input[j];
      }
      k = input[0];
      for(j=0; j<size-1; j++)  {
         input[j] = input[j+1];
      }
      input[size-1] = k;
   }

   //  algorithm for finding the normal form 
   while (1) {
      min = 12;
      for (i=0; i<rowCount; i++) {
         dist = matrix[i][colCount-1] - matrix[i][0];
         if (dist < 0) {
            dist += 12;
         }

         if (dist < min) { 
            duplicateCount = 0;
            min = dist;
            minIndex = i;
            memcpy(duplicates[duplicateCount++], matrix[i], sizeof(matrix[i]));
            continue;
         }
         if (dist == min) {
            memcpy(duplicates[duplicateCount++], matrix[i], sizeof(matrix[i]));
         }
      }

      if (duplicateCount == 0) {          // this means the array has length 1, 
         break;                          // so there is only 1 minimum
      }
      if (duplicateCount == size - 1)  { // symmetrical set 
         break;
      }

      rowCount = duplicateCount;
      duplicateCount = 0;
      j = sizeof(duplicates);
      memcpy(matrix, duplicates, j);
      colCount--;
   }

   if (duplicateCount != 0)  {
      // find smallest first value and that is the normal 
      for (i=0; i<size; i++)  {
         first[i] = matrix[i][0];
      }
      smallest = returnSmallest(first, size);
   } else {
      // minIndex is the normal
      smallest = minIndex;
   }

   // store the normal form from the correct row
   for (i=0; i<size; i++)  {
      output[i] = matrix[smallest][i];
   }
    
   return size;
}



//////////////////////////////
//
// calculatePrimeForm -- the input is the normal form vector.
//   The output will be the prime form vector.  Returns the
//   size of the prime-form vector which should be the
//   same as the normal-form vector.
//

int calculatePrimeForm(long *output, long *input, long size) {
   long i;
   long dist;
   long normalTotal = 0;
   long invertedTotal = 0;
   long inverted[12];
   long newinverted[12];

   // transpose the set so the first element is 0
   dist = input[0];
   for (i=0; i<size; i++)  {
      output[i] = input[i] - dist;
      if (output[i] < 0)  {
         output[i] += 12;
      }
      normalTotal += output[i];
   }

   // invert the set 
   for (i=0; i<size; i++)  {
      inverted[i] = 0 - (input[i]);
      if (inverted[i] < 0)  {
         inverted[i] += 12;
      }
   }

   // find the normal form of the inverted set 
   qsort(inverted, size, sizeof(long), longcompare);
   calculateNormalForm(newinverted, inverted, size);
   for (i=0; i<size; i++) {
      inverted[i] = newinverted[i];
   }

   // transpose the inverted set to start on 0 
   dist = inverted[0];
   for (i=0; i<size; i++)  {
      inverted[i] = inverted[i] - dist;
      if (inverted[i] < 0)  {
         inverted[i] += 12;
      }
      invertedTotal += inverted[i];
   }

   // figure out which one (the inverted or the regular) set
   // is more packed to the left.  i think we can do this by
   // simply adding up the members.  I hope that works.
   // anyway, we've been keeping track of the total above.

   if (invertedTotal < normalTotal)  {
      memcpy(output, inverted, sizeof(inverted));
   }

   return size;
}



//////////////////////////////
//
// returnSmallest -- This function will return the index of the smallest
//  element in the given array, -1 if there are any ties.
//

long returnSmallest(long *list, long count)  {
   long i;
   long tiny;
   long loc;
   long tim;

   loc  = 0;
   tiny = list[0];
   tim  = 13;  

   for (i=1; i<count; i++)  {
     if (list[i] < tiny)  {
        loc = i;
     }
     if (list[i] == tiny)  {
        tim = list[i];
     }
   }       

   if (tiny == tim)  {
      return -1;
   }

   return loc;
}



//////////////////////////////
//
// calculateICV -- 
//    calculate the Interval Class Vector for the incoming set
//    we do this by measuring the intervals between each member
//    and all the members greater than it, and keep a tally
//    of how many of each interval class there is.  Returns
//    the size of the output vector (which is always 6).
//

long calculateICV(long *output, long *input, long size) {
   int i = 0;
   int j = 0;
   int k = 0;
    
   for (i=0; i<6; i++) {          //  clear the ICV array 
      output[i] = 0;
   }

   for (i=0; i<size-1; i++)  {
      for (j=i+1; j<size; j++)  {
         k = getIntervalClass(input[j] - input[i]);
         if (k >= 0) {
            output[k]++;
         }
      }
   }

   return 6;
} 



//////////////////////////////
//
// getIntervalClass -- put the interval value in to the range from 0 to 5.
//

long getIntervalClass(long input) {
   long output = (input + 12000) % 12;
   if (output > 6) {
      return 12 - output - 1;
   } else {
      return output - 1;
   }
}



//////////////////////////////
//
// longcompare -- compare two integers for ordering
//

int longcompare(const void* a, const void* b) {
   if (*((long*)a) < *((long*)b)) {
      return -1;
   } else if (*((long*)a) > *((long*)b)) {
      return 1;
   } else {
      return 0;
   }
}