//
// Programmer:    Craig Stuart Sapp <craig@ccrma.stanford.edu>
// Creation Date: Sat Jul  1 02:48:55 PDT 2006
// Last Modified: Tue Jul  4 06:53:48 PDT 2006
// Last Modified: Mon Aug 21 04:06:54 PDT 2006 (added polycorrelation plots)
// Last Modified: Sat May 12 04:49:15 PDT 2007 (added gradient coloring)
// Last Modified: Sat May 12 09:23:11 PDT 2007 (added smoothing)
// Last Modified: Sun May 13 00:07:45 PDT 2007 (added data plotting)
// Last Modified: Sat Aug 18 21:54:01 PDT 2007 (added bi-feature plotting)
// Filename:      ...sig/examples/all/hicor.cpp
// Web Address:   http://sig.sapp.org/examples/museinfo/hicor/hicor.cpp
// Syntax:        C++; museinfo
//
// Description:   Generate a hierarchical correlation plot from a
//                pair of number sequences.  
//

#include "Array.h"
#include "Options.h"
#include "PixelColor.h"

#include <string.h>
#include <stdio.h>
#include <math.h>

#ifndef OLDCPP
   using namespace std;
   #include <iostream>
   #include <fstream>
#else
   #include <iostream.h>
   #include <fstream.h>
#endif


#define COLOR_BW  0
#define COLOR_HUE 1

#define SHAPE_TRIANGLE  0
#define SHAPE_RECTANGLE 1


// function declarations:
void     checkOptions          (Options& opts, int argc, char** argv);
void     example               (void);
void     usage                 (const char* command);
void     calculateAverage      (Array<Array<double> >& correlations, 
                                Array<double>& a);
void     calculateCorrelations (Array<Array<double> >& correlations,
                                Array<double>& a, 
                                Array<double>& b);
double   pearsonCorrelation    (int size, double* a, double* b);
double   getMean               (int size, double* a);
void     printImageTriangle    (Array<Array<double> >& correlation);
void     printImageRectangle   (Array<Array<double> >& correlation);
void     printCorrelationPixel (double value, int style);
double   average               (int size, double* data);
void     normalizeAverageData  (Array<Array<double> >& correlations);
double   getStandardDeviation  (double mean, Array<double>& data);
int      logScale              (int level, int maxlevel, double factor);
void     calculateAverageDifference(Array<Array<double> >& correlations, 
                                Array<double>& a, 
                                Array<double>& b);

void     makeSineArch          (Array<double>& arch, int size);
void     makeRamp              (Array<double>& ramp, int size);
void     calculateArchCorrelations(Array<Array<double> >& correlations, 
                                Array<double>& a);
void     calculateRampCorrelations(Array<Array<double> >& correlations, 
                                Array<double>& a);
void     calculateGradient      (Array<Array<double> >& gradient,
                                 Array<Array<double> >& correlations);
int      getGradient            (Array<Array<double> >& scape, 
                                 int row, int col);
void     smoothSequence         (Array<double>& sequence, double gain);
void     unsmoothSequence       (Array<double>& sequence, double gain);
void     printDataPlot          (Array<double>& a, Array<double>& b, 
                                 int height, int gridlines);
void     printBiCorrelationPixel(double value, double valueB, int style);
void     printBiImageTriangle   (Array<Array<double> >& correlation,
                                 Array<Array<double> >& correlationB);
void     printBiImageRectangle  (Array<Array<double> >& correlation,
                                 Array<Array<double> >& correlationB);

// User interface variables:
Options   options;
int       dataq     = 0;         // debug: display input data only with -d
int       correlq   = 0;         // debug: display correlations with -c
int       colortype = COLOR_BW;  // used with various color options
int       plotshape = SHAPE_TRIANGLE;  // used with -r option
int       averageq  = 0;         // used with -a option
int       diffq     = 0;         // display average differences with -d option
int       condiffq  = 0;         // display contoured differences, --cd option
int       column    = 0;         // used with -n option (for averaging)
int       logq      = 0;         // used with -g option
double    logfactor = 1500;      // used with -f option
int       flipq     = 0;         // flip the data from tempo to duration
int       polyq     = 0;         // used with -p option
int       archq     = 0;         // used with -A option
int       rampq     = 0;         // used with -R option
int       maxheight = -1;        // used with --height option
int       gradientq = 0;         // used with -G option
double    smooth    = 0.0;       // used with -S option
int       hline     = 0;         // used with -H option
int       plotq     = 0;         // used with -P option
int       unsmoothq = 0;         // used with -U option
int       bifeature = 0;         // used with -b option
int       bgval     = 80;        // no interface to user yet

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

int main(int argc, char** argv) {
   // process the command-line options
   checkOptions(options, argc, argv);

   double aval;
   double bval;

   Array<double> a;
   Array<double> b;

   a.setSize(100000);
   a.setGrowth(100000);
   a.setSize(0);

   b.setSize(100000);
   b.setGrowth(100000);
   b.setSize(0);

   ifstream infile;
   string filename = "";

   if (options.getArgCount() <= 0) {
      filename = "<STDIN>";
   
   } else {
      filename = options.getArg(1);

      infile.open(filename.c_str(), ios::in);
   }

   if (!infile.is_open()) {
      cerr << "Cannot open file for reading: " << filename << endl;
      exit(1);
   }

   char templine[4096];

   if (averageq) {

      while (!infile.eof()) {
         infile.getline(templine, 4096, '\n');
         if (infile.eof() && (strcmp(templine, "") == 0)) {
            break;
         } else if (sscanf(templine, "%lf\t%lf", &aval, &bval) == 2) {
            if (flipq) {
               aval = 60.0 / aval;
               bval = 60.0 / bval;
            }
            a.append(aval);
            b.append(bval);
         } else if (sscanf(templine, "%lf", &aval) == 1) {
            if (flipq) {
               aval = 60.0 / aval;
            }
            a.append(aval);
            b.append(aval);
         }
      }

   } else {

      while (!infile.eof()) {
         infile.getline(templine, 4096, '\n');
         if (infile.eof() && (strcmp(templine, "") == 0)) {
            break;
         } else {
            if (sscanf(templine, "%lf\t%lf", &aval, &bval) == 2) {
               if (flipq) {
                  aval = 60.0 / aval;
                  bval = 60.0 / bval;
               }
               a.append(aval);
               b.append(bval);
            } else if (sscanf(templine, "%lf", &aval) == 1) {
               if (flipq) {
                  aval = 60.0 / aval;
               }
               a.append(aval);
               b.append(aval);
            }
         }
      }

   }

   if (smooth != 0.0) {
      if (!unsmoothq) {
         smoothSequence(a, smooth);
         smoothSequence(b, smooth);
      } else {
         unsmoothSequence(a, smooth);
         unsmoothSequence(b, smooth);
      }
   }

   int i, j;
   if (dataq) {
      for (i=0; i<a.getSize(); i++) {
         cout << a[i] << "\t" << b[i] << endl;
      }
      exit(0);
   } 

   Array<Array<double> > correlations;
   Array<Array<double> > correlationsB;

   correlations.setSize(a.getSize());
   correlationsB.setSize(a.getSize());
   for (i=0; i<a.getSize(); i++) {
      correlations[i].setSize(i+1);
      correlations[i].allowGrowth(0);
      correlationsB[i].setSize(i+1);
      correlationsB[i].allowGrowth(0);
   }

   if (averageq) {
      if (bifeature) {
         calculateAverage(correlations,  a);
         calculateAverage(correlationsB, b);
      } else if (column == 0) {
         calculateAverage(correlations, a);
      } else {
         calculateAverage(correlations, b);
      }
   } else if (archq) {
      calculateArchCorrelations(correlations, a);
      if (bifeature) {
         calculateArchCorrelations(correlationsB, b);
      }
   } else if (rampq) {
      calculateRampCorrelations(correlations, a);
      if (bifeature) {
         calculateRampCorrelations(correlationsB, b);
      }
   } else if (condiffq || diffq) {
      calculateAverageDifference(correlations, a, b);
   } else {
      calculateCorrelations(correlations, a, b);
   }

   if (correlq) {
      for (i=0; i<correlations.getSize(); i++) {
         for (j=0; j<correlations[i].getSize(); j++) {
            cout << correlations[i][j];
            if (j < correlations[i].getSize()-1) {
               cout << '\t';
            }
         }
         cout << '\n';
      }
      exit(0);
   }

   Array<Array<double> > gradient;
   if (gradientq) {
      calculateGradient(gradient, correlations);

      switch (plotshape) {
         case SHAPE_RECTANGLE:
            printImageRectangle(gradient);
            break;
         case SHAPE_TRIANGLE:
         default:
            printImageTriangle(gradient);
      }
   } else {
      if (bifeature) {
         switch (plotshape) {
            case SHAPE_RECTANGLE:
               printBiImageRectangle(correlations, correlationsB);
               break;
            case SHAPE_TRIANGLE:
            default:
               printBiImageTriangle(correlations, correlationsB);
         }
      } else {
         switch (plotshape) {
            case SHAPE_RECTANGLE:
               printImageRectangle(correlations);
               break;
            case SHAPE_TRIANGLE:
            default:
               printImageTriangle(correlations);
         }
      }
   }

   if (plotq) {
      printDataPlot(a, b, plotq, hline);
   }

   return 0;
}


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


//////////////////////////////
//
// printDataPlot -- print an input data plot underneath the scape.
//

void printDataPlot(Array<double>& a, Array<double>& b, int height,
      int gridlines) {
   PixelColor pc;
   pc.setColor(255,0,0);

   int i, j;
   int rows = height;
   int cols = a.getSize() * 2;

   Array<Array<PixelColor> > plotregion;
   plotregion.setSize(rows);
   plotregion.allowGrowth(0);
   for (i=0; i<rows; i++) {
      plotregion[i].setSize(cols);
      plotregion[i].allowGrowth(0);
      for (j=0; j<cols; j++) {
         plotregion[i][j].setColor(255,255,255);
      }
   }

   // draw vertical lines if gridlines > 0:
   if (gridlines > 0) {
      for (j=0; j<cols/2; j++) {
         if ((j+1) % gridlines == 0) {
            for (i=0; i<rows; i++) {
               plotregion[i][j*2].setColor(bgval,bgval,bgval);
               plotregion[i][j*2+1].setColor(bgval,bgval,bgval);
            }
         }
      }
   }
   

   double minn = a[0];
   double maxx = a[0];
   for (i=1; i<a.getSize(); i++) {
      if (a[i] < minn) { minn = a[i]; }
      if (a[i] > maxx) { maxx = a[i]; }
   }

   double i1;
   double i2;
   int x1;
   int x2;
   int k;
   int xlow;
   int xhi;
   for (j=0; j<a.getSize(); j++) {
      i1 = (a[j] - minn)/(maxx - minn) * rows;
      i = rows - int(i1) - 1;
      if (i < 0) { i = 0; }
      if (i >= rows) { i = rows-1; }
      plotregion[i][2*j].setColor(0,0,255);
      if (j < a.getSize()-1) {
         i2 = (a[j+1] - minn)/(maxx - minn) * rows;
         x1 = rows - int(i1) - 1;
         x2 = rows - int(i2) - 1;
         if (x1 < 0) { x1 = 0; }
         if (x1 >= rows) { x1 = rows-1; }
         if (x2 < 0) { x2 = 0; }
         if (x2 >= rows) { x2 = rows-1; }
         if (x1 < x2) {
            xlow = x1;
            xhi  = x2;
         } else {
            xlow = x2;
            xhi  = x1;
         }
         for (k=xlow+1; k<xhi; k++) {
            plotregion[k][2*j+1].setColor(175,175,255);
         }
	 if (xlow == xhi) {
            plotregion[xlow][2*j+1].setColor(175,175,255);
         }
	 if (xhi-xlow == 1) {
            plotregion[x1][2*j+1].setColor(175,175,255);
         }
      }
   }

   for (i=0; i<rows; i++) {
      for (j=0; j<cols; j++) {
         cout << " ";
         plotregion[i][j].writePpm3(cout);
      }
      cout << endl;
   }

}



//////////////////////////////
//
// printBiImageTriangle --
//

void printBiImageTriangle(Array<Array<double> >& correlation,
      Array<Array<double> >& correlationB) {
   string background = options.getString("background");

   int maxcolumn = correlation.getSize() * 2;
   int maxrow    = (correlation.getSize()-1) * 2;

   cout << "P3\n";
   cout << maxcolumn << " " << maxrow << "\n";
   cout << "255\n";

   int i, j;
   int ii;
   int datawidth;
   int bgwidth;
   for (i=0; i<correlation.getSize()-1; i++) {
      ii = i;
      if (logq) {
         ii = logScale(i, correlation.getSize(), logfactor);
      }
      datawidth = correlation[ii].getSize();
      bgwidth = (maxcolumn - datawidth * 2)/2;

      for (j=0; j<bgwidth; j++) {
         cout << " " << background << " ";
      }
      for (j=0; j<datawidth; j++) {
         printBiCorrelationPixel(correlation[ii][j], 
              correlationB[ii][j], colortype);
      }
      for (j=0; j<bgwidth; j++) {
         cout << " " << background << " ";
      }
      cout << "\n";

      // repeat the row to double it (to make image square)
      for (j=0; j<bgwidth; j++) {
         cout << " " << background << " ";
      }
      for (j=0; j<datawidth; j++) {
         printBiCorrelationPixel(correlation[ii][j], 
               correlationB[ii][j], colortype);
      }
      for (j=0; j<bgwidth; j++) {
         cout << " " << background << " ";
      }
      cout << "\n";
   }
}



//////////////////////////////
//
// printImageTriangle --
//

void printImageTriangle(Array<Array<double> >& correlation) {
   string background = options.getString("background");

   int maxcolumn = correlation.getSize() * 2;
   int maxrow    = (correlation.getSize()-1) * 2;

   cout << "P3\n";
   cout << maxcolumn << " " << maxrow << "\n";
   cout << "255\n";

   int i, j;
   int ii;
   int datawidth;
   int bgwidth;
   for (i=0; i<correlation.getSize()-1; i++) {
      ii = i;
      if (logq) {
         ii = logScale(i, correlation.getSize(), logfactor);
      }
      datawidth = correlation[ii].getSize();
      bgwidth = (maxcolumn - datawidth * 2)/2;

      for (j=0; j<bgwidth; j++) {
         cout << " " << background << " ";
      }
      for (j=0; j<datawidth; j++) {
         printCorrelationPixel(correlation[ii][j], colortype);
      }
      for (j=0; j<bgwidth; j++) {
         cout << " " << background << " ";
      }
      cout << "\n";

      // repeat the row to double it (to make image square)
      for (j=0; j<bgwidth; j++) {
         cout << " " << background << " ";
      }
      for (j=0; j<datawidth; j++) {
         printCorrelationPixel(correlation[ii][j], colortype);
      }
      for (j=0; j<bgwidth; j++) {
         cout << " " << background << " ";
      }
      cout << "\n";
   }
}


/*  Rebuild plot by pre-calulating the image and then print it:


//////////////////////////////
//
// printImageTriangle --
//

void printImageTriangle(Array<Array<double> >& correlation) {
   PixelColor background;
   background.setColor(options.getString("background"));

   int maxcolumn = correlation.getSize() * 2;
   int maxrow    = (correlation.getSize()) * 2;

   int startrow = 0;
   if (maxheight > 0 && maxheight < correlation.getSize()) {
      startrow = correlation.getSize() - maxheight;
      maxrow = (maxheight) * 2;
   }

   cout << "P3\n";
   cout << maxcolumn << " " << maxrow + plotq << "\n";
   cout << "255\n";

   int i, j;

   Array<Array<PixelColor> > mainplot;
   mainplot.setSize(maxrow);
   mainplot.allowGrowth(0);
   for (i=0; i<mainplot.getSize(); i++) {
      mainplot[i].setSize(maxcolumn);
      mainplot[i].allowGrowth(0);
      for (j=0; j<mainplot[i].getSize(); j++) {
         mainplot[i][j] = background;
      }
   }

   int ii;
   int datawidth;
   int bgwidth;
   int bgline = 0;
   for (i=startrow; i<correlation.getSize(); i++) {
      if (hline > 0) {
         if ((correlation.getSize() - i - 1 + 1) % hline == 0) {
            bgline = 1;
         } else {
            bgline = 0;
         }
      }

      ii = i;
      if (logq) {
         ii = logScale(i, correlation.getSize(), logfactor);
      }
      datawidth = correlation[ii].getSize();
      bgwidth = (maxcolumn - datawidth * 2)/2;

xxx
      for (j=0; j<bgwidth; j++) {
         if (bgline == 1) mainplot[i][j].setColor(bgval,bgval,bgval);
      }
      for (j=0+bgwidth; j<datawidth+bgwidth; j++) {
         printCorrelationPixel(correlation[ii][j], colortype);
      }
      for (j=0+bgwidth+datawidth; j<bgwidth+datawidth+bgwidth; j++) {
         if (bgline == 1) mainplot[ii][j].setColor(bgval,bgval,bgval);
      }

      // repeat the row to double it (to make image square)      
      for (j=0; j<bgwidth; j++) {
         if (bgline == 0) cout << " " << background << " ";
         if (bgline == 1) mainplot[ii][j].setColor(bgval,bgval,bgval);
      }
      for (j=0; j<datawidth; j++) {
         printCorrelationPixel(correlation[ii][j], colortype);
      }
      for (j=0; j<bgwidth; j++) {
         if (bgline == 0) cout << " " << background << " ";
         if (bgline == 1) mainplot[ii][j].setColor(bgval,bgval,bgval);
      }
   }

}

*/


//////////////////////////////
//
// logScale -- 
//

int logScale(int level, int maxlevel, double factor) {
   level = maxlevel - level - 1;
   double dummy = factor;
   double minval = log10(10.0);
   double valrange = log10(10.0 + dummy) - minval;
   int    templevel = maxlevel - level;
   double tempval = (log10(10.0 + dummy * templevel/maxlevel)-minval)/valrange;
   int    newlevel = (int)(tempval * maxlevel + 0.5);
   // newlevel = maxlevel - newlevel;
   newlevel = newlevel - 1;
   return newlevel;
}



//////////////////////////////
//
// printBiImageRectangle --
//

void printBiImageRectangle(Array<Array<double> >& correlation,
      Array<Array<double> >& correlationB) {

   int maxcolumn = correlation.getSize() * 2;
   int maxrow    = (correlation.getSize()-1) * 2;
   if (maxheight > 0 && maxheight < maxrow) {
      maxrow = maxheight;
   }

   cout << "P3\n";
   cout << maxcolumn << " " << maxrow << "\n";
   cout << "255\n";

   int i, j;
   int datawidth;
   int jstretch;

   for (i=0; i<correlation.getSize()-1; i++) {
      datawidth = correlation[i].getSize();

      for (j=0; j<maxcolumn/2; j++) {
         jstretch = int(double(j)/maxcolumn*2*datawidth);
         printBiCorrelationPixel(correlation[i][jstretch], 
             correlationB[i][jstretch], colortype);
      }
      cout << "\n";

      // repeat for double pixel size
      for (j=0; j<maxcolumn/2; j++) {
         jstretch = int(double(j)/maxcolumn*2*datawidth);
         printBiCorrelationPixel(correlation[i][jstretch], 
             correlationB[i][jstretch], colortype);
      }
      cout << "\n";

   }
}



//////////////////////////////
//
// printImageRectangle --
//

void printImageRectangle(Array<Array<double> >& correlation) {

   int maxcolumn = correlation.getSize() * 2;
   int maxrow    = (correlation.getSize()-1) * 2;
   if (maxheight > 0 && maxheight < maxrow) {
      maxrow = maxheight;
   }

   cout << "P3\n";
   cout << maxcolumn << " " << maxrow << "\n";
   cout << "255\n";

   int i, j;
   int datawidth;
   int jstretch;

   for (i=0; i<correlation.getSize()-1; i++) {
      datawidth = correlation[i].getSize();

      for (j=0; j<maxcolumn/2; j++) {
         jstretch = int(double(j)/maxcolumn*2*datawidth);
         printCorrelationPixel(correlation[i][jstretch], colortype);
      }
      cout << "\n";

      // repeat for double pixel size
      for (j=0; j<maxcolumn/2; j++) {
         jstretch = int(double(j)/maxcolumn*2*datawidth);
         printCorrelationPixel(correlation[i][jstretch], colortype);
      }
      cout << "\n";

   }
}



 
//////////////////////////////
//
// printBiCorrelationPixel --
//

void printBiCorrelationPixel(double value, double valueB, int style) {
   // input values are in the range from -1.0 to +1.0
   value = (value+1.0)/2.0;
   valueB = (valueB+1.0)/2.0;

   // brighten the colors a bit:
   value = sqrt(value);
   valueB = sqrt(valueB);

   PixelColor pc;

   pc.setGreen(0);   // green channel is not used (red/green color blind people)

   // gradient style is ignored for now.
   // black and white style is ignored.

   pc.setRedF(value);
   pc.setBlueF(valueB);

   cout << " ";
   pc.writePpm3(cout);
   cout << " ";
	   
   // repeat the cell for a square image
   cout << " ";
   pc.writePpm3(cout);
   cout << " ";

}

 
//////////////////////////////
//
// printCorrelationPixel --
//

void printCorrelationPixel(double value, int style) {

   PixelColor pc;


   if (!gradientq) {
      switch (style) {
         case COLOR_HUE:
            pc.setHue(-((value+1)/2.0)*0.80 - 0.20); // -1=purple, 
	                                             // 0=green +1=red
            //pc.setHue(value);	 // -1=purple, 0=green +1=red
            break;
         case COLOR_BW:
         default:
            pc.setGrayNormalized((value + 1.0)/2.0);
      }
   
      if (diffq) {
         if (value > 0) {
            pc.setColor(255,0,0);
         } else if (value < 0) {
            pc.setColor(0,0,255);
         } else {
            pc.setColor(255,255,255);
         }
      } else if (condiffq) {
         if      (value >   0.210)  {  pc.setColor(252,41,41);   } // red
         else if (value >   0.150)  {  pc.setColor(252,41,210);  } // pink
         else if (value >   0.090)  {  pc.setColor(255,144,41);  } // orange
         else if (value >   0.030)  {  pc.setColor(223,250,35);  } // yellow
         else if (value >  -0.030)  {  pc.setColor(239,239,239); } // whiteish
         else if (value >  -0.090)  {  pc.setColor(148,247,138); } // green
         else if (value >  -0.150)  {  pc.setColor(82,246,221);  } // lt blue
         else if (value >  -0.210)  {  pc.setColor(36,59,227);   } // dk blue
         else if (value <= -0.210)  {  pc.setColor(80,4,184);    } // violet
      }
   } else {
      // print gradient colors
      switch (int(value+0.1)) {
         case 6:  pc.setColor(80,4,184);    break; // violet
         case 5:  pc.setColor(36,59,227);   break; // dk blue
         case 4:  pc.setColor(82,246,221);  break; // lt blue
         case 3:  pc.setColor(148,247,138); break; // green
         case 2:  pc.setColor(223,250,35);  break; // yellow
         case 1:  pc.setColor(255,144,41);  break; // orange
         case 0:  pc.setColor(252,41,41);   break; // red
         default: pc.setColor(255,255,255); break; // white (invalid option)
      }
      pc.invert();
      pc = pc / 4.0;
      pc.invert();
      if (int(value+0.1) == 0) { pc.setColor(255,0,0); }
   }

   cout << " ";
   pc.writePpm3(cout);
   cout << " ";
	   
   // repeat the cell for a square image
   cout << " ";
   pc.writePpm3(cout);
   cout << " ";

}



//////////////////////////////
//
// calculateAverage --
//

void calculateAverage(Array<Array<double> >& correlations, Array<double>& a) {

   int i, j;
   double *aa;
   int rowsize;
   int corelsize;

   for (i=0; i<correlations.getSize()-1; i++) {
      rowsize = correlations[i].getSize();
      corelsize = a.getSize() - i;
      for (j=0; j<rowsize; j++) {
         aa = a.getBase() + j;
         correlations[i][j] = average(corelsize, aa);
      }
   }

   // set the bottom row to data
   i = correlations.getSize()-1;
   for (j=0; j<correlations[i].getSize(); j++) {
      correlations[i][j] = a[j];
   }

   // normalized the data into the range from -1 to 1 so that
   // the colorization of the pixels will work.

   normalizeAverageData(correlations);

}



//////////////////////////////
//
// normalizeAverageData -- transform the data range into -1 to +1 
// using the sigmoid function positions at the mean of the bottom row
// (i.e. the top item), and transition with equal to 4 standard 
// deviations (-2s to +2s) of the bottom row.
//

#define SIGMOID(x,c,w) (1.0/(1.0+exp(-(((x)-(c))/((w)/8)))))

void normalizeAverageData(Array<Array<double> >& correlations) {

   int i, j;
   double mean = 0.0;
   for (i=0; i<correlations[correlations.getSize()-1].getSize(); i++) {
      mean += correlations[correlations.getSize()-1][i];
   }
   mean = mean / correlations.getSize();

   double sd = getStandardDeviation(mean,correlations[correlations.getSize()-1]);

   double& center = mean;
   double  width  = 4 * sd;

   for (i=0; i<correlations.getSize(); i++) {
      for (j=0; j<correlations[i].getSize(); j++) {
         correlations[i][j] = SIGMOID(correlations[i][j], center, width)*2-1;
      }
   } 
}


//////////////////////////////
//
// getStandardDeviation --
//

double getStandardDeviation(double mean, Array<double>& data) {
   double sum = 0.0;
   double value;
   int i;
   for (i=0; i<data.getSize(); i++) {
      value = data[i] - mean;
      sum += value * value;
   }
   
   return sqrt(sum / data.getSize());
}



//////////////////////////////
//
// calculateAverageDifference --
//

void calculateAverageDifference(Array<Array<double> >& correlations, 
      Array<double>& a, Array<double>& b) {

   int i, j;
   double *aa;
   double *bb;
   double aval = 1.0;
   double bval = 1.0;
   int rowsize;
   int corelsize;

   for (i=0; i<correlations.getSize(); i++) {
      rowsize = correlations[i].getSize();
      corelsize = a.getSize() - i;
      for (j=0; j<rowsize; j++) {
         aa = a.getBase() + j;
         bb = b.getBase() + j;
         aval = average(corelsize, aa);
         bval = average(corelsize, bb);
         correlations[i][j] = bval - aval;    // switch to b-a later

         if (condiffq) {
            // segment by tempo bands.
            correlations[i][j] = correlations[i][j] / aval;
         }
      }

   }


   // normalized the data into the range from -1 to 1 so that
   // the colorization of the pixels will work.

   //normalizeAverageData(correlations);


}



//////////////////////////////
//
// calculateCorrelations --
//

void calculateCorrelations(Array<Array<double> >& correlations,
      Array<double>& a, Array<double>& b) {

   int i, j;
   double *aa;
   double *bb;
   int rowsize;
   int corelsize;

   for (i=0; i<correlations.getSize()-1; i++) {
      rowsize = correlations[i].getSize();
      corelsize = a.getSize() - i;
      for (j=0; j<rowsize; j++) {
         aa = a.getBase() + j;
         bb = b.getBase() + j;
         correlations[i][j] = pearsonCorrelation(corelsize, aa, bb);
      }
   }

   // set the bottom row to 1.0 correlations
   i = correlations.getSize()-1;
   for (j=0; j<correlations[i].getSize(); j++) {
      correlations[i][j] = 1.0;
   }

}



//////////////////////////////
//
// average -- calculate the average of the input values.
//

double average(int size, double* data) {
   double sum = 0.0;

   int i;
   for (i=0; i<size; i++) {
      sum += data[i];
   }

   return sum / size;
}



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

double pearsonCorrelation(int size, double* a, double* b) {
   double meana = getMean(size, a);
   double meanb = getMean(size, b);

   double topsum     = 0.0;
   double bottomasum = 0.0;
   double bottombsum = 0.0;

   int i;
   for (i=0; i<size; i++) {
      topsum += (a[i] - meana) * (b[i] - meanb);
      bottomasum += (a[i] - meana) * (a[i] - meana);
      bottombsum += (b[i] - meanb) * (b[i] - meanb);
   }

   if (bottomasum == 0.0 || bottombsum == 0.0) {
      return 0.0;
   }
   return topsum / sqrt(bottomasum * bottombsum);
}



//////////////////////////////
//
// getMean --
//

double getMean(int size, double* a) {
   if (size <= 0) {
      return 0.0;
   }

   int i;
   double sum = 0.0;
   for (i=0; i<size; i++) {
      sum += a[i];
   }

   return sum / size;
}



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

void checkOptions(Options& opts, int argc, char* argv[]) {
   opts.define("g|log=b", "display vertical axis on a log scale");
   opts.define("flip=b", "flip the input data from tempo to duration");
   opts.define("d|difference=b", "display average difference plots");
   opts.define("A|arch=b", "calculate arch correlations");
   opts.define("R|ramp=b", "calculate ramp correlations");
   opts.define("cd|condiff=b", "display contoured difference plots");
   opts.define("f|logfactor=d:1500.0", "factor for strenth of log scaling");
   opts.define("n|column=i:1", "which column to average");
   opts.define("a|average=b", "display data average rather than correlation");
   opts.define("x|data=b",    "display input data and then quit");
   opts.define("c|correlation=b", "display correlation data and then quit");
   opts.define("h|hue=b", "display correlations in hue coloring");
   opts.define("r|rectangle=b", "display correlations in rectangular plot");
   opts.define("p|poly=b", "poly correlation plot");
   opts.define("b|bg|background=s:255 255 255", "background color");
   opts.define("bi=b", "Binary feature correlation");
   opts.define("height=i:-1", "Maximum height scope");
   opts.define("G|gradient=b", "Gradient analysis (of A or R option)");
   opts.define("S|smooth=d:0.0", "Smoothing gain (0.0 = don't smooth)");
   opts.define("U|unsmooth=b", "apply unsmoothing");
   opts.define("H|horizontal=i:0", "hoizontal line step");
   opts.define("P|plot=i:0", "display original data underneath");

   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, Jul 2006" << endl;
      exit(0);
   } else if (opts.getBoolean("version")) {
      cout << argv[0] << ", version: 13 May 2007" << endl;
      cout << "compiled: " << __DATE__ << endl;
      exit(0);
   } else if (opts.getBoolean("help")) {
      usage(opts.getCommand().c_str());
      exit(0);
   } else if (opts.getBoolean("example")) {
      example();
      exit(0);
   }

   logq      = opts.getBoolean("log");
   logfactor = opts.getDouble("logfactor");
   smooth    = opts.getDouble("smooth");
   unsmoothq = opts.getBoolean("unsmooth");
   hline     = opts.getInteger("horizontal");
   flipq     = opts.getBoolean("flip");
   averageq  = opts.getBoolean("average");
   diffq     = opts.getBoolean("difference");
   condiffq  = opts.getBoolean("condiff");
   column    = opts.getInteger("column") - 1;
   polyq     = opts.getBoolean("poly");
   archq     = opts.getBoolean("arch");
   maxheight = opts.getInteger("height");
   rampq     = opts.getBoolean("ramp");
   plotq     = opts.getInteger("plot");
   bifeature = opts.getBoolean("bi");
   if (plotq < 0) {
      plotq = 0;
   }
   if (column < 0) {
      column = 0;
   } else if (column > 1) {
      column = 1;
   }
   dataq     = opts.getBoolean("data");
   correlq   = opts.getBoolean("correlation");
   if (opts.getBoolean("hue")) {
      colortype = COLOR_HUE;
   }

   if ((rampq || archq) && opts.getBoolean("gradient")) {
      gradientq = 1;
   }
   
   if (opts.getBoolean("rectangle"))  {
      plotshape = SHAPE_RECTANGLE;
   }
}



//////////////////////////////
//
// makeSineArch --  generate a sqrt sine half-cycle arch of the
//     specified length.
//

void makeSineArch(Array<double>& arch, int size) {
   arch.setSize(size);
   arch.allowGrowth(0);

   if (size == 1) {
      arch[0] = 0;
      return;
   }

   int i;
   double freq = 2.0 * M_PI / ((size-1) * 2.0);
   for (i=0; i<size-1; i++) {
      arch[i] = sin(freq * i);
   }
   arch[size-1] = 0.0;

}



//////////////////////////////
//
// makeRamp --
//

void makeRamp(Array<double>& ramp, int size) {
   ramp.setSize(size);
   ramp.allowGrowth(0);

   if (size == 1) {
      ramp[0] = 0;
   }

   double factor = 1.0 / (size - 1.0);
   int i;
   for (i=0; i<size; i++) {
      ramp[i] = factor * i;
   }

}



//////////////////////////////
//
// calculateArchCorrelations --
//

void calculateArchCorrelations(Array<Array<double> >& correlations, 
      Array<double>& a) {

   Array<double> arch;

   int i, j;
   double *aa;
   double *bb;
   int rowsize;
   int corelsize;

   for (i=0; i<correlations.getSize()-1; i++) {
      makeSineArch(arch, i+1);
      rowsize = correlations[i].getSize();
      corelsize = a.getSize() - i;
      makeSineArch(arch, corelsize);
      for (j=0; j<rowsize; j++) {
         aa = a.getBase() + j;
         bb = arch.getBase();
         correlations[i][j] = pearsonCorrelation(corelsize, aa, bb);
      }
   }

   // set the bottom row to 1.0 correlations
   i = correlations.getSize()-1;
   for (j=0; j<correlations[i].getSize(); j++) {
      correlations[i][j] = 1.0;
   }

}



//////////////////////////////
//
// calculateGradient --
//

void calculateGradient(Array<Array<double> >& gradient, 
      Array<Array<double> >& correlations) {

   int i, j;
   int colsize;

   gradient.setSize(correlations.getSize());
   gradient.allowGrowth(0);

   for (i=0; i<correlations.getSize(); i++) {
      colsize = correlations[i].getSize();
      gradient[i].setSize(colsize);
      gradient[i].allowGrowth(0);

      for (j=0; j<colsize; j++) {
         gradient[i][j] = getGradient(correlations, i, j);
      }
   }
}



//////////////////////////////
//
// getGradient -- search for the uphill direction in the neighborhood
//     of a pixel in a scape plot. 
//   0 = pixel is at top of surface (red)
//   1 = right-hand pixel is upward slope (orange)
//   2 = top-right-hand pixel is upward slope (yellow)
//   3 = top-left-hand pixel is upward slope (green)
//   4 = left-hand pixel is upward slope (light blue)
//   5 = bottom-left hand pixel is upward slope (dark blue)
//   6 = bottom-right hand pixel is upward slope (violet)
//

int getGradient(Array<Array<double> >& scape, int row, int col) {
   int rows = scape.getSize();
   int cols = scape[row].getSize();

   int output = 0;

   double ax = scape[row][col];

/* When the rows are indexedfrom the bottom:
   // to the right
   if (col + 1 < cols) {
      if (ax < scape[row][col+1]) { output = 1; }
   }
   // to the upper right
   if (row + 1 < rows && col < cols-1) { 
      if (ax < scape[row+1][col]) { output = 2; }
   }
   // to the upper left
   if (row + 1 < rows && col - 1 >= 0) {  
      if (ax < scape[row+1][col-1]) { output = 3; }
   }
   // to the left
   if (col - 1 >= 0) {  
      if (ax < scape[row][col-1]) { output = 4; }
   }
   // to the lower left
   if (row - 1 >= 0) {  
      if (ax < scape[row-1][col]) { output = 5; }
   }
   // to the lower right
   if (row - 1 >= 0 && col + 1 < cols+1) {  
      if (ax < scape[row-1][col]+1) { output = 6; }
   }
*/

   // when the rows are indexed from the top:

   // to the right
   if (col + 1 < cols) {
      if (ax < scape[row][col+1]) { output = 1; ax = scape[row][col+1]; }
   }
   // to the upper right
   if (row - 1 >= 0 && col < cols-1) { 
      if (ax < scape[row-1][col]) { output = 2; ax = scape[row-1][col]; }
   }
   // to the upper left
   if (row - 1 >= 0 && col - 1 >= 0) {  
      if (ax < scape[row-1][col-1]) { output = 3; ax = scape[row-1][col-1]; }
   }
   // to the left
   if (col - 1 >= 0) {  
      if (ax < scape[row][col-1]) { output = 4; ax = scape[row][col-1]; }
   }
   // to the lower left
   if (row + 1 < rows) {  
      if (ax < scape[row+1][col]) { output = 5; ax = scape[row+1][col]; }
   }
   // to the lower right
   if (row + 1 < rows && col + 1 < cols+1) {  
      if (ax < scape[row+1][col+1]) { output = 6; ax = scape[row+1][col+1]; }
   }

   return output;
}



//////////////////////////////
//
// smoothSequence -- smooth the sequence with a
//    symmetric exponential smoothing filter (applied in the forward
//    and reverse directions with the specified input gain.
//
//    Difference equation for smoothing: y[n] = k * x[n] + (1-k) * y[n-1]
//

void smoothSequence(Array<double>& sequence, double gain) {
   double oneminusgain = 1.0 - gain;
   int i;
   int ssize = sequence.getSize();

   // reverse filtering first 
   for (i=ssize-2; i>=0; i--) {
      sequence[i] = gain*sequence[i] + oneminusgain*sequence[i+1];
   }

   // then forward filtering
   for (i=1; i<ssize; i++) {
      sequence[i] = gain*sequence[i] + oneminusgain*sequence[i-1];
   }

}

//////////////////////////////
//
// unsmoothSequence -- removed smoothed sequence from original sequence.
//

void unsmoothSequence(Array<double>& sequence, double gain) {
   Array<double> smoothed = sequence;
   smoothSequence(smoothed, gain);
   int i;
   for (i=0; i<sequence.getSize(); i++) {
      sequence[i] -= smoothed[i];
   } 
}



//////////////////////////////
//
// calculateRampCorrelations --
//

void calculateRampCorrelations(Array<Array<double> >& correlations, 
      Array<double>& a) {

   Array<double> ramp;

   int i, j;
   double *aa;
   double *bb;
   int rowsize;
   int corelsize;

   for (i=0; i<correlations.getSize()-1; i++) {
      rowsize = correlations[i].getSize();
      corelsize = a.getSize() - i;
      makeRamp(ramp, corelsize);
      for (j=0; j<rowsize; j++) {
         aa = a.getBase() + j;
         bb = ramp.getBase();
         correlations[i][j] = pearsonCorrelation(corelsize, aa, bb);
      }
   }

   // set the bottom row to 1.0 correlations
   i = correlations.getSize()-1;
   for (j=0; j<correlations[i].getSize(); j++) {
      correlations[i][j] = 1.0;
   }
}




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

void example(void) {


}



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

void usage(const char* command) {

}



// md5sum: cd1143e80ebe077a00e0469fa901b449 hicor.cpp [20160320]