working on windows paralellization, added trimOligos class to be used by trim.flows...

[mothur.git] / shhhercommand.h

diff --git a/shhhercommand.h b/shhhercommand.h
index bdc850362b4b08e6c2e14655360e455353e74793..00bd41abef486338da1d33463b0698bd895266a1 100644 (file)
--- a/shhhercommand.h
+++ b/shhhercommand.h
@@ -13,6 +13,15 @@
 #include "mothur.h"
 #include "command.hpp"
 
+//**********************************************************************************************************************
+
+#define NUMBINS 1000
+#define HOMOPS 10
+#define MIN_COUNT 0.1
+#define MIN_WEIGHT 0.1
+#define MIN_TAU 0.0001
+#define MIN_ITER 10
+//**********************************************************************************************************************
 
 class ShhherCommand : public Command {
        
@@ -116,6 +125,125 @@ private:
        
 };
 
+/**************************************************************************************************/
+//custom data structure for threads to use.
+// This is passed by void pointer so it can be any data type
+// that can be passed using a single void pointer (LPVOID).
+typedef struct flowDistParentForkData {
+       string distFileName; 
+       vector<int> mapUniqueToSeq;
+       vector<int> mapSeqToUnique;
+       vector<int> lengths;
+       vector<short> flowDataIntI;
+       vector<double> flowDataPrI;
+       vector<double> jointLookUp;
+       MothurOut* m;
+       int threadID, startSeq, stopSeq, numFlowCells;
+       float cutoff;
+       
+       flowDistParentForkData(){}
+       flowDistParentForkData(string d, vector<int> mapU, vector<int> mapS, vector<int> l, vector<short> flowD, vector<double> flowDa, vector<double> j, MothurOut* mout, int st, int sp, int n, float cut, int tid) {
+               distFileName = d;
+               mapUniqueToSeq = mapU;
+               mapSeqToUnique = mapS;
+               lengths = l;
+               flowDataIntI = flowD;
+               flowDataPrI = flowDa;
+               jointLookUp = j;
+               m = mout;
+               startSeq = st;
+               stopSeq = sp;
+               numFlowCells = n;
+               cutoff= cut;
+               threadID = tid;
+       }
+};
+
+/**************************************************************************************************/
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+#else
+static DWORD WINAPI MyflowDistParentForkThreadFunction(LPVOID lpParam){ 
+       flowDistParentForkData* pDataArray;
+       pDataArray = (flowDistParentForkData*)lpParam;
+       
+       try {
+               ostringstream outStream;
+               outStream.setf(ios::fixed, ios::floatfield);
+               outStream.setf(ios::dec, ios::basefield);
+               outStream.setf(ios::showpoint);
+               outStream.precision(6);
+               
+               int begTime = time(NULL);
+               double begClock = clock();
+               string tempOut = "start and end = " + toString(pDataArray->startSeq) +'\t' + toString(pDataArray->stopSeq) + "-";
+               cout << tempOut << endl;
+       
+               for(int i=pDataArray->startSeq;i<pDataArray->stopSeq;i++){
+                       
+                       if (pDataArray->m->control_pressed) { break; }
+                       cout << "thread i = " << i << endl;
+                       for(int j=0;j<i;j++){
+                               
+                               cout << "thread j = " << j << endl;
+                               float flowDistance = 0.0;
+                               ////////////////// calcPairwiseDist ///////////////////
+                               //needed because this is a static global function that can't see the classes internal functions
+                               
+                               int minLength = pDataArray->lengths[pDataArray->mapSeqToUnique[pDataArray->mapUniqueToSeq[i]]];
+                               if(pDataArray->lengths[pDataArray->mapUniqueToSeq[j]] < minLength){     minLength = pDataArray->lengths[pDataArray->mapSeqToUnique[pDataArray->mapUniqueToSeq[j]]];     }
+                               
+                               int ANumFlowCells = pDataArray->mapUniqueToSeq[i] * pDataArray->numFlowCells;
+                               int BNumFlowCells = pDataArray->mapUniqueToSeq[j] * pDataArray->numFlowCells;
+                               
+                               for(int k=0;k<minLength;k++){
+                                       
+                                       if (pDataArray->m->control_pressed) { break; }
+                                       
+                                       int flowAIntI = pDataArray->flowDataIntI[ANumFlowCells + k];
+                                       float flowAPrI = pDataArray->flowDataPrI[ANumFlowCells + k];
+                                       
+                                       int flowBIntI = pDataArray->flowDataIntI[BNumFlowCells + k];
+                                       float flowBPrI = pDataArray->flowDataPrI[BNumFlowCells + k];
+                                       flowDistance += pDataArray->jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
+                               }
+                               
+                               flowDistance /= (float) minLength;
+                               //cout << flowDistance << endl;
+                               ////////////////// end of calcPairwiseDist ///////////////////
+                                                               
+                               if(flowDistance < 1e-6){
+                                       outStream << pDataArray->mapUniqueToSeq[i] << '\t' << pDataArray->mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
+                               }
+                               else if(flowDistance <= pDataArray->cutoff){
+                                       outStream << pDataArray->mapUniqueToSeq[i] << '\t' << pDataArray->mapUniqueToSeq[j] << '\t' << flowDistance << endl;
+                               }
+                       }
+                       if(i % 100 == 0){
+                               pDataArray->m->mothurOut(toString(i) + "\t" + toString(time(NULL) - begTime));
+                               pDataArray->m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
+                               pDataArray->m->mothurOutEndLine();
+                       }
+               }
+               
+               ofstream distFile(pDataArray->distFileName.c_str());
+               distFile << outStream.str();            
+               distFile.close();
+               
+               if (pDataArray->m->control_pressed) {}
+               else {
+                       pDataArray->m->mothurOut(toString(pDataArray->stopSeq-1) + "\t" + toString(time(NULL) - begTime));
+                       pDataArray->m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
+                       pDataArray->m->mothurOutEndLine();
+               }               
+               
+       }
+       catch(exception& e) {
+               pDataArray->m->errorOut(e, "ShhherCommand", "MyflowDistParentForkThreadFunction");
+               exit(1);
+       }
+} 
+#endif
+
 
 #endif