1 // $Id: maximizer.h,v 1.36 2012/04/17 19:03:20 ewalkup Exp $ 2 3 /* 4 Copyright 2002 Peter Beerli, Mary Kuhner, Jon Yamato and Joseph Felsenstein 5 6 This software is distributed free of charge for non-commercial use 7 and is copyrighted. Of course, we do not guarantee that the software 8 works, and are not responsible for any damage you may cause or have. 9 */ 10 11 // Maximizer Class ---------------------------------------------- 12 // 13 // Maximizer for likelihood function (see PostLike::Calculate()) 14 // using the Broyden-Fletcher-Shanno-Goldfarb method 15 // other methods (e.g. Newton-Raphson) can be easily incorporated 16 // by adding another Maximizer::Calculate() 17 // 18 // 19 // calculates the maximum of the likelihood function and returns 20 // the parameters as a vector at the ML. 21 // 22 // Call sequence and existence through run 23 // foreach(locus) 24 // { 25 // likelihood = SinglePostLike(forces) 26 // Maximizer(likelihood) 27 // foreach(replicate) [independent replicates, if any] 28 // { 29 // foreach(single_chains) 30 // Calculate(data) uses data and MLE 31 // foreach(long_chain) 32 // Calculate(data) 33 // } 34 // ~Maximizer() 35 // likelihood = ReplicatePostLike(forces) 36 // Maximizer(likelihood) 37 // ~Maximizer() 38 // } 39 // if(locus>1) 40 // { 41 // likelihood = RegionPostLike(forces) 42 // Maximizer(likelihood) 43 // ~Maximizer() 44 // likelihood = GammaRegionPostLike(forces) 45 // Maximizer(likelihood) 46 // ~Maximizer() 47 // } 48 // STILL TODO 49 // Line() is not yet implemented 50 // 51 52 #ifndef MAXIMIZER_H 53 #define MAXIMIZER_H 54 55 #include <algorithm> 56 #include <cmath> 57 #include <functional> 58 #include <iostream> 59 #include <numeric> 60 #include <vector> 61 62 #include "definitions.h" 63 #include "vectorx.h" 64 #include "likelihood.h" 65 #include "runreport.h" 66 #include "paramstat.h" 67 68 using 69 std::vector; 70 71 // declarations -------------------------------------------- 72 73 class 74 Maximizer 75 { 76 public: 77 Maximizer (long int thisNparam); 78 ~Maximizer (); 79 80 void Initialize(long int thisNparam); 81 void SetLikelihood (PostLike * thispostlike); 82 bool Calculate (DoubleVec1d & thisparam, double& lnL, string& message); GetPostLike()83 PostLike *GetPostLike() 84 { 85 return m_pPostLike; // Analyzer and others need access to this. 86 }; 87 double GetLastNorm(); // Debugging function. 88 void SetMLEs(DoubleVec1d newMLEs); // Need for sumfile reading. 89 90 // Interface for performing constrained maximization-- 91 // e.g., maximization under the constraint that forward and backward 92 // migration rates are equal (M12 == M21). 93 bool SetConstraints(const vector<vector<unsigned long int> > & constraints); ClearConstraints(void)94 void ClearConstraints(void) { m_constraints.clear(); }; 95 void AppendConstraintOnAlpha(ParamStatus alphaStatus); // Append to the gradient guide. 96 97 // Gradient guide interface, moved here from the PostLike class, 98 // because it's better for the maximizer to control it. 99 // (The gradient guide still is, and will be, used by PostLike.) 100 void GradientGuideSetup (const ParamVector &thisToDolist, PostLike* pPL); 101 void ProfileGuideRestore(); 102 void ProfileGuideFix(long int guide); 103 void ProfileGuideFixAll(); 104 void ProfileGuideRestore(long int guide); GetPostlikeTag()105 likelihoodtype GetPostlikeTag() { return m_pPostLike->GetTag(); }; 106 107 private: Maximizer()108 Maximizer () 109 { 110 }; 111 double m_lastnormg; 112 113 protected: 114 PostLike *m_pPostLike; // pointer to PostLike object 115 unsigned long int m_nparam; // number of parameters (m_param.size()) 116 long int m_nloci; // number of loci 117 long int m_nrep; // number of replicates 118 DoubleVec1d m_param; // parameters 119 DoubleVec1d m_lparam; // log parameters 120 DoubleVec1d m_oldlparam; // old log parameters 121 DoubleVec1d m_gradient; // gradient 122 DoubleVec1d m_oldgradient; // old gradient 123 DoubleVec1d m_paramdelta; // m_newparam - m_oldparam 124 DoubleVec1d m_gradientdelta; // m_gradient - m_oldgradient 125 DoubleVec1d m_direction; // direction 126 vector<DoubleVec1d> m_second; // approx second derivative 127 128 vector<unsigned long int> m_isLinearParam; // track which parameters get treated 129 // logarithmically and which linearly 130 bool m_dataHasLinearParam; 131 bool m_dataHasGamma; 132 bool m_constraintsVectorHasSlotForAlpha; 133 double m_maxAllowedValueForAlpha; 134 135 // MARY -- temporary storage, kept as class variables for speed 136 DoubleVec1d m_newparam; 137 DoubleVec1d m_newlparam; 138 DoubleVec1d m_minparamvalues; 139 DoubleVec1d m_minlparamvalues; 140 DoubleVec1d m_maxparamvalues; 141 DoubleVec1d m_maxlparamvalues; 142 DoubleVec1d m_temp; 143 DoubleVec2d m_dd; 144 145 vector<vector<unsigned long int> > m_constraints; // user can constrain params 146 // to equal one another 147 bool m_constrained; // assume this is faster than evaluating m_constraints.empty() 148 DoubleVec1d m_constraintratio; // for multiplicative constraints 149 150 double Norm (const DoubleVec1d&d); // utility function: norm of vector 151 // (i.e., its length) 152 153 bool CalculateSteepest(double&, string& message); 154 bool CalculateByParts(double&, string& message); 155 bool CalculateBroyden(unsigned long int maxtrials, double &, string & message); 156 157 // calculates the new parameter 158 double SetParam(const double &lambda, const DoubleVec1d& direction); 159 bool SetParam1d(const double& step, const vector<unsigned long int> * pWhichparam); 160 161 void CoutNewParams() const; // for debugging 162 void CoutCurParams() const; // for debugging 163 void CoutByLinOrNoMult(double mult, const DoubleVec1d&) const; // for debugging 164 165 // Used by CalculateSteepest(). 166 bool BracketTheMaximumAndFindIt(const vector<unsigned long int> 167 * pWhichparam, const double epsilon, 168 double oldlike, double & newlike, 169 string & message); 170 171 // Handles constrained parameters when applicable. 172 bool DCalculate(const DoubleVec1d & param, DoubleVec1d & gradient); 173 174 void ResetSecond (); //resets the second derivative matrix 175 void CalcSecond (); // calculate approximative second derivatives 176 void CalcDirection (); // calculate the direction of change 177 SetGradGuide(const vector<ParamStatus> & thisguide)178 void SetGradGuide(const vector<ParamStatus>& thisguide) 179 { 180 std::copy(thisguide.begin(), thisguide.end(), 181 m_pPostLike->m_working_pstatusguides.begin()); 182 }; 183 184 void ExplainParamChange(long index, double oldVal, double newVal, std::string kind); 185 void AdjustExtremeLinearParam(long i); 186 void AdjustExtremeLogParam(long i); 187 }; 188 189 // helper functions 190 191 // simply returns result = one - two 192 void CalcDelta (const DoubleVec1d & one, const DoubleVec1d & two, DoubleVec1d & result); 193 194 inline bool SlopesHaveSameSign(const double & oldgrad, const double & newgrad, const double & tolerance); 195 inline double sign(const double & x); 196 inline double tolerance(double number); 197 198 #endif // MAXIMIZER_H 199 200 //____________________________________________________________________________________ 201