//
// 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<double>& a, Array<double>& 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<double>& 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<double>& 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().c_str());
      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: c88c206a3ead0beda9745cb938faa150 noisesequence.cpp [20160320]