1 // $Id: dlcalc.h,v 1.44 2011/10/11 16:42:16 mkkuhner 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 #ifndef DLCALCULATOR_H 12 #define DLCALCULATOR_H 13 14 #include <cmath> 15 #include <vector> 16 17 #include "vectorx.h" 18 #include "constants.h" 19 #include "rangex.h" 20 #include "dlcell.h" 21 #include "branch.h" 22 23 /********************************************************************* 24 DLCalculator manages the calculation of data likelihoods on 25 a given locus specified at constructor time. Call 26 the member function Calculate() to perform the actual calculation. 27 Calculate() makes use of DataModel member functions 28 ComputeExponentials(), ComputeSiteDLs() and ComputeTreeDL(), 29 and returns the likelihood. 30 31 The class is polymorphic on data type (not data model). 32 33 We assume that we start with a correct tree with all update flags set 34 appropriately. We end with a tree with correct likelihoods. This 35 class does NOT reset the update flags after updating the likelihood, 36 as the flags may be needed by subsequent calculations. 37 38 This file also holds the small helper class ChildInfo used to 39 simplify data likelihood calculation code. 40 41 Written by Jon Yamato, revised by Jim Sloan, revised by Jon Yamato 42 2002/01/28 added microsatellite support -- Mary Kuhner 43 moved markerweights to base class 44 2002/07/08 added k-allele model -- Mary 45 2002/07/21 refactoring to remove duplicate code -- Mary 46 2004/07/13 added trait-locus code -- Mary 47 2005/03/01 moved aliases to tree -- Mary 48 2005/09/22 deleted markerweights entirely. 49 50 **********************************************************************/ 51 52 // The following are used in the .cpp code: 53 // #include "datapack.h" for access to 54 // datamodel, GetDataLength(), GetNTips() 55 // #include "dlmodel.h" for access to 56 // Finalize(), ComputeExponentials(), ComputeSiteDLs(), 57 // ComputeTreeDLs() 58 // #include "tree.h" for access to 59 // timelist, SetDLValue() 60 // #include "timelist.h" for access to 61 // BeginBranch(), GetBranch(), BeginBody(), InRange() 62 63 class Locus; 64 class Tree; 65 class DataModel; 66 class NucModel; 67 class ChildInfo; 68 class TimeList; 69 70 //------------------------------------------------------------------------------------ 71 //------------------------------------------------------------------------------------ 72 73 class DLCalculator 74 { 75 private: 76 DLCalculator(); // undefined 77 DLCalculator& operator=(const DLCalculator&); // undefined 78 79 protected: 80 DataModel& m_datamodel; 81 LongVec1d m_markerpos; 82 83 DLCalculator(const DLCalculator& src); // internal use 84 85 rangepair SitePairToMarkerPair(rangepair sites); 86 ChildInfo GetChildInfo(Branch_ptr branch, long locus, long childindex, 87 long cellindex, long posn, bool moving) const; 88 ChildInfo NullChildInfo() const; 89 90 public: 91 DLCalculator(const Locus& locus); ~DLCalculator()92 virtual ~DLCalculator() {}; 93 virtual DLCalculator* Clone() const = 0; 94 virtual double Calculate(Tree& tree, const Locus& locus, bool moving) = 0; 95 virtual void SimulateData(Tree& tree, Locus& locus); 96 virtual void CopyDataFrom(Locus& destloc, Locus& origloc, Tree& tree); 97 virtual void Randomize(Locus& destloc, rangeset rset, Tree& tree); 98 virtual void MarkPanelBranches(Tree& tree, const Locus& locus); RecalculateAliases(const Tree &,const Locus &)99 virtual LongVec1d RecalculateAliases(const Tree&, const Locus&) const 100 { LongVec1d empty; return empty; }; 101 102 }; 103 104 //------------------------------------------------------------------------------------ 105 //------------------------------------------------------------------------------------ 106 107 class NucCalculator : public DLCalculator 108 { 109 private: 110 LongVec1d CalculateAliasesFromData(const std::vector<DoubleVec2d>& data) const; 111 LongVec1d SetupAliases(const Locus& locus) const; 112 113 protected: 114 typedef double** SiteDL; 115 116 NucCalculator(const NucCalculator& src); // internal use 117 virtual void CalculateSite(Cell_ptr child1, Cell_ptr child2, 118 Cell_ptr newcell, long pos, long alias); 119 void Breakalias(LongVec1d& aliases, const rangevector& subtrees); 120 121 public: 122 NucCalculator(const Locus& locus); ~NucCalculator()123 virtual ~NucCalculator() {}; 124 virtual LongVec1d RecalculateAliases(const Tree& tree, const Locus& locus) const; 125 126 }; 127 128 //------------------------------------------------------------------------------------ 129 //------------------------------------------------------------------------------------ 130 131 class DNACalculator : public NucCalculator 132 { 133 private: 134 135 protected: DNACalculator(const DNACalculator & src)136 DNACalculator(const DNACalculator& src) : NucCalculator(src) {}; 137 138 public: DNACalculator(const Locus & locus)139 DNACalculator(const Locus& locus) : NucCalculator(locus) {}; ~DNACalculator()140 virtual ~DNACalculator() {}; 141 virtual DLCalculator* Clone() const; 142 virtual double Calculate(Tree& tree, const Locus& locus, bool moving); 143 }; 144 145 //------------------------------------------------------------------------------------ 146 //------------------------------------------------------------------------------------ 147 148 class SNPCalculator : public NucCalculator 149 { 150 private: 151 NucModel* m_invarmodel; 152 153 void CalculateInvarSite(Cell_ptr child1, Cell_ptr child2, Cell_ptr newcell, long pos); 154 void MaskPanelTips(TimeList& timeList, long loci, long first, long last, long i); 155 156 protected: 157 SNPCalculator(const SNPCalculator& src); 158 159 // no-panel pathway 160 double CalcNoPanel(Tree& tree, long loc, std::pair<long,long> sitespan, LongVec1d aliases); 161 double SumNoPanel(Cell_ptr varcell, Cell_ptr invarcell, std::pair<long,long> sitespan); 162 163 // panel pathway 164 Branch_ptr CalcPanel(Tree& tree, long loc, std::pair<long,long> sitespan, LongVec1d aliases); 165 Branch_ptr CalcPanelInvariants(Tree& tree, long loc, std::pair<long,long> sitespan); 166 double SumPanel(Cell_ptr varcell, std::pair<long,long> sitespan, Cell_ptr sumcell); 167 168 public: 169 SNPCalculator(const Locus& locus); 170 virtual ~SNPCalculator(); 171 virtual DLCalculator* Clone() const; 172 virtual double Calculate(Tree& tree, const Locus& locus, bool moving); 173 virtual void SimulateData(Tree& tree, Locus& locus); 174 175 }; 176 177 //------------------------------------------------------------------------------------ 178 //------------------------------------------------------------------------------------ 179 // This subclass is suitable for any type of allelic data, such as 180 // microsats, electrophoretic alleles or the k-allele data type. 181 // It used to have subclasses, but they weren't needed. 182 183 class AlleleCalculator : public DLCalculator 184 { 185 protected: 186 AlleleCalculator(const AlleleCalculator& src); // internal use 187 188 public: 189 AlleleCalculator(const Locus& locus); ~AlleleCalculator()190 virtual ~AlleleCalculator() {}; Clone()191 virtual DLCalculator* Clone() const { return new AlleleCalculator(*this); }; 192 virtual double Calculate(Tree& tree, const Locus& locus, bool moving); 193 194 }; 195 196 //------------------------------------------------------------------------------------ 197 //------------------------------------------------------------------------------------ 198 // This tiny class saves some information for use by data likelihood 199 // calculations. It is similar to std::pair. 200 201 class ChildInfo 202 { 203 public: 204 Cell_ptr m_cell; 205 double m_length; 206 ChildInfo(Cell_ptr cell,double len)207 ChildInfo(Cell_ptr cell, double len) : m_cell(cell), m_length(len) {}; ChildInfo()208 ChildInfo() : m_length(0) {}; 209 }; 210 211 #endif // DLCALCULATOR_H 212 213 //____________________________________________________________________________________ 214