//
// Programmer:    Craig Stuart Sapp <craig@ccrma.stanford.edu>
// Creation Date: Sun Feb  3 16:03:55 PST 2008
// Last Modified: Sun Feb  3 16:03:58 PST 2008
// Filename:      ...sig/examples/all/mkeyscape.cpp
// Web Address:   http://sig.sapp.org/examples/museinfo/humdrum/mkeyscape.cpp
// Syntax:        C++; museinfo
//
// Description:   Create Fast Keyscapes using MIDI file or Humdrum file.
//

#include "humdrum.h"
#include "MidiFile.h"

#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>             /* for qsort and bsearch functions */

#ifndef OLDCPP
   #include <iostream>
   #include <fstream>
#else
   #include <iostream.h>
   #include <fstream.h>
#endif

#define HISTTYPE double

#define UNKNOWNFILE 0
#define HUMDRUMFILE 1 
#define MIDIFILE 2 

// function declarations:
void     checkOptions(Options& opts, int argc, char** argv);
void     example(void);
void     usage(const char* command);
double   loadHistogramFromHumdrumFile
                                (Array<Array<HISTTYPE> >& histogram,
                                 const char* filename, int segments);
double   loadHistogramFromMidiFile
                                (Array<Array<HISTTYPE> >& histogram,
                                 const char* filename, int segments);
void     printBaseHistogram(Array<Array<HISTTYPE> >& histogram);
void     printBaseHistogramHumdrumStyle     
                                (Array<Array<HISTTYPE> >& histogram, 
				 double width);
void     printNormalizedHistogram(Array<Array<HISTTYPE> >& histogram);
void     addToHistogramInteger(Array<Array<int> >& histogram, int pc,
                                 double start, double dur, double tdur,
                                 int segments);
void     addToHistogramDouble(Array<Array<double> >& histogram, int pc,
                                 double start, double dur, double tdur,
                                 int segments);
void     fillAllHistograms(Array<Array<Array<HISTTYPE> > >& histograms);
void     printBest(Array<Array<Array<HISTTYPE> > >& histogram);
void     calculateBestKeys(Array<Array<Array<HISTTYPE> > >& histograms);
void     identifyKey(Array<HISTTYPE>& histogram);
void     displayRawAnalysis(Array<Array<Array<HISTTYPE> > >& histogram);
void     displayAnalysisHistogram(Array<Array<Array<HISTTYPE> > >& histograms);
void     identifyKeyDouble(Array<HISTTYPE>& histogram);
void     printPPM(Array<Array<Array<HISTTYPE> > >& histograms);
double   pearsonCorrelation(int size, double* x, double* y);
void     setFilterOptions(Array<int>& channelfilter, 
                                 const char* exclude);
void     processColorFile(const char* filename, HumdrumFile& cfile);
int      ranksort(const void* A, const void* B);
void     fillWeightsWithAardenEssen
                                (Array<HISTTYPE>& maj, Array<HISTTYPE>& min);
void     fillWeightsWithBellmanBudge
                                (Array<HISTTYPE>& maj, Array<HISTTYPE>& min);
void     fillWeightsWithKostkaPayne
                                (Array<HISTTYPE>& maj, Array<HISTTYPE>& min);
void     fillWeightsWithKrumhanslKessler
                                (Array<HISTTYPE>& maj, Array<HISTTYPE>& min);
void     fillWeightsWithSimpleValues
                                (Array<HISTTYPE>& maj, Array<HISTTYPE>& min);
void     processWeights(const char* filename, Array<HISTTYPE>& major, 
                                 Array<HISTTYPE>& minor);
void     printColorMap(Array<const char*>& colorindex);
void     printWeights(Array<HISTTYPE>& maj, Array<HISTTYPE>& min);

// User interface variables:
Options   options;
int       segments     = 300;   // used with -s option
int       rawQ         = 0;     // used with -r option
int       histQ        = 0;     // used with --hist option
int       khistQ       = 0;     // used with --khist option
int       transpose    = 0;     // used with -t option

Array<int> channelfilter;       // used with -x option
Array<const char*> colorindex;  // used with -c option
HumdrumFile colorfile;          // used with -c option (needs to be global)
Array<HISTTYPE> majorweights;
Array<HISTTYPE> minorweights;

//////////////////////////////////////////////////////////////////////////

int main(int argc, char** argv) {
   channelfilter.setSize(16);
   channelfilter.allowGrowth(0);
   channelfilter.setAll(1);

   majorweights.setSize(12);
   majorweights.allowGrowth(0);
   majorweights.zero(0);

   minorweights.setSize(12);
   minorweights.allowGrowth(0);
   minorweights.zero(0);

   fillWeightsWithBellmanBudge(majorweights, minorweights);

   colorindex.setSize(26);
   colorindex.allowGrowth(0);
   colorindex[0]  = "0 255 0";		// C major
   colorindex[1]  = "18 237 73";	// C-sharp major
   colorindex[2]  = "63 95 255";	// D major
   colorindex[3]  = "218 9 73";		// E-flat major
   colorindex[4]  = "255 0 0";		// E major
   colorindex[5]  = "255 255 0";	// F major
   colorindex[6]  = "182 255 0";	// F-sharp major
   colorindex[7]  = "63 191 255";	// G major
   colorindex[8]  = "109 50 255";	// A-flat major
   colorindex[9]  = "127 31 255";	// A major
   colorindex[10] = "255 132 0";	// B-flat major
   colorindex[11] = "255 127 0";	// B major
   colorindex[12] = "0 191 0";		// C minor
   colorindex[13] = "14 178 55";	// C-sharp minor
   colorindex[14] = "47 71 191";	// D minor
   colorindex[15] = "164 7 55";		// E-flat minor
   colorindex[16] = "191 0 0";		// E minor
   colorindex[17] = "191 191 0";	// F minor
   colorindex[18] = "137 191 0";	// F-sharp minor
   colorindex[19] = "47 143 191";	// G minor
   colorindex[20] = "47 92 191";	// G-sharp minor
   colorindex[21] = "95 23 191";	// A minor
   colorindex[22] = "191 99 0";		// B-flat minor
   colorindex[23] = "191 95 0";		// B minor
   colorindex[24] = "0 0 0";		// silence
   colorindex[25] = "255 255 255";	// background

   // process the command-line options
   checkOptions(options, argc, argv);

   Array<Array<Array<HISTTYPE> > > histograms;
   histograms.setSize(segments);
   histograms.allowGrowth(0);
   int i;
   int j;
   for (i=0; i<segments; i++) {
      histograms[i].setSize(i+1);
      histograms[i].allowGrowth(0);
      for (j=0; j<i+1; j++) {
         histograms[i][j].setSize(13);   // last spot is for key result
         histograms[i][j].allowGrowth(0);
         histograms[i][j].zero();
      }
   }

   double totalduration = 0;
   const char* filename = "";
   if (options.getArgCount() > 0) {
      filename = options.getArg(1);
   }
   int fnlength = strlen(filename);
   int filetype = UNKNOWNFILE;
   if (fnlength != 0) {
      if (strcmp(filename + (fnlength - 4), ".mid") == 0) {
         filetype = MIDIFILE;
         totalduration = loadHistogramFromMidiFile(
               histograms[histograms.getSize()-1], filename, segments);
      } else {
         filetype = HUMDRUMFILE;
         totalduration = loadHistogramFromHumdrumFile(
               histograms[histograms.getSize()-1], filename, segments);
      }
   } else {
      filetype = HUMDRUMFILE;
      totalduration = loadHistogramFromHumdrumFile(
               histograms[histograms.getSize()-1], filename, segments);
   }

   fillAllHistograms(histograms);

   if (histQ) {
      // printBaseHistogram(histograms[histograms.getSize()-1]);
      printBaseHistogramHumdrumStyle(histograms[histograms.getSize()-1],
            totalduration);
      exit(0);
   }
   //printNormalizedHistogram(histograms[0]);
   //printBaseHistogram(histograms[1]);
   //cout << "========================" << endl;
   //printBaseHistogram(histograms[histograms.getSize()-1]);

   //identifyKey(histograms[0][0]);
   //exit(0);

   calculateBestKeys(histograms);

   //printBest(histograms);

   if (rawQ) {
      displayRawAnalysis(histograms);
      exit(0);
   } else if (khistQ) {
      displayAnalysisHistogram(histograms);
      exit(0);
   }
	    
	    

   printPPM(histograms);

   return 0;
}

//////////////////////////////////////////////////////////////////////////

//////////////////////////////
//
// printWeights --
//

void printWeights(Array& maj, Array& min) {
   cout << "**kern	**weight\n";
   cout << "!!major weights (using C as the tonic)\n";
   cout << "C	" << maj[0]  << "\n";
   cout << "C#	" << maj[1]  << "\n";
   cout << "D	" << maj[2]  << "\n";
   cout << "E-	" << maj[3]  << "\n";
   cout << "E	" << maj[4]  << "\n";
   cout << "F	" << maj[5]  << "\n";
   cout << "F#	" << maj[6]  << "\n";
   cout << "G	" << maj[7]  << "\n";
   cout << "G#	" << maj[8]  << "\n";
   cout << "A	" << maj[9]  << "\n";
   cout << "B-	" << maj[10] << "\n";
   cout << "B	" << maj[11] << "\n";
   cout << "!!minor weights (using c as the tonic)\n";
   cout << "c	" << min[0]  << "\n";
   cout << "c#	" << min[1]  << "\n";
   cout << "d	" << min[2]  << "\n";
   cout << "e-	" << min[3]  << "\n";
   cout << "e	" << min[4]  << "\n";
   cout << "f	" << min[5]  << "\n";
   cout << "f#	" << min[6]  << "\n";
   cout << "g	" << min[7]  << "\n";
   cout << "g#	" << min[8]  << "\n";
   cout << "a	" << min[9]  << "\n";
   cout << "b-	" << min[10] << "\n";
   cout << "b	" << min[11] << "\n";
   cout << "*-	*-\n";
}



//////////////////////////////
//
// printColorMap --
//

void printColorMap(Array& ci) {
   cout << "**index	**pixel	**comment\n";
   cout << "!! Minor Keys\n";
   cout << "0	" << ci[0]	<< "\tC major\n";
   cout << "1	" << ci[1]	<< "\tC-sharp major\n";
   cout << "2	" << ci[2]	<< "\tD major\n";
   cout << "3	" << ci[3]	<< "\tE-flat major\n";
   cout << "4	" << ci[4]	<< "\tE major\n";
   cout << "5	" << ci[5]	<< "\tF major\n";
   cout << "6	" << ci[6]	<< "\tF-sharp major\n";
   cout << "7	" << ci[7]	<< "\tG major\n";
   cout << "8	" << ci[8]	<< "\tA-flat major\n";
   cout << "9	" << ci[9]	<< "\tA major\n";
   cout << "10	" << ci[10]	<< "\tB-flat major\n";
   cout << "11	" << ci[11]	<< "\tB major\n";
   cout << "!! Minor Keys\n";
   cout << "12	" << ci[12]	<< "\tC minor\n";
   cout << "13	" << ci[13]	<< "\tC-sharp minor\n";
   cout << "14	" << ci[14]	<< "\tD minor\n";
   cout << "15	" << ci[15]	<< "\tE-flat minor\n";
   cout << "16	" << ci[16]	<< "\tE minor\n";
   cout << "17	" << ci[17]	<< "\tF minor\n";
   cout << "18	" << ci[18]	<< "\tF-sharp minor\n";
   cout << "19	" << ci[19]	<< "\tG minor\n";
   cout << "20	" << ci[20]	<< "\tG-sharp minor\n";
   cout << "21	" << ci[21]	<< "\tA minor\n";
   cout << "22	" << ci[22]	<< "\tB-flat minor\n";
   cout << "23	" << ci[23]	<< "\tB minor\n";
   cout << "!! Other Colors\n";
   cout << "24	" << ci[24]	<< "\tsilence\n";
   cout << "25	" << ci[25]	<< "\tbackground\n";
   cout << "*-	*-	*-\n";
}


//////////////////////////////
//
// printPPM --
//

void printPPM(Array > >& histograms) {

   // Pitch to Color translations

   int width = histograms.getSize();

   #define BGCOLOR colorindex[25]

   cout << "P3\n";
   cout << width * 2 << " " << width  << "\n";
   cout << "255\n";

   int blankcells;
   int i;
   int j;
   const char* color;
   for (i=0; i<width; i++) {
      blankcells = width - histograms[i].getSize();
      for (j=0; j<blankcells; j++) {
         cout << ' ' << BGCOLOR;
      }
      for (j=0; j<histograms[i].getSize(); j++) {
         color = colorindex[(int)histograms[i][j][12]];
         cout << ' ' << color << ' ' << color;
      }
      for (j=0; j<blankcells; j++) {
         cout << ' ' << BGCOLOR;
      }
      cout << "\n";
   }
}



//////////////////////////////
//
// calculateBestKeys --
//

void calculateBestKeys(Array > >& histograms) {
   int i;
   int j;

   for (i=0; i<histograms.getSize(); i++) {
      for (j=0; j<histograms[i].getSize(); j++) {
          identifyKeyDouble(histograms[i][j]);
      }
   }
}



////////////////////////////////////////
//
// identifyKeyDouble --
//

void identifyKeyDouble(Array& histogram) {
   int i;

   double h[24];
   for (i=0; i<12; i++) {
      h[i] = histogram[i];
      h[i+12] = h[i];
   }

   double keysum[24];

   double testsum = 0.0;
   for (i=0; i<12; i++) {
      testsum += h[i];
      keysum[i]    = pearsonCorrelation(12, majorweights.getBase(), h+i);
      keysum[i+12] = pearsonCorrelation(12, minorweights.getBase(), h+i);
   }

   if (testsum == 0.0) {
      histogram[12] = 24;  // empty histogram, so going to display black
      return;
   }

   // find max value
   int besti = 0;
   HISTTYPE bestsum = keysum[0];
   for (i=1; i<24; i++) {
      if (keysum[i] > bestsum) {
         besti = i;
         bestsum = keysum[i];
      }
   }

   histogram[12] = besti;
}



////////////////////////////////////////
//
// identifyKey --
//

void identifyKey(Array& histogram) {
   int i;

   int h[24];
   int scaled[24];
   int mean = 0;
   for (i=0; i<12; i++) {
      h[i] = (int)histogram[i];
      h[i+12] = h[i];
      scaled[i] = int(h[i] * 0.75 + 0.5);
      scaled[i+12] = scaled[i];
      mean += h[i];
   }
   mean = (int)(mean / 12.0 + 0.5);
   int keysum[24];

   //double majorw[12] = {2,0,1,0,1,1,0,2,0,1,0,1};
   //double minorw[12] = {2,0,1,1,0,1,0,2,1,0,1,0};

   for (i=0; i<12; i++) {
      // keysums[i]    = pearsonCorrelation(12, majorw, h+i);
      // keysums[i+12] = pearsonCorrelation(12, minorw, h+i);

      keysum[i] += ((h[i]-mean) << 1) - scaled[i];
      keysum[i] += -scaled[i+1];
      keysum[i] += h[i+2] - mean - scaled[i+2];
      keysum[i] += -scaled[i+3];
      keysum[i] += h[i+4] - mean - scaled[i+4];
      keysum[i] += h[i+5] - mean - scaled[i+5];
      keysum[i] += -scaled[i+6];
      keysum[i] += ((h[i]-mean) << 1) - scaled[i+7];
      keysum[i] += -scaled[i+8];
      keysum[i] += h[i+9] - mean - scaled[i+9];
      keysum[i] += -scaled[i+10];
      keysum[i] += h[i+11] - mean - scaled[i+11];

      keysum[i+12] += ((h[i]-mean) << 1) - scaled[i];
      keysum[i+12] += -scaled[i+1];
      keysum[i+12] += h[i+2] - mean - scaled[i+2];
      keysum[i+12] += h[i+3] - mean - scaled[i+3];
      keysum[i+12] += -scaled[i+4];
      keysum[i+12] += h[i+5] - mean - scaled[i+5];
      keysum[i+12] += -scaled[i+6];
      keysum[i+12] += ((h[i]-mean) << 1) - scaled[i+7];
      keysum[i+12] += h[i+8] - mean - scaled[i+8];
      keysum[i+12] += -scaled[i+9];
      keysum[i+12] += h[i+10] - mean - scaled[i+10];
      keysum[i+12] += -scaled[i+11];
   }

   // find max value
   int besti = 0;
   HISTTYPE bestsum = keysum[0];
   for (i=1; i<24; i++) {
// cout << i << ":\t" << keysum[i] << endl;
      if (keysum[i] > bestsum) {
         besti = i;
         bestsum = keysum[i];
      }
   }

   histogram[12] = besti;
}



//////////////////////////////
//
// fillAllHistograms --
//

void fillAllHistograms(Array > >& histograms) {
   int size = histograms.getSize();
   int i, j, k;
   for (i=size-2; i>=0; i--) {
      for (j=0; j<histograms[i].getSize(); j++) {
         for (k=0; k<12; k++) {
            histograms[i][j][k] = histograms[i+1][j][k] +
                                  histograms[size-1][size-1-i+j][k];
         }
      }
   }
}



//////////////////////////////
//
// displayAnalysisHistogram --
//

void displayAnalysisHistogram(Array > >& histograms) {
   int i, j;
   int key = 0;
   int size = histograms.getSize();
   int total = size * (size + 1) / 2;   // triangle number
   Array<int> counts;
   counts.setSize(25);
   counts.allowGrowth(0);
   counts.zero();
   for (i=0; i<histograms.getSize(); i++) {
      for (j=0; j<histograms[i].getSize(); j++) {
         key = (int)histograms[i][j][12];
         counts[key]++;
      }
   }

   Array<int*> rank;
   rank.setSize(25);
   rank.allowGrowth(0);
   Array<int> cc;
   cc.setSize(25);
   for (i=0; i<rank.getSize(); i++) {
      cc[i] = counts[i];
      rank[i] = &(cc[i]);
     
   }
   qsort(rank.getBase(), rank.getSize(), sizeof(int*), ranksort);
   for (i=0; i<cc.getSize(); i++) {
      *(rank[i]) = i+1;
   }

   cout << "!! Count Total = " << total << "\n";
   cout << "**key\t**rank\t**count\t**fraction\n";
   for (i=0; i<counts.getSize(); i++) {
      key = i;
      if (key < 12) {
         key = key+1;
      } else if (key < 24) {
         key = -(key - 11);
      } else {
         key = 0;
      }
      cout << key   << "\t";
      cout << cc[i] << "\t";
      cout << counts[i] << "\t";
      cout << (double)counts[i] / total << "\n";
   }

   int majorsum = 0;
   int minorsum = 0;
   for (i=0; i<12; i++) {
      majorsum += counts[i];
   }
   for (i=12; i<24; i++) {
      minorsum += counts[i];
   }
   cout << "*-\t*-\t*-\t*-\n";
   cout << "!! Major fraction: " << (double)majorsum / total << "\n";
   cout << "!! Minor fraction: " << (double)minorsum / total << "\n";
   if (majorsum > minorsum) {
      cout << "!! The music is happy\n";
   } else {
      cout << "!! The music is sad\n";
   }

}



//////////////////////////////
//
// ranksort -- sort counts by largest first
//

int ranksort(const void* A, const void* B) {
   int& a = *(*((int**)A));
   int& b = *(*((int**)B));
   if (a < b) {
      return +1;
   } else if (a > b) {
      return -1;
   } else {
      return 0;
   }
}



//////////////////////////////
//
// displayRawAnalysis --
//

void displayRawAnalysis(Array > >& histogram) {
   int i, j;
   int key;
   for (i=0; i<histogram.getSize(); i++) {
      for (j=0; j<histogram[i].getSize(); j++) {
         if (j > 0) {
            cout << '\t';
         } 
         key = (int)histogram[i][j][12];
         if (key < 12) {
            cout << key+1;             // major key C=1 C#=2 D=3, etc.
         } else if (key < 24) {
            cout << -(key - 11);       // minor key C=-1 C#=-2 D=-3, etc.
         } else {
            cout << 0;                 // silence
         }
      }
      cout << '\n';
   }
}



//////////////////////////////
//
// printBaseHistogramHumdrumStyle --
//

void printBaseHistogramHumdrumStyle(Array<Array<HISTTYPE> >& histogram,
      double totalduration) {
   int i;
   int j;
   cout << "!! Total Duration of Music:\t" << totalduration << endl;
   cout << "!! Duration Units per Frame:\t" 
	<< totalduration / histogram.getSize() << endl;
   cout << "**frame\t**bin00\t**bin01\t**bin02\t**bin03"
        << "\t**bin04\t**bin05\t**bin06\t**bin07\t**bin08\t**bin09"
        << "\t**bin10\t**bin11\n";
   cout << "!\t!C\t!C#\t!D\t!D#\t!E\t!F\t!F#\t!G\t!G#\t!A\t!A#\t!B\n";
   for (i=0; i<histogram.getSize(); i++) {
      cout << i << ':';
      for (j=0; j<histogram[i].getSize()-1; j++) {
         // subtracting one from limit becuase last item
         // is storage for key analysis
         cout << '\t' << histogram[i][j];
      }
      cout << '\n';
   }
   cout << "*-\t*-\t*-\t*-\t*-\t*-\t*-\t*-\t*-\t*-\t*-\t*-\t*-\n";
}

//////////////////////////////
//
// printBaseHistogram --
//

void printBaseHistogram(Array >& histogram) {
   int i;
   int j;
   for (i=0; i<histogram.getSize(); i++) {
      cout << i << ':';
      for (j=0; j<histogram[i].getSize(); j++) {
         cout << '\t' << histogram[i][j];
      }
      cout << '\n';
   }
}



//////////////////////////////
//
// printNormalizedHistogram --
//

void printNormalizedHistogram(Array >& histogram) {
   int i;
   int j;
   double sum = 0;
   for (i=0; i<histogram.getSize(); i++) {
      cout << i << ':';
      sum = 0;
      for (j=0; j<histogram[i].getSize(); j++) {
         sum += histogram[i][j];
      }
      for (j=0; j<histogram[i].getSize(); j++) {
         cout << '\t' << (double)histogram[i][j] / sum;
      }
      cout << '\n';
   }
}



//////////////////////////////
//
// printBest --
//

void printBest(Array > >& histogram) {
   int i;
   int j;
   for (i=0; i<histogram.getSize(); i++) {
      cout << i << ':';
      for (j=0; j<histogram[i].getSize(); j++) {
         cout << '\t' << histogram[i][j][12];
      }
      cout << '\n';
   }
}



//////////////////////////////
//
// loadHistogramFromHumdrumFile --
//

double loadHistogramFromHumdrumFile(Array<Array<HISTTYPE> >& histograms,
   const char* filename, int segments) {

   HumdrumFile infile;
   if (strcmp(filename, "") == 0) {
      infile.read(cin);
   } else {
      infile.read(filename);
   }
   infile.analyzeRhythm("4");
   double totalduration = infile.getTotalDuration();

   double duration;
   int i;
   int j;
   int k;
   char buffer[10000] = {0};
   int pitch;
   double start;
   int tokencount;
   for (i=0; i<infile.getNumLines(); i++) {
      if (infile[i].getType() !=  E_humrec_data) {
         continue;
      }
      start = infile[i].getAbsBeat();
      for (j=0; j<infile[i].getFieldCount(); j++) {
         if (strcmp(infile[i].getExInterp(j), "**kern") != 0) {
            continue;
         }
         if (strcmp(infile[i][j], ".") == 0) {
            continue; // ignore null tokens
         }
         tokencount = infile[i].getTokenCount(j);
         for (k=0; k<tokencount; k++) {
            infile[i].getToken(buffer, j, k);
            if (strcmp(buffer, ".") == 0) {
               continue;  // ignore illegal inline null tokens
            }
            pitch = Convert::kernToMidiNoteNumber(buffer);
            if (pitch < 0) {
               continue;  // ignore rests or strange objects
            }
            pitch = pitch % 12;  // convert to chromatic pitch-class
            duration = Convert::kernToDuration(buffer);
            if (duration <= 0.0) {
               continue;   // ignore grace notes and strange objects
            }
            
            addToHistogramDouble(histograms, pitch,
                  start, duration, totalduration, segments);
         }
      }
   }

   return totalduration;
}



//////////////////////////////
//
// addToHistogramDouble -- fill the histogram in the right spots.
//

void addToHistogramDouble(Array<Array<double> >& histogram, int pc, 
      double start, double dur, double tdur, int segments) {

   pc = (pc + transpose + 12) % 12;

   double startseg = start / tdur * segments;
   double startfrac = startseg - (int)startseg;

   double segdur = dur / tdur * segments;

   if (segdur <= 1.0 - startfrac) {
      histogram[(int)startseg][pc] += segdur;
      return;
   } else if (1.0 - startfrac > 0.0) {
      histogram[(int)startseg][pc] += (1.0 - startfrac);
      segdur -= (1.0 - startfrac);
   }

   int i = (int)(startseg + 1);
   while (segdur > 0.0 && i < histogram.getSize()) {
      if (segdur < 1.0) {
         histogram[i][pc] += segdur;
         segdur = 0.0;
      } else {
         histogram[i][pc] += 1.0;
         segdur -= 1.0;
      }
      i++;
   }
}



//////////////////////////////
//
// addToHistogramInteger -- fill the histogram in the right spots.
//

#define WID 100000

void addToHistogramInteger(Array<Array<int> >& histogram, int pc, 
      double start, double dur, double tdur, int segments) {

   double startseg = start / tdur * segments;
   double startfrac = startseg - (int)startseg;

   double segdur = dur / tdur * segments;

   if (segdur <= 1.0 - startfrac) {
      histogram[(int)startseg][pc] += int(segdur * WID + 0.5);
      return;
   } else if (1.0 - startfrac > 0.0) {
      histogram[(int)startseg][pc] += int((1.0 - startfrac) * WID + 0.5);
      segdur -= (1.0 - startfrac);
   }

   int i = (int)(startseg + 1);
   while (segdur > 0.0 && i < histogram.getSize()) {
      if (segdur < 1.0) {
         histogram[i][pc] += int(segdur * WID + 0.5);
         segdur = 0.0;
      } else {
         histogram[i][pc] += WID;
         segdur -= 1.0;
      }
      i++;
   }
}



//////////////////////////////
//
// loadHistogramFromMidiFile --
//

double loadHistogramFromMidiFile(Array<Array<HISTTYPE> >& histogram,
   const char* filename, int segments) {

   MidiFile midifile(filename);
   midifile.absoluteTime();
   midifile.joinTracks();

   Array<int> ontimes(128*16);
   Array<int> onvelocities(128*16);
   ontimes.setAll(-1);
   onvelocities.zero();

   int i;
   int command = 0;
   int totalduration = midifile.getEvent(0,midifile.getNumEvents(0)-1).time;
   // using the last item as the total duration is not so great because
   // some MIDI files have some junk messages long after the music 
   // has stopped (2_ase.mid is an example), so search backwards
   // through the midifile for the last note-on or note-off message.
   for (i=midifile.getNumEvents(0)-1; i>=0; i--) {
      command = midifile.getEvent(0,i).data[0] & 0xf0;
      if (command == 0x90 || command == 0x80) {
         totalduration = midifile.getEvent(0, i).time;
         break;
      }
   }
   
   int key;
   int channel;
   int duration;
   int ontime;
   for (i=0; i<midifile.getNumEvents(0); i++) {
      command = midifile.getEvent(0,i).data[0] & 0xf0;
      channel = midifile.getEvent(0,i).data[0] & 0x0f;
      if (channelfilter[channel] == 0) {
         // ignore events on this channel
         continue;
      }
      if (command == 0x90 && midifile.getEvent(0,i).data[2] != 0) {
         // store note-on velocity and time.
         key = midifile.getEvent(0,i).data[1];
         ontime = midifile.getEvent(0,i).time;
         if (ontimes[key * channel] > -1) {
            // the previous note was not turned off, to turn
            // it off now and store that note in the histogram
            duration = ontime - ontimes[key * channel];
            addToHistogramDouble(histogram, key % 12, ontimes[key*channel], 
                  duration, totalduration, segments);
            ontimes[key * channel] = ontime;
         } else {
            // no note exists in the slot, to store for later
            ontimes[key * channel] = ontime;
         }
      } else if (command == 0x90 || command == 0x80) {
         // process a note-off command
         key = midifile.getEvent(0,i).data[1];
         ontime = midifile.getEvent(0,i).time;
         if (ontimes[key * channel] > -1) {
            // process the note which has been waiting
            duration = ontime - ontimes[key * channel];
            addToHistogramDouble(histogram, key % 12, ontimes[key*channel], 
                  duration, totalduration, segments);
         }
         ontimes[key * channel] = -1;
      }
   }

   return totalduration;
}



//////////////////////////////
//
// pearsonCorrelation --
//

double pearsonCorrelation(int size, double* x, double* y) {

   double sumx  = 0.0;
   double sumy  = 0.0;
   double sumco = 0.0;
   double meanx = x[0];
   double meany = y[0];
   double sweep;
   double deltax;
   double deltay;

   int i;
   for (i=2; i<=size; i++) {
      sweep = (i-1.0) / i;
      deltax = x[i-1] - meanx;
      deltay = y[i-1] - meany;
      sumx  += deltax * deltax * sweep;
      sumy  += deltay * deltay * sweep;
      sumco += deltax * deltay * sweep;
      meanx += deltax / i;
      meany += deltay / i;
   }

   double popsdx = sqrt(sumx / size);
   double popsdy = sqrt(sumy / size);
   double covxy  = sumco / size;

   return covxy / (popsdx * popsdy);
}




//////////////////////////////
//
// checkOptions --
//

void checkOptions(Options& opts, int argc, char* argv[]) {
   opts.define("t|transpose=i:0", "Transposition by half-steps");
   opts.define("s|segments=i:300", "The height in pixel of the output plot");
   opts.define("r|raw=b", "Display the analyzed keys");
   opts.define("c|colorfile=s", "key to color mapping specification");
   opts.define("C|printcolors=b", "print key to color mapping");
   opts.define("D|no-drum=b", "Don't analyze the General MIDI drum track");
   opts.define("x|exclude=s", "exclude MIDI channel notes");
   opts.define("hist|histograms|histogram=b", "display the baseline histograms");
   opts.define("w|weights=s", "arbitrary set of pitch weights");
   opts.define("W|printweights=b", "display weights which will be used");

   opts.define("aa|aarden=b",        "load Aarden-Essen weights");
   opts.define("bb|bellman|budge=b", "load Bellman-Budge weights");
   opts.define("kk|krumhansl=b",     "load Krumhansl-Kessler weights");
   opts.define("kp|kostkapayne=b",   "load Kostka-Payne weights");
   opts.define("ss|simple|sapp=b",   "load Simple weights");

   opts.define("khist=b", "display the analysis key histogram");
   opts.define("author=b",  "author of program");
   opts.define("version=b", "compilation info");
   opts.define("example=b", "example usages");
   opts.define("help=b",  "short description");
   opts.process(argc, argv);

   // handle basic options:
   if (opts.getBoolean("author")) {
      cout << "Written by Craig Stuart Sapp, "
           << "craig@ccrma.stanford.edu, Jan 2008" << endl;
      exit(0);
   } else if (opts.getBoolean("version")) {
      cout << argv[0] << ", version: 30 Jan 2008" << endl;
      cout << "compiled: " << __DATE__ << endl;
      cout << MUSEINFO_VERSION << endl;
      exit(0);
   } else if (opts.getBoolean("help")) {
      usage(opts.getCommand());
      exit(0);
   } else if (opts.getBoolean("example")) {
      example();
      exit(0);
   }

   segments  = opts.getInteger("segments");
   histQ     = opts.getBoolean("histogram");
   khistQ    = opts.getBoolean("khist");
   rawQ      = opts.getBoolean("raw");
   transpose = opts.getInteger("transpose");
   if (transpose < -12 || transpose > 12) {
      cerr << "Error: transposition is out of range: " << transpose << endl;
      exit(1);
   } 

   setFilterOptions(channelfilter, opts.getString("exclude"));
   if (opts.getBoolean("no-drum")) {
      // turn off MIDI channel 10
      channelfilter[9] = 0;
   }

   if (opts.getBoolean("colorfile")) {
      processColorFile(opts.getString("colorfile"), colorfile);
   }
   if (opts.getBoolean("printcolors")) {
      printColorMap(colorindex);
      exit(0);
   }
	    
   
   if (opts.getBoolean("aarden")) {
      fillWeightsWithAardenEssen(majorweights, minorweights);
   } else if (opts.getBoolean("bellman")) {
      fillWeightsWithBellmanBudge(majorweights, minorweights);
   } else if (opts.getBoolean("krumhansl")) {
      fillWeightsWithKrumhanslKessler(majorweights, minorweights);
   } else if (opts.getBoolean("kostkapayne")) {
      fillWeightsWithKostkaPayne(majorweights, minorweights);
   } else if (opts.getBoolean("simple")) {
      fillWeightsWithSimpleValues(majorweights, minorweights);
   }

   if (opts.getBoolean("weights")) {
      processWeights(opts.getString("weights"), majorweights, minorweights);
   }

   if (opts.getBoolean("printweights")) {
      printWeights(majorweights, minorweights);  
      exit(0);
   }

}



//////////////////////////////
//
// fillWeightsWith* -- Set key weights to specified default values.
//

void fillWeightsWithAardenEssen(Array& maj, Array& min) {
   // as found in Bret Aarden's dissertation
   maj[0]  = 17.7661  ;  // C major weights
   maj[1]  =  0.145624;  // C#
   maj[2]  = 14.9265  ;  // D
   maj[3]  =  0.160186;  // D#
   maj[4]  = 19.8049  ;  // E
   maj[5]  = 11.3587  ;  // F
   maj[6]  =  0.291248;  // F#
   maj[7]  = 22.062   ;  // G
   maj[8]  =  0.145624;  // G#
   maj[9]  =  8.15494 ;  // A
   maj[10] =  0.232998;  // A#
   maj[11] =  4.95122 ;  // B
   min[0]  = 18.2648  ;  // c minor weights
   min[1]  =  0.737619;  // c#
   min[2]  = 14.0499  ;  // d
   min[3]  = 16.8599  ;  // d#
   min[4]  =  0.702494;  // e
   min[5]  = 14.4362  ;  // f
   min[6]  =  0.702494;  // f#
   min[7]  = 18.6161  ;  // g
   min[8]  =  4.56621 ;  // g#
   min[9]  =  1.93186 ;  // a
   min[10] =  7.37619 ;  // a#
   min[11] =  1.75623 ;  // b
} 


void fillWeightsWithBellmanBudge(Array& maj, Array& min) {
   // as found in Bellman CMMR 2005 paper, derived from Budge 1943 dissertation
   maj[0]  = 16.80;	// C major weights
   maj[1]  =  0.86;	// C#
   maj[2]  = 12.95;	// D
   maj[3]  =  1.41;	// D#
   maj[4]  = 13.49;	// E
   maj[5]  = 11.93;	// F
   maj[6]  =  1.25;	// F#
   maj[7]  = 20.28;	// G
   maj[8]  =  1.80;	// G#
   maj[9]  =  8.04;	// A
   maj[10] =  0.62;	// A#
   maj[11] = 10.57;	// B
   min[0]  = 18.16;	// c minor weights
   min[1]  =  0.69;	// c#
   min[2]  = 12.99;	// d
   min[3]  = 13.34;	// d#
   min[4]  =  1.07;	// e
   min[5]  = 11.15;	// f
   min[6]  =  1.38;	// f#
   min[7]  = 21.07;	// g
   min[8]  =  7.49;	// g#
   min[9]  =  1.53;	// a
   min[10] =  0.92;	// a#
   min[11] = 10.21;	// b
}


void fillWeightsWithKrumhanslKessler(Array& maj, Array& min) {
   // as found in Carol Krumhansl's 1990 book
   maj[0]  = 6.35; 	// C major weights
   maj[1]  = 2.23; 	// C#
   maj[2]  = 3.48; 	// D
   maj[3]  = 2.33; 	// D#
   maj[4]  = 4.38; 	// E
   maj[5]  = 4.09; 	// F
   maj[6]  = 2.52; 	// F#
   maj[7]  = 5.19; 	// G
   maj[8]  = 2.39; 	// G#
   maj[9]  = 3.66; 	// A
   maj[10] = 2.29; 	// A#
   maj[11] = 2.88;	// B
   min[0]  = 6.33; 	// c minor weights
   min[1]  = 2.68; 	// c#
   min[2]  = 3.52; 	// d
   min[3]  = 5.38; 	// d#
   min[4]  = 2.60; 	// e
   min[5]  = 3.53; 	// f
   min[6]  = 2.54; 	// f#
   min[7]  = 4.75; 	// g
   min[8]  = 3.98; 	// g#
   min[9]  = 2.69; 	// a
   min[10] = 3.34; 	// a#
   min[11] = 3.17;	// b
}


void fillWeightsWithKostkaPayne(Array& maj, Array& min) {
   // found in David Temperley: Music and Probability 2006
   maj[0]  = 0.748;	// C major weights
   maj[1]  = 0.060;	// C#
   maj[2]  = 0.488;	// D
   maj[3]  = 0.082;	// D#
   maj[4]  = 0.670;	// E
   maj[5]  = 0.460;	// F
   maj[6]  = 0.096;	// F#
   maj[7]  = 0.715;	// G
   maj[8]  = 0.104;	// G#
   maj[9]  = 0.366;	// A
   maj[10] = 0.057;	// A#
   maj[11] = 0.400;	// B
   min[0]  = 0.712;	// c minor weights
   min[1]  = 0.084;	// c#
   min[2]  = 0.474;	// d
   min[3]  = 0.618;	// d#
   min[4]  = 0.049;	// e
   min[5]  = 0.460;	// f
   min[6]  = 0.105;	// f#
   min[7]  = 0.747;	// g
   min[8]  = 0.404;	// g#
   min[9]  = 0.067;	// a
   min[10] = 0.133;	// a#
   min[11] = 0.330;	// b
}


void fillWeightsWithSimpleValues(Array& maj, Array& min) {
   maj[0]  = 2;		// C major weights
   maj[1]  = 0;		// C#
   maj[2]  = 1;		// D
   maj[3]  = 0;		// D#
   maj[4]  = 1;		// E
   maj[5]  = 1;		// F
   maj[6]  = 0;		// F#
   maj[7]  = 2;		// G
   maj[8]  = 0;		// G#
   maj[9]  = 1;		// A
   maj[10] = 0;		// A#
   maj[11] = 1;		// B
   min[0]  = 2;		// c minor weights
   min[1]  = 0;		// c#
   min[2]  = 1;		// d
   min[3]  = 1;		// d#
   min[4]  = 0;		// e
   min[5]  = 1;		// f
   min[6]  = 0;		// f#
   min[7]  = 2;		// g
   min[8]  = 1;		// g#
   min[9]  = 0;		// a
   min[10] = 1;		// a#
   min[11] = 0;		// b
}



//////////////////////////////
//
// processWeights --
//

void processWeights(const char* filename, Array<HISTTYPE>& major, 
      Array<HISTTYPE>& minor) {

   HumdrumFile wfile;
   wfile.read(filename);
   int i, j;
   int key;
   double value;
   for (i=0; i<wfile.getNumLines(); i++) {
      if (wfile[i].getType() != E_humrec_data) {
         continue;
      }
      key = -1;
      value = -1.0;
      for (j=0; j<wfile[i].getFieldCount(); j++) {
         if (strcmp(wfile[i].getExInterp(j), "**kern") == 0) {
            key = Convert::kernToMidiNoteNumber(wfile[i][j]) % 12;
            if (islower(wfile[i][j][0])) {
               key += 12;
            }
         } else if (strcmp(wfile[i].getExInterp(j), "**weight") == 0) {
            sscanf(wfile[i][j], "%lf", &value);
         }
      }
      if ((key >= 0) && (key < 24) && (value != -1.0)) {
         if (key < 12) {
            major[key] = value;
         } else {
            minor[key-12] = value;
         }
      }
   }
}



//////////////////////////////
//
// processColorFile --
//

void processColorFile(const char* filename, HumdrumFile& cfile) {
   cfile.read(filename);
   int i;
   int j;
   int index;
   int test1;
   int test2;
   int test3;
   const char* string;

   for (i=0; i<cfile.getNumLines(); i++) {
      if (cfile[i].getType() != E_humrec_data) {
         continue;
      }
      index = -1;
      string = NULL;
      for (j=0; j<cfile[i].getFieldCount(); j++) {
         if (strcmp(cfile[i].getExInterp(j), "**index") == 0) {
            sscanf(cfile[i][j], "%d", &index);
            if (index < 0 || index > 25) {
               index = -1;  // ignore out-of range indices
            }
         }
         if (strcmp(cfile[i].getExInterp(j), "**pixel") == 0) {
            if (sscanf(cfile[i][j], "%d %d %d", &test1, &test2, &test3) == 3) {
               string = cfile[i][j]; 
            } else {
               cerr << "Error in color file on line " << i + 1 
                    << ":" << cfile[i] << endl;
               exit(1);
            }
         }
      }
      if (index >= 0) {
         colorindex[index] = string;
         // note that the "string" can not be deallocated until the
         // program is finished.
      }
   }
}



//////////////////////////////
//
// setFilterOptions --
//

void setFilterOptions(Array& channelfilter, const char* exclude) {
   int length = strlen(exclude);
   int character;
   int i;
   int value;

   for (i=0; i<length; i++) {
      character = toupper(exclude[i]);
      if (!isxdigit(character)) {
         continue;
      }
      if (isdigit(character)) {
         value = character - '0';
      } else {
         value = character - 'A' + 10;
      }
      if (value >= 0 && value <= 15) {
         channelfilter[value] = 0;
      }
   }
}



//////////////////////////////
//
// example --
//

void example(void) {


}



//////////////////////////////
//
// usage --
//

void usage(const char* command) {

}



// md5sum: 9e18757243079343f0a83e7267595491 mkeyscape.cpp [20080227]