shhhercommand.cpp

   1 /*
   2  *  shhher.cpp
   3  *  Mothur
   4  *
   5  *  Created by Pat Schloss on 12/27/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "shhhercommand.h"
  11
  12 #include "readcolumn.h"
  13 #include "readmatrix.hpp"
  14 #include "rabundvector.hpp"
  15 #include "sabundvector.hpp"
  16 #include "listvector.hpp"
  17 #include "cluster.hpp"
  18 #include "sparsematrix.hpp"
  19 #include <cfloat>
  20
  21 //**********************************************************************************************************************
  22
  23 #define NUMBINS 1000
  24 #define HOMOPS 10
  25 #define MIN_COUNT 0.1
  26 #define MIN_WEIGHT 0.1
  27 #define MIN_TAU 0.0001
  28 #define MIN_ITER 10
  29
  30 //**********************************************************************************************************************
  31
  32 vector<string> ShhherCommand::getValidParameters(){
  33         try {
  34                 string Array[] =  {
  35                         "file", "flow", "lookup", "cutoff", "sigma", "outputdir","inputdir", "processors", "maxiter", "mindelta"
  36                 };
  37
  38                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  39                 return myArray;
  40         }
  41         catch(exception& e) {
  42                 m->errorOut(e, "ShhherCommand", "getValidParameters");
  43                 exit(1);
  44         }
  45 }
  46
  47 //**********************************************************************************************************************
  48
  49 ShhherCommand::ShhherCommand(){
  50         try {
  51                 abort = true; calledHelp = true;
  52
  53                 //initialize outputTypes
  54                 vector<string> tempOutNames;
  55                 outputTypes["pn.dist"] = tempOutNames;
  56
  57         }
  58         catch(exception& e) {
  59                 m->errorOut(e, "ShhherCommand", "ShhherCommand");
  60                 exit(1);
  61         }
  62 }
  63
  64 //**********************************************************************************************************************
  65
  66 vector<string> ShhherCommand::getRequiredParameters(){
  67         try {
  68                 string Array[] =  {"flow"};
  69                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  70                 return myArray;
  71         }
  72         catch(exception& e) {
  73                 m->errorOut(e, "ShhherCommand", "getRequiredParameters");
  74                 exit(1);
  75         }
  76 }
  77
  78 //**********************************************************************************************************************
  79
  80 vector<string> ShhherCommand::getRequiredFiles(){
  81         try {
  82                 vector<string> myArray;
  83                 return myArray;
  84         }
  85         catch(exception& e) {
  86                 m->errorOut(e, "ShhherCommand", "getRequiredFiles");
  87                 exit(1);
  88         }
  89 }
  90
  91 //**********************************************************************************************************************
  92
  93 ShhherCommand::ShhherCommand(string option) {
  94         try {
  95
  96 #ifdef USE_MPI
  97                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
  98                 MPI_Comm_size(MPI_COMM_WORLD, &ncpus);
  99
 100                 if(pid == 0){
 101 #endif
 102
 103
 104                 abort = false; calledHelp = false;
 105
 106
 107                 //allow user to run help
 108                 if(option == "help") { help(); abort = true; calledHelp = true; }
 109
 110                 else {
 111
 112                         //valid paramters for this command
 113                         string AlignArray[] =  {
 114                                 "file", "flow", "lookup", "cutoff", "sigma", "outputdir","inputdir", "processors", "maxiter", "mindelta"
 115                         };
 116
 117                         vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
 118
 119                         OptionParser parser(option);
 120                         map<string,string> parameters = parser.getParameters();
 121
 122                         ValidParameters validParameter;
 123                         map<string,string>::iterator it;
 124
 125                         //check to make sure all parameters are valid for command
 126                         for (it = parameters.begin(); it != parameters.end(); it++) {
 127                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 128                         }
 129
 130                         //initialize outputTypes
 131                         vector<string> tempOutNames;
 132                         outputTypes["pn.dist"] = tempOutNames;
 133                         //                      outputTypes["fasta"] = tempOutNames;
 134
 135                         //if the user changes the input directory command factory will send this info to us in the output parameter
 136                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 137                         if (inputDir == "not found"){   inputDir = "";          }
 138                         else {
 139                                 string path;
 140                                 it = parameters.find("flow");
 141                                 //user has given a template file
 142                                 if(it != parameters.end()){
 143                                         path = m->hasPath(it->second);
 144                                         //if the user has not given a path then, add inputdir. else leave path alone.
 145                                         if (path == "") {       parameters["flow"] = inputDir + it->second;             }
 146                                 }
 147
 148                                 it = parameters.find("lookup");
 149                                 //user has given a template file
 150                                 if(it != parameters.end()){
 151                                         path = m->hasPath(it->second);
 152                                         //if the user has not given a path then, add inputdir. else leave path alone.
 153                                         if (path == "") {       parameters["lookup"] = inputDir + it->second;           }
 154                                 }
 155
 156                                 it = parameters.find("file");
 157                                 //user has given a template file
 158                                 if(it != parameters.end()){
 159                                         path = m->hasPath(it->second);
 160                                         //if the user has not given a path then, add inputdir. else leave path alone.
 161                                         if (path == "") {       parameters["file"] = inputDir + it->second;             }
 162                                 }
 163                         }
 164
 165
 166                         //check for required parameters
 167                         flowFileName = validParameter.validFile(parameters, "flow", true);
 168                         flowFilesFileName = validParameter.validFile(parameters, "file", true);
 169                         if (flowFileName == "not found" && flowFilesFileName == "not found") {
 170                                 m->mothurOut("values for either flow or file must be provided for the shhh.seqs command.");
 171                                 m->mothurOutEndLine();
 172                                 abort = true;
 173                         }
 174                         else if (flowFileName == "not open" || flowFilesFileName == "not open") { abort = true; }
 175
 176                         if(flowFileName != "not found"){        compositeFASTAFileName = "";    }
 177                         else{
 178                                 compositeFASTAFileName = flowFilesFileName.substr(0, flowFilesFileName.length()-10) + "pn.fasta";
 179                                 ofstream temp;
 180                                 m->openOutputFile(compositeFASTAFileName, temp);
 181                                 temp.close();
 182                         }
 183
 184                         //if the user changes the output directory command factory will send this info to us in the output parameter
 185                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 186                                 outputDir = "";
 187                                 outputDir += m->hasPath(flowFileName); //if user entered a file with a path then preserve it
 188                         }
 189
 190
 191                         //check for optional parameter and set defaults
 192                         // ...at some point should added some additional type checking...
 193                         string temp;
 194                         temp = validParameter.validFile(parameters, "lookup", true);
 195                         if (temp == "not found")        {       lookupFileName = "LookUp_Titanium.pat"; }
 196                         else if(temp == "not open")     {       abort = true;                   }
 197                         else                                            {       lookupFileName = temp;  }
 198
 199                         temp = validParameter.validFile(parameters, "processors", false);if (temp == "not found"){      temp = "1";                     }
 200                         convert(temp, processors);
 201
 202                         temp = validParameter.validFile(parameters, "cutoff", false);   if (temp == "not found"){       temp = "0.01";          }
 203                         convert(temp, cutoff);
 204
 205                         temp = validParameter.validFile(parameters, "mindelta", false); if (temp == "not found"){       temp = "0.000001";      }
 206                         convert(temp, minDelta);
 207
 208                         temp = validParameter.validFile(parameters, "maxiter", false);  if (temp == "not found"){       temp = "1000";          }
 209                         convert(temp, maxIters);
 210
 211                         temp = validParameter.validFile(parameters, "sigma", false);if (temp == "not found")    {       temp = "60";            }
 212                         convert(temp, sigma);
 213
 214                         globaldata = GlobalData::getInstance();
 215                 }
 216
 217 #ifdef USE_MPI
 218                 }
 219 #endif
 220
 221         }
 222         catch(exception& e) {
 223                 m->errorOut(e, "ShhherCommand", "ShhherCommand");
 224                 exit(1);
 225         }
 226 }
 227
 228 //**********************************************************************************************************************
 229
 230 ShhherCommand::~ShhherCommand(){}
 231
 232 //**********************************************************************************************************************
 233
 234 void ShhherCommand::help(){
 235         try {
 236                 m->mothurOut("The shhher command reads a file containing flowgrams and creates a file of corrected sequences.\n");
 237         }
 238         catch(exception& e) {
 239                 m->errorOut(e, "ShhherCommand", "help");
 240                 exit(1);
 241         }
 242 }
 243
 244 //**********************************************************************************************************************
 245 #ifdef USE_MPI
 246 int ShhherCommand::execute(){
 247         try {
 248                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 249
 250                 int tag = 1976;
 251                 MPI_Status status;
 252
 253                 if(pid == 0){
 254
 255                         for(int i=1;i<ncpus;i++){
 256                                 MPI_Send(&abort, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 257                         }
 258                         if(abort == 1){ return 0;       }
 259
 260                         processors = ncpus;
 261
 262                         m->mothurOut("\nGetting preliminary data...\n");
 263                         getSingleLookUp();
 264                         getJointLookUp();
 265
 266                         vector<string> flowFileVector;
 267                         if(flowFilesFileName != "not found"){
 268                                 string fName;
 269
 270                                 ifstream flowFilesFile;
 271                                 m->openInputFile(flowFilesFileName, flowFilesFile);
 272                                 while(flowFilesFile){
 273                                         flowFilesFile >> fName;
 274                                         flowFileVector.push_back(fName);
 275                                         m->gobble(flowFilesFile);
 276                                 }
 277                         }
 278                         else{
 279                                 flowFileVector.push_back(flowFileName);
 280                         }
 281                         int numFiles = flowFileVector.size();
 282
 283                         for(int i=1;i<ncpus;i++){
 284                                 MPI_Send(&numFiles, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 285                         }
 286
 287                         for(int i=0;i<numFiles;i++){
 288                                 double begClock = clock();
 289                                 unsigned long int begTime = time(NULL);
 290
 291                                 flowFileName = flowFileVector[i];
 292
 293                                 m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(numFiles) + ")\t<<<<<\n");
 294                                 m->mothurOut("Reading flowgrams...\n");
 295                                 getFlowData();
 296
 297                                 m->mothurOut("Identifying unique flowgrams...\n");
 298                                 getUniques();
 299
 300                                 m->mothurOut("Calculating distances between flowgrams...\n");
 301                                 char fileName[1024];
 302                                 strcpy(fileName, flowFileName.c_str());
 303
 304                                 for(int i=1;i<ncpus;i++){
 305                                         MPI_Send(&fileName[0], 1024, MPI_CHAR, i, tag, MPI_COMM_WORLD);
 306
 307                                         MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 308                                         MPI_Send(&numUniques, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 309                                         MPI_Send(&numFlowCells, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 310                                         MPI_Send(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, i, tag, MPI_COMM_WORLD);
 311                                         MPI_Send(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
 312                                         MPI_Send(&mapUniqueToSeq[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
 313                                         MPI_Send(&mapSeqToUnique[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
 314                                         MPI_Send(&lengths[0], numSeqs, MPI_INT, i, tag, MPI_COMM_WORLD);
 315                                         MPI_Send(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
 316                                         MPI_Send(&cutoff, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
 317                                 }
 318
 319                                 string distFileName = flowDistMPI(0, int(sqrt(1.0/float(ncpus)) * numUniques));
 320
 321                                 int done;
 322                                 for(int i=1;i<ncpus;i++){
 323                                         MPI_Recv(&done, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
 324
 325                                         m->appendFiles((distFileName + ".temp." + toString(i)), distFileName);
 326                                         remove((distFileName + ".temp." + toString(i)).c_str());
 327                                 }
 328
 329                                 string namesFileName = createNamesFile();
 330
 331                                 m->mothurOut("\nClustering flowgrams...\n");
 332                                 string listFileName = cluster(distFileName, namesFileName);
 333
 334                                 getOTUData(listFileName);
 335                                 initPyroCluster();
 336
 337                                 for(int i=1;i<ncpus;i++){
 338                                         MPI_Send(&numOTUs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 339                                         MPI_Send(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
 340                                         MPI_Send(&uniqueFlowgrams[0], numFlowCells * numUniques, MPI_SHORT, i, tag, MPI_COMM_WORLD);
 341                                         MPI_Send(&sigma, 1, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
 342                                 }
 343
 344
 345                                 double maxDelta = 0;
 346                                 int iter = 0;
 347
 348                                 int numOTUsOnCPU = numOTUs / ncpus;
 349                                 int numSeqsOnCPU = numSeqs / ncpus;
 350                                 m->mothurOut("\nDenoising flowgrams...\n");
 351                                 m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
 352
 353                                 while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
 354
 355                                         double cycClock = clock();
 356                                         unsigned long int cycTime = time(NULL);
 357                                         fill();
 358
 359                                         int total = singleTau.size();
 360                                         for(int i=1;i<ncpus;i++){
 361                                                 MPI_Send(&total, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 362                                                 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
 363                                                 MPI_Send(&centroids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
 364
 365                                                 MPI_Send(&singleTau[0], total, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
 366                                                 MPI_Send(&seqNumber[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
 367                                                 MPI_Send(&seqIndex[0], total, MPI_INT, i, tag, MPI_COMM_WORLD);
 368                                                 MPI_Send(&nSeqsPerOTU[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
 369                                                 MPI_Send(&cumNumSeqs[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
 370                                         }
 371
 372                                         calcCentroidsDriver(0, numOTUsOnCPU);
 373
 374                                         for(int i=1;i<ncpus;i++){
 375                                                 int otuStart = i * numOTUs / ncpus;
 376                                                 int otuStop = (i + 1) * numOTUs / ncpus;
 377
 378                                                 vector<int> tempCentroids(numOTUs, 0);
 379                                                 vector<short> tempChange(numOTUs, 0);
 380
 381                                                 MPI_Recv(&tempCentroids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
 382                                                 MPI_Recv(&tempChange[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD, &status);
 383
 384                                                 for(int j=otuStart;j<otuStop;j++){
 385                                                         centroids[j] = tempCentroids[j];
 386                                                         change[j] = tempChange[j];
 387                                                 }
 388                                         }
 389
 390                                         maxDelta = getNewWeights();
 391                                         double nLL = getLikelihood();
 392                                         checkCentroids();
 393
 394                                         for(int i=1;i<ncpus;i++){
 395                                                 MPI_Send(&centroids[0], numOTUs, MPI_INT, i, tag, MPI_COMM_WORLD);
 396                                                 MPI_Send(&weight[0], numOTUs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
 397                                                 MPI_Send(&change[0], numOTUs, MPI_SHORT, i, tag, MPI_COMM_WORLD);
 398                                         }
 399
 400                                         calcNewDistancesParent(0, numSeqsOnCPU);
 401
 402                                         total = singleTau.size();
 403
 404                                         for(int i=1;i<ncpus;i++){
 405                                                 int childTotal;
 406                                                 int seqStart = i * numSeqs / ncpus;
 407                                                 int seqStop = (i + 1) * numSeqs / ncpus;
 408
 409                                                 MPI_Recv(&childTotal, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
 410
 411                                                 vector<int> childSeqIndex(childTotal, 0);
 412                                                 vector<double> childSingleTau(childTotal, 0);
 413                                                 vector<double> childDist(numSeqs * numOTUs, 0);
 414                                                 vector<int> otuIndex(childTotal, 0);
 415
 416                                                 MPI_Recv(&childSeqIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
 417                                                 MPI_Recv(&childSingleTau[0], childTotal, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
 418                                                 MPI_Recv(&childDist[0], numOTUs * numSeqs, MPI_DOUBLE, i, tag, MPI_COMM_WORLD, &status);
 419                                                 MPI_Recv(&otuIndex[0], childTotal, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
 420
 421                                                 int oldTotal = total;
 422                                                 total += childTotal;
 423                                                 singleTau.resize(total, 0);
 424                                                 seqIndex.resize(total, 0);
 425                                                 seqNumber.resize(total, 0);
 426
 427                                                 int childIndex = 0;
 428
 429                                                 for(int j=oldTotal;j<total;j++){
 430                                                         int otuI = otuIndex[childIndex];
 431                                                         int seqI = childSeqIndex[childIndex];
 432
 433                                                         singleTau[j] = childSingleTau[childIndex];
 434
 435                                                         aaP[otuI][nSeqsPerOTU[otuI]] = j;
 436                                                         aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
 437                                                         nSeqsPerOTU[otuI]++;
 438                                                         childIndex++;
 439                                                 }
 440
 441                                                 int index = seqStart * numOTUs;
 442                                                 for(int j=seqStart;j<seqStop;j++){
 443                                                         for(int k=0;k<numOTUs;k++){
 444                                                                 dist[index] = childDist[index];
 445                                                                 index++;
 446                                                         }
 447                                                 }
 448                                         }
 449
 450                                         iter++;
 451
 452                                         m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
 453
 454                                         if((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
 455                                                 int live = 1;
 456                                                 for(int i=1;i<ncpus;i++){
 457                                                         MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 458                                                 }
 459                                         }
 460                                         else{
 461                                                 int live = 0;
 462                                                 for(int i=1;i<ncpus;i++){
 463                                                         MPI_Send(&live, 1, MPI_INT, i, tag, MPI_COMM_WORLD); //send kill command
 464                                                 }
 465                                         }
 466
 467                                 }
 468
 469                                 m->mothurOut("\nFinalizing...\n");
 470                                 fill();
 471                                 setOTUs();
 472                                 vector<int> otuCounts(numOTUs, 0);
 473                                 for(int i=0;i<numSeqs;i++)      {       otuCounts[otuData[i]]++;        }
 474                                 calcCentroidsDriver(0, numOTUs);
 475                                 writeQualities(otuCounts);
 476                                 writeSequences(otuCounts);
 477                                 writeNames(otuCounts);
 478                                 writeClusters(otuCounts);
 479                                 writeGroups();
 480
 481                                 remove(distFileName.c_str());
 482                                 remove(namesFileName.c_str());
 483                                 remove(listFileName.c_str());
 484
 485                                 m->mothurOut("Total time to process " + toString(flowFileName) + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
 486                         }
 487                 }
 488                 else{
 489                         int abort = 1;
 490
 491                         MPI_Recv(&abort, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 492                         if(abort){      return 0;       }
 493
 494                         int numFiles;
 495                         MPI_Recv(&numFiles, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 496
 497                         for(int i=0;i<numFiles;i++){
 498                                 //Now into the pyrodist part
 499                                 bool live = 1;
 500
 501                                 char fileName[1024];
 502                                 MPI_Recv(&fileName, 1024, MPI_CHAR, 0, tag, MPI_COMM_WORLD, &status);
 503                                 MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 504                                 MPI_Recv(&numUniques, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 505                                 MPI_Recv(&numFlowCells, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 506
 507                                 flowDataIntI.resize(numSeqs * numFlowCells);
 508                                 flowDataPrI.resize(numSeqs * numFlowCells);
 509                                 mapUniqueToSeq.resize(numSeqs);
 510                                 mapSeqToUnique.resize(numSeqs);
 511                                 lengths.resize(numSeqs);
 512                                 jointLookUp.resize(NUMBINS * NUMBINS);
 513
 514                                 MPI_Recv(&flowDataIntI[0], numSeqs * numFlowCells, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
 515                                 MPI_Recv(&flowDataPrI[0], numSeqs * numFlowCells, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
 516                                 MPI_Recv(&mapUniqueToSeq[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 517                                 MPI_Recv(&mapSeqToUnique[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 518                                 MPI_Recv(&lengths[0], numSeqs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 519                                 MPI_Recv(&jointLookUp[0], NUMBINS * NUMBINS, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
 520                                 MPI_Recv(&cutoff, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
 521
 522                                 flowFileName = string(fileName);
 523                                 int flowDistStart = int(sqrt(float(pid)/float(ncpus)) * numUniques);
 524                                 int flowDistEnd = int(sqrt(float(pid+1)/float(ncpus)) * numUniques);
 525
 526                                 string distanceStringChild = flowDistMPI(flowDistStart, flowDistEnd);
 527
 528                                 int done = 1;
 529                                 MPI_Send(&done, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
 530
 531                                 //Now into the pyronoise part
 532                                 MPI_Recv(&numOTUs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 533
 534                                 singleLookUp.resize(HOMOPS * NUMBINS);
 535                                 uniqueFlowgrams.resize(numUniques * numFlowCells);
 536                                 weight.resize(numOTUs);
 537                                 centroids.resize(numOTUs);
 538                                 change.resize(numOTUs);
 539                                 dist.assign(numOTUs * numSeqs, 0);
 540                                 nSeqsPerOTU.resize(numOTUs);
 541                                 cumNumSeqs.resize(numOTUs);
 542
 543                                 MPI_Recv(&singleLookUp[0], singleLookUp.size(), MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
 544                                 MPI_Recv(&uniqueFlowgrams[0], uniqueFlowgrams.size(), MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
 545                                 MPI_Recv(&sigma, 1, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
 546
 547                                 int startOTU = pid * numOTUs / ncpus;
 548                                 int endOTU = (pid + 1) * numOTUs / ncpus;
 549
 550                                 int startSeq = pid * numSeqs / ncpus;
 551                                 int endSeq = (pid + 1) * numSeqs /ncpus;
 552
 553                                 int total;
 554
 555                                 while(live){
 556
 557                                         MPI_Recv(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 558                                         singleTau.assign(total, 0.0000);
 559                                         seqNumber.assign(total, 0);
 560                                         seqIndex.assign(total, 0);
 561
 562                                         MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
 563                                         MPI_Recv(&centroids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 564                                         MPI_Recv(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
 565                                         MPI_Recv(&seqNumber[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
 566                                         MPI_Recv(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 567                                         MPI_Recv(&nSeqsPerOTU[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 568                                         MPI_Recv(&cumNumSeqs[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 569
 570                                         calcCentroidsDriver(startOTU, endOTU);
 571
 572                                         MPI_Send(&centroids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD);
 573                                         MPI_Send(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD);
 574
 575                                         MPI_Recv(&centroids[0], numOTUs, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 576                                         MPI_Recv(&weight[0], numOTUs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD, &status);
 577                                         MPI_Recv(&change[0], numOTUs, MPI_SHORT, 0, tag, MPI_COMM_WORLD, &status);
 578
 579                                         vector<int> otuIndex(total, 0);
 580                                         calcNewDistancesChildMPI(startSeq, endSeq, otuIndex);
 581                                         total = otuIndex.size();
 582
 583                                         MPI_Send(&total, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
 584                                         MPI_Send(&seqIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
 585                                         MPI_Send(&singleTau[0], total, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
 586                                         MPI_Send(&dist[0], numOTUs * numSeqs, MPI_DOUBLE, 0, tag, MPI_COMM_WORLD);
 587                                         MPI_Send(&otuIndex[0], total, MPI_INT, 0, tag, MPI_COMM_WORLD);
 588
 589                                         MPI_Recv(&live, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 590                                 }
 591                         }
 592                 }
 593                 MPI_Barrier(MPI_COMM_WORLD);
 594
 595                 return 0;
 596
 597         }
 598         catch(exception& e) {
 599                 m->errorOut(e, "ShhherCommand", "execute");
 600                 exit(1);
 601         }
 602 }
 603
 604 /**************************************************************************************************/
 605
 606 string ShhherCommand::flowDistMPI(int startSeq, int stopSeq){
 607         try{
 608                 ostringstream outStream;
 609                 outStream.setf(ios::fixed, ios::floatfield);
 610                 outStream.setf(ios::dec, ios::basefield);
 611                 outStream.setf(ios::showpoint);
 612                 outStream.precision(6);
 613
 614                 int begTime = time(NULL);
 615                 double begClock = clock();
 616
 617                 for(int i=startSeq;i<stopSeq;i++){
 618                         for(int j=0;j<i;j++){
 619                                 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
 620
 621                                 if(flowDistance < 1e-6){
 622                                         outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
 623                                 }
 624                                 else if(flowDistance <= cutoff){
 625                                         outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
 626                                 }
 627                         }
 628                         if(i % 100 == 0){
 629                                 m->mothurOut(toString(i) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
 630                         }
 631                 }
 632
 633                 m->mothurOut(toString(stopSeq) + '\t' + toString(time(NULL) - begTime) + '\t' + toString((clock()-begClock)/CLOCKS_PER_SEC) + '\n');
 634
 635                 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.dist";
 636                 if(pid != 0){   fDistFileName += ".temp." + toString(pid);      }
 637
 638                 ofstream distFile(fDistFileName.c_str());
 639                 distFile << outStream.str();
 640                 distFile.close();
 641
 642                 return fDistFileName;
 643         }
 644         catch(exception& e) {
 645                 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
 646                 exit(1);
 647         }
 648 }
 649
 650 #else
 651 //**********************************************************************************************************************
 652
 653 int ShhherCommand::execute(){
 654         try {
 655                 if (abort == true) { return 0; }
 656
 657                 getSingleLookUp();
 658                 getJointLookUp();
 659
 660                 vector<string> flowFileVector;
 661                 if(flowFilesFileName != "not found"){
 662                         string fName;
 663
 664                         ifstream flowFilesFile;
 665                         m->openInputFile(flowFilesFileName, flowFilesFile);
 666                         while(flowFilesFile){
 667                                 flowFilesFile >> fName;
 668                                 flowFileVector.push_back(fName);
 669                                 m->gobble(flowFilesFile);
 670                         }
 671                 }
 672                 else{
 673                         flowFileVector.push_back(flowFileName);
 674                 }
 675                 int numFiles = flowFileVector.size();
 676
 677
 678                 for(int i=0;i<numFiles;i++){
 679                         flowFileName = flowFileVector[i];
 680
 681                         m->mothurOut("\n>>>>>\tProcessing " + flowFileName + " (file " + toString(i+1) + " of " + toString(numFiles) + ")\t<<<<<\n");
 682                         m->mothurOut("Reading flowgrams...\n");
 683                         getFlowData();
 684
 685                         m->mothurOut("Identifying unique flowgrams...\n");
 686                         getUniques();
 687
 688
 689                         m->mothurOut("Calculating distances between flowgrams...\n");
 690                         string distFileName = createDistFile(processors);
 691                         string namesFileName = createNamesFile();
 692
 693                         m->mothurOut("\nClustering flowgrams...\n");
 694                         string listFileName = cluster(distFileName, namesFileName);
 695                         getOTUData(listFileName);
 696
 697                         initPyroCluster();
 698
 699                         double maxDelta = 0;
 700                         int iter = 0;
 701
 702                         double begClock = clock();
 703                         unsigned long int begTime = time(NULL);
 704
 705                         m->mothurOut("\nDenoising flowgrams...\n");
 706                         m->mothurOut("iter\tmaxDelta\tnLL\t\tcycletime\n");
 707
 708                         while((maxIters == 0 && maxDelta > minDelta) || iter < MIN_ITER || (maxDelta > minDelta && iter < maxIters)){
 709
 710                                 double cycClock = clock();
 711                                 unsigned long int cycTime = time(NULL);
 712                                 fill();
 713
 714                                 calcCentroids();
 715
 716                                 maxDelta = getNewWeights();
 717                                 double nLL = getLikelihood();
 718                                 checkCentroids();
 719
 720                                 calcNewDistances();
 721
 722                                 iter++;
 723
 724                                 m->mothurOut(toString(iter) + '\t' + toString(maxDelta) + '\t' + toString(nLL) + '\t' + toString(time(NULL) - cycTime) + '\t' + toString((clock() - cycClock)/(double)CLOCKS_PER_SEC) + '\n');
 725
 726                         }
 727
 728                         m->mothurOut("\nFinalizing...\n");
 729                         fill();
 730                         setOTUs();
 731
 732                         vector<int> otuCounts(numOTUs, 0);
 733                         for(int i=0;i<numSeqs;i++)      {       otuCounts[otuData[i]]++;        }
 734
 735                         calcCentroidsDriver(0, numOTUs);
 736                         writeQualities(otuCounts);
 737                         writeSequences(otuCounts);
 738                         writeNames(otuCounts);
 739                         writeClusters(otuCounts);
 740                         writeGroups();
 741
 742                         remove(distFileName.c_str());
 743                         remove(namesFileName.c_str());
 744                         remove(listFileName.c_str());
 745
 746                         m->mothurOut("Total time to process " + flowFileName + ":\t" + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/(double)CLOCKS_PER_SEC) + '\n');
 747                 }
 748                 return 0;
 749         }
 750         catch(exception& e) {
 751                 m->errorOut(e, "ShhherCommand", "execute");
 752                 exit(1);
 753         }
 754 }
 755 #endif
 756 /**************************************************************************************************/
 757
 758 void ShhherCommand::getFlowData(){
 759         try{
 760                 ifstream flowFile;
 761                 m->openInputFile(flowFileName, flowFile);
 762
 763                 string seqName;
 764                 seqNameVector.clear();
 765                 lengths.clear();
 766                 flowDataIntI.clear();
 767                 nameMap.clear();
 768
 769
 770                 int currentNumFlowCells;
 771
 772                 float intensity;
 773
 774                 flowFile >> numFlowCells;
 775                 int index = 0;//pcluster
 776                 while(!flowFile.eof()){
 777                         flowFile >> seqName >> currentNumFlowCells;
 778                         lengths.push_back(currentNumFlowCells);
 779
 780                         seqNameVector.push_back(seqName);
 781                         nameMap[seqName] = index++;//pcluster
 782
 783                         for(int i=0;i<numFlowCells;i++){
 784                                 flowFile >> intensity;
 785                                 if(intensity > 9.99)    {       intensity = 9.99;       }
 786                                 int intI = int(100 * intensity + 0.0001);
 787                                 flowDataIntI.push_back(intI);
 788                         }
 789                         m->gobble(flowFile);
 790                 }
 791                 flowFile.close();
 792
 793                 numSeqs = seqNameVector.size();
 794
 795                 for(int i=0;i<numSeqs;i++){
 796                         int iNumFlowCells = i * numFlowCells;
 797                         for(int j=lengths[i];j<numFlowCells;j++){
 798                                 flowDataIntI[iNumFlowCells + j] = 0;
 799                         }
 800                 }
 801
 802         }
 803         catch(exception& e) {
 804                 m->errorOut(e, "ShhherCommand", "getFlowData");
 805                 exit(1);
 806         }
 807 }
 808
 809 /**************************************************************************************************/
 810
 811 void ShhherCommand::getSingleLookUp(){
 812         try{
 813                 //      these are the -log probabilities that a signal corresponds to a particular homopolymer length
 814                 singleLookUp.assign(HOMOPS * NUMBINS, 0);
 815
 816                 int index = 0;
 817                 ifstream lookUpFile;
 818                 m->openInputFile(lookupFileName, lookUpFile);
 819
 820                 for(int i=0;i<HOMOPS;i++){
 821                         float logFracFreq;
 822                         lookUpFile >> logFracFreq;
 823
 824                         for(int j=0;j<NUMBINS;j++)      {
 825                                 lookUpFile >> singleLookUp[index];
 826                                 index++;
 827                         }
 828                 }
 829                 lookUpFile.close();
 830         }
 831         catch(exception& e) {
 832                 m->errorOut(e, "ShhherCommand", "getSingleLookUp");
 833                 exit(1);
 834         }
 835 }
 836
 837 /**************************************************************************************************/
 838
 839 void ShhherCommand::getJointLookUp(){
 840         try{
 841
 842                 //      the most likely joint probability (-log) that two intenities have the same polymer length
 843                 jointLookUp.resize(NUMBINS * NUMBINS, 0);
 844
 845                 for(int i=0;i<NUMBINS;i++){
 846                         for(int j=0;j<NUMBINS;j++){
 847
 848                                 double minSum = 100000000;
 849
 850                                 for(int k=0;k<HOMOPS;k++){
 851                                         double sum = singleLookUp[k * NUMBINS + i] + singleLookUp[k * NUMBINS + j];
 852
 853                                         if(sum < minSum)        {       minSum = sum;           }
 854                                 }
 855                                 jointLookUp[i * NUMBINS + j] = minSum;
 856                         }
 857                 }
 858         }
 859         catch(exception& e) {
 860                 m->errorOut(e, "ShhherCommand", "getJointLookUp");
 861                 exit(1);
 862         }
 863 }
 864
 865 /**************************************************************************************************/
 866
 867 double ShhherCommand::getProbIntensity(int intIntensity){
 868         try{
 869
 870                 double minNegLogProb = 100000000;
 871
 872
 873                 for(int i=0;i<HOMOPS;i++){//loop signal strength
 874                         float negLogProb = singleLookUp[i * NUMBINS + intIntensity];
 875                         if(negLogProb < minNegLogProb)  {       minNegLogProb = negLogProb; }
 876                 }
 877
 878                 return minNegLogProb;
 879         }
 880         catch(exception& e) {
 881                 m->errorOut(e, "ShhherCommand", "getProbIntensity");
 882                 exit(1);
 883         }
 884 }
 885
 886 /**************************************************************************************************/
 887
 888 void ShhherCommand::getUniques(){
 889         try{
 890
 891
 892                 numUniques = 0;
 893                 uniqueFlowgrams.assign(numFlowCells * numSeqs, -1);
 894                 uniqueCount.assign(numSeqs, 0);                                                 //      anWeights
 895                 uniqueLengths.assign(numSeqs, 0);
 896                 mapSeqToUnique.assign(numSeqs, -1);
 897                 mapUniqueToSeq.assign(numSeqs, -1);
 898
 899                 vector<short> uniqueFlowDataIntI(numFlowCells * numSeqs, -1);
 900
 901                 for(int i=0;i<numSeqs;i++){
 902                         int index = 0;
 903
 904                         vector<short> current(numFlowCells);
 905                         for(int j=0;j<numFlowCells;j++){
 906                                 current[j] = short(((flowDataIntI[i * numFlowCells + j] + 50.0)/100.0));
 907                         }
 908
 909                         for(int j=0;j<numUniques;j++){
 910                                 int offset = j * numFlowCells;
 911                                 bool toEnd = 1;
 912
 913                                 int shorterLength;
 914                                 if(lengths[i] < uniqueLengths[j])       {       shorterLength = lengths[i];                     }
 915                                 else                                                            {       shorterLength = uniqueLengths[j];       }
 916
 917                                 for(int k=0;k<shorterLength;k++){
 918                                         if(current[k] != uniqueFlowgrams[offset + k]){
 919                                                 toEnd = 0;
 920                                                 break;
 921                                         }
 922                                 }
 923
 924                                 if(toEnd){
 925                                         mapSeqToUnique[i] = j;
 926                                         uniqueCount[j]++;
 927                                         index = j;
 928                                         if(lengths[i] > uniqueLengths[j])       {       uniqueLengths[j] = lengths[i];  }
 929                                         break;
 930                                 }
 931                                 index++;
 932                         }
 933
 934                         if(index == numUniques){
 935                                 uniqueLengths[numUniques] = lengths[i];
 936                                 uniqueCount[numUniques] = 1;
 937                                 mapSeqToUnique[i] = numUniques;//anMap
 938                                 mapUniqueToSeq[numUniques] = i;//anF
 939
 940                                 for(int k=0;k<numFlowCells;k++){
 941                                         uniqueFlowgrams[numUniques * numFlowCells + k] = current[k];
 942                                         uniqueFlowDataIntI[numUniques * numFlowCells + k] = flowDataIntI[i * numFlowCells + k];
 943                                 }
 944
 945                                 numUniques++;
 946                         }
 947                 }
 948                 uniqueFlowDataIntI.resize(numFlowCells * numUniques);
 949                 uniqueLengths.resize(numUniques);
 950
 951                 flowDataPrI.resize(numSeqs * numFlowCells, 0);
 952                 for(int i=0;i<flowDataPrI.size();i++)   {       flowDataPrI[i] = getProbIntensity(flowDataIntI[i]);             }
 953         }
 954         catch(exception& e) {
 955                 m->errorOut(e, "ShhherCommand", "getUniques");
 956                 exit(1);
 957         }
 958 }
 959
 960 /**************************************************************************************************/
 961
 962 float ShhherCommand::calcPairwiseDist(int seqA, int seqB){
 963         try{
 964                 int minLength = lengths[mapSeqToUnique[seqA]];
 965                 if(lengths[seqB] < minLength){  minLength = lengths[mapSeqToUnique[seqB]];      }
 966
 967                 int ANumFlowCells = seqA * numFlowCells;
 968                 int BNumFlowCells = seqB * numFlowCells;
 969
 970                 float dist = 0;
 971
 972                 for(int i=0;i<minLength;i++){
 973                         int flowAIntI = flowDataIntI[ANumFlowCells + i];
 974                         float flowAPrI = flowDataPrI[ANumFlowCells + i];
 975
 976                         int flowBIntI = flowDataIntI[BNumFlowCells + i];
 977                         float flowBPrI = flowDataPrI[BNumFlowCells + i];
 978                         dist += jointLookUp[flowAIntI * NUMBINS + flowBIntI] - flowAPrI - flowBPrI;
 979                 }
 980
 981                 dist /= (float) minLength;
 982                 return dist;
 983         }
 984         catch(exception& e) {
 985                 m->errorOut(e, "ShhherCommand", "calcPairwiseDist");
 986                 exit(1);
 987         }
 988 }
 989
 990 /**************************************************************************************************/
 991
 992 void ShhherCommand::flowDistParentFork(string distFileName, int startSeq, int stopSeq){
 993         try{
 994
 995                 ostringstream outStream;
 996                 outStream.setf(ios::fixed, ios::floatfield);
 997                 outStream.setf(ios::dec, ios::basefield);
 998                 outStream.setf(ios::showpoint);
 999                 outStream.precision(6);
1000
1001                 int begTime = time(NULL);
1002                 double begClock = clock();
1003
1004                 for(int i=startSeq;i<stopSeq;i++){
1005                         for(int j=0;j<i;j++){
1006                                 float flowDistance = calcPairwiseDist(mapUniqueToSeq[i], mapUniqueToSeq[j]);
1007
1008                                 if(flowDistance < 1e-6){
1009                                         outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << 0.000000 << endl;
1010                                 }
1011                                 else if(flowDistance <= cutoff){
1012                                         outStream << mapUniqueToSeq[i] << '\t' << mapUniqueToSeq[j] << '\t' << flowDistance << endl;
1013                                 }
1014                         }
1015                         if(i % 100 == 0){
1016                                 m->mothurOut(toString(i) + "\t" + toString(time(NULL) - begTime));
1017                                 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
1018                                 m->mothurOutEndLine();
1019                         }
1020                 }
1021                 m->mothurOut(toString(stopSeq-1) + "\t" + toString(time(NULL) - begTime));
1022                 m->mothurOut("\t" + toString((clock()-begClock)/CLOCKS_PER_SEC));
1023                 m->mothurOutEndLine();
1024
1025                 ofstream distFile(distFileName.c_str());
1026                 distFile << outStream.str();
1027                 distFile.close();
1028         }
1029         catch(exception& e) {
1030                 m->errorOut(e, "ShhherCommand", "flowDistParentFork");
1031                 exit(1);
1032         }
1033 }
1034
1035 /**************************************************************************************************/
1036
1037 string ShhherCommand::createDistFile(int processors){
1038         try{
1039                 string fDistFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.dist";
1040
1041                 unsigned long int begTime = time(NULL);
1042                 double begClock = clock();
1043
1044                 vector<int> start;
1045                 vector<int> end;
1046
1047 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1048                 if(processors == 1)     {       flowDistParentFork(fDistFileName, 0, numUniques);               }
1049                 else{ //you have multiple processors
1050
1051                         if (numSeqs < processors){      processors = 1; }
1052
1053                         vector<int> start(processors, 0);
1054                         vector<int> end(processors, 0);
1055
1056                         for (int i = 0; i < processors; i++) {
1057                                 start[i] = int(sqrt(float(i)/float(processors)) * numUniques);
1058                                 end[i] = int(sqrt(float(i+1)/float(processors)) * numUniques);
1059                         }
1060
1061                         int process = 1;
1062                         vector<int> processIDs;
1063
1064                         //loop through and create all the processes you want
1065                         while (process != processors) {
1066                                 int pid = fork();
1067
1068                                 if (pid > 0) {
1069                                         processIDs.push_back(pid);  //create map from line number to pid so you can append files in correct order later
1070                                         process++;
1071                                 }else if (pid == 0){
1072                                         flowDistParentFork(fDistFileName + toString(getpid()) + ".temp", start[process], end[process]);
1073                                         exit(0);
1074                                 }else {
1075                                         m->mothurOut("[ERROR]: unable to spawn the necessary processes. Error code: " + toString(pid)); m->mothurOutEndLine();
1076                                         perror(" : ");
1077                                         for (int i=0;i<processIDs.size();i++) {  int temp = processIDs[i]; kill (temp, SIGINT); }
1078                                         exit(0);
1079                                 }
1080                         }
1081
1082                         //parent does its part
1083                         flowDistParentFork(fDistFileName, start[0], end[0]);
1084
1085                         //force parent to wait until all the processes are done
1086                         for (int i=0;i<processIDs.size();i++) {
1087                                 int temp = processIDs[i];
1088                                 wait(&temp);
1089                         }
1090
1091                         //append and remove temp files
1092                         for (int i=0;i<processIDs.size();i++) {
1093                                 m->appendFiles((fDistFileName + toString(processIDs[i]) + ".temp"), fDistFileName);
1094                                 remove((fDistFileName + toString(processIDs[i]) + ".temp").c_str());
1095                         }
1096
1097                 }
1098
1099 #else
1100                 flowDistParentFork(fDistFileName, 0, numUniques);
1101 #endif
1102
1103                 m->mothurOutEndLine();
1104
1105                 m->mothurOut("Total time: " + toString(time(NULL) - begTime) + '\t' + toString((clock() - begClock)/CLOCKS_PER_SEC) + '\n');
1106
1107
1108                 return fDistFileName;
1109         }
1110         catch(exception& e) {
1111                 m->errorOut(e, "ShhherCommand", "createDistFile");
1112                 exit(1);
1113         }
1114
1115 }
1116
1117 /**************************************************************************************************/
1118
1119 string ShhherCommand::createNamesFile(){
1120         try{
1121
1122                 vector<string> duplicateNames(numUniques, "");
1123                 for(int i=0;i<numSeqs;i++){
1124                         duplicateNames[mapSeqToUnique[i]] += seqNameVector[i] + ',';
1125                 }
1126
1127                 string nameFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.names";
1128
1129                 ofstream nameFile;
1130                 m->openOutputFile(nameFileName, nameFile);
1131
1132                 for(int i=0;i<numUniques;i++){
1133 //                      nameFile << seqNameVector[mapUniqueToSeq[i]] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
1134                         nameFile << mapUniqueToSeq[i] << '\t' << duplicateNames[i].substr(0, duplicateNames[i].find_last_of(',')) << endl;
1135                 }
1136
1137                 nameFile.close();
1138                 return  nameFileName;
1139         }
1140         catch(exception& e) {
1141                 m->errorOut(e, "ShhherCommand", "createNamesFile");
1142                 exit(1);
1143         }
1144 }
1145
1146 //**********************************************************************************************************************
1147
1148 string ShhherCommand::cluster(string distFileName, string namesFileName){
1149         try {
1150
1151
1152                 globaldata->setNameFile(namesFileName);
1153                 globaldata->setColumnFile(distFileName);
1154                 globaldata->setFormat("column");
1155
1156                 ReadMatrix* read = new ReadColumnMatrix(distFileName);
1157                 read->setCutoff(cutoff);
1158
1159                 NameAssignment* clusterNameMap = new NameAssignment(namesFileName);
1160                 clusterNameMap->readMap();
1161                 read->read(clusterNameMap);
1162
1163                 ListVector* list = read->getListVector();
1164                 SparseMatrix* matrix = read->getMatrix();
1165
1166                 delete read;
1167                 delete clusterNameMap;
1168
1169                 RAbundVector* rabund = new RAbundVector(list->getRAbundVector());
1170
1171                 Cluster* cluster = new CompleteLinkage(rabund, list, matrix, cutoff, "furthest");
1172                 string tag = cluster->getTag();
1173
1174                 double clusterCutoff = cutoff;
1175                 while (matrix->getSmallDist() <= clusterCutoff && matrix->getNNodes() > 0){
1176                         cluster->update(clusterCutoff);
1177                 }
1178
1179                 list->setLabel(toString(cutoff));
1180
1181                 string listFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.list";
1182                 ofstream listFile;
1183                 m->openOutputFile(listFileName, listFile);
1184                 list->print(listFile);
1185                 listFile.close();
1186
1187                 delete matrix;  delete cluster; delete rabund; delete list;
1188
1189                 return listFileName;
1190         }
1191         catch(exception& e) {
1192                 m->errorOut(e, "ShhherCommand", "cluster");
1193                 exit(1);
1194         }
1195 }
1196
1197 /**************************************************************************************************/
1198
1199 void ShhherCommand::getOTUData(string listFileName){
1200         try {
1201
1202                 ifstream listFile;
1203                 m->openInputFile(listFileName, listFile);
1204                 string label;
1205
1206                 listFile >> label >> numOTUs;
1207
1208                 otuData.assign(numSeqs, 0);
1209                 cumNumSeqs.assign(numOTUs, 0);
1210                 nSeqsPerOTU.assign(numOTUs, 0);
1211                 aaP.clear();aaP.resize(numOTUs);
1212
1213                 seqNumber.clear();
1214                 aaI.clear();
1215                 seqIndex.clear();
1216
1217                 string singleOTU = "";
1218
1219                 for(int i=0;i<numOTUs;i++){
1220
1221                         listFile >> singleOTU;
1222
1223                         istringstream otuString(singleOTU);
1224
1225                         while(otuString){
1226
1227                                 string seqName = "";
1228
1229                                 for(int j=0;j<singleOTU.length();j++){
1230                                         char letter = otuString.get();
1231
1232                                         if(letter != ','){
1233                                                 seqName += letter;
1234                                         }
1235                                         else{
1236                                                 map<string,int>::iterator nmIt = nameMap.find(seqName);
1237                                                 int index = nmIt->second;
1238
1239                                                 nameMap.erase(nmIt);
1240
1241                                                 otuData[index] = i;
1242                                                 nSeqsPerOTU[i]++;
1243                                                 aaP[i].push_back(index);
1244                                                 seqName = "";
1245                                         }
1246                                 }
1247
1248                                 map<string,int>::iterator nmIt = nameMap.find(seqName);
1249
1250                                 int index = nmIt->second;
1251                                 nameMap.erase(nmIt);
1252
1253                                 otuData[index] = i;
1254                                 nSeqsPerOTU[i]++;
1255                                 aaP[i].push_back(index);
1256
1257                                 otuString.get();
1258                         }
1259
1260                         sort(aaP[i].begin(), aaP[i].end());
1261                         for(int j=0;j<nSeqsPerOTU[i];j++){
1262                                 seqNumber.push_back(aaP[i][j]);
1263                         }
1264                         for(int j=nSeqsPerOTU[i];j<numSeqs;j++){
1265                                 aaP[i].push_back(0);
1266                         }
1267
1268
1269                 }
1270
1271                 for(int i=1;i<numOTUs;i++){
1272                         cumNumSeqs[i] = cumNumSeqs[i-1] + nSeqsPerOTU[i-1];
1273                 }
1274                 aaI = aaP;
1275                 seqIndex = seqNumber;
1276
1277                 listFile.close();
1278
1279         }
1280         catch(exception& e) {
1281                 m->errorOut(e, "ShhherCommand", "getOTUData");
1282                 exit(1);
1283         }
1284 }
1285
1286 /**************************************************************************************************/
1287
1288 void ShhherCommand::initPyroCluster(){
1289         try{
1290                 dist.assign(numSeqs * numOTUs, 0);
1291                 change.assign(numOTUs, 1);
1292                 centroids.assign(numOTUs, -1);
1293                 weight.assign(numOTUs, 0);
1294                 singleTau.assign(numSeqs, 1.0);
1295
1296                 nSeqsBreaks.assign(processors+1, 0);
1297                 nOTUsBreaks.assign(processors+1, 0);
1298
1299                 nSeqsBreaks[0] = 0;
1300                 for(int i=0;i<processors;i++){
1301                         nSeqsBreaks[i+1] = nSeqsBreaks[i] + (int)((double) numSeqs / (double) processors);
1302                         nOTUsBreaks[i+1] = nOTUsBreaks[i] + (int)((double) numOTUs / (double) processors);
1303                 }
1304                 nSeqsBreaks[processors] = numSeqs;
1305                 nOTUsBreaks[processors] = numOTUs;
1306         }
1307         catch(exception& e) {
1308                 m->errorOut(e, "ShhherCommand", "initPyroCluster");
1309                 exit(1);
1310         }
1311 }
1312
1313 /**************************************************************************************************/
1314
1315 void ShhherCommand::fill(){
1316         try {
1317                 int index = 0;
1318                 for(int i=0;i<numOTUs;i++){
1319                         cumNumSeqs[i] = index;
1320                         for(int j=0;j<nSeqsPerOTU[i];j++){
1321                                 seqNumber[index] = aaP[i][j];
1322                                 seqIndex[index] = aaI[i][j];
1323
1324                                 index++;
1325                         }
1326                 }
1327         }
1328         catch(exception& e) {
1329                 m->errorOut(e, "ShhherCommand", "fill");
1330                 exit(1);
1331         }
1332 }
1333
1334 /**************************************************************************************************/
1335
1336 void ShhherCommand::calcCentroids(){
1337         try{
1338
1339 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1340
1341                 if(processors == 1)     {
1342                         calcCentroidsDriver(0, numOTUs);
1343                 }
1344                 else{ //you have multiple processors
1345                         if (numOTUs < processors){      processors = 1; }
1346
1347                         int process = 1;
1348                         vector<int> processIDs;
1349
1350                         //loop through and create all the processes you want
1351                         while (process != processors) {
1352                                 int pid = vfork();
1353
1354                                 if (pid > 0) {
1355                                         processIDs.push_back(pid);  //create map from line number to pid so you can append files in correct order later
1356                                         process++;
1357                                 }else if (pid == 0){
1358                                         calcCentroidsDriver(nOTUsBreaks[process], nOTUsBreaks[process+1]);
1359                                         exit(0);
1360                                 }else {
1361                                         m->mothurOut("[ERROR]: unable to spawn the necessary processes. Error code: " + toString(pid)); m->mothurOutEndLine();
1362                                         perror(" : ");
1363                                         for (int i=0;i<processIDs.size();i++) {  int temp = processIDs[i]; kill (temp, SIGINT); }
1364                                         exit(0);
1365                                 }
1366                         }
1367
1368                         //parent does its part
1369                         calcCentroidsDriver(nOTUsBreaks[0], nOTUsBreaks[1]);
1370
1371                         //force parent to wait until all the processes are done
1372                         for (int i=0;i<processIDs.size();i++) {
1373                                 int temp = processIDs[i];
1374                                 wait(&temp);
1375                         }
1376                 }
1377
1378 #else
1379                 calcCentroidsDriver(0, numOTUs);
1380 #endif
1381         }
1382         catch(exception& e) {
1383                 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1384                 exit(1);
1385         }
1386 }
1387
1388 /**************************************************************************************************/
1389
1390 void ShhherCommand::calcCentroidsDriver(int start, int finish){
1391
1392         //this function gets the most likely homopolymer length at a flow position for a group of sequences
1393         //within an otu
1394
1395         try{
1396
1397
1398                 for(int i=start;i<finish;i++){
1399
1400                         double count = 0;
1401                         int position = 0;
1402                         int minFlowGram = 100000000;
1403                         double minFlowValue = 1e8;
1404                         change[i] = 0; //FALSE
1405
1406                         for(int j=0;j<nSeqsPerOTU[i];j++){
1407                                 count += singleTau[seqNumber[cumNumSeqs[i] + j]];
1408                         }
1409
1410                         if(nSeqsPerOTU[i] > 0 && count > MIN_COUNT){
1411                                 vector<double> adF(nSeqsPerOTU[i]);
1412                                 vector<int> anL(nSeqsPerOTU[i]);
1413
1414                                 for(int j=0;j<nSeqsPerOTU[i];j++){
1415                                         int index = cumNumSeqs[i] + j;
1416                                         int nI = seqIndex[index];
1417                                         int nIU = mapSeqToUnique[nI];
1418
1419                                         int k;
1420                                         for(k=0;k<position;k++){
1421                                                 if(nIU == anL[k]){
1422                                                         break;
1423                                                 }
1424                                         }
1425                                         if(k == position){
1426                                                 anL[position] = nIU;
1427                                                 adF[position] = 0.0000;
1428                                                 position++;
1429                                         }
1430                                 }
1431
1432                                 for(int j=0;j<nSeqsPerOTU[i];j++){
1433                                         int index = cumNumSeqs[i] + j;
1434                                         int nI = seqIndex[index];
1435
1436                                         double tauValue = singleTau[seqNumber[index]];
1437
1438                                         for(int k=0;k<position;k++){
1439                                                 double dist = getDistToCentroid(anL[k], nI, lengths[nI]);
1440                                                 adF[k] += dist * tauValue;
1441                                         }
1442                                 }
1443
1444                                 for(int j=0;j<position;j++){
1445                                         if(adF[j] < minFlowValue){
1446                                                 minFlowGram = j;
1447                                                 minFlowValue = adF[j];
1448                                         }
1449                                 }
1450
1451                                 if(centroids[i] != anL[minFlowGram]){
1452                                         change[i] = 1;
1453                                         centroids[i] = anL[minFlowGram];
1454                                 }
1455                         }
1456                         else if(centroids[i] != -1){
1457                                 change[i] = 1;
1458                                 centroids[i] = -1;
1459                         }
1460                 }
1461         }
1462         catch(exception& e) {
1463                 m->errorOut(e, "ShhherCommand", "calcCentroidsDriver");
1464                 exit(1);
1465         }
1466 }
1467
1468 /**************************************************************************************************/
1469
1470 double ShhherCommand::getDistToCentroid(int cent, int flow, int length){
1471         try{
1472
1473                 int flowAValue = cent * numFlowCells;
1474                 int flowBValue = flow * numFlowCells;
1475
1476                 double dist = 0;
1477
1478                 for(int i=0;i<length;i++){
1479                         dist += singleLookUp[uniqueFlowgrams[flowAValue] * NUMBINS + flowDataIntI[flowBValue]];
1480                         flowAValue++;
1481                         flowBValue++;
1482                 }
1483
1484                 return dist / (double)length;
1485         }
1486         catch(exception& e) {
1487                 m->errorOut(e, "ShhherCommand", "getDistToCentroid");
1488                 exit(1);
1489         }
1490 }
1491
1492 /**************************************************************************************************/
1493
1494 double ShhherCommand::getNewWeights(){
1495         try{
1496
1497                 double maxChange = 0;
1498
1499                 for(int i=0;i<numOTUs;i++){
1500
1501                         double difference = weight[i];
1502                         weight[i] = 0;
1503
1504                         for(int j=0;j<nSeqsPerOTU[i];j++){
1505                                 int index = cumNumSeqs[i] + j;
1506                                 double tauValue = singleTau[seqNumber[index]];
1507                                 weight[i] += tauValue;
1508                         }
1509
1510                         difference = fabs(weight[i] - difference);
1511                         if(difference > maxChange){     maxChange = difference; }
1512                 }
1513                 return maxChange;
1514         }
1515         catch(exception& e) {
1516                 m->errorOut(e, "ShhherCommand", "getNewWeights");
1517                 exit(1);
1518         }
1519 }
1520
1521 /**************************************************************************************************/
1522
1523 double ShhherCommand::getLikelihood(){
1524
1525         try{
1526
1527                 vector<long double> P(numSeqs, 0);
1528                 int effNumOTUs = 0;
1529
1530                 for(int i=0;i<numOTUs;i++){
1531                         if(weight[i] > MIN_WEIGHT){
1532                                 effNumOTUs++;
1533                         }
1534                 }
1535
1536                 string hold;
1537                 for(int i=0;i<numOTUs;i++){
1538                         for(int j=0;j<nSeqsPerOTU[i];j++){
1539                                 int index = cumNumSeqs[i] + j;
1540                                 int nI = seqIndex[index];
1541                                 double singleDist = dist[seqNumber[index]];
1542
1543                                 P[nI] += weight[i] * exp(-singleDist * sigma);
1544                         }
1545                 }
1546                 double nLL = 0.00;
1547                 for(int i=0;i<numSeqs;i++){
1548                         if(P[i] == 0){  P[i] = DBL_EPSILON;     }
1549
1550                         nLL += -log(P[i]);
1551                 }
1552
1553                 nLL = nLL -(double)numSeqs * log(sigma);
1554
1555                 return nLL;
1556         }
1557         catch(exception& e) {
1558                 m->errorOut(e, "ShhherCommand", "getNewWeights");
1559                 exit(1);
1560         }
1561 }
1562
1563 /**************************************************************************************************/
1564
1565 void ShhherCommand::checkCentroids(){
1566         try{
1567                 vector<int> unique(numOTUs, 1);
1568
1569                 for(int i=0;i<numOTUs;i++){
1570                         if(centroids[i] == -1 || weight[i] < MIN_WEIGHT){
1571                                 unique[i] = -1;
1572                         }
1573                 }
1574
1575                 for(int i=0;i<numOTUs;i++){
1576                         if(unique[i] == 1){
1577                                 for(int j=i+1;j<numOTUs;j++){
1578                                         if(unique[j] == 1){
1579
1580                                                 if(centroids[j] == centroids[i]){
1581                                                         unique[j] = 0;
1582                                                         centroids[j] = -1;
1583
1584                                                         weight[i] += weight[j];
1585                                                         weight[j] = 0.0;
1586                                                 }
1587                                         }
1588                                 }
1589                         }
1590                 }
1591         }
1592         catch(exception& e) {
1593                 m->errorOut(e, "ShhherCommand", "checkCentroids");
1594                 exit(1);
1595         }
1596 }
1597
1598 /**************************************************************************************************/
1599
1600 void ShhherCommand::calcNewDistances(){
1601         try{
1602
1603 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1604
1605                 if(processors == 1)     {
1606                         calcNewDistancesParent(0, numSeqs);
1607                 }
1608                 else{ //you have multiple processors
1609                         if (numSeqs < processors){      processors = 1; }
1610
1611                         vector<vector<int> > child_otuIndex(processors);
1612                         vector<vector<int> > child_seqIndex(processors);
1613                         vector<vector<double> > child_singleTau(processors);
1614                         vector<int> totals(processors);
1615
1616                         int process = 1;
1617                         vector<int> processIDs;
1618
1619                         //loop through and create all the processes you want
1620                         while (process != processors) {
1621                                 int pid = vfork();
1622
1623                                 if (pid > 0) {
1624                                         processIDs.push_back(pid);  //create map from line number to pid so you can append files in correct order later
1625                                         process++;
1626                                 }else if (pid == 0){
1627                                         calcNewDistancesChild(nSeqsBreaks[process], nSeqsBreaks[process+1], child_otuIndex[process], child_seqIndex[process], child_singleTau[process]);
1628                                         totals[process] = child_otuIndex[process].size();
1629
1630                                         exit(0);
1631                                 }else {
1632                                         m->mothurOut("[ERROR]: unable to spawn the necessary processes. Error code: " + toString(pid)); m->mothurOutEndLine();
1633                                         perror(" : ");
1634                                         for (int i=0;i<processIDs.size();i++) {  int temp = processIDs[i]; kill (temp, SIGINT); }
1635                                         exit(0);
1636                                 }
1637                         }
1638
1639                         //parent does its part
1640                         calcNewDistancesParent(nSeqsBreaks[0], nSeqsBreaks[1]);
1641                         int total = seqIndex.size();
1642
1643                         //force parent to wait until all the processes are done
1644                         for (int i=0;i<processIDs.size();i++) {
1645                                 int temp = processIDs[i];
1646                                 wait(&temp);
1647                         }
1648
1649                         for(int i=1;i<processors;i++){
1650                                 int oldTotal = total;
1651                                 total += totals[i];
1652
1653                                 singleTau.resize(total, 0);
1654                                 seqIndex.resize(total, 0);
1655                                 seqNumber.resize(total, 0);
1656
1657                                 int childIndex = 0;
1658
1659                                 for(int j=oldTotal;j<total;j++){
1660                                         int otuI = child_otuIndex[i][childIndex];
1661                                         int seqI = child_seqIndex[i][childIndex];
1662
1663                                         singleTau[j] = child_singleTau[i][childIndex];
1664                                         aaP[otuI][nSeqsPerOTU[otuI]] = j;
1665                                         aaI[otuI][nSeqsPerOTU[otuI]] = seqI;
1666                                         nSeqsPerOTU[otuI]++;
1667
1668                                         childIndex++;
1669                                 }
1670                         }
1671                 }
1672 #else
1673                 calcNewDistancesParent(0, numSeqs);
1674 #endif
1675         }
1676         catch(exception& e) {
1677                 m->errorOut(e, "ShhherCommand", "calcNewDistances");
1678                 exit(1);
1679         }
1680 }
1681
1682 /**************************************************************************************************/
1683 #ifdef USE_MPI
1684 void ShhherCommand::calcNewDistancesChildMPI(int startSeq, int stopSeq, vector<int>& otuIndex){
1685
1686         try{
1687                 vector<double> newTau(numOTUs,0);
1688                 vector<double> norms(numSeqs, 0);
1689                 otuIndex.clear();
1690                 seqIndex.clear();
1691                 singleTau.clear();
1692
1693
1694
1695                 for(int i=startSeq;i<stopSeq;i++){
1696                         double offset = 1e8;
1697                         int indexOffset = i * numOTUs;
1698
1699                         for(int j=0;j<numOTUs;j++){
1700
1701                                 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1702                                         dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1703                                 }
1704                                 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1705                                         offset = dist[indexOffset + j];
1706                                 }
1707                         }
1708
1709                         for(int j=0;j<numOTUs;j++){
1710                                 if(weight[j] > MIN_WEIGHT){
1711                                         newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1712                                         norms[i] += newTau[j];
1713                                 }
1714                                 else{
1715                                         newTau[j] = 0.0;
1716                                 }
1717                         }
1718
1719                         for(int j=0;j<numOTUs;j++){
1720
1721                                 newTau[j] /= norms[i];
1722
1723                                 if(newTau[j] > MIN_TAU){
1724                                         otuIndex.push_back(j);
1725                                         seqIndex.push_back(i);
1726                                         singleTau.push_back(newTau[j]);
1727                                 }
1728                         }
1729
1730                 }
1731         }
1732         catch(exception& e) {
1733                 m->errorOut(e, "ShhherCommand", "calcNewDistancesChildMPI");
1734                 exit(1);
1735         }
1736 }
1737 #endif
1738 /**************************************************************************************************/
1739
1740 void ShhherCommand::calcNewDistancesChild(int startSeq, int stopSeq, vector<int>& child_otuIndex, vector<int>& child_seqIndex, vector<double>& child_singleTau){
1741
1742         try{
1743                 vector<double> newTau(numOTUs,0);
1744                 vector<double> norms(numSeqs, 0);
1745                 child_otuIndex.resize(0);
1746                 child_seqIndex.resize(0);
1747                 child_singleTau.resize(0);
1748
1749                 for(int i=startSeq;i<stopSeq;i++){
1750                         double offset = 1e8;
1751                         int indexOffset = i * numOTUs;
1752
1753
1754                         for(int j=0;j<numOTUs;j++){
1755                                 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1756                                         dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1757                                 }
1758
1759                                 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1760                                         offset = dist[indexOffset + j];
1761                                 }
1762                         }
1763
1764                         for(int j=0;j<numOTUs;j++){
1765                                 if(weight[j] > MIN_WEIGHT){
1766                                         newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1767                                         norms[i] += newTau[j];
1768                                 }
1769                                 else{
1770                                         newTau[j] = 0.0;
1771                                 }
1772                         }
1773
1774                         for(int j=0;j<numOTUs;j++){
1775                                 newTau[j] /= norms[i];
1776
1777                                 if(newTau[j] > MIN_TAU){
1778                                         child_otuIndex.push_back(j);
1779                                         child_seqIndex.push_back(i);
1780                                         child_singleTau.push_back(newTau[j]);
1781                                 }
1782                         }
1783                 }
1784         }
1785         catch(exception& e) {
1786                 m->errorOut(e, "ShhherCommand", "calcNewDistancesChild");
1787                 exit(1);
1788         }
1789 }
1790
1791 /**************************************************************************************************/
1792
1793 void ShhherCommand::calcNewDistancesParent(int startSeq, int stopSeq){
1794
1795         try{
1796
1797                 int total = 0;
1798                 vector<double> newTau(numOTUs,0);
1799                 vector<double> norms(numSeqs, 0);
1800                 nSeqsPerOTU.assign(numOTUs, 0);
1801
1802                 for(int i=startSeq;i<stopSeq;i++){
1803                         int indexOffset = i * numOTUs;
1804
1805                         double offset = 1e8;
1806
1807                         for(int j=0;j<numOTUs;j++){
1808                                 if(weight[j] > MIN_WEIGHT && change[j] == 1){
1809                                         dist[indexOffset + j] = getDistToCentroid(centroids[j], i, lengths[i]);
1810                                 }
1811
1812                                 if(weight[j] > MIN_WEIGHT && dist[indexOffset + j] < offset){
1813                                         offset = dist[indexOffset + j];
1814                                 }
1815                         }
1816
1817                         for(int j=0;j<numOTUs;j++){
1818                                 if(weight[j] > MIN_WEIGHT){
1819                                         newTau[j] = exp(sigma * (-dist[indexOffset + j] + offset)) * weight[j];
1820                                         norms[i] += newTau[j];
1821                                 }
1822                                 else{
1823                                         newTau[j] = 0.0;
1824                                 }
1825                         }
1826
1827                         for(int j=0;j<numOTUs;j++){
1828                                 newTau[j] /= norms[i];
1829                         }
1830
1831                         for(int j=0;j<numOTUs;j++){
1832                                 if(newTau[j] > MIN_TAU){
1833
1834                                         int oldTotal = total;
1835
1836                                         total++;
1837
1838                                         singleTau.resize(total, 0);
1839                                         seqNumber.resize(total, 0);
1840                                         seqIndex.resize(total, 0);
1841
1842                                         singleTau[oldTotal] = newTau[j];
1843
1844                                         aaP[j][nSeqsPerOTU[j]] = oldTotal;
1845                                         aaI[j][nSeqsPerOTU[j]] = i;
1846                                         nSeqsPerOTU[j]++;
1847                                 }
1848                         }
1849                 }
1850         }
1851         catch(exception& e) {
1852                 m->errorOut(e, "ShhherCommand", "calcNewDistancesParent");
1853                 exit(1);
1854         }
1855 }
1856
1857 /**************************************************************************************************/
1858
1859 void ShhherCommand::setOTUs(){
1860
1861         try {
1862                 vector<double> bigTauMatrix(numOTUs * numSeqs, 0.0000);
1863
1864                 for(int i=0;i<numOTUs;i++){
1865                         for(int j=0;j<nSeqsPerOTU[i];j++){
1866                                 int index = cumNumSeqs[i] + j;
1867                                 double tauValue = singleTau[seqNumber[index]];
1868                                 int sIndex = seqIndex[index];
1869                                 bigTauMatrix[sIndex * numOTUs + i] = tauValue;
1870                         }
1871                 }
1872
1873                 for(int i=0;i<numSeqs;i++){
1874                         double maxTau = -1.0000;
1875                         int maxOTU = -1;
1876                         for(int j=0;j<numOTUs;j++){
1877                                 if(bigTauMatrix[i * numOTUs + j] > maxTau){
1878                                         maxTau = bigTauMatrix[i * numOTUs + j];
1879                                         maxOTU = j;
1880                                 }
1881                         }
1882
1883                         otuData[i] = maxOTU;
1884                 }
1885
1886                 nSeqsPerOTU.assign(numOTUs, 0);
1887
1888                 for(int i=0;i<numSeqs;i++){
1889                         int index = otuData[i];
1890
1891                         singleTau[i] = 1.0000;
1892                         dist[i] = 0.0000;
1893
1894                         aaP[index][nSeqsPerOTU[index]] = i;
1895                         aaI[index][nSeqsPerOTU[index]] = i;
1896
1897                         nSeqsPerOTU[index]++;
1898                 }
1899                 fill();
1900         }
1901         catch(exception& e) {
1902                 m->errorOut(e, "ShhherCommand", "calcNewDistances");
1903                 exit(1);
1904         }
1905 }
1906
1907 /**************************************************************************************************/
1908
1909 void ShhherCommand::writeQualities(vector<int> otuCounts){
1910
1911         try {
1912                 string qualityFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.qual";
1913
1914                 ofstream qualityFile;
1915                 m->openOutputFile(qualityFileName, qualityFile);
1916
1917                 qualityFile.setf(ios::fixed, ios::floatfield);
1918                 qualityFile.setf(ios::showpoint);
1919                 qualityFile << setprecision(6);
1920
1921                 vector<vector<int> > qualities(numOTUs);
1922                 vector<double> pr(HOMOPS, 0);
1923
1924
1925                 for(int i=0;i<numOTUs;i++){
1926                         int index = 0;
1927                         int base = 0;
1928
1929                         if(nSeqsPerOTU[i] > 0){
1930                                 qualities[i].assign(1024, -1);
1931
1932                                 while(index < numFlowCells){
1933                                         double maxPrValue = 1e8;
1934                                         short maxPrIndex = -1;
1935                                         double count = 0.0000;
1936
1937                                         pr.assign(HOMOPS, 0);
1938
1939                                         for(int j=0;j<nSeqsPerOTU[i];j++){
1940                                                 int lIndex = cumNumSeqs[i] + j;
1941                                                 double tauValue = singleTau[seqNumber[lIndex]];
1942                                                 int sequenceIndex = aaI[i][j];
1943                                                 short intensity = flowDataIntI[sequenceIndex * numFlowCells + index];
1944
1945                                                 count += tauValue;
1946
1947                                                 for(int s=0;s<HOMOPS;s++){
1948                                                         pr[s] += tauValue * singleLookUp[s * NUMBINS + intensity];
1949                                                 }
1950                                         }
1951
1952                                         maxPrIndex = uniqueFlowgrams[centroids[i] * numFlowCells + index];
1953                                         maxPrValue = pr[maxPrIndex];
1954
1955                                         if(count > MIN_COUNT){
1956                                                 double U = 0.0000;
1957                                                 double norm = 0.0000;
1958
1959                                                 for(int s=0;s<HOMOPS;s++){
1960                                                         norm += exp(-(pr[s] - maxPrValue));
1961                                                 }
1962
1963                                                 for(int s=1;s<=maxPrIndex;s++){
1964                                                         int value = 0;
1965                                                         double temp = 0.0000;
1966
1967                                                         U += exp(-(pr[s-1]-maxPrValue))/norm;
1968
1969                                                         if(U>0.00){
1970                                                                 temp = log10(U);
1971                                                         }
1972                                                         else{
1973                                                                 temp = -10.1;
1974                                                         }
1975                                                         temp = floor(-10 * temp);
1976                                                         value = (int)floor(temp);
1977                                                         if(value > 100){        value = 100;    }
1978
1979                                                         qualities[i][base] = (int)value;
1980                                                         base++;
1981                                                 }
1982                                         }
1983
1984                                         index++;
1985                                 }
1986                         }
1987
1988
1989                         if(otuCounts[i] > 0){
1990                                 qualityFile << '>' << seqNameVector[mapUniqueToSeq[i]] << endl;
1991
1992                                 int j=4;        //need to get past the first four bases
1993                                 while(qualities[i][j] != -1){
1994                                         qualityFile << qualities[i][j] << ' ';
1995                                         j++;
1996                                 }
1997                                 qualityFile << endl;
1998                         }
1999                 }
2000                 qualityFile.close();
2001
2002         }
2003         catch(exception& e) {
2004                 m->errorOut(e, "ShhherCommand", "writeQualities");
2005                 exit(1);
2006         }
2007 }
2008
2009 /**************************************************************************************************/
2010
2011 void ShhherCommand::writeSequences(vector<int> otuCounts){
2012         try {
2013                 string bases = "TACG";
2014
2015                 string fastaFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.fasta";
2016                 ofstream fastaFile;
2017                 m->openOutputFile(fastaFileName, fastaFile);
2018
2019                 vector<string> names(numOTUs, "");
2020
2021                 for(int i=0;i<numOTUs;i++){
2022                         int index = centroids[i];
2023
2024                         if(otuCounts[i] > 0){
2025                                 fastaFile << '>' << seqNameVector[aaI[i][0]] << endl;
2026
2027                                 for(int j=8;j<numFlowCells;j++){
2028
2029                                         char base = bases[j % 4];
2030                                         for(int k=0;k<uniqueFlowgrams[index * numFlowCells + j];k++){
2031                                                 fastaFile << base;
2032                                         }
2033                                 }
2034                                 fastaFile << endl;
2035                         }
2036                 }
2037                 fastaFile.close();
2038
2039                 if(compositeFASTAFileName != ""){
2040                         m->appendFiles(fastaFileName, compositeFASTAFileName);
2041                 }
2042         }
2043         catch(exception& e) {
2044                 m->errorOut(e, "ShhherCommand", "writeSequences");
2045                 exit(1);
2046         }
2047 }
2048
2049 /**************************************************************************************************/
2050
2051 void ShhherCommand::writeNames(vector<int> otuCounts){
2052         try {
2053                 string nameFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.final.names";
2054                 ofstream nameFile;
2055                 m->openOutputFile(nameFileName, nameFile);
2056
2057                 for(int i=0;i<numOTUs;i++){
2058                         if(otuCounts[i] > 0){
2059                                 nameFile << seqNameVector[aaI[i][0]] << '\t' << seqNameVector[aaI[i][0]];
2060
2061                                 for(int j=1;j<nSeqsPerOTU[i];j++){
2062                                         nameFile << ',' << seqNameVector[aaI[i][j]];
2063                                 }
2064
2065                                 nameFile << endl;
2066                         }
2067                 }
2068                 nameFile.close();
2069         }
2070         catch(exception& e) {
2071                 m->errorOut(e, "ShhherCommand", "writeNames");
2072                 exit(1);
2073         }
2074 }
2075
2076 /**************************************************************************************************/
2077
2078 void ShhherCommand::writeGroups(){
2079         try {
2080                 string fileRoot = flowFileName.substr(0,flowFileName.find_last_of('.'));
2081                 string groupFileName = fileRoot + ".pn.groups";
2082                 ofstream groupFile;
2083                 m->openOutputFile(groupFileName, groupFile);
2084
2085                 for(int i=0;i<numSeqs;i++){
2086                         groupFile << seqNameVector[i] << '\t' << fileRoot << endl;
2087                 }
2088                 groupFile.close();
2089         }
2090         catch(exception& e) {
2091                 m->errorOut(e, "ShhherCommand", "writeGroups");
2092                 exit(1);
2093         }
2094 }
2095
2096 /**************************************************************************************************/
2097
2098 void ShhherCommand::writeClusters(vector<int> otuCounts){
2099         try {
2100                 string otuCountsFileName = flowFileName.substr(0,flowFileName.find_last_of('.')) + ".pn.counts";
2101                 ofstream otuCountsFile;
2102                 m->openOutputFile(otuCountsFileName, otuCountsFile);
2103
2104                 string bases = "TACG";
2105
2106                 for(int i=0;i<numOTUs;i++){
2107                         //output the translated version of the centroid sequence for the otu
2108                         if(otuCounts[i] > 0){
2109                                 int index = centroids[i];
2110
2111                                 otuCountsFile << "ideal\t";
2112                                 for(int j=8;j<numFlowCells;j++){
2113                                         char base = bases[j % 4];
2114                                         for(int s=0;s<uniqueFlowgrams[index * numFlowCells + j];s++){
2115                                                 otuCountsFile << base;
2116                                         }
2117                                 }
2118                                 otuCountsFile << endl;
2119
2120                                 for(int j=0;j<nSeqsPerOTU[i];j++){
2121                                         int sequence = aaI[i][j];
2122                                         otuCountsFile << seqNameVector[sequence] << '\t';
2123
2124                                         for(int k=8;k<lengths[sequence];k++){
2125                                                 char base = bases[k % 4];
2126                                                 int freq = int(0.01 * (double)flowDataIntI[sequence * numFlowCells + k] + 0.5);
2127
2128                                                 for(int s=0;s<freq;s++){
2129                                                         otuCountsFile << base;
2130                                                 }
2131                                         }
2132                                         otuCountsFile << endl;
2133                                 }
2134                                 otuCountsFile << endl;
2135                         }
2136                 }
2137                 otuCountsFile.close();
2138         }
2139         catch(exception& e) {
2140                 m->errorOut(e, "ShhherCommand", "writeClusters");
2141                 exit(1);
2142         }
2143 }
2144
2145 //**********************************************************************************************************************