1 // $Id: sumfilehandler.cpp,v 1.18 2011/03/07 06:08:47 bobgian Exp $
2
3 /*
4 Copyright 2002 Peter Beerli, Mary Kuhner, Jon Yamato and Joseph Felsenstein
5
6 This software is distributed free of charge for non-commercial use
7 and is copyrighted. Of course, we do not guarantee that the software
8 works, and are not responsible for any damage you may cause or have.
9 */
10
11 #include "sumfilehandler.h"
12 #include "forceparam.h"
13 #include "chainout.h"
14 #include "chainpack.h"
15 #include "collmanager.h"
16 #include "registry.h"
17 #include "region.h"
18 #include "runreport.h"
19 #include "stringx.h"
20 #include "timex.h"
21 #include "treesum.h"
22 #include "xmlsum_strings.h" // for xml sumfile handling
23
24 using namespace std;
25
SumFileHandler()26 SumFileHandler::SumFileHandler()
27 : m_sumin(), m_sumout(), m_lastRegion(0), m_lastReplicate(0),
28 m_lastChain(-1),m_lastRegionChainSum(-1), m_lastReplicateChainSum(-1),
29 m_lastReplicateSummary(0), m_regionSummary(false)
30 {
31 }
32
33 /****************************************************************************
34 *
35 * Sumfile Reading functions
36 *
37 ****************************************************************************/
38
39 // should be only reading fxn to read EOF, all others should not
ReadInSumFile(ChainPack & chainpack,CollectionManager & collectionmanager,long int numchains)40 void SumFileHandler::ReadInSumFile(ChainPack& chainpack,
41 CollectionManager& collectionmanager,
42 long int numchains)
43 {
44 string infilename = registry.GetUserParameters().GetTreeSumInFileName();
45 m_sumin.open(infilename.c_str(), ios::in );
46 if (!m_sumin)
47 {
48 string err_string = "Could not open \"" + infilename + "\". "
49 + "Please check that this file:\n"
50 + " 1) exists in directory that lamarc is being run from,\n"
51 + " 2) is read enabled, and\n"
52 + " 3) is not in use by another program.\n";
53 throw file_error( err_string );
54 }
55
56 string tag;
57 m_sumin >> tag;
58 ReadInCheckFileFormat("ReadInSumFile", xmlsum::SUMFILE_START, tag );
59 m_sumin >> tag;
60
61 while ( (!m_sumin.eof()) && (tag != xmlsum::SUMFILE_END))
62 {
63 if( tag == xmlsum::COMMENT_START )
64 {
65 SkipComments();
66 }
67 else if( tag == xmlsum::CHAINPACK_START )
68 {
69 ReadInChainPack(chainpack);
70 }
71 else if( tag == xmlsum::CHAINSUM_START )
72 {
73 ReadInChainSum(chainpack, collectionmanager, numchains);
74 }
75 else if( tag == xmlsum::REPLICATE_SUMMARY_START )
76 {
77 ReadInReplicateSummary(chainpack);
78 }
79 else if( tag == xmlsum::END_REGION_START )
80 {
81 ReadInEndRegion(chainpack); //This no longer does anything.
82 }
83 else if( tag == xmlsum::REGION_SUMMARY_START )
84 {
85 ReadInRegionSummary(chainpack);
86 }
87 else {//The file is formatted incorrectly.
88 ReadInCheckFileFormat("ReadInSumFile", "a top-level sumfile tag", tag );
89 }
90
91 m_sumin >> tag;
92 }
93 m_sumin.close();
94 } // SumFileHandler::ReadInSumFile
95
96 //------------------------------------------------------------------------------------
97 //------------------------------------------------------------------------------------
98
99 // does not handle nested comments, but should eventually
100 // precondition: last tag read in from stream in: xmlsum::COMMENT_START
101 // postcondition: last string read in from stream in: xmlsum::COMMENT_END
SkipComments()102 void SumFileHandler::SkipComments()
103 {
104 string comment;
105 m_sumin >> comment;
106 while ( !m_sumin.eof() && comment != xmlsum::COMMENT_END )
107 {
108 m_sumin >> comment;
109 if (comment == xmlsum::COMMENT_START) //for nested comments.
110 {
111 SkipComments();
112 }
113 }
114 }
115
116 //------------------------------------------------------------------------------------
117 //------------------------------------------------------------------------------------
118
119 // populate ChainManager's chainpack with a chain summary.
120 // precondition: last tag read in from stream in: xmlsum::CHAINPACK_START
121 // postcondition: last tag read in from stream in: xmlsum::CHAINPACK_END
ReadInChainPack(ChainPack & chainpack)122 void SumFileHandler::ReadInChainPack(ChainPack& chainpack)
123 {
124 string tag;
125 m_sumin >> tag;
126 ReadInCheckFileFormat("ReadInChainPack", xmlsum::NUMBER_START, tag );
127 m_sumin >> m_lastRegion;
128 m_sumin >> m_lastReplicate;
129 m_sumin >> m_lastChain;
130 m_sumin >> tag;
131 ReadInCheckFileFormat("ReadInChainPack", xmlsum::NUMBER_END, tag );
132 m_sumin >> tag;
133 ReadInCheckFileFormat("ReadInChainPack", xmlsum::CHAINOUT_START, tag );
134
135 ChainOut co;
136 ReadInChainOut(co);
137 chainpack.SetChain(co, m_lastRegion, m_lastReplicate, m_lastChain);
138
139 m_sumin >> tag;
140 while (tag==xmlsum::ALPHA_START1)
141 {
142 m_sumin >> tag;
143 ReadInCheckFileFormat( "ReadInChainPack", xmlsum::ALPHA_START2, tag );
144 long int loc;
145 m_sumin >> loc;
146 m_sumin >> tag;
147 ReadInCheckFileFormat( "ReadInChainPack", xmlsum::ALPHA_START3, tag );
148 double alpha;
149 m_sumin >> alpha;
150 registry.GetDataPack().GetRegion(m_lastRegion).GetLocus(loc).GetDataModel()->
151 SetAlpha(alpha, m_lastReplicate, m_lastChain+1);
152 m_sumin >> tag;
153 ReadInCheckFileFormat( "ReadInChainPack", xmlsum::ALPHA_END, tag );
154 m_sumin >> tag;
155 }
156 ReadInCheckFileFormat( "ReadInChainPack", xmlsum::CHAINPACK_END, tag );
157 }
158
159 //------------------------------------------------------------------------------------
160 //------------------------------------------------------------------------------------
161
ReadInReplicateSummary(ChainPack & chainpack)162 void SumFileHandler::ReadInReplicateSummary(ChainPack& chainpack)
163 {
164 m_lastReplicateSummary++;
165 string tag;
166 ForceParameters fp(global_region);
167
168 m_sumin>>tag;
169 ReadInCheckFileFormat("ReadInReplicateSummary", xmlsum::ESTIMATES_START, tag );
170 ReadInForceParameters(fp);
171
172 double maxlike;
173 m_sumin >> tag;
174 ReadInCheckFileFormat("ReadInReplicateSummary", xmlsum::MAXLIKE_START, tag );
175 m_sumin >> maxlike;
176 m_sumin >> tag;
177 ReadInCheckFileFormat("ReadInReplicateSummary", xmlsum::MAXLIKE_END, tag );
178 m_sumin >> tag;
179 ReadInCheckFileFormat("ReadInReplicateSummary", xmlsum::REPLICATE_SUMMARY_END, tag );
180
181 ChainOut repout;
182 repout.SetEstimates(fp);
183 repout.SetLlikemle(maxlike);
184 chainpack.SetSummaryOverReps(repout);
185
186 } // ReadInReplicateSummary
187
ReadInEndRegion(ChainPack & chainpack)188 void SumFileHandler::ReadInEndRegion(ChainPack& chainpack)
189 {
190 //This tag is only present when there is no Replicate Summary information.
191 string tag;
192 m_sumin >> tag;
193 ReadInCheckFileFormat("ReadInSumFile", xmlsum::END_REGION_END, tag );
194 }
195
ReadInRegionSummary(ChainPack & chainpack)196 void SumFileHandler::ReadInRegionSummary(ChainPack& chainpack)
197 {
198 m_regionSummary = true;
199 string tag;
200 ForceParameters fp(global_region);
201
202 m_sumin >> tag;
203 ReadInCheckFileFormat("ReadInRegionSummary", xmlsum::ESTIMATES_START, tag );
204 ReadInForceParameters(fp);
205
206 double maxlike;
207
208 m_sumin >> tag;
209 ReadInCheckFileFormat("ReadInRegionSummary", xmlsum::MAXLIKE_START, tag );
210 m_sumin >> maxlike;
211 m_sumin >> tag;
212 ReadInCheckFileFormat("ReadInRegionSummary", xmlsum::MAXLIKE_END, tag );
213 m_sumin >> tag;
214 ReadInCheckFileFormat("ReadInRegionSummary", xmlsum::REGION_SUMMARY_END, tag );
215
216 ForceSummary& forcesum = registry.GetForceSummary();
217 ChainOut regionout;
218
219 regionout.SetEstimates(fp);
220 regionout.SetLlikemle(maxlike);
221 chainpack.SetSummaryOverRegions(regionout);
222 forcesum.SetOverallMLE(regionout);
223 } // ReadInRegionSummary
224
225 // populate chainout object c with data from the sumfile
226 // precondition: last string read in from stream in: xmlsum::CHAINOUT_START
227 // postcondition: last string read in from stream in: xmlsum::CHAINOUT_END
ReadInChainOut(ChainOut & c)228 void SumFileHandler::ReadInChainOut(ChainOut &c)
229 {
230 string tag, srate;
231 long int p1, p2;
232 long int badtrees, stretchedtrees, tinytrees, zerodltrees;
233 double accrate, llikemle, llikedata;
234 time_t starttime, endtime;
235 ForceParameters fp(global_region);
236
237 do{
238 m_sumin >> tag;
239 if ( tag == xmlsum::BADTREES_START )
240 {
241 m_sumin >> badtrees;
242 c.SetNumBadTrees( badtrees );
243 m_sumin >> tag;
244 ReadInCheckFileFormat( "ReadInChainOut", xmlsum::BADTREES_END, tag );
245 }
246 if ( tag == xmlsum::STRETCHEDTREES_START )
247 {
248 m_sumin >> stretchedtrees;
249 c.SetNumStretchedTrees( stretchedtrees );
250 m_sumin >> tag;
251 ReadInCheckFileFormat( "ReadInChainOut", xmlsum::STRETCHEDTREES_END, tag );
252 }
253 if ( tag == xmlsum::TINYTREES_START )
254 {
255 m_sumin >> tinytrees;
256 c.SetNumTinyPopTrees( tinytrees );
257 m_sumin >> tag;
258 ReadInCheckFileFormat( "ReadInChainOut", xmlsum::TINYTREES_END, tag );
259 }
260 if ( tag == xmlsum::ZERODLTREES_START )
261 {
262 m_sumin >> zerodltrees;
263 c.SetNumZeroDLTrees( zerodltrees );
264 m_sumin >> tag;
265 ReadInCheckFileFormat( "ReadInChainOut", xmlsum::ZERODLTREES_END, tag );
266 }
267 if ( tag == xmlsum::ACCRATE_START )
268 {
269 m_sumin >> accrate;
270 c.SetAccrate( accrate );
271 m_sumin >> tag;
272 ReadInCheckFileFormat( "ReadInChainOut", xmlsum::ACCRATE_END, tag );
273 }
274 if ( tag == xmlsum::LLIKEMLE_START )
275 {
276 m_sumin >> llikemle;
277 c.SetLlikemle( llikemle );
278 m_sumin >> tag;
279 ReadInCheckFileFormat( "ReadInChainOut", xmlsum::LLIKEMLE_END, tag );
280 }
281 if ( tag == xmlsum::LLIKEDATA_START )
282 {
283 m_sumin >> llikedata;
284 c.SetLlikedata( llikedata );
285 m_sumin >> tag;
286 ReadInCheckFileFormat( "ReadInChainOut", xmlsum::LLIKEDATA_END, tag );
287 }
288 if ( tag == xmlsum::STARTTIME_START )
289 {
290 m_sumin >> starttime;
291 c.SetStarttime( starttime );
292 m_sumin >> tag;
293 ReadInCheckFileFormat( "ReadInChainOut", xmlsum::STARTTIME_END, tag );
294 }
295 if ( tag == xmlsum::ENDTIME_START )
296 {
297 m_sumin >> endtime;
298 c.SetEndtime( endtime );
299 m_sumin >> tag;
300 ReadInCheckFileFormat( "ReadInChainOut", xmlsum::ENDTIME_END, tag );
301 }
302 if ( tag == xmlsum::RATES_START )
303 {
304 m_sumin >> tag;
305 ratemap r;
306 while ( tag == xmlsum::MAP_START )
307 {
308 m_sumin >> srate;
309 m_sumin >> p1;
310 m_sumin >> p2;
311 m_sumin >> tag;
312 pair<long int, long int> rpair(p1, p2);
313 r.insert( make_pair(srate, rpair) );
314 c.SetRates( r );
315 ReadInCheckFileFormat( "ReadInChainOut", xmlsum::MAP_END, tag );
316 m_sumin >> tag;
317 }
318 ReadInCheckFileFormat( "ReadInChainOut", xmlsum::RATES_END, tag );
319 }
320 if ( tag == xmlsum::ESTIMATES_START )
321 {
322 ReadInForceParameters(fp);
323 c.SetEstimates(fp);
324 }
325 if ( tag == xmlsum::TEMPERATURES_START)
326 {
327 DoubleVec1d temperatures;
328 if ( ReadInVec1D( temperatures, xmlsum::TEMPERATURES_END ) )
329 {
330 c.SetTemperatures( temperatures );
331 c.SetNumtemps(temperatures.size());
332 }
333 }
334 if ( tag == xmlsum::SWAPRATES_START)
335 {
336 DoubleVec1d swaprates;
337 if ( ReadInVec1D( swaprates, xmlsum::SWAPRATES_END ) )
338 c.SetSwaprates( swaprates );
339 }
340 if ( tag == xmlsum::BAYESUNIQUE_START)
341 {
342 LongVec1d bayesunique;
343 if ( ReadInVec1D( bayesunique, xmlsum::BAYESUNIQUE_END ) )
344 c.SetBayesUnique( bayesunique );
345 }
346
347 } while (tag != xmlsum::CHAINOUT_END);
348 } // ReadInChainOut
349
350 //------------------------------------------------------------------------------------
351 //------------------------------------------------------------------------------------
352
353 // populate forceparameter object fp with data from the sumfile
354 // precondition: last string read in: xmlsum::ESTIMATES_START
355 // postcondition: last string read in: xmlsum::ESTIMATES_END
ReadInForceParameters(ForceParameters & fp)356 void SumFileHandler::ReadInForceParameters( ForceParameters& fp )
357 {
358 vector<double> vd;
359 string tag;
360
361 m_sumin >> tag;
362 while ( tag != xmlsum::ESTIMATES_END )
363 {
364 if ( tag == xmlsum::THETAS_START )
365 {
366 if ( ReadInVec1D( vd, xmlsum::THETAS_END ) )
367 fp.SetGlobalThetas( vd );
368 }
369 if ( tag == xmlsum::MIGRATES_START )
370 {
371 if ( ReadInVec1D( vd, xmlsum::MIGRATES_END ) )
372 fp.SetMigRates( vd );
373 }
374 if ( tag == xmlsum::RECRATES_START )
375 {
376 if ( ReadInVec1D( vd, xmlsum::RECRATES_END ) )
377 fp.SetRecRates( vd );
378 }
379 if ( tag == xmlsum::GROWTHRATES_START )
380 {
381 if ( ReadInVec1D( vd, xmlsum::GROWTHRATES_END ) )
382 fp.SetGrowthRates( vd );
383 }
384 if ( tag == xmlsum::LOGISTICSELECTION_START )
385 {
386 if ( ReadInVec1D( vd, xmlsum::LOGISTICSELECTION_END ) )
387 fp.SetLogisticSelectionCoefficient( vd );
388 }
389 if ( tag == xmlsum::DISEASERATES_START )
390 {
391 if ( ReadInVec1D( vd, xmlsum::DISEASERATES_END ) )
392 fp.SetDiseaseRates( vd );
393 }
394 if ( tag == xmlsum::GAMMAOVERREGIONS_START )
395 {
396 if ( ReadInVec1D( vd, xmlsum::GAMMAOVERREGIONS_END ) )
397 {
398 if ( !registry.GetRegionGammaInfo() )
399 {
400 string msg = "Warning! You are reading summary-file data ";
401 msg += "from a run in which you allowed the mutation rate ";
402 msg += "to vary over genomic regions according to a gamma ";
403 msg += "distribution, and in which the scaled shape parameter ";
404 msg += "of this distribution was estimated to be ";
405 msg += ToString(vd[0]);
406 msg += ". In the present analysis, you have the gamma ";
407 msg += "distribution turned OFF, which will produce different ";
408 msg += "results. You may wish to re-start LAMARC, and ";
409 msg += "turn the gamma distribution on, using the evolutionary ";
410 msg += "forces menu.";
411 registry.GetRunReport().ReportUrgent(msg);
412 }
413 }
414 }
415 m_sumin >> tag;
416 }
417 registry.GetForceSummary().ValidateForceParamOrBarf(fp);
418 ReadInCheckFileFormat( "SumFileHandler::ReadInForceParameters", xmlsum::ESTIMATES_END, tag);
419 } // ReadInForceParameters
420
421 //------------------------------------------------------------------------------------
422 //------------------------------------------------------------------------------------
423
424 // populate vector vd with data from sumfile
425 // clears vd before adding doubles to it
426 // if vd not empty, return true, else return false
ReadInVec1D(vector<double> & vd,string endtag)427 bool SumFileHandler::ReadInVec1D( vector<double>& vd, string endtag )
428 {
429 string sdouble;
430 m_sumin >> sdouble;
431 vd.clear();
432 while ( sdouble != endtag )
433 {
434 vd.push_back( atof(sdouble.c_str()) );
435 m_sumin >> sdouble;
436 }
437
438 if (vd.size() > 0)
439 return true;
440 else
441 return false;
442 } // ReadInVec1D
443
ReadInVec1D(vector<long int> & vd,string endtag)444 bool SumFileHandler::ReadInVec1D( vector<long int> & vd, string endtag )
445 {
446 string slong;
447 m_sumin >> slong;
448 vd.clear();
449 while ( slong != endtag )
450 {
451 vd.push_back( atol(slong.c_str()) );
452 m_sumin >> slong;
453 }
454
455 if (vd.size() > 0)
456 return true;
457 else
458 return false;
459 } // ReadInVec1D
460
461 //------------------------------------------------------------------------------------
462
ReadInChainSum(ChainPack & chainpack,CollectionManager & collectionmanager,long int numchains)463 void SumFileHandler::ReadInChainSum(ChainPack& chainpack, CollectionManager& collectionmanager, long int numchains)
464 {
465 long int region, replicate;
466 string tag;
467 m_sumin >> tag;
468 ReadInCheckFileFormat("ReadInChainSum", xmlsum::REG_REP_START, tag );
469 m_sumin >> region;
470 m_sumin >> replicate;
471 m_sumin >> tag;
472 ReadInCheckFileFormat("ReadInChainSum", xmlsum::REG_REP_END, tag );
473 m_sumin >> tag;
474
475 m_lastRegionChainSum = region;
476 m_lastReplicateChainSum = replicate;
477
478 collectionmanager.StartChain(region, replicate, true);
479 while (((tag == xmlsum::TREESUM_START) ||
480 (tag == xmlsum::PARAM_SUMMARY_START))
481 && (!m_sumin.eof()))
482 {
483 if (tag==xmlsum::TREESUM_START)
484 {
485 TreeSummary* ts = registry.GetProtoTreeSummary().Clone();
486 ts->ReadInTreeSummary(m_sumin);
487 collectionmanager.GetTreeColl(region, replicate)->AddTreeSummary(ts);
488 }
489 else if (tag==xmlsum::PARAM_SUMMARY_START)
490 {
491 m_sumin >> tag;
492 ReadInCheckFileFormat("ReadInChainSum", xmlsum::NCOPY_START, tag );
493 long int ncopy;
494 m_sumin >> ncopy;
495 m_sumin >> tag;
496 ReadInCheckFileFormat("ReadInChainSum", xmlsum::NCOPY_END, tag );
497 m_sumin >> tag;
498 ReadInCheckFileFormat("ReadInChainSum", xmlsum::ESTIMATES_START, tag );
499 ForceParameters fp(region);
500 ReadInForceParameters(fp);
501 collectionmanager.GetParamColl(region, replicate)->AddParamSummary(fp, ncopy);
502 m_sumin >> tag;
503 ReadInCheckFileFormat("ReadInChainSum", xmlsum::PARAM_SUMMARY_END, tag );
504
505 }
506 else
507 {
508 string eitheror = "either " + xmlsum::TREESUM_START + " or "
509 + xmlsum::PARAM_SUMMARY_START;
510 ReadInCheckFileFormat("ReadInChainSum", eitheror, tag);
511 }
512 m_sumin >> tag;
513 }
514
515 if (collectionmanager.GetSampleTrees())
516 {
517 ForceParameters fp(unknown_region);
518 if (numchains >=2 )
519 {
520 fp = chainpack.GetChain(region, replicate, numchains-2).GetEstimates();
521 //-2 is -1 for size->index, and -1 for next-to-last.
522 }
523 else if (numchains == 1)
524 {
525 ForceParameters fpstart(registry.GetForceSummary().GetStartParameters(), region);
526 fp = fpstart;
527 }
528 else
529 {
530 //numchains is negative or zero?
531 throw implementation_error
532 ("The number of chains seems to be zero or negative. This should be impossible.");
533 }
534 collectionmanager.GetTreeColl(region, replicate)->SetStartParameters(fp);
535
536 //The old format had estimates here, but now they're optional.
537 if (tag == xmlsum::ESTIMATES_START)
538 {
539 ForceParameters fp(region);
540 ReadInForceParameters(fp);
541 collectionmanager.GetTreeColl(region, replicate)->SetStartParameters(fp);
542 m_sumin >> tag;
543 }
544 }
545
546 ReadInCheckFileFormat("ReadInChainSum", xmlsum::CHAINSUM_END, tag );
547 //We don't need AdjustSummaries since it was called if needed originally.
548
549 } // ReadInChainSum
550
551 // ReadInCheckFileFormat: checks the file format one tag at a time
552 // args:
553 // callingfxn - should include class calling fxn belongs to
554 // expectedtag - expected xml tag in sumfile
555 // readtag - xml tag just read in by the calling fxn
ReadInCheckFileFormat(string callingfxn,string expectedtag,string readtag)556 void SumFileHandler::ReadInCheckFileFormat(string callingfxn, string expectedtag, string readtag)
557 {
558 if (readtag != expectedtag)
559 {
560 string sumfile = registry.GetUserParameters().GetTreeSumInFileName();
561 string err_string = callingfxn + " - " + sumfile + " has syntactic error.";
562 err_string += " Expected " + expectedtag + " but read " + readtag + "\n";
563 throw incorrect_xml(err_string);
564 }
565 } // ReadInCheckFileFormat
566
567 /****************************************************************************
568 *
569 * Sumfile Writing functions
570 *
571 ****************************************************************************/
572
WriteSumFileStart()573 void SumFileHandler::WriteSumFileStart()
574 {
575 // open file
576 m_sumout.open(registry.GetUserParameters().GetTreeSumOutFileName().c_str(), ios::out );
577 if (!m_sumout)
578 {
579 HandleSumOutFileFormatError("SumFileHandler::WriteSumFileStart");
580 }
581
582 // write preamble
583 m_sumout << xmlsum::SUMFILE_START << endl;
584 m_sumout << xmlsum::COMMENT_START << " Lamarc v. "
585 << VERSION << endl
586 << " Please do not modify. " << xmlsum::COMMENT_END << endl;
587
588 // set the precision of the data (should be greater than that of outfile)
589 m_sumout.precision(SUMFILE_PRECISION);
590 }
591
WriteSumFileEnd(const ChainPack & chainpack)592 void SumFileHandler::WriteSumFileEnd(const ChainPack& chainpack)
593 {
594 if ( m_sumout.is_open() )
595 {
596 // finish up
597 m_sumout << xmlsum::COMMENT_START << " End summary file" << endl
598 << "\t Generated from run that started at: "
599 << PrintTime(chainpack.GetStartTime(), "%c") << endl
600 << "\t and ended at: " << PrintTime(chainpack.GetEndTime(), "%c")
601 << " " << xmlsum::COMMENT_END << endl;
602 m_sumout << xmlsum::SUMFILE_END << endl;
603 m_sumout.close();
604 }
605 else
606 HandleSumOutFileFormatError("WriteSumFileEnd");
607 } // WriteSumFileStart
608
WriteLastChain(const ChainPack & chainpack)609 void SumFileHandler::WriteLastChain(const ChainPack& chainpack)
610 {
611 chainpack.WriteLastChain(m_sumout);
612 }
613
614 //------------------------------------------------------------------------------------
615 //------------------------------------------------------------------------------------
616
617 // currently called after tree generation/processing done
618 // currently outputs necessary data to skip tree gen to calculate MLE's
619
WriteChainSumStart(long int whichregion,long int whichreplicate,const CollectionManager & collectionmanager)620 void SumFileHandler::WriteChainSumStart( long int whichregion, long int whichreplicate,
621 const CollectionManager& collectionmanager )
622 {
623 if ( m_sumout.is_open() )
624 {
625 m_sumout << xmlsum::CHAINSUM_START << endl
626 << "\t" << xmlsum::REG_REP_START << " " << whichregion << " "
627 << whichreplicate << " " << xmlsum::REG_REP_END << endl;
628
629 collectionmanager.WriteThisChainsCollections(&m_sumout);
630 }
631 else
632 HandleSumOutFileFormatError("WriteChainSumStart");
633 } // WriteChainSumsStart
634
WriteChainSumEnd(const CollectionManager & collectionmanager)635 void SumFileHandler::WriteChainSumEnd(const CollectionManager& collectionmanager )
636 {
637 if (m_sumout.is_open())
638 {
639 collectionmanager.WriteLastSummaries();
640 m_sumout << xmlsum::CHAINSUM_END << endl;
641 }
642 else
643 HandleSumOutFileFormatError("WriteChainSumEnd");
644 }
645
646 //Write summary info if more than one region is summarized.
WriteRegionSummary(ForceParameters & fp,double maxlike)647 void SumFileHandler::WriteRegionSummary(ForceParameters& fp, double maxlike)
648 {
649 if ( m_sumout.is_open() )
650 {
651 m_sumout << xmlsum::REGION_SUMMARY_START << endl;
652 fp.WriteForceParameters(m_sumout, 1);
653 m_sumout << "\t" << xmlsum::MAXLIKE_START
654 << " " << maxlike << " "
655 << xmlsum::MAXLIKE_END << endl;
656 m_sumout << xmlsum::REGION_SUMMARY_END << endl;
657 }
658 else
659 HandleSumOutFileFormatError("WriteRegionSummary");
660 } // WriteRegionSummary
661
662 //Write summary info if more than one replicate is summarized.
WriteReplicateSummary(ForceParameters & fp,double maxlike,const ChainPack & chainpack)663 void SumFileHandler::WriteReplicateSummary( ForceParameters& fp, double maxlike, const ChainPack& chainpack)
664 {
665 if ( m_sumout.is_open() )
666 {
667 m_sumout << xmlsum::REPLICATE_SUMMARY_START << endl;
668 fp.WriteForceParameters(m_sumout, 1);
669 m_sumout << "\t" << xmlsum::MAXLIKE_START
670 << " " << maxlike << " "
671 << xmlsum::MAXLIKE_END << endl;
672 m_sumout << xmlsum::REPLICATE_SUMMARY_END << endl;
673 }
674 else
675 HandleSumOutFileFormatError("WriteReplicateSummary");
676 } // WriteReplicateSummary
677
678 //------------------------------------------------------------------------------------
679 //------------------------------------------------------------------------------------
680
681 // purpose: called if sumfile is closed when a summarizing fxn expected it
682 // to be open.
683 // causes: user could have done something to sumfile, or lamarc was
684 // unable to open it in WriteSumFileStart now also called by
685 // non-Summarize fxns, ie TreeSummary::ReadInTreeSummary
HandleSumOutFileFormatError(string callingfxn)686 void SumFileHandler::HandleSumOutFileFormatError(string callingfxn)
687 {
688 if (registry.GetUserParameters().GetWriteSumFile())
689 {
690 string sumfile = registry.GetUserParameters().GetTreeSumOutFileName();
691 string err_string = "In function, " + callingfxn + ", file," + sumfile
692 + " unexpectedly closed or could not be opened and is no longer valid."
693 + " Possible causes of this problem include, but are not limited to: "
694 + "1) read permissions for this file are not or are no longer enabled."
695 + "2) this file was unexpectedly moved, renamed, or deleted. "
696 + "Continuing, but without writing to the summary file any more.\n";
697 registry.GetRunReport().ReportUrgent(err_string);
698 registry.GetUserParameters().SetWriteSumFile(false);
699 }
700 } // HandleSumOutFileFormatError
701
702 // WriteWhatWasRead is called when the user is both reading and writing to a
703 // summary file. This function writes everything at once, instead of doing
704 // it piecemeal within DoChain/DoReplicate/etc. as it would normally.
705 //
706 // This function relies on ReadInRecover working properly and setting
707 // the recover* member variables if the read summary file was not complete.
WriteWhatWasRead(bool recoversumfile,long int recover_region,long int recover_replicate,long int recover_chaintype,long int recover_chain,bool recover_redochain,bool recover_redomaximization,long int nregions,long int nreplicates,const ChainPack & chainpack,const CollectionManager & collectionmanager)708 void SumFileHandler::WriteWhatWasRead(bool recoversumfile,
709 long int recover_region,
710 long int recover_replicate,
711 long int recover_chaintype,
712 long int recover_chain,
713 bool recover_redochain,
714 bool recover_redomaximization,
715 long int nregions,
716 long int nreplicates,
717 const ChainPack& chainpack,
718 const CollectionManager& collectionmanager)
719 {
720 if (m_sumout.is_open())
721 {
722 m_sumout << xmlsum::COMMENT_START
723 << " This summary file should match the input summary file "
724 << registry.GetUserParameters().GetTreeSumInFileName()
725 << ",\n";
726 if (recover_redochain)
727 {
728 m_sumout << " with the exception of the final chainpack, which was re-created. ";
729 }
730 else
731 m_sumout << " up until that file's end. ";
732 m_sumout << xmlsum::COMMENT_END << endl;
733
734 long int total_chains = 0 ;
735 for (int i=0; i<NCHAINTYPES; i++ )
736 {
737 total_chains += registry.GetChainParameters().GetNChains(i);
738 }
739
740 long int last_region = recover_region + 1;
741 for (long int rw_region = 0; rw_region < last_region; rw_region++ )
742 {
743 long int last_replicate = nreplicates;
744 if (rw_region == last_region-1)
745 {
746 //This is the last region, so the number of replicates might be
747 // different than the full number.
748 last_replicate = recover_replicate+1;
749 }
750 for (long int rw_replicate = 0; rw_replicate < last_replicate; rw_replicate++ )
751 {
752 long int last_chain = total_chains;
753 if ((rw_replicate == last_replicate-1) &&
754 (rw_region == last_region-1))
755 {
756 //This is the last replicate, so the number of chains might be
757 // different than the full number.
758 last_chain = 0;
759 for (int i = 0; i < recover_chaintype; i++ )
760 {
761 last_chain += registry.GetChainParameters().GetNChains(i);
762 }
763 last_chain += recover_chain;
764 if (recover_redomaximization)
765 {
766 //We have chain summaries we need to write
767 last_chain++;
768 }
769 }
770 for (long int rw_chain = 0; rw_chain < last_chain; rw_chain++)
771 {
772 bool last_loop = false;
773 if ((rw_region == last_region-1) &&
774 (rw_replicate == last_replicate-1) &&
775 (rw_chain == last_chain-1))
776 {
777 last_loop = true;
778 }
779 if ((rw_chain == total_chains-1) &&
780 (!last_loop || !recover_redochain))
781 {
782 //Write out the chain summary.
783 WriteChainSumStart(rw_region, rw_replicate, collectionmanager);
784 collectionmanager.WriteAllSummaries(rw_region, rw_replicate, m_sumout);
785 m_sumout << xmlsum::CHAINSUM_END << endl;
786 }
787 if (!last_loop || !recover_redomaximization)
788 {
789 //Write out the chainpack.
790 chainpack.WriteChain(m_sumout, rw_region, rw_replicate, rw_chain);
791 }
792 }
793 }
794 if (chainpack.GetLenRegionsVec() - 1 >= rw_region)
795 {
796 //Write the summary over replicates
797 ChainOut rw_chainout = chainpack.GetRegion(rw_region);
798 double rw_maxlike = rw_chainout.GetLlikemle();
799 ForceParameters rw_fp = rw_chainout.GetEstimates();
800 WriteReplicateSummary(rw_fp, rw_maxlike, chainpack);
801 }
802 }
803 if (chainpack.GetLenOverallVec() != 0)
804 {
805 //Write the summary over regions
806 ChainOut rw_chainout = chainpack.GetOverall();
807 double rw_maxlike = rw_chainout.GetLlikemle();
808 ForceParameters rw_fp = rw_chainout.GetEstimates();
809 WriteRegionSummary(rw_fp, rw_maxlike);
810 }
811
812 if (recoversumfile)
813 {
814 m_sumout << xmlsum::COMMENT_START
815 << " New information past this point. "
816 << xmlsum::COMMENT_END << endl;
817 }
818 else
819 WriteSumFileEnd(chainpack);
820 }
821 else
822 HandleSumOutFileFormatError("SumFileHandler::WriteWhatWasRead");
823 } // WriteWhatWasRead
824
WriteVec1D(ofstream & sumout,vector<double> & vd)825 void SumFileHandler::WriteVec1D( ofstream& sumout, vector<double> &vd)
826 {
827 vector<double>::iterator itstart = vd.begin();
828 vector<double>::iterator itend = vd.end();
829 if ( sumout.is_open() )
830 {
831 for( ; itstart != itend; ++itstart )
832 sumout << *itstart << " ";
833 }
834 else
835 SumFileHandler::HandleSumOutFileFormatError("WriteVec1D");
836 } // WriteVec1D
837
WriteVec1D(ofstream & sumout,vector<long int> & vd)838 void SumFileHandler::WriteVec1D( ofstream& sumout, vector<long int> & vd)
839 {
840 vector<long int>::iterator itstart = vd.begin();
841 vector<long int>::iterator itend = vd.end();
842 if ( sumout.is_open() )
843 {
844 for( ; itstart != itend; ++itstart )
845 sumout << *itstart << " ";
846 }
847 else
848 SumFileHandler::HandleSumOutFileFormatError("WriteVec1D");
849 } // WriteVec1D
850
CloseSumOut()851 void SumFileHandler::CloseSumOut()
852 {
853 m_sumout.close();
854 }
855
856 //____________________________________________________________________________________
857