+ string name;
+ //Name Length Overlap_Length Overlap_Start Overlap_End MisMatches Num_Ns
+ int length, OLength, thisOStart, thisOEnd, numMisMatches, numns;
+
+ ifstream in;
+ m->openInputFile(contigsreport, in);
+
+ in.seekg(filePos.start);
+ if (filePos.start == 0) { //read headers
+ m->getline(in); m->gobble(in);
+ }
+
+ bool done = false;
+ int count = 0;
+
+ while (!done) {
+
+ if (m->control_pressed) { in.close(); return 1; }
+
+ //seqname start end nbases ambigs polymer numSeqs
+ in >> name >> length >> OLength >> thisOStart >> thisOEnd >> numMisMatches >> numns; m->gobble(in);
+
+ int num = 1;
+ if ((namefile != "") || (countfile !="")){
+ //make sure this sequence is in the namefile, else error
+ map<string, int>::iterator it = nameMap.find(name);
+
+ if (it == nameMap.end()) { m->mothurOut("[ERROR]: " + name + " is not in your namefile, please correct."); m->mothurOutEndLine(); m->control_pressed = true; }
+ else { num = it->second; }
+ }
+
+ //for each sequence this sequence represents
+ for (int i = 0; i < num; i++) {
+ ostartPosition.push_back(thisOStart);
+ oendPosition.push_back(thisOEnd);
+ oLength.push_back(OLength);
+ omismatches.push_back(numMisMatches);
+ numNs.push_back(numns);
+ }
+
+ count++;
+
+ //if((count) % 100 == 0){ m->mothurOut("Optimizing sequence: " + toString(count)); m->mothurOutEndLine(); }
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+ unsigned long long pos = in.tellg();
+ if ((pos == -1) || (pos >= filePos.end)) { break; }
+#else
+ if (in.eof()) { break; }
+#endif
+ }
+
+ in.close();
+
+ return count;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "driverContigsSummary");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+int ScreenSeqsCommand::createProcessesContigsSummary(vector<int>& oLength, vector<int>& ostartPosition, vector<int>& oendPosition, vector<int>& omismatches, vector<int>& numNs, vector<linePair> contigsLines) {
+ try {
+
+ int process = 1;
+ int num = 0;
+ vector<int> processIDS;
+
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+
+ //loop through and create all the processes you want
+ while (process != processors) {
+ int pid = fork();
+
+ if (pid > 0) {
+ processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
+ process++;
+ }else if (pid == 0){
+ num = driverContigsSummary(oLength, ostartPosition, oendPosition, omismatches, numNs, contigsLines[process]);
+
+ //pass numSeqs to parent
+ ofstream out;
+ string tempFile = contigsreport + toString(getpid()) + ".num.temp";
+ m->openOutputFile(tempFile, out);
+
+ out << num << endl;
+ out << ostartPosition.size() << endl;
+ for (int k = 0; k < ostartPosition.size(); k++) { out << ostartPosition[k] << '\t'; } out << endl;
+ for (int k = 0; k < oendPosition.size(); k++) { out << oendPosition[k] << '\t'; } out << endl;
+ for (int k = 0; k < oLength.size(); k++) { out << oLength[k] << '\t'; } out << endl;
+ for (int k = 0; k < omismatches.size(); k++) { out << omismatches[k] << '\t'; } out << endl;
+ for (int k = 0; k < numNs.size(); k++) { out << numNs[k] << '\t'; } out << endl;
+
+ out.close();
+
+ exit(0);
+ }else {
+ m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
+ for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
+ exit(0);
+ }
+ }
+
+ num = driverContigsSummary(oLength, ostartPosition, oendPosition, omismatches, numNs, contigsLines[0]);
+
+ //force parent to wait until all the processes are done
+ for (int i=0;i<processIDS.size();i++) {
+ int temp = processIDS[i];
+ wait(&temp);
+ }
+
+ //parent reads in and combine Filter info
+ for (int i = 0; i < processIDS.size(); i++) {
+ string tempFilename = contigsreport + toString(processIDS[i]) + ".num.temp";
+ ifstream in;
+ m->openInputFile(tempFilename, in);
+
+ int temp, tempNum;
+ in >> tempNum; m->gobble(in); num += tempNum;
+ in >> tempNum; m->gobble(in);
+ for (int k = 0; k < tempNum; k++) { in >> temp; ostartPosition.push_back(temp); } m->gobble(in);
+ for (int k = 0; k < tempNum; k++) { in >> temp; oendPosition.push_back(temp); } m->gobble(in);
+ for (int k = 0; k < tempNum; k++) { in >> temp; oLength.push_back(temp); } m->gobble(in);
+ for (int k = 0; k < tempNum; k++) { in >> temp; omismatches.push_back(temp); } m->gobble(in);
+ for (int k = 0; k < tempNum; k++) { in >> temp; numNs.push_back(temp); } m->gobble(in);
+
+ in.close();
+ m->mothurRemove(tempFilename);
+ }
+
+
+#else
+ //////////////////////////////////////////////////////////////////////////////////////////////////////
+ //Windows version shared memory, so be careful when passing variables through the seqSumData struct.
+ //Above fork() will clone, so memory is separate, but that's not the case with windows,
+ //Taking advantage of shared memory to allow both threads to add info to vectors.
+ //////////////////////////////////////////////////////////////////////////////////////////////////////
+ /*
+ vector<contigsSumData*> pDataArray;
+ DWORD dwThreadIdArray[processors-1];
+ HANDLE hThreadArray[processors-1];
+
+ //Create processor worker threads.
+ for( int i=0; i<processors-1; i++ ){
+
+ // Allocate memory for thread data.
+ contigsSumData* tempSum = new contigsSumData(contigsreport, m, contigsLines[i].start, contigsLines[i].end, namefile, countfile, nameMap);
+ pDataArray.push_back(tempSum);
+
+ //MySeqSumThreadFunction is in header. It must be global or static to work with the threads.
+ //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
+ hThreadArray[i] = CreateThread(NULL, 0, MyContigsSumThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
+ }
+ */
+ contigsLines[processors-1].start = 0;
+ //do your part
+ num = driverContigsSummary(oLength, ostartPosition, oendPosition, omismatches, numNs, contigsLines[processors-1]);
+ /*
+ //Wait until all threads have terminated.
+ WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
+
+ //Close all thread handles and free memory allocations.
+ for(int i=0; i < pDataArray.size(); i++){
+ num += pDataArray[i]->count;
+ if (pDataArray[i]->count != pDataArray[i]->end) {
+ m->mothurOut("[ERROR]: process " + toString(i) + " only processed " + toString(pDataArray[i]->count) + " of " + toString(pDataArray[i]->end) + " sequences assigned to it, quitting. \n"); m->control_pressed = true;
+ }
+ for (int k = 0; k < pDataArray[i]->ostartPosition.size(); k++) { ostartPosition.push_back(pDataArray[i]->ostartPosition[k]); }
+ for (int k = 0; k < pDataArray[i]->oendPosition.size(); k++) { oendPosition.push_back(pDataArray[i]->oendPosition[k]); }
+ for (int k = 0; k < pDataArray[i]->oLength.size(); k++) { oLength.push_back(pDataArray[i]->oLength[k]); }
+ for (int k = 0; k < pDataArray[i]->omismatches.size(); k++) { omismatches.push_back(pDataArray[i]->omismatches[k]); }
+ for (int k = 0; k < pDataArray[i]->numNs.size(); k++) { numNs.push_back(pDataArray[i]->numNs[k]); }
+ CloseHandle(hThreadArray[i]);
+ delete pDataArray[i];
+ }
+ */
+#endif
+ return num;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "createProcessesContigsSummary");
+ exit(1);
+ }
+}
+//***************************************************************************************************************
+int ScreenSeqsCommand::optimizeAlign(){
+ try {
+
+ vector<float> sims;
+ vector<float> scores;
+ vector<int> inserts;
+
+ vector<unsigned long long> positions;
+ vector<linePair> alignLines;
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+ positions = m->divideFilePerLine(alignreport, processors);
+ for (int i = 0; i < (positions.size()-1); i++) { alignLines.push_back(linePair(positions[i], positions[(i+1)])); }
+#else
+ if(processors == 1){ alignLines.push_back(linePair(0, 1000)); }
+ else {
+ int numAlignSeqs = 0;
+ positions = m->setFilePosEachLine(alignreport, numAlignSeqs);
+ if (positions.size() < processors) { processors = positions.size(); }
+
+ //figure out how many sequences you have to process
+ int numSeqsPerProcessor = numAlignSeqs / processors;
+ for (int i = 0; i < processors; i++) {
+ int startIndex = i * numSeqsPerProcessor;
+ if(i == (processors - 1)){ numSeqsPerProcessor = numAlignSeqs - i * numSeqsPerProcessor; }
+ alignLines.push_back(linePair(positions[startIndex], numSeqsPerProcessor));
+ }
+ }
+#endif
+
+#ifdef USE_MPI
+ int pid;
+ MPI_Comm_rank(MPI_COMM_WORLD, &pid);
+
+ if (pid == 0) {
+ driverAlignSummary(sims, scores, inserts, alignLines[0]);
+#else
+ createProcessesAlignSummary(sims, scores, inserts, alignLines);
+
+ if (m->control_pressed) { return 0; }
+#endif
+ sort(sims.begin(), sims.end());
+ sort(scores.begin(), scores.end());
+ sort(inserts.begin(), inserts.end());
+
+ //numSeqs is the number of unique seqs, startPosition.size() is the total number of seqs, we want to optimize using all seqs
+ int criteriaPercentile = int(sims.size() * (criteria / (float) 100));
+
+ for (int i = 0; i < optimize.size(); i++) {
+ if (optimize[i] == "minsim") { int mincriteriaPercentile = int(sims.size() * ((100 - criteria) / (float) 100)); minSim = sims[mincriteriaPercentile]; m->mothurOut("Optimizing minsim to " + toString(minSim) + "."); m->mothurOutEndLine();}
+ else if (optimize[i] == "minscore") { int mincriteriaPercentile = int(scores.size() * ((100 - criteria) / (float) 100)); minScore = scores[mincriteriaPercentile]; m->mothurOut("Optimizing minscore to " + toString(minScore) + "."); m->mothurOutEndLine(); }
+ else if (optimize[i] == "maxinsert") { maxInsert = inserts[criteriaPercentile]; m->mothurOut("Optimizing maxinsert to " + toString(maxInsert) + "."); m->mothurOutEndLine(); }
+ }
+
+#ifdef USE_MPI
+ }
+
+ MPI_Status status;
+ MPI_Comm_rank(MPI_COMM_WORLD, &pid);
+ MPI_Comm_size(MPI_COMM_WORLD, &processors);
+
+ if (pid == 0) {
+ //send file positions to all processes
+ for(int i = 1; i < processors; i++) {
+ MPI_Send(&minSim, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
+ MPI_Send(&minScore, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
+ MPI_Send(&maxInsert, 1, MPI_INT, i, 2001, MPI_COMM_WORLD);
+ }
+ }else {
+ MPI_Recv(&minSim, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
+ MPI_Recv(&minScore, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
+ MPI_Recv(&maxInsert, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
+ }
+ MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
+#endif
+ return 0;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "optimizeContigs");
+ exit(1);
+ }
+}
+/**************************************************************************************/
+int ScreenSeqsCommand::driverAlignSummary(vector<float>& sims, vector<float>& scores, vector<int>& inserts, linePair filePos) {
+ try {
+
+ string name, TemplateName, SearchMethod, AlignmentMethod;
+ //QueryName QueryLength TemplateName TemplateLength SearchMethod SearchScore AlignmentMethod QueryStart QueryEnd TemplateStart TemplateEnd PairwiseAlignmentLength GapsInQuery GapsInTemplate LongestInsert SimBtwnQuery&Template
+ //checking for minScore, maxInsert, minSim
+ int length, TemplateLength, QueryStart, QueryEnd, TemplateStart, TemplateEnd, PairwiseAlignmentLength, GapsInQuery, GapsInTemplate, LongestInsert;
+ float SearchScore, SimBtwnQueryTemplate;
+
+ ifstream in;
+ m->openInputFile(alignreport, in);
+
+ in.seekg(filePos.start);
+ if (filePos.start == 0) { //read headers
+ m->getline(in); m->gobble(in);
+ }
+
+ bool done = false;
+ int count = 0;
+
+ while (!done) {
+
+ if (m->control_pressed) { in.close(); return 1; }
+
+ in >> name >> length >> TemplateName >> TemplateLength >> SearchMethod >> SearchScore >> AlignmentMethod >> QueryStart >> QueryEnd >> TemplateStart >> TemplateEnd >> PairwiseAlignmentLength >> GapsInQuery >> GapsInTemplate >> LongestInsert >> SimBtwnQueryTemplate; m->gobble(in);
+
+ int num = 1;
+ if ((namefile != "") || (countfile !="")){
+ //make sure this sequence is in the namefile, else error
+ map<string, int>::iterator it = nameMap.find(name);
+
+ if (it == nameMap.end()) { m->mothurOut("[ERROR]: " + name + " is not in your namefile, please correct."); m->mothurOutEndLine(); m->control_pressed = true; }
+ else { num = it->second; }
+ }
+
+ //for each sequence this sequence represents
+ for (int i = 0; i < num; i++) {
+ sims.push_back(SimBtwnQueryTemplate);
+ scores.push_back(SearchScore);
+ inserts.push_back(LongestInsert);
+ }
+
+ count++;
+
+ //if((count) % 100 == 0){ m->mothurOut("Optimizing sequence: " + toString(count)); m->mothurOutEndLine(); }
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+ unsigned long long pos = in.tellg();
+ if ((pos == -1) || (pos >= filePos.end)) { break; }
+#else
+ if (in.eof()) { break; }
+#endif
+ }
+
+ in.close();
+
+ return count;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "driverAlignSummary");
+ exit(1);
+ }
+}
+
+/**************************************************************************************************/
+int ScreenSeqsCommand::createProcessesAlignSummary(vector<float>& sims, vector<float>& scores, vector<int>& inserts, vector<linePair> alignLines) {
+ try {
+
+ int process = 1;
+ int num = 0;
+ vector<int> processIDS;
+
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+
+ //loop through and create all the processes you want
+ while (process != processors) {
+ int pid = fork();
+
+ if (pid > 0) {
+ processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
+ process++;
+ }else if (pid == 0){
+ num = driverAlignSummary(sims, scores, inserts, alignLines[process]);
+
+ //pass numSeqs to parent
+ ofstream out;
+ string tempFile = alignreport + toString(getpid()) + ".num.temp";
+ m->openOutputFile(tempFile, out);
+
+ out << num << endl;
+ out << sims.size() << endl;
+ for (int k = 0; k < sims.size(); k++) { out << sims[k] << '\t'; } out << endl;
+ for (int k = 0; k < scores.size(); k++) { out << scores[k] << '\t'; } out << endl;
+ for (int k = 0; k < inserts.size(); k++) { out << inserts[k] << '\t'; } out << endl;
+
+ out.close();
+
+ exit(0);
+ }else {
+ m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
+ for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
+ exit(0);
+ }
+ }
+
+ num = driverAlignSummary(sims, scores, inserts, alignLines[0]);
+
+ //force parent to wait until all the processes are done
+ for (int i=0;i<processIDS.size();i++) {
+ int temp = processIDS[i];
+ wait(&temp);
+ }
+
+ //parent reads in and combine Filter info
+ for (int i = 0; i < processIDS.size(); i++) {
+ string tempFilename = alignreport + toString(processIDS[i]) + ".num.temp";
+ ifstream in;
+ m->openInputFile(tempFilename, in);
+
+ int temp, tempNum;
+ float temp2;
+ in >> tempNum; m->gobble(in); num += tempNum;
+ in >> tempNum; m->gobble(in);
+ for (int k = 0; k < tempNum; k++) { in >> temp2; sims.push_back(temp2); } m->gobble(in);
+ for (int k = 0; k < tempNum; k++) { in >> temp2; scores.push_back(temp2); } m->gobble(in);
+ for (int k = 0; k < tempNum; k++) { in >> temp; inserts.push_back(temp); } m->gobble(in);
+
+ in.close();
+ m->mothurRemove(tempFilename);
+ }
+
+
+#else
+ //////////////////////////////////////////////////////////////////////////////////////////////////////
+ //Windows version shared memory, so be careful when passing variables through the seqSumData struct.
+ //Above fork() will clone, so memory is separate, but that's not the case with windows,
+ //Taking advantage of shared memory to allow both threads to add info to vectors.
+ //////////////////////////////////////////////////////////////////////////////////////////////////////
+ /*
+ vector<alignsData*> pDataArray;
+ DWORD dwThreadIdArray[processors-1];
+ HANDLE hThreadArray[processors-1];
+
+ //Create processor worker threads.
+ for( int i=0; i<processors-1; i++ ){
+
+ // Allocate memory for thread data.
+ alignsData* tempSum = new alignsData(alignreport, m, alignLines[i].start, alignLines[i].end, namefile, countfile, nameMap);
+ pDataArray.push_back(tempSum);
+
+ //MySeqSumThreadFunction is in header. It must be global or static to work with the threads.
+ //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
+ hThreadArray[i] = CreateThread(NULL, 0, MyAlignsThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
+ }*/
+ alignLines[processors-1].start = 0;
+ //do your part
+ num = driverAlignSummary(sims, scores, inserts, alignLines[processors-1]);
+ /*
+ //Wait until all threads have terminated.
+ WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
+
+ //Close all thread handles and free memory allocations.
+ for(int i=0; i < pDataArray.size(); i++){
+ num += pDataArray[i]->count;
+ if (pDataArray[i]->count != pDataArray[i]->end) {
+ m->mothurOut("[ERROR]: process " + toString(i) + " only processed " + toString(pDataArray[i]->count) + " of " + toString(pDataArray[i]->end) + " sequences assigned to it, quitting. \n"); m->control_pressed = true;
+ }
+ for (int k = 0; k < pDataArray[i]->sims.size(); k++) { sims.push_back(pDataArray[i]->sims[k]); }
+ for (int k = 0; k < pDataArray[i]->scores.size(); k++) { scores.push_back(pDataArray[i]->scores[k]); }
+ for (int k = 0; k < pDataArray[i]->inserts.size(); k++) { inserts.push_back(pDataArray[i]->inserts[k]); }
+ CloseHandle(hThreadArray[i]);
+ delete pDataArray[i];
+ }
+ */
+#endif
+ return num;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ScreenSeqsCommand", "createProcessesAlignSummary");
+ exit(1);
+ }
+}
+//***************************************************************************************************************
+int ScreenSeqsCommand::getSummary(vector<unsigned long long>& positions){
+ try {
+
+ vector<int> startPosition;
+ vector<int> endPosition;
+ vector<int> seqLength;
+ vector<int> ambigBases;
+ vector<int> longHomoPolymer;
+ vector<int> numNs;
+
+ vector<unsigned long long> positions;
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
+ positions = m->divideFile(fastafile, processors);
+ for (int i = 0; i < (positions.size()-1); i++) { lines.push_back(linePair(positions[i], positions[(i+1)])); }
+#else
+ if(processors == 1){ lines.push_back(linePair(0, 1000)); }
+ else {
+ int numFastaSeqs = 0;
+ positions = m->setFilePosFasta(fastafile, numFastaSeqs);
+ if (positions.size() < processors) { processors = positions.size(); }
+
+ //figure out how many sequences you have to process
+ int numSeqsPerProcessor = numFastaSeqs / processors;
+ for (int i = 0; i < processors; i++) {
+ int startIndex = i * numSeqsPerProcessor;
+ if(i == (processors - 1)){ numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor; }
+ lines.push_back(linePair(positions[startIndex], numSeqsPerProcessor));
+ }
+ }
+#endif
+
+#ifdef USE_MPI
+ int pid;
+ MPI_Comm_rank(MPI_COMM_WORLD, &pid);
+
+ if (pid == 0) {
+ driverCreateSummary(startPosition, endPosition, seqLength, ambigBases, longHomoPolymer, numNs, fastafile, lines[0]);
+#else