X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=treegroupscommand.cpp;h=4b4f21c86366000b8efb195bf94232059a78d03c;hb=aca78ed4a47dff8672ea8fd93cef0dfbaf0f7495;hp=90d336da38f57f70fa0f5ddb59125cc2099ec346;hpb=1a5c2356c1b955c6ec024b2baf9f46377ee7c72e;p=mothur.git diff --git a/treegroupscommand.cpp b/treegroupscommand.cpp index 90d336d..4b4f21c 100644 --- a/treegroupscommand.cpp +++ b/treegroupscommand.cpp @@ -25,7 +25,7 @@ vector TreeGroupCommand::setParameters(){ CommandParameter pprecision("precision", "Number", "", "100", "", "", "","",false,false); parameters.push_back(pprecision); CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel); CommandParameter pgroups("groups", "String", "", "", "", "", "","",false,false); parameters.push_back(pgroups); - CommandParameter pcalc("calc", "Multiple", "sharedsobs-sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan-kstest-sharednseqs-ochiai-anderberg-kulczynski-kulczynskicody-lennon-morisitahorn-braycurtis-whittaker-odum-canberra-structeuclidean-structchord-hellinger-manhattan-structpearson-soergel-spearman-structkulczynski-speciesprofile-hamming-structchi2-gower-memchi2-memchord-memeuclidean-mempearson", "jclass-thetayc", "", "", "","",true,false,true); parameters.push_back(pcalc); + CommandParameter pcalc("calc", "Multiple", "sharedsobs-sharedchao-sharedace-jabund-sorabund-jclass-sorclass-jest-sorest-thetayc-thetan-kstest-sharednseqs-ochiai-anderberg-kulczynski-kulczynskicody-lennon-morisitahorn-braycurtis-whittaker-odum-canberra-structeuclidean-structchord-hellinger-manhattan-structpearson-soergel-spearman-structkulczynski-speciesprofile-hamming-structchi2-gower-memchi2-memchord-memeuclidean-mempearson-jsd-rjsd", "jclass-thetayc", "", "", "","",true,false,true); parameters.push_back(pcalc); CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors); //CommandParameter poutput("output", "Multiple", "lt-square", "lt", "", "", "",false,false); parameters.push_back(poutput); @@ -397,7 +397,12 @@ int TreeGroupCommand::execute(){ treeCalculators.push_back(new MemEuclidean()); }else if (Estimators[i] == "mempearson") { treeCalculators.push_back(new MemPearson()); - } + }else if (Estimators[i] == "jsd") { + treeCalculators.push_back(new JSD()); + }else if (Estimators[i] == "rjsd") { + treeCalculators.push_back(new RJSD()); + } + } } @@ -429,7 +434,10 @@ int TreeGroupCommand::execute(){ m->Treenames.clear(); //fills globaldatas tree names - m->Treenames = m->getGroups(); + //m->Treenames = m->getGroups(); + for (int k = 0; k < lookup.size(); k++) { + m->Treenames.push_back(lookup[k]->getGroup()); + } if (m->control_pressed) { return 0; } @@ -455,7 +463,7 @@ int TreeGroupCommand::execute(){ readMatrix->read(nameMap); }else if (countfile != "") { ct = new CountTable(); - ct->readTable(countfile); + ct->readTable(countfile, true, false); readMatrix->read(ct); }else { readMatrix->read(nameMap); @@ -463,7 +471,10 @@ int TreeGroupCommand::execute(){ list = readMatrix->getListVector(); SparseDistanceMatrix* dMatrix = readMatrix->getDMatrix(); - + + //clear globaldatas old tree names if any + m->Treenames.clear(); + //make treemap if (ct != NULL) { delete ct; } ct = new CountTable(); @@ -475,17 +486,12 @@ int TreeGroupCommand::execute(){ nameMap.insert(bin); gps.insert(bin); groupMap[bin] = bin; + m->Treenames.push_back(bin); } ct->createTable(nameMap, groupMap, gps); vector namesGroups = ct->getNamesOfGroups(); m->setGroups(namesGroups); - - //clear globaldatas old tree names if any - m->Treenames.clear(); - - //fills globaldatas tree names - m->Treenames = m->getGroups(); //used in tree constructor m->runParse = false; @@ -643,6 +649,7 @@ int TreeGroupCommand::makeSimsShared() { }else { m->clearGroups(); Groups.clear(); + m->Treenames.clear(); vector temp; for (int i = 0; i < lookup.size(); i++) { if (lookup[i]->getNumSeqs() < subsampleSize) { @@ -651,6 +658,7 @@ int TreeGroupCommand::makeSimsShared() { }else { Groups.push_back(lookup[i]->getGroup()); temp.push_back(lookup[i]); + m->Treenames.push_back(lookup[i]->getGroup()); } } lookup = temp; @@ -659,8 +667,18 @@ int TreeGroupCommand::makeSimsShared() { if (lookup.size() < 2) { m->mothurOut("You have not provided enough valid groups. I cannot run the command."); m->mothurOutEndLine(); m->control_pressed = true; return 0; } } - numGroups = lookup.size(); + + //sanity check to make sure processors < numComparisions + int numDists = 0; + for(int i=0;i processors) { break; } + } + } + if (numDists < processors) { processors = numDists; } + lines.resize(processors); for (int i = 0; i < processors; i++) { lines[i].start = int (sqrt(float(i)/float(processors)) * numGroups); @@ -782,7 +800,7 @@ int TreeGroupCommand::process(vector thisLookup) { #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix) //loop through and create all the processes you want while (process != processors) { - int pid = fork(); + pid_t pid = fork(); if (pid > 0) { processIDS.push_back(pid); @@ -791,7 +809,7 @@ int TreeGroupCommand::process(vector thisLookup) { driver(thisItersLookup, lines[process].start, lines[process].end, calcDists); - string tempdistFileName = m->getRootName(m->getSimpleName(sharedfile)) + toString(getpid()) + ".dist"; + string tempdistFileName = m->getRootName(m->getSimpleName(sharedfile)) + m->mothurGetpid(process) + ".dist"; ofstream outtemp; m->openOutputFile(tempdistFileName, outtemp); @@ -889,6 +907,9 @@ int TreeGroupCommand::process(vector thisLookup) { //Close all thread handles and free memory allocations. for(int i=0; i < pDataArray.size(); i++){ + if (pDataArray[i]->count != (pDataArray[i]->end-pDataArray[i]->start)) { + m->mothurOut("[ERROR]: process " + toString(i) + " only processed " + toString(pDataArray[i]->count) + " of " + toString(pDataArray[i]->end-pDataArray[i]->start) + " groups assigned to it, quitting. \n"); m->control_pressed = true; + } for (int j = 0; j < pDataArray[i]->thisLookup.size(); j++) { delete pDataArray[i]->thisLookup[j]; } for (int k = 0; k < calcDists.size(); k++) { @@ -912,35 +933,17 @@ int TreeGroupCommand::process(vector thisLookup) { thisItersLookup.clear(); for (int i = 0; i < calcDists.size(); i++) { calcDists[i].clear(); } } + + if (m->debug) { m->mothurOut("[DEBUG]: iter = " + toString(thisIter) + ".\n"); } } - + + if (m->debug) { m->mothurOut("[DEBUG]: done with iters.\n"); } + if (iters != 1) { //we need to find the average distance and standard deviation for each groups distance + vector< vector > calcAverages = m->getAverages(calcDistsTotals); - vector< vector > calcAverages; calcAverages.resize(treeCalculators.size()); - for (int i = 0; i < calcAverages.size(); i++) { //initialize sums to zero. - calcAverages[i].resize(calcDistsTotals[0][i].size()); - - for (int j = 0; j < calcAverages[i].size(); j++) { - calcAverages[i][j].seq1 = calcDists[i][j].seq1; - calcAverages[i][j].seq2 = calcDists[i][j].seq2; - calcAverages[i][j].dist = 0.0; - } - } - - for (int thisIter = 0; thisIter < iters; thisIter++) { //sum all groups dists for each calculator - for (int i = 0; i < calcAverages.size(); i++) { //initialize sums to zero. - for (int j = 0; j < calcAverages[i].size(); j++) { - calcAverages[i][j].dist += calcDistsTotals[thisIter][i][j].dist; - } - } - } - - for (int i = 0; i < calcAverages.size(); i++) { //finds average. - for (int j = 0; j < calcAverages[i].size(); j++) { - calcAverages[i][j].dist /= (float) iters; - } - } + if (m->debug) { m->mothurOut("[DEBUG]: found averages.\n"); } //create average tree for each calc for (int i = 0; i < calcDists.size(); i++) { @@ -971,6 +974,8 @@ int TreeGroupCommand::process(vector thisLookup) { if (newTree != NULL) { writeTree(outputFile, newTree); } } + if (m->debug) { m->mothurOut("[DEBUG]: done averages trees.\n"); } + //create all trees for each calc and find their consensus tree for (int i = 0; i < calcDists.size(); i++) { if (m->control_pressed) { break; } @@ -1002,7 +1007,7 @@ int TreeGroupCommand::process(vector thisLookup) { int row = calcDistsTotals[myIter][i][j].seq1; int column = calcDistsTotals[myIter][i][j].seq2; double dist = calcDistsTotals[myIter][i][j].dist; - + matrix[row][column] = dist; matrix[column][row] = dist; } @@ -1017,11 +1022,15 @@ int TreeGroupCommand::process(vector thisLookup) { outAll.close(); if (m->control_pressed) { for (int k = 0; k < trees.size(); k++) { delete trees[k]; } } + if (m->debug) { m->mothurOut("[DEBUG]: done all trees.\n"); } + Consensus consensus; //clear old tree names if any m->Treenames.clear(); m->Treenames = m->getGroups(); //may have changed if subsample eliminated groups Tree* conTree = consensus.getTree(trees); + if (m->debug) { m->mothurOut("[DEBUG]: done cons tree.\n"); } + //create a new filename variables["[tag]"] = "cons"; string conFile = getOutputFileName("tree",variables);