added amova command

[mothur.git] / shhhercommand.cpp

diff --git a/shhhercommand.cpp b/shhhercommand.cpp
index 24301a2c18b3decebd1a9c5a3954ffc96d57b90b..ffddb8cfc9de4dd92feb4966dd8d31477e1363fa 100644 (file)
--- a/shhhercommand.cpp
+++ b/shhhercommand.cpp
@@ -16,6 +16,7 @@
 #include "listvector.hpp"
 #include "cluster.hpp"
 #include "sparsematrix.hpp"
+#include <cfloat>
 
 //**********************************************************************************************************************
 
@@ -183,7 +184,7 @@ ShhherCommand::ShhherCommand(string option) {
                        // ...at some point should added some additional type checking...
                        string temp;
                        temp = validParameter.validFile(parameters, "lookup", true);
-                       if (temp == "not found")        {       lookupFileName = "LookUp_E123.pat";     }
+                       if (temp == "not found")        {       lookupFileName = "LookUp_Titanium.pat"; }
                        else if(temp == "not open")     {       abort = true;                   } 
                        else                                            {       lookupFileName = temp;  }
                        
@@ -259,7 +260,7 @@ int ShhherCommand::execute(){
                        cout << "\nGetting preliminary data..." << endl;
                        getSingleLookUp();
                        getJointLookUp();
-                                               
+                       
                        vector<string> flowFileVector;
                        if(flowFilesFileName != "not found"){
                                string fName;
@@ -655,8 +656,7 @@ int ShhherCommand::execute(){
                
                getSingleLookUp();
                getJointLookUp();
-               
-               
+                               
                vector<string> flowFileVector;
                if(flowFilesFileName != "not found"){
                        string fName;
@@ -838,10 +838,10 @@ void ShhherCommand::getJointLookUp(){
                for(int i=0;i<NUMBINS;i++){
                        for(int j=0;j<NUMBINS;j++){             
                                
-                               float minSum = 10000000000;
+                               double minSum = 100000000;
                                
                                for(int k=0;k<HOMOPS;k++){
-                                       float sum = singleLookUp[k * NUMBINS + i] + singleLookUp[k * NUMBINS + j];
+                                       double sum = singleLookUp[k * NUMBINS + i] + singleLookUp[k * NUMBINS + j];
                                        
                                        if(sum < minSum)        {       minSum = sum;           }
                                }       
@@ -857,9 +857,11 @@ void ShhherCommand::getJointLookUp(){
 
 /**************************************************************************************************/
 
-float ShhherCommand::getProbIntensity(int intIntensity){                          
+double ShhherCommand::getProbIntensity(int intIntensity){                          
        try{
-               float minNegLogProb = 10000000000; 
+
+               double minNegLogProb = 100000000; 
+
                
                for(int i=0;i<HOMOPS;i++){//loop signal strength
                        float negLogProb = singleLookUp[i * NUMBINS + intIntensity];
@@ -1006,7 +1008,7 @@ void ShhherCommand::flowDistParentFork(string distFileName, int startSeq, int st
                m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
                m->mothurOutEndLine();
                
-               ofstream distFile((distFileName + "temp." + toString(pid)).c_str());
+               ofstream distFile(distFileName.c_str());
                distFile << outStream.str();            
                distFile.close();
        }
@@ -1160,7 +1162,6 @@ string ShhherCommand::cluster(string distFileName, string namesFileName){
                double clusterCutoff = cutoff;
                while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
                        cluster->update(clusterCutoff);
-                       float dist = matrix->getSmallDist();
                }
                
                list->setLabel(toString(cutoff));
@@ -1411,7 +1412,6 @@ void ShhherCommand::calcCentroidsDriver(int start, int finish){
                                for(int j=0;j<nSeqsPerOTU[i];j++){
                                        int index = cumNumSeqs[i] + j;
                                        int nI = seqIndex[index];
-                                       int nIU = mapSeqToUnique[nI];
                                        
                                        double tauValue = singleTau[seqNumber[index]];
                                        
@@ -1483,7 +1483,6 @@ double ShhherCommand::getNewWeights(){
                        
                        for(int j=0;j<nSeqsPerOTU[i];j++){
                                int index = cumNumSeqs[i] + j;
-                               int nI = seqIndex[index];
                                double tauValue = singleTau[seqNumber[index]];
                                weight[i] += tauValue;
                        }
@@ -1505,7 +1504,7 @@ double ShhherCommand::getLikelihood(){
        
        try{
                
-               vector<double> P(numSeqs, 0);
+               vector<long double> P(numSeqs, 0);
                int effNumOTUs = 0;
                
                for(int i=0;i<numOTUs;i++){
@@ -1514,6 +1513,7 @@ double ShhherCommand::getLikelihood(){
                        }
                }
                
+               string hold;
                for(int i=0;i<numOTUs;i++){
                        for(int j=0;j<nSeqsPerOTU[i];j++){
                                int index = cumNumSeqs[i] + j;
@@ -1523,14 +1523,15 @@ double ShhherCommand::getLikelihood(){
                                P[nI] += weight[i] * exp(-singleDist * sigma);
                        }
                }
-               
                double nLL = 0.00;
                for(int i=0;i<numSeqs;i++){
+                       if(P[i] == 0){  P[i] = DBL_EPSILON;     }
+
                        nLL += -log(P[i]);
                }
                
                nLL = nLL -(double)numSeqs * log(sigma);
-               
+
                return nLL; 
        }
        catch(exception& e) {
@@ -1900,7 +1901,6 @@ void ShhherCommand::writeQualities(vector<int> otuCounts){
                vector<vector<int> > qualities(numOTUs);
                vector<double> pr(HOMOPS, 0);
                
-               int index = 0;
                
                for(int i=0;i<numOTUs;i++){
                        int index = 0;