pintail.cpp

   1 /*
   2  *  pintail.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 7/9/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "pintail.h"
  11 #include "ignoregaps.h"
  12 #include "eachgapdist.h"
  13
  14 //********************************************************************************************************************
  15 //sorts lowest to highest
  16 inline bool compareQuanMembers(quanMember left, quanMember right){
  17         return (left.score < right.score);
  18 }
  19 //***************************************************************************************************************
  20
  21 Pintail::Pintail(string filename, string temp, bool f, int p, string mask, string cons, string q, int win, int inc, string o) : Chimera() {
  22         try {
  23
  24                 fastafile = filename;
  25                 templateFileName = temp; templateSeqs = readSeqs(temp);
  26                 filter = f;
  27                 processors = p;
  28                 setMask(mask);
  29                 consfile = cons;
  30                 quanfile = q;
  31                 window = win;
  32                 increment = inc;
  33                 outputDir = o;
  34
  35                 distcalculator = new eachGapDist();
  36                 decalc = new DeCalculator();
  37
  38                 doPrep();
  39         }
  40         catch(exception& e) {
  41                 m->errorOut(e, "Pintail", "Pintail");
  42                 exit(1);
  43         }
  44
  45 }
  46 //***************************************************************************************************************
  47
  48 Pintail::~Pintail() {
  49         try {
  50
  51                 delete distcalculator;
  52                 delete decalc;
  53         }
  54         catch(exception& e) {
  55                 m->errorOut(e, "Pintail", "~Pintail");
  56                 exit(1);
  57         }
  58 }
  59 //***************************************************************************************************************
  60 int Pintail::doPrep() {
  61         try {
  62
  63                 mergedFilterString = "";
  64                 windowSizesTemplate.resize(templateSeqs.size(), window);
  65                 quantiles.resize(100);  //one for every percent mismatch
  66                 quantilesMembers.resize(100);  //one for every percent mismatch
  67
  68                 //if the user does not enter a mask then you want to keep all the spots in the alignment
  69                 if (seqMask.length() == 0)      {       decalc->setAlignmentLength(templateSeqs[0]->getAligned().length());     }
  70                 else                                            {       decalc->setAlignmentLength(seqMask.length());                                           }
  71
  72                 decalc->setMask(seqMask);
  73
  74         #ifdef USE_MPI
  75                 //do nothing
  76         #else
  77                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
  78                         //find breakup of templatefile for quantiles
  79                         if (processors == 1) {   templateLines.push_back(new linePair(0, templateSeqs.size()));  }
  80                         else {
  81                                 for (int i = 0; i < processors; i++) {
  82                                         templateLines.push_back(new linePair());
  83                                         templateLines[i]->start = int (sqrt(float(i)/float(processors)) * templateSeqs.size());
  84                                         templateLines[i]->end = int (sqrt(float(i+1)/float(processors)) * templateSeqs.size());
  85                                 }
  86                         }
  87                 #else
  88                         templateLines.push_back(new linePair(0, templateSeqs.size()));
  89                 #endif
  90         #endif
  91
  92                 m->mothurOut("Getting conservation... "); cout.flush();
  93                 if (consfile == "") {
  94                         m->mothurOut("Calculating probability of conservation for your template sequences.  This can take a while...  I will output the frequency of the highest base in each position to a .freq file so that you can input them using the conservation parameter next time you run this command.  Providing the .freq file will improve speed.    "); cout.flush();
  95                         probabilityProfile = decalc->calcFreq(templateSeqs, templateFileName);
  96                         if (m->control_pressed) {  return 0;  }
  97                         m->mothurOut("Done."); m->mothurOutEndLine();
  98                 }else                           {   probabilityProfile = readFreq();    m->mothurOut("Done.");            }
  99                 m->mothurOutEndLine();
 100
 101                 //make P into Q
 102                 for (int i = 0; i < probabilityProfile.size(); i++)  { probabilityProfile[i] = 1 - probabilityProfile[i];  }  //
 103
 104                 bool reRead = false;
 105                 //create filter if needed for later
 106                 if (filter) {
 107
 108                         //read in all query seqs
 109                         vector<Sequence*> tempQuerySeqs = readSeqs(fastafile);
 110
 111                         vector<Sequence*> temp;
 112                         //merge query seqs and template seqs
 113                         temp = templateSeqs;
 114                         for (int i = 0; i < tempQuerySeqs.size(); i++) {  temp.push_back(tempQuerySeqs[i]);  }
 115
 116                         if (seqMask != "") {
 117                             reRead = true;
 118                                 //mask templates
 119                                 for (int i = 0; i < temp.size(); i++) {
 120                                         if (m->control_pressed) {
 121                                                 for (int i = 0; i < tempQuerySeqs.size(); i++) { delete tempQuerySeqs[i];  }
 122                                                 return 0;
 123                                         }
 124                                         decalc->runMask(temp[i]);
 125                                 }
 126                         }
 127
 128                         mergedFilterString = createFilter(temp, 0.5);
 129
 130                         if (m->control_pressed) {
 131                                 for (int i = 0; i < tempQuerySeqs.size(); i++) { delete tempQuerySeqs[i];  }
 132                                 return 0;
 133                         }
 134
 135                         //reread template seqs
 136                         for (int i = 0; i < tempQuerySeqs.size(); i++) { delete tempQuerySeqs[i];  }
 137                 }
 138
 139
 140                 //quantiles are used to determine whether the de values found indicate a chimera
 141                 //if you have to calculate them, its time intensive because you are finding the de and deviation values for each
 142                 //combination of sequences in the template
 143                 if (quanfile != "") {
 144                         quantiles = readQuantiles();
 145                 }else {
 146                         if ((!filter) && (seqMask != "")) { //if you didn't filter but you want to mask. if you filtered then you did mask first above.
 147                                 reRead = true;
 148                                 //mask templates
 149                                 for (int i = 0; i < templateSeqs.size(); i++) {
 150                                         if (m->control_pressed) {  return 0;  }
 151                                         decalc->runMask(templateSeqs[i]);
 152                                 }
 153                         }
 154
 155                         if (filter) {
 156                                 reRead = true;
 157                                 for (int i = 0; i < templateSeqs.size(); i++) {
 158                                         if (m->control_pressed) {  return 0;  }
 159                                         runFilter(templateSeqs[i]);
 160                                 }
 161                         }
 162
 163                         m->mothurOut("Calculating quantiles for your template.  This can take a while...  I will output the quantiles to a .quan file that you can input them using the quantiles parameter next time you run this command.  Providing the .quan file will dramatically improve speed.    "); cout.flush();
 164                         if (processors == 1) {
 165                                 quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, 0, templateSeqs.size());
 166                         }else {         createProcessesQuan();          }
 167
 168                         if (m->control_pressed) {  return 0;  }
 169
 170                         string noOutliers, outliers;
 171
 172                         if ((!filter) && (seqMask == "")) {
 173                                 noOutliers = m->getRootName(m->getSimpleName(templateFileName)) + "pintail.quan";
 174                         }else if ((!filter) && (seqMask != "")) {
 175                                 noOutliers =m->getRootName(m->getSimpleName(templateFileName)) + "pintail.masked.quan";
 176                         }else if ((filter) && (seqMask != "")) {
 177                                 noOutliers = m->getRootName(m->getSimpleName(templateFileName)) + "pintail.filtered." + m->getSimpleName(m->getRootName(fastafile)) + "masked.quan";
 178                         }else if ((filter) && (seqMask == "")) {
 179                                 noOutliers = m->getRootName(m->getSimpleName(templateFileName)) + "pintail.filtered." + m->getSimpleName(m->getRootName(fastafile)) + "quan";
 180                         }
 181
 182                         decalc->removeObviousOutliers(quantilesMembers, templateSeqs.size());
 183
 184                         if (m->control_pressed) {  return 0;  }
 185
 186                         string outputString = "#" + m->getVersion() + "\n";
 187
 188                         //adjust quantiles
 189                         for (int i = 0; i < quantilesMembers.size(); i++) {
 190                                 vector<float> temp;
 191
 192                                 if (quantilesMembers[i].size() == 0) {
 193                                         //in case this is not a distance found in your template files
 194                                         for (int g = 0; g < 6; g++) {
 195                                                 temp.push_back(0.0);
 196                                         }
 197                                 }else{
 198
 199                                         sort(quantilesMembers[i].begin(), quantilesMembers[i].end());
 200
 201                                         //save 10%
 202                                         temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.10)]);
 203                                         //save 25%
 204                                         temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.25)]);
 205                                         //save 50%
 206                                         temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.5)]);
 207                                         //save 75%
 208                                         temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.75)]);
 209                                         //save 95%
 210                                         temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.95)]);
 211                                         //save 99%
 212                                         temp.push_back(quantilesMembers[i][int(quantilesMembers[i].size() * 0.99)]);
 213
 214                                 }
 215
 216                                 //output quan value
 217                                 outputString += toString(i+1) + "\t";
 218                                 for (int u = 0; u < temp.size(); u++) {   outputString += toString(temp[u]) + "\t"; }
 219                                 outputString += "\n";
 220
 221                                 quantiles[i] = temp;
 222
 223                         }
 224
 225                         printQuanFile(noOutliers, outputString);
 226
 227                         //free memory
 228                         quantilesMembers.clear();
 229
 230                         m->mothurOut("Done."); m->mothurOutEndLine();
 231                 }
 232
 233                 if (reRead) {
 234                         for (int i = 0; i < templateSeqs.size(); i++) { delete templateSeqs[i];  }
 235                         templateSeqs.clear();
 236                         templateSeqs = readSeqs(templateFileName);
 237                 }
 238
 239
 240                 //free memory
 241                 for (int i = 0; i < templateLines.size(); i++) { delete templateLines[i];  }
 242
 243                 return 0;
 244
 245         }
 246         catch(exception& e) {
 247                 m->errorOut(e, "Pintail", "doPrep");
 248                 exit(1);
 249         }
 250 }
 251 //***************************************************************************************************************
 252 Sequence Pintail::print(ostream& out, ostream& outAcc) {
 253         try {
 254
 255                 int index = ceil(deviation);
 256
 257                 //is your DE value higher than the 95%
 258                 string chimera;
 259                 if (index != 0) {  //if index is 0 then its an exact match to a template seq
 260                         if (quantiles[index][4] == 0.0) {
 261                                 chimera = "Your template does not include sequences that provide quantile values at distance " + toString(index);
 262                         }else {
 263                                 if (DE > quantiles[index][4])           {       chimera = "Yes";        }
 264                                 else                                                            {       chimera = "No";         }
 265                         }
 266                 }else{ chimera = "No";          }
 267
 268                 out << querySeq->getName() << '\t' << "div: " << deviation << "\tstDev: " << DE << "\tchimera flag: " << chimera << endl;
 269                 if (chimera == "Yes") {
 270                         m->mothurOut(querySeq->getName() + "\tdiv: " + toString(deviation) + "\tstDev: " + toString(DE) + "\tchimera flag: " + chimera); m->mothurOutEndLine();
 271                         outAcc << querySeq->getName() << endl;
 272                 }
 273                 out << "Observed\t";
 274
 275                 for (int j = 0; j < obsDistance.size(); j++) {  out << obsDistance[j] << '\t';  }
 276                 out << endl;
 277
 278                 out << "Expected\t";
 279
 280                 for (int m = 0; m < expectedDistance.size(); m++) {  out << expectedDistance[m] << '\t';  }
 281                 out << endl;
 282
 283                 return *querySeq;
 284
 285         }
 286         catch(exception& e) {
 287                 m->errorOut(e, "Pintail", "print");
 288                 exit(1);
 289         }
 290 }
 291 #ifdef USE_MPI
 292 //***************************************************************************************************************
 293 Sequence Pintail::print(MPI_File& out, MPI_File& outAcc) {
 294         try {
 295
 296                 string outputString = "";
 297                 int index = ceil(deviation);
 298
 299                 //is your DE value higher than the 95%
 300                 string chimera;
 301                 if (index != 0) {  //if index is 0 then its an exact match to a template seq
 302                         if (quantiles[index][4] == 0.0) {
 303                                 chimera = "Your template does not include sequences that provide quantile values at distance " + toString(index);
 304                         }else {
 305                                 if (DE > quantiles[index][4])           {       chimera = "Yes";        }
 306                                 else                                                            {       chimera = "No";         }
 307                         }
 308                 }else{ chimera = "No";          }
 309
 310                 outputString += querySeq->getName() + "\tdiv: " + toString(deviation) + "\tstDev: " + toString(DE) + "\tchimera flag: " + chimera + "\n";
 311                 if (chimera == "Yes") {
 312                         cout << querySeq->getName() << "\tdiv: " << toString(deviation) << "\tstDev: " << toString(DE) << "\tchimera flag: " << chimera << endl;
 313                         string outAccString = querySeq->getName() + "\n";
 314
 315                         MPI_Status statusAcc;
 316                         int length = outAccString.length();
 317                         char* buf = new char[length];
 318                         memcpy(buf, outAccString.c_str(), length);
 319
 320                         MPI_File_write_shared(outAcc, buf, length, MPI_CHAR, &statusAcc);
 321                         delete buf;
 322
 323                         return *querySeq;
 324                 }
 325                 outputString += "Observed\t";
 326
 327                 for (int j = 0; j < obsDistance.size(); j++) {  outputString += toString(obsDistance[j]) + "\t";  }
 328                 outputString += "\n";
 329
 330                 outputString += "Expected\t";
 331
 332                 for (int m = 0; m < expectedDistance.size(); m++) {  outputString += toString(expectedDistance[m]) + "\t";  }
 333                 outputString += "\n";
 334
 335                 MPI_Status status;
 336                 int length = outputString.length();
 337                 char* buf2 = new char[length];
 338                 memcpy(buf2, outputString.c_str(), length);
 339
 340                 MPI_File_write_shared(out, buf2, length, MPI_CHAR, &status);
 341                 delete buf2;
 342
 343                 return *querySeq;
 344         }
 345         catch(exception& e) {
 346                 m->errorOut(e, "Pintail", "print");
 347                 exit(1);
 348         }
 349 }
 350 #endif
 351 //***************************************************************************************************************
 352 int Pintail::getChimeras(Sequence* query) {
 353         try {
 354                 querySeq = query;
 355                 trimmed.clear();
 356                 windowSizes = window;
 357
 358                 //find pairs has to be done before a mask
 359                 bestfit = findPairs(query);
 360
 361                 if (m->control_pressed) {  return 0; }
 362
 363                 //if they mask
 364                 if (seqMask != "") {
 365                         decalc->runMask(query);
 366                         decalc->runMask(bestfit);
 367                 }
 368
 369                 if (filter) { //must be done after a mask
 370                         runFilter(query);
 371                         runFilter(bestfit);
 372                 }
 373
 374
 375                 //trim seq
 376                 decalc->trimSeqs(query, bestfit, trimmed);
 377
 378                 //find windows
 379                 it = trimmed.begin();
 380                 windowsForeachQuery = decalc->findWindows(query, it->first, it->second, windowSizes, increment);
 381
 382                 //find observed distance
 383                 obsDistance = decalc->calcObserved(query, bestfit, windowsForeachQuery, windowSizes);
 384
 385                 if (m->control_pressed) {  return 0; }
 386
 387                 Qav = decalc->findQav(windowsForeachQuery, windowSizes, probabilityProfile);
 388
 389                 if (m->control_pressed) {  return 0; }
 390
 391                 //find alpha
 392                 seqCoef = decalc->getCoef(obsDistance, Qav);
 393
 394                 //calculating expected distance
 395                 expectedDistance = decalc->calcExpected(Qav, seqCoef);
 396
 397                 if (m->control_pressed) {  return 0; }
 398
 399                 //finding de
 400                 DE = decalc->calcDE(obsDistance, expectedDistance);
 401
 402                 if (m->control_pressed) {  return 0; }
 403
 404                 //find distance between query and closest match
 405                 it = trimmed.begin();
 406                 deviation = decalc->calcDist(query, bestfit, it->first, it->second);
 407
 408                 delete bestfit;
 409
 410                 return 0;
 411         }
 412         catch(exception& e) {
 413                 m->errorOut(e, "Pintail", "getChimeras");
 414                 exit(1);
 415         }
 416 }
 417
 418 //***************************************************************************************************************
 419
 420 vector<float> Pintail::readFreq() {
 421         try {
 422                 //read in probabilities and store in vector
 423                 int pos; float num;
 424
 425                 vector<float> prob;
 426                 set<int> h = decalc->getPos();  //positions of bases in masking sequence
 427
 428         #ifdef USE_MPI
 429
 430                 MPI_File inMPI;
 431                 MPI_Offset size;
 432                 MPI_Status status;
 433
 434                 //char* inFileName = new char[consfile.length()];
 435                 //memcpy(inFileName, consfile.c_str(), consfile.length());
 436
 437                 char inFileName[1024];
 438                 strcpy(inFileName, consfile.c_str());
 439
 440                 MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);
 441                 MPI_File_get_size(inMPI, &size);
 442                 //delete inFileName;
 443
 444                 char* buffer = new char[size];
 445                 MPI_File_read(inMPI, buffer, size, MPI_CHAR, &status);
 446
 447                 string tempBuf = buffer;
 448                 delete buffer;
 449
 450                 if (tempBuf.length() > size) { tempBuf = tempBuf.substr(0, size);  }
 451                 istringstream iss (tempBuf,istringstream::in);
 452
 453                 //read version
 454                 string line = m->getline(iss); m->gobble(iss);
 455
 456                 while(!iss.eof()) {
 457                         iss >> pos >> num;
 458
 459                         if (h.count(pos) > 0) {
 460                                 float Pi;
 461                                 Pi =  (num - 0.25) / 0.75;
 462
 463                                 //cannot have probability less than 0.
 464                                 if (Pi < 0) { Pi = 0.0; }
 465
 466                                 //do you want this spot
 467                                 prob.push_back(Pi);
 468                         }
 469
 470                         m->gobble(iss);
 471                 }
 472
 473                 MPI_File_close(&inMPI);
 474
 475         #else
 476
 477                 ifstream in;
 478                 m->openInputFile(consfile, in);
 479
 480                 //read version
 481                 string line = m->getline(in); m->gobble(in);
 482
 483                 while(!in.eof()){
 484
 485                         in >> pos >> num;
 486
 487                         if (h.count(pos) > 0) {
 488                                 float Pi;
 489                                 Pi =  (num - 0.25) / 0.75;
 490
 491                                 //cannot have probability less than 0.
 492                                 if (Pi < 0) { Pi = 0.0; }
 493
 494                                 //do you want this spot
 495                                 prob.push_back(Pi);
 496                         }
 497
 498                         m->gobble(in);
 499                 }
 500                 in.close();
 501
 502         #endif
 503
 504                 return prob;
 505
 506         }
 507         catch(exception& e) {
 508                 m->errorOut(e, "Pintail", "readFreq");
 509                 exit(1);
 510         }
 511 }
 512
 513 //***************************************************************************************************************
 514 //calculate the distances from each query sequence to all sequences in the template to find the closest sequence
 515 Sequence* Pintail::findPairs(Sequence* q) {
 516         try {
 517
 518                 Sequence* seqsMatches;
 519
 520                 seqsMatches = decalc->findClosest(q, templateSeqs);
 521                 return seqsMatches;
 522
 523         }
 524         catch(exception& e) {
 525                 m->errorOut(e, "Pintail", "findPairs");
 526                 exit(1);
 527         }
 528 }
 529 //**************************************************************************************************
 530 void Pintail::createProcessesQuan() {
 531         try {
 532 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 533                 int process = 1;
 534                 vector<int> processIDS;
 535
 536                 //loop through and create all the processes you want
 537                 while (process != processors) {
 538                         int pid = fork();
 539
 540                         if (pid > 0) {
 541                                 processIDS.push_back(pid);
 542                                 process++;
 543                         }else if (pid == 0){
 544
 545                                 quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, templateLines[process]->start, templateLines[process]->end);
 546
 547                                 //write out data to file so parent can read it
 548                                 ofstream out;
 549                                 string s = toString(getpid()) + ".temp";
 550                                 m->openOutputFile(s, out);
 551
 552                                 //output observed distances
 553                                 for (int i = 0; i < quantilesMembers.size(); i++) {
 554                                         out << quantilesMembers[i].size() << '\t';
 555                                         for (int j = 0; j < quantilesMembers[i].size(); j++) {
 556                                                 out << quantilesMembers[i][j] << '\t';
 557                                         }
 558                                         out << endl;
 559                                 }
 560
 561                                 out.close();
 562
 563                                 exit(0);
 564                         }else {
 565                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 566                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 567                                 exit(0);
 568                         }
 569                 }
 570
 571                 //parent does its part
 572                 quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, templateLines[0]->start, templateLines[0]->end);
 573
 574                 //force parent to wait until all the processes are done
 575                 for (int i=0;i<(processors-1);i++) {
 576                         int temp = processIDS[i];
 577                         wait(&temp);
 578                 }
 579
 580                 //get data created by processes
 581                 for (int i=0;i<(processors-1);i++) {
 582                         ifstream in;
 583                         string s = toString(processIDS[i]) + ".temp";
 584                         m->openInputFile(s, in);
 585
 586                         vector< vector<float> > quan;
 587                         quan.resize(100);
 588
 589                         //get quantiles
 590                         for (int h = 0; h < quan.size(); h++) {
 591                                 int num;
 592                                 in >> num;
 593
 594                                 m->gobble(in);
 595
 596                                 vector<float> q;  float w;
 597                                 for (int j = 0; j < num; j++) {
 598                                         in >> w;
 599                                         q.push_back(w);
 600                                 }
 601
 602                                 quan[h] = q;
 603                                 m->gobble(in);
 604                         }
 605
 606
 607                         //save quan in quantiles
 608                         for (int j = 0; j < quan.size(); j++) {
 609                                 //put all values of q[i] into quan[i]
 610                                 for (int l = 0; l < quan[j].size(); l++) {  quantilesMembers[j].push_back(quan[j][l]);   }
 611                                 //quantilesMembers[j].insert(quantilesMembers[j].begin(), quan[j].begin(), quan[j].end());
 612                         }
 613
 614                         in.close();
 615                         m->mothurRemove(s);
 616                 }
 617
 618 #else
 619                 quantilesMembers = decalc->getQuantiles(templateSeqs, windowSizesTemplate, window, probabilityProfile, increment, 0, templateSeqs.size());
 620 #endif
 621         }
 622         catch(exception& e) {
 623                 m->errorOut(e, "Pintail", "createProcessesQuan");
 624                 exit(1);
 625         }
 626 }
 627 //***************************************************************************************************************
 628 vector< vector<float> > Pintail::readQuantiles() {
 629         try {
 630                 int num;
 631                 float ten, twentyfive, fifty, seventyfive, ninetyfive, ninetynine;
 632
 633                 vector< vector<float> > quan;
 634                 vector <float> temp; temp.resize(6, 0);
 635
 636                 //to fill 0
 637                 quan.push_back(temp);
 638
 639         #ifdef USE_MPI
 640
 641                 MPI_File inMPI;
 642                 MPI_Offset size;
 643                 MPI_Status status;
 644
 645                 //char* inFileName = new char[quanfile.length()];
 646                 //memcpy(inFileName, quanfile.c_str(), quanfile.length());
 647
 648                 char inFileName[1024];
 649                 strcpy(inFileName, quanfile.c_str());
 650
 651                 MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);
 652                 MPI_File_get_size(inMPI, &size);
 653                 //delete inFileName;
 654
 655
 656                 char* buffer = new char[size];
 657                 MPI_File_read(inMPI, buffer, size, MPI_CHAR, &status);
 658
 659                 string tempBuf = buffer;
 660                 if (tempBuf.length() > size) { tempBuf = tempBuf.substr(0, size);  }
 661                 istringstream iss (tempBuf,istringstream::in);
 662                 delete buffer;
 663
 664                 //read version
 665                 string line = m->getline(iss); m->gobble(iss);
 666
 667                 while(!iss.eof()) {
 668                         iss >> num >> ten >> twentyfive >> fifty >> seventyfive >> ninetyfive >> ninetynine;
 669
 670                         temp.clear();
 671
 672                         temp.push_back(ten);
 673                         temp.push_back(twentyfive);
 674                         temp.push_back(fifty);
 675                         temp.push_back(seventyfive);
 676                         temp.push_back(ninetyfive);
 677                         temp.push_back(ninetynine);
 678
 679                         quan.push_back(temp);
 680
 681                         m->gobble(iss);
 682                 }
 683
 684                 MPI_File_close(&inMPI);
 685
 686         #else
 687
 688                 ifstream in;
 689                 m->openInputFile(quanfile, in);
 690
 691                 //read version
 692                 string line = m->getline(in); m->gobble(in);
 693
 694                 while(!in.eof()){
 695
 696                         in >> num >> ten >> twentyfive >> fifty >> seventyfive >> ninetyfive >> ninetynine;
 697
 698                         temp.clear();
 699
 700                         temp.push_back(ten);
 701                         temp.push_back(twentyfive);
 702                         temp.push_back(fifty);
 703                         temp.push_back(seventyfive);
 704                         temp.push_back(ninetyfive);
 705                         temp.push_back(ninetynine);
 706
 707                         quan.push_back(temp);
 708
 709                         m->gobble(in);
 710                 }
 711                 in.close();
 712         #endif
 713
 714                 return quan;
 715
 716         }
 717         catch(exception& e) {
 718                 m->errorOut(e, "Pintail", "readQuantiles");
 719                 exit(1);
 720         }
 721 }
 722 //***************************************************************************************************************/
 723
 724 void Pintail::printQuanFile(string file, string outputString) {
 725         try {
 726
 727                 #ifdef USE_MPI
 728
 729                         MPI_File outQuan;
 730                         MPI_Status status;
 731
 732                         int pid;
 733                         MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 734
 735                         int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
 736
 737                         //char* FileName = new char[file.length()];
 738                         //memcpy(FileName, file.c_str(), file.length());
 739
 740                         char FileName[1024];
 741                         strcpy(FileName, file.c_str());
 742
 743                         if (pid == 0) {
 744                                 MPI_File_open(MPI_COMM_SELF, FileName, outMode, MPI_INFO_NULL, &outQuan);  //comm, filename, mode, info, filepointer
 745
 746                                 int length = outputString.length();
 747                                 char* buf = new char[length];
 748                                 memcpy(buf, outputString.c_str(), length);
 749
 750                                 MPI_File_write(outQuan, buf, length, MPI_CHAR, &status);
 751                                 delete buf;
 752
 753                                 MPI_File_close(&outQuan);
 754                         }
 755
 756                         //delete FileName;
 757                 #else
 758                         ofstream outQuan;
 759                         m->openOutputFile(file, outQuan);
 760
 761                         outQuan << outputString;
 762
 763                         outQuan.close();
 764                 #endif
 765         }
 766         catch(exception& e) {
 767                 m->errorOut(e, "Pintail", "printQuanFile");
 768                 exit(1);
 769         }
 770 }
 771
 772 //***************************************************************************************************************/
 773
 774
 775