FivePointsNumericalDerivative.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1

#include "FivePointsNumericalDerivative.h"

using namespace bpp;
using namespace std;

void FivePointsNumericalDerivative::updateDerivatives(const ParameterList& parameters)
{
  if (computeD1_ && variables_.size() > 0)
  {
    if (function1_)
      function1_->enableFirstOrderDerivatives(false);
    if (function2_)
      function2_->enableSecondOrderDerivatives(false);
    function_->setParameters(parameters);
    f3_ = function_->getValue();
    string lastVar;
    bool functionChanged = false;
    ParameterList p;
    bool start = true;
    for (size_t i = 0; i < variables_.size(); ++i)
    {
      string var = variables_[i];
      if (!parameters.hasParameter(var))
        continue;
      if (!start)
      {
        vector<string> vars(2);
        vars[0] = var;
        vars[1] = lastVar;
        p = parameters.createSubList(vars);
      }
      else
      {
        p = parameters.createSubList(var);
        start = false;
      }
      lastVar = var;
      functionChanged = true;
      double value = function_->getParameterValue(var);
      double h = (1. + std::abs(value)) * h_;
      // Compute four other points:
      try
      {
        p[0].setValue(value - 2 * h);
        function_->setParameters(p);
        f1_ = function_->getValue();
        try
        {
          p[0].setValue(value + 2 * h);
          function_->setParameters(p);
          f5_ = function_->getValue();
          // No limit raised, use central approximation:
          p[0].setValue(value - h);
          function_->setParameters(p);
          f2_ = function_->getValue();
          p[0].setValue(value + h);
          function_->setParameters(p);
          f4_ = function_->getValue();
          der1_[i] = (f1_ - 8. * f2_ + 8. * f4_ - f5_) / (12. * h);
          der2_[i] = (-f1_ + 16. * f2_ - 30. * f3_ + 16. * f4_ - f5_) / (12. * h * h);
        }
        catch (ConstraintException& ce)
        {
          // Right limit raised, use backward approximation:
          p[0].setValue(value - h);
          function_->setParameters(p);
          f2_ = function_->getValue();
          p[0].setValue(value - 2 * h);
          function_->setParameters(p);
          f1_ = function_->getValue();
          der1_[i] = (f3_ - f2_) / h;
          der2_[i] = (f3_ - 2. * f2_ + f1_) / (h * h);
        }
      }
      catch (ConstraintException& ce)
      {
        // Left limit raised, use forward approximation:
        p[0].setValue(value + h);
        function_->setParameters(p);
        f4_ = function_->getValue();
        p[0].setValue(value + 2 * h);
        function_->setParameters(p);
        f5_ = function_->getValue();
        der1_[i] = (f4_ - f3_) / h;
        der2_[i] = (f5_ - 2. * f4_ + f3_) / (h * h);
      }
    }
    // Reset last parameter and compute analytical derivatives if any.
    if (function1_)
      function1_->enableFirstOrderDerivatives(computeD1_);
    if (function2_)
      function2_->enableSecondOrderDerivatives(computeD2_);
    if (functionChanged)
      function_->setParameters(parameters.createSubList(lastVar));
  }
  else
  {
    // Reset initial value and compute analytical derivatives if any.
    if (function1_)
      function1_->enableFirstOrderDerivatives(computeD1_);
    if (function2_)
      function2_->enableSecondOrderDerivatives(computeD2_);
    function_->setParameters(parameters);
    // Just in  case derivatives are not computed:
    f3_ = function_->getValue();
  }
}