//
// Programmer: Craig Stuart Sapp <craig@ccrma.stanford.edu>
// Creation Date: Thu Dec 25 22:45:39 PST 2008
// Last Modified: Thu Dec 25 22:45:43 PST 2008
// Filename: ...sig/examples/all/noisesequence.cpp
// Web Address: http://sig.sapp.org/examples/museinfo/humdrum/noisesequence.cpp
// Syntax: C++; museinfo
//
// Description: Create a dissimilarity plot between two sequences.
//
#include "humdrum.h"
#include <stdlib.h> /* for qsort and bsearch functions */
#include <time.h> /* for random number seeding */
#include <math.h>
// function declarations:
void checkOptions(Options& opts, int argc, char** argv);
void example(void);
void usage(const char* command);
double getMean(Array<double>& data);
double getSampleSD(double mean, Array<double>& data);
double pearsonCorrelation(int size, double* x, double* y);
void printData(Array<double>& a, Array<double>& b,
int counter);
void getSequences(Array<double>& a, Array<double>& b,
int seqlength, double noiseadd);
// User interface variables:
Options options;
int seqlength = 100; // used with -l option
double corrtarget = 0.0; // used with -r option
double corrtolerance = 0.00001; // used with -t option
int limit = 1000; // used with -L option
int randomseed = 0; // used with -s option
double noiseadd = 0.0; // used witn -n option
//////////////////////////////////////////////////////////////////////////
int main(int argc, char** argv) {
// process the command-line options
checkOptions(options, argc, argv);
if (randomseed <= 0) {
randomseed = time(NULL); // time in seconds
}
srand48(randomseed);
Array<double> a;
Array<double> b;
getSequences(a, b, seqlength, noiseadd);
double corr = pearsonCorrelation(seqlength, a.getBase(), b.getBase());
int counter = 1;
while (counter < limit) {
if (fabs(corr - corrtarget) < corrtolerance) {
break;
}
getSequences(a, b, seqlength, noiseadd);
corr = pearsonCorrelation(seqlength, a.getBase(), b.getBase());
counter++;
}
printData(a, b, counter);
return 0;
}
//////////////////////////////////////////////////////////////////////////
//////////////////////////////
//
// getSequence --
//
void getSequences(Array<double>& a, Array<double>& b, int seqlength,
double noiseadd) {
a.setSize(seqlength);
a.allowGrowth(0);
b.setSize(seqlength);
b.allowGrowth(0);
int i;
for (i=0; i<seqlength; i++) {
a[i] = drand48();
b[i] = drand48();
}
if (noiseadd != 0.0) {
for (i=0; i<seqlength; i++) {
b[i] += noiseadd * a[i];
}
}
}
//////////////////////////////
//
// printData --
//
void printData(Array& a, Array& b, int counter) {
int i;
int size = a.getSize();
double corr = pearsonCorrelation(seqlength, a.getBase(), b.getBase());
cout << "!!!correlation:\t" << corr << "\n";
cout << "!!!randomseed:\t" << randomseed << "\n";
cout << "!!!target:\t" << corrtarget << "\n";
if (fabs(corr - corrtarget) > corrtolerance) {
cout << "!!!status:\tFAILED to match tolerance\n";
}
cout << "!!!tolerance:\t" << corrtolerance << "\n";
cout << "!!!attempts:\t" << counter << "\n";
cout << "!!!length:\t" << seqlength << "\n";
if (noiseadd != 0.0) {
cout << "!!!noiseadd:\t" << noiseadd << "\n";
}
cout << "**seq\t**seq\n";
for (i=0; i<size; i++) {
cout << a[i] << "\t" << b[i] << "\n";
}
cout << "*-\t*-\n";
}
//////////////////////////////
//
// getSampleSD --
//
double getSampleSD(double mean, Array& data) {
int size = data.getSize();
double sum = 0.0;
double value;
int i;
for (i=0; i<size; i++) {
value = data[i] - mean;
sum += value * value;
}
return sqrt(sum / (size - 1.0));
}
//////////////////////////////
//
// getMean --
//
double getMean(Array& data) {
int size = data.getSize();
if (size <= 0) {
return 0.0;
}
int i;
double sum = 0.0;
for (i=0; i<size; i++) {
sum += data[i];
}
return sum / size;
}
//////////////////////////////
//
// checkOptions --
//
void checkOptions(Options& opts, int argc, char* argv[]) {
opts.define("l|length=i:100", "length of output sequences");
opts.define("r|correlation=d:0.0", "target correlation");
opts.define("t|tolerance=d:0.00001", "tolerance of correlation target");
opts.define("L|limit=i:1000000", "limit to the number of attempts");
opts.define("s|seed=i:0", "random number generator seed value");
opts.define("n|noiseadd=d:0.0", "add fraction of seq A to seq B");
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, Dec 2008" << endl;
exit(0);
} else if (opts.getBoolean("version")) {
cout << argv[0] << ", version: 25 Dec 2008" << endl;
cout << "compiled: " << __DATE__ << endl;
cout << MUSEINFO_VERSION << endl;
exit(0);
} else if (opts.getBoolean("help")) {
usage(opts.getCommand());
exit(0);
} else if (opts.getBoolean("example")) {
example();
exit(0);
}
seqlength = opts.getInteger("length");
corrtarget = opts.getDouble("correlation");
corrtolerance = opts.getDouble("tolerance");
limit = opts.getInteger("limit");
randomseed = opts.getInteger("seed");
noiseadd = opts.getDouble("noiseadd");
}
//////////////////////////////
//
// pearsonCorrelation --
//
double pearsonCorrelation(int size, double* x, double* y) {
double sumx = 0.0;
double sumy = 0.0;
double sumco = 0.0;
double meanx = x[0];
double meany = y[0];
double sweep;
double deltax;
double deltay;
int i;
for (i=2; i<=size; i++) {
sweep = (i-1.0) / i;
deltax = x[i-1] - meanx;
deltay = y[i-1] - meany;
sumx += deltax * deltax * sweep;
sumy += deltay * deltay * sweep;
sumco += deltax * deltay * sweep;
meanx += deltax / i;
meany += deltay / i;
}
double popsdx = sqrt(sumx / size);
double popsdy = sqrt(sumy / size);
double covxy = sumco / size;
return covxy / (popsdx * popsdy);
}
//////////////////////////////
//
// example --
//
void example(void) {
}
//////////////////////////////
//
// usage --
//
void usage(const char* command) {
}
// md5sum: 512b899a94c1e175683ce9b5b055bad5 noisesequence.cpp [20081229]