```//
// Programmer:    Craig Stuart Sapp <craig@ccrma.stanford.edu>
// Creation Date: Sun Jul 28 21:48:39 PDT 2002
// Filename:      ...sig/examples/all/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>

#include <iomanip>

#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);
void   printMmaPlot(Array<double>& iw, Array<double>& min,
Array<double>& max,
Array<ArrayInt> pitches,
Array<int> root,
Array<double> count);
double getMidRange(Array<double>& min, Array<double>& max,
Array<ArrayInt>& pitches, Array<int>& root,
Array<double>& count);
void   printLineInfo(int cc, int ii, Array<double>& iw,
double minval, double maxval,
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
int          mmaQ       = 0;     // used with -M optnion

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

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 if (mmaQ) {
printMmaPlot(initialweights, min, max, pitches, root, count);
} 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);
}

//////////////////////////////
//
// printMmaPlot --
//

void printMmaPlot(Array<double>& iw, Array<double>& min,
Array<double>& max, Array<ArrayInt> pitches,
Array<int> root, Array<double> count) {

int cc = 100;
double fifthname = 0;
char buffer[100] = {0};

cout << "(* Weights: line-of-fifths interval, base-40 pitch class, ";
cout << "interval abbreviation, weight *)\n";
cout << "weights = {\n";
int i;
for (i=0; i<iw.getSize(); i++) {
fifthname = Convert::base40IntervalToLineOfFifths((i-2+40)%40);
if (fifthname > 40) {
continue;
}
Convert::base40ToIntervalAbbr2(buffer, (i-2+40)%40);
cout << "{";
cout << fifthname;
cout << ", ";
cout << i;
cout << ", \"";
cout << buffer;
cout << "\", ";
cout << iw[i];
cout << "}";
if (i < iw.getSize()-1) {
cout << ",";
}
cout << "\n";
}
cout << "};\n\n";

cout << "(* Data: Line-of-Fifths interval, interval abbreviation, ";
cout << "interval weight, minimum, maximum, {improvement, base errors, ";
cout << "total chords}, {{weight, relative error},{...}} *)\n";
cout << "data = {\n";

for (i=0; i<iw.getSize(); i++) {
fifthname = Convert::base40IntervalToLineOfFifths((i-2+40)%40);
if (fifthname > 40) {
continue;
}
Convert::base40ToIntervalAbbr2(buffer, (i-2+40)%40);
cout << "{";
cout << fifthname << ", ";
cout << "\"" << buffer << "\", ";
cout << fixed << iw[i] << ", ";
cout << fixed << min[i] << ", ";
cout << fixed << max[i] << ", ";
printLineInfo(cc, i, iw, min[i], max[i], pitches, root, count);
cout << "}";
if (i < iw.getSize()-1) {
cout << ",\n";
}
}
cout << "};\n";
}

//////////////////////////////
//
// printLineInfo --
//

void printLineInfo(int cc, int ii, Array<double>& iw, double minval,
double maxval, Array<ArrayInt> pitches, Array<int> root,
Array<double> count) {

double sum = 0.0;
int i;
for (i=0; i<count.getSize(); i++) {
sum += count[i];
}

double baseerror = getErrors(iw, pitches, root, count);
double localerror;

Array<double> tempweights = iw;
double minimumerror = baseerror;

Array<double> results;
results.setSize(cc+1);
results.allowGrowth(0);

for (i=0; i<=cc; i++) {
tempweights[ii] = minval + i * (maxval - minval) / cc;
localerror = getErrors(tempweights, pitches, root, count);
results[i] = localerror;
if (localerror < minimumerror) {
minimumerror = localerror;
}
}

// print error rates and improvements
cout << "{";
cout << int(baseerror - minimumerror + 0.5);
cout << ", ";
cout << int(baseerror + 0.5);
cout << ", ";
cout << int (sum + 0.5);
cout << "}, ";

cout << "{";
for (i=0; i<=cc; i++) {
if (i==0 || i==cc || results[i] != results[i-1] || results[i] != results[i+1]) {
cout << "{";
cout << fixed << minval + i * (maxval - minval) / cc;
cout << ", ";
if (baseerror - minimumerror > 0.0) {
cout << (results[i] - minimumerror) / (baseerror - minimumerror);
} else {
cout << 1;
}
cout << "}";
if (i<cc) {
cout << ", ";
}
}
}
cout << "}";

}

//////////////////////////////
//
// 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
};
-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;
cout << "@START: POINTS\n";
for (i=0; i<40; i++) {
if (i==5||i==11||i==22||i==28||i==34) {
continue;
}
}
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];
while ((fabs(step) > MINIM) && (fabs(tune - starttune) < limit*2)) {
newerrors = getErrors(weights, pitches, root, count);
bestweights = weights;
}
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& a, Array& 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& 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& 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 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 ((datafile[i].getType() == E_humrec_none) ||
(datafile[i].getType() == E_humrec_empty) ||
(datafile[i].getType() == E_humrec_global_comment) ||
(datafile[i].getType() == E_humrec_bibliography) ) {
continue;
}

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("M|mma=b", "Mathematica plotting data");

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");
mmaQ       = opts.getBoolean("mma");

}

//////////////////////////////
//
// 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]
```