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