CostFuncFitting.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/CostFuncFitting.h"
#include "MantidAPI/IConstraint.h"
#include "MantidCurveFitting/EigenJacobian.h"
#include "MantidCurveFitting/GSLFunctions.h"
#include "MantidCurveFitting/SeqDomain.h"
#include "MantidKernel/Exception.h"
 
#include <gsl/gsl_multifit_nlin.h>
#include <limits>
#include <utility>
 
namespace Mantid::CurveFitting::CostFunctions {
 
CostFuncFitting::CostFuncFitting()
    : m_dirtyVal(true), m_dirtyDeriv(true), m_dirtyHessian(true), m_includePenalty(true), m_value(0), m_pushed(false),
      m_pushedValue(false), m_ignoreInvalidData(false) {}
 
void CostFuncFitting::setDirty() {
  m_dirtyVal = true;
  m_dirtyDeriv = true;
  m_dirtyHessian = true;
}
 
double CostFuncFitting::getParameter(size_t i) const {
  checkValidity();
  return m_function->activeParameter(m_indexMap[i]);
}
 
void CostFuncFitting::setParameter(size_t i, const double &value) {
  checkValidity();
  m_function->setActiveParameter(m_indexMap[i], value);
  setDirty();
}
 
std::string CostFuncFitting::parameterName(size_t i) const {
  checkValidity();
  return m_function->parameterName(m_indexMap[i]);
}
 
size_t CostFuncFitting::nParams() const {
  checkValidity();
  return m_indexMap.size();
}
 
void CostFuncFitting::setFittingFunction(API::IFunction_sptr function, API::FunctionDomain_sptr domain,
                                         API::FunctionValues_sptr values) {
  m_function = std::move(function);
  m_domain = std::move(domain);
  m_values = std::move(values);
  updateValidateFitWeights();
  reset();
}
 
bool CostFuncFitting::isValid() const {
  if (m_function == API::IFunction_sptr()) {
    return false;
  }
  if (m_function->nParams() != m_numberFunParams) {
    reset();
  }
  auto nActive = m_indexMap.size();
  for (size_t i = 0, j = 0; i < m_numberFunParams; ++i) {
    if (m_function->isActive(i)) {
      if (j >= nActive || m_indexMap[j] != i) {
        reset();
        break;
      }
      ++j;
    }
  }
  return true;
}
 
void CostFuncFitting::checkValidity() const {
  if (!isValid()) {
    throw std::runtime_error("Fitting cost function isn't set");
  }
}
 
void CostFuncFitting::calActiveCovarianceMatrix(EigenMatrix &covar, double epsrel) {
  // construct the jacobian
  EigenJacobian J(*m_function, m_values->size());
  size_t na = this->nParams(); // number of active parameters
  const auto &j = J.getJ();
  assert(static_cast<size_t>(j.cols()) == na);
  covar.resize(na, na);
 
  // calculate the derivatives
  m_function->functionDeriv(*m_domain, J);
 
  // let the GSL to compute the covariance matrix
  EigenMatrix covarTr(covar.tr());
  EigenMatrix tempJTr;
  tempJTr = j.transpose();
  gsl_matrix_const_view tempJTrGSL = getGSLMatrixView_const(tempJTr.inspector());
  gsl_matrix_view covarTrGSL = getGSLMatrixView(covarTr.mutator());
 
  gsl_multifit_covar(&tempJTrGSL.matrix, epsrel, &covarTrGSL.matrix);
  covar = covarTr.tr();
}
 
void CostFuncFitting::calCovarianceMatrix(EigenMatrix &covar, double epsrel) {
  checkValidity();
  EigenMatrix c;
  calActiveCovarianceMatrix(c, epsrel);
 
  size_t np = m_function->nParams();
 
  bool isTransformationIdentity = true;
  for (size_t i = 0; i < np; ++i) {
    if (!m_function->isActive(i))
      continue;
    isTransformationIdentity =
        isTransformationIdentity && (m_function->activeParameter(i) == m_function->getParameter(i));
  }
 
  if (isTransformationIdentity) {
    // if the transformation is identity simply copy the matrix
    covar = c;
  } else {
    // else do the transformation
    EigenMatrix tm;
    calTransformationMatrixNumerically(tm);
    covar = tm.tr() * c * tm;
  }
}
 
void CostFuncFitting::calFittingErrors(const EigenMatrix &covar, double chi2) {
  checkValidity();
  size_t np = m_function->nParams();
  auto covarMatrix = std::shared_ptr<Kernel::Matrix<double>>(new Kernel::Matrix<double>(np, np));
  m_function->setCovarianceMatrix(covarMatrix);
  size_t ia = 0;
  for (size_t i = 0; i < np; ++i) {
    if (!m_function->isActive(i)) {
      m_function->setError(i, 0);
    } else {
      size_t ja = 0;
      for (size_t j = 0; j < np; ++j) {
        if (m_function->isActive(j)) {
          (*covarMatrix)[i][j] = covar.get(ia, ja);
          ++ja;
 
          if (ja >= covar.size2())
            break;
        }
      }
      double err = sqrt(covar.get(ia, ia));
      m_function->setError(i, err);
      ++ia;
 
      if (ia >= covar.size1())
        break;
    }
  }
  m_function->setCovarianceMatrix(covarMatrix);
  m_function->setReducedChiSquared(chi2 / static_cast<double>((m_values->size() - np)));
}
 
void CostFuncFitting::calTransformationMatrixNumerically(EigenMatrix &tm) {
  const double epsilon = std::numeric_limits<double>::epsilon() * 100;
  size_t np = m_function->nParams();
  size_t na = nParams();
  tm.resize(na, na);
  size_t ia = 0;
  for (size_t i = 0; i < np; ++i) {
    if (!m_function->isActive(i))
      continue;
    double p0 = m_function->getParameter(i);
    for (size_t j = 0; j < na; ++j) {
      double ap = getParameter(j);
      double step = ap == 0.0 ? epsilon : ap * epsilon;
      setParameter(j, ap + step);
      tm.set(ia, j, (m_function->getParameter(i) - p0) / step);
      setParameter(j, ap);
    }
    ++ia;
  }
}
 
void CostFuncFitting::applyTies() {
  if (m_function) {
    m_function->applyTies();
  }
}
 
void CostFuncFitting::reset() const {
  m_numberFunParams = m_function->nParams();
  m_indexMap.clear();
  for (size_t i = 0; i < m_numberFunParams; ++i) {
    if (m_function->isActive(i)) {
      m_indexMap.emplace_back(i);
    }
    API::IConstraint *c = m_function->getConstraint(i);
    if (c) {
      c->setParamToSatisfyConstraint();
    }
  }
  m_dirtyDeriv = true;
  m_dirtyHessian = true;
}
 
void CostFuncFitting::setParameters(const EigenVector &params) {
  auto np = nParams();
  if (np != params.size()) {
    throw Kernel::Exception::FitSizeWarning(params.size(), np);
  }
  for (size_t i = 0; i < np; ++i) {
    setParameter(i, params.get(i));
  }
  m_function->applyTies();
}
 
void CostFuncFitting::getParameters(EigenVector &params) const {
  auto np = nParams();
  if (params.size() != np) {
    params.resize(np);
  }
  for (size_t i = 0; i < np; ++i) {
    params.set(i, getParameter(i));
  }
}
 
double CostFuncFitting::val() const {
  if (!m_dirtyVal) {
    return m_value;
  }
 
  checkValidity();
 
  m_value = 0.0;
 
  auto seqDomain = std::dynamic_pointer_cast<SeqDomain>(m_domain);
 
  if (seqDomain) {
    seqDomain->additiveCostFunctionVal(*this);
  } else {
    if (!m_values) {
      throw std::runtime_error("CostFunction: undefined FunctionValues.");
    }
    addVal(m_domain, m_values);
  }
 
  // add penalty
  if (m_includePenalty) {
    for (size_t i = 0; i < m_function->nParams(); ++i) {
      if (!m_function->isActive(i))
        continue;
      API::IConstraint *c = m_function->getConstraint(i);
      if (c) {
        m_value += c->check();
      }
    }
  }
 
  m_dirtyVal = false;
  return m_value;
}
 
void CostFuncFitting::deriv(std::vector<double> &der) const {
  valDerivHessian(true, false);
 
  if (der.size() != nParams()) {
    der.resize(nParams());
  }
  for (size_t i = 0; i < nParams(); ++i) {
    der[i] = m_der.get(i);
  }
}
 
double CostFuncFitting::valAndDeriv(std::vector<double> &der) const {
  valDerivHessian(true, false);
 
  if (der.size() != nParams()) {
    der.resize(nParams());
  }
  for (size_t i = 0; i < nParams(); ++i) {
    der[i] = m_der.get(i);
  }
  return m_value;
}
 
double CostFuncFitting::valDerivHessian(bool evalDeriv, bool evalHessian) const {
  if (m_pushed || !evalDeriv) {
    return val();
  }
 
  if (!m_dirtyVal && !m_dirtyDeriv && !m_dirtyHessian) {
    return m_value;
  }
 
  const size_t numParams = nParams();
 
  checkValidity();
 
  m_value = 0.0;
 
  m_der.resize(numParams);
  m_der.zero();
  if (evalHessian) {
    m_hessian.resize(numParams, numParams);
    m_hessian.zero();
  }
 
  auto seqDomain = std::dynamic_pointer_cast<SeqDomain>(m_domain);
 
  if (seqDomain) {
    seqDomain->additiveCostFunctionValDerivHessian(*this, evalDeriv, evalHessian);
  } else {
    if (!m_values) {
      throw std::runtime_error("CostFunction: undefined FunctionValues.");
    }
    addValDerivHessian(m_function, m_domain, m_values, evalDeriv, evalHessian);
  }
 
  // Add constraints penalty
  size_t np = m_function->nParams();
  if (m_includePenalty) {
    for (size_t i = 0; i < np; ++i) {
      API::IConstraint *c = m_function->getConstraint(i);
      if (c) {
        m_value += c->check();
      }
    }
  }
  m_dirtyVal = false;
 
  if (evalDeriv) {
    if (m_includePenalty) {
      size_t i = 0;
      for (size_t ip = 0; ip < np; ++ip) {
        if (!m_function->isActive(ip))
          continue;
        API::IConstraint *c = m_function->getConstraint(ip);
        if (c) {
          double d = m_der.get(i) + c->checkDeriv();
          m_der.set(i, d);
        }
        ++i;
      }
    }
    m_dirtyDeriv = false;
 
    if (m_includePenalty) {
      size_t i = 0;
      for (size_t ip = 0; ip < np; ++ip) {
        if (!m_function->isActive(ip))
          continue;
        API::IConstraint *c = m_function->getConstraint(ip);
        if (c && !m_hessian.isEmpty()) {
          double d = m_hessian.get(i, i) + c->checkDeriv2();
          m_hessian.set(i, i, d);
        }
        ++i;
      }
    }
    // clear the dirty flag if hessian was actually calculated
    m_dirtyHessian = m_hessian.isEmpty();
  }
 
  return m_value;
}
 
const EigenVector &CostFuncFitting::getDeriv() const {
  if (m_pushed) {
    return m_der;
  }
  if (m_dirtyVal || m_dirtyDeriv || m_dirtyHessian) {
    valDerivHessian();
  }
  return m_der;
}
 
const EigenMatrix &CostFuncFitting::getHessian() const {
  if (m_pushed) {
    return m_hessian;
  }
  if (m_dirtyVal || m_dirtyDeriv || m_dirtyHessian) {
    valDerivHessian();
  }
  return m_hessian;
}
 
void CostFuncFitting::push() {
  if (m_pushed) {
    throw std::runtime_error("Cost Function: double push.");
  }
  // make sure we are not dirty
  m_pushedValue = valDerivHessian();
  getParameters(m_pushedParams);
  m_pushed = true;
}
 
void CostFuncFitting::pop() {
  if (!m_pushed) {
    throw std::runtime_error("Cost Function: empty stack.");
  }
  setParameters(m_pushedParams);
  m_value = m_pushedValue;
  m_pushed = false;
  m_dirtyVal = false;
  m_dirtyDeriv = false;
  m_dirtyHessian = false;
}
 
void CostFuncFitting::drop() {
  if (!m_pushed) {
    throw std::runtime_error("Cost Function: empty stack.");
  }
  m_pushed = false;
  setDirty();
}
 
void CostFuncFitting::validateNegativeFitWeights() {
  if (m_ignoreInvalidData) {
    return;
  }
 
  for (size_t i = 0; i < m_values->size(); i++) {
    if (m_values->getFitWeight(i) < 0) {
      throw std::runtime_error("Invalid data found at point=" + std::to_string(i) + " in fit weight.");
    }
  }
}
 
} // namespace Mantid::CurveFitting::CostFunctions