+ it = trim.begin();
+
+ string refQuery = query->getAligned();
+
+ //j<i, so you don't find distances from i to j and then j to i.
+ for (int j = 0; j < ref.size(); j++) {
+
+ string refJ = ref[j].seq->getAligned();
+
+ for (int k = 0; k < windows.size(); k++) {
+
+ string QueryWindowk, refJWindowk;
+
+ if (k < windows.size()-1) {
+ //get window strings
+ QueryWindowk = refQuery.substr(windows[k], windowSizes);
+ refJWindowk = refJ.substr(windows[k], windowSizes);
+ }else { //last window may be smaller than rest - see findwindows
+ //get window strings
+ QueryWindowk = refQuery.substr(windows[k], (it->second-windows[k]));
+ refJWindowk = refJ.substr(windows[k], (it->second-windows[k]));
+ }
+
+ //find differences
+ int diff = getDiff(QueryWindowk, refJWindowk);
+
+ //save differences
+ diffs[j][k] = diff;
+
+ }//k
+ }//j
+
+
+ //initialize sumRef for this query
+ sumQuery.resize(windows.size(), 0);
+ sumSquaredQuery.resize(windows.size(), 0);
+ averageQuery.resize(windows.size(), 0);
+
+ //find the sum of the differences
+ for (int j = 0; j < diffs.size(); j++) {
+ for (int k = 0; k < diffs[j].size(); k++) {
+ sumQuery[k] += diffs[j][k];
+ sumSquaredQuery[k] += (diffs[j][k]*diffs[j][k]);
+ }//k
+ }//j
+
+
+ //find the average of the differences for the references for each window
+ for (int i = 0; i < windows.size(); i++) {
+ averageQuery[i] = sumQuery[i] / (float) ref.size();
+ }
+ }
+ catch(exception& e) {
+ m->errorOut(e, "Ccode", "getAverageQuery");
+ exit(1);
+ }
+}
+/**************************************************************************************************/
+void Ccode::findVarianceRef() {
+ try {
+
+ varRef.resize(windows.size(), 0);
+ sdRef.resize(windows.size(), 0);
+
+ //for each window
+ for (int i = 0; i < windows.size(); i++) {
+ varRef[i] = (sumSquaredRef[i] - ((sumRef[i]*sumRef[i])/(float)refCombo)) / (float)(refCombo-1);
+ sdRef[i] = sqrt(varRef[i]);
+
+ //set minimum error rate to 0.001 - to avoid potential divide by zero - not sure if this is necessary but it follows ccode implementation
+ if (averageRef[i] < 0.001) { averageRef[i] = 0.001; }
+ if (sumRef[i] < 0.001) { sumRef[i] = 0.001; }
+ if (varRef[i] < 0.001) { varRef[i] = 0.001; }
+ if (sumSquaredRef[i] < 0.001) { sumSquaredRef[i] = 0.001; }
+ if (sdRef[i] < 0.001) { sdRef[i] = 0.001; }