1 // $Id: gc_structures.cpp,v 1.65 2012/05/02 21:39:13 jmcgill Exp $
2
3 /*
4 Copyright 2002 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 "gc_creation_info.h"
14 #include "gc_default.h"
15 #include "gc_datastore.h"
16 #include "gc_errhandling.h"
17 #include "gc_file.h"
18 #include "gc_genotype_resolution.h"
19 #include "gc_loci_match.h"
20 #include "gc_locus_err.h"
21 #include "gc_parse.h"
22 #include "gc_parse_block.h"
23 #include "gc_pop_match.h"
24 #include "gc_strings.h"
25 #include "gc_strings_infile.h"
26 #include "gc_strings_structures.h"
27 #include "gc_structures.h"
28 #include "gc_structures_err.h"
29 #include "gc_trait_err.h"
30 #include "wx/log.h"
31
32 //------------------------------------------------------------------------------------
33
34 const wxString &
infoString() const35 gcNameSet::infoString() const
36 {
37 return gcstr_structures::objDict;
38 }
39
gcNameSet()40 gcNameSet::gcNameSet()
41 {
42 }
43
~gcNameSet()44 gcNameSet::~gcNameSet()
45 {
46 }
47
48 //------------------------------------------------------------------------------------
49
50 const wxString &
infoString() const51 gcTraitNameSet::infoString() const
52 {
53 return gcstr_structures::traitDict;
54 }
55
gcTraitNameSet()56 gcTraitNameSet::gcTraitNameSet()
57 {
58 }
59
~gcTraitNameSet()60 gcTraitNameSet::~gcTraitNameSet()
61 {
62 }
63
64 //------------------------------------------------------------------------------------
65
GCStructures(const GCDataStore * dsp)66 GCStructures::GCStructures(const GCDataStore * dsp)
67 :
68 m_dataStoreP(dsp),
69 m_divergenceState(false),
70 //m_divMigMatrixDefined(false),
71 //m_migMatrixDefined(false)
72 m_panelsState(false)
73 {
74 }
75
~GCStructures()76 GCStructures::~GCStructures()
77 {
78 }
79
80 gcIdSet
GetBlockIds(size_t popId,size_t locusId) const81 GCStructures::GetBlockIds(size_t popId, size_t locusId) const
82 {
83 gcIdSet outSet;
84 gcPopLocusIdPair p(popId,locusId);
85 gcBlockSetMap::const_iterator iter = m_blocks.find(p);
86 if(iter != m_blocks.end())
87 {
88 outSet = (*iter).second;
89 }
90 return outSet;
91 }
92
93 gcIdSet
GetBlocksForLocus(size_t locusId) const94 GCStructures::GetBlocksForLocus(size_t locusId) const
95 {
96 gcIdSet blockIds;
97
98 for(gcPopMap::const_iterator i=m_pops.begin();
99 i != m_pops.end(); i++)
100 {
101 size_t popId = (*i).first;
102 gcPopLocusIdPair p(popId,locusId);
103 gcBlockSetMap::const_iterator biter = m_blocks.find(p);
104 if(biter != m_blocks.end())
105 {
106 const gcIdSet & blocks = (*biter).second;
107 if(!blocks.empty())
108 {
109 blockIds.insert(blocks.begin(),blocks.end());
110 }
111 }
112 }
113 return blockIds;
114 }
115
116 bool
HaveBlocksForRegion(size_t regionId) const117 GCStructures::HaveBlocksForRegion(size_t regionId) const
118 {
119 gcIdVec locusIds = GetLocusIdsForRegionByCreation(regionId);
120
121 if(locusIds.size() == 0)
122 {
123 return false;
124 }
125 else
126 {
127 for(gcIdVec::const_iterator liter=locusIds.begin(); liter != locusIds.end(); liter++)
128 {
129 size_t locusId = (*liter);
130 if(HaveBlocksForLocus(locusId))
131 {
132 return true;
133 }
134 }
135 }
136 return false;
137 }
138
139 bool
HaveBlocksForLocus(size_t locusId) const140 GCStructures::HaveBlocksForLocus(size_t locusId) const
141 {
142 for(gcPopMap::const_iterator i=m_pops.begin();
143 i != m_pops.end(); i++)
144 {
145 size_t popId = (*i).first;
146 gcPopLocusIdPair p(popId,locusId);
147 gcBlockSetMap::const_iterator biter = m_blocks.find(p);
148 if(biter != m_blocks.end())
149 {
150 const gcIdSet & blocks = (*biter).second;
151 if(!blocks.empty())
152 {
153 return true;
154 }
155 }
156 }
157 return false;
158 }
159
160 bool
HaveBlocksForPop(size_t popId) const161 GCStructures::HaveBlocksForPop(size_t popId) const
162 {
163 // EWFIX.P5 SPEED -- could make faster with additional structures
164 for(gcLocusMap::const_iterator i = m_loci.begin(); i != m_loci.end(); i++)
165 {
166 size_t locusId = (*i).first;
167 gcPopLocusIdPair p(popId,locusId);
168 gcBlockSetMap::const_iterator biter = m_blocks.find(p);
169 if(biter != m_blocks.end())
170 {
171 const gcIdSet & blocks = (*biter).second;
172 if(!blocks.empty())
173 {
174 return true;
175 }
176 }
177 }
178 return false;
179 }
180
181 bool
IsBlessedRegion(size_t regionId) const182 GCStructures::IsBlessedRegion(size_t regionId) const
183 {
184 const gcRegion & region = GetRegion(regionId);
185 if(region.GetBlessed()) return true;
186 gcIdVec loci = GetLocusIdsForRegionByCreation(regionId);
187 for(gcIdVec::const_iterator iter=loci.begin(); iter!=loci.end();iter++)
188 {
189 const gcLocus & locus = GetLocus(*iter);
190 if(locus.GetBlessed()) return true;
191 }
192 return false;
193 }
194
195 bool
IsBlessedLocus(size_t locusId) const196 GCStructures::IsBlessedLocus(size_t locusId) const
197 {
198 return GetLocus(locusId).GetBlessed();
199 }
200
201 bool
IsBlessedPop(size_t popId) const202 GCStructures::IsBlessedPop(size_t popId) const
203 {
204 return GetPop(popId).GetBlessed();
205 }
206
207 bool
GetDivergenceState() const208 GCStructures::GetDivergenceState() const
209 {
210 return m_divergenceState;
211 }
212
213 bool
GetDivergenceState()214 GCStructures::GetDivergenceState()
215 {
216 return m_divergenceState;
217 }
218
219 void
SetDivergenceState(bool state)220 GCStructures::SetDivergenceState(bool state)
221 {
222 m_divergenceState = state;
223 }
224
225 bool
GetPanelsState() const226 GCStructures::GetPanelsState() const
227 {
228 return m_panelsState;
229 }
230
231 bool
GetPanelsState()232 GCStructures::GetPanelsState()
233 {
234 return m_panelsState;
235 }
236
237 void
SetPanelsState(bool state)238 GCStructures::SetPanelsState(bool state)
239 {
240 m_panelsState = state;
241 }
242
243
244 bool
GetDivMigMatrixDefined() const245 GCStructures::GetDivMigMatrixDefined() const
246 {
247 return m_divMigMatrixDefined;
248 }
249
250 bool
GetDivMigMatrixDefined()251 GCStructures::GetDivMigMatrixDefined()
252 {
253 return m_divMigMatrixDefined;
254 }
255
256 void
SetDivMigMatrixDefined(bool state)257 GCStructures::SetDivMigMatrixDefined(bool state)
258 {
259 m_divMigMatrixDefined = state;
260 }
261
262 bool
GetMigMatrixDefined() const263 GCStructures::GetMigMatrixDefined() const
264 {
265 return m_migMatrixDefined;
266 }
267
268 bool
GetMigMatrixDefined()269 GCStructures::GetMigMatrixDefined()
270 {
271 return m_migMatrixDefined;
272 }
273
274 void
SetMigMatrixDefined(bool state)275 GCStructures::SetMigMatrixDefined(bool state)
276 {
277 m_migMatrixDefined = state;
278 }
279
280 constObjVector
GetConstDisplayableRegions() const281 GCStructures::GetConstDisplayableRegions() const
282 {
283 constObjVector toReturn;
284 for(gcDisplayOrder::const_iterator i=m_regionDisplay.begin();
285 i != m_regionDisplay.end();
286 i++)
287 {
288 size_t regionId = *i;
289 if(HaveBlocksForRegion(regionId) || IsBlessedRegion(regionId))
290 {
291 toReturn.push_back(&(GetRegion(regionId)));
292 }
293 }
294 return toReturn;
295 }
296
297 objVector
GetDisplayableRegions()298 GCStructures::GetDisplayableRegions()
299 {
300 objVector toReturn;
301 for(gcDisplayOrder::const_iterator i=m_regionDisplay.begin();
302 i != m_regionDisplay.end();
303 i++)
304 {
305 size_t regionId = *i;
306 if(HaveBlocksForRegion(regionId) || IsBlessedRegion(regionId))
307 {
308 toReturn.push_back(&(GetRegion(regionId)));
309 }
310 }
311 return toReturn;
312 }
313
314 gcDisplayOrder
GetDisplayableRegionIds() const315 GCStructures::GetDisplayableRegionIds() const
316 {
317 gcDisplayOrder toReturn;
318 for(gcDisplayOrder::const_iterator i=m_regionDisplay.begin();
319 i != m_regionDisplay.end();
320 i++)
321 {
322 size_t regionId = *i;
323 if(HaveBlocksForRegion(regionId) || IsBlessedRegion(regionId))
324 {
325 toReturn.push_back(regionId);
326 }
327 }
328 return toReturn;
329 }
330
331 constObjVector
GetConstDisplayableLoci() const332 GCStructures::GetConstDisplayableLoci() const
333 {
334 constObjVector regionsInOrder = GetConstDisplayableRegions();
335
336 constObjVector toReturn;
337 for(gcRegionMap::const_iterator giter=m_regions.begin(); giter != m_regions.end(); giter++)
338 {
339 const gcRegion & region = (*giter).second;
340 gcIdVec loci = GetLocusIdsForRegionByMapPosition(region.GetId());
341 for(gcIdVec::const_iterator liter=loci.begin(); liter != loci.end(); liter++)
342 {
343 size_t locusId = *liter;
344 if(HaveBlocksForLocus(locusId) || IsBlessedLocus(locusId))
345 {
346 toReturn.push_back(&(GetLocus(locusId)));
347 }
348 }
349 }
350 return toReturn;
351 }
352
353 objVector
GetDisplayableLoci()354 GCStructures::GetDisplayableLoci()
355 {
356 objVector regionsInOrder = GetDisplayableRegions();
357
358 objVector toReturn;
359 for(gcRegionMap::const_iterator giter=m_regions.begin(); giter != m_regions.end(); giter++)
360 {
361 const gcRegion & region = (*giter).second;
362 gcIdVec loci = GetLocusIdsForRegionByMapPosition(region.GetId());
363 for(gcIdVec::const_iterator liter=loci.begin(); liter != loci.end(); liter++)
364 {
365 size_t locusId = *liter;
366 if(HaveBlocksForLocus(locusId) || IsBlessedLocus(locusId))
367 {
368 toReturn.push_back(&(GetLocus(locusId)));
369 }
370 }
371 }
372 return toReturn;
373 }
374
375 gcDisplayOrder
GetDisplayableLocusIds() const376 GCStructures::GetDisplayableLocusIds() const
377 {
378 gcDisplayOrder toReturn;
379
380 for(gcRegionMap::const_iterator giter=m_regions.begin(); giter != m_regions.end(); giter++)
381 {
382 const gcRegion & region = (*giter).second;
383 gcIdVec loci = GetLocusIdsForRegionByMapPosition(region.GetId());
384 for(gcIdVec::const_iterator liter=loci.begin(); liter != loci.end(); liter++)
385 {
386 size_t locusId = *liter;
387 if(HaveBlocksForLocus(locusId) || IsBlessedLocus(locusId))
388 {
389 toReturn.push_back(locusId);
390 }
391 }
392 }
393 return toReturn;
394 }
395
396 #if 0 // Potentially DEAD CODE (bobgian, Feb 2010)
397 constObjVector
398 GCStructures::GetConstDisplayableLociFor(size_t regionId) const
399 {
400 constObjVector toReturn;
401 const gcRegion & region = GetRegion(regionId);
402 gcIdVec loci = GetLocusIdsdForRegionByMapPosition(region.GetId());
403 for(gcIdVec::const_iterator liter=loci.begin(); liter != loci.end(); liter++)
404 {
405 size_t locusId = *liter;
406 if(HaveBlocksForLocus(locusId) || IsBlessedLocus(locusId))
407 {
408 toReturn.push_back(&(GetLocus(locusId)));
409 }
410 }
411 return toReturn;
412 }
413 #endif
414
415 constObjVector
GetConstDisplayableLociInMapOrderFor(size_t regionId) const416 GCStructures::GetConstDisplayableLociInMapOrderFor(size_t regionId) const
417 {
418 constObjVector toReturn;
419 const gcRegion & region = GetRegion(regionId);
420 gcIdVec loci = GetLocusIdsForRegionByMapPosition(region.GetId());
421 for(gcIdVec::const_iterator liter=loci.begin(); liter != loci.end(); liter++)
422 {
423 size_t locusId = *liter;
424 if(HaveBlocksForLocus(locusId) || IsBlessedLocus(locusId))
425 {
426 toReturn.push_back(&(GetLocus(locusId)));
427 }
428 }
429 return toReturn;
430 }
431
432 constObjVector
GetConstDisplayableLinkedLociInMapOrderFor(size_t regionId) const433 GCStructures::GetConstDisplayableLinkedLociInMapOrderFor(size_t regionId) const
434 {
435 constObjVector toReturn;
436 const gcRegion & region = GetRegion(regionId);
437 gcIdVec loci = GetLocusIdsForRegionByMapPosition(region.GetId());
438 for(gcIdVec::const_iterator liter=loci.begin(); liter != loci.end(); liter++)
439 {
440 size_t locusId = *liter;
441 if(HaveBlocksForLocus(locusId) || IsBlessedLocus(locusId))
442 {
443 const gcLocus & locusRef = GetLocus(locusId);
444 if(locusRef.GetLinked())
445 {
446 toReturn.push_back(&locusRef);
447 }
448 }
449 }
450 return toReturn;
451 }
452
453 objVector
GetDisplayableLociFor(size_t regionId)454 GCStructures::GetDisplayableLociFor(size_t regionId)
455 {
456
457 objVector toReturn;
458 const gcRegion & region = GetRegion(regionId);
459 const GCLocusInfoMap & loci = region.GetLocusInfoMap();
460 for(GCLocusInfoMap::const_iterator liter=loci.begin(); liter != loci.end(); liter++)
461 {
462 size_t locusId = (*liter).first;
463 if(HaveBlocksForLocus(locusId) || IsBlessedLocus(locusId))
464 {
465 toReturn.push_back(&(GetLocus(locusId)));
466 }
467 }
468 return toReturn;
469 }
470
471 constObjVector
GetConstDisplayablePops() const472 GCStructures::GetConstDisplayablePops() const
473 {
474 constObjVector toReturn;
475 for(gcDisplayOrder::const_iterator i=m_popDisplay.begin();
476 i != m_popDisplay.end();
477 i++)
478 {
479 size_t popId = *i;
480 if(HaveBlocksForPop(popId) || IsBlessedPop(popId))
481 {
482 toReturn.push_back(&(GetPop(popId)));
483 }
484 }
485 return toReturn;
486 }
487
488 gcDisplayOrder
GetDisplayablePopIds() const489 GCStructures::GetDisplayablePopIds() const
490 {
491 gcDisplayOrder toReturn;
492 for(gcDisplayOrder::const_iterator i=m_popDisplay.begin();
493 i != m_popDisplay.end();
494 i++)
495 {
496 size_t popId = *i;
497 if(HaveBlocksForPop(popId) || IsBlessedPop(popId))
498 {
499 toReturn.push_back(popId);
500 }
501 }
502 return toReturn;
503 }
504
505 objVector
GetDisplayablePops()506 GCStructures::GetDisplayablePops()
507 {
508 objVector toReturn;
509 for(gcDisplayOrder::const_iterator i=m_popDisplay.begin();
510 i != m_popDisplay.end();
511 i++)
512 {
513 size_t popId = *i;
514 if(HaveBlocksForPop(popId) || IsBlessedPop(popId))
515 {
516 toReturn.push_back(&(GetPop(popId)));
517 }
518 }
519 return toReturn;
520 }
521
522 gcDisplayOrder
GetParentIds() const523 GCStructures::GetParentIds() const
524 {
525 gcDisplayOrder toReturn;
526 for(gcParentMap::const_iterator i=m_parents.begin();
527 i != m_parents.end();
528 i++)
529 {
530 size_t parentId = (*i).first;
531 toReturn.push_back(parentId);
532 }
533 return toReturn;
534 }
535
536 objVector
GetParents()537 GCStructures::GetParents()
538 {
539 objVector toReturn;
540 for(gcParentMap::const_iterator i=m_parents.begin();
541 i != m_parents.end();
542 i++)
543 {
544 size_t parentId = (*i).first;
545 toReturn.push_back(&(GetParent(parentId)));
546 }
547 return toReturn;
548 }
549
550 constObjVector
GetConstParents() const551 GCStructures::GetConstParents() const
552 {
553 constObjVector toReturn;
554 for(gcParentMap::const_iterator i=m_parents.begin();
555 i != m_parents.end();
556 i++)
557 {
558 size_t parentId = (*i).first;
559 toReturn.push_back(&(GetParent(parentId)));
560 }
561 return toReturn;
562 }
563
564 constObjVector
GetConstTraits() const565 GCStructures::GetConstTraits() const
566 {
567 constObjVector toReturn;
568 for(gcTraitMap::const_iterator i=m_traitClasses.begin();
569 i != m_traitClasses.end();
570 i++)
571 {
572 const gcTraitInfo & trait = (*i).second;
573 toReturn.push_back(&trait);
574 }
575 return toReturn;
576 }
577
578 gcTraitAllele &
GetAllele(size_t id)579 GCStructures::GetAllele(size_t id)
580 {
581 gcAlleleMap::iterator iter = m_alleles.find(id);
582 assert(iter != m_alleles.end());
583 return (*iter).second;
584 }
585
586 const gcTraitAllele &
GetAllele(size_t id) const587 GCStructures::GetAllele(size_t id) const
588 {
589 gcAlleleMap::const_iterator iter = m_alleles.find(id);
590 assert(iter != m_alleles.end());
591 return (*iter).second;
592 }
593
594 gcRegion &
GetRegion(size_t id)595 GCStructures::GetRegion(size_t id)
596 {
597 gcRegionMap::iterator iter = m_regions.find(id);
598 assert(iter != m_regions.end());
599 return (*iter).second;
600 }
601
602 const gcRegion &
GetRegion(size_t id) const603 GCStructures::GetRegion(size_t id) const
604 {
605 gcRegionMap::const_iterator iter = m_regions.find(id);
606 assert(iter != m_regions.end());
607 return (*iter).second;
608 }
609
610 gcLocus &
GetLocus(size_t id)611 GCStructures::GetLocus(size_t id)
612 {
613 gcLocusMap::iterator iter = m_loci.find(id);
614 assert(iter != m_loci.end());
615 return (*iter).second;
616 }
617
618 const gcLocus &
GetLocus(size_t id) const619 GCStructures::GetLocus(size_t id) const
620 {
621 gcLocusMap::const_iterator iter = m_loci.find(id);
622 assert(iter != m_loci.end());
623 return (*iter).second;
624 }
625
626 gcPanel &
GetPanel(size_t id)627 GCStructures::GetPanel(size_t id)
628 {
629 for (gcPanelMap::iterator i=m_panels.begin();i != m_panels.end();i++)
630 {
631 if((*i).second.GetId() == id)
632 {
633 return (*i).second;
634 }
635 }
636 wxString badId = wxString::Format("%i",(int)id);
637 throw missing_panel_id(badId);
638 }
639
640 const gcPanel &
GetPanel(size_t id) const641 GCStructures::GetPanel(size_t id) const
642 {
643 for (gcPanelMap::const_iterator i=m_panels.begin();i != m_panels.end();i++)
644 {
645 if((*i).second.GetId() == id)
646 {
647 return (*i).second;
648 }
649 }
650 wxString badId = wxString::Format("%i",(int)id);
651 throw missing_panel_id(badId);
652 }
653
654 gcParent &
GetParent(size_t id)655 GCStructures::GetParent(size_t id)
656 {
657 for (gcParentMap::iterator i=m_parents.begin();i != m_parents.end();i++)
658 {
659 if((*i).second.GetId() == id)
660 {
661 return (*i).second;
662 }
663 }
664 wxString badId = wxString::Format("%i",(int)id);
665 throw missing_parent_id(badId);
666 }
667
668 const gcParent &
GetParent(size_t id) const669 GCStructures::GetParent(size_t id) const
670 {
671 for (gcParentMap::const_iterator i=m_parents.begin();i != m_parents.end();i++)
672 {
673 if((*i).second.GetId() == id)
674 {
675 return (*i).second;
676 }
677 }
678 wxString badId = wxString::Format("%i",(int)id);
679 throw missing_parent_id(badId);
680 }
681
682 gcPhenotype &
GetPhenotype(size_t id)683 GCStructures::GetPhenotype(size_t id)
684 {
685 gcPhenoMap::iterator iter = m_phenotypes.find(id);
686 assert(iter != m_phenotypes.end());
687 return (*iter).second;
688 }
689
690 const gcPhenotype &
GetPhenotype(size_t id) const691 GCStructures::GetPhenotype(size_t id) const
692 {
693 gcPhenoMap::const_iterator iter = m_phenotypes.find(id);
694 assert(iter != m_phenotypes.end());
695 return (*iter).second;
696 }
697
698 gcPopulation &
GetPop(size_t id)699 GCStructures::GetPop(size_t id)
700 {
701 gcPopMap::iterator iter = m_pops.find(id);
702 assert(iter != m_pops.end());
703 return (*iter).second;
704 }
705
706 const gcPopulation &
GetPop(size_t id) const707 GCStructures::GetPop(size_t id) const
708 {
709 gcPopMap::const_iterator iter = m_pops.find(id);
710 assert(iter != m_pops.end());
711 return (*iter).second;
712 }
713
714 gcTraitInfo &
GetTrait(size_t id)715 GCStructures::GetTrait(size_t id)
716 {
717 assert(m_traitClasses.find(id) != m_traitClasses.end());
718 return m_traitClasses[id];
719 }
720
721 const gcTraitInfo &
GetTrait(size_t id) const722 GCStructures::GetTrait(size_t id) const
723 {
724 gcTraitMap::const_iterator iter = m_traitClasses.find(id);
725 assert(iter != m_traitClasses.end());
726 return (*iter).second;
727 }
728
729 void
AssignAllele(gcTraitAllele & allele,gcTraitInfo & trait)730 GCStructures::AssignAllele(gcTraitAllele & allele, gcTraitInfo & trait)
731 {
732 if(allele.HasTraitId())
733 {
734 size_t oldTraitId = allele.GetTraitId();
735 gcTraitInfo & oldTrait = GetTrait(oldTraitId);
736 oldTrait.RemoveAllele(allele);
737 }
738 allele.SetTraitId(trait.GetId());
739 trait.AddAllele(allele);
740 }
741
742 void
AssignLocus(size_t locusId,size_t regionId)743 GCStructures::AssignLocus(size_t locusId, size_t regionId)
744 {
745 gcLocus & locus = GetLocus(locusId);
746 gcRegion & newRegion = GetRegion(regionId);
747 return AssignLocus(locus,newRegion);
748 }
749
750 void
AssignLocus(gcLocus & locus,gcRegion & newRegion)751 GCStructures::AssignLocus(gcLocus & locus, gcRegion & newRegion)
752 {
753 if(locus.HasRegion())
754 {
755 size_t regionId = locus.GetRegionId();
756 gcRegion & regionRef = GetRegion(regionId);
757 regionRef.RemoveLocusId(locus.GetId());
758 }
759
760 locus.SetRegionId(newRegion.GetId());
761 newRegion.AddLocus(locus);
762 }
763
764 void
AssignPhenotype(gcPhenotype & phenotype,gcTraitInfo & trait)765 GCStructures::AssignPhenotype(gcPhenotype & phenotype, gcTraitInfo & trait)
766 {
767 if(phenotype.HasTraitId())
768 {
769 size_t oldTraitId = phenotype.GetTraitId();
770 gcTraitInfo & oldTrait = GetTrait(oldTraitId);
771 oldTrait.RemovePhenotype(phenotype);
772 }
773 phenotype.SetTraitId(trait.GetId());
774 trait.AddPhenotype(phenotype);
775 }
776
777 void
AssignTrait(size_t traitId,size_t regionId)778 GCStructures::AssignTrait(size_t traitId, size_t regionId)
779 {
780 gcTraitInfo & trait = GetTrait(traitId);
781 gcRegion & newRegion = GetRegion(regionId);
782 AssignTrait(trait,newRegion);
783 }
784
785 void
AssignTrait(gcTraitInfo & trait,gcRegion & newRegion)786 GCStructures::AssignTrait(gcTraitInfo & trait, gcRegion & newRegion)
787 {
788 if(trait.HasRegionId())
789 {
790 size_t oldRegionId = trait.GetRegionId();
791 gcRegion & oldRegion = GetRegion(oldRegionId);
792 oldRegion.RemoveTraitId(trait.GetId());
793 }
794 newRegion.AddTraitId(trait.GetId());
795 }
796
797 size_t
GetPopForBlock(size_t blockId) const798 GCStructures::GetPopForBlock(size_t blockId) const
799 {
800 for(gcBlockSetMap::const_iterator i=m_blocks.begin(); i != m_blocks.end(); i++)
801 {
802 gcPopLocusIdPair popLoc = (*i).first;
803 const gcIdSet & blocks = (*i).second;
804 if(blocks.find(blockId) != blocks.end())
805 {
806 return popLoc.first;
807 }
808 }
809 assert(false);
810 return gcdefault::badIndex;
811 }
812
813 size_t
GetLocusForBlock(size_t blockId) const814 GCStructures::GetLocusForBlock(size_t blockId) const
815 {
816 for(gcBlockSetMap::const_iterator i=m_blocks.begin(); i != m_blocks.end(); i++)
817 {
818 gcPopLocusIdPair popLoc = (*i).first;
819 const gcIdSet & blocks = (*i).second;
820 if(blocks.find(blockId) != blocks.end())
821 {
822 return popLoc.second;
823 }
824 }
825 assert(false);
826 return gcdefault::badIndex;
827 }
828
829 void
AssignBlockToPop(size_t blockId,size_t popId)830 GCStructures::AssignBlockToPop(size_t blockId, size_t popId)
831 {
832 size_t oldLocus = GetLocusForBlock(blockId);
833 AssignBlock(blockId,popId,oldLocus);
834 }
835
836 void
AssignBlockToLocus(size_t blockId,size_t locusId)837 GCStructures::AssignBlockToLocus(size_t blockId, size_t locusId)
838 {
839 size_t oldPop = GetPopForBlock(blockId);
840 AssignBlock(blockId,oldPop,locusId);
841 }
842
843 void
AssignBlock(size_t blockId,size_t popId,size_t locusId)844 GCStructures::AssignBlock(size_t blockId, size_t popId, size_t locusId)
845 {
846 RemoveBlockAssignment(blockId);
847
848 gcPopLocusIdPair plPair(popId,locusId);
849 gcBlockSetMap::iterator iter = m_blocks.find(plPair);
850 if(iter == m_blocks.end())
851 // we need to insert a new, empty set
852 {
853 m_blocks[plPair] = gcIdSet();
854 iter = m_blocks.find(plPair);
855 }
856 gcIdSet & blockSet = (*iter).second;
857 blockSet.insert(blockId);
858 }
859
860 void
RemoveBlockAssignment(size_t blockId)861 GCStructures::RemoveBlockAssignment(size_t blockId)
862 {
863 size_t foundCount = 0;
864 for(gcBlockSetMap::iterator i=m_blocks.begin(); i != m_blocks.end(); i++)
865 {
866 gcIdSet & blockSetRef = (*i).second;
867 gcIdSet::iterator j = blockSetRef.find(blockId);
868 if(j != blockSetRef.end())
869 {
870 foundCount++;
871 blockSetRef.erase(j);
872 }
873 }
874 assert(foundCount < 2);
875 }
876
877 gcTraitAllele &
FetchOrMakeAllele(gcTraitInfo & trait,wxString name)878 GCStructures::FetchOrMakeAllele(gcTraitInfo& trait,wxString name)
879 {
880 if(HasAllele(name))
881 {
882 return GetAllele(trait,name);
883 }
884 return MakeAllele(name);
885 }
886
887 gcLocus &
FetchOrMakeLocus(gcRegion & region,wxString name,const gcCreationInfo & createInfo)888 GCStructures::FetchOrMakeLocus(gcRegion & region, wxString name, const gcCreationInfo & createInfo)
889 {
890 if(HasLocus(name))
891 {
892 return GetLocus(region,name);
893 }
894 return MakeLocus(region,name,true, createInfo); // true == blessed
895 }
896
897 gcPopulation &
FetchOrMakePop(wxString name)898 GCStructures::FetchOrMakePop(wxString name)
899 {
900 if(HasPop(name))
901 {
902 return GetPop(name);
903 }
904 return MakePop(name,true); // true == blessed
905 }
906
907 gcRegion &
FetchOrMakeRegion(wxString name)908 GCStructures::FetchOrMakeRegion(wxString name)
909 {
910 wxLogVerbose("****FetchOrMakeRegion name: %s", name.c_str()); // JMDBG
911 if(HasRegion(name))
912 {
913 return GetRegion(name);
914 }
915 return MakeRegion(name,true); // true == blessed
916 }
917
918 gcTraitInfo &
FetchOrMakeTrait(wxString name)919 GCStructures::FetchOrMakeTrait(wxString name)
920 {
921 if(HasTrait(name))
922 {
923 return GetTrait(name);
924 }
925 return MakeTrait(name);
926 }
927
928 gcTraitAllele &
MakeAllele(wxString name)929 GCStructures::MakeAllele(wxString name)
930 {
931 wxString nameToUse = m_traitNames.ReserveOrMakeName(name,gcstr::allele);
932 gcTraitAllele newAllele;
933 newAllele.SetName(nameToUse);
934 size_t newId = newAllele.GetId();
935 // don't do any setting after this step --
936 // it won't propagate
937 m_alleles[newId] = newAllele;
938 return m_alleles[newId];
939 }
940
941 gcLocus &
MakeLocus(size_t regionId,wxString name,bool blessed,const gcCreationInfo & createInfo)942 GCStructures::MakeLocus(size_t regionId, wxString name, bool blessed, const gcCreationInfo & createInfo)
943 {
944 gcRegion & region = GetRegion(regionId);
945 return MakeLocus(region,name,blessed,createInfo);
946 }
947
948 gcLocus &
MakeLocus(gcRegion & region,wxString name,bool blessed,const gcCreationInfo & creationInfo)949 GCStructures::MakeLocus(gcRegion & region, wxString name, bool blessed, const gcCreationInfo & creationInfo)
950 {
951 wxString nameToUse = m_names.ReserveOrMakeName(name,gcstr::locus);
952 gcLocus newLocus;
953 newLocus.SetName(nameToUse);
954 newLocus.SetBlessed(blessed);
955 newLocus.SetCreationInfo(creationInfo);
956 size_t newId = newLocus.GetId();
957 // don't make any direct changes to newLocus after
958 // this step -- it won't propagate
959 m_loci[newId] = newLocus;
960 AssignLocus(m_loci[newId],region);
961 return m_loci[newId];
962 }
963
964 gcPopulation &
MakePop(wxString name,bool blessed)965 GCStructures::MakePop(wxString name,bool blessed)
966 {
967 wxString nameToUse = m_names.ReserveOrMakeName(name,gcstr::population);
968 gcPopulation newPop;
969 newPop.SetName(nameToUse);
970 newPop.SetBlessed(blessed);
971 size_t newId = newPop.GetId();
972 // don't do any setting after this step --
973 // it won't propagate
974 m_pops[newId] = newPop;
975 m_popDisplay.push_back(newId);
976 return m_pops[newId];
977 }
978
979 gcPhenotype &
MakePhenotype(wxString name)980 GCStructures::MakePhenotype(wxString name)
981 {
982 wxString nameToUse = m_traitNames.ReserveOrMakeName(name,gcstr::phenotype);
983 gcPhenotype newPhenotype;
984 newPhenotype.SetName(nameToUse);
985 size_t newId = newPhenotype.GetId();
986 // don't do any setting after this step --
987 // it won't propagate
988 m_phenotypes[newId] = newPhenotype;
989 return m_phenotypes[newId];
990 }
991
992 gcRegion &
MakeRegion(wxString name,bool blessed)993 GCStructures::MakeRegion(wxString name,bool blessed)
994 {
995 wxLogVerbose("****MakeRegion name: %s", name.c_str()); // JMDBG
996 wxString nameToUse = m_names.ReserveOrMakeName(name,gcstr::region);
997 gcRegion newRegion;
998 newRegion.SetName(nameToUse);
999 newRegion.SetBlessed(blessed);
1000 wxLogVerbose("****MakeRegion nameToUse: %s blessed: %i", nameToUse.c_str(), (int)blessed); // JMDBG
1001 size_t newId = newRegion.GetId();
1002
1003 // create panels for each region / population pair
1004 for(gcPopMap::const_iterator i=m_pops.begin(); i != m_pops.end(); i++)
1005 {
1006 size_t popId = (*i).first;
1007 wxString popname = (*i).second.GetName();
1008 wxString pname;
1009 pname.Printf(wxT("panel:%s:%s"), newRegion.GetName().c_str(), (*i).second.GetName().c_str());
1010 wxString panelNameToUse = m_names.ReserveOrMakeName(pname,gcstr::panel);
1011 wxLogVerbose("***MakeRegion create Panel name: %s", panelNameToUse.c_str()); // JMDBG
1012 gcPanel newPanel = MakePanel(panelNameToUse, blessed, newId, popId);
1013 }
1014
1015 // don't do any setting after this step --
1016 // it won't propagate
1017 m_regions[newId] = newRegion;
1018 m_regionDisplay.push_back(newId);
1019 return m_regions[newId];
1020 }
1021
1022 gcPanel &
CreatePanel(size_t regionId,size_t popId)1023 GCStructures::CreatePanel(size_t regionId, size_t popId )
1024 {
1025 wxString pname;
1026 pname.Printf(wxT("panel:%s:%s"), GetRegion(regionId).GetName().c_str(), GetPop(popId).GetName().c_str());
1027 wxString panelNameToUse = m_names.ReserveOrMakeName(pname,gcstr::panel);
1028 wxLogVerbose("***CreatePanel name: %s", panelNameToUse.c_str()); // JMDBG
1029 return MakePanel(panelNameToUse, false, regionId, popId);
1030 }
1031
1032 gcPanel &
MakePanel(wxString name,bool blessed,size_t regionId,size_t popId)1033 GCStructures::MakePanel(wxString name, bool blessed, size_t regionId, size_t popId )
1034 {
1035 gcPanel newPanel;
1036 newPanel.SetName(name);
1037 newPanel.SetBlessed(blessed);
1038 newPanel.SetRegionId(regionId);
1039 newPanel.SetPopId(popId);
1040 wxLogVerbose("****MakePanel name: %s", name.c_str()); // JMDBG
1041 size_t newId = newPanel.GetId();
1042
1043 gcRegionPopIdPair regPopPair = gcRegionPopIdPair(regionId, popId);
1044
1045 // don't do any setting after this step --
1046 // it won't propagate
1047 m_panels[regPopPair] = newPanel;
1048 ////m_panelDisplay.push_back(newId);
1049 return m_panels[regPopPair];
1050 }
1051
1052 gcParent &
MakeParent(wxString name)1053 GCStructures::MakeParent(wxString name)
1054 {
1055 wxString nameToUse = m_names.ReserveOrMakeName(name,gcstr::parent);
1056 gcParent newParent;
1057 newParent.SetName(nameToUse);
1058 newParent.SetDispLevel(0);
1059 wxLogVerbose("****MakeParent name: %s", name.c_str()); // JMDBG
1060 size_t newId = newParent.GetId();
1061 // don't do any setting after this step --
1062 // it won't propagate
1063 m_parents[newId] = newParent;
1064 return m_parents[newParent.GetId()];
1065 }
1066
1067 gcMigration &
MakeMigration(bool blessed,size_t fromId,size_t toId)1068 GCStructures::MakeMigration(bool blessed, size_t fromId, size_t toId )
1069 {
1070 gcMigration newMigration;
1071
1072 newMigration.m_blessed = true;
1073 newMigration.m_hasFrom = true;
1074 newMigration.m_fromId = fromId;
1075 newMigration.m_hasTo = true;
1076 newMigration.m_toId = toId;
1077 newMigration.m_startValue = 50.0;
1078 newMigration.m_method = migmethod_USER;
1079 newMigration.m_profile = migprofile_NONE;
1080 newMigration.m_constraint = migconstraint_UNCONSTRAINED;
1081
1082 newMigration.SetName(wxString::Format(_T("cell%i:%i "), (int)fromId, (int)toId)); //JMDBG
1083
1084 wxLogVerbose("****MakeMigration from: %i to: %i", (int)fromId, (int)toId ); // JMDBG
1085 size_t newId = newMigration.GetId();
1086
1087 gcMigrationPair fromToPair = gcMigrationPair(fromId, toId);
1088
1089 // don't do any setting after this step --
1090 // it won't propagate
1091 m_migrations[fromToPair] = newMigration;
1092 return m_migrations[fromToPair];
1093 }
1094
1095 gcTraitInfo &
MakeTrait(wxString name)1096 GCStructures::MakeTrait(wxString name)
1097 {
1098 wxString nameToUse = m_traitNames.ReserveOrMakeName(name,gcstr::trait);
1099 gcTraitInfo newTrait;
1100 newTrait.SetName(nameToUse);
1101 // don't do any setting after this step --
1102 // it won't propagate
1103 assert(m_traitClasses.find(newTrait.GetId()) == m_traitClasses.end());
1104 m_traitClasses[newTrait.GetId()] = newTrait;
1105 return m_traitClasses[newTrait.GetId()];
1106 }
1107
1108 void
Rename(GCQuantum & object,wxString newName)1109 GCStructures::Rename(GCQuantum & object, wxString newName)
1110 {
1111 wxString oldName = object.GetName();
1112 wxLogVerbose("old name: %s \nnew name: %s", oldName.c_str(), newName.c_str()); // JMDBG
1113 if(newName.Cmp(oldName) == 0) return;
1114
1115 assert(m_names.HasName(oldName) || m_traitNames.HasName(oldName));
1116 assert(! (m_names.HasName(oldName) && m_traitNames.HasName(oldName)));
1117
1118 wxLogVerbose("checking m_names"); // JMDBG
1119 if(m_names.HasName(oldName))
1120 {
1121 wxLogVerbose("in m_names"); // JMDBG
1122 m_names.ReserveName(newName);
1123 object.SetName(newName);
1124 m_names.FreeName(oldName);
1125 }
1126
1127 wxLogVerbose("checking m_traitNames"); // JMDBG
1128 if(m_traitNames.HasName(oldName))
1129 {
1130 wxLogVerbose("in m_traitNames"); // JMDBG
1131 m_traitNames.ReserveName(newName);
1132 object.SetName(newName);
1133 m_traitNames.FreeName(oldName);
1134 }
1135 }
1136
1137 void
DebugDump(wxString prefix) const1138 GCStructures::DebugDump(wxString prefix) const
1139 {
1140 wxLogDebug(gcstr::structureDump,prefix.c_str()); // EWDUMPOK
1141 m_pops.DebugDump((prefix+gcstr::indent).c_str());
1142 m_regions.DebugDump((prefix+gcstr::indent).c_str());
1143 m_loci.DebugDump((prefix+gcstr::indent).c_str());
1144 m_traitClasses.DebugDump((prefix+gcstr::indent).c_str());
1145 m_blocks.DebugDump((prefix+gcstr::indent).c_str());
1146 m_files.DebugDump((prefix+gcstr::indent).c_str());
1147 }
1148
1149 gcRegion &
GetRegion(wxString name)1150 GCStructures::GetRegion(wxString name)
1151 {
1152 for(gcRegionMap::iterator i=m_regions.begin();
1153 i != m_regions.end();
1154 i++)
1155 {
1156 if((*i).second.GetName().CmpNoCase(name) == 0)
1157 {
1158 return (*i).second;
1159 }
1160 }
1161 missing_region e(wxString::Format(gcerr::missingRegion,name.c_str()).c_str());
1162 throw e;
1163 }
1164
1165 const gcRegion &
GetRegion(wxString name) const1166 GCStructures::GetRegion(wxString name) const
1167 {
1168 for(gcRegionMap::const_iterator i=m_regions.begin();
1169 i != m_regions.end();
1170 i++)
1171 {
1172 if((*i).second.GetName().CmpNoCase(name) == 0)
1173 {
1174 return (*i).second;
1175 }
1176 }
1177 missing_region e(wxString::Format(gcerr::missingRegion,name.c_str()).c_str());
1178 throw e;
1179 }
1180
1181 bool
HasRegion(wxString name) const1182 GCStructures::HasRegion(wxString name) const
1183 {
1184 for(gcRegionMap::const_iterator i=m_regions.begin();
1185 i != m_regions.end();
1186 i++)
1187 {
1188 if((*i).second.GetName().CmpNoCase(name) == 0)
1189 {
1190 return true;
1191 }
1192 }
1193 return false;
1194 }
1195
1196 gcTraitAllele &
GetAllele(wxString name)1197 GCStructures::GetAllele(wxString name)
1198 {
1199 for(gcAlleleMap::iterator i=m_alleles.begin();
1200 i != m_alleles.end();
1201 i++)
1202 {
1203 if((*i).second.GetName().Cmp(name) == 0) // alleles case sensitive
1204 {
1205 return (*i).second;
1206 }
1207 }
1208 throw gc_missing_allele(name);
1209 }
1210
1211 const gcTraitAllele &
GetAllele(wxString name) const1212 GCStructures::GetAllele(wxString name) const
1213 {
1214 for(gcAlleleMap::const_iterator i=m_alleles.begin();
1215 i != m_alleles.end();
1216 i++)
1217 {
1218 if((*i).second.GetName().Cmp(name) == 0) // alleles case sensitive
1219 {
1220 return (*i).second;
1221 }
1222 }
1223 throw gc_missing_allele(name);
1224 }
1225
1226 gcTraitAllele &
GetAllele(gcTraitInfo & trait,wxString alleleName)1227 GCStructures::GetAllele(gcTraitInfo& trait, wxString alleleName)
1228 {
1229 gcTraitAllele & allele = GetAllele(alleleName);
1230 assert(allele.HasTraitId());
1231 if(allele.GetTraitId() != trait.GetId())
1232 {
1233 const gcTraitInfo & oldTrait = GetTrait(allele.GetTraitId());
1234 throw gc_allele_trait_mismatch(alleleName,trait.GetName(),oldTrait.GetName());
1235 }
1236 return allele;
1237 }
1238
1239 bool
HasAllele(wxString name) const1240 GCStructures::HasAllele(wxString name) const
1241 {
1242 for(gcAlleleMap::const_iterator i=m_alleles.begin();
1243 i != m_alleles.end();
1244 i++)
1245 {
1246 if((*i).second.GetName().Cmp(name) == 0) // alleles case sensitive
1247 {
1248 return true;
1249 }
1250 }
1251 return false;
1252 }
1253
1254 gcLocus &
GetLocus(wxString name)1255 GCStructures::GetLocus(wxString name)
1256 {
1257 for(gcLocusMap::iterator i=m_loci.begin();
1258 i != m_loci.end();
1259 i++)
1260 {
1261 if((*i).second.GetName().CmpNoCase(name) == 0)
1262 {
1263 return (*i).second;
1264 }
1265 }
1266 throw gc_missing_locus(name);
1267 }
1268
1269 const gcLocus &
GetLocus(wxString name) const1270 GCStructures::GetLocus(wxString name) const
1271 {
1272 for(gcLocusMap::const_iterator i=m_loci.begin();
1273 i != m_loci.end();
1274 i++)
1275 {
1276 if((*i).second.GetName().CmpNoCase(name) == 0)
1277 {
1278 return (*i).second;
1279 }
1280 }
1281 throw gc_missing_locus(name);
1282 }
1283
1284 gcLocus &
GetLocus(gcRegion & region,wxString locusName)1285 GCStructures::GetLocus(gcRegion& region, wxString locusName)
1286 {
1287 gcLocus & locus = GetLocus(locusName);
1288 assert(locus.HasRegion());
1289 if(locus.GetRegionId() != region.GetId())
1290 {
1291 const gcRegion & oldRegion = GetRegion(locus.GetRegionId());
1292 throw gc_locus_region_mismatch(locusName,region.GetName(),oldRegion.GetName());
1293 }
1294 return locus;
1295 }
1296
1297 gcPanel &
GetPanel(wxString name)1298 GCStructures::GetPanel(wxString name)
1299 {
1300 for(gcPanelMap::iterator i=m_panels.begin();i != m_panels.end();i++)
1301 {
1302 if((*i).second.GetName().CmpNoCase(name) == 0)
1303 {
1304 return (*i).second;
1305 }
1306 }
1307 throw missing_panel(name);
1308 }
1309
1310 gcMigration &
GetMigration(size_t id)1311 GCStructures::GetMigration(size_t id)
1312 {
1313 for (gcMigrationMap::iterator i=m_migrations.begin();i != m_migrations.end();i++)
1314 {
1315 if((*i).second.GetId() == id)
1316 {
1317 return (*i).second;
1318 }
1319 }
1320 wxString badId = wxString::Format("%i",(int)id);
1321 throw missing_migration_id(badId);
1322 }
1323
1324 const gcMigration &
GetMigration(size_t id) const1325 GCStructures::GetMigration(size_t id) const
1326 {
1327 for (gcMigrationMap::const_iterator i=m_migrations.begin();i != m_migrations.end();i++)
1328 {
1329 if((*i).second.GetId() == id)
1330 {
1331 return (*i).second;
1332 }
1333 }
1334 wxString badId = wxString::Format("%i",(int)id);
1335 throw missing_migration_id(badId);
1336 }
1337
1338 bool
HasMigration(size_t fromId,size_t toId) const1339 GCStructures::HasMigration(size_t fromId, size_t toId) const
1340 {
1341 gcMigrationPair p(fromId, toId);
1342 gcMigrationMap::const_iterator iter = m_migrations.find(p);
1343 if(iter != m_migrations.end())
1344 {
1345 return true;
1346 }
1347 return false;
1348 }
1349
1350 gcMigration &
GetMigration(size_t fromId,size_t toId)1351 GCStructures::GetMigration(size_t fromId, size_t toId)
1352 {
1353 gcMigrationPair p(fromId, toId);
1354 gcMigrationMap::iterator iter = m_migrations.find(p);
1355 if(iter != m_migrations.end())
1356 {
1357 return (*iter).second;
1358 }
1359
1360 wxString fromName;
1361 if (IsPop(fromId))
1362 {
1363 fromName = GetPop(fromId).GetName();
1364 }
1365 else
1366 {
1367 if(IsParent(fromId))
1368 {
1369 fromName = GetParent(fromId).GetName();
1370 }
1371 else
1372 {
1373 fromName = wxString("Bad from ID:%s", fromId);
1374 }
1375 }
1376
1377 wxString toName;
1378 if (IsPop(toId))
1379 {
1380 toName = GetPop(toId).GetName();
1381 }
1382 else
1383 {
1384 if(IsParent(toId))
1385 {
1386 toName = GetParent(toId).GetName();
1387 }
1388 else
1389 {
1390 toName = wxString("Bad to ID:%s", toId);
1391 }
1392 }
1393 throw missing_migration(fromName, toName);
1394 }
1395
1396 const gcMigration &
GetMigration(size_t fromId,size_t toId) const1397 GCStructures::GetMigration(size_t fromId, size_t toId) const
1398 {
1399 gcMigrationPair p(fromId, toId);
1400 gcMigrationMap::const_iterator iter = m_migrations.find(p);
1401 if(iter != m_migrations.end())
1402 {
1403 return (*iter).second;
1404 }
1405
1406 wxString fromName;
1407 if (IsPop(fromId))
1408 {
1409 fromName = GetPop(fromId).GetName();
1410 }
1411 else
1412 {
1413 if(IsParent(fromId))
1414 {
1415 fromName = GetParent(fromId).GetName();
1416 }
1417 else
1418 {
1419 fromName = wxString("Bad from ID:%s", fromId);
1420 }
1421 }
1422
1423 wxString toName;
1424 if (IsPop(toId))
1425 {
1426 toName = GetPop(toId).GetName();
1427 }
1428 else
1429 {
1430 if(IsParent(toId))
1431 {
1432 toName = GetParent(toId).GetName();
1433 }
1434 else
1435 {
1436 toName = wxString("Bad to ID:%s", toId);
1437 }
1438 }
1439 throw missing_migration(fromName, toName);
1440 }
1441
1442 const gcPanel &
GetPanel(wxString name) const1443 GCStructures::GetPanel(wxString name) const
1444 {
1445 for(gcPanelMap::const_iterator i=m_panels.begin();i != m_panels.end();i++)
1446 {
1447 if((*i).second.GetName().CmpNoCase(name) == 0)
1448 {
1449 return (*i).second;
1450 }
1451 }
1452 throw missing_panel(name);
1453 }
1454
1455 bool
HasPanel(size_t regionId,size_t popId) const1456 GCStructures::HasPanel(size_t regionId, size_t popId) const
1457 {
1458 gcRegionPopIdPair p(regionId, popId);
1459 gcPanelMap::const_iterator iter = m_panels.find(p);
1460 if(iter == m_panels.end())
1461 {
1462 return false;
1463 }
1464 return true;
1465 }
1466
1467 gcPanel &
GetPanel(size_t regionId,size_t popId)1468 GCStructures::GetPanel(size_t regionId, size_t popId)
1469 {
1470 gcRegionPopIdPair p(regionId, popId);
1471 gcPanelMap::iterator iter = m_panels.find(p);
1472 if(iter != m_panels.end())
1473 {
1474 return (*iter).second;
1475 }
1476 wxString regionName = GetRegion(regionId).GetName();
1477 wxString popName = GetPop(popId).GetName();
1478 throw missing_panel(regionName, popName);
1479 }
1480
1481 const gcPanel &
GetPanel(size_t regionId,size_t popId) const1482 GCStructures::GetPanel(size_t regionId, size_t popId) const
1483 {
1484 gcRegionPopIdPair p(regionId, popId);
1485 gcPanelMap::const_iterator iter = m_panels.find(p);
1486 if(iter != m_panels.end())
1487 {
1488 return (*iter).second;
1489 }
1490 wxString regionName = GetRegion(regionId).GetName();
1491 wxString popName = GetPop(popId).GetName();
1492 throw missing_panel(regionName, popName);
1493 }
1494
1495 gcParent &
GetParent(wxString name)1496 GCStructures::GetParent(wxString name)
1497 {
1498 for(gcParentMap::iterator i=m_parents.begin();i != m_parents.end();i++)
1499 {
1500 if((*i).second.GetName().CmpNoCase(name) == 0)
1501 {
1502 return (*i).second;
1503 }
1504 }
1505 throw missing_parent(name);
1506 }
1507
1508 const gcParent &
GetParent(wxString name) const1509 GCStructures::GetParent(wxString name) const
1510 {
1511 for(gcParentMap::const_iterator i=m_parents.begin();i != m_parents.end();i++)
1512 {
1513 if((*i).second.GetName().CmpNoCase(name) == 0)
1514 {
1515 return (*i).second;
1516 }
1517 }
1518 throw missing_parent(name);
1519 }
1520
1521 gcPhenotype &
GetPhenotype(wxString name)1522 GCStructures::GetPhenotype(wxString name)
1523 {
1524 for(gcPhenoMap::iterator i=m_phenotypes.begin();
1525 i != m_phenotypes.end();
1526 i++)
1527 {
1528 if((*i).second.GetName().Cmp(name) == 0) // phenotypes case sensitive
1529 {
1530 return (*i).second;
1531 }
1532 }
1533 throw gc_missing_phenotype(name);
1534 }
1535
1536 const gcPhenotype &
GetPhenotype(wxString name) const1537 GCStructures::GetPhenotype(wxString name) const
1538 {
1539 for(gcPhenoMap::const_iterator i=m_phenotypes.begin();
1540 i != m_phenotypes.end();
1541 i++)
1542 {
1543 if((*i).second.GetName().Cmp(name) == 0) // phenotypes case sensitive
1544 {
1545 return (*i).second;
1546 }
1547 }
1548 throw gc_missing_phenotype(name);
1549 }
1550
1551 bool
HasLocus(wxString name) const1552 GCStructures::HasLocus(wxString name) const
1553 {
1554 for(gcLocusMap::const_iterator i=m_loci.begin();
1555 i != m_loci.end();
1556 i++)
1557 {
1558 if((*i).second.GetName().CmpNoCase(name) == 0)
1559 {
1560 return true;
1561 }
1562 }
1563 return false;
1564 }
1565
1566 gcPopulation &
GetPop(wxString name)1567 GCStructures::GetPop(wxString name)
1568 {
1569 for(gcPopMap::iterator i=m_pops.begin();
1570 i != m_pops.end();
1571 i++)
1572 {
1573 if((*i).second.GetName().CmpNoCase(name) == 0)
1574 {
1575 return (*i).second;
1576 }
1577 }
1578 throw gc_missing_population(name);
1579 }
1580
1581 const gcPopulation &
GetPop(wxString name) const1582 GCStructures::GetPop(wxString name) const
1583 {
1584 for(gcPopMap::const_iterator i=m_pops.begin();
1585 i != m_pops.end();
1586 i++)
1587 {
1588 if((*i).second.GetName().CmpNoCase(name) == 0)
1589 {
1590 return (*i).second;
1591 }
1592 }
1593 throw gc_missing_population(name);
1594 }
1595
1596 bool
HasPop(wxString name) const1597 GCStructures::HasPop(wxString name) const
1598 {
1599 for(gcPopMap::const_iterator i=m_pops.begin();
1600 i != m_pops.end();
1601 i++)
1602 {
1603 if((*i).second.GetName().CmpNoCase(name) == 0)
1604 {
1605 return true;
1606 }
1607 }
1608 return false;
1609 }
1610
1611 gcTraitInfo &
GetTrait(wxString name)1612 GCStructures::GetTrait(wxString name)
1613 {
1614 for(gcTraitMap::iterator i=m_traitClasses.begin();
1615 i != m_traitClasses.end();
1616 i++)
1617 {
1618 if((*i).second.GetName().CmpNoCase(name) == 0)
1619 {
1620 return (*i).second;
1621 }
1622 }
1623 missing_trait e(wxString::Format(gcerr::missingTrait,name.c_str()).c_str());
1624 throw e;
1625 }
1626
1627 const gcTraitInfo &
GetTrait(wxString name) const1628 GCStructures::GetTrait(wxString name) const
1629 {
1630 for(gcTraitMap::const_iterator i=m_traitClasses.begin();
1631 i != m_traitClasses.end();
1632 i++)
1633 {
1634 if((*i).second.GetName().CmpNoCase(name) == 0)
1635 {
1636 return (*i).second;
1637 }
1638 }
1639 missing_trait e(wxString::Format(gcerr::missingTrait,name.c_str()).c_str());
1640 throw e;
1641 }
1642
1643 bool
HasTrait(wxString name) const1644 GCStructures::HasTrait(wxString name) const
1645 {
1646 for(gcTraitMap::const_iterator i=m_traitClasses.begin();
1647 i != m_traitClasses.end();
1648 i++)
1649 {
1650 if((*i).second.GetName().CmpNoCase(name) == 0)
1651 {
1652 return true;
1653 }
1654 }
1655 return false;
1656 }
1657
1658 long
GetPopDisplayIndexOf(size_t popId) const1659 GCStructures::GetPopDisplayIndexOf(size_t popId) const
1660 {
1661 long popIndex = 0;
1662 for( gcDisplayOrder::const_iterator i = m_popDisplay.begin();
1663 i != m_popDisplay.end();
1664 i++, popIndex++)
1665 {
1666 if((*i) == popId) return popIndex;
1667 }
1668 return gcdefault::badDisplayIndex;
1669 }
1670
1671 long
GetLocusDisplayIndexOf(size_t locusId) const1672 GCStructures::GetLocusDisplayIndexOf(size_t locusId) const
1673 {
1674 long locusIndex = 0;
1675 for( gcDisplayOrder::const_iterator i = m_regionDisplay.begin();
1676 i != m_regionDisplay.end();
1677 i++)
1678 {
1679 size_t regionId = (*i);
1680 const gcRegion & region = GetRegion(regionId);
1681
1682 const GCLocusInfoMap & loci = region.GetLocusInfoMap();
1683
1684 if(loci.size() == 0)
1685 {
1686 locusIndex++;
1687 }
1688 else
1689 {
1690
1691 for(GCLocusInfoMap::const_iterator liter=loci.begin();
1692 liter != loci.end(); liter++)
1693 {
1694 if(locusId == ((*liter).first))
1695 {
1696 return locusIndex;
1697 }
1698 locusIndex++;
1699 }
1700 }
1701
1702 }
1703 return gcdefault::badDisplayIndex;
1704 }
1705
1706 bool
RegionHasAnyLinkedLoci(size_t regionId) const1707 GCStructures::RegionHasAnyLinkedLoci(size_t regionId) const
1708 {
1709 const gcRegion & region = GetRegion(regionId);
1710 const GCLocusInfoMap & loci = region.GetLocusInfoMap();
1711 for(GCLocusInfoMap::const_iterator i=loci.begin(); i != loci.end(); i++)
1712 {
1713 size_t locusId = (*i).first;
1714 const gcLocus & locus = GetLocus(locusId);
1715 if(locus.GetLinked())
1716 {
1717 return true;
1718 }
1719 }
1720 return false;
1721 }
1722
1723 bool
RegionHasAnyUnLinkedLoci(size_t regionId) const1724 GCStructures::RegionHasAnyUnLinkedLoci(size_t regionId) const
1725 {
1726 const gcRegion & region = GetRegion(regionId);
1727 const GCLocusInfoMap & loci = region.GetLocusInfoMap();
1728 for(GCLocusInfoMap::const_iterator i=loci.begin(); i != loci.end(); i++)
1729 {
1730 size_t locusId = (*i).first;
1731 const gcLocus & locus = GetLocus(locusId);
1732 if(!(locus.GetLinked()))
1733 {
1734 return true;
1735 }
1736 }
1737 return false;
1738 }
1739
1740 void
FragmentRegion(size_t regionId)1741 GCStructures::FragmentRegion(size_t regionId)
1742 {
1743 gcRegion & regionRef = GetRegion(regionId);
1744 gcIdVec loci = GetLocusIdsForRegionByCreation(regionId);
1745 for(gcIdVec::const_iterator i=loci.begin(); i != loci.end(); i++)
1746 {
1747 size_t locusId = *i;
1748 gcLocus & locusRef = GetLocus(locusId);
1749
1750 wxLogVerbose("****FragmentRegion "); // JMDBG
1751 gcRegion & newRegion = MakeRegion("",false); // false == not blessed
1752 if(regionRef.HasEffectivePopulationSize())
1753 {
1754 newRegion.SetEffectivePopulationSize(regionRef.GetEffectivePopulationSize());
1755 }
1756
1757 // propogate panel size out to all the new panels
1758 for(gcPopMap::const_iterator i=m_pops.begin(); i != m_pops.end(); i++)
1759 {
1760 size_t popId = (*i).first;
1761 if (HasBlock(locusId, popId))
1762 {
1763 GetPanel(newRegion.GetId(),popId).SetNumPanels(GetPanel(regionId,popId).GetNumPanels());
1764 GetPanel(newRegion.GetId(),popId).SetBlessed(true);
1765 }
1766 }
1767
1768 AssignLocus(locusRef,newRegion);
1769 }
1770
1771 // unbless old region
1772 regionRef.SetBlessed(false);
1773 }
1774
1775 void
LocusToOwnRegion(size_t locusId)1776 GCStructures::LocusToOwnRegion(size_t locusId)
1777 {
1778 gcLocus & locusRef = GetLocus(locusId);
1779
1780 wxLogVerbose("****LocusToOwnRegion "); // JMDBG
1781 gcRegion & newRegion = MakeRegion("",false); // false == not blessed
1782 if(locusRef.HasRegion())
1783 {
1784 gcRegion & regionRef = GetRegion(locusRef.GetRegionId());
1785 if(regionRef.HasEffectivePopulationSize())
1786 {
1787 newRegion.SetEffectivePopulationSize(regionRef.GetEffectivePopulationSize());
1788 }
1789 }
1790 AssignLocus(locusRef,newRegion);
1791 }
1792
1793 void
MergeRegions(gcIdVec regions)1794 GCStructures::MergeRegions(gcIdVec regions)
1795 {
1796 if(regions.size() < 2) return;
1797
1798 // testing mergeability
1799
1800 // population size set
1801 bool haveEffPopSize = false;
1802 double effPopSize = 1;
1803
1804 for(gcIdVec::iterator iter=regions.begin(); iter != regions.end(); iter++)
1805 {
1806 size_t regionId = *iter;
1807 gcRegion * regionP = &GetRegion(regionId);
1808
1809 // effective population size
1810 if(regionP->HasEffectivePopulationSize())
1811 {
1812 double pSize = regionP->GetEffectivePopulationSize();
1813 if(haveEffPopSize)
1814 {
1815 if(effPopSize != pSize)
1816 {
1817 throw effective_pop_size_clash(effPopSize,pSize);
1818 }
1819 }
1820 else
1821 {
1822 haveEffPopSize = true;
1823 effPopSize = pSize;
1824 }
1825 }
1826 }
1827
1828 // test panel size and blessed state
1829 gcDisplayOrder popIds = GetDisplayablePopIds();
1830
1831 for(gcDisplayOrder::iterator piter=popIds.begin(); piter !=popIds.end(); piter++)
1832 {
1833 size_t popId = *piter;
1834 gcPopulation * curPop = &GetPop(popId);
1835 gcRegion * region1;
1836 gcRegion * region2;
1837
1838 bool isBlessed = false;
1839 long panelSize = 0;
1840 //size_t lastId = 0;
1841
1842 for(gcIdVec::iterator iter=regions.begin(); iter != regions.end(); iter++)
1843 {
1844 size_t regionId = *iter;
1845 wxLogVerbose("In MergeRegions testing region: %i pop: %i", regionId, popId);
1846 if (HasPanel(regionId, popId))
1847 {
1848 const gcPanel &panelId = GetPanel(regionId, popId);
1849 long numP = panelId.GetNumPanels();
1850 wxLogVerbose("In MergeRegions region: %i pop: %i HasPanel numP: %i", regionId, popId, numP);
1851 if (panelId.GetBlessed())
1852 {
1853 if(!isBlessed)
1854 {
1855 // blessed overides unblessed, no matter what the value
1856 panelSize = numP;
1857 isBlessed = true;
1858 region1 = &GetRegion(regionId);
1859 }
1860 else
1861 {
1862 if (panelSize != numP)
1863 {
1864 wxLogVerbose("In MergeRegions panelSize: %i numP: %i", panelSize, numP);
1865 // can't merge panels that are both blessed and have different number of members
1866 region2 = &GetRegion(regionId);
1867 throw panel_size_clash(curPop->GetName(),region1->GetName(), region2->GetName());
1868 }
1869 }
1870 }
1871 else
1872 {
1873 if (numP > 0)
1874 {
1875 // this is a bug - you can't have >0 panel size that isn't blessed
1876 region1 = &GetRegion(regionId);
1877 throw panel_blessed_error(region1->GetName(), curPop->GetName());
1878 }
1879 }
1880 }
1881 }
1882 }
1883
1884 // everything worked so combine panels
1885 for(gcDisplayOrder::iterator piter=popIds.begin(); piter !=popIds.end(); piter++)
1886 {
1887 size_t popId = *piter;
1888
1889 bool isFirstRegion = true;
1890 gcPanel * firstPanel;
1891 gcPanel * curPanel;
1892 for(gcIdVec::iterator iter=regions.begin(); iter != regions.end(); iter++)
1893 {
1894 size_t regionId = *iter;
1895 wxLogVerbose("In MergeRegions - regionId: %i popId: %i", regionId, popId);
1896 if (isFirstRegion)
1897 {
1898 // everything propogates to the first region selected
1899 isFirstRegion = false;
1900 if (HasPanel(regionId, popId))
1901 {
1902 firstPanel = &GetPanel(regionId, popId);
1903 wxLogVerbose(" first panel regionId: %i popId: %i", regionId, popId);
1904 firstPanel->SetBlessed(true);
1905 }
1906 else
1907 {
1908 // corner case - first region did not contain this pop so create it
1909 CreatePanel(regionId, popId);
1910 }
1911 }
1912 else
1913 {
1914 // largest panel size wins
1915 if (HasPanel(regionId, popId))
1916 {
1917 curPanel = &GetPanel(regionId, popId);
1918 }
1919 if (curPanel->GetNumPanels() > firstPanel->GetNumPanels())
1920 {
1921 firstPanel->SetNumPanels(curPanel->GetNumPanels());
1922 }
1923 }
1924 }
1925 }
1926
1927 // combine regions
1928 bool firstItem = true;
1929 gcRegion * firstRegionPointer = NULL;
1930 for(gcIdVec::iterator iter=regions.begin(); iter != regions.end(); iter++)
1931 {
1932 size_t regionId = *iter;
1933 gcRegion * regionP = &GetRegion(regionId);
1934 assert(regionP != NULL);
1935 if(firstItem)
1936 {
1937 firstRegionPointer = regionP;
1938 }
1939 else
1940 {
1941 gcRegion & regionRef = *regionP;
1942 const GCLocusInfoMap & lmap = regionRef.GetLocusInfoMap();
1943 gcIdSet locusIds;
1944 for(GCLocusInfoMap::const_iterator miter=lmap.begin(); miter != lmap.end(); miter++)
1945 {
1946 locusIds.insert((*miter).first);
1947 }
1948 for(gcIdSet::const_iterator liter=locusIds.begin(); liter != locusIds.end(); liter++)
1949 {
1950 size_t locusId = *liter;
1951 AssignLocus(locusId,firstRegionPointer->GetId()); // EWFIX.P5 SPEED
1952 }
1953
1954 // unbless region that's now empty
1955 regionRef.SetBlessed(false);
1956 }
1957 firstItem = false;
1958 }
1959 if(haveEffPopSize)
1960 {
1961 firstRegionPointer->SetEffectivePopulationSize(effPopSize);
1962 }
1963 }
1964
1965 void
MergeLoci(gcIdVec locusIds)1966 GCStructures::MergeLoci(gcIdVec locusIds)
1967 // this operation can fail part way through. If that happens
1968 // in the batch converter, the thrown error will result in
1969 // aborting the program (cleanly, one would hope). If we're
1970 // in the gui, the thrown error should be caught in
1971 // gcUpdatingDialog::OnButton 's call to DoUpdateData
1972 {
1973 // find all blocks that are assigned to one of these loci and
1974 // assign them all to the first
1975 if(locusIds.size() < 2) return;
1976
1977 bool firstItem = true;
1978 gcLocus * firstLocusPointer = NULL;
1979 for(gcIdVec::iterator iter = locusIds.begin(); iter != locusIds.end(); iter++)
1980 {
1981 size_t locusId = *iter;
1982 gcLocus * locusP = &(GetLocus(locusId));
1983 assert(locusP != NULL);
1984 if(firstItem)
1985 {
1986 firstLocusPointer = locusP;
1987 }
1988 else
1989 {
1990
1991 firstLocusPointer->MergeWith(*locusP);
1992
1993 for(gcPopMap::const_iterator piter=m_pops.begin(); piter != m_pops.end(); piter++)
1994 {
1995 size_t popId = (*piter).first;
1996 gcPopLocusIdPair oldPair(popId,locusP->GetId());
1997 gcBlockSetMap::iterator biter = m_blocks.find(oldPair);
1998 if(biter != m_blocks.end())
1999 {
2000 gcIdSet blocks = (*biter).second;
2001 for(gcIdSet::iterator bsetIter=blocks.begin(); bsetIter != blocks.end(); bsetIter++)
2002 {
2003 size_t blockId = *bsetIter;
2004 AssignBlock(blockId,popId,firstLocusPointer->GetId()); // EWFIX.P5 SPEED
2005 }
2006 }
2007 }
2008 // unbless locus that's now empty
2009 locusP->SetBlessed(false);
2010 }
2011 firstItem = false;
2012 }
2013
2014 }
2015
2016 void
MergePops(gcIdVec popIds)2017 GCStructures::MergePops(gcIdVec popIds)
2018 {
2019 // find all blocks that are assigned to one of these pops and
2020 // assign them all to the first
2021 if(popIds.size() < 2) return;
2022
2023 bool firstItem = true;
2024 gcPopulation * firstPopP = NULL;
2025 for(gcIdVec::iterator iter = popIds.begin(); iter != popIds.end(); iter++)
2026 {
2027 size_t popId = *iter;
2028 gcPopulation * popP = &(GetPop(popId));
2029 assert(popP != NULL);
2030 if (popP->GetBlessed())
2031 {
2032 wxLogVerbose("****in GCStructures::MergePops popP: %i Blessed: true before merge", popId); // JMDBG
2033 }
2034 else
2035 {
2036 wxLogVerbose("****in GCStructures::MergePops popP: %i Blessed: false before merge", popId); // JMDBG
2037 }
2038
2039 if(firstItem)
2040 {
2041 firstPopP = popP;
2042 }
2043 else
2044 {
2045 for(gcLocusMap::const_iterator liter=m_loci.begin(); liter != m_loci.end(); liter++)
2046 {
2047 size_t locusId = (*liter).first;
2048 gcPopLocusIdPair oldPair(popP->GetId(),locusId);
2049 gcBlockSetMap::iterator biter = m_blocks.find(oldPair);
2050 if(biter != m_blocks.end())
2051 {
2052 gcIdSet blocks = (*biter).second;
2053 for(gcIdSet::iterator bsetIter=blocks.begin(); bsetIter != blocks.end(); bsetIter++)
2054 {
2055 size_t blockId = *bsetIter;
2056 AssignBlock(blockId,firstPopP->GetId(),locusId);
2057 }
2058 }
2059 }
2060 // unbless pop that's now empty
2061 popP->SetBlessed(false);
2062 }
2063 firstItem = false;
2064 if (popP->GetBlessed())
2065 {
2066 wxLogVerbose("****in GCStructures::MergePops popP: %i Blessed: true after merge", popId); // JMDBG
2067 }
2068 else
2069 {
2070 wxLogVerbose("****in GCStructures::MergePops popP: %i Blessed: false after merge", popId); // JMDBG
2071 }
2072 }
2073
2074 }
2075
2076 bool
HasBlock(size_t locusId,size_t popId) const2077 GCStructures::HasBlock(size_t locusId, size_t popId) const
2078 {
2079 gcPopLocusIdPair p(popId,locusId);
2080 gcBlockSetMap::const_iterator iter = m_blocks.find(p);
2081 if(iter != m_blocks.end())
2082 {
2083 return true;
2084 }
2085 return false;
2086 }
2087
2088 void
RemoveBlocks(gcIdSet blocksInFile)2089 GCStructures::RemoveBlocks(gcIdSet blocksInFile)
2090 {
2091
2092 // check in structures
2093 for(gcBlockSetMap::iterator iter=m_blocks.begin(); iter != m_blocks.end(); )
2094 {
2095 gcIdSet & bSet = (*iter).second;
2096 for(gcIdSet::iterator biter=bSet.begin(); biter != bSet.end(); )
2097 {
2098 size_t blockId = *biter;
2099 if(blocksInFile.find(blockId) != blocksInFile.end())
2100 {
2101 bSet.erase(biter++);
2102 }
2103 else
2104 {
2105 biter++;
2106 }
2107 }
2108 if(bSet.empty())
2109 {
2110 m_blocks.erase(iter++);
2111 }
2112 else
2113 {
2114 iter++;
2115 }
2116 }
2117 }
2118
2119 void
RemoveBlocksForLocus(size_t locusId)2120 GCStructures::RemoveBlocksForLocus(size_t locusId)
2121 {
2122 for(gcBlockSetMap::iterator iter=m_blocks.begin(); iter != m_blocks.end(); )
2123 {
2124 gcPopLocusIdPair p = (*iter).first;
2125 if(p.second == locusId)
2126 {
2127 m_blocks.erase(iter++);
2128 }
2129 else
2130 {
2131 iter++;
2132 }
2133 }
2134 }
2135
2136 void
RemoveBlocksForPop(size_t popId)2137 GCStructures::RemoveBlocksForPop(size_t popId)
2138 {
2139 for(gcBlockSetMap::iterator iter=m_blocks.begin(); iter != m_blocks.end(); )
2140 {
2141 gcPopLocusIdPair p = (*iter).first;
2142 if(p.first == popId)
2143 {
2144 m_blocks.erase(iter++);
2145 }
2146 else
2147 {
2148 iter++;
2149 }
2150 }
2151 }
2152
2153 gcIdVec
GetLocusIdsForRegionByCreation(size_t regionId) const2154 GCStructures::GetLocusIdsForRegionByCreation(size_t regionId) const
2155 {
2156 const gcRegion & region = GetRegion(regionId);
2157 const GCLocusInfoMap & loci = region.GetLocusInfoMap();
2158 gcIdVec vec;
2159 for(GCLocusInfoMap::const_iterator iter=loci.begin(); iter != loci.end(); iter++)
2160 {
2161 size_t locusId = (*iter).first;
2162 if(iter != loci.begin())
2163 {
2164 assert(locusId > vec.back());
2165 }
2166 vec.push_back(locusId);
2167 }
2168 return vec;
2169 }
2170
2171 gcIdVec
GetLocusIdsForRegionByMapPosition(size_t regionId) const2172 GCStructures::GetLocusIdsForRegionByMapPosition(size_t regionId) const
2173 {
2174 const gcRegion & region = GetRegion(regionId);
2175 const GCLocusInfoMap & loci = region.GetLocusInfoMap();
2176
2177 std::multimap<long,size_t> mapByPos;
2178
2179 for(GCLocusInfoMap::const_iterator iter=loci.begin(); iter != loci.end(); iter++)
2180 {
2181 size_t locusId = (*iter).first;
2182
2183 const gcLocus & locus = GetLocus(locusId);
2184 long pos = gcdefault::badMapPosition;
2185 if(locus.HasMapPosition())
2186 {
2187 pos = locus.GetMapPosition();
2188 }
2189 mapByPos.insert(std::pair<long, size_t>(pos,locusId));
2190 }
2191
2192 gcIdVec vec;
2193 for(std::multimap<long,size_t>::iterator i=mapByPos.begin(); i!=mapByPos.end(); i++)
2194 {
2195 size_t id = (*i).second;
2196 vec.push_back(id);
2197 }
2198 return vec;
2199 }
2200
2201 void
RemoveRegion(gcRegion & region)2202 GCStructures::RemoveRegion(gcRegion & region)
2203 {
2204
2205 // remove all blocks for this region
2206 gcIdVec locusIds = GetLocusIdsForRegionByCreation(region.GetId());
2207 for(gcIdVec::const_iterator iter=locusIds.begin(); iter != locusIds.end(); iter++)
2208 {
2209 size_t locusId = (*iter);
2210 RemoveLocus(GetLocus(locusId));
2211 }
2212
2213 // remove from display and master map
2214 // no need to delete, it lives there
2215 gcDisplayOrder::iterator displayIter =
2216 find(m_regionDisplay.begin(),m_regionDisplay.end(),region.GetId());
2217 m_regionDisplay.erase(displayIter);
2218 gcRegionMap::iterator masterIter = m_regions.find(region.GetId());
2219 m_regions.erase(masterIter);
2220
2221 }
2222
2223 void
RemoveLocus(gcLocus & locus)2224 GCStructures::RemoveLocus(gcLocus & locus)
2225 {
2226 // we do everything by id in this class
2227 size_t locId = locus.GetId();
2228
2229 // remove all blocks for this locus
2230 RemoveBlocksForLocus(locId);
2231
2232 // remove this locus from its region
2233 if(locus.HasRegion())
2234 {
2235 size_t oldRegionId = locus.GetRegionId();
2236 gcRegion & oldRegion = GetRegion(oldRegionId);
2237 oldRegion.RemoveLocusId(locId);
2238 }
2239
2240 // remove from master map -- no need to delete, it
2241 // lives there
2242 gcLocusMap::iterator masterIter = m_loci.find(locId);
2243 m_loci.erase(masterIter);
2244
2245 }
2246
2247 void
RemoveMigration(size_t fromId,size_t toId)2248 GCStructures::RemoveMigration(size_t fromId, size_t toId)
2249 {
2250 gcMigrationPair p(fromId, toId);
2251 gcMigrationMap::iterator masterIter = m_migrations.find(p);
2252 m_migrations.erase(masterIter);
2253 }
2254
2255 void
RemovePop(gcPopulation & pop)2256 GCStructures::RemovePop(gcPopulation & pop)
2257 {
2258 // we do everything by id in this class
2259 size_t popId = pop.GetId();
2260
2261 // remove all blocks for this pop
2262 RemoveBlocksForPop(popId);
2263
2264 // remove from display and master map
2265 gcDisplayOrder::iterator displayIter =
2266 find(m_popDisplay.begin(),m_popDisplay.end(),popId);
2267 m_popDisplay.erase(displayIter);
2268 gcPopMap::iterator masterIter = m_pops.find(popId);
2269 m_pops.erase(masterIter);
2270
2271 }
2272
2273 void
RemovePanel(size_t regionId,size_t popId)2274 GCStructures::RemovePanel(size_t regionId, size_t popId)
2275 {
2276 gcRegionPopIdPair p(regionId, popId);
2277 // remove from master map
2278 gcPanelMap::iterator masterIter = m_panels.find(p);
2279 m_panels.erase(masterIter);
2280 }
2281
2282 void
RemoveParent(size_t parentId)2283 GCStructures::RemoveParent(size_t parentId)
2284 {
2285 const gcParent & parRef = GetParent(parentId);
2286 m_names.FreeName(parRef.GetName());
2287 gcParentMap::iterator masterIter = m_parents.find(parentId);
2288 m_parents.erase(masterIter);
2289 }
2290
2291
2292 void
RemoveRegions(objVector regions)2293 GCStructures::RemoveRegions(objVector regions)
2294 {
2295 for(objVector::iterator i=regions.begin(); i != regions.end(); i++)
2296 {
2297 GCQuantum * object = *i;
2298 gcRegion * regionP = dynamic_cast<gcRegion *>(object);
2299 assert(regionP != NULL);
2300 RemoveRegion(*regionP);
2301 }
2302 }
2303
2304 void
RemoveLoci(objVector loci)2305 GCStructures::RemoveLoci(objVector loci)
2306 {
2307 for(objVector::iterator iter = loci.begin(); iter != loci.end(); iter++)
2308 {
2309 GCQuantum * object = *iter;
2310 gcLocus * locP = dynamic_cast<gcLocus *>(object);
2311 assert(locP != NULL);
2312 RemoveLocus(*locP);
2313 }
2314 }
2315
2316 void
RemovePops(objVector pops)2317 GCStructures::RemovePops(objVector pops)
2318 {
2319 for(objVector::iterator iter = pops.begin(); iter != pops.end(); iter++)
2320 {
2321 GCQuantum * object = *iter;
2322 gcPopulation * popP = dynamic_cast<gcPopulation *>(object);
2323 assert(popP != NULL);
2324 RemovePop(*popP);
2325 }
2326 }
2327
2328 void
RemoveParents()2329 GCStructures::RemoveParents()
2330 {
2331 for (gcParentMap::iterator i=m_parents.begin();i != m_parents.end();i++)
2332 {
2333 RemoveParent((*i).second.GetId());
2334 }
2335
2336 // clean up the populations
2337 for (gcPopMap::const_iterator i=m_pops.begin();i != m_pops.end();i++)
2338 {
2339 GetPop((*i).second.GetId()).ClearParent();
2340 }
2341 }
2342
2343 bool
GetFileSelection(size_t fileId) const2344 GCStructures::GetFileSelection(size_t fileId) const
2345 {
2346 gcFileMap::const_iterator iter = m_files.find(fileId);
2347 assert(iter != m_files.end());
2348 const gcFileInfo & fileInfo = (*iter).second;
2349 return fileInfo.GetSelected();
2350 }
2351
2352 void
SetFileSelection(size_t fileId,bool value)2353 GCStructures::SetFileSelection(size_t fileId, bool value)
2354 {
2355 gcFileMap::iterator iter = m_files.find(fileId);
2356 assert(iter != m_files.end());
2357 gcFileInfo & fileInfo = (*iter).second;
2358 return fileInfo.SetSelected(value);
2359 }
2360
2361 size_t
SelectedFileCount() const2362 GCStructures::SelectedFileCount() const
2363 {
2364 size_t count = 0;
2365 for(gcFileMap::const_iterator iter=m_files.begin(); iter != m_files.end(); iter++)
2366 {
2367 const gcFileInfo & fileInfo = (*iter).second;
2368 if(fileInfo.GetSelected())
2369 {
2370 count++;
2371 }
2372 }
2373 return count;
2374 }
2375
2376 void
AllFileSelectionsTo(bool value)2377 GCStructures::AllFileSelectionsTo(bool value)
2378 {
2379 for(gcFileMap::iterator iter=m_files.begin(); iter != m_files.end(); iter++)
2380 {
2381 gcFileInfo & fileInfo = (*iter).second;
2382 fileInfo.SetSelected(value);
2383 }
2384 }
2385
2386 const GCLocusMatcher &
GetLocusMatcher(const GCFile & fileRef) const2387 GCStructures::GetLocusMatcher(const GCFile & fileRef) const
2388 {
2389 gcFileMap::const_iterator iter = m_files.find(fileRef.GetId());
2390 assert(iter != m_files.end());
2391 const gcFileInfo & fileInfo = (*iter).second;
2392 return fileInfo.GetLocMatcher();
2393 }
2394
2395 const GCPopMatcher &
GetPopMatcher(const GCFile & fileRef) const2396 GCStructures::GetPopMatcher(const GCFile & fileRef) const
2397 {
2398 gcFileMap::const_iterator iter = m_files.find(fileRef.GetId());
2399 assert(iter != m_files.end());
2400 const gcFileInfo & fileInfo = (*iter).second;
2401 return fileInfo.GetPopMatcher();
2402 }
2403
2404 void
SetLocusMatcher(const GCFile & fileRef,const GCLocusMatcher & l)2405 GCStructures::SetLocusMatcher(const GCFile & fileRef, const GCLocusMatcher & l)
2406 {
2407 gcFileMap::iterator iter = m_files.find(fileRef.GetId());
2408 assert(iter != m_files.end());
2409 gcFileInfo & fileInfo = (*iter).second;
2410 fileInfo.SetLocMatcher(l);
2411 }
2412
2413 void
SetPopMatcher(const GCFile & fileRef,const GCPopMatcher & p)2414 GCStructures::SetPopMatcher(const GCFile & fileRef, const GCPopMatcher & p)
2415 {
2416 gcFileMap::iterator iter = m_files.find(fileRef.GetId());
2417 assert(iter != m_files.end());
2418 gcFileInfo & fileInfo = (*iter).second;
2419 fileInfo.SetPopMatcher(p);
2420 }
2421
2422 const gcPhenoMap &
GetPhenotypeMap() const2423 GCStructures::GetPhenotypeMap() const
2424 {
2425 return m_phenotypes;
2426 }
2427
2428 // EWFIX.P3 -- this should share code with ::GetParse(const GCFile&) const
2429 const GCParse &
GetParse(size_t fileId) const2430 GCStructures::GetParse(size_t fileId) const
2431 {
2432 gcFileMap::const_iterator iter = m_files.find(fileId);
2433 assert(iter != m_files.end());
2434 const gcFileInfo & fileInfo = (*iter).second;
2435 size_t parseId = fileInfo.GetParse();
2436 for(size_t i=0; i < m_dataStoreP->GetDataFile(fileId).GetParseCount(); i++)
2437 {
2438 const GCParse & parse = m_dataStoreP->GetDataFile(fileId).GetParse(i);
2439 if(parse.GetId() == parseId)
2440 {
2441 return parse;
2442 }
2443 }
2444 wxString msg = wxString::Format(gcerr::missingParse,m_dataStoreP->GetDataFile(fileId).GetName().c_str());
2445 gc_implementation_error e(msg.c_str());
2446 // implementation error -- you should have checked HasParse first!
2447 throw e;
2448 }
2449
2450 const GCParse &
GetParse(const GCFile & file) const2451 GCStructures::GetParse(const GCFile & file) const
2452 {
2453 gcFileMap::const_iterator iter = m_files.find(file.GetId());
2454 assert(iter != m_files.end());
2455 const gcFileInfo & fileInfo = (*iter).second;
2456 size_t parseId = fileInfo.GetParse();
2457 for(size_t i=0; i < file.GetParseCount(); i++)
2458 {
2459 const GCParse & parse = file.GetParse(i);
2460 if(parse.GetId() == parseId)
2461 {
2462 return parse;
2463 }
2464 }
2465 wxString msg = wxString::Format(gcerr::missingParse,file.GetName().c_str());
2466 gc_implementation_error e(msg.c_str());
2467 // implementation error -- you should have checked HasParse first!
2468 throw e;
2469 }
2470
2471 bool
HasParse(const GCFile & file) const2472 GCStructures::HasParse(const GCFile & file) const
2473 {
2474 gcFileMap::const_iterator iter = m_files.find(file.GetId());
2475 assert(iter != m_files.end());
2476 const gcFileInfo & fileInfo = (*iter).second;
2477 return fileInfo.HasParse();
2478 }
2479
2480 bool
HasUnparsedFiles() const2481 GCStructures::HasUnparsedFiles() const
2482 {
2483 for(gcFileMap::const_iterator iter=m_files.begin(); iter != m_files.end(); iter++)
2484 {
2485 const gcFileInfo & fileInfo = (*iter).second;
2486 if(!(fileInfo.HasParse())) return true;
2487 }
2488 return false;
2489 }
2490
2491 void
SetParse(const GCParse & parse)2492 GCStructures::SetParse(const GCParse & parse)
2493 {
2494 const GCFile & fileRef = parse.GetFileRef();
2495 gcFileMap::iterator iter = m_files.find(fileRef.GetId());
2496 assert(iter != m_files.end());
2497 gcFileInfo & fileInfo = (*iter).second;
2498 fileInfo.SetParse(parse.GetId());
2499
2500 // EWFIX.BUG.588
2501 if(parse.GetHasSpacesInNames())
2502 {
2503 m_dataStoreP->GCWarning(wxString::Format(gcerr_infile::shortSampleName,fileRef.GetShortName().c_str()));
2504 }
2505 }
2506
2507 void
UnsetParse(const GCParse & parse)2508 GCStructures::UnsetParse(const GCParse & parse)
2509 {
2510 const GCFile & fileRef = parse.GetFileRef();
2511 gcFileMap::iterator iter = m_files.find(fileRef.GetId());
2512 assert(iter != m_files.end());
2513 gcFileInfo & fileInfo = (*iter).second;
2514 if(fileInfo.HasParse())
2515 {
2516 if(fileInfo.GetParse() == parse.GetId())
2517 {
2518 fileInfo.UnsetParse();
2519 }
2520 }
2521
2522 for(size_t popIndex = 0; popIndex < parse.GetPopCount(); popIndex++)
2523 {
2524 for(size_t locIndex = 0; locIndex < parse.GetLociCount(); locIndex++)
2525 {
2526 const GCParseBlock & block = parse.GetBlock(popIndex,locIndex);
2527 RemoveBlockAssignment(block.GetId());
2528 }
2529 }
2530 }
2531
2532 void
AddFile(const GCFile & file)2533 GCStructures::AddFile(const GCFile & file)
2534 {
2535 gcFileInfo g;
2536 g.UnsetParse();
2537 g.SetSelected(false);
2538 m_files[file.GetId()] = g;
2539 }
2540
2541 void
RemoveFile(size_t fileId)2542 GCStructures::RemoveFile(size_t fileId)
2543 {
2544 gcFileMap::iterator iter = m_files.find(fileId);
2545 assert(iter != m_files.end());
2546 m_files.erase(iter);
2547 }
2548
2549 size_t
GetHapFileAdjacent(size_t fileId) const2550 GCStructures::GetHapFileAdjacent(size_t fileId) const
2551 {
2552 assert(HasHapFileAdjacent(fileId));
2553
2554 gcFileMap::const_iterator iter = m_files.find(fileId);
2555 assert(iter != m_files.end());
2556 const gcFileInfo & fileInfo = (*iter).second;
2557
2558 return fileInfo.GetAdjacentHapAssignment();
2559 }
2560
2561 void
SetHapFileAdjacent(size_t fileId,size_t numHaps)2562 GCStructures::SetHapFileAdjacent(size_t fileId, size_t numHaps)
2563 {
2564 gcFileMap::iterator iter = m_files.find(fileId);
2565 assert(iter != m_files.end());
2566 gcFileInfo & fileInfo = (*iter).second;
2567
2568 // EWFIX.P4
2569 // this part is here to throw an error if this would cause
2570 // a problem
2571 const GCParse & parse = GetParse(fileId);
2572 gcPhaseInfo * info = parse.GetPhaseRecordsForAdjacency(numHaps);
2573 delete info;
2574
2575 fileInfo.SetAdjacentHapAssignment(numHaps);
2576 }
2577
2578 void
UnsetHapFileAdjacent(size_t fileId)2579 GCStructures::UnsetHapFileAdjacent(size_t fileId)
2580 {
2581 gcFileMap::iterator iter = m_files.find(fileId);
2582 assert(iter != m_files.end());
2583 gcFileInfo & fileInfo = (*iter).second;
2584
2585 fileInfo.UnsetAdjacentHapAssignment();
2586 }
2587
2588 bool
HasHapFileAdjacent(size_t fileId) const2589 GCStructures::HasHapFileAdjacent(size_t fileId) const
2590 {
2591 gcFileMap::const_iterator iter = m_files.find(fileId);
2592 assert(iter != m_files.end());
2593 const gcFileInfo & fileInfo = (*iter).second;
2594 return fileInfo.HasAdjacentHapAssignment();
2595 }
2596
2597 void
UnsetGenoFile(size_t fileId)2598 GCStructures::UnsetGenoFile(size_t fileId)
2599 {
2600 gcFileMap::iterator iter = m_files.find(fileId);
2601 assert(iter != m_files.end());
2602 gcFileInfo & fileInfo = (*iter).second;
2603
2604 fileInfo.UnsetGenoFile();
2605 }
2606
2607 bool
HasGenoFile(size_t fileId) const2608 GCStructures::HasGenoFile(size_t fileId) const
2609 {
2610 gcFileMap::const_iterator iter = m_files.find(fileId);
2611 assert(iter != m_files.end());
2612 const gcFileInfo & fileInfo = (*iter).second;
2613 return fileInfo.HasGenoFile();
2614 }
2615
2616 void
VerifyLocusSeparations(const gcRegion & regionRef) const2617 GCStructures::VerifyLocusSeparations(const gcRegion& regionRef) const
2618 {
2619 constObjVector loci =
2620 GetConstDisplayableLociInMapOrderFor(regionRef.GetId());
2621
2622 bool havePrevious = false;
2623 long lastExtent = 0;
2624 size_t lastLocusId = 0;
2625
2626 for(constObjVector::const_iterator i=loci.begin(); i != loci.end(); i++)
2627 {
2628 const gcLocus * locusP = dynamic_cast<const gcLocus *>(*i);
2629 assert(locusP != NULL);
2630
2631 if(locusP->GetLinked())
2632 // no need to check unlinked loci, they are their own region
2633 {
2634 long start = 0;
2635 long stop = 0;
2636 if(locusP->HasMapPosition())
2637 // EWFIX.P3 -- should we insist ?
2638 {
2639 start += locusP->GetMapPosition();
2640 stop += locusP->GetMapPosition();
2641 }
2642 if(locusP->HasOffset())
2643 // EWFIX.P3 -- should we insist ?
2644 {
2645 start += locusP->GetOffset();
2646 stop += locusP->GetOffset();
2647 }
2648 stop += locusP->GetLength();
2649 stop -= 1; // EWFIX explain why
2650
2651 if(havePrevious)
2652 {
2653 if(start <= lastExtent)
2654 {
2655 long end = locusP->GetMapPosition()+locusP->GetLength()-1;
2656 long lastStart = GetLocus(lastLocusId).GetMapPosition();
2657 throw gc_locus_overlap(locusP->GetName(),start,end,
2658 GetLocus(lastLocusId).GetName(),lastStart,lastExtent);
2659 }
2660 }
2661
2662 havePrevious = true;
2663 lastExtent = stop;
2664 lastLocusId = locusP->GetId();
2665 }
2666 }
2667 }
2668
2669 bool
AnyZeroes() const2670 GCStructures::AnyZeroes() const
2671 {
2672 constObjVector loci = GetConstDisplayableLoci();
2673
2674 for(constObjVector::const_iterator i = loci.begin(); i != loci.end(); i++)
2675 {
2676 const GCQuantum * q = *i;
2677 const gcLocus * locusP = dynamic_cast<const gcLocus*>(q);
2678 assert(locusP != NULL);
2679
2680 if(locusP->HasMapPosition())
2681 {
2682 if (locusP->GetMapPosition() == 0) return true;
2683 }
2684
2685 if(locusP->HasOffset())
2686 {
2687 if (locusP->GetOffset() == 0) return true;
2688 }
2689
2690 if(locusP->HasLocationZero()) return true;
2691
2692 if (locusP->HasUnphasedMarkers())
2693 {
2694 const gcUnphasedMarkers * phaseInfo = locusP->GetUnphasedMarkers();
2695 if(phaseInfo != NULL)
2696 {
2697 if(phaseInfo->HasZero()) return true;
2698 }
2699 }
2700
2701 }
2702
2703 assert(m_dataStoreP != NULL);
2704 return m_dataStoreP->PhaseInfoHasAnyZeroes();
2705 }
2706
2707 int
GetPopCount()2708 GCStructures::GetPopCount()
2709 {
2710 return m_pops.size();
2711 }
2712
2713 int
GetPopCount() const2714 GCStructures::GetPopCount() const
2715 {
2716 return m_pops.size();
2717 }
2718
2719 int
GetParentCount()2720 GCStructures::GetParentCount()
2721 {
2722 return m_parents.size();
2723 }
2724
2725 int
GetParentCount() const2726 GCStructures::GetParentCount() const
2727 {
2728 return m_parents.size();
2729 }
2730
2731 bool
IsPop(size_t id)2732 GCStructures::IsPop(size_t id)
2733 {
2734 gcPopMap::iterator iter = m_pops.find(id);
2735 if (iter == m_pops.end())
2736 {
2737 return false;
2738 }
2739 return true;
2740 }
2741
2742 bool
IsPop(size_t id) const2743 GCStructures::IsPop(size_t id) const
2744 {
2745 gcPopMap::const_iterator iter = m_pops.find(id);
2746 if (iter == m_pops.end())
2747 {
2748 return false;
2749 }
2750 return true;
2751 }
2752 bool
IsPop(wxString name)2753 GCStructures::IsPop(wxString name)
2754 {
2755 for(gcPopMap::iterator i=m_pops.begin();
2756 i != m_pops.end();
2757 i++)
2758 {
2759 if((*i).second.GetName().CmpNoCase(name) == 0)
2760 {
2761 return true;
2762 }
2763 }
2764 return false;
2765 }
2766
2767 bool
IsPop(wxString name) const2768 GCStructures::IsPop(wxString name) const
2769 {
2770 for(gcPopMap::const_iterator i=m_pops.begin();
2771 i != m_pops.end();
2772 i++)
2773 {
2774 if((*i).second.GetName().CmpNoCase(name) == 0)
2775 {
2776 return true;
2777 }
2778 }
2779 return false;
2780 }
2781
2782 bool
IsParent(size_t id)2783 GCStructures::IsParent(size_t id)
2784 {
2785 gcParentMap::iterator iter = m_parents.find(id);
2786 if (iter == m_parents.end())
2787 {
2788 return false;
2789 }
2790 return true;
2791 }
2792
2793 bool
IsParent(size_t id) const2794 GCStructures::IsParent(size_t id) const
2795 {
2796 gcParentMap::const_iterator iter = m_parents.find(id);
2797 if (iter == m_parents.end())
2798 {
2799 return false;
2800 }
2801 return true;
2802 }
2803
2804 int
GetUnusedPopCount()2805 GCStructures::GetUnusedPopCount()
2806 {
2807 int count = 0;
2808 gcDisplayOrder ids = GetDisplayablePopIds();
2809 for(gcDisplayOrder::iterator iter=ids.begin(); iter != ids.end(); iter++)
2810 {
2811 size_t id = *iter;
2812 wxLogVerbose("GCStructures::GetUnusedPopCount pop id:%i name:%s", id, GetPop(id).GetName().c_str()); // JMDBG
2813 if (!GetPop(id).HasParent())
2814 {
2815 count++;
2816 wxLogVerbose(" has no parent");
2817 }
2818 else
2819 {
2820 wxLogVerbose(" has parent: %i", (int)GetPop(id).GetParentId());
2821 }
2822 }
2823 wxLogVerbose("GCStructures::GetUnusedPopCount count: %i", count);
2824 return count;
2825 }
2826
2827 int
GetUnusedPopCount() const2828 GCStructures::GetUnusedPopCount() const
2829 {
2830 int count = 0;
2831 gcDisplayOrder ids = GetDisplayablePopIds();
2832 for(gcDisplayOrder::iterator iter=ids.begin(); iter != ids.end(); iter++)
2833 {
2834 size_t id = *iter;
2835 if (!GetPop(id).HasParent())
2836 {
2837 count++;
2838 wxLogVerbose("GCStructures::GetUnusedPopCount const pop id:%i name:%s has no parent", id, GetPop(id).GetName().c_str()); // JMDBG
2839 }
2840 else
2841 {
2842 wxLogVerbose("GCStructures::GetUnusedPopCount const pop id:%i name:%s has parent: %i", id, GetPop(id).GetName().c_str(), (int)GetPop(id).GetParentId());// JMDBG
2843 }
2844 }
2845 wxLogVerbose("GCStructures::GetUnusedPopCount count: %i", count);
2846 return count;
2847 }
2848
2849 int
GetUnusedParentCount()2850 GCStructures::GetUnusedParentCount()
2851 {
2852 int count = 0;
2853 gcDisplayOrder ids = GetParentIds();
2854 for(gcDisplayOrder::iterator iter=ids.begin(); iter != ids.end(); iter++)
2855 {
2856 size_t id = *iter;
2857 if (!GetParent(id).HasParent())
2858 {
2859 count++;
2860 wxLogVerbose("GCStructures::GetUnusedParentCount par id: %i name:%s has no parent", id, GetParent(id).GetName().c_str()); // JMDBG
2861 }
2862 else
2863 {
2864 wxLogVerbose("GCStructures::GetUnusedParentCount par id: %i name:%s has parent: %s", id, GetParent(id).GetName().c_str(), (int)GetParent(id).GetParentId()); // JMDBG
2865 }
2866 }
2867 wxLogVerbose("GCStructures::GetUnusedParentCount count: %i", count);
2868 return count;
2869 }
2870
2871 int
GetUnusedParentCount() const2872 GCStructures::GetUnusedParentCount() const
2873 {
2874 int count = 0;
2875 gcDisplayOrder ids = GetParentIds();
2876 for(gcDisplayOrder::iterator iter=ids.begin(); iter != ids.end(); iter++)
2877 {
2878 size_t id = *iter;
2879 wxLogVerbose("GCStructures::GetUnusedParentCount const par id: %i name:%s", id, GetParent(id).GetName().c_str()); // JMDBG
2880 if (!GetParent(id).HasParent())
2881 {
2882 count++;
2883 wxLogVerbose(" has no parent");
2884 }
2885 else
2886 {
2887 wxLogVerbose(" has parent: &s", (int)GetParent(id).GetParentId());
2888 }
2889 }
2890 wxLogVerbose("GCStructures::GetUnusedParentCount count: %i", count);
2891 return count;
2892 }
2893
2894 void
ClearPopDisplayOrder()2895 GCStructures::ClearPopDisplayOrder()
2896 {
2897 gcDisplayOrder ids = GetDisplayablePopIds() ;
2898 for(gcDisplayOrder::iterator iter=ids.begin(); iter != ids.end(); iter++)
2899 {
2900 GetPop(*iter).SetDispOrder(0);
2901 }
2902 }
2903
2904 size_t
FindTopParent()2905 GCStructures::FindTopParent()
2906 {
2907 size_t retid = gcdefault::badIndex;
2908 gcDisplayOrder ids = GetParentIds();
2909 for(gcDisplayOrder::iterator iter=ids.begin(); iter != ids.end(); iter++)
2910 {
2911 if (!GetParent(*iter).HasParent())
2912 {
2913 retid = *iter;
2914 }
2915 }
2916 assert(retid != gcdefault::badIndex);
2917 return retid;
2918 }
2919
2920 bool
HasParents() const2921 GCStructures::HasParents() const
2922 {
2923 if (m_parents.size() > 0)
2924 {
2925 return true;
2926 }
2927 else
2928 {
2929 return false;
2930 }
2931 }
2932
2933 void
MakeMigrationMatrix()2934 GCStructures::MakeMigrationMatrix()
2935 {
2936 // convenience to make the code easier to read
2937 gcDisplayOrder popids = GetDisplayablePopIds();
2938 gcDisplayOrder parids = GetParentIds();
2939
2940 int matrixDim = popids.size() + parids.size();
2941
2942 assert(matrixDim > 0);
2943
2944 size_t idex;
2945 size_t jdex;
2946
2947 // build whole matrix
2948 for (int i=0; i<matrixDim; i++)
2949 {
2950 for (int j=0; j<matrixDim; j++)
2951 {
2952 if (j!=i)
2953 {
2954 if (i < (int)popids.size())
2955 {
2956 int count = 0;
2957 for(gcDisplayOrder::const_iterator iter = popids.begin(); iter != popids.end(); iter++)
2958 {
2959 if (count == i)
2960 {
2961 idex = *iter;
2962 }
2963 count ++;
2964 }
2965 }
2966 else
2967 {
2968 int ioff = i - popids.size();
2969 int count = 0;
2970 for(gcDisplayOrder::const_iterator iter = parids.begin(); iter != parids.end(); iter++)
2971 {
2972 if (count == ioff)
2973 {
2974 idex = *iter;
2975 }
2976 count ++;
2977 }
2978 }
2979
2980 if (j < (int) popids.size())
2981 {
2982 int count = 0;
2983 for(gcDisplayOrder::const_iterator jter = popids.begin(); jter != popids.end(); jter++)
2984 {
2985 if (count == j)
2986 {
2987 jdex = *jter;
2988 }
2989 count ++;
2990 }
2991 }
2992 else
2993 {
2994 int joff = j - popids.size();
2995 int count = 0;
2996 for(gcDisplayOrder::const_iterator jter = parids.begin(); jter != parids.end(); jter++)
2997 {
2998 if (count == joff)
2999 {
3000 jdex = *jter;
3001 }
3002 count ++;
3003 }
3004 }
3005 if (!HasMigration(idex, jdex))
3006 {
3007 gcMigration Mig1 = MakeMigration(true, idex, jdex );
3008 }
3009 if (!HasMigration(jdex, idex))
3010 {
3011 gcMigration Mig2 = MakeMigration(true, jdex, idex );
3012 }
3013 }
3014 }
3015 }
3016
3017 if (GetUnusedParentCount() < 2)
3018 {
3019 // divergence linkages finished so need to remove some of the matrix members
3020
3021 // creation order list - used to decide which children have been used
3022 // and are thus no longer available for migration
3023 size_t* orderids = new size_t[matrixDim];
3024 for (int i=0; i<matrixDim; i++)
3025 {
3026 orderids[i] = gcdefault::badIndex;
3027 }
3028
3029 // populations can always migrate between each other before parents are declared
3030 int mdx = 0;
3031 for(gcDisplayOrder::const_iterator iter = popids.begin(); iter != popids.end(); iter++)
3032 {
3033 orderids[mdx] = *iter;
3034 mdx++;
3035 }
3036
3037 // now remove those that disappear as parents are declared
3038 for(gcDisplayOrder::const_iterator iter = parids.begin(); iter != parids.end(); iter++)
3039 {
3040 orderids[mdx] = *iter;
3041 const gcParent & parent = GetParent(*iter);
3042
3043 if (!parent.HasParent())
3044 {
3045 // no migrations possible at top
3046 for(gcDisplayOrder::const_iterator jter = popids.begin(); jter != popids.end(); jter++)
3047 {
3048 RemoveMigration(*iter, *jter );
3049 RemoveMigration(*jter, *iter );
3050 }
3051
3052 for(gcDisplayOrder::const_iterator jter = parids.begin(); jter != parids.end(); jter++)
3053 {
3054 if (*jter != *iter)
3055 {
3056 RemoveMigration(*iter, *jter );
3057 RemoveMigration(*jter, *iter );
3058 }
3059 }
3060 }
3061 else
3062 {
3063 size_t child1Id = parent.GetChild1Id();
3064 size_t child2Id = parent.GetChild2Id();
3065 for (int i=0; i<matrixDim; i++)
3066 {
3067 if (child1Id == orderids[i])
3068 {
3069 RemoveMigration(*iter, child1Id );
3070 RemoveMigration(child1Id, *iter );
3071 }
3072 if (child2Id == orderids[i])
3073 {
3074 RemoveMigration(*iter, child2Id );
3075 RemoveMigration(child2Id, *iter );
3076 }
3077 }
3078 }
3079 mdx ++;
3080 }
3081 }
3082 }
3083
3084 //____________________________________________________________________________________
3085