//
// Programmer:    Craig Stuart Sapp <craig@ccrma.stanford.edu>
// Creation Date: Sun Sep 30 23:06:34 PDT 2007
// Last Modified: Sun Sep 30 23:06:40 PDT 2007
// Last Modified: Sat Apr 12 21:00:21 PDT 2008 (added -a and -A options)
// Filename:      ...sig/examples/all/smoother.cpp
// Web Address:   http://sig.sapp.org/examples/museinfo/humdrum/smoother.cpp
// Syntax:        C++; museinfo
//
// Description:   Smoothes data with exponential smoothing filter.
//

#include "humdrum.h"
#include <stdlib.h>

#ifndef OLDCPP
   using namespace std;
#endif

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

// function declarations
void      checkOptions       (Options& opts, int argc, char* argv[]);
void      example            (void);
void      usage              (const char* command);
void      processFileBasic   (HumdrumFile& infile);
void      processFileSmoothAll(HumdrumFile& infile);
void      processFileDesmoothAll(HumdrumFile& infile);
void      fillDesmooth       (Array<double>& sequence, 
                              int spine, HumdrumFile& infile);
void      fillSmooth         (Array<double>& sequence, 
                              int spine, HumdrumFile& infile);

// global variables
Options   options;            // database for command-line arguments
double    sgain = 0.5;        // used with -s option (smoothing amount)
double    allsmoothQ  = 0;    // used with -a option
double    allresidueQ = 0;    // used with -A option

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

int main(int argc, char* argv[]) {
   HumdrumFile infile;

   // process the command-line options
   checkOptions(options, argc, argv);

   string filename;
   infile.clear();
   // if no command-line arguments read data file from standard input
   int numinputs = options.getArgCount();
   if (numinputs < 1) {
      infile.read(cin);
   } else {
      filename = options.getArg(1);
      infile.read(filename.c_str());
   }
   if (allsmoothQ) {
      processFileSmoothAll(infile);
   } else if (allresidueQ) {
      processFileDesmoothAll(infile);
   } else {
      processFileBasic(infile);
   }

}

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


///////////////////////////////
//
// processFileSmoothAll --
//

void processFileSmoothAll(HumdrumFile& infile) {
   Array<Array<double> > spines;
   spines.setSize(infile.getMaxTracks());
   int i;
   for (i=0; i<spines.getSize(); i++) {
      fillSmooth(spines[i], i, infile);
   }

   Array<int> pointers;
   pointers.setSize(spines.getSize());
   pointers.allowGrowth(0);
   pointers.setAll(0);   
   int jsize;

   int j;
   for (i=0; i<infile.getNumLines(); i++) {
      if (infile[i].getType() == E_humrec_data) {
         jsize = infile[i].getFieldCount();
         for (j=0; j<jsize; j++) {
            if (strcmp(infile[i][j], ".") == 0) {
               cout << infile[i][j];
               if (j < jsize-1) {
                  cout << "\t";
               }
               continue;
            }
            if (infile[i][j][0] == '=') {
               cout << infile[i][j];
               if (j < jsize-1) {
                  cout << "\t";
               }
               continue;
            }
            cout << spines[j][pointers[j]++];
            if (j < jsize-1) {
               cout << "\t";
            }
            
         }
         cout << "\n";
      } else {
         cout << infile[i] << "\n";
      }
   }

}



///////////////////////////////
//
// processFileDesmoothAll --
//

void processFileDesmoothAll(HumdrumFile& infile) {
   Array<Array<double> > spines;
   spines.setSize(infile.getMaxTracks());
   int i;
   for (i=0; i<spines.getSize(); i++) {
      fillDesmooth(spines[i], i, infile);
   }

   Array<int> pointers;
   pointers.setSize(spines.getSize());
   pointers.allowGrowth(0);
   pointers.setAll(0);   
   int jsize;

   int j;
   for (i=0; i<infile.getNumLines(); i++) {
      if (infile[i].getType() == E_humrec_data) {
         jsize = infile[i].getFieldCount();
         for (j=0; j<jsize; j++) {
            if (strcmp(infile[i][j], ".") == 0) {
               cout << infile[i][j];
               if (j < jsize-1) {
                  cout << "\t";
               }
               continue;
            }
            if (infile[i][j][0] == '=') {
               cout << infile[i][j];
               if (j < jsize-1) {
                  cout << "\t";
               }
               continue;
            }
            cout << spines[j][pointers[j]++];
            if (j < jsize-1) {
               cout << "\t";
            }
            
         }
         cout << "\n";
      } else {
         cout << infile[i] << "\n";
      }
   }

}



//////////////////////////////
//
// fillSmooth --
//

void fillSmooth(Array<double>& sequence, int spine, HumdrumFile& infile) {

   sequence.setSize(infile.getNumLines());
   sequence.setSize(0);
   double value;

   int i;
   for (i=0; i<infile.getNumLines(); i++) {
      if (infile[i].getType() == E_humrec_data) {
         if (strcmp(infile[i][spine], ".") == 0) {
            continue;
         }
         if (infile[i][spine][0] == '=') {
            continue;
         }
         value = atof(infile[i][spine]);
         sequence.append(value);
      }
   }

   int thesize = sequence.getSize();
   Array<double> smoothed;
   smoothed.setSize(thesize);

   double filterk = sgain;
   double oneminusk = 1.0 - filterk;

   // reverse filtering first
   smoothed[thesize-1] = sequence[thesize-1];
   for (i=thesize-2; i>=0; i--) {
      smoothed[i] = filterk*sequence[i] + oneminusk*smoothed[i+1];
   }

   // then forward filtering
   for (i=1; i<thesize; i++) {
      smoothed[i] = filterk*smoothed[i] + oneminusk*smoothed[i-1];
   }

   for (i=0; i<thesize; i++) {
      sequence[i] = smoothed[i];
   }
}



//////////////////////////////
//
// fillDesmooth --
//

void fillDesmooth(Array<double>& sequence, int spine, HumdrumFile& infile) {
   sequence.setSize(infile.getNumLines());
   sequence.setSize(0);
   double value;

   int i;
   for (i=0; i<infile.getNumLines(); i++) {
      if (infile[i].getType() == E_humrec_data) {
         if (strcmp(infile[i][spine], ".") == 0) {
            continue;
         }
         if (infile[i][spine][0] == '=') {
            continue;
         }
         value = atof(infile[i][spine]);
         sequence.append(value);
      }
   }

   int thesize = sequence.getSize();
   Array<double> smoothed;
   smoothed.setSize(thesize);

   double filterk = sgain;
   double oneminusk = 1.0 - filterk;

   // reverse filtering first
   smoothed[thesize-1] = sequence[thesize-1];
   for (i=thesize-2; i>=0; i--) {
      smoothed[i] = filterk*sequence[i] + oneminusk*smoothed[i+1];
   }

   // then forward filtering
   for (i=1; i<thesize; i++) {
      smoothed[i] = filterk*smoothed[i] + oneminusk*smoothed[i-1];
   }

   for (i=0; i<thesize; i++) {
      sequence[i] = sequence[i] - smoothed[i];
   }
}



//////////////////////////////
//
// processFileBasic -- only read the first row for processing for now...
//

void processFileBasic(HumdrumFile& infile) {
   Array<double> sequence;
   sequence.setSize(infile.getNumLines());
   sequence.setSize(0);

   double value;
   int i;
   for (i=0; i<infile.getNumLines(); i++) {
      if (infile[i].getType() == E_humrec_data) {
         value = atof(infile[i][0]);
         sequence.append(value);
      }
   }

   int thesize = sequence.getSize();
   Array<double> smoothed;
   smoothed.setSize(thesize);

   double filterk = sgain;
   double oneminusk = 1.0 - filterk;
   

   // reverse filtering first
   smoothed[thesize-1] = sequence[thesize-1];
   for (i=thesize-2; i>=0; i--) {
      smoothed[i] = filterk*sequence[i] + oneminusk*smoothed[i+1];
   }

   // then forward filtering
   for (i=1; i<thesize; i++) {
      smoothed[i] = filterk*smoothed[i] + oneminusk*smoothed[i-1];
   }



   // print original, smoothed, and unsmoothed data in columns:

   for (i=0; i<thesize; i++) {
      cout << sequence[i] << "\t" << smoothed[i] << "\t";
      cout << sequence[i] - smoothed[i];
      cout << "\n";
   }

}



//////////////////////////////
//
// checkOptions -- validate and process command-line options.
//

void checkOptions(Options& opts, int argc, char* argv[]) {
   opts.define("s|smooth=d:0.5");       // smoothing factor
   opts.define("a|all|allsmooth=b");    // smoothing factor
   opts.define("A|All|allresidue=b");   // smoothing factor
   opts.define("debug=b");              // determine bad input line num
   opts.define("author=b");             // author of program
   opts.define("version=b");            // compilation info
   opts.define("example=b");            // example usages
   opts.define("h|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, Oct 2007" << endl;
      exit(0);
   } else if (opts.getBoolean("version")) {
      cout << argv[0] << ", version: 1 Oct 2007" << 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);
   }

   sgain = opts.getDouble("smooth");
   allsmoothQ  = opts.getBoolean("allsmooth");
   allresidueQ = opts.getBoolean("allresidue");
}



//////////////////////////////
//
// example -- example usage of the program
//

void example(void) {
   cout <<
   "                                                                         \n"
   << endl;
}



//////////////////////////////
//
// usage -- gives the usage statement for the program
//

void usage(const char* command) {
   cout <<
   "                                                                         \n"
   << endl;
}



// md5sum: 57147157ce5389485452772f4f9a9730 smoother.cpp [20160320]