shhhercommand.cpp

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