IMWDomainCreator.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/IMWDomainCreator.h"
#include "MantidCurveFitting/Functions/Convolution.h"
#include "MantidCurveFitting/Jacobian.h"
#include "MantidCurveFitting/ParameterEstimator.h"
#include "MantidCurveFitting/SeqDomain.h"
 
#include "MantidAPI/CompositeFunction.h"
#include "MantidAPI/FunctionDomain1D.h"
#include "MantidAPI/FunctionProperty.h"
#include "MantidAPI/FunctionValues.h"
#include "MantidAPI/IEventWorkspace.h"
#include "MantidAPI/IFunctionMW.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceProperty.h"
 
#include "MantidAPI/TextAxis.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/EmptyValues.h"
#include "MantidKernel/Matrix.h"
 
#include <algorithm>
#include <boost/math/distributions/students_t.hpp>
 
namespace Mantid::CurveFitting {
 
namespace {
bool greaterIsLess(double x1, double x2) { return x1 > x2; }
} // namespace
 
using namespace Kernel;
using API::MatrixWorkspace;
using API::Workspace;
 
IMWDomainCreator::IMWDomainCreator(Kernel::IPropertyManager *fit, const std::string &workspacePropertyName,
                                   IMWDomainCreator::DomainType domainType)
    : API::IDomainCreator(fit, std::vector<std::string>(1, workspacePropertyName), domainType), m_workspaceIndex(-1),
      m_startX(EMPTY_DBL()), m_endX(EMPTY_DBL()), m_startIndex(0) {
  if (m_workspacePropertyNames.empty()) {
    throw std::runtime_error("Cannot create FitMW: no workspace given");
  }
  m_workspacePropertyName = m_workspacePropertyNames[0];
}
 
void IMWDomainCreator::setParameters() const {
  // if property manager is set overwrite any set parameters
  if (m_manager) {
 
    // get the workspace
    API::Workspace_sptr ws = m_manager->getProperty(m_workspacePropertyName);
    m_matrixWorkspace = std::dynamic_pointer_cast<API::MatrixWorkspace>(ws);
    if (!m_matrixWorkspace) {
      throw std::invalid_argument("InputWorkspace must be a MatrixWorkspace.");
    }
 
    int index = m_manager->getProperty(m_workspaceIndexPropertyName);
    m_workspaceIndex = static_cast<size_t>(index);
    m_startX = m_manager->getProperty(m_startXPropertyName);
    m_endX = m_manager->getProperty(m_endXPropertyName);
  }
}
 
void IMWDomainCreator::declareDatasetProperties(const std::string &suffix, bool addProp) {
  m_workspaceIndexPropertyName = "WorkspaceIndex" + suffix;
  m_startXPropertyName = "StartX" + suffix;
  m_endXPropertyName = "EndX" + suffix;
 
  if (addProp && !m_manager->existsProperty(m_workspaceIndexPropertyName)) {
    auto mustBePositive = std::make_shared<BoundedValidator<int>>();
    mustBePositive->setLower(0);
    declareProperty(new PropertyWithValue<int>(m_workspaceIndexPropertyName, 0, mustBePositive),
                    "The Workspace Index to fit in the input workspace");
    declareProperty(new PropertyWithValue<double>(m_startXPropertyName, EMPTY_DBL()),
                    "A value of x in, or on the low x boundary of, the first bin to "
                    "include in\n"
                    "the fit (default lowest value of x)");
    declareProperty(new PropertyWithValue<double>(m_endXPropertyName, EMPTY_DBL()),
                    "A value in, or on the high x boundary of, the last bin the fitting "
                    "range\n"
                    "(default the highest value of x)");
  }
}
 
std::pair<size_t, size_t> IMWDomainCreator::getXInterval() const {
 
  const auto &X = m_matrixWorkspace->x(m_workspaceIndex);
  if (X.empty()) {
    throw std::runtime_error("Workspace contains no data.");
  }
 
  setParameters();
 
  // From points to the first occurrence of StartX in the workspace interval.
  // End points to the last occurrence of EndX in the workspace interval.
  // Find the fitting interval: from -> to
  Mantid::MantidVec::const_iterator from;
  Mantid::MantidVec::const_iterator to;
 
  bool isXAscending = X.front() < X.back();
 
  if (m_startX == EMPTY_DBL() && m_endX == EMPTY_DBL()) {
    m_startX = X.front();
    from = X.begin();
    m_endX = X.back();
    to = X.end();
  } else if (m_startX == EMPTY_DBL() || m_endX == EMPTY_DBL()) {
    throw std::invalid_argument("Both StartX and EndX must be given to set fitting interval.");
  } else if (isXAscending) {
    if (m_startX > m_endX) {
      std::swap(m_startX, m_endX);
    }
    from = std::lower_bound(X.begin(), X.end(), m_startX);
    to = std::upper_bound(from, X.end(), m_endX);
  } else { // x is descending
    if (m_startX < m_endX) {
      std::swap(m_startX, m_endX);
    }
    from = std::lower_bound(X.begin(), X.end(), m_startX, greaterIsLess);
    to = std::upper_bound(from, X.end(), m_endX, greaterIsLess);
  }
 
  // Check whether the fitting interval defined by StartX and EndX is 0.
  // This occurs when StartX and EndX are both less than the minimum workspace
  // x-value or greater than the maximum workspace x-value.
  if (to - from == 0) {
    throw std::invalid_argument("StartX and EndX values do not capture a range "
                                "within the workspace interval.");
  }
 
  if (m_matrixWorkspace->isHistogramData()) {
    if (X.end() == to) {
      to = X.end() - 1;
    }
  }
  return std::make_pair(std::distance(X.begin(), from), std::distance(X.begin(), to));
}
 
size_t IMWDomainCreator::getDomainSize() const {
  setParameters();
  size_t startIndex, endIndex;
  std::tie(startIndex, endIndex) = getXInterval();
  return endIndex - startIndex;
}
 
void IMWDomainCreator::initFunction(API::IFunction_sptr function) {
  setParameters();
  if (!function) {
    throw std::runtime_error("Cannot initialize empty function.");
  }
  function->setWorkspace(m_matrixWorkspace);
  function->setMatrixWorkspace(m_matrixWorkspace, m_workspaceIndex, m_startX, m_endX);
  setInitialValues(*function);
}
 
API::MatrixWorkspace_sptr IMWDomainCreator::createEmptyResultWS(const size_t nhistograms, const size_t nyvalues) {
  size_t nxvalues(nyvalues);
  if (m_matrixWorkspace->isHistogramData())
    nxvalues += 1;
  API::MatrixWorkspace_sptr ws =
      API::WorkspaceFactory::Instance().create("Workspace2D", nhistograms, nxvalues, nyvalues);
  ws->setTitle("");
  ws->setYUnitLabel(m_matrixWorkspace->YUnitLabel());
  ws->setYUnit(m_matrixWorkspace->YUnit());
  ws->getAxis(0)->unit() = m_matrixWorkspace->getAxis(0)->unit();
  auto tAxis = std::make_unique<API::TextAxis>(nhistograms);
  ws->replaceAxis(1, std::move(tAxis));
 
  auto &inputX = m_matrixWorkspace->x(m_workspaceIndex);
  auto &inputY = m_matrixWorkspace->y(m_workspaceIndex);
  auto &inputE = m_matrixWorkspace->e(m_workspaceIndex);
  // X values for all
  for (size_t i = 0; i < nhistograms; i++) {
    ws->mutableX(i).assign(inputX.begin() + m_startIndex, inputX.begin() + m_startIndex + nxvalues);
  }
  // Data values for the first histogram
  ws->mutableY(0).assign(inputY.begin() + m_startIndex, inputY.begin() + m_startIndex + nyvalues);
  ws->mutableE(0).assign(inputE.begin() + m_startIndex, inputE.begin() + m_startIndex + nyvalues);
 
  return ws;
}
 
void IMWDomainCreator::setInitialValues(API::IFunction &function) {
  auto domain = m_domain.lock();
  auto values = m_values.lock();
  if (domain && values) {
    ParameterEstimator::estimate(function, *domain, *values);
  }
}
 
std::shared_ptr<API::Workspace> IMWDomainCreator::createOutputWorkspace(
    const std::string &baseName, API::IFunction_sptr function, std::shared_ptr<API::FunctionDomain> domain,
    std::shared_ptr<API::FunctionValues> values, const std::string &outputWorkspacePropertyName) {
  if (!values) {
    throw std::logic_error("FunctionValues expected");
  }
 
  // Compile list of functions to output. The top-level one is first
  std::list<API::IFunction_sptr> functionsToDisplay(1, function);
  if (m_outputCompositeMembers) {
    appendCompositeFunctionMembers(functionsToDisplay, function);
  }
 
  // Nhist = Data histogram, Difference Histogram + nfunctions
  const size_t nhistograms = functionsToDisplay.size() + 2;
  const size_t nyvalues = values->size();
  auto ws = createEmptyResultWS(nhistograms, nyvalues);
  // The workspace was constructed with a TextAxis
  API::TextAxis *textAxis = static_cast<API::TextAxis *>(ws->getAxis(1));
  textAxis->setLabel(0, "Data");
  textAxis->setLabel(1, "Calc");
  textAxis->setLabel(2, "Diff");
 
  // Add each calculated function
  auto iend = functionsToDisplay.end();
  size_t wsIndex(1); // Zero reserved for data
  for (auto it = functionsToDisplay.begin(); it != iend; ++it) {
    if (wsIndex > 2)
      textAxis->setLabel(wsIndex, (*it)->name());
    addFunctionValuesToWS(*it, ws, wsIndex, domain, values);
    if (it == functionsToDisplay.begin())
      wsIndex += 2; // Skip difference histogram for now
    else
      ++wsIndex;
  }
 
  // Set the difference spectrum
  const auto &Ycal = ws->y(1);
  auto &Diff = ws->mutableY(2);
  const size_t nData = values->size();
  for (size_t i = 0; i < nData; ++i) {
    if (values->getFitWeight(i) != 0.0) {
      Diff[i] = values->getFitData(i) - Ycal[i];
    } else {
      Diff[i] = 0.0;
    }
  }
 
  if (!outputWorkspacePropertyName.empty()) {
    declareProperty(new API::WorkspaceProperty<MatrixWorkspace>(outputWorkspacePropertyName, "", Direction::Output),
                    "Name of the output Workspace holding resulting simulated spectrum");
    m_manager->setPropertyValue(outputWorkspacePropertyName, baseName + "Workspace");
    m_manager->setProperty(outputWorkspacePropertyName, ws);
  }
 
  // If the input is a not an EventWorkspace and is a distrubution, then convert
  // the output also to a distribution
  if (!std::dynamic_pointer_cast<Mantid::API::IEventWorkspace>(m_matrixWorkspace)) {
    if (m_matrixWorkspace->isDistribution()) {
      ws->setDistribution(true);
    }
  }
 
  return ws;
}
 
void IMWDomainCreator::appendCompositeFunctionMembers(std::list<API::IFunction_sptr> &functionList,
                                                      const API::IFunction_sptr &function) const {
  // if function is a Convolution then output of convolved model's mebers may be
  // required
  if (m_convolutionCompositeMembers && std::dynamic_pointer_cast<Functions::Convolution>(function)) {
    appendConvolvedCompositeFunctionMembers(functionList, function);
  } else {
    const auto compositeFn = std::dynamic_pointer_cast<API::CompositeFunction>(function);
    if (!compositeFn)
      return;
 
    const size_t nlocals = compositeFn->nFunctions();
    for (size_t i = 0; i < nlocals; ++i) {
      auto localFunction = compositeFn->getFunction(i);
      auto localComposite = std::dynamic_pointer_cast<API::CompositeFunction>(localFunction);
      if (localComposite)
        appendCompositeFunctionMembers(functionList, localComposite);
      else
        functionList.insert(functionList.end(), localFunction);
    }
  }
}
 
void IMWDomainCreator::appendConvolvedCompositeFunctionMembers(std::list<API::IFunction_sptr> &functionList,
                                                               const API::IFunction_sptr &function) const {
  std::shared_ptr<Functions::Convolution> convolution = std::dynamic_pointer_cast<Functions::Convolution>(function);
 
  const auto compositeFn = std::dynamic_pointer_cast<API::CompositeFunction>(convolution->getFunction(1));
  if (!compositeFn) {
    functionList.insert(functionList.end(), convolution);
  } else {
    auto resolution = convolution->getFunction(0);
    const size_t nlocals = compositeFn->nFunctions();
    for (size_t i = 0; i < nlocals; ++i) {
      auto localFunction = compositeFn->getFunction(i);
      std::shared_ptr<Functions::Convolution> localConvolution = std::make_shared<Functions::Convolution>();
      localConvolution->addFunction(resolution);
      localConvolution->addFunction(localFunction);
      functionList.insert(functionList.end(), localConvolution);
    }
  }
}
 
void IMWDomainCreator::addFunctionValuesToWS(const API::IFunction_sptr &function,
                                             std::shared_ptr<API::MatrixWorkspace> &ws, const size_t wsIndex,
                                             const std::shared_ptr<API::FunctionDomain> &domain,
                                             const std::shared_ptr<API::FunctionValues> &resultValues) const {
  const size_t nData = resultValues->size();
  resultValues->zeroCalculated();
  // Confidence bands are calculated based on the example in
  // www.astro.rug.nl/software/kapteyn/kmpfittutorial.html#confidence-and-prediction-intervals,
  // which references J.Wolberg, Data Analysis Using the Method of Least Squares, 2006, Springer.
  // Here we assume a confidence band of 1 sigma
  double sigma = 1;
  double prob = std::erf(sigma / sqrt(2));
  // critical value for t distribution
  double alpha = (1 + prob) / 2;
 
  // Function value
  function->function(*domain, *resultValues);
 
  size_t nParams = function->nParams();
 
  // the function should contain the parameter's covariance matrix
  auto covar = function->getCovarianceMatrix();
  bool hasErrors = false;
  if (!covar) {
    for (size_t j = 0; j < nParams; ++j) {
      if (function->getError(j) != 0.0) {
        hasErrors = true;
        break;
      }
    }
  }
 
  if (covar || hasErrors) {
    // and errors
    Jacobian J(nData, nParams);
    try {
      function->functionDeriv(*domain, J);
    } catch (...) {
      function->calNumericalDeriv(*domain, J);
    }
    if (covar) {
      // if the function has a covariance matrix attached - use it for the
      // errors
      const Kernel::Matrix<double> &C = *covar;
      // The formula is E = J * C * J^T
      // We don't do full 3-matrix multiplication because we only need the
      // diagonals of E
      std::vector<double> E(nData);
      for (size_t k = 0; k < nData; ++k) {
        double s = 0.0;
        for (size_t i = 0; i < nParams; ++i) {
          double tmp = J.get(k, i);
          s += C[i][i] * tmp * tmp;
          for (size_t j = i + 1; j < nParams; ++j) {
            s += J.get(k, i) * C[i][j] * J.get(k, j) * 2;
          }
        }
        E[k] = s;
      }
 
      size_t dof = nData - nParams;
      auto &yValues = ws->mutableY(wsIndex);
      auto &eValues = ws->mutableE(wsIndex);
      double T = 1.0;
      if (dof != 0) {
        boost::math::students_t dist(static_cast<double>(dof));
        T = boost::math::quantile(dist, alpha);
      }
      for (size_t i = 0; i < nData; i++) {
        yValues[i] = resultValues->getCalculated(i);
        eValues[i] = T * std::sqrt(E[i]);
      }
    } else {
      // otherwise use the parameter errors which is OK for uncorrelated
      // parameters
      auto &yValues = ws->mutableY(wsIndex);
      auto &eValues = ws->mutableE(wsIndex);
      for (size_t i = 0; i < nData; i++) {
        yValues[i] = resultValues->getCalculated(i);
        double err = 0.0;
        for (size_t j = 0; j < nParams; ++j) {
          double d = J.get(i, j) * function->getError(j);
          err += d * d;
        }
        eValues[i] = std::sqrt(err);
      }
    }
  } else {
    // No errors
    auto &yValues = ws->mutableY(wsIndex);
    for (size_t i = 0; i < nData; i++) {
      yValues[i] = resultValues->getCalculated(i);
    }
  }
}
 
} // namespace Mantid::CurveFitting