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