42 #include "../TreeTemplateTools.h" 46 #include <Bpp/Numeric/VectorTools.h> 47 #include <Bpp/Numeric/Function/ConjugateGradientMultiDimensions.h> 48 #include <Bpp/Text/TextTools.h> 49 #include <Bpp/App/ApplicationTools.h> 57 return NumTools::abs<double>(
59 dynamic_cast<const PseudoNewtonOptimizer*>(optimizer_)->
previousValue_);
66 return getCurrentTolerance() < tolerance_;
72 AbstractOptimizer(function),
80 setDefaultStopCondition_(
new FunctionStopCondition(
this));
81 setStopCondition(*getDefaultStopCondition());
88 n_ = getParameters().size();
89 params_ = getParameters().getParameterNames();
90 getFunction()->enableSecondOrderDerivatives(
true);
91 getFunction()->setParameters(getParameters());
98 ParameterList* bckPoint = 0;
99 if (updateParameters()) bckPoint =
new ParameterList(
getFunction()->getParameters());
102 std::vector<double> movements(
n_);
103 ParameterList newPoint = getParameters();
104 for (
size_t i = 0; i <
n_; i++)
108 if (secondOrderDerivative == 0)
112 else if (secondOrderDerivative < 0)
114 printMessage(
"!!! Second order derivative is negative for parameter " +
params_[i] +
"(" + TextTools::toString(getParameters()[i].getValue()) +
"). Moving in the other direction.");
117 movements[i] = -firstOrderDerivative / secondOrderDerivative;
119 else movements[i] = firstOrderDerivative / secondOrderDerivative;
120 if (std::isnan(movements[i]))
122 printMessage(
"!!! Non derivable point at " +
params_[i] +
". No move performed. (f=" + TextTools::toString(currentValue_) +
", d1=" + TextTools::toString(firstOrderDerivative) +
", d2=" + TextTools::toString(secondOrderDerivative) +
").");
127 newPoint[i].setValue(getParameters()[i].getValue() - movements[i]);
129 movements[i] = getParameters()[i].getValue() - newPoint[i].getValue();
134 unsigned int count = 0;
135 while ((count <
maxCorrection_) && ((newValue > currentValue_ + getStopCondition()->getTolerance()) || std::isnan(newValue))) {
137 if ((count==0) && updateParameters())
getFunction()->setParameters(*bckPoint);
139 if (!(
useCG_ && (count==3))){
140 printMessage(
"!!! Function at new point is greater than at current point: " + TextTools::toString(newValue) +
">" + TextTools::toString(currentValue_) +
". Applying Felsenstein-Churchill correction: " + TextTools::toString(count));
142 for (
size_t i = 0; i < movements.size(); i++) {
143 movements[i] = movements[i] / 2;
144 newPoint[i].setValue(getParameters()[i].getValue() - movements[i]);
149 printMessage(
"!!! Felsenstein-Churchill correction applied too many times.");
150 printMessage(
"Use conjugate gradients optimization.");
151 getFunction()->enableSecondOrderDerivatives(
false);
152 ConjugateGradientMultiDimensions opt(
getFunction());
153 opt.setConstraintPolicy(getConstraintPolicy());
154 opt.setProfiler(getProfiler());
155 opt.setMessageHandler(getMessageHandler());
157 double tol = std::max(getStopCondition()->getCurrentTolerance() / 2., getStopCondition()->getTolerance());
158 opt.getStopCondition()->setTolerance(tol);
159 opt.setMaximumNumberOfEvaluations(nbEvalMax_);
161 opt.init(getParameters());
163 newPoint = opt.getParameters();
164 newValue = opt.getFunctionValue();
166 if (newValue > currentValue_ + tol) {
167 printMessage(
"!!! Conjugate gradient method failed to improve likelihood.");
168 printMessage(
"Back to Felsenstein-Churchill method.");
174 if (newValue > currentValue_ + getStopCondition()->getTolerance()){
175 printMessage(
"PseudoNewtonOptimizer::doStep. Value could not be ameliorated!");
176 newValue = currentValue_;
185 getParameters_() = newPoint;
188 if (updateParameters())
delete bckPoint;
PseudoNewtonOptimizer(DerivableSecondOrder *function)
ParameterList previousPoint_
std::vector< std::string > params_
void doInit(const ParameterList ¶ms)
bool isToleranceReached() const
This Optimizer implements Newton's algorithm for finding a minimum of a function. This is in fact a m...
unsigned int maxCorrection_
double getCurrentTolerance() const
const DerivableSecondOrder * getFunction() const