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