+/**********************************************************************************************************************
+ vector< vector<double> > NMDSCommand::calculateStressGradientVector(vector<seqDist>& eDists, vector<seqDist>& D, double rawStress, double stress, vector< vector<double> >& axes) {
+ try {
+ vector< vector<double> > gradient; gradient.resize(dimension);
+ for (int i = 0; i < gradient.size(); i++) { gradient[i].resize(axes[0].size(), 0.0); }
+
+ double sumDij = 0.0;
+ for (int i = 0; i < eDists.size(); i++) { sumDij += (eDists[i].dist * eDists[i].dist); }
+
+ for (int i = 0; i < eDists.size(); i++) {
+
+ for (int j = 0; j < dimension; j++) {
+
+ if (m->control_pressed) { return gradient; }
+
+ double firstTerm1 = (stress / rawStress) * (eDists[i].dist - D[i].dist);
+ double firstTerm2 = (stress / sumDij) * eDists[i].dist;
+ double firstTerm = firstTerm1 - firstTerm2;
+
+ float r = (dimension-1.0);
+ double temp = 1.0 / (pow(eDists[i].dist, r));
+ float absTemp = abs(axes[j][eDists[i].seq1] - axes[j][eDists[i].seq2]);
+ double secondTerm = pow(absTemp, r) * temp;
+
+ double sigNum = 1.0;
+ if ((axes[j][eDists[i].seq1] - axes[j][eDists[i].seq2]) == 0) { sigNum = 0.0; }
+ else if ((axes[j][eDists[i].seq1] - axes[j][eDists[i].seq2]) < 0) { sigNum = -1.0; }
+
+ double results = (firstTerm * secondTerm * sigNum);
+ cout << i << '\t' << j << '\t' << "results = " << results << endl;
+ gradient[j][eDists[i].seq1] += results;
+ gradient[j][eDists[i].seq2] -= results;
+ }
+ }
+
+ return gradient;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "NMDSCommand", "calculateStressGradientVector");
+ exit(1);
+ }
+ }
+ //**********************************************************************************************************************
+ double NMDSCommand::calculateMagnitude(vector< vector<double> >& gradient) {
+ try {
+ double magnitude = 0.0;
+
+ double sum = 0.0;
+ for (int i = 0; i < gradient.size(); i++) {
+ for (int j = 0; j < gradient[i].size(); j++) {
+ sum += (gradient[i][j] * gradient[i][j]);
+ }
+ }
+
+ magnitude = sqrt(((1.0/(float)gradient[0].size()) * sum));
+
+ return magnitude;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "NMDSCommand", "calculateMagnitude");
+ exit(1);
+ }
+ }
+ //**********************************************************************************************************************
+ //described in Kruskal paper page 121 + 122
+ double NMDSCommand::calculateStep(vector< vector<double> >& prevGrad, vector< vector<double> >& grad, vector<double>& prevStress) {
+ try {
+ double newStep = step;
+
+ //calc the cos theta
+ double sumNum = 0.0;
+ double sumDenom1 = 0.0;
+ double sumDenom2 = 0.0;
+ for (int i = 0; i < prevGrad.size(); i++) {
+ for (int j = 0; j < prevGrad[i].size(); j++) {
+ sumDenom1 += (grad[i][j] * grad[i][j]);
+ sumDenom2 += (prevGrad[i][j] * prevGrad[i][j]);
+ sumNum += (grad[i][j] * prevGrad[i][j]);
+ }
+ }
+
+ double cosTheta = sumNum / (sqrt(sumDenom1) * sqrt(sumDenom2));
+ cosTheta *= cosTheta;
+
+ //calc angle factor
+ double angle = pow(4.0, cosTheta);
+
+ //calc 5 step ratio
+ double currentStress = prevStress[prevStress.size()-1];
+ double lastStress = prevStress[0];
+ if (prevStress.size() > 1) { lastStress = prevStress[prevStress.size()-2]; }
+ double fivePrevStress = prevStress[0];
+ if (prevStress.size() > 5) { fivePrevStress = prevStress[prevStress.size()-6]; }
+
+ double fiveStepRatio = min(1.0, (currentStress / fivePrevStress));
+
+ //calc relaxation factor
+ double relaxation = 1.3 / (1.0 + pow(fiveStepRatio, 5.0));
+
+ //calc good luck factor
+ double goodLuck = min(1.0, (currentStress / lastStress));
+
+ //calc newStep
+ //cout << "\ncos = " << cosTheta << " step = " << step << " angle = " << angle << " relaxation = " << relaxation << " goodluck = " << goodLuck << endl;
+ newStep = step * angle * relaxation * goodLuck;
+
+ return newStep;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "NMDSCommand", "calculateStep");
+ exit(1);
+ }
+ }
+ //**********************************************************************************************************************
+ vector< vector<double> > NMDSCommand::calculateNewConfiguration(double magnitude, vector< vector<double> >& axes, vector< vector<double> >& gradient) {
+ try {
+
+ vector< vector<double> > newAxes = axes;
+
+ for (int i = 0; i < newAxes.size(); i++) {
+
+ if (m->control_pressed) { return newAxes; }
+
+ for (int j = 0; j < newAxes[i].size(); j++) {
+ newAxes[i][j] = axes[i][j] + ((step / magnitude) * gradient[i][j]);
+ }
+ }
+
+ return newAxes;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "NMDSCommand", "calculateNewConfiguration");
+ exit(1);
+ }
+ }*/
+/**********************************************************************************************************************
+ //adjust eDists so that it creates monotonically increasing series of succesive values that increase or stay the same, but never decrease
+ vector<seqDist> NMDSCommand::satisfyMonotonicity(vector<seqDist> eDists, vector<int> partitions) {
+ try {
+
+ //find averages of each partitions
+ vector<double> sums; sums.resize(partitions.size(), 0.0);
+ vector<int> sizes; sizes.resize(partitions.size(), 0);
+
+ for (int i = 0; i < partitions.size(); i++) {
+ //i is not the last one
+ int start = partitions[i];
+ int end;
+ if (i != (partitions.size()-1)) { end = partitions[i+1]; }
+ else{ end = eDists.size(); }
+
+ for (int j = start; j < end; j++) { sums[i] += eDists[j].dist; }
+
+ sizes[i] = (end - start);
+ }
+
+
+ vector<seqDist> D = eDists;
+
+ //i represents the "active block"
+ int i = 0;
+ while (i < partitions.size()) {
+
+ if (m->control_pressed) { return D; }
+
+ bool upActive = true;
+ bool upSatisfied = false;
+ bool downSatisfied = false;
+
+ //while we are not done with this block
+ while ((!upSatisfied) || (!downSatisfied)) {
+
+ if (upActive) {
+
+ //are we are upSatisfied? - is the average of the next block greater than mine?
+ if (i != (partitions.size()-1)) { //if we are the last guy then we are upsatisfied
+ if ((sums[i+1]/(float)sizes[i+1]) >= (sums[i]/(float)sizes[i])) {
+ upSatisfied = true;
+ upActive = false;
+ }else {
+ //find new weighted average
+ double newSum = sums[i] + sums[i+1];
+
+ //merge blocks - putting everything in i
+ sums[i] = newSum;
+ sizes[i] += sizes[i+1];
+ partitions[i] = partitions[i+1];
+
+ sums.erase(sums.begin()+(i+1));
+ sizes.erase(sizes.begin()+(i+1));
+ partitions.erase(partitions.begin()+(i+1));
+
+ upActive = false;
+ }
+ }else { upSatisfied = true; upActive = false; }
+
+ }else { //downActive
+
+ //are we are DownSatisfied? - is the average of the previous block less than mine?
+ if (i != 0) { //if we are the first guy then we are downSatisfied
+ if ((sums[i-1]/(float)sizes[i-1]) <= (sums[i]/(float)sizes[i])) {
+ downSatisfied = true;
+ upActive = true;
+ }else {
+ //find new weighted average
+ double newSum = sums[i] + sums[i-1];;
+
+ //merge blocks - putting everything in i-1
+ sums[i-1] = newSum;
+ sizes[i-1] += sizes[i];
+
+ sums.erase(sums.begin()+i);
+ sizes.erase(sizes.begin()+i);
+ partitions.erase(partitions.begin()+i);
+ i--;
+
+ upActive = true;
+ }
+ }else { downSatisfied = true; upActive = true; }
+ }
+ }
+
+ i++; // go to next block
+ }
+
+ //sanity check - for rounding errors
+ vector<double> averages; averages.resize(sums.size(), 0.0);
+ for (int i = 0; i < sums.size(); i++) { averages[i] = sums[i] / (float) sizes[i]; }
+ for (int i = 0; i < averages.size(); i++) { if (averages[i+1] < averages[i]) { averages[i+1] = averages[i]; } }
+
+ //fill D
+ int placeHolder = 0;
+ for (int i = 0; i < averages.size(); i++) {
+ for (int j = 0; j < sizes[i]; j++) {
+ D[placeHolder].dist = averages[i];
+ placeHolder++;
+ }
+ }
+
+ return D;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "NMDSCommand", "satisfyMonotonicity");
+ exit(1);
+ }
+ }*/
+