frg.left = start; frg.right = start + len;
assert (mol != fastasizes.begin());  vector<long>::iterator prev = mol - 1; frg.left = start - *prev; frg.right = frg.left + len;
int main(int argc, char ** argv){  if (!GetOptions(argc, argv)){    cerr << "Command line parsing failed, use -h option for usage info" 	 << endl;    //printHelpText();    exit(1);  }  if (globals.find("conffile") == globals.end()){    cerr << "A configuration file must be provided.  see -h" << endl;    exit(1);  }  if (globals.find("input") == globals.end()){    cerr << "An input file must be provided.  see -h" << endl;    exit(1);  }  ifstream conf(globals["conffile"].c_str());  if (! conf.is_open()){    cerr << "Could not open conffile" << globals["conffile"] << endl;    exit(1);  }  parseConfFile(conf);  ifstream in(globals["input"].c_str());  if (! in.is_open()){    cerr << "Coult not open input " << globals["input"] << endl;    exit(1);  }  parseInFile(in);  vector<vector<Fragment_t> > fragments;  fragments.reserve(fastasizes.size());  fragments.resize(fastasizes.size());  // set up the random number generator  gsl_rng * rng = gsl_rng_alloc(gsl_rng_mt19937);  for (int i = 0; i < libraries.size(); i++){    if (libraries[i].coverage == 0)      continue;    long target = (long) libraries[i].coverage * genomeSize;    long totLen = 0;    long mean, stdev;     if (libraries[i].mean == 0){      mean = libraries[i].readmean;      stdev = libraries[i].readstdev;    }    while (totLen < target){      unsigned long start = gsl_rng_uniform_int(rng, genomeSize);      unsigned long len = (long) (mean + gsl_ran_gaussian(rng, stdev));            // find molecule that contains start      vector<long>::iterator mol = 	upper_bound(fastasizes.begin(), fastasizes.end(), start);            if (len < libraries[i].readmin)	continue;               // inserts should be bigger than the				// smallest read      if (mol == fastasizes.end())	if (start + len > genomeSize) // beyond end of genome	  continue;      if (mol != fastasizes.end())	if (start + len > *mol) // beyond end of molecule	  continue;      int idx = mol - fastasizes.begin() - 1; // where to insert fragment      Fragment_t frg;      frg.left = start;      frg.right = start + len;      frg.libid = i;      // now the reads      if (libraries[i].mean == 0){ // unpaired read	frg.forwlen = len;	frg.revlen = 0;	totLen += len;      } else {	int forwlen = (int) (libraries[i].readmean 			     + gsl_ran_gaussian(rng, libraries[i].readstdev));	int revlen = (int) (libraries[i].readmean 			    + gsl_ran_gaussian(rng, libraries[i].readstdev));	// make sure reads are between min and max	while (forwlen > libraries[i].readmax || forwlen < libraries[i].readmin)	  forwlen = (int) (libraries[i].readmean 			     + gsl_ran_gaussian(rng, libraries[i].readstdev));	while (revlen > libraries[i].readmax || revlen < libraries[i].readmin)	  revlen = (int) (libraries[i].readmean 			     + gsl_ran_gaussian(rng, libraries[i].readstdev));	if (forwlen > len) forwlen = len;	if (revlen > len) revlen = len;	frg.forwlen = forwlen;	frg.revlen = revlen;	totLen += forwlen + revlen;      }      fragments[idx].push_back(frg);    } // while coverage not attained  } // for each library  // output libraries    for (int i = 0; i < libraries.size(); i++){    if (libraries[i].coverage == 0)      continue;    if (libraries[i].mean == 0)      continue;      cout << "{LIB\n";    cout << "iid:" << i + 1 << "\n";    cout << "eid:" << libraries[i].name << "\n";    cout << "{DST\n";    cout << "mea:" << libraries[i].mean << "\n";    cout << "std:" << libraries[i].stdev << "\n";    cout << "}\n";    cout << "}\n";  } // for each library again  // now the fragments are all built, we just need to write out the files  in.clear();  // clear flags  in.seekg(0, ios_base::beg); // re-process input file    int seq = 0;  string sequence;  int numseq = 0;  while ((sequence = parseFasta(in)) != ""){    for (int f = 0; f < fragments[seq].size(); f++){      int frgid = 0;      string forw = sequence.substr(fragments[seq][f].left, fragments[seq][f].forwlen);//this one always exists      int doRead = 0;  // do both      double mates = gsl_rng_uniform(rng);      if (fragments[seq][f].revlen == 0 || 	  libraries[fragments[seq][f].libid].permates == 1 ||	  mates < libraries[fragments[seq][f].libid].permates)	doRead = 0;      else	doRead = 1;      mates = gsl_rng_uniform(rng);      if (mates < 0.5)	doRead = -doRead; // do just reverse      if (fragments[seq][f].revlen != 0) { // the whole shebang	string rev = revcompl(sequence.substr(fragments[seq][f].right - fragments[seq][f].revlen, fragments[seq][f].revlen));	frgid = ++numseq;	cout << "{FRG\n";	cout << "iid:" << frgid << "\n";	cout << "lib:" << fragments[seq][f].libid << "\n";	cout << "typ:I\n";	cout << "}\n";	if (doRead <= 0)	  printRED(rev, ++numseq, frgid);      }      if (doRead >= 0)	printRED(forw, ++numseq, frgid);      if (frgid != 0 && doRead == 0){	cout << "{MTP\n";	cout << "rd1:" << numseq - 1 << "\n";	cout << "rd2:" << numseq << "\n";	cout << "typ:E\n";	cout << "}\n";      }    } // for each fragment in this sequence        seq++;  }  return(0);} // main