CostFuncPoisson.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2019 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidCurveFitting/CostFunctions/CostFuncPoisson.h"
#include "MantidAPI/CompositeDomain.h"
#include "MantidAPI/FunctionDomain.h"
#include "MantidAPI/FunctionValues.h"
#include "MantidAPI/IConstraint.h"
#include "MantidCurveFitting/Jacobian.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/MultiThreaded.h"
 
#include <cmath>
#include <limits>
 
using namespace Mantid::API;
 
namespace {
// predicted < 0 is forbidden as it causes inf cost
constexpr double absoluteCutOff = 0.0;
constexpr double effectiveCutOff = 0.0001;
constexpr double epsilon = 1e-10;
 
double calculatePoissonLoss(double observedCounts, double predicted) {
  double retVal = (predicted - observedCounts);
  retVal += observedCounts * (log(observedCounts) - log(predicted));
  return retVal;
}
 
} // namespace
 
namespace Mantid::CurveFitting::CostFunctions {
 
namespace PoissonLossLM {
template MANTID_CURVEFITTING_DLL int sgn<double>(double val);
 
double calculatePoissonResidualLM(double observedCounts, double predicted) {
  observedCounts = std::max(absoluteCutOff, observedCounts);
  double retVal = (predicted - observedCounts);
  if (predicted <= absoluteCutOff) {
    return std::numeric_limits<double>::max();
  }
  // at observed = 0 the Poisson function reduces to predicted - observed
  if (observedCounts != 0) {
    if (const auto x = retVal / observedCounts; std::abs(x) < effectiveCutOff) {
      // taylor around small
      return observedCounts * (x * x) * (0.5 - (x / 3.0) + (x * x) / 4.0);
    }
 
    retVal += observedCounts * (log(observedCounts) - log(predicted));
  }
  return retVal;
}
double calculateJacobianScaleFactor(double observedCounts, double predicted) {
  if (predicted <= absoluteCutOff) {
    return std::numeric_limits<double>::max();
  }
  const double delta = predicted - observedCounts;
  const double signDelta = sgn(delta);
  // If observed is zero
  if (observedCounts == absoluteCutOff) {
    return predicted > epsilon ? signDelta / std::sqrt(2.0 * predicted) : signDelta / epsilon;
  }
  // taylor
  if (const auto x = delta / observedCounts; std::abs(x) < effectiveCutOff) {
    return (1.0 / std::sqrt(observedCounts)) * (1.0 - (2.0 / 3.0) * x + (7.0 / 12.0) * (x * x));
  }
 
  return signDelta * (1 - observedCounts / predicted) * 1 /
         std::sqrt(2 * calculatePoissonResidualLM(observedCounts, predicted));
}
 
// calculate loss for the LM minimizer as needed for gsl minimizer, taking the root of the residual.
double calculatePoissonLossLM(double observedCounts, double predicted) {
  return sgn(predicted - observedCounts) * std::sqrt(2 * calculatePoissonResidualLM(observedCounts, predicted));
}
 
} // namespace PoissonLossLM
 
DECLARE_COSTFUNCTION(CostFuncPoisson, Poisson)
//----------------------------------------------------------------------------------------------
CostFuncPoisson::CostFuncPoisson() : CostFuncFitting() {}
 
void CostFuncPoisson::addVal(API::FunctionDomain_sptr domain, API::FunctionValues_sptr values) const {
  m_function->function(*domain, *values);
  size_t ny = values->size();
 
  double retVal = 0.0;
 
  for (size_t i = 0; i < ny; i++) {
    const double predicted = values->getCalculated(i);
 
    if (predicted <= absoluteCutOff) {
      retVal = std::numeric_limits<double>::infinity();
      break;
    }
 
    const double observed = values->getFitData(i);
    if (predicted <= effectiveCutOff) {
      retVal += (effectiveCutOff - predicted) / predicted;
    } else if (observed == 0.0) {
      // at observed = 0 the Poisson function reduces to simply adding predicted
      retVal += predicted;
    } else {
      retVal += calculatePoissonLoss(observed, predicted);
    }
  }
 
  PARALLEL_ATOMIC
  m_value += 2.0 * retVal;
}
 
void CostFuncPoisson::addValDerivHessian(API::IFunction_sptr function, API::FunctionDomain_sptr domain,
                                         API::FunctionValues_sptr values, bool evalDeriv, bool evalHessian) const {
  const size_t numParams = nParams();
 
  if (evalDeriv) {
    m_der.resize(numParams);
    m_der.zero();
    calculateDerivative(*function, *domain, *values);
  }
 
  if (evalHessian) {
    m_hessian.resize(numParams, numParams);
    m_hessian.zero();
    calculateHessian(*function, *domain, *values);
  }
}
 
void CostFuncPoisson::calculateDerivative(API::IFunction &function, const FunctionDomain &domain,
                                          FunctionValues &values) const {
  const size_t numParams = function.nParams();
  const size_t numDataPoints = domain.size();
 
  Jacobian jacobian(numDataPoints, numParams);
  function.function(domain, values);
  function.functionDeriv(domain, jacobian);
 
  size_t activeParamIndex = 0;
  double costVal = 0.0;
 
  for (size_t paramIndex = 0; paramIndex < numParams; ++paramIndex) {
    if (!function.isActive(paramIndex))
      continue;
 
    double determinant = 0.0;
    for (size_t i = 0; i < numDataPoints; ++i) {
      double calc = values.getCalculated(i);
      double obs = values.getFitData(i);
 
      if (calc <= absoluteCutOff) {
        if (activeParamIndex == 0) {
          costVal += std::numeric_limits<double>::infinity();
        }
        determinant += std::numeric_limits<double>::infinity();
        continue;
      }
 
      if (calc <= effectiveCutOff) {
        if (activeParamIndex == 0) {
          costVal += (effectiveCutOff - calc) / (calc - absoluteCutOff);
        }
        double tmp = calc - absoluteCutOff;
        determinant += jacobian.get(i, paramIndex) * (absoluteCutOff - effectiveCutOff) / (tmp * tmp);
 
      } else if (obs == 0.0) {
        if (activeParamIndex == 0) {
          costVal += calc;
        }
        determinant += jacobian.get(i, paramIndex);
 
      } else {
        if (activeParamIndex == 0) {
          costVal += calculatePoissonLoss(obs, calc);
        }
        determinant += jacobian.get(i, paramIndex) * (1.0 - obs / calc);
      }
    }
    PARALLEL_CRITICAL(der_set) {
      double der = m_der.get(activeParamIndex);
      m_der.set(activeParamIndex, der + determinant);
    }
    ++activeParamIndex;
  }
 
  PARALLEL_ATOMIC
  m_value += 2.0 * costVal;
}
 
void CostFuncPoisson::calculateHessian(API::IFunction &function, const API::FunctionDomain &domain,
                                       const API::FunctionValues &values) const {
  size_t numParams = function.nParams(); // number of parameters
  size_t numDataPoints = domain.size();  // number of data points
 
  Jacobian jacobian(numDataPoints, numParams);
  function.functionDeriv(domain, jacobian);
 
  size_t activeParamFirstIndex = 0; // The params are split into two halves and iterated through
  for (size_t paramIndex = 0; paramIndex < numParams; ++paramIndex) {
 
    if (!function.isActive(paramIndex))
      continue;
    size_t activeParamSecondIndex = 0; // The counterpart index
    double parameter = function.getParameter(paramIndex);
 
    double scalingFactor = 1e-4;
    if (parameter != 0.0) {
      scalingFactor *= parameter;
    }
 
    function.setParameter(paramIndex, parameter + scalingFactor);
    Jacobian jacobian2(numDataPoints, numParams);
    function.functionDeriv(domain, jacobian2);
    function.setParameter(paramIndex, parameter);
 
    for (size_t j = 0; j <= paramIndex; ++j) // over ~ half of parameters
    {
      if (!function.isActive(j))
        continue;
      double d = 0.0;
      for (size_t k = 0; k < numDataPoints; ++k) // over fitting data
      {
        double d2 = (jacobian2.get(k, j) - jacobian.get(k, j)) / scalingFactor;
 
        double calc = values.getCalculated(k);
        double obs = values.getFitData(k);
        if (calc <= absoluteCutOff) {
          d += std::numeric_limits<double>::infinity();
        } else {
          if (calc <= effectiveCutOff) {
            double constrainedCalc = calc - absoluteCutOff;
            d += d2 * (absoluteCutOff - effectiveCutOff) / (constrainedCalc * constrainedCalc);
            d += jacobian.get(k, paramIndex) * jacobian.get(k, j) * (effectiveCutOff - absoluteCutOff) * 2 /
                 (constrainedCalc * constrainedCalc * constrainedCalc);
          } else if (obs == 0.0) {
            d += d2;
          } else {
            d += d2 * (1.0 - obs / calc);
            d += jacobian.get(k, paramIndex) * jacobian.get(k, j) * obs / (calc * calc);
          }
        }
      }
      PARALLEL_CRITICAL(hessian_set) {
        double h = m_hessian.get(activeParamFirstIndex, activeParamSecondIndex) + d;
        m_hessian.set(activeParamFirstIndex, activeParamSecondIndex, h);
        if (activeParamFirstIndex != activeParamSecondIndex) {
          m_hessian.set(activeParamSecondIndex, activeParamFirstIndex, h);
        }
      }
      ++activeParamSecondIndex;
    }
    ++activeParamFirstIndex;
  }
}
 
} // namespace Mantid::CurveFitting::CostFunctions