1 // $Id: forceparam.cpp,v 1.40 2011/03/07 06:08: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
13 #include "local_build.h"
14
15 #include "forceparam.h"
16 #include "forcesummary.h"
17 #include "force.h"
18 #include "mathx.h"
19 #include "region.h"
20 #include "sumfilehandler.h"
21 #include "vectorx.h" // for LongSquareRootOfLong
22 #include "xmlsum_strings.h" // for xml sumfile handling
23
24 #ifdef DMALLOC_FUNC_CHECK
25 #include "/usr/local/include/dmalloc.h"
26 #endif
27
28 using namespace std;
29
30 //------------------------------------------------------------------------------------
31
ForceParameters(long region)32 ForceParameters::ForceParameters(long region)
33 : m_region(region),
34 m_space(known_region),
35 m_forceTags(registry.GetForceSummary().GetForceTags()),
36 m_forceSizes(registry.GetForceSummary().GetForceSizes()),
37 m_epochptr(registry.GetForceSummary().GetEpochs())
38 {
39 assert (region >= 0);
40 }
41
42 //------------------------------------------------------------------------------------
43
ForceParameters(param_space space)44 ForceParameters::ForceParameters(param_space space)
45 : m_region(FLAGLONG),
46 m_space(space),
47 m_forceTags(registry.GetForceSummary().GetForceTags()),
48 m_forceSizes(registry.GetForceSummary().GetForceSizes())
49 {
50 assert (m_space != known_region); //can be either global or unknown
51 }
52
53 //------------------------------------------------------------------------------------
54 //Need this constructor for the registry, when we don't have a forcesummary yet.
55
ForceParameters(param_space space,ForceTypeVec1d types,LongVec1d sizes)56 ForceParameters::ForceParameters(param_space space, ForceTypeVec1d types, LongVec1d sizes)
57 : m_region(FLAGLONG),
58 m_space(space),
59 m_forceTags(types),
60 m_forceSizes(sizes)
61 {
62 assert (m_space != known_region); // Can be either global or unknown.
63 }
64
65 //------------------------------------------------------------------------------------
66
ForceParameters(const ForceParameters & fp,long region)67 ForceParameters::ForceParameters(const ForceParameters &fp, long region)
68 : m_region(region),
69 m_space(known_region),
70 m_forceTags(fp.m_forceTags),
71 m_forceSizes(fp.m_forceSizes)
72 {
73 assert (region >= 0);
74 switch (fp.GetParamSpace())
75 {
76 case unknown_region:
77 CopyRegionalMembers(fp);
78 FillGlobalParamsFromRegionalParams();
79 break;
80 case known_region:
81 if (fp.m_region != region)
82 {
83 assert(false); //Why are we assigning fps for different regions to each other?
84 }
85 //else fall through to:
86 case global_region:
87 CopyGlobalMembers(fp);
88 FillRegionalParamsFromGlobalParams();
89 break;
90 }
91 }
92
93 //------------------------------------------------------------------------------------
94
CopyGlobalMembers(const ForceParameters & fp)95 void ForceParameters::CopyGlobalMembers(const ForceParameters& fp)
96 {
97 m_globalThetas = fp.m_globalThetas;
98 m_migrates = fp.m_migrates;
99 m_recrates = fp.m_recrates;
100 m_growths = fp.m_growths;
101 m_diseases = fp.m_diseases;
102 m_logisticSelectionCoefficient = fp.m_logisticSelectionCoefficient;
103 m_epochtimes = fp.m_epochtimes;
104 } // ForceParameters::CopyGlobalMembers
105
106 //------------------------------------------------------------------------------------
107
CopyRegionalMembers(const ForceParameters & fp)108 void ForceParameters::CopyRegionalMembers(const ForceParameters& fp)
109 {
110 m_regionalThetas = fp.m_regionalThetas;
111 m_migrates = fp.m_migrates;
112 m_recrates = fp.m_recrates;
113 m_growths = fp.m_growths;
114 m_diseases = fp.m_diseases;
115 m_logisticSelectionCoefficient = fp.m_logisticSelectionCoefficient;
116 m_epochtimes = fp.m_epochtimes;
117 } // ForceParameters::CopyRegionalMembers
118
119 //------------------------------------------------------------------------------------
120
SetGlobalParametersByTag(force_type tag,const DoubleVec1d & v)121 void ForceParameters::SetGlobalParametersByTag(force_type tag, const DoubleVec1d& v)
122 {
123 assert (m_space != unknown_region);
124 string msg = "ForceParameters::SetGlobalParametersByTag() received tag ";
125 switch (tag)
126 {
127 case force_COAL:
128 SetGlobalThetas(v);
129 return;
130 case force_MIG:
131 case force_DIVMIG:
132 SetMigRates(v);
133 return;
134 case force_REC:
135 SetRecRates(v);
136 return;
137 case force_EXPGROWSTICK:
138 case force_GROW:
139 SetGrowthRates(v);
140 return;
141 case force_LOGSELECTSTICK:
142 case force_LOGISTICSELECTION:
143 SetLogisticSelectionCoefficient(v);
144 return;
145 case force_DISEASE:
146 SetDiseaseRates(v);
147 return;
148 case force_DIVERGENCE:
149 SetEpochTimes(v);
150 return;
151 case force_REGION_GAMMA:
152 msg += "\"" + ToString(tag) + ",\" which should never happen. (This is a ";
153 msg += "pseudo-force that\'s only treated as a force within certain ";
154 msg += "contexts of the program.)";
155 throw implementation_error(msg);
156 case force_NONE:
157 msg += "\"" + ToString(tag) + ",\" which should never happen.";
158 throw implementation_error(msg);
159 }
160
161 msg += "\"" + ToString(tag) + ",\" but it does not know how to use it.";
162 throw implementation_error(msg);
163 } // SetGlobalParametersByTag
164
165 //------------------------------------------------------------------------------------
166
SetRegionalParametersByTag(force_type tag,const DoubleVec1d & v)167 void ForceParameters::SetRegionalParametersByTag(force_type tag, const DoubleVec1d& v)
168 {
169 assert (m_space != global_region);
170 string msg = "ForceParameters::SetRegionalParametersByTag() received tag ";
171 switch (tag)
172 {
173 case force_COAL:
174 SetRegionalThetas(v);
175 return;
176 case force_MIG:
177 case force_DIVMIG:
178 SetMigRates(v);
179 return;
180 case force_REC:
181 SetRecRates(v);
182 return;
183 case force_EXPGROWSTICK:
184 case force_GROW:
185 SetGrowthRates(v);
186 return;
187 case force_LOGSELECTSTICK:
188 case force_LOGISTICSELECTION:
189 SetLogisticSelectionCoefficient(v);
190 return;
191 case force_DISEASE:
192 SetDiseaseRates(v);
193 return;
194 case force_DIVERGENCE:
195 SetEpochTimes(v);
196 return;
197 case force_REGION_GAMMA:
198 msg += "\"" + ToString(tag) + ",\" which should never happen. (This is a ";
199 msg += "pseudo-force that\'s only treated as a force within certain ";
200 msg += "contexts of the program.)";
201 throw implementation_error(msg);
202 case force_NONE:
203 msg += "\"" + ToString(tag) + ",\" which should never happen.";
204 throw implementation_error(msg);
205 }
206
207 msg += "\"" + ToString(tag) + ",\" but it does not know how to use it.";
208 throw implementation_error(msg);
209 } // SetRegionalParametersByTag
210
211 //------------------------------------------------------------------------------------
212
SetGlobalThetas(const DoubleVec1d & v)213 void ForceParameters::SetGlobalThetas(const DoubleVec1d& v)
214 {
215 assert(m_space != unknown_region);
216 m_globalThetas = v;
217 if (m_space == known_region)
218 {
219 FillRegionalParamsFromGlobalParams();
220 }
221 } // SetThetas
222
223 //------------------------------------------------------------------------------------
224
SetRegionalThetas(const DoubleVec1d & v)225 void ForceParameters::SetRegionalThetas(const DoubleVec1d& v)
226 {
227 assert(m_space != global_region);
228 m_regionalThetas = v;
229 if (m_space == known_region)
230 {
231 FillGlobalParamsFromRegionalParams();
232 }
233 } // SetThetas
234
235 //------------------------------------------------------------------------------------
236
SetMigRates(const DoubleVec1d & v)237 void ForceParameters::SetMigRates(const DoubleVec1d& v)
238 {
239 long rowsize = LongSquareRootOfLong(v.size());
240
241 // put the given values into migrates, but never put a
242 // non-zero value into the diagonal entries!
243
244 long i, j;
245 long index = 0;
246 m_migrates.clear();
247 for (i = 0; i < rowsize; ++i)
248 {
249 for (j = 0; j < rowsize; ++j)
250 {
251 if (i == j)
252 {
253 m_migrates.push_back(0.0);
254 assert(v[index] == 0.0);
255 }
256 else m_migrates.push_back(v[index]);
257 ++index;
258 }
259 }
260 } // SetMigRates
261
262 //------------------------------------------------------------------------------------
263
SetRecRates(const DoubleVec1d & v)264 void ForceParameters::SetRecRates(const DoubleVec1d& v)
265 {
266 assert(static_cast<long> (v.size()) == 1);// program supports only 1 rec rate
267 m_recrates = v;
268 //LS NOTE: if this gets split into SetRegionalRecRates and SetGlobalRecRates
269 // we will need to also fill the appropriate parallel vector if m_space is
270 // known_region (see SetGlobalThetas)
271 } // SetRecRates
272
273 //------------------------------------------------------------------------------------
274
SetGrowthRates(const DoubleVec1d & v)275 void ForceParameters::SetGrowthRates(const DoubleVec1d& v)
276 {
277 m_growths = v;
278 } // SetGrowths
279
280 //------------------------------------------------------------------------------------
281
SetLogisticSelectionCoefficient(const DoubleVec1d & v)282 void ForceParameters::SetLogisticSelectionCoefficient(const DoubleVec1d& v)
283 {
284 if (1 != v.size())
285 {
286 string msg = "ForceParameters::SetLogisticSelectionCoefficient() received ";
287 msg += "a vector of " + ToString(v.size()) + " elements. This vector must ";
288 msg += "contain exactly one element, corresponding to the selection coefficient.";
289 throw implementation_error(msg);
290 }
291
292 m_logisticSelectionCoefficient = v;
293 } // SetLogisticSelectionCoefficient
294
295 //------------------------------------------------------------------------------------
296
SetDiseaseRates(const DoubleVec1d & v)297 void ForceParameters::SetDiseaseRates(const DoubleVec1d& v)
298 {
299 long rowsize = LongSquareRootOfLong(v.size());
300
301 // put the given values into disease rates, but never put a
302 // non-zero value into the diagonal entries!
303
304 long i, j;
305 long index = 0;
306 m_diseases.clear();
307 for (i = 0; i < rowsize; ++i)
308 {
309 for (j = 0; j < rowsize; ++j)
310 {
311 if (i == j)
312 {
313 m_diseases.push_back(0.0);
314 assert(v[index] == 0.0);
315 }
316 else
317 {
318 m_diseases.push_back(v[index]);
319 }
320 ++index;
321 }
322 }
323 } // SetDiseaseRates
324
325 //------------------------------------------------------------------------------------
326
SetEpochTimes(const DoubleVec1d & v)327 void ForceParameters::SetEpochTimes(const DoubleVec1d& v)
328 {
329 m_epochtimes = v;
330 }
331
332 //------------------------------------------------------------------------------------
333
SetGlobalParameters(const DoubleVec1d & v)334 void ForceParameters::SetGlobalParameters(const DoubleVec1d& v)
335 {
336 assert(m_space != unknown_region);
337 SetParameters(v, true);
338 } // SetGlobalParameters
339
340 //------------------------------------------------------------------------------------
341
SetRegionalParameters(const DoubleVec1d & v)342 void ForceParameters::SetRegionalParameters(const DoubleVec1d& v)
343 {
344 assert(m_space != global_region);
345 SetParameters(v, false);
346 } // SetRegionalParameters
347
348 //------------------------------------------------------------------------------------
349 // SetParameters is private; use the two above functions.
350
SetParameters(const DoubleVec1d & v,bool isGlobal)351 void ForceParameters::SetParameters(const DoubleVec1d& v, bool isGlobal)
352 {
353 for (unsigned long fnum = 0, pnum = 0; fnum<m_forceTags.size(); fnum++)
354 {
355 DoubleVec1d oneForceVec;
356 for (long fpnum = 0; fpnum < m_forceSizes[fnum]; pnum++, fpnum++)
357 {
358 assert (v.size() > static_cast<unsigned long>(pnum));
359 oneForceVec.push_back(v[pnum]);
360 }
361 if (isGlobal)
362 {
363 SetGlobalParametersByTag(m_forceTags[fnum], oneForceVec);
364 }
365 else
366 {
367 SetRegionalParametersByTag(m_forceTags[fnum], oneForceVec);
368 }
369 }
370 }
371
372 //------------------------------------------------------------------------------------
373
GetGlobalThetas() const374 const DoubleVec1d& ForceParameters::GetGlobalThetas() const
375 {
376 assert(m_space != unknown_region);
377 return m_globalThetas;
378 }
379
380 //------------------------------------------------------------------------------------
381
GetRegionalThetas() const382 const DoubleVec1d& ForceParameters::GetRegionalThetas()const
383 {
384 assert(m_space != global_region);
385 return m_regionalThetas;
386 }
387
388 //------------------------------------------------------------------------------------
389
GetGlobalLogParameters() const390 DoubleVec1d ForceParameters::GetGlobalLogParameters() const
391 {
392 assert(m_space != unknown_region);
393 return GetLogParameters(true);
394 }
395
396 //------------------------------------------------------------------------------------
397
GetRegionalLogParameters() const398 DoubleVec1d ForceParameters::GetRegionalLogParameters() const
399 {
400 assert(m_space != global_region);
401 return GetLogParameters(false);
402 }
403
404 //------------------------------------------------------------------------------------
405
GetOnlyDiseaseRates() const406 DoubleVec1d ForceParameters::GetOnlyDiseaseRates() const
407 {
408 DoubleVec1d returnvec;
409 long nallparams(m_diseases.size());
410 long ncontigparams(sqrt(nallparams));
411 long ind(1); // skip the first element as it is always invalid
412 while (ind < nallparams)
413 {
414 long startind(ind);
415 for(; ind < startind+ncontigparams; ++ind)
416 returnvec.push_back(m_diseases[ind]);
417 ++ind; // there's always 1 invalid param between sets of valid
418 // params
419 }
420 assert(static_cast<long>(returnvec.size()) == ncontigparams * (ncontigparams-1));
421 return returnvec;
422 }
423
424 //------------------------------------------------------------------------------------
425
426 //GetLogParameters is private; use the above two calls in code.
GetLogParameters(bool isGlobal) const427 DoubleVec1d ForceParameters::GetLogParameters(bool isGlobal) const
428 {
429 DoubleVec1d tempvec;
430 DoubleVec1d resultvec;
431
432 for (unsigned long fnum=0; fnum<m_forceTags.size(); fnum++)
433 {
434 if (isGlobal)
435 {
436 tempvec = GetGlobalParametersByTag(m_forceTags[fnum]);
437 }
438 else
439 {
440 tempvec = GetRegionalParametersByTag(m_forceTags[fnum]);
441 }
442 if (m_forceTags[fnum] != force_GROW)
443 {
444 LogVec0(tempvec, tempvec);
445 }
446 resultvec.insert(resultvec.end(), tempvec.begin(), tempvec.end());
447 }
448
449 return resultvec;
450 } // GetLogParameters
451
452 //------------------------------------------------------------------------------------
453
GetRegionalLogParameter(long pnum) const454 double ForceParameters::GetRegionalLogParameter(long pnum) const
455 {
456 DoubleVec1d tempvec = GetRegionalParameters();
457 const ParamVector paramvec(true);
458 if (paramvec[pnum].IsForce(force_GROW))
459 {
460 return tempvec[pnum];
461 }
462 else
463 {
464 return SafeLog(tempvec[pnum]);
465 }
466 }
467
468 //------------------------------------------------------------------------------------
469 //------------------------------------------------------------------------------------
470
GetGlobalParameters() const471 DoubleVec1d ForceParameters::GetGlobalParameters() const
472 {
473 assert(m_space != unknown_region);
474 return GetParameters(true);
475 }
476
477 //------------------------------------------------------------------------------------
478
GetRegionalParameters() const479 DoubleVec1d ForceParameters::GetRegionalParameters() const
480 {
481 assert(m_space != global_region);
482 return GetParameters(false);
483 }
484
485 //------------------------------------------------------------------------------------
486
GetParameters(bool isGlobal) const487 DoubleVec1d ForceParameters::GetParameters(bool isGlobal) const
488 {
489 DoubleVec1d tempvec;
490 DoubleVec1d resultvec;
491
492 for (unsigned long fnum=0; fnum<m_forceTags.size(); fnum++)
493 {
494 if (isGlobal)
495 {
496 tempvec = GetGlobalParametersByTag(m_forceTags[fnum]);
497 }
498 else
499 {
500 tempvec = GetRegionalParametersByTag(m_forceTags[fnum]);
501 }
502 resultvec.insert(resultvec.end(), tempvec.begin(), tempvec.end());
503 }
504
505 return resultvec;
506 } // GetParameters
507
508 //------------------------------------------------------------------------------------
509 //------------------------------------------------------------------------------------
510
GetRegionalParametersByTag(force_type tag) const511 const DoubleVec1d& ForceParameters::GetRegionalParametersByTag(force_type tag) const
512 {
513 assert(m_space != global_region);
514 string msg = "ForceParameters::GetRegionalParametersByTag() received tag ";
515 switch (tag)
516 {
517 case force_COAL:
518 return GetRegionalThetas();
519 case force_MIG:
520 case force_DIVMIG:
521 return GetMigRates();
522 case force_REC:
523 return GetRecRates();
524 case force_EXPGROWSTICK:
525 case force_GROW:
526 return GetGrowthRates();
527 case force_DISEASE:
528 return GetDiseaseRates();
529 case force_LOGSELECTSTICK:
530 case force_LOGISTICSELECTION:
531 return GetLogisticSelectionCoefficient();
532 case force_DIVERGENCE:
533 return GetEpochTimes();
534 case force_REGION_GAMMA:
535 msg += "\"" + ToString(tag) + ",\" which should never happen. (This is a ";
536 msg += "pseudo-force that\'s only treated as a force within certain ";
537 msg += "contexts of the program.)";
538 throw implementation_error(msg);
539 case force_NONE:
540 msg += "\"" + ToString(tag) + ",\" which should never happen.";
541 throw implementation_error(msg);
542 }
543
544 msg += "\"" + ToString(tag) + ",\" but it does not know how to use it.";
545 throw implementation_error(msg);
546 } // GetParametersByTag
547
548 //------------------------------------------------------------------------------------
549
GetGlobalParametersByTag(force_type tag) const550 const DoubleVec1d& ForceParameters::GetGlobalParametersByTag(force_type tag) const
551 {
552 assert(m_space != unknown_region);
553 string msg = "ForceParameters::GetGlobalParametersByTag() received tag ";
554 switch (tag)
555 {
556 case force_COAL:
557 return GetGlobalThetas();
558 case force_MIG:
559 case force_DIVMIG:
560 return GetMigRates();
561 case force_REC:
562 return GetRecRates();
563 case force_EXPGROWSTICK:
564 case force_GROW:
565 return GetGrowthRates();
566 case force_LOGSELECTSTICK:
567 case force_LOGISTICSELECTION:
568 return GetLogisticSelectionCoefficient();
569 case force_DISEASE:
570 return GetDiseaseRates();
571 case force_DIVERGENCE:
572 return GetEpochTimes();
573 case force_REGION_GAMMA:
574 msg += "\"" + ToString(tag) + ",\" which should never happen. (This is a ";
575 msg += "pseudo-force that\'s only treated as a force within certain ";
576 msg += "contexts of the program.)";
577 throw implementation_error(msg);
578 case force_NONE:
579 msg += "\"" + ToString(tag) + ",\" which should never happen.";
580 throw implementation_error(msg);
581 }
582
583 msg += "\"" + ToString(tag) + ",\" but it does not know how to use it.";
584 throw implementation_error(msg);
585 } // GetParametersByTag
586
587 //------------------------------------------------------------------------------------
588
GetGlobal2dRates(force_type tag) const589 DoubleVec2d ForceParameters::GetGlobal2dRates(force_type tag) const
590 {
591 assert(tag == force_DISEASE || tag == force_MIG || tag == force_DIVMIG);
592 return SquareOffVector(GetGlobalParametersByTag(tag));
593 } // Get2dRates
594
595 //------------------------------------------------------------------------------------
596
GetRegional2dRates(force_type tag) const597 DoubleVec2d ForceParameters::GetRegional2dRates(force_type tag) const
598 {
599 assert(tag == force_DISEASE || tag == force_MIG || tag == force_DIVMIG);
600 return SquareOffVector(GetRegionalParametersByTag(tag));
601 } // Get2dRates
602
603 //------------------------------------------------------------------------------------
604
FillRegionalParamsFromGlobalParams()605 void ForceParameters::FillRegionalParamsFromGlobalParams()
606 {
607 if (m_region == FLAGLONG)
608 {
609 assert(false);
610 throw implementation_error("Unable to convert global parameters to regional parameters"
611 " because we don't have a region to convert to.");
612 }
613 double effPopSize = registry.GetDataPack().GetRegion(m_region).GetEffectivePopSize();
614
615 m_regionalThetas = m_globalThetas;
616 std::transform(m_globalThetas.begin(),
617 m_globalThetas.end(),
618 m_regionalThetas.begin(),
619 std::bind2nd(std::multiplies<double>(),effPopSize));
620 }
621
622 //------------------------------------------------------------------------------------
623
FillGlobalParamsFromRegionalParams()624 void ForceParameters::FillGlobalParamsFromRegionalParams()
625 {
626 if (m_region == FLAGLONG)
627 {
628 assert(false);
629 throw implementation_error("Unable to convert regional parameters to global parameters"
630 " because we don't know what region we're converting from.");
631 }
632
633 double effPopSize = registry.GetDataPack().GetRegion(m_region).GetEffectivePopSize();
634
635 m_globalThetas = m_regionalThetas;
636 std::transform(m_regionalThetas.begin(),
637 m_regionalThetas.end(),
638 m_globalThetas.begin(),
639 std::bind2nd(std::divides<double>(),effPopSize));
640 }
641
642 //------------------------------------------------------------------------------------
643
GetRegionalScalars() const644 DoubleVec1d ForceParameters::GetRegionalScalars() const
645 {
646 DoubleVec1d resultvec;
647 for (unsigned long force=0; force<m_forceTags.size(); force++)
648 {
649 double val = 1.0;
650 if (m_forceTags[force] == force_COAL)
651 {
652 val = registry.GetDataPack().GetRegion(m_region).GetEffectivePopSize();
653 }
654 if (force_REGION_GAMMA == m_forceTags[force])
655 continue; // should never hit this, but let's put it here for safety
656 DoubleVec1d tempvec(m_forceSizes[force], val);
657 resultvec.insert(resultvec.end(), tempvec.begin(), tempvec.end());
658 }
659
660 return resultvec;
661 } // GetRegionalScalars
662
663 //------------------------------------------------------------------------------------
664
WriteForceParameters(ofstream & sumout,long numtabs) const665 void ForceParameters::WriteForceParameters( ofstream& sumout, long numtabs) const
666 {
667 string tabs(numtabs, '\t');
668 if ( sumout.is_open() )
669 {
670 vector<double> vd;
671
672 sumout << tabs << xmlsum::ESTIMATES_START << endl;
673
674 vd = GetGlobalThetas();
675 if (vd.size() > 0)
676 {
677 sumout << tabs << "\t" << xmlsum::THETAS_START << " ";
678 SumFileHandler::WriteVec1D( sumout, vd );
679 sumout << xmlsum::THETAS_END << endl;
680 }
681
682 vd = GetMigRates();
683 if (vd.size() > 0)
684 {
685 sumout << tabs << "\t" << xmlsum::MIGRATES_START << " ";
686 SumFileHandler::WriteVec1D( sumout, vd );
687 sumout << xmlsum::MIGRATES_END << endl;
688 }
689
690 vd = GetRecRates();
691 if (vd.size() > 0)
692 {
693 sumout << tabs << "\t" << xmlsum::RECRATES_START << " ";
694 SumFileHandler::WriteVec1D( sumout, vd );
695 sumout << xmlsum::RECRATES_END << endl;
696 }
697 vd = GetGrowthRates();
698 if (vd.size() > 0)
699 {
700 sumout << tabs << "\t" << xmlsum::GROWTHRATES_START << " ";
701 SumFileHandler::WriteVec1D( sumout, vd );
702 sumout << xmlsum::GROWTHRATES_END << endl;
703 }
704 vd = GetLogisticSelectionCoefficient();
705 if (vd.size() > 0)
706 {
707 sumout << tabs << "\t" << xmlsum::LOGISTICSELECTION_START << " ";
708 SumFileHandler::WriteVec1D( sumout, vd );
709 sumout << xmlsum::LOGISTICSELECTION_END << endl;
710 }
711 vd = GetDiseaseRates();
712 if (vd.size() > 0)
713 {
714 sumout << tabs << "\t" << xmlsum::DISEASERATES_START << " ";
715 SumFileHandler::WriteVec1D( sumout, vd );
716 sumout << xmlsum::DISEASERATES_END << endl;
717 }
718 vd = GetEpochTimes();
719 if (vd.size() > 0)
720 {
721 sumout << tabs << "\t" << xmlsum::EPOCHTIMES_START << " ";
722 SumFileHandler::WriteVec1D( sumout, vd );
723 sumout << xmlsum::EPOCHTIMES_END << endl;
724 }
725
726 if (global_region == m_space)
727 {
728 const RegionGammaInfo *pRegionGammaInfo = registry.GetRegionGammaInfo();
729 if (pRegionGammaInfo)
730 {
731 vd.clear();
732 vd.push_back(pRegionGammaInfo->GetMLE());
733 sumout << tabs << "\t" << xmlsum::GAMMAOVERREGIONS_START << " ";
734 SumFileHandler::WriteVec1D( sumout, vd );
735 sumout << xmlsum::GAMMAOVERREGIONS_END << endl;
736 }
737 }
738 sumout << tabs << xmlsum::ESTIMATES_END << endl;
739 }
740 else
741 SumFileHandler::HandleSumOutFileFormatError("WriteForceParameters");
742 } // WriteForceParameters
743
744 //____________________________________________________________________________________
745