shhhercommand.cpp

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