```//
// 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 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);
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;
const char* filename = "";

if (options.getArgCount() <= 0) {
filename = "<STDIN>";

} else {
filename = options.getArg(1);

#ifndef OLDCPP
infile.open(filename, ios::in);
#else
infile.open(filename, ios::in | ios::nocreate);
#endif
}

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);
}

switch (plotshape) {
case SHAPE_RECTANGLE:
break;
case SHAPE_TRIANGLE:
default:
}
} 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) {
const char* 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 >& correlation) {
const char* 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 >& 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 >& 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;

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 >& correlations, Array& 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 >& 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& 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;
}
}

//////////////////////////////
//
// 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());
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")) {
}

if (opts.getBoolean("rectangle"))  {
plotshape = SHAPE_RECTANGLE;
}
}

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

void makeSineArch(Array& 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& 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;
}

}

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

Array<Array<double> >& correlations) {

int i, j;
int colsize;

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

for (j=0; j<colsize; 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 >& 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& 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& 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: 7ce7a3769a2f96495b317272ac15c3ca hicor.cpp [20080518]
```