//
// Programmer:    Craig Stuart Sapp <craig@ccrma.stanford.edu>
// Creation Date: Sun Jul 28 21:48:39 PDT 2002
// Last Modified: Mon Jul 29 13:54:33 PDT 2002
// Filename:      ...sig/examples/all/iwrange.cpp
// Web Address:   http://sig.sapp.org/examples/museinfo/humdrum/iwrange.cpp
// Syntax:        C++; museinfo
//
// Description:   Find the independent variation in interval weights
//	which yeild as good or better results than the current value.
// 

#include "humdrum.h"

#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define MINIM 0.00000001

typedef Array<double> ArrayDouble;
typedef Array<int>    ArrayInt;

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

// function declarations
void   checkOptions             (Options& opts, int argc, char* argv[]);
void   example                  (void);
void   usage                    (const char* command);
void   getStartingWeights       (Array<double>& weights, 
                                 HumdrumFile& weightfile);
void   getChordInformation      (Array<ArrayInt>& pitches, Array<int>& root, 
                                 Array<double>& count, HumdrumFile& datafile);
double getErrors                (Array<double>& weights, 
                                 Array<ArrayInt>& pitches, Array<int>& root, 
                                 Array<double>& count);
void   makeRandomWeights        (Array<ArrayDouble>& stepweights, double range);
void   printWeights             (Array<double> weights);
double getDistance              (Array<double>& a, Array<double>& b);
void   findRange                (Array<ArrayInt>& pitches, Array<int>& root, 
                                 Array<double>& count, 
                                 Array<double>& initialweights, 
                                 int index, double& tune, double step);
void   printKern                (Array<double>& initialweights, 
                                 Array<double>& min, Array<double>& max);
void   printPlot                (Array<double>& iw, Array<double>& min,
                                 Array<double>& max);
double getMidRange              (Array<double>& min, Array<double>& max,
                                 Array<ArrayInt>& pitches, Array<int>& root, 
                                 Array<double>& count);


// global variables
Options      options;            // database for command-line arguments
int          debugQ     = 0;     // used with --debug option
int          verboseQ   = 0;     // used with -v option
double       decay      = 0.5;   // decay of range between each step
double       limit      = 10.0;  // maximum distance which can be searched
double       step       = 1.0;   // maximum distance which can be searched
int          plotQ      = 0;     // used with -p option
int          bestQ      = 0;     // used with the -b option

double       baseerror  = 0;
double       minerror   = -1;
double       miderror   = 0;
Array<double> bestweights;


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

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

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

   datafile.read(options.getArg(1));
   weightfile.read(options.getArg(2));

   Array<double> count;          // the frequency of the chords occurence
   Array<int> root;              // the root of the chord
   Array<ArrayInt> pitches;      // the pitch class set of the chord
   Array<double> initialweights; // starting weights of the search

   getStartingWeights(initialweights, weightfile);
   getChordInformation(pitches, root, count, datafile);
   bestweights = initialweights;

   int i;
   Array<double> min(initialweights.getSize());
   Array<double> max(initialweights.getSize());
   for (i=0; i<initialweights.getSize(); i++) {
      if (i==5||i==11||i==22||i==28||i==34) {
         continue;
      }
      findRange(pitches, root, count, initialweights, i, min[i], -step);
      findRange(pitches, root, count, initialweights, i, max[i], step);
   }
     

   miderror = getMidRange(min, max, pitches, root, count);
   if (plotQ) {
      printPlot(initialweights, min, max);
   } else {
      printKern(initialweights, min, max);
   }

   return 0;
}



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


//////////////////////////////
//
// getMidRange --
//

double getMidRange(Array<double>& min, Array<double>& max,
      Array<ArrayInt>& pitches, Array<int>& root, Array<double>& count) {
   Array<double> weights(40);
   weights.setAll(10000);
   int i;
   for (i=0; i<40; i++) {
      if (i==5||i==11||i==22||i==28||i==34) {
         continue;
      }
      weights[i] = (max[i] - min[i])/2.0 + min[i];
   }
   return getErrors(weights, pitches, root, count);
}



//////////////////////////////
//
// printPlot --
//

void printPlot(Array<double>& iw, Array<double>& min,
      Array<double>& max) {

   printKern(iw, min,max);
   cout << "@START: FIGURE\n\n";
   int yvals[40] = {
        -14, -7, 0, 7, 14,  // C
        1000,
        -12, -5, 2, 9, 16,  // D
        1000,
        -10, -3, 4, 11, 18,  // E
        -15, -8, -1, 6, 13,  // F
        1000,
        -13, -6, 1, 8, 15,  // G
        1000,
        -11, -4, 3, 10, 17,  // A
        1000,
        -9, -2, 5, 12, 19  // B
   };
   double radius[40] = {
        -0.4, -0.4, -0.4, -0.4, -0.4, 
        0.0,
        0.4, 0.4, 0.4, 0.4, 0.4, 
        0.0,
        -0.3, -0.3, -0.3, -0.3, -0.3, 
        0.3, 0.3, 0.3, 0.3, 0.3, 
        0.0,
        -0.2, -0.2, -0.2, -0.2, -0.2, 
        0.0,
        0.2, 0.2, 0.2, 0.2, 0.2, 
        0.0,
        -0.1, -0.1, -0.1, -0.1, -0.1
   };

   int i;
   for (i=0; i<40; i++) {
      if (i==5||i==11||i==22||i==28||i==34) {
         continue;
      }
      cout << "@START: PLOTDATA\n";
      cout << min[i] << "\t" << yvals[i] << "\n";
      cout << max[i] << "\t" << yvals[i] << "\n";
      cout << "@END: PLOTDATA\n\n";
   }

   double curstyle = -1;
   double laststyle = -1;
   double currad = -1;
   cout << "@START: POINTS\n";
   for (i=0; i<40; i++) {
      if (i==5||i==11||i==22||i==28||i==34) {
         continue;
      }
      if (currad != fabs(radius[i])) {
         currad = fabs(radius[i]);
         cout << "@RADIUS: " << currad << "\n";
      }
      laststyle = curstyle;
      curstyle = radius[i] < 0 ? 0 : 1;
      if (curstyle != laststyle) {
         if (curstyle == 1) {
            cout << "@STYLE: circle\n";
         } else {
            cout << "@STYLE: opencircle\n";
         }
      }
      cout << iw[i] << "\t" << yvals[i] << "\n";
   }

   cout <<"@END: FIGURE\n\n";
}



//////////////////////////////
//
// printKern --
//

void printKern(Array<double>& initialweights, Array<double>& min,
      Array<double>& max) {
   char buffer[1024] = {0};

   cout << "!! Base Errors:\t"    << baseerror << endl;
   cout << "!! Best Errors:\t"     << minerror << endl;
   cout << "!! MidRange Errors:\t" << miderror << endl;
   cout << "**kern\t**weight\t**min\t**max\t**range\t**midrange\n";
   int i;

   for (i=0; i<initialweights.getSize(); i++) {
      if (i==5||i==11||i==22||i==28||i==34) {
         continue;
      }
      cout << Convert::base40ToKern(buffer, i+3*40) << "\t"
           << initialweights[i] << "\t"
           << min[i] << "\t" << max[i] << "\t"
           << max[i] - min[i] << "\t"
           << (max[i] - min[i])/2 + min[i] << "\n";
   }

   cout << "*-\t*-\t*-\t*-\t*-\t*-\n";

   if (!bestQ) {
      return;
   }

   if (bestweights == initialweights) {
      cout << "!! no better weights found\n";
      return;
   }
   cout << "\n";
   cout << "!! best weight set found:\n";
   cout << "**kern\t**weight\n";

   for (i=0; i<bestweights.getSize(); i++) {
      if (i==5||i==11||i==22||i==28||i==34) {
         continue;
      }
      cout << Convert::base40ToKern(buffer, i+3*40) << "\t"
           << bestweights[i] << "\n";
   }

   cout << "*-\t*-\n";

}




//////////////////////////////
//
// findRange --
//

void findRange(Array<ArrayInt>& pitches, Array<int>& root, 
      Array<double>& count, Array<double>& initialweights, 
      int index, double& tune, double step) {

   Array<double> weights;
   weights = initialweights;
   baseerror = getErrors(weights, pitches, root, count);
   if (minerror < 0) {
      minerror = baseerror;
   }

   double starttune = tune;
   tune = weights[index];
   double newerrors = 0.0;
   while ((fabs(step) > MINIM) && (fabs(tune - starttune) < limit*2)) {
      newerrors = getErrors(weights, pitches, root, count);
      if (newerrors < minerror) {
         minerror = newerrors;
         bestweights = weights;
      }
      if (newerrors <= baseerror) {
         tune += step;
      } else {
         tune -= step;
         step *= decay;
         if (verboseQ) {
            cerr << "Step set to " << step << endl;
         }
      }
      weights[index] = tune;
      if (verboseQ) {
         cerr << "weight = " << tune << endl;
      }
   }

   char buffer[32] = {0};
   if (verboseQ) {
      cerr << Convert::base40ToKern(buffer, index+3*40) 
           << "=====================================" << endl;
   }

}



//////////////////////////////
//
// getErrors -- try all chords and count how many root errors occured.
//

double getErrors(Array<double>& weights, Array<ArrayInt>& pitches,
      Array<int>& root, Array<double>& count) {

   Array<double> rootscores;
   rootscores.setSize(40);
   rootscores.allowGrowth(0);
   int i, j, m;
   double errors = 0;
   int min;
   for (m=0; m<root.getSize(); m++) {
      rootscores.setAll(0.0);
      for (i=0; i<rootscores.getSize(); i++) {
         for (j=0; j<pitches[m].getSize(); j++) {
            rootscores[i] += weights[(pitches[m][j]-i+400)%40];
         }
      }
      min = 0;
      for (i=0; i<rootscores.getSize(); i++) {
         if (rootscores[min] > rootscores[i]) {  
            min = i;
         }
      }
      if (root[m] != min+2) {
         if (debugQ) {
            cout << "Error: root=" << root[m] 
                 << "\tbut measured: " << min << endl;
         }
         errors += count[m];
      }
   }

   return errors;
}



//////////////////////////////
//
// getDistance -- find the Euclidian distance between two vectors.
//

double getDistance(Array<double>& a, Array<double>& b)  {
   Array<double> c(40);
   int i;
   double sum = 0.0;
   for (i=0; i<c.getSize(); i++) {
      if (i==5 || i==11 || i==22 || i==28 || i == 34) {
         continue;
      }
      c[i] = a[i] - b[i];
      sum += c[i] * c[i];
   }
   
   return sqrt(sum);
}



//////////////////////////////
//
// makeRandomWeights --
//

void makeRandomWeights(Array<ArrayDouble>& stepweights, double range) {
   int i, j;
   double delta = 0.0;
   for (i=1; i<stepweights.getSize(); i++) {
      for (j=0; j<stepweights[i].getSize(); j++) {
         delta = drand48() * 2 * range - range;
         stepweights[i][j] = stepweights[0][j] + delta;
      }
   }
}



//////////////////////////////
//
// getStartingWeights --
//

void getStartingWeights(Array<double>& weights, HumdrumFile& weightfile) {
   weights.setSize(40);
   weights.setAll(100000);
   weights.allowGrowth(0);

   int i, j;
   int root;
   double weight;
   for (i=0; i<weightfile.getNumLines(); i++) {
      root = -1;
      weight = 10000;
      if (!weightfile[i].isData()) {
         continue;
      }
      for (j=0; j<weightfile[i].getFieldCount(); j++) {
         if (strcmp("**kern", weightfile[i].getExInterp(j)) == 0) {
            root = Convert::kernToBase40(weightfile[i][j]) % 40;
         } else if (strcmp("**weight", weightfile[i].getExInterp(j)) == 0) {
            weight = strtod(weightfile[i][j], NULL);
         }
      }
      if (root >= 0) {
         weights[root] = weight;
      }
   }

   if (debugQ) {
      printWeights(weights);
   }
}


//////////////////////////////
//
// printWeights --
//

void printWeights(Array<double> weights) {
   char buffer[128] = {0};
   cout << "**kern\t**weight\n";
   int i;
   for (i=0; i<weights.getSize(); i++) {
      if (i==5 || i==11 || i==22 || i==28 || i==34) {
         continue;   // invalid interval
      }
      cout << Convert::base40ToKern(buffer, i+4*40);
      cout << "\t" << weights[i] << "\n";
   }
   cout << "*-\t*-\n";
}



//////////////////////////////
//
// getChordInformation --
//

void getChordInformation(Array<ArrayInt>& pitches, Array<int>& root, 
     Array<double>& count, HumdrumFile& datafile) {

   count.setSize(datafile.getNumLines());
   count.setSize(0);
   root.setSize(datafile.getNumLines());
   root.setSize(0);
   pitches.setSize(datafile.getNumLines());
   pitches.setSize(0);

   char buffer[1024] = {0};
   int i, j, k;
   int troot = -1;
   int tpitch = -1;
   double tcount = 0.0;
   Array<int> tpitches;
   tpitches.setSize(100);
   tpitches.setSize(0);
   for (i=0; i<datafile.getNumLines(); i++) {
      tpitches.setSize(0);
      troot = -1;
      tcount = 0.0;
      for (j=0; j<datafile[i].getFieldCount(); j++) {
         if (strcmp("**root", datafile[i].getExInterp(j)) == 0) {
            troot = Convert::kernToBase40(datafile[i][j]) % 40;
         } else if (strcmp("**count", datafile[i].getExInterp(j)) == 0) {
            tcount = strtod(datafile[i][j], NULL);
         } else if (strcmp("**kern", datafile[i].getExInterp(j)) == 0) {
            int notecount = datafile[i].getTokenCount(j);
            for (k=0; k<notecount; k++) { 
               datafile[i].getToken(buffer, j, k);
               tpitch = Convert::kernToBase40(buffer);
               tpitches.append(tpitch);
            }
         }
      }

      if (troot < 0 || tcount <= 0.0 || tpitches.getSize() == 0) {
         continue;
      }
      root.append(troot);
      count.append(tcount);
      pitches.append(tpitches);
   }


   if (debugQ) {
      cout << "**count\t**root\t**kern\n";
      for (i=0; i<count.getSize(); i++) {
         cout << count[i] << "\t";
         cout << Convert::base40ToKern(buffer, root[i]+3*40) << "\t";
         for (j=0; j<pitches[i].getSize(); j++) {
            cout << Convert::base40ToKern(buffer, pitches[i][j]);
            if (j < pitches[i].getSize() - 1) {
               cout << " ";
            }
         }
         cout << "\n";
      }
      cout << "*-\t*-\t*-\n";
   }

}



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

void checkOptions(Options& opts, int argc, char* argv[]) {
   opts.define("d|decay=d:0.5",   "decay of range between each step");
   opts.define("l|limit=d:10.0",  "number of trials with same errors");
   opts.define("s|step=d:1.0",    "number of trials with same errors");
   opts.define("v|verbose=b",     "monitor status");   
   opts.define("p|plot=b",        "output format as a plot");   
   opts.define("b|best=b", "output best single parameter modified weights");

   opts.define("debug=b",         "trace input parsing");   
   opts.define("author=b",        "author of the program");   
   opts.define("version=b",       "compilation information"); 
   opts.define("example=b",       "example usage"); 
   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, July 2002" << endl;
      exit(0);
   } else if (opts.getBoolean("version")) {
      cout << argv[0] << ", version: 28 July 2002" << 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);
   }

   decay      = opts.getDouble("decay");
   limit      = opts.getDouble("limit");
   step       = opts.getDouble("step");
   debugQ     = opts.getBoolean("debug");
   verboseQ   = opts.getBoolean("verbose");
   plotQ      = opts.getBoolean("plot");
   bestQ      = opts.getBoolean("best");

}



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

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



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

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



// md5sum: ffafec806053581db9fe92bcf9d46cdb iwrange.cpp [20050403]