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