1 // $Id: forcesummary.cpp,v 1.61 2011/04/23 02:02:49 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 <cassert>
12 #include <cmath>
13
14 #include "local_build.h"
15
16 #include "chainout.h"
17 #include "datapack.h" // for access to Region functions in SummarizeData()
18 #include "force.h" // to create Force objects
19 #include "forceparam.h"
20 #include "forcesummary.h"
21 #include "plotstat.h"
22 #include "region.h"
23 #include "runreport.h"
24 #include "summary.h"
25 #include "tree.h" //to create Tree objects
26 #include "ui_vars.h"
27 #include "vector_constants.h"
28 #include "xml_strings.h" // for ToXML()
29 #include "event.h" // to create Event objects
30 #include "regiongammainfo.h"
31 #include "timemanager.h"
32
33 #ifdef DMALLOC_FUNC_CHECK
34 #include "/usr/local/include/dmalloc.h"
35 #endif
36
37 using namespace std;
38
39 //------------------------------------------------------------------------------------
40
ForceSummary(ForceVec fvec,ForceParameters fparams,DataPack & dpack)41 ForceSummary::ForceSummary(ForceVec fvec, ForceParameters fparams, DataPack& dpack)
42 : m_forcevec(fvec),
43 m_startParameters(fparams),
44 m_llikemles(dpack.GetNRegions()), // likelihoods at the MLEs per region.
45 m_overallLlikemle()
46 {
47 // first set the partition info
48 LongVec1d partitioncounts;
49 long ind = 0;
50 ForceVec::iterator fit;
51 for(fit = m_forcevec.begin(); fit != m_forcevec.end(); ++fit)
52 {
53 if ((*fit)->SetPartIndex(ind))
54 {
55 ++ind;
56 long partCount = (*fit)->GetNPartitions();
57 partitioncounts.push_back(partCount);
58 }
59 //Set up the identical groups;
60 ULongVec2d newgroups = (*fit)->GetIdenticalGroupedParams();
61 for (unsigned long gnum=0; gnum<newgroups.size(); gnum++)
62 {
63 m_identicalGroups.push_back(newgroups[gnum]);
64 }
65 //Set up the multiplicative groups;
66 newgroups = (*fit)->GetMultiplicativeGroupedParams();
67 for (unsigned long gnum=0; gnum<newgroups.size(); gnum++)
68 {
69 m_multiplicativeGroups.push_back(newgroups[gnum]);
70 }
71 //Set up the multiplicative multipliers;
72 DoubleVec1d newmults = (*fit)->GetParamGroupMultipliers();
73 m_multgroupmultiplier.insert(m_multgroupmultiplier.end(),newmults.begin(),newmults.end());
74 }
75 dpack.SetFinalPartitionCounts(partitioncounts);
76
77 } // ForceSummary constructor
78
79 //------------------------------------------------------------------------------------
80
~ForceSummary()81 ForceSummary::~ForceSummary()
82 {
83 ForceVec::iterator fit;
84 for (fit = m_forcevec.begin(); fit != m_forcevec.end(); ++fit)
85 delete(*fit);
86 m_forcevec.clear();
87 } // ForceSummary destructor
88
89 //------------------------------------------------------------------------------------
90
GetAllParameters()91 vector<Parameter> ForceSummary::GetAllParameters()
92 // a private function; public users should create a ParamVector
93 {
94 ForceVec::iterator force;
95 vector<Parameter> result;
96
97 for (force = m_forcevec.begin(); force != m_forcevec.end(); ++force)
98 {
99 vector<Parameter>& tempvec = (*force)->GetParameters();
100 result.insert(result.end(), tempvec.begin(), tempvec.end());
101 }
102
103 return result;
104
105 } // GetAllParameters
106
107 //------------------------------------------------------------------------------------
108
SetAllParameters(const vector<Parameter> & src)109 void ForceSummary::SetAllParameters(const vector<Parameter>& src)
110 // a private function; public users should create a ParamVector
111 {
112 ForceVec::iterator force;
113 vector<Parameter>::const_iterator startparam = src.begin();
114 vector<Parameter>::const_iterator endparam;
115 for (force = m_forcevec.begin(); force != m_forcevec.end(); ++force)
116 {
117 endparam = startparam + (*force)->GetNParams();
118 vector<Parameter> result(startparam, endparam);
119 (*force)->SetParameters(result);
120 startparam = endparam;
121 }
122 //LS NOTE: We might not be at the end of the param vector if there's a gamma, but
123 // we throw that information away. Or keep it somewhere else. Or who knows what.
124
125 } // SetAllParameters
126
127 //------------------------------------------------------------------------------------
128
CreateProtoTree() const129 Tree *ForceSummary::CreateProtoTree() const
130 {
131 Tree *tree;
132
133 if (CheckForce(force_REC))
134 {
135 tree = new RecTree();
136 }
137
138 else
139 {
140 tree = new PlainTree();
141 }
142
143 tree->SetTreeTimeManager(CreateTimeManager());
144
145 return tree;
146 }
147
148 //------------------------------------------------------------------------------------
149
CreateProtoTreeSummary() const150 TreeSummary* ForceSummary::CreateProtoTreeSummary() const
151 {
152 bool shortform = true;
153
154 // if growth is in effect, short form summaries cannot be used
155 if (HasGrowth() || HasLogisticSelection())
156 {
157 shortform = false;
158 }
159
160 // if recombination is in effect, a recombinant summary must be
161 // created
162 TreeSummary* treesum;
163 if (CheckForce(force_REC))
164 {
165 treesum = new RecTreeSummary();
166 }
167 else
168 {
169 treesum = new TreeSummary();
170 }
171
172 IntervalData& interval = treesum->GetIntervalData();
173
174 ForceVec::const_iterator fit = m_forcevec.begin();
175 for ( ; fit != m_forcevec.end(); ++fit)
176 {
177 Summary* sum;
178 if (force_REC == (*fit)->GetTag() && !shortform)
179 {
180 if (CheckForce(force_DISEASE))
181 sum = (*fit)->CreateSummary(interval, false);
182 else
183 sum = (*fit)->CreateSummary(interval, true);
184 }
185 else
186 sum = (*fit)->CreateSummary(interval, shortform);
187 if (sum)
188 {
189 treesum->AddSummary((*fit)->GetTag(), sum);
190 }
191 }
192
193 return treesum;
194
195 } // CreateProtoTreeSummary
196
197 //------------------------------------------------------------------------------------
198
CreateEventVec() const199 vector<Event*> ForceSummary::CreateEventVec() const
200 {
201 vector<Event*> events;
202 ForceVec::const_iterator fit(m_forcevec.begin());
203 for( ; fit != m_forcevec.end(); ++fit)
204 {
205 vector<Event*> ev((*fit)->MakeEvents(*this));
206 events.insert(events.end(),ev.begin(),ev.end());
207 }
208
209 for(fit = m_forcevec.begin(); fit != m_forcevec.end(); ++fit)
210 (*fit)->ModifyEvents(*this,events);
211
212 return events;
213 }
214
215 //------------------------------------------------------------------------------------
216
GetNLocalPartitionForces() const217 long ForceSummary::GetNLocalPartitionForces() const
218 {
219 long nforces(0L);
220 ForceVec::const_iterator fit(m_forcevec.begin());
221 for( ; fit != m_forcevec.end(); ++fit)
222 {
223 if ((*fit)->IsLocalPartitionForce()) ++nforces;
224 }
225 return nforces;
226 }
227
228 //------------------------------------------------------------------------------------
229
GetForceByTag(force_type tag) const230 ForceVec::const_iterator ForceSummary::GetForceByTag(force_type tag) const
231 {
232 // returns an iterator to the element, or end() if none is found
233 ForceVec::const_iterator it;
234 for (it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
235 {
236 if ((*it)->GetTag() == tag) return (it);
237 }
238 it = m_forcevec.end();
239 return (it);
240
241 } // GetForceByTag
242
243 //------------------------------------------------------------------------------------
244
GetStickForce() const245 const StickForce& ForceSummary::GetStickForce() const
246 {
247 ForceVec::const_iterator it;
248 for (it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
249 {
250 if ((*it)->IsStickForce())
251 {
252 StickForce* fp(dynamic_cast<StickForce*>(*it));
253 return *fp;
254 }
255 }
256
257 assert(false); // asked for a stickforce when there wasn't one to be had!!
258
259 return *(dynamic_cast<StickForce*>(*it)); // this is probably illegal,
260 // but we should never be here
261 // see assert above
262
263 } // GetStickForce
264
265 //------------------------------------------------------------------------------------
266
GetLocalPartitionIndexes() const267 LongVec1d ForceSummary::GetLocalPartitionIndexes() const
268 {
269 LongVec1d lpindexes;
270
271 ForceVec::const_iterator it;
272 for (it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
273 {
274 if ((*it)->IsLocalPartitionForce())
275 {
276 lpindexes.push_back(dynamic_cast<PartitionForce*>(*it)->GetPartIndex());
277 }
278 }
279
280 return lpindexes;
281
282 } // GetLocalPartitionIndexes
283
284 //------------------------------------------------------------------------------------
285
GetNonLocalPartitionIndexes() const286 LongVec1d ForceSummary::GetNonLocalPartitionIndexes() const
287 {
288 LongVec1d lpindexes;
289
290 ForceVec::const_iterator it;
291 for (it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
292 {
293 if ((*it)->IsPartitionForce() && (!(*it)->IsLocalPartitionForce()))
294 {
295 lpindexes.push_back(dynamic_cast<PartitionForce*>(*it)->GetPartIndex());
296 }
297 }
298
299 return lpindexes;
300
301 } // GetNonLocalPartitionIndexes
302
303 //------------------------------------------------------------------------------------
304
GetNonLocalPartitionForceTag() const305 force_type ForceSummary::GetNonLocalPartitionForceTag() const
306 {
307 if (CheckForce(force_MIG)) return force_MIG;
308 if (CheckForce(force_DIVMIG)) return force_DIVMIG;
309 assert(false);
310 return force_NONE; // should never happen
311 } // GetNonLocalPartitionForceTag
312
313 //------------------------------------------------------------------------------------
314
IsMissingForce(force_type tag) const315 bool ForceSummary::IsMissingForce(force_type tag) const
316 {
317 ForceVec::const_iterator it = GetForceByTag(tag);
318 return (it == m_forcevec.end());
319 } // IsMissingForce
320
321 //------------------------------------------------------------------------------------
322
AnyForcesOtherThanGrowCoal() const323 bool ForceSummary::AnyForcesOtherThanGrowCoal() const
324 {
325 ForceVec::const_iterator it;
326 for (it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
327 {
328 if ((*it)->GetTag() != force_COAL && (*it)->GetTag() != force_GROW)
329 return true;
330 }
331
332 return false;
333
334 } // AnyForcesOtherThanGrowCoal
335
336 //------------------------------------------------------------------------------------
337
UsingStick() const338 bool ForceSummary::UsingStick() const
339 {
340 ForceVec::const_iterator it;
341 for (it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
342 {
343 if ((*it)->IsStickForce()) return true;
344 }
345
346 return false;
347 }
348
349 //------------------------------------------------------------------------------------
350
CheckForce(force_type tag) const351 bool ForceSummary::CheckForce(force_type tag) const
352 {
353 if (GetForceByTag(tag) == m_forcevec.end()) return false;
354 return true;
355 }
356
357 //------------------------------------------------------------------------------------
358
GetMaxEvents(force_type tag) const359 long ForceSummary::GetMaxEvents(force_type tag) const
360 {
361 ForceVec::const_iterator it = GetForceByTag(tag);
362 assert(it != m_forcevec.end());
363 return((*it)->GetMaxEvents());
364 } // GetMaxEvents
365
366 //------------------------------------------------------------------------------------
367
SetRegionMLEs(const ChainOut & chout,long region)368 void ForceSummary::SetRegionMLEs(const ChainOut& chout, long region)
369 {
370 // we do this by way of the ParamVector linearized form; it's easier
371
372 ParamVector paramvec(false);
373 vector<double> estimates = chout.GetEstimates().GetGlobalParameters();
374
375 assert(estimates.size() == paramvec.size()); // need one estimate per parameter
376
377 ParamVector::iterator param;
378 vector<double>::const_iterator est = estimates.begin();
379
380 for (param = paramvec.begin(); param != paramvec.end(); ++param, ++est)
381 {
382 if (param->IsValidParameter()) param->AddMLE(*est, region);
383 }
384
385 m_llikemles[region] = chout.GetLlikemle();
386
387 } // SetRegionMLEs
388
389 //------------------------------------------------------------------------------------
390
SetOverallMLE(const ChainOut & chout)391 void ForceSummary::SetOverallMLE(const ChainOut& chout)
392 {
393 ParamVector paramvec(false);
394 vector<double> estimates = chout.GetEstimates().GetGlobalParameters();
395
396 assert(estimates.size() == paramvec.size()); // need one estimate per parameter
397
398 ParamVector::iterator param;
399 vector<double>::const_iterator est = estimates.begin();
400
401 for (param = paramvec.begin(); param != paramvec.end(); ++param, ++est)
402 {
403 if (param->IsValidParameter()) param->AddOverallMLE(*est);
404 }
405
406 m_overallLlikemle = chout.GetLlikemle();
407
408 } // SetOverallMLEs
409
410 //------------------------------------------------------------------------------------
411
GetMethods(force_type tag) const412 MethodTypeVec1d ForceSummary::GetMethods(force_type tag) const
413 {
414 ForceVec::const_iterator it = GetForceByTag(tag);
415 assert(it != m_forcevec.end());
416 return (*it)->GetMethods();
417
418 } // GetMethods
419
420 //------------------------------------------------------------------------------------
421
GetEpochs() const422 const vector<Epoch>* ForceSummary::GetEpochs() const
423 {
424 if (!CheckForce(force_DIVERGENCE)) return NULL;
425 ForceVec::const_iterator it = GetForceByTag(force_DIVERGENCE);
426 const DivForce* divforce = dynamic_cast<const DivForce*>(*it);
427 return divforce->GetEpochs();
428 } // GetEpochs
429
430 //------------------------------------------------------------------------------------
431
ToXML(unsigned long nspaces) const432 StringVec1d ForceSummary::ToXML(unsigned long nspaces) const
433 {
434 StringVec1d xmllines;
435
436 string line = MakeIndent(MakeTag(xmlstr::XML_TAG_FORCES),nspaces);
437 xmllines.push_back(line);
438 nspaces += INDENT_DEPTH;
439 ForceVec::const_iterator force;
440 for(force = m_forcevec.begin(); force != m_forcevec.end(); ++force)
441 {
442 StringVec1d forcexml((*force)->ToXML(nspaces));
443 xmllines.insert(xmllines.end(),forcexml.begin(),forcexml.end());
444 }
445 if (registry.GetRegionGammaInfo()) //checks a pointer value
446 {
447 StringVec1d gammaforcexml((*registry.GetRegionGammaInfo()).ToXML(nspaces));
448 xmllines.insert(xmllines.end(),gammaforcexml.begin(),gammaforcexml.end());
449 }
450 nspaces -= INDENT_DEPTH;
451 line = MakeIndent(MakeCloseTag(xmlstr::XML_TAG_FORCES),nspaces);
452 xmllines.push_back(line);
453
454 return xmllines;
455
456 } // ToXML
457
458 //------------------------------------------------------------------------------------
459
ConstrainParameterValues(ForceParameters & fp) const460 bool ForceSummary::ConstrainParameterValues(ForceParameters& fp) const
461 {
462 bool parameters_constrained = false;
463 const ParamVector pvec(true);
464 ForceVec::const_iterator fit;
465 DoubleVec1d params = fp.GetGlobalParameters();
466 bool islog = false;
467 for (unsigned long pnum = 0; pnum< pvec.size(); pnum++)
468 {
469 if (pvec[pnum].IsVariable())
470 {
471 force_type ftype = pvec[pnum].WhichForce();
472 fit = GetForceByTag(ftype);
473 double trimvalue = (*fit)->Truncate(params[pnum]);
474 if (trimvalue != params[pnum])
475 {
476 string msg = "The parameter \"" + pvec[pnum].GetName() +
477 "\" maximized at " + ToString(params[pnum]) +
478 " but is being constrained to " + ToString(trimvalue) + ".\n";
479 RunReport& runreport = registry.GetRunReport();
480 runreport.ReportNormal(msg);
481 SetParamWithConstraints(pnum, trimvalue, params, islog);
482 //params[param] = trimvalue;
483 parameters_constrained = true;
484 }
485 }
486 }
487 fp.SetGlobalParameters(params);
488
489 return parameters_constrained;
490 } // ConstrainParameterValues
491
492 //------------------------------------------------------------------------------------
493
Get2DMethods(force_type tag) const494 MethodTypeVec2d ForceSummary::Get2DMethods(force_type tag) const
495 {
496 ForceVec::const_iterator it = GetForceByTag(tag);
497 assert(it != m_forcevec.end());
498 MethodTypeVec2d result = SquareOffVector((*it)->GetMethods());
499 return result;
500 } // Get2DMethods
501
502 //------------------------------------------------------------------------------------
503
GetPartitionForces() const504 const ForceVec ForceSummary::GetPartitionForces() const
505 {
506 ForceVec partitions;
507 ForceVec::const_iterator it;
508 for(it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
509 if ((*it)->IsPartitionForce()) partitions.push_back(*it);
510
511 return partitions;
512 } // GetPartitionForces
513
514 //------------------------------------------------------------------------------------
515
GetLocalPartitionForces() const516 const ForceVec ForceSummary::GetLocalPartitionForces() const
517 {
518 ForceVec partitions;
519 ForceVec::const_iterator it;
520 for(it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
521 if ((*it)->IsLocalPartitionForce()) partitions.push_back(*it);
522
523 return partitions;
524 } // GetLocalPartitionForces
525
526 //------------------------------------------------------------------------------------
527
GetParamsIdenticalGroupedWith(unsigned long pindex) const528 ULongVec1d ForceSummary::GetParamsIdenticalGroupedWith(unsigned long pindex) const
529 {
530 for (unsigned long gnum=0; gnum<m_identicalGroups.size(); gnum++)
531 {
532 for (unsigned long gpnum=0; gpnum<m_identicalGroups[gnum].size(); gpnum++)
533 {
534 if (m_identicalGroups[gnum][gpnum] == pindex)
535 {
536 return m_identicalGroups[gnum];
537 }
538 }
539 }
540 assert(false); //It's probably not a good idea to check this for *all*
541 // pindices.
542 ULongVec1d emptyvec;
543 emptyvec.push_back(pindex);
544 return emptyvec;
545 }
546
547 //------------------------------------------------------------------------------------
548
GetParamsMultiplicativelyGroupedWith(unsigned long pindex) const549 ULongVec1d ForceSummary::GetParamsMultiplicativelyGroupedWith(unsigned long pindex) const
550 {
551 for (unsigned long gnum=0; gnum<m_multiplicativeGroups.size(); gnum++)
552 {
553 for (unsigned long gpnum=0; gpnum<m_multiplicativeGroups[gnum].size(); gpnum++)
554 {
555 if (m_multiplicativeGroups[gnum][gpnum] == pindex)
556 {
557 return m_multiplicativeGroups[gnum];
558 }
559 }
560 }
561 assert(false); //It's probably not a good idea to check this for *all*
562 // pindices.
563 ULongVec1d emptyvec;
564 emptyvec.push_back(pindex);
565 return emptyvec;
566 }
567
568 //------------------------------------------------------------------------------------
569
GetGroupMultiplier(unsigned long pindex) const570 double ForceSummary::GetGroupMultiplier(unsigned long pindex) const
571 {
572 for (unsigned long gnum=0; gnum<m_multiplicativeGroups.size(); gnum++)
573 {
574 for (unsigned long gpnum=0; gpnum<m_multiplicativeGroups[gnum].size(); gpnum++)
575 {
576 if (m_multiplicativeGroups[gnum][gpnum] == pindex)
577 {
578 return m_multgroupmultiplier[gnum];
579 }
580 }
581 }
582 assert(false); //It's probably not a good idea to check this for *all*
583 // pindices.
584 return FLAGDOUBLE;
585 }
586
587 //------------------------------------------------------------------------------------
588
GetForceTags() const589 vector<force_type> ForceSummary::GetForceTags() const
590 {
591 vector<force_type> retVec;
592 for (unsigned long fnum = 0; fnum<m_forcevec.size(); fnum++)
593 {
594 retVec.push_back(m_forcevec[fnum]->GetTag());
595 }
596 return retVec;
597 }
598
599 //------------------------------------------------------------------------------------
600
GetForceSizes() const601 vector<long> ForceSummary::GetForceSizes() const
602 {
603 vector<long> retVec;
604 for (unsigned long fnum = 0; fnum<m_forcevec.size(); fnum++)
605 {
606 retVec.push_back(m_forcevec[fnum]->GetNParams());
607 }
608 return retVec;
609 }
610
611 //------------------------------------------------------------------------------------
612
SetParamWithConstraints(long pindex,double newval,DoubleVec1d & pvecnums,bool islog) const613 void ForceSummary::SetParamWithConstraints(long pindex, double newval,
614 DoubleVec1d& pvecnums, bool islog) const
615 {
616 //Note that this function is sometimes called with 'pvecnums' equal to
617 // *log* parameters, not always 'normal' parameters.
618
619 const ParamVector pvec(true);
620 // MDEBUG a good refactor would be to unify identical-group and multiplicative-group logic here
621 ParamStatus mystatus = pvec[pindex].GetStatus();
622 ULongVec1d groupmembers;
623 double multiplier = 1.0;
624 if (mystatus.Status() == pstat_identical_head)
625 {
626 groupmembers = GetParamsIdenticalGroupedWith(pindex);
627 }
628 if (mystatus.Status() == pstat_multiplicative_head)
629 {
630 groupmembers = GetParamsMultiplicativelyGroupedWith(pindex);
631 multiplier = GetGroupMultiplier(pindex);
632 if (islog) multiplier = log(multiplier);
633 }
634 mystatus.SetWithConstraints(pindex, newval, pvecnums, groupmembers, multiplier);
635 }
636
637 //------------------------------------------------------------------------------------
638
GetNParameters(force_type tag) const639 long ForceSummary::GetNParameters(force_type tag) const
640 {
641 ForceVec::const_iterator it = GetForceByTag(tag);
642 assert(it != m_forcevec.end());
643 return (*it)->GetNParams();
644 } // GetNParameters
645
646 //------------------------------------------------------------------------------------
647
GetAllNParameters() const648 long ForceSummary::GetAllNParameters() const
649 {
650 long total = 0;
651 ForceVec::const_iterator it;
652 for(it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
653 total += (*it)->GetNParams();
654 return total;
655 } // ForceSummary::GetAllNParameters
656
657 //------------------------------------------------------------------------------------
658
GetModifiers(long posinparamvec) const659 vector<double> ForceSummary::GetModifiers(long posinparamvec) const
660 {
661 const ParamVector paramvec(true); // read-only copy
662 vector<double> results;
663
664 // diagnose whether we have fixed or percentile profiles.
665
666 proftype profiletype = paramvec[posinparamvec].GetProfileType();
667
668 //Note: The profiling type must be consistent within a force, but can
669 // change from force to force. Also within a force, profiling can be turned
670 // on or off.
671
672 if (profiletype == profile_NONE) return results; // no profiles, hence no modifiers
673
674 verbosity_type verbosity = registry.GetUserParameters().GetVerbosity();
675
676 if (profiletype == profile_PERCENTILE) // percentile
677 {
678 if (verbosity == CONCISE || verbosity == NONE)
679 results = vecconst::percentiles_short;
680 else
681 results = vecconst::percentiles;
682 return results;
683 }
684 else
685 { // fixed
686 if (verbosity == CONCISE || verbosity == NONE)
687 {
688 results = vecconst::multipliers_short; //growth or not, just do this one.
689 }
690 else
691 {
692 if(paramvec[posinparamvec].IsForce(force_GROW))
693 {
694 results = vecconst::growthmultipliers;
695 DoubleVec1d gfixed = vecconst::growthfixed;
696 for (unsigned long i=0; i<gfixed.size(); i++)
697 results.push_back(gfixed[i]);
698 //There are two components to this vector--the first the normal
699 // multipliers of growth, but the second actual values for growth
700 // which we want to test. We'll have to catch this later in
701 // Analyzer::CalcProfileFixed (in profile.cpp)
702 }
703 else
704 results = vecconst::multipliers;
705 }
706 return results;
707 }
708 } // GetModifiers
709
710 //------------------------------------------------------------------------------------
711
GetModifiersByForce(force_type tag) const712 DoubleVec1d ForceSummary::GetModifiersByForce(force_type tag) const
713 {
714 const ParamVector paramvec(true); // read-only copy
715 for (unsigned long i=0; i<paramvec.size(); ++i)
716 if (paramvec[i].IsValidParameter() && paramvec[i].IsForce(tag))
717 if (paramvec[i].GetProfileType() != profile_NONE)
718 return GetModifiers(i);
719 return GetModifiers(0);
720 //This assumes that all parameters of a given force have either the
721 // same profiling type or none at all.
722 } // GetModifiersByForce
723
724 //------------------------------------------------------------------------------------
725
GetOverallProfileType() const726 proftype ForceSummary::GetOverallProfileType() const
727 {
728 unsigned long i;
729 for (i = 0; i < m_forcevec.size(); ++i)
730 {
731 proftype mytype = m_forcevec[i]->SummarizeProfTypes();
732 if (mytype != profile_NONE) return mytype;
733 }
734
735 return profile_NONE;
736
737 } // GetOverallProfileType
738
739 //------------------------------------------------------------------------------------
740
GetLowParameter(long posinparamvec) const741 double ForceSummary::GetLowParameter(long posinparamvec) const
742 {
743 const ParamVector paramvec(true); // read-only copy
744 force_type tag = paramvec[posinparamvec].WhichForce();
745 if (force_REGION_GAMMA == tag)
746 {
747 const RegionGammaInfo *pRegionGammaInfo = registry.GetRegionGammaInfo();
748 if (!pRegionGammaInfo)
749 {
750 string msg = "ForceSummary::GetLowParameter() attempted to get the low ";
751 msg += "parameter value for force \"" + ToString(tag);
752 msg += ",\" but the necessary RegionGammaInfo was not found in the registry.";
753 throw implementation_error(msg);
754 }
755 return pRegionGammaInfo->GetLowValue();
756 }
757 ForceVec::const_iterator fit = GetForceByTag(tag);
758 if (m_forcevec.end() == fit)
759 {
760 string msg = "ForceSummary::GetLowParameter() failed to get a force object ";
761 msg += "corresponding to tag \"" + ToString(tag) + ".\"";
762 throw implementation_error(msg);
763 }
764 return (*fit)->GetLowVal();
765 }
766
767 //------------------------------------------------------------------------------------
768
GetHighParameter(long posinparamvec) const769 double ForceSummary::GetHighParameter(long posinparamvec) const
770 {
771 const ParamVector paramvec(true); // read-only copy
772 force_type tag = paramvec[posinparamvec].WhichForce();
773 if (force_REGION_GAMMA == tag)
774 {
775 const RegionGammaInfo *pRegionGammaInfo = registry.GetRegionGammaInfo();
776 if (!pRegionGammaInfo)
777 {
778 string msg = "ForceSummary::GetHighParameter() attempted to get the high ";
779 msg += "parameter value for force \"" + ToString(tag);
780 msg += ",\" but the necessary RegionGammaInfo was not found in the registry.";
781 throw implementation_error(msg);
782 }
783 return pRegionGammaInfo->GetHighValue();
784 }
785 ForceVec::const_iterator fit = GetForceByTag(tag);
786 if (m_forcevec.end() == fit)
787 {
788 string msg = "ForceSummary::GetHighParameter() failed to get a force object ";
789 msg += "corresponding to tag \"" + ToString(tag) + ".\"";
790 throw implementation_error(msg);
791 }
792 return (*fit)->GetHighVal();
793 }
794
795 //------------------------------------------------------------------------------------
796
GetLowMult(long posinparamvec) const797 double ForceSummary::GetLowMult(long posinparamvec) const
798 {
799 const ParamVector paramvec(true); // read-only copy
800 force_type tag = paramvec[posinparamvec].WhichForce();
801 if (force_REGION_GAMMA == tag)
802 {
803 const RegionGammaInfo *pRegionGammaInfo = registry.GetRegionGammaInfo();
804 if (!pRegionGammaInfo)
805 {
806 string msg = "ForceSummary::GetLowMult() attempted to get the low ";
807 msg += "multiplier value for force \"" + ToString(tag);
808 msg += ",\" but the necessary RegionGammaInfo was not found in the registry.";
809 throw implementation_error(msg);
810 }
811 return pRegionGammaInfo->GetLowMultiplier();
812 }
813 ForceVec::const_iterator fit = GetForceByTag(tag);
814 if (m_forcevec.end() == fit)
815 {
816 string msg = "ForceSummary::GetLowMult() failed to get a force object ";
817 msg += "corresponding to tag \"" + ToString(tag) + ".\"";
818 throw implementation_error(msg);
819 }
820 return (*fit)->GetLowMult();
821 }
822
823 //------------------------------------------------------------------------------------
824
GetHighMult(long posinparamvec) const825 double ForceSummary::GetHighMult(long posinparamvec) const
826 {
827 const ParamVector paramvec(true); // read-only copy
828 force_type tag = paramvec[posinparamvec].WhichForce();
829 if (force_REGION_GAMMA == tag)
830 {
831 const RegionGammaInfo *pRegionGammaInfo = registry.GetRegionGammaInfo();
832 if (!pRegionGammaInfo)
833 {
834 string msg = "ForceSummary::GetHighMult() attempted to get the high ";
835 msg += "multiplier value for force \"" + ToString(tag);
836 msg += ",\" but the necessary RegionGammaInfo was not found in the registry.";
837 throw implementation_error(msg);
838 }
839 return pRegionGammaInfo->GetHighMultiplier();
840 }
841 ForceVec::const_iterator fit = GetForceByTag(tag);
842 if (m_forcevec.end() == fit)
843 {
844 string msg = "ForceSummary::GetHighMult() failed to get a force object ";
845 msg += "corresponding to tag \"" + ToString(tag) + ".\"";
846 throw implementation_error(msg);
847 }
848 return (*fit)->GetHighMult();
849 }
850
851 //------------------------------------------------------------------------------------
852
GetPartIndex(force_type forcename) const853 long ForceSummary::GetPartIndex(force_type forcename) const
854 {
855 assert(dynamic_cast<PartitionForce*>(*GetForceByTag(forcename)));
856 return dynamic_cast<PartitionForce*>(*GetForceByTag(forcename))->
857 GetPartIndex();
858 } // ForceSummary::GetPartIndex
859
860 //------------------------------------------------------------------------------------
861
CreateTimeManager() const862 TimeManager* ForceSummary::CreateTimeManager() const
863 {
864 if (UsingStick())
865 {
866 if (HasGrowth())
867 {
868 if (HasSelection())
869 {
870 string msg = "ForceSummary::CreateTimeManager(), detected the ";
871 msg += "simultaneous presence of the growth and selection ";
872 msg += "forces. If LAMARC has been updated to allow both of these ";
873 msg += "forces to act at the same time, then ";
874 msg += "ForceSummary::CreateTimeManager() needs to be updated ";
875 msg += "accordingly.";
876 throw implementation_error(msg);
877 }
878 return new ExpGrowStickTimeManager();
879 }
880
881 if (HasSelection())
882 return new SelectStickTimeManager();
883
884 }
885 else
886 {
887 if (HasGrowth())
888 {
889 if (HasSelection())
890 {
891 string msg = "ForceSummary::CreateTimeManager(), detected the ";
892 msg += "simultaneous presence of the growth and selection ";
893 msg += "forces. If LAMARC has been updated to allow both of these ";
894 msg += "forces to act at the same time, then ";
895 msg += "ForceSummary::CreateTimeManager() needs to be updated ";
896 msg += "accordingly.";
897 throw implementation_error(msg);
898 }
899 return new ExpGrowTimeManager();
900 }
901
902 if (HasLogisticSelection())
903 return new LogSelectTimeManager();
904 }
905
906 return new ConstantTimeManager();
907
908 } // ForceSummary::CreateTimeManager
909
910 //------------------------------------------------------------------------------------
911 // WARNING warning--butt ugly code! maybe add new classes of
912 // forces LocalPartitionForces and SexualForces (REC && CONVERSION)
913
GetSelectedSites() const914 LongVec1d ForceSummary::GetSelectedSites() const
915 {
916 bool foundsex(false);
917 ForceVec::const_iterator force;
918 for(force = m_forcevec.begin(); force != m_forcevec.end(); ++force)
919 {
920 if ((*force)->IsSexualForce())
921 {
922 foundsex = true;
923 break;
924 }
925 }
926
927 LongVec1d selsites;
928
929 if (foundsex)
930 {
931 for(force = m_forcevec.begin(); force != m_forcevec.end(); ++force)
932 {
933 if ((*force)->IsLocalPartitionForce())
934 selsites.push_back(dynamic_cast<LocalPartitionForce*>(*force)->GetLocalSite());
935 }
936 }
937
938 return selsites;
939
940 } // ForceSummary::GetSelectedSites
941
942 //------------------------------------------------------------------------------------
943
IsValidForceSummary() const944 bool ForceSummary::IsValidForceSummary() const
945 {
946 // Some parameters must be present!
947 if (m_forcevec.size() == 0) return false;
948
949 // unsigned long npops = nparams;
950
951 unsigned long ncoal = 0, nmig = 0, ndisease = 0, nrec = 0, ngrow = 0, ngamma = 0;
952 unsigned long nlselect = 0, nepoch = 0;
953
954 if (!CheckForce(force_COAL)) return false;
955 ncoal = m_startParameters.GetGlobalThetas().size();
956
957 if (CheckForce(force_MIG))
958 {
959 nmig = m_startParameters.GetMigRates().size();
960 if (nmig < 2) return false;
961 }
962
963 if (CheckForce(force_DIVMIG))
964 {
965 nmig = m_startParameters.GetMigRates().size();
966 if (nmig < 2) return false;
967 }
968
969 if (CheckForce(force_DISEASE))
970 {
971 ndisease = m_startParameters.GetDiseaseRates().size();
972 if (ndisease < 2) return false;
973 }
974
975 //nmig and ndisease are set to 1 because I later take the square root to
976 // determine the number of partitions for each.
977
978 if (CheckForce(force_DIVERGENCE))
979 {
980 nepoch = m_startParameters.GetEpochTimes().size();
981 unsigned long nmigparts = static_cast<long>(sqrt(static_cast<double>(nmig)));
982 if (nepoch != (nmigparts - 1)/2) return false; // truncation okay
983 }
984
985 if (CheckForce(force_REC))
986 {
987 nrec = m_startParameters.GetRecRates().size();
988 if (nrec != 1) return false;
989 // Until we implement distinct rec rates per population, at least.
990 }
991
992 if (CheckForce(force_GROW))
993 {
994 ngrow = m_startParameters.GetGrowthRates().size();
995 if (ngrow != ncoal) return false;
996 }
997
998 if (CheckForce(force_LOGISTICSELECTION))
999 {
1000 nlselect = m_startParameters.GetLogisticSelectionCoefficient().size();
1001 if (1 != nlselect) return false;
1002 }
1003
1004 if (CheckForce(force_EXPGROWSTICK))
1005 {
1006 ngrow = m_startParameters.GetGrowthRates().size();
1007 if (ngrow != ncoal) return false;
1008 }
1009
1010 unsigned long nmigparts;
1011 if (nmig)
1012 {
1013 nmigparts = static_cast<long>(sqrt(static_cast<double>(nmig)));
1014 if (nmigparts*nmigparts != nmig) return false;
1015 }
1016 else nmigparts = 1;
1017
1018 unsigned long ndiseaseparts;
1019 if (ndisease)
1020 {
1021 ndiseaseparts = static_cast<long>(sqrt(static_cast<double>(ndisease)));
1022 if (ndiseaseparts*ndiseaseparts != ndisease) return false;
1023 }
1024 else ndiseaseparts = 1;
1025
1026 if (ncoal != nmigparts * ndiseaseparts) return false;
1027
1028 //Now check out a parameter vector and check to be sure the number of total
1029 // parameters is accurate, and that each parameter's index value is right.
1030 const ParamVector testvec(true);
1031
1032 if (testvec.size() != ncoal + nmig + ndisease + nrec + nepoch + ngrow + ngamma + nlselect)
1033 return false;
1034
1035 for (unsigned long index=0; index < testvec.size(); index++)
1036 {
1037 if (index != testvec[index].GetParamVecIndex()) return false;
1038 }
1039
1040 // WARNING DEBUG should check validity of the parameterlist
1041 // (but I don't know how) -- Mary
1042
1043 return true;
1044 } // IsValidForceSummary
1045
1046 //------------------------------------------------------------------------------------
1047
ValidateForceParamOrBarf(const ForceParameters & fp)1048 void ForceSummary::ValidateForceParamOrBarf(const ForceParameters& fp)
1049 {
1050 DoubleVec1d params;
1051 vector<force_type> allforces;
1052 allforces.push_back(force_COAL);
1053 allforces.push_back(force_MIG);
1054 allforces.push_back(force_DIVMIG);
1055 allforces.push_back(force_DISEASE);
1056 allforces.push_back(force_REC);
1057 allforces.push_back(force_GROW);
1058 allforces.push_back(force_LOGISTICSELECTION);
1059 allforces.push_back(force_DIVERGENCE);
1060
1061 for (vector<force_type>::iterator force = allforces.begin();
1062 force != allforces.end(); force++)
1063 {
1064 params = fp.GetGlobalParametersByTag((*force));
1065 if (CheckForce((*force)))
1066 {
1067 //This force is active, and should have a non-zero-sized vector in the forceparameters object.
1068 if (params.size() == 0)
1069 {
1070 string msg = "Error: The force " + ToString((*force)) +
1071 " was turned on in this run of LAMARC, but\n" +
1072 " was not on in the LAMARC run that created this summary file.\n\n"
1073 + "Re-run LAMARC with " + ToString((*force)) +
1074 " turned off to accurately continue the old run.";
1075 throw data_error(msg);
1076 }
1077 }
1078 else
1079 {
1080 //This force is inactive, and should have a zero-sized vector in the forceparameters object.
1081 if (params.size() > 0)
1082 {
1083 string msg = "Error: The force " + ToString((*force)) +
1084 " was turned off in this run of LAMARC, but\n" +
1085 " was on in the LAMARC run that created this summary file.\n\n" +
1086 "Re-run LAMARC with " + ToString((*force)) +
1087 " turned on to accurately continue the old run.";
1088 throw data_error(msg);
1089 }
1090 }
1091 }
1092 }
1093
1094 //____________________________________________________________________________________
1095