1 // $Id: parsetreetosettings.cpp,v 1.83 2012/05/15 06:21:48 ewalkup 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 "constants.h"
14 #include "errhandling.h"
15 #include "local_build.h"
16 #include "mathx.h"
17 #include "parsetreetosettings.h"
18 #include "parsetreewalker.h"
19 #include "stringx.h"
20 #include "ui_constants.h" // for uiconst::GLOBAL_ID
21 #include "ui_interface.h"
22 #include "ui_regid.h"
23 #include "ui_strings.h"
24 #include "ui_vars.h"
25 #include "xml.h"
26 #include "xml_strings.h"
27
28 using std::string;
29
30 //------------------------------------------------------------------------------------
31
ParseTreeToSettings(XmlParser & parser,UIInterface & ui)32 ParseTreeToSettings::ParseTreeToSettings(XmlParser& parser, UIInterface& ui)
33 :
34 ParseTreeWalker(parser),
35 uiInterface(ui)
36 {
37 }
38
~ParseTreeToSettings()39 ParseTreeToSettings::~ParseTreeToSettings()
40 {
41 }
42
43 void
ProcessFileSettings()44 ParseTreeToSettings::ProcessFileSettings()
45 {
46 TiXmlElement * docElement = m_parser.GetRootElement();
47
48 const char * tagValue = docElement->Value();
49 string tagString(tagValue);
50 bool matches = CaselessStrCmp(xmlstr::XML_TAG_LAMARC,tagValue);
51 if(!matches)
52 {
53 string msg = m_parser.GetFileName() + ": " + xmlstr::XML_ERR_NOT_LAMARC;
54 m_parser.ThrowDataError(msg);
55 }
56
57 vector<TiXmlElement *> globalModelElements
58 = getAllOptionalDescendantElements(docElement,xmlstr::XML_TAG_MODEL);
59 for (unsigned long i=0; i<globalModelElements.size(); ++i)
60 {
61 DoDLModel(globalModelElements[i],uiconst::GLOBAL_ID);
62 }
63
64 DoDataModels(singleRequiredChild(docElement,xmlstr::XML_TAG_DATA));
65 DoTraits(singleRequiredChild(docElement,xmlstr::XML_TAG_DATA));
66
67 TiXmlElement * chainsElement
68 = singleOptionalChild(docElement,xmlstr::XML_TAG_CHAINS);
69 if(chainsElement != NULL)
70 {
71 DoChainParams(chainsElement);
72 }
73
74 TiXmlElement * formatElement
75 = singleOptionalChild(docElement,xmlstr::XML_TAG_FORMAT);
76 if(formatElement != NULL)
77 {
78 DoUserParams(formatElement);
79 }
80
81 // This *must* be after DoChainParams because DIVMIG requires
82 // Bayesian and that the Epoch Arranger be on
83 TiXmlElement * forcesElement
84 = singleOptionalChild(docElement,xmlstr::XML_TAG_FORCES);
85 if(forcesElement != NULL)
86 {
87 DoForces(forcesElement);
88 }
89 }
90
91 void
DoOptionalElement(TiXmlElement * ancestor,string tagName,string uiLabel,UIId id)92 ParseTreeToSettings::DoOptionalElement(TiXmlElement* ancestor, string tagName,string uiLabel, UIId id)
93 {
94 TiXmlElement * node = singleOptionalChild(ancestor,tagName);
95 if(node == NULL)
96 {
97 return;
98 }
99 try
100 {
101 uiInterface.doSet(uiLabel,getNodeText(node),id);
102 }
103 catch (data_error & d)
104 {
105 m_parser.ThrowDataError(d.what(), node->Row());
106 }
107 }
108
109 void
DoRequiredElement(TiXmlElement * ancestor,string tagName,string uiLabel,UIId id)110 ParseTreeToSettings::DoRequiredElement(TiXmlElement* ancestor, string tagName,string uiLabel, UIId id)
111 {
112 TiXmlElement * node = singleRequiredChild(ancestor,tagName);
113 if(node == NULL)
114 {
115 m_parser.ThrowDataError(xmlstr::XML_ERR_MISSING_TAG_HIER_0
116 + tagName
117 + xmlstr::XML_ERR_MISSING_TAG_HIER_1
118 + ancestor->Value()
119 + xmlstr::XML_ERR_MISSING_TAG_HIER_2,
120 node->Row()
121 );
122 }
123 uiInterface.doSet(uiLabel,getNodeText(node),id);
124 }
125
126 void
DoDataModels(TiXmlElement * dataElement)127 ParseTreeToSettings::DoDataModels(TiXmlElement * dataElement)
128 {
129 TiXmlNode * child = NULL;
130
131 long regionNumber = 0;
132
133 while((child = dataElement->IterateChildren(xmlstr::XML_TAG_REGION,child)))
134 {
135 TiXmlElement * regionElement = child->ToElement();
136 DoRegionDataModels(regionElement,regionNumber);
137 regionNumber++;
138 }
139 }
140
141 void
DoRegionDataModels(TiXmlElement * regionElement,long regionId)142 ParseTreeToSettings::DoRegionDataModels(TiXmlElement * regionElement, long regionId)
143 {
144 TiXmlNode * child = NULL;
145
146 long maxLoci = uiInterface.GetCurrentVars().datapackplus.GetMaxLoci();
147 long locusId = 0;
148
149 while((child = regionElement->IterateChildren(xmlstr::XML_TAG_MODEL,child)))
150 {
151 TiXmlElement * modelElement = child->ToElement();
152 // EWFIX.P4 LOCUS.HACK.HACK.HACK
153 // probably making more than we need
154 // also need to handle per-locus models later
155 // probably does something awful for mismatched types
156 //
157 //JDEBUG--evil hack
158 // we now handle multiple loci, but...ick...sequential implied
159 // ordering
160 if (locusId >= maxLoci) // too many loci?
161 {
162 std::string errstring("Warning: More datamodels found");
163 errstring += " than segments, only the first ";
164 errstring += ToString(locusId);
165 errstring += " datamodels were used.\n";
166 uiInterface.AddWarning(errstring);
167 return;
168 }
169
170 DoDLModel(modelElement,regionId,locusId);
171 ++locusId;
172 }
173 // JDEBUG
174 // ideally we would now fill out the remaining loci (if any) with
175 // default datamodels for their datatype...but this may conflict
176 // with registry::InstallDataModel/registry::CreateDataModel and
177 // it's use of the single "regional" datamodel
178 }
179
180 void
DoDLModel(TiXmlElement * modelElement,long regionId,long locusId)181 ParseTreeToSettings::DoDLModel(TiXmlElement * modelElement, long regionId, long locusId)
182 {
183 string modelName = getNodeAttributeValue(modelElement,xmlstr::XML_ATTRTYPE_NAME);
184 //Here is where we change the regionId to the correct flag if it had been
185 // set to GLOBAL_ID. It must be recursive for the KAllele case (for now).
186 if (regionId == uiconst::GLOBAL_ID)
187 {
188 switch (ProduceModelTypeOrBarf(modelName))
189 {
190 case F84:
191 case GTR:
192 DoDLModel(modelElement, uiconst::GLOBAL_DATAMODEL_NUC_ID);
193 return;
194 case Brownian:
195 case Stepwise:
196 case MixedKS:
197 DoDLModel(modelElement, uiconst::GLOBAL_DATAMODEL_MSAT_ID);
198 return;
199 case KAllele:
200 //Here we have to set the global model for *both* the microsat *and*
201 // 'kallele' data.
202 DoDLModel(modelElement, uiconst::GLOBAL_DATAMODEL_MSAT_ID);
203 DoDLModel(modelElement, uiconst::GLOBAL_DATAMODEL_KALLELE_ID);
204 return;
205 }
206 assert(false); //Uncaught model type.
207 regionId = uiconst::GLOBAL_DATAMODEL_NUC_ID;
208 }
209 assert(regionId != uiconst::GLOBAL_ID);
210
211 UIId locusUIId(regionId,locusId);
212
213 uiInterface.doSet(uistr::dataModel,modelName,locusUIId);
214
215 // comments below show which models should allow the setting
216
217 // base frequencies -- F84, GTR (but don't allow calculated)
218 TiXmlElement * baseFrequenciesElement
219 = singleOptionalChild(modelElement,xmlstr::XML_TAG_BASE_FREQS);
220 if(baseFrequenciesElement != NULL)
221 {
222 DoBaseFrequencies(baseFrequenciesElement,locusUIId);
223 }
224 // ttratio -- F84
225 DoOptionalElement(modelElement,xmlstr::XML_TAG_TTRATIO,uistr::TTRatio,locusUIId);
226 // normalize -- F84 GTR Stepwise KAllele
227 DoOptionalElement(modelElement,xmlstr::XML_TAG_NORMALIZE,uistr::normalization,locusUIId);
228 // gtrrates -- GTR
229 TiXmlElement * gtrElement
230 = singleOptionalChild(modelElement,xmlstr::XML_TAG_GTRRATES);
231 if(gtrElement != NULL)
232 {
233 DoGTR(gtrElement,locusUIId);
234 }
235 // error rate -- GTR and F84
236 DoOptionalElement(modelElement,xmlstr::XML_TAG_PER_BASE_ERROR_RATE,uistr::perBaseErrorRate,locusUIId);
237 // alpha -- mixedKS
238 DoOptionalElement(modelElement,xmlstr::XML_TAG_ALPHA,uistr::alpha,locusUIId);
239 // optimize -- mixedKS
240 DoOptionalElement(modelElement,xmlstr::XML_TAG_ISOPT,uistr::optimizeAlpha,locusUIId);
241
242 // categories F84 GTR Stepwise Brownian KAllele
243 TiXmlElement * categoriesElement
244 = singleOptionalChild(modelElement,xmlstr::XML_TAG_CATEGORIES);
245 if(categoriesElement != NULL)
246 {
247 DoCategories(categoriesElement,locusUIId);
248 }
249
250 // relative mutation rate -- all
251 DoOptionalElement(modelElement,xmlstr::XML_TAG_RELATIVE_MURATE,uistr::relativeMuRate,locusUIId);
252 } // DoDLModel
253
254 void
DoCategories(TiXmlElement * categoriesElement,UIId locusUIId)255 ParseTreeToSettings::DoCategories(TiXmlElement * categoriesElement, UIId locusUIId)
256 {
257 DoOptionalElement(categoriesElement,xmlstr::XML_TAG_NUM_CATEGORIES,uistr::categoryCount,locusUIId);
258 DoOptionalElement(categoriesElement,xmlstr::XML_TAG_AUTOCORRELATION,uistr::autoCorrelation,locusUIId);
259
260 TiXmlElement * catRates = singleOptionalChild(categoriesElement,xmlstr::XML_TAG_RATES);
261 TiXmlElement * catProbs = singleOptionalChild(categoriesElement,xmlstr::XML_TAG_PROBABILITIES);
262
263 if(catRates != NULL)
264 {
265 StringVec1d strings = getNodeTextSplitOnWhitespace(catRates);
266 long catCount = uiInterface.doGetLong(uistr::categoryCount,locusUIId);
267 long numToSet = catCount;
268 if((long)strings.size() > catCount)
269 {
270 uiInterface.AddWarning("Warning: The number of supplied category rates and/or probabilities is greater than the number of categories. Discarding the extras.");
271 }
272 else
273 {
274 if((long)strings.size() < catCount)
275 {
276 uiInterface.AddWarning("Warning: The number of supplied category rates and/or probabilities is less than the number of categories. Using defaults for the extras.");
277 numToSet = strings.size();
278 }
279 }
280
281 for(long index=0;index < numToSet; index++)
282 {
283 UIId innerId(locusUIId.GetIndex1(),locusUIId.GetIndex2(),index);
284 uiInterface.doSet(uistr::categoryRate,strings[index],innerId);
285 }
286 }
287 if(catProbs != NULL)
288 {
289 StringVec1d strings = getNodeTextSplitOnWhitespace(catProbs);
290 long catCount = uiInterface.doGetLong(uistr::categoryCount,locusUIId);
291 long numToSet = catCount;
292 if((long)strings.size() > catCount)
293 {
294 uiInterface.AddWarning("Warning: The number of supplied category rates and/or probabilities is greater than the number of categories. Discarding the extras.");
295 }
296 else
297 {
298 if((long)strings.size() < catCount)
299 {
300 uiInterface.AddWarning("Warning: The number of supplied category rates and/or probabilities is less than the number of categories. Using defaults for the extras.");
301 numToSet = strings.size();
302 }
303 }
304
305 for(long index=0;index < numToSet; index++)
306 {
307 UIId innerId(locusUIId.GetIndex1(),locusUIId.GetIndex2(),index);
308 uiInterface.doSet(uistr::categoryProbability,strings[index],innerId);
309 }
310 }
311 }
312
313 void
DoBaseFrequencies(TiXmlElement * baseFreqsElem,UIId thisId)314 ParseTreeToSettings::DoBaseFrequencies(TiXmlElement * baseFreqsElem, UIId thisId)
315 {
316
317 // kinda grotty, but we want to ignore case differences
318 string frequencies = getNodeText(baseFreqsElem);
319 string tag = xmlstr::XML_TAG_CALCULATED;
320 LowerCase(frequencies);
321 LowerCase(tag);
322 string::size_type index = frequencies.find(tag);
323
324 if(index != std::string::npos)
325 {
326 uiInterface.doSet(uistr::freqsFromData,"true",thisId);
327 }
328 else
329 {
330 StringVec1d strings;
331 bool gotStrings = FromString(frequencies,strings);
332 if(gotStrings)
333 {
334 if(strings.size() == 4)
335 {
336 uiInterface.doSet(uistr::baseFrequencyA,strings[0],thisId);
337 uiInterface.doSet(uistr::baseFrequencyC,strings[1],thisId);
338 uiInterface.doSet(uistr::baseFrequencyG,strings[2],thisId);
339 uiInterface.doSet(uistr::baseFrequencyT,strings[3],thisId);
340 }
341 else if (strings.size() == 0)
342 {
343 uiInterface.AddWarning("Warning: no supplied base frequencies; using a set of defaults.");
344
345 }
346 else
347 {
348 uiInterface.AddWarning("Warning: wrong number of supplied base frequencies; using a set of defaults.");
349 }
350 }
351 else
352 {
353 uiInterface.AddWarning("Warning: no supplied base frequencies; using a set of defaults.");
354 }
355 }
356 } // DoBaseFrequencies
357
358 void
DoGTR(TiXmlElement * gtrElem,UIId thisId)359 ParseTreeToSettings::DoGTR(TiXmlElement * gtrElem, UIId thisId)
360 {
361 StringVec1d strings = getNodeTextSplitOnWhitespace(gtrElem);
362 if(strings.size() == 6)
363 {
364 uiInterface.doSet(uistr::gtrRateAC,strings[0],thisId);
365 uiInterface.doSet(uistr::gtrRateAG,strings[1],thisId);
366 uiInterface.doSet(uistr::gtrRateAT,strings[2],thisId);
367 uiInterface.doSet(uistr::gtrRateCG,strings[3],thisId);
368 uiInterface.doSet(uistr::gtrRateCT,strings[4],thisId);
369 uiInterface.doSet(uistr::gtrRateGT,strings[5],thisId);
370 }
371 else if (strings.size() == 0)
372 {
373 uiInterface.AddWarning("Warning: no GTR rates supplied: using a set of defaults.");
374 }
375 else
376 {
377 uiInterface.AddWarning("Warning: incorrect number of GTR rates supplied: using a set of defaults.");
378 }
379 } // DoGTR
380
381 void
DoChainParams(TiXmlElement * chainsElement)382 ParseTreeToSettings::DoChainParams(TiXmlElement * chainsElement)
383 {
384 DoOptionalElement(chainsElement,xmlstr::XML_TAG_REPLICATES,uistr::replicates);
385
386 TiXmlElement * initialChainElement
387 = singleOptionalChild(chainsElement,xmlstr::XML_TAG_INITIAL);
388 if(initialChainElement != NULL)
389 {
390 DoOptionalElement(initialChainElement,xmlstr::XML_TAG_NUMBER,uistr::initialChains);
391 DoOptionalElement(initialChainElement,xmlstr::XML_TAG_SAMPLES,uistr::initialSamples);
392 DoOptionalElement(initialChainElement,xmlstr::XML_TAG_INTERVAL,uistr::initialInterval);
393 DoOptionalElement(initialChainElement,xmlstr::XML_TAG_DISCARD,uistr::initialDiscard);
394 }
395
396 TiXmlElement * finalChainElement
397 = singleOptionalChild(chainsElement,xmlstr::XML_TAG_FINAL);
398 if(finalChainElement != NULL)
399 {
400 DoOptionalElement(finalChainElement,xmlstr::XML_TAG_NUMBER,uistr::finalChains);
401 DoOptionalElement(finalChainElement,xmlstr::XML_TAG_SAMPLES,uistr::finalSamples);
402 DoOptionalElement(finalChainElement,xmlstr::XML_TAG_INTERVAL,uistr::finalInterval);
403 DoOptionalElement(finalChainElement,xmlstr::XML_TAG_DISCARD,uistr::finalDiscard);
404 }
405
406 TiXmlElement * heatingElement
407 = singleOptionalChild(chainsElement,xmlstr::XML_TAG_HEATING);
408 if(heatingElement != NULL)
409 {
410 TiXmlElement * temperatures
411 = singleOptionalChild(heatingElement,xmlstr::XML_TAG_TEMPERATURES);
412 if(temperatures != NULL) DoTemperatures(getNodeTextSplitOnWhitespace(temperatures));
413
414 TiXmlElement * intervals
415 = singleOptionalChild(heatingElement,xmlstr::XML_TAG_SWAP_INTERVAL);
416 if(intervals != NULL) DoSwapInterval(intervals);
417
418 DoOptionalElement(heatingElement,xmlstr::XML_TAG_HEATING_STRATEGY,uistr::tempAdapt);
419 }
420
421 // must do processing for XML_TAG_BAYESIAN_ANALYSIS before processing for
422 // the XML_TAG_BAYESIAN strategy element because if they are inconsistent
423 // we need to throw--the BAYESIAN_ANALYSIS tag takes precedence.
424 DoOptionalElement(chainsElement,xmlstr::XML_TAG_BAYESIAN_ANALYSIS,uistr::bayesian);
425 TiXmlElement * strategyElement
426 = singleOptionalChild(chainsElement,xmlstr::XML_TAG_STRATEGY);
427 if(strategyElement != NULL)
428 {
429 DoOptionalElement(strategyElement,xmlstr::XML_TAG_RESIMULATING,uistr::dropArranger);
430 uiInterface.GetCurrentVars().chains.RescaleDefaultSizeArranger();
431 DoOptionalElement(strategyElement,xmlstr::XML_TAG_TREESIZE,uistr::sizeArranger);
432 DoOptionalElement(strategyElement,xmlstr::XML_TAG_STAIRARRANGER,uistr::stairArranger);
433 DoOptionalElement(strategyElement,xmlstr::XML_TAG_HAPLOTYPING,uistr::hapArranger);
434 DoOptionalElement(strategyElement,xmlstr::XML_TAG_BAYESIAN,uistr::bayesArranger);
435 DoOptionalElement(strategyElement,xmlstr::XML_TAG_EPOCHSIZEARRANGER,uistr::epochSizeArranger);
436 }
437
438 } // DoChainParams
439
440 //------------------------------------------------------------------------------------
441 //------------------------------------------------------------------------------------
442
443 void
DoTemperatures(StringVec1d temperatures)444 ParseTreeToSettings::DoTemperatures(StringVec1d temperatures)
445 {
446 if(temperatures.empty())
447 {
448 uiInterface.AddWarning("Warning: No supplied temperatures for heated chains: using a set of defaults.");
449 }
450 else
451 {
452 uiInterface.doSet(uistr::heatedChainCount,ToString(temperatures.size()));
453 for(long index = 0; index < (long)temperatures.size(); index++)
454 {
455 uiInterface.doSet(uistr::heatedChain,temperatures[index],UIId(index));
456 }
457 }
458 }
459
460 //------------------------------------------------------------------------------------
461
462 void
DoSwapInterval(TiXmlElement * intervals)463 ParseTreeToSettings::DoSwapInterval(TiXmlElement * intervals)
464 {
465 StringVec1d values = getNodeTextSplitOnWhitespace(intervals);
466 if(values.size() > 0)
467 {
468 uiInterface.doSet(uistr::tempInterval,values[values.size()-1]);
469 }
470 if(values.empty())
471 {
472 uiInterface.AddWarning("Warning: empty <" + xmlstr::XML_TAG_SWAP_INTERVAL + "> tag found; using defaults.");
473 }
474
475 } // DoSwapInterval
476
477 //------------------------------------------------------------------------------------
478
479 void
DoSeedParams(TiXmlElement * formatElement)480 ParseTreeToSettings::DoSeedParams(TiXmlElement * formatElement)
481 {
482 TiXmlElement * node1 = singleOptionalChild(formatElement,xmlstr::XML_TAG_SEED);
483 TiXmlElement * node2 = singleOptionalChild(formatElement,xmlstr::XML_TAG_SEED_FROM_CLOCK);
484 if(node1 != NULL && node2 != NULL)
485 {
486 m_parser.ThrowDataError(xmlstr::XML_ERR_BOTH_SEED_TYPES_0
487 + ToString(node1->Row())
488 + xmlstr::XML_ERR_BOTH_SEED_TYPES_1
489 + ToString(node2->Row())
490 + xmlstr::XML_ERR_BOTH_SEED_TYPES_2
491 + m_parser.GetFileName()
492 + xmlstr::XML_ERR_BOTH_SEED_TYPES_3
493 );
494 }
495
496 DoOptionalElement(formatElement,xmlstr::XML_TAG_SEED,uistr::randomSeed);
497 DoOptionalElement(formatElement,xmlstr::XML_TAG_SEED_FROM_CLOCK,uistr::setOldSeedFromClock);
498 }
499
500 void
DoUserParams(TiXmlElement * formatElement)501 ParseTreeToSettings::DoUserParams(TiXmlElement * formatElement)
502 {
503 DoOptionalElement(formatElement,xmlstr::XML_TAG_VERBOSITY,uistr::verbosity);
504 DoOptionalElement(formatElement,xmlstr::XML_TAG_PROGRESS_REPORTS,uistr::progress);
505
506 TiXmlElement * plottingElement
507 = singleOptionalChild(formatElement,xmlstr::XML_TAG_PLOTTING);
508 if(plottingElement != NULL)
509 {
510 DoOptionalElement(plottingElement,xmlstr::XML_TAG_POSTERIOR,uistr::plotPost);
511 }
512
513 DoSeedParams(formatElement);
514 DoOptionalElement(formatElement,xmlstr::XML_TAG_RESULTS_FILE,uistr::resultsFileName);
515 //LS NOTE: These are off by default, so we only warn that there's a
516 // problem if we turn them off, meaning we should set the filename first.
517 DoOptionalElement(formatElement,xmlstr::XML_TAG_OLD_SUMMARY_FILE,uistr::treeSumInFileName);
518 DoOptionalElement(formatElement,xmlstr::XML_TAG_IN_SUMMARY_FILE,uistr::treeSumInFileName);
519 DoOptionalElement(formatElement,xmlstr::XML_TAG_OUT_SUMMARY_FILE,uistr::treeSumOutFileName);
520 DoOptionalElement(formatElement,xmlstr::XML_TAG_USE_IN_SUMMARY,uistr::treeSumInFileEnabled);
521 DoOptionalElement(formatElement,xmlstr::XML_TAG_USE_OUT_SUMMARY,uistr::treeSumOutFileEnabled);
522 #ifdef LAMARC_QA_TREE_DUMP
523 DoOptionalElement(formatElement,xmlstr::XML_TAG_ARGFILE_PREFIX,uistr::argFilePrefix);
524 DoOptionalElement(formatElement,xmlstr::XML_TAG_USE_ARGFILES,uistr::useArgFiles);
525 DoOptionalElement(formatElement,xmlstr::XML_TAG_MANY_ARGFILES,uistr::manyArgFiles);
526 #endif // LAMARC_QA_TREE_DUMP
527 DoOptionalElement(formatElement,xmlstr::XML_TAG_NEWICKTREEFILE_PREFIX,uistr::newickTreeFilePrefix);
528 DoOptionalElement(formatElement,xmlstr::XML_TAG_USE_NEWICKTREEFILE,uistr::useNewickTreeFiles);
529 //LS NOTE: And these two are on by default, so we warn only if it's still
530 // on when we change the name.
531 DoOptionalElement(formatElement,xmlstr::XML_TAG_USE_CURVEFILES,uistr::useCurveFiles);
532 DoOptionalElement(formatElement,xmlstr::XML_TAG_CURVEFILE_PREFIX,uistr::curveFilePrefix);
533 DoOptionalElement(formatElement,xmlstr::XML_TAG_PROFILE_PREFIX,uistr::profileprefix);
534 DoOptionalElement(formatElement,xmlstr::XML_TAG_USE_RECLOCFILE,uistr::useReclocFiles);
535 DoOptionalElement(formatElement,xmlstr::XML_TAG_RECLOCFILE_PREFIX,uistr::reclocFilePrefix);
536 DoOptionalElement(formatElement,xmlstr::XML_TAG_USE_TRACEFILE,uistr::useTraceFiles);
537 DoOptionalElement(formatElement,xmlstr::XML_TAG_TRACEFILE_PREFIX,uistr::traceFilePrefix);
538 //LS NOTE: And this one we always write, so there are no warnings about
539 // the setting.
540 DoOptionalElement(formatElement,xmlstr::XML_TAG_OUT_XML_FILE,uistr::xmlOutFileName);
541 DoOptionalElement(formatElement,xmlstr::XML_TAG_REPORT_XML_FILE,uistr::xmlReportFileName);
542
543 } // DoUserParams
544
545 //------------------------------------------------------------------------------------
546 // Force parameter functions
547 //------------------------------------------------------------------------------------
548
549 void
DoForces(TiXmlElement * forcesElement)550 ParseTreeToSettings::DoForces(TiXmlElement * forcesElement)
551 {
552 DoForceIfPresent(forcesElement,force_DIVERGENCE,xmlstr::XML_TAG_DIVERGENCE);
553 DoForceIfPresent(forcesElement,force_COAL,xmlstr::XML_TAG_COALESCENCE);
554 DoForceIfPresent(forcesElement,force_MIG,xmlstr::XML_TAG_MIGRATION);
555 DoForceIfPresent(forcesElement,force_REC,xmlstr::XML_TAG_RECOMBINATION);
556 DoForceIfPresent(forcesElement,force_GROW,xmlstr::XML_TAG_GROWTH);
557 DoForceIfPresent(forcesElement,force_LOGISTICSELECTION,
558 xmlstr::XML_TAG_LOGISTICSELECTION);
559 DoForceIfPresent(forcesElement,force_LOGSELECTSTICK,
560 xmlstr::XML_TAG_STOCHASTICSELECTION);
561 DoForceIfPresent(forcesElement,force_DISEASE,xmlstr::XML_TAG_DISEASE);
562 DoForceIfPresent(forcesElement,force_REGION_GAMMA,
563 xmlstr::XML_TAG_REGION_GAMMA);
564 DoForceIfPresent(forcesElement,force_DIVMIG,xmlstr::XML_TAG_DIVMIG);
565
566 } // DoForces
567
568 //------------------------------------------------------------------------------------
569
570 void
DoForceIfPresent(TiXmlElement * forcesElement,force_type forcetype,const string & forcetag)571 ParseTreeToSettings::DoForceIfPresent(
572 TiXmlElement * forcesElement,
573 force_type forcetype,
574 const string& forcetag)
575 {
576 TiXmlElement * forceElement = singleOptionalChild(forcesElement,forcetag);
577
578 if(forceElement != NULL)
579 {
580 DoForce(forceElement,forcetype);
581 }
582
583 }
584
585 //------------------------------------------------------------------------------------
586
587 void
DoForce(TiXmlElement * forceElement,force_type forcetype)588 ParseTreeToSettings::DoForce(TiXmlElement* forceElement, force_type forcetype)
589 {
590 UIId forceId(forcetype);
591
592 string forceOnOff
593 = getNodeAttributeValue(forceElement,xmlstr::XML_ATTRTYPE_VALUE);
594 if(!forceOnOff.empty())
595 {
596 uiInterface.doSet(uistr::forceOnOff,forceOnOff,forceId);
597 }
598 else
599 {
600 uiInterface.doSet(uistr::forceOnOff,"on",forceId);
601 }
602
603 if(forcetype == force_GROW)
604 {
605 string growthType = getNodeAttributeValue(forceElement,xmlstr::XML_ATTRTYPE_TYPE);
606 if(!growthType.empty())
607 {
608 uiInterface.doSet(uistr::growthType,growthType);
609 }
610 }
611
612 if(forcetype == force_LOGSELECTSTICK)
613 {
614 string stype(lamarcstrings::shortSSelectionName);
615 uiInterface.doSet(uistr::selectType,stype);
616 }
617
618 if(forcetype == force_DIVERGENCE)
619 {
620 // parse the population-tree hierarchy
621 TiXmlElement* popElement =
622 singleRequiredChild(forceElement,xmlstr::XML_TAG_POPTREE);
623 TiXmlNode * child = NULL;
624
625 // dreadful quick-n-dirty structures to make sure we have a legal tree
626 std::map<std::string,std::string> popToAnc;
627 std::set<std::string> ancSet;
628 std::map<std::string,std::string>::iterator mit;
629 std::set<std::string>::iterator sit;
630
631 while((child = popElement->IterateChildren(xmlstr::XML_TAG_EPOCH_BOUNDARY,child)))
632 {
633 TiXmlElement * boundaryElement = child->ToElement();
634 TiXmlElement* newpopElement =
635 singleRequiredChild(boundaryElement,xmlstr::XML_TAG_NEWPOP);
636 StringVec1d newpops = getNodeTextSplitOnWhitespace(newpopElement);
637
638 if (newpops.size() != 2)
639 {
640 m_parser.ThrowDataError(xmlstr::XML_ERR_NEWPOP, newpopElement->Row());
641 }
642
643 // JNOTE: not using the SetGet machinery because of obscure "unsolved problems"
644 // fix this right! at some point
645 uiInterface.GetCurrentVars().forces.AddNewPops(newpops);
646
647 TiXmlElement* ancestorElement =
648 singleRequiredChild(boundaryElement,xmlstr::XML_TAG_ANCESTOR);
649 string ancname = getNodeText(ancestorElement);
650 if(ancname.empty())
651 {
652 m_parser.ThrowDataError(xmlstr::XML_ERR_EMPTY_ANCESTOR, ancestorElement->Row());
653 }
654 else
655 {
656 // JNOTE: not using the SetGet machinery because of obscure "unsolved problems"
657 // fix this right! at some point
658 // uiInterface.doSet(uistr::divergenceEpochAncestor,ancname);
659 uiInterface.GetCurrentVars().forces.AddAncestor(ancname);
660 }
661
662
663
664 // check newpops[0] not already a child
665 mit = popToAnc.find(newpops[0]);
666 if(mit != popToAnc.end())
667 {
668 std::string oldAnc = (*mit).second;
669 m_parser.ThrowDataError(xmlstr::XML_ERR_MULTIPLE_ANCESTORS_0
670 + newpops[0]
671 + xmlstr::XML_ERR_MULTIPLE_ANCESTORS_1
672 + oldAnc
673 + xmlstr::XML_ERR_MULTIPLE_ANCESTORS_2
674 + ancname
675 + xmlstr::XML_ERR_MULTIPLE_ANCESTORS_3,
676 ancestorElement->Row()
677 );
678 }
679
680 // check newpops[1] not already a child
681 mit = popToAnc.find(newpops[1]);
682 if(mit != popToAnc.end())
683 {
684 std::string oldAnc = (*mit).second;
685 m_parser.ThrowDataError(xmlstr::XML_ERR_MULTIPLE_ANCESTORS_0
686 + newpops[1]
687 + xmlstr::XML_ERR_MULTIPLE_ANCESTORS_1
688 + oldAnc
689 + xmlstr::XML_ERR_MULTIPLE_ANCESTORS_2
690 + ancname
691 + xmlstr::XML_ERR_MULTIPLE_ANCESTORS_3,
692 ancestorElement->Row()
693 );
694 }
695
696 // check ancname not already used
697 sit = ancSet.find(ancname);
698 if(sit != ancSet.end())
699 {
700 m_parser.ThrowDataError(xmlstr::XML_ERR_DUPLICATE_ANCESTOR_0
701 + ancname
702 + xmlstr::XML_ERR_DUPLICATE_ANCESTOR_1,
703 ancestorElement->Row()
704 );
705
706 }
707
708 // Everything OK! Let's add data -- edit only one of the three at your peril
709 popToAnc[newpops[0]] = ancname;
710 popToAnc[newpops[1]] = ancname;
711 ancSet.insert(ancname);
712
713
714 }
715
716
717 // Now, check that we have the following:
718 // * number of keys mapping pops to ancestors in (2*pops) - 2
719 // num pops + num internal nodes - one root
720 // * exactly one ancestor not a key in map from pops to ancestors
721 // * every pop name is a key in map from pops to ancestors
722 // * no pop name a key in set of ancestors
723
724 size_t numNotRoot = popToAnc.size();
725 size_t numInternal = ancSet.size();
726
727 std::set<std::string> popNamesSpecifiedInData;
728 TiXmlElement * docElement = m_parser.GetRootElement();
729 TiXmlElement * dataElem = singleRequiredChild(docElement,xmlstr::XML_TAG_DATA);
730
731 TiXmlNode * regChild = NULL;
732 while((regChild = dataElem->IterateChildren(xmlstr::XML_TAG_REGION,regChild)))
733 {
734 TiXmlElement * regionElement = regChild->ToElement();
735 TiXmlNode * popChild = NULL;
736 while((popChild = regionElement->IterateChildren(xmlstr::XML_TAG_POPULATION,popChild)))
737 {
738 TiXmlElement * populationElement = popChild->ToElement();
739 std::string thisPopName = getNodeAttributeValue(populationElement,xmlstr::XML_ATTRTYPE_NAME);
740 popNamesSpecifiedInData.insert(thisPopName);
741 }
742 }
743
744
745
746
747 size_t popCount = popNamesSpecifiedInData.size();
748
749
750 for(std::set<std::string>::const_iterator si = popNamesSpecifiedInData.begin(); si != popNamesSpecifiedInData.end(); si++)
751 {
752 std::string popName = *si;
753 mit = popToAnc.find(popName);
754 if (mit == popToAnc.end())
755 {
756 m_parser.ThrowDataError(xmlstr::XML_ERR_MISSING_NEWPOP_0
757 + popName
758 + xmlstr::XML_ERR_MISSING_NEWPOP_1,
759 popElement->Row()
760 );
761 }
762
763 sit = ancSet.find(popName);
764 if(sit != ancSet.end())
765 {
766 m_parser.ThrowDataError(xmlstr::XML_ERR_ANCESTOR_TRUEPOP_0
767 + popName
768 + xmlstr::XML_ERR_ANCESTOR_TRUEPOP_1,
769 popElement->Row()
770 );
771 }
772
773 }
774
775 if(numNotRoot != (popNamesSpecifiedInData.size() * 2) - 2 )
776 {
777 m_parser.ThrowDataError(xmlstr::XML_ERR_BAD_ANCESTOR_TREE,popElement->Row());
778 }
779
780 size_t nonRootAncestors = 0;
781 for(sit = ancSet.begin(); sit != ancSet.end(); sit++)
782 {
783 mit = popToAnc.find(*sit);
784 if (mit != popToAnc.end())
785 {
786 nonRootAncestors++;
787 }
788 }
789
790 if(nonRootAncestors != ancSet.size() - 1)
791 {
792 m_parser.ThrowDataError(xmlstr::XML_ERR_BAD_ANCESTOR_TREE,popElement->Row());
793 }
794
795 // UGH -- this is awful -- we must change to bayesian and force
796 // epoch arranger to be non-zero
797 if(!uiInterface.doGetBool(uistr::bayesian))
798 {
799 // this should do a decent enough job of choosing a value for the arranger frequency
800 uiInterface.doSet(uistr::bayesian,"true");
801 }
802 double bayesFreq = uiInterface.doGetDouble(uistr::bayesArranger);
803 double resimFreq = uiInterface.doGetDouble(uistr::dropArranger);
804 double epochFreq = (bayesFreq > resimFreq) ? resimFreq : bayesFreq;
805 uiInterface.doSet(uistr::epochSizeArranger,ToString(epochFreq));
806
807
808
809 #if 0
810 TiXmlElement* boundaryElement =
811 singleRequiredChild(popElement,xmlstr::XML_TAG_EPOCH_BOUNDARY);
812 do
813 {
814 TiXmlElement* newpopElement =
815 singleRequiredChild(boundaryElement,xmlstr::XML_TAG_NEWPOP);
816 StringVec1d newpops = getNodeTextSplitOnWhitespace(newpopElement);
817 // JNOTE: not using the SetGet machinery because of obscure "unsolved problems"
818 // fix this right! at some point
819 uiInterface.GetCurrentVars().forces.AddNewPops(newpops);
820
821 TiXmlElement* ancestorElement =
822 singleRequiredChild(boundaryElement,xmlstr::XML_TAG_ANCESTOR);
823 string ancname = getNodeText(ancestorElement);
824 // JNOTE: not using the SetGet machinery because of obscure "unsolved problems"
825 // fix this right! at some point
826 uiInterface.GetCurrentVars().forces.AddAncestor(ancname);
827
828 boundaryElement = singleOptionalChild(popElement,xmlstr::XML_TAG_EPOCH_BOUNDARY);
829 } while (boundaryElement != NULL);
830 #endif
831
832 }
833
834 if(uiInterface.doGetBool(uistr::forceLegal,forceId))
835 {
836 DoOptionalElement(forceElement,xmlstr::XML_TAG_MAX_EVENTS,uistr::maxEvents,forceId);
837 DoOptionalElement(forceElement,xmlstr::XML_TAG_DISEASELOCATION,uistr::diseaseLocation,forceId);
838
839 //LS NOTE: Order matters in this next section. First, we do profiles,
840 // because if they're done after the groups, the 'none' entries
841 // override any (correct) 'on' entries for grouped parameters. Next
842 // we do constraints, then groups, because the groups tags need
843 // to override the constraints tags.
844 //
845 // Then we do the priors, because they will be affected by the groups
846 // and constraints.
847 //
848 // Then we do start values, because they will be affected by the
849 // constraints, groups, priors, and methods. Sadly, we have to do a bit
850 // of munging in DoStartValuesAndMethods because the Methods are, in
851 // turn, affected by the process of setting the start values.
852 DoProfiles(forceElement,forceId);
853 DoConstraints(forceElement,forceId);
854 DoGroups(forceElement, forceId);
855 DoPriors(forceElement, forceId);
856 DoStartValuesAndMethods(forceElement,forceId);
857 DoTrueValues(forceElement,forceId);
858
859 // LS DEBUG -- we can at least check the zeroes for validity.
860 // Eventually the method to do this will change, so don't copy this
861 // technique elsewhere.
862 if (!uiInterface.GetCurrentVars().forces.GetForceZeroesValidity(forcetype))
863 {
864 string err = "Invalid settings for force ";
865 err += ToString(forcetype) + ". Too many parameters are set invalid or have a start value of 0.0.";
866 throw data_error(err);
867 }
868 uiInterface.GetCurrentVars().forces.FixGroups(forcetype);
869 }
870 }
871
872 void
DoStartValuesAndMethods(TiXmlElement * forceElement,UIId forceId)873 ParseTreeToSettings::DoStartValuesAndMethods(TiXmlElement* forceElement, UIId forceId)
874 {
875 TiXmlElement * methodsElement =
876 singleOptionalChild(forceElement,xmlstr::XML_TAG_METHOD);
877 TiXmlElement * startValuesElement =
878 singleOptionalChild(forceElement,xmlstr::XML_TAG_START_VALUES);
879
880 if(methodsElement == NULL && startValuesElement == NULL)
881 // simply return and take the default values
882 {
883 return;
884 }
885
886 force_type thisForce = forceId.GetForceType();
887 long expectedNumParameters =
888 uiInterface.GetCurrentVars().forces.GetNumParameters(thisForce);
889 StringVec1d values;
890 StringVec1d methods;
891
892 if(startValuesElement != NULL)
893 {
894 values = getNodeTextSplitOnWhitespace(startValuesElement);
895 if(static_cast<long>(values.size()) != expectedNumParameters)
896 {
897
898 m_parser.ThrowDataError(xmlstr::XML_ERR_START_VALUE_COUNT_0
899 + ToString(expectedNumParameters)
900 + xmlstr::XML_ERR_START_VALUE_COUNT_1
901 + ToString(thisForce)
902 + xmlstr::XML_ERR_START_VALUE_COUNT_2
903 + ToString(values.size())
904 + xmlstr::XML_ERR_START_VALUE_COUNT_3,
905 startValuesElement->Row()
906
907 );
908 }
909 }
910
911 if(methodsElement != NULL)
912 {
913 methods = getNodeTextSplitOnWhitespace(methodsElement);
914 if(static_cast<long>(methods.size()) != expectedNumParameters)
915 {
916 m_parser.ThrowDataError(xmlstr::XML_ERR_METHOD_TYPE_COUNT_0
917 + ToString(expectedNumParameters)
918 + xmlstr::XML_ERR_METHOD_TYPE_COUNT_1
919 + ToString(thisForce)
920 + xmlstr::XML_ERR_METHOD_TYPE_COUNT_2
921 + ToString(methods.size())
922 + xmlstr::XML_ERR_METHOD_TYPE_COUNT_2,
923 methodsElement->Row()
924 );
925 }
926 }
927
928 // at this point, we either have a full set of method types or none
929 // and we have either a full set of start values or none
930
931 if(methodsElement == NULL)
932 // an easy case -- all methods become type USER
933 {
934 for(long index=0; index < (long)values.size(); index++)
935 {
936 UIId id(forceId.GetForceType(),index);
937 uiInterface.doSet(uistr::startValue,values[index],id);
938 uiInterface.doSet(uistr::startValueMethod,ToString(method_USER),id);
939 }
940 return;
941 }
942
943 if(startValuesElement == NULL)
944 // OK as long as method isn't USER
945 {
946 for(long index=0; index < (long)methods.size(); index++)
947 {
948 UIId id(forceId.GetForceType(),index);
949 uiInterface.doSet(uistr::startValueMethod,methods[index],id);
950
951 if(StringMatchesMethodType(methods[index],method_USER))
952 {
953 double defaultStartVal
954 = uiInterface.GetCurrentVars().forces.GetStartValue(forceId.GetForceType(),index);
955 uiInterface.AddWarning(xmlstr::XML_ERR_METHOD_USER_WITHOUT_VALUE_0
956 +ToString(forceId.GetForceType())
957 +xmlstr::XML_ERR_METHOD_USER_WITHOUT_VALUE_1
958 +ToString(defaultStartVal)
959 +xmlstr::XML_ERR_METHOD_USER_WITHOUT_VALUE_2);
960 }
961 }
962 return;
963 }
964
965 // OK. Now we should have a full set of method types and start values
966 assert((long)methods.size() == expectedNumParameters);
967 assert((long)values.size() == expectedNumParameters);
968
969 for(long index=0; index < (long)values.size(); index++)
970 {
971 UIId id(forceId.GetForceType(),index);
972 uiInterface.doSet(uistr::startValue,values[index],id);
973 }
974 DoubleVec1d usersStartValues = uiInterface.GetCurrentVars().forces.GetStartValues(forceId.GetForceType());
975
976 for(long index=0; index < (long)values.size(); index++)
977 {
978 UIId id(forceId.GetForceType(),index);
979 uiInterface.doSet(uistr::startValueMethod,methods[index],id);
980 }
981 DoubleVec1d overriddenStartValues = uiInterface.GetCurrentVars().forces.GetStartValues(forceId.GetForceType());
982
983 assert(usersStartValues.size() == overriddenStartValues.size());
984 assert(usersStartValues.size() == methods.size());
985
986 for(long index=0; index < (long)methods.size(); index++)
987 {
988 double userVal = usersStartValues[index];
989 double overriddenValue = overriddenStartValues[index];
990
991 if (fabs(userVal - overriddenValue) > 0.00001)
992 {
993 method_type culprit
994 = uiInterface.doGetMethodType(uistr::startValueMethod,UIId(thisForce,index));
995 uiInterface.AddWarning("Warning: setting the start method for a "
996 + ToString(forceId.GetForceType())
997 + " parameter to "
998 + ToString(culprit, true)
999 + " will override the value set in the <"
1000 + xmlstr::XML_TAG_START_VALUES
1001 + "> tag.");
1002 }
1003 }
1004
1005 #if 0 // Potentially DEAD CODE (bobgian, Feb 2010)
1006 bool someStartValues = false;
1007 // set start values first since these calls set methods
1008 if(startValuesElement != NULL)
1009 {
1010 someStartValues = true;
1011 StringVec1d values = getNodeTextSplitOnWhitespace(startValuesElement);
1012 for(long index=0; index < (long)values.size(); index++)
1013 {
1014 UIId id(forceId.GetForceType(),index);
1015 uiInterface.doSet(uistr::startValue,values[index],id);
1016 }
1017 if(values.empty())
1018 {
1019 someStartValues = false;
1020 uiInterface.AddWarning("Warning: empty <" + xmlstr::XML_TAG_START_VALUES + "> tag found; using defaults.");
1021 }
1022
1023 }
1024 else
1025 {
1026 uiInterface.AddWarning("Warning: missing <" + xmlstr::XML_TAG_START_VALUES + "> tag ; using defaults.");
1027 }
1028
1029 DoubleVec1d originalStartValues = uiInterface.GetCurrentVars().forces.GetStartValues(forceId.GetForceType());
1030 if(methodsElement != NULL)
1031 {
1032 StringVec1d methods = getNodeTextSplitOnWhitespace(methodsElement);
1033 for(long index=0; index < (long)methods.size(); index++)
1034 {
1035 UIId id(forceId.GetForceType(),index);
1036 uiInterface.doSet(uistr::startValueMethod,methods[index],id);
1037 }
1038 if(methods.empty())
1039 {
1040 uiInterface.AddWarning("Warning: empty <" + xmlstr::XML_TAG_METHOD + "> tag found; using defaults");
1041 }
1042 else
1043 {
1044 if (someStartValues)
1045 {
1046 DoubleVec1d newStartValues = uiInterface.GetCurrentVars().forces.GetStartValues(forceId.GetForceType());
1047 for (unsigned long i=0; i<originalStartValues.size(); i++)
1048 {
1049 if (fabs(originalStartValues[i] - newStartValues[i]) > 0.00001)
1050 {
1051 method_type culprit = uiInterface.doGetMethodType(uistr::startValueMethod,UIId(forceId.GetForceType(), i));
1052 uiInterface.AddWarning("Warning: setting the start method for a "
1053 + ToString(forceId.GetForceType())
1054 + " parameter to "
1055 + ToString(culprit, true)
1056 + " will override the value set in the <"
1057 + xmlstr::XML_TAG_START_VALUES + "> tag.");
1058 }
1059 }
1060 }
1061 }
1062 }
1063 #endif
1064 }
1065
1066 //------------------------------------------------------------------------------------
1067
1068 void
DoTrueValues(TiXmlElement * forceElement,UIId forceId)1069 ParseTreeToSettings::DoTrueValues(TiXmlElement* forceElement, UIId forceId)
1070 {
1071 TiXmlElement * trueValuesElement =
1072 singleOptionalChild(forceElement,xmlstr::XML_TAG_TRUEVALUE);
1073 if(trueValuesElement != NULL)
1074 {
1075 StringVec1d values = getNodeTextSplitOnWhitespace(trueValuesElement);
1076 for(unsigned long index=0; index < values.size(); index++)
1077 {
1078 UIId id(forceId.GetForceType(),index);
1079 uiInterface.doSet(uistr::trueValue,values[index],id);
1080 }
1081 if(values.empty())
1082 {
1083 uiInterface.AddWarning("Warning: empty <" + xmlstr::XML_TAG_TRUEVALUE + "> tag found; using defaults");
1084 }
1085 }
1086 }
1087
1088 //------------------------------------------------------------------------------------
1089
1090 void
DoProfiles(TiXmlElement * forceElement,UIId forceId)1091 ParseTreeToSettings::DoProfiles(TiXmlElement * forceElement, UIId forceId)
1092 {
1093 TiXmlElement * profilesElement =
1094 singleOptionalChild(forceElement,xmlstr::XML_TAG_PROFILES);
1095
1096 if(profilesElement != NULL)
1097 {
1098 string profilesString = getNodeText(profilesElement);
1099 ProftypeVec1d profiles = ProduceProftypeVec1dOrBarf(profilesString);
1100 size_t index;
1101 bool hasFix = false;
1102 bool hasPerc= false;
1103 for(index=0; index < profiles.size(); index++)
1104 {
1105 UIId localId(forceId.GetForceType(),index);
1106 switch(profiles[index])
1107 {
1108 case profile_PERCENTILE:
1109 hasPerc = true;
1110 uiInterface.doSet(uistr::profileByID,"true",localId);
1111 break;
1112 case profile_FIX:
1113 hasFix = true;
1114 uiInterface.doSet(uistr::profileByID,"true",localId);
1115 break;
1116 case profile_NONE:
1117 uiInterface.doSet(uistr::profileByID,"false",localId);
1118 break;
1119 }
1120 }
1121 if(hasFix && !hasPerc)
1122 {
1123 uiInterface.doSet(uistr::profileByForce,"fixed",forceId);
1124 }
1125 else
1126 {
1127 uiInterface.doSet(uistr::profileByForce,"percentile",forceId);
1128 if(hasFix)
1129 {
1130 // add warning
1131 std::string errstring("Warning: <profiles> element near line ");
1132 errstring += ToString(profilesElement->Row());
1133 errstring += " contains both fixed and percentile profiling.";
1134 errstring += " We only allow one type per force and are";
1135 errstring += " changing the type to percentile profiling.";
1136 uiInterface.AddWarning(errstring);
1137 }
1138 }
1139 }
1140 }
1141
1142 void
DoConstraints(TiXmlElement * forceElement,UIId forceId)1143 ParseTreeToSettings::DoConstraints(TiXmlElement * forceElement, UIId forceId)
1144 {
1145 TiXmlElement * constraintsElement =
1146 singleOptionalChild(forceElement,xmlstr::XML_TAG_CONSTRAINTS);
1147
1148 if(constraintsElement != NULL)
1149 {
1150 string constraintsString = getNodeText(constraintsElement);
1151 vector < ParamStatus > constraints
1152 = ProduceParamstatusVec1dOrBarf(constraintsString);
1153 size_t index;
1154 for(index=0; index < constraints.size(); index++)
1155 {
1156 UIId localID(forceId.GetForceType(),index);
1157 uiInterface.doSet(uistr::constraintType,ToString(constraints[index].Status()) ,localID);
1158 }
1159 }
1160 }
1161
1162 void
DoGroups(TiXmlElement * forceElement,UIId forceId)1163 ParseTreeToSettings::DoGroups(TiXmlElement * forceElement, UIId forceId)
1164 {
1165 TiXmlNode * child = NULL;
1166 size_t index = 0;
1167 while((child = forceElement->IterateChildren(xmlstr::XML_TAG_GROUP,child)))
1168 {
1169 TiXmlElement * groupElement = child->ToElement();
1170 DoGroup(groupElement, forceId, index);
1171 index++;
1172 }
1173 }
1174
1175 void
DoGroup(TiXmlElement * groupElement,UIId forceId,size_t index)1176 ParseTreeToSettings::DoGroup(TiXmlElement * groupElement, UIId forceId,
1177 size_t index)
1178 {
1179 string constraintType
1180 = getNodeAttributeValue(groupElement,xmlstr::XML_ATTRTYPE_CONSTRAINT);
1181 string indicesString = getNodeText(groupElement);
1182
1183 UIId localID(forceId.GetForceType(),index);
1184 uiInterface.doSet(uistr::groupParamList, indicesString, localID);
1185 uiInterface.doSet(uistr::groupConstraintType, constraintType, localID);
1186 }
1187
1188 void
DoPriors(TiXmlElement * forceElement,UIId forceId)1189 ParseTreeToSettings::DoPriors(TiXmlElement * forceElement, UIId forceId)
1190 {
1191 TiXmlNode * child = NULL;
1192 while((child = forceElement->IterateChildren(xmlstr::XML_TAG_PRIOR,child)))
1193 {
1194 TiXmlElement * priorElement = child->ToElement();
1195 DoPrior(priorElement, forceId);
1196 }
1197 }
1198
1199 void
DoPrior(TiXmlElement * priorElement,UIId forceId)1200 ParseTreeToSettings::DoPrior(TiXmlElement * priorElement, UIId forceId)
1201 {
1202 TiXmlElement * parameterElement =
1203 singleOptionalChild(priorElement,xmlstr::XML_TAG_PARAMINDEX);
1204
1205 long paramID = uiconst::GLOBAL_ID;
1206
1207 if (parameterElement)
1208 {
1209 string paramName = getNodeText(parameterElement);
1210 try
1211 {
1212 paramID = ProduceLongOrBarf(paramName);
1213 paramID--;
1214 }
1215 catch (const data_error& e)
1216 {
1217 if ((!CaselessStrCmp(paramName, uistr::defaultStr)) && (!CaselessStrCmp(paramName, uistr::allStr)))
1218 {
1219 m_parser.ThrowDataError("Parameter index must be an integer or 'all'", priorElement->Row());
1220 }
1221 paramID = uiconst::GLOBAL_ID;
1222 }
1223 }
1224 UIId localId(forceId.GetForceType(), paramID);
1225
1226 string priorTypeString
1227 = getNodeAttributeValue(priorElement,xmlstr::XML_ATTRTYPE_TYPE);
1228 uiInterface.doSet(uistr::priorType, priorTypeString, localId);
1229
1230 try
1231 {
1232 TiXmlElement * lowerBoundElement =
1233 singleRequiredChild(priorElement,xmlstr::XML_TAG_PRIORLOWERBOUND);
1234 string lowerString = getNodeText(lowerBoundElement);
1235 uiInterface.doSet(uistr::priorLowerBound, lowerString, localId);
1236
1237 TiXmlElement * upperBoundElement =
1238 singleRequiredChild(priorElement,xmlstr::XML_TAG_PRIORUPPERBOUND);
1239 string upperString = getNodeText(upperBoundElement);
1240 uiInterface.doSet(uistr::priorUpperBound, upperString, localId);
1241 }
1242 catch (const data_error&)
1243 {
1244 //Presumably, the lower bound was higher than the *default* upper
1245 // bound--try them in the reverse order instead.
1246 TiXmlElement * upperBoundElement =
1247 singleRequiredChild(priorElement,xmlstr::XML_TAG_PRIORUPPERBOUND);
1248 string upperString = getNodeText(upperBoundElement);
1249 uiInterface.doSet(uistr::priorUpperBound, upperString, localId);
1250
1251 TiXmlElement * lowerBoundElement =
1252 singleRequiredChild(priorElement,xmlstr::XML_TAG_PRIORLOWERBOUND);
1253 string lowerString = getNodeText(lowerBoundElement);
1254 uiInterface.doSet(uistr::priorLowerBound, lowerString, localId);
1255
1256 }
1257
1258 #ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
1259 TiXmlElement * sampleRateElement =
1260 singleOptionalChild(priorElement,xmlstr::XML_TAG_RELATIVE_SAMPLE_RATE);
1261 if (sampleRateElement)
1262 {
1263 string rateString = getNodeText(sampleRateElement);
1264 uiInterface.doSet(uistr::relativeSampleRate,rateString,localId);
1265 }
1266 #endif
1267 }
1268
DoTraits(TiXmlElement * dataElement)1269 void ParseTreeToSettings::DoTraits(TiXmlElement * dataElement)
1270 {
1271 TiXmlNode * child = NULL;
1272
1273 long regionNumber = 0;
1274
1275 while((child = dataElement->IterateChildren(xmlstr::XML_TAG_REGION,child)))
1276 {
1277 TiXmlElement * traitsElement = singleOptionalChild(child->ToElement(), xmlstr::XML_TAG_TRAITS);
1278
1279 DoRegionTraits(traitsElement,regionNumber);
1280 regionNumber++;
1281 }
1282 }
1283
DoRegionTraits(TiXmlElement * traitsElement,long regionNumber)1284 void ParseTreeToSettings::DoRegionTraits(TiXmlElement * traitsElement, long regionNumber)
1285 {
1286 if (traitsElement != NULL)
1287 {
1288 TiXmlNode * traitNode = NULL;
1289 while((traitNode = traitsElement->IterateChildren(xmlstr::XML_TAG_TRAIT,traitNode)))
1290 {
1291 TiXmlElement * traitElement = traitNode->ToElement();
1292 DoRegionTrait(traitElement, regionNumber);
1293 }
1294 }
1295 }
1296
DoRegionTrait(TiXmlElement * traitElement,long regionNumber)1297 void ParseTreeToSettings::DoRegionTrait(TiXmlElement * traitElement, long regionNumber)
1298 {
1299 TiXmlElement * nameElement = singleRequiredChild(traitElement, xmlstr::XML_TAG_NAME);
1300 string tname = getNodeText(nameElement);
1301 if (!uiInterface.GetCurrentVars().datapackplus.HasLocus(regionNumber, tname))
1302 {
1303 throw data_error("Error: trait settings found for '"
1304 + tname
1305 + "', but no data was found for this trait. Check spelling, or delete these settings.");
1306 }
1307 long locusNumber = uiInterface.GetCurrentVars().datapackplus.GetLocusIndex(regionNumber, tname);
1308
1309 TiXmlElement * modelElement = singleOptionalChild(traitElement, xmlstr::XML_TAG_MODEL);
1310 if (modelElement != NULL)
1311 {
1312 DoDLModel(modelElement,regionNumber,locusNumber);
1313 }
1314
1315 TiXmlElement * locationsElement = singleOptionalChild(traitElement, xmlstr::XML_TAG_POSSIBLE_LOCATIONS);
1316 if (locationsElement != NULL)
1317 {
1318 DoLocations(locationsElement, regionNumber, locusNumber);
1319 }
1320
1321 // Do this after setting the range (locations) so we can warn/throw appropriately
1322 TiXmlElement * analysisElement = singleOptionalChild(traitElement, xmlstr::XML_TAG_ANALYSIS);
1323 if (analysisElement != NULL)
1324 {
1325 DoAnalysis(analysisElement, regionNumber, locusNumber);
1326 }
1327
1328 TiXmlElement * phenotypesElement = singleOptionalChild(traitElement, xmlstr::XML_TAG_PHENOTYPES);
1329 if (phenotypesElement != NULL)
1330 {
1331 DoPhenotypes(phenotypesElement, regionNumber, locusNumber);
1332 }
1333 }
1334
DoLocations(TiXmlElement * locationsElement,long regionNumber,long locusNumber)1335 void ParseTreeToSettings::DoLocations(TiXmlElement* locationsElement, long regionNumber, long locusNumber)
1336 {
1337 UIRegId locId(regionNumber,locusNumber, uiInterface.GetCurrentVars());
1338 rangeset rset;
1339
1340 //Collect the whole range, then set it afterwards.
1341 TiXmlNode* rangeNode = NULL;
1342 while ((rangeNode = locationsElement->IterateChildren(xmlstr::XML_TAG_RANGE, rangeNode)))
1343 {
1344 TiXmlElement * rangeElement = rangeNode->ToElement();
1345 TiXmlElement * startElement = singleRequiredChild(rangeElement, xmlstr::XML_TAG_START);
1346 TiXmlElement * endElement = singleRequiredChild(rangeElement, xmlstr::XML_TAG_END);
1347 long start = ProduceLongOrBarf(getNodeText(startElement));
1348 long end = ProduceLongOrBarf(getNodeText(endElement));
1349 rangepair range = std::make_pair(start, end);
1350 range.second++;
1351 rset = AddPairToRange(range, rset);
1352 }
1353 uiInterface.GetCurrentVars().traitmodels.SetRange(locId, rset);
1354 }
1355
DoAnalysis(TiXmlElement * analysisElement,long regionNumber,long locusNumber)1356 void ParseTreeToSettings::DoAnalysis(TiXmlElement* analysisElement, long regionNumber, long locusNumber)
1357 {
1358 UIRegId locId(regionNumber,locusNumber, uiInterface.GetCurrentVars());
1359 string analysisString = getNodeText(analysisElement);
1360 mloc_type analysis = ProduceMlocTypeOrBarf(analysisString);
1361 uiInterface.GetCurrentVars().traitmodels.SetAnalysisType(locId, analysis);
1362 }
1363
DoPhenotypes(TiXmlElement * phenotypesElement,long regionNumber,long locusNumber)1364 void ParseTreeToSettings::DoPhenotypes(TiXmlElement* phenotypesElement, long regionNumber, long locusNumber)
1365 {
1366 StringVec2d alleles = uiInterface.GetCurrentVars().datapackplus.GetUniqueAlleles(regionNumber, locusNumber);
1367 size_t nalleles = alleles[0].size(); //Assumes only one marker
1368
1369 std::set<long> ploidies = uiInterface.GetCurrentVars().datapackplus.GetPloidies(regionNumber);
1370 //Now calculate how many genotypes we need. For a given number
1371 // of alleles A, and ploidy P, we need:
1372 //
1373 // (A + P - 1)!
1374 // ------------
1375 // P!(A - 1)!
1376 //
1377 // which is the formula for 'combination with repetition'. See, for example
1378 // http://en.wikipedia.org/wiki/Combinatorics
1379 //
1380 // That's gotta be some sort of milestone--a URL in a comment! Wow.
1381 //
1382 // At any rate, we need to calculate this formula for each possible ploidy,
1383 // and add them up. Normally, there'd be only one, but there are exceptions
1384 // like if our region is on the X chromosome.
1385
1386 size_t targetCombs = 0;
1387 for (std::set<long>::iterator ploid = ploidies.begin(); ploid != ploidies.end(); ploid++)
1388 {
1389 targetCombs += static_cast<long>(factorial(nalleles + (*ploid) - 1)) /
1390 (static_cast<long>(factorial(*ploid)) * static_cast<long>(factorial(nalleles - 1)));
1391 }
1392
1393 vector<TiXmlElement *> genotypeElements = getAllOptionalDescendantElements(phenotypesElement, xmlstr::XML_TAG_GENOTYPE);
1394 if (genotypeElements.size() < targetCombs)
1395 {
1396 string msg = "Error: not enough genotypes listed in your phenotype list. You must provide a set of phenotypes"
1397 " for every possible combination of alleles (for which you have "
1398 + ToString(nalleles) + ")";
1399 if (ploidies.size() > 1)
1400 {
1401 msg += ". Also note that you have "
1402 + ToString(ploidies.size())
1403 + " different ploidies in your list of individuals, meaning that you must provide a complete"
1404 " list of combinations of alleles for all different possible numbers of genes, too.";
1405 }
1406 else
1407 {
1408 msg += ", for the ploidy of your individuals (which are all " + ToString(*ploidies.begin()) + ".)";
1409 }
1410 msg += " All told, you should have "
1411 + ToString(targetCombs)
1412 + " genotypes, while you currently have "
1413 + ToString(genotypeElements.size()) + ".";
1414
1415 throw data_error(msg);
1416 }
1417 if (genotypeElements.size() > targetCombs)
1418 {
1419 string msg = "Error: You have provided too many genotypes for the number of alleles for your trait ("
1420 + ToString(nalleles)
1421 + ") for the ploidies of your sampled individuals. You should have "
1422 + ToString(targetCombs)
1423 + ", but currently have "
1424 + ToString(genotypeElements.size())
1425 + ". Remember that order doesn't matter for genotypes--being heterozygous should give you the same phenotype regardless of which allele is which.";
1426
1427 throw data_error(msg);
1428 }
1429
1430 for (vector<TiXmlElement* >::iterator genotypeElement = genotypeElements.begin(); genotypeElement != genotypeElements.end(); genotypeElement++)
1431 {
1432 //Get the set of alleles.
1433 TiXmlElement* alleleElement = singleRequiredChild(*genotypeElement, xmlstr::XML_TAG_ALLELES);
1434 StringVec1d alleles = getNodeTextSplitOnWhitespace(alleleElement);
1435
1436 //Get the (potentially many) associated phenotypes.
1437 TiXmlNode * phenotypeNode = NULL;
1438 while ((phenotypeNode = (*genotypeElement)->IterateChildren(xmlstr::XML_TAG_PHENOTYPE, phenotypeNode)))
1439 {
1440 TiXmlElement* phenotypeElement = phenotypeNode->ToElement();
1441 TiXmlElement* phenNameElement = singleRequiredChild(phenotypeElement, xmlstr::XML_TAG_PHENOTYPE_NAME);
1442 TiXmlElement* penetranceElement = singleRequiredChild(phenotypeElement, xmlstr::XML_TAG_PENETRANCE);
1443
1444 string phenName = getNodeText(phenNameElement);
1445 double penetrance = ProduceDoubleOrBarf(getNodeText(penetranceElement));
1446 UIRegId locId(regionNumber,locusNumber, uiInterface.GetCurrentVars());
1447 uiInterface.GetCurrentVars().traitmodels.AddPhenotype(locId, alleles, phenName, penetrance);
1448 }
1449 }
1450 }
1451
1452 //____________________________________________________________________________________
1453