GSLFunctions.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 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/GSLFunctions.h"
#include "MantidAPI/IConstraint.h"
#include "MantidAPI/ICostFunction.h"
 
#include <MantidCurveFitting/CostFunctions/CostFuncPoisson.h>
 
namespace Mantid::CurveFitting {
 
namespace {
double leastSquaresLoss(const std::shared_ptr<API::FunctionValues> &values, const size_t index) {
  return (values->getCalculated(index) - values->getFitData(index)) * values->getFitWeight(index);
}
 
double poissonLoss(const std::shared_ptr<API::FunctionValues> &values, const size_t index) {
  return CostFunctions::PoissonLossLM::calculatePoissonLossLM(values->getFitData(index), values->getCalculated(index));
}
 
double leastSquaresWeight(const std::shared_ptr<API::FunctionValues> &values, const size_t index) {
  return values->getFitWeight(index);
}
 
double poissonWeight(const std::shared_ptr<API::FunctionValues> &values, const size_t index) {
  return CostFunctions::PoissonLossLM::calculateJacobianScaleFactor(values->getFitData(index),
                                                                    values->getCalculated(index));
}
} // namespace
 
int gsl_f(const gsl_vector *x, void *params, gsl_vector *f) {
  assert(x->data);
  auto *p = reinterpret_cast<struct GSL_FitData *>(params);
 
  // update function parameters
  size_t ia = 0;
  for (size_t i = 0; i < p->function->nParams(); ++i) {
    if (p->function->isActive(i)) {
      if (ia < x->size) {
        p->function->setActiveParameter(i, x->data[ia]);
        ++ia;
      } else {
        // The number of active parameters now exceeds the space
        // originally allocated
        throw Kernel::Exception::FitSizeWarning(x->size);
      }
    }
  }
  p->function->applyTies();
 
  auto values = std::dynamic_pointer_cast<API::FunctionValues>(p->costFunction->getValues());
  if (!values) {
    throw std::invalid_argument("FunctionValues expected");
  }
  p->function->function(*p->costFunction->getDomain(), *values);
 
  // Add penalty
  double penalty = 0.;
  for (size_t i = 0; i < p->function->nParams(); ++i) {
    API::IConstraint *c = p->function->getConstraint(i);
    if (c) {
      penalty += c->checkDeriv();
    }
  }
 
  size_t n = values->size() - 1;
  // add penalty to first and last point and every 10th point in between
  if (penalty != 0.0) {
    values->addToCalculated(0, penalty);
    values->addToCalculated(n, penalty);
 
    for (size_t i = 9; i < n; i += 10) {
      values->addToCalculated(i, penalty);
    }
  }
 
  // function() return calculated data values. Need to convert this values into
  // calculated-observed divided by error values used by GSL
 
  for (size_t i = 0; i < p->n; i++) {
    f->data[i] = p->m_loss(values, i);
  }
  return GSL_SUCCESS;
}
 
int gsl_df(const gsl_vector *x, void *params, gsl_matrix *J) {
 
  auto *p = reinterpret_cast<struct GSL_FitData *>(params);
  std::unique_ptr<gsl_matrix, decltype(&gsl_matrix_free)> J_tr(gsl_matrix_calloc(J->size2, J->size1), gsl_matrix_free);
  gsl_matrix_transpose_memcpy(J_tr.get(), J);
  EigenMatrix m(J_tr->size2, J_tr->size1);
  p->J.setJ(&m);
 
  // update function parameters
  if (x->data) {
    size_t ia = 0;
    for (size_t i = 0; i < p->function->nParams(); ++i) {
      if (p->function->isActive(i)) {
        p->function->setActiveParameter(i, x->data[ia]);
        ++ia;
      }
    }
  }
  p->function->applyTies();
 
  // calculate the Jacobian
  p->function->functionDeriv(*p->costFunction->getDomain(), p->J);
 
  // p->function->addPenaltyDeriv(&p->J);
  // add penalty
  size_t n = p->costFunction->getValues()->size() - 1;
  size_t ia = 0;
  for (size_t i = 0; i < p->function->nParams(); ++i) {
    if (!p->function->isActive(i))
      continue;
    API::IConstraint *c = p->function->getConstraint(i);
    if (c) {
      double penalty = c->checkDeriv2();
      if (penalty != 0.0) {
        double deriv = p->J.get(0, ia);
        p->J.set(0, ia, deriv + penalty);
        deriv = p->J.get(n, ia);
        p->J.set(n, ia, deriv + penalty);
 
        for (size_t j = 9; j < n; j += 10) {
          deriv = p->J.get(j, ia);
          p->J.set(j, ia, deriv + penalty);
        }
      }
    } // if (c)
    ++ia;
  }
 
  // functionDeriv() return derivatives of calculated data values. Need to
  // convert this values into
  // derivatives of calculated-observed divided by error values used by GSL
  auto values = std::dynamic_pointer_cast<API::FunctionValues>(p->costFunction->getValues());
  if (!values) {
    throw std::invalid_argument("FunctionValues expected");
  }
 
  EigenMatrix m_tr = m.tr();
  std::copy(&m_tr.mutator().data()[0], &m_tr.mutator().data()[J_tr->size1 * J_tr->size2], &J->data[0]);
  for (size_t iY = 0; iY < p->n; iY++) {
    const double weight = p->m_scaleFactor(values, iY);
    for (size_t iP = 0; iP < p->p; iP++) {
      J->data[iY * p->p + iP] *= weight;
    }
  }
 
  return GSL_SUCCESS;
}
 
int gsl_fdf(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J) {
  gsl_f(x, params, f);
  gsl_df(x, params, J);
  return GSL_SUCCESS;
}
 
GSL_FitData::GSL_FitData(const std::shared_ptr<CostFunctions::CostFuncFitting> &cf)
    : function(cf->getFittingFunction()), costFunction(cf) {
  gsl_set_error_handler_off();
 
  if (std::dynamic_pointer_cast<CostFunctions::CostFuncPoisson>(cf)) {
    this->m_loss = &poissonLoss;
    this->m_scaleFactor = &poissonWeight;
  } else {
    this->m_loss = &leastSquaresLoss;
    this->m_scaleFactor = &leastSquaresWeight;
  }
 
  // number of active parameters
  p = 0;
  for (size_t i = 0; i < function->nParams(); ++i) {
    if (function->isActive(i))
      ++p;
  }
 
  // number of fitting data
  n = cf->getValues()->size();
 
  bool functionFixed = false;
  if (p == 0) {
    p = 1;
    functionFixed = true;
  }
 
  // holdCalculatedJacobian =  gsl_matrix_alloc (n, p);
 
  initFuncParams = gsl_vector_alloc(p);
 
  if (functionFixed) {
    gsl_vector_set(initFuncParams, 0, 0.0);
  } else {
    size_t ia = 0;
    for (size_t i = 0; i < function->nParams(); ++i) {
      if (function->isActive(i)) {
        gsl_vector_set(initFuncParams, ia, function->activeParameter(i));
        ++ia;
      }
    }
  }
 
  int j = 0;
  for (size_t i = 0; i < function->nParams(); ++i) {
    if (function->isActive(i)) {
      J.m_index.emplace_back(j);
      j++;
    } else
      J.m_index.emplace_back(-1);
  }
}
 
GSL_FitData::~GSL_FitData() {
  // gsl_matrix_free(holdCalculatedJacobian);
  gsl_vector_free(initFuncParams);
}
 
} // namespace Mantid::CurveFitting