ConvolutionFit.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 +
#include "MantidAPI/AlgorithmManager.h"
#include "MantidAPI/CompositeFunction.h"
#include "MantidAPI/FunctionDomain1D.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/ITableWorkspace.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidAPI/Progress.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/TextAxis.h"
#include "MantidAPI/WorkspaceGroup.h"
 
#include "MantidCurveFitting/Algorithms/ConvolutionFit.h"
#include "MantidCurveFitting/Algorithms/QENSFitSequential.h"
#include "MantidCurveFitting/Algorithms/QENSFitSimultaneous.h"
 
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/MandatoryValidator.h"
#include "MantidKernel/StringContainsValidator.h"
#include "MantidKernel/VectorHelper.h"
 
#include <algorithm>
#include <cmath>
#include <utility>
 
namespace {
using namespace Mantid::API;
using namespace Mantid::Kernel;
using Mantid::MantidVec;
 
std::size_t numberOfFunctions(const IFunction_sptr &function, const std::string &functionName);
 
std::size_t numberOfFunctions(const CompositeFunction_sptr &composite, const std::string &functionName) {
  std::size_t count = 0;
  for (auto i = 0u; i < composite->nFunctions(); ++i)
    count += numberOfFunctions(composite->getFunction(i), functionName);
  return count;
}
 
std::size_t numberOfFunctions(const IFunction_sptr &function, const std::string &functionName) {
  const auto composite = std::dynamic_pointer_cast<CompositeFunction>(function);
  if (composite)
    return numberOfFunctions(composite, functionName);
  return function->name() == functionName ? 1 : 0;
}
 
bool containsFunction(const IFunction_sptr &function, const std::string &functionName);
 
bool containsFunction(const CompositeFunction_sptr &composite, const std::string &functionName) {
  for (auto i = 0u; i < composite->nFunctions(); ++i) {
    if (containsFunction(composite->getFunction(i), functionName))
      return true;
  }
  return false;
}
 
bool containsFunction(const IFunction_sptr &function, const std::string &functionName) {
  const auto composite = std::dynamic_pointer_cast<CompositeFunction>(function);
  if (function->name() == functionName)
    return true;
  else if (composite)
    return containsFunction(composite, functionName);
  return false;
}
 
template <typename T, typename F, typename... Ts>
std::vector<T, Ts...> transformVector(const std::vector<T, Ts...> &vec, F const &functor) {
  auto target = std::vector<T, Ts...>();
  target.reserve(vec.size());
  std::transform(vec.begin(), vec.end(), std::back_inserter(target), functor);
  return target;
}
 
template <typename T, typename F, typename... Ts>
std::vector<T, Ts...> combineVectors(const std::vector<T, Ts...> &vec, const std::vector<T, Ts...> &vec2,
                                     F const &combinator) {
  auto combined = std::vector<T, Ts...>();
  combined.reserve(vec.size());
  std::transform(vec.begin(), vec.end(), vec2.begin(), std::back_inserter(combined), combinator);
  return combined;
}
 
template <typename T, typename... Ts>
std::vector<T, Ts...> divideVectors(const std::vector<T, Ts...> &dividend, const std::vector<T, Ts...> &divisor) {
  return combineVectors(dividend, divisor, std::divides<T>());
}
 
template <typename T, typename... Ts>
std::vector<T, Ts...> addVectors(const std::vector<T, Ts...> &vec, const std::vector<T, Ts...> &vec2) {
  return combineVectors(vec, vec2, std::plus<T>());
}
 
template <typename T, typename... Ts>
std::vector<T, Ts...> multiplyVectors(const std::vector<T, Ts...> &vec, const std::vector<T, Ts...> &vec2) {
  return combineVectors(vec, vec2, std::multiplies<T>());
}
 
template <typename T, typename... Ts> std::vector<T, Ts...> squareVector(const std::vector<T, Ts...> &vec) {
  return transformVector(vec, VectorHelper::Squares<T>());
}
 
template <typename T, typename... Ts> std::vector<T, Ts...> squareRootVector(const std::vector<T, Ts...> &vec) {
  return transformVector(vec, static_cast<T (*)(T)>(sqrt));
}
 
IFunction_sptr extractFirstBackground(IFunction_sptr function);
 
IFunction_sptr extractFirstBackground(const CompositeFunction_sptr &composite) {
  for (auto i = 0u; i < composite->nFunctions(); ++i) {
    auto background = extractFirstBackground(composite->getFunction(i));
    if (background)
      return background;
  }
  return nullptr;
}
 
IFunction_sptr extractFirstBackground(IFunction_sptr function) {
  auto composite = std::dynamic_pointer_cast<CompositeFunction>(function);
 
  if (composite)
    return extractFirstBackground(composite);
  else if (function->category() == "Background")
    return function;
  return nullptr;
}
 
std::string extractBackgroundType(IFunction_sptr function) {
  auto background = extractFirstBackground(std::move(function));
  if (!background)
    return "None";
 
  auto backgroundType = background->name();
  auto position = backgroundType.rfind("Background");
 
  if (position != std::string::npos)
    backgroundType = backgroundType.substr(0, position);
 
  if (background->isFixed(0))
    backgroundType = "Fixed " + backgroundType;
  else
    backgroundType = "Fit " + backgroundType;
  return backgroundType;
}
 
std::vector<std::size_t> searchForFitParameters(const std::string &suffix, const ITableWorkspace_sptr &tableWorkspace) {
  auto indices = std::vector<std::size_t>();
 
  for (auto i = 0u; i < tableWorkspace->columnCount(); ++i) {
    auto name = tableWorkspace->getColumn(i)->name();
    auto position = name.rfind(suffix);
    if (position != std::string::npos && position + suffix.size() == name.size())
      indices.emplace_back(i);
  }
  return indices;
}
 
std::pair<MantidVec, MantidVec> calculateEISFAndError(const MantidVec &height, const MantidVec &heightError,
                                                      const MantidVec &amplitude, const MantidVec &amplitudeError) {
  auto total = addVectors(height, amplitude);
  auto eisfY = divideVectors(height, total);
 
  auto heightESq = squareVector(heightError);
 
  auto ampErrSq = squareVector(amplitudeError);
  auto totalErr = addVectors(heightESq, ampErrSq);
 
  auto heightYSq = squareVector(height);
  auto totalSq = squareVector(total);
 
  auto errOverTotalSq = divideVectors(totalErr, totalSq);
  auto heightESqOverYSq = divideVectors(heightESq, heightYSq);
 
  auto sqrtESqOverYSq = squareRootVector(heightESqOverYSq);
  auto eisfYSumRoot = multiplyVectors(eisfY, sqrtESqOverYSq);
 
  return {eisfY, addVectors(eisfYSumRoot, errOverTotalSq)};
}
 
void addEISFToTable(ITableWorkspace_sptr &tableWs) {
  // Get height data from parameter table
  const auto height = searchForFitParameters("Height", tableWs).at(0);
  const auto heightErr = searchForFitParameters("Height_Err", tableWs).at(0);
  const auto heightY = tableWs->getColumn(height)->numeric_fill();
  const auto heightE = tableWs->getColumn(heightErr)->numeric_fill();
 
  // Get amplitude column names
  const auto ampIndices = searchForFitParameters("Amplitude", tableWs);
  const auto ampErrorIndices = searchForFitParameters("Amplitude_Err", tableWs);
 
  // For each lorentzian, calculate EISF
  auto maxSize = ampIndices.size();
  if (ampErrorIndices.size() > maxSize)
    maxSize = ampErrorIndices.size();
 
  for (auto i = 0u; i < maxSize; ++i) {
    // Get amplitude from column in table workspace
    auto ampY = tableWs->getColumn(ampIndices[i])->numeric_fill();
    auto ampErr = tableWs->getColumn(ampErrorIndices[i])->numeric_fill();
    auto eisfAndError = calculateEISFAndError(heightY, heightE, ampY, ampErr);
 
    // Append the calculated values to the table workspace
    auto ampName = tableWs->getColumn(ampIndices[i])->name();
    auto ampErrorName = tableWs->getColumn(ampErrorIndices[i])->name();
    auto columnName = ampName.substr(0, (ampName.size() - std::string("Amplitude").size()));
    columnName += "EISF";
    auto errorColumnName = ampErrorName.substr(0, (ampErrorName.size() - std::string("Amplitude_Err").size()));
    errorColumnName += "EISF_Err";
 
    tableWs->addColumn("double", columnName);
    tableWs->addColumn("double", errorColumnName);
    auto maxEisf = eisfAndError.first.size();
    if (eisfAndError.second.size() > maxEisf) {
      maxEisf = eisfAndError.second.size();
    }
 
    auto col = tableWs->getColumn(columnName);
    auto errCol = tableWs->getColumn(errorColumnName);
    for (auto j = 0u; j < maxEisf; j++) {
      col->cell<double>(j) = eisfAndError.first.at(j);
      errCol->cell<double>(j) = eisfAndError.second.at(j);
    }
  }
}
} // namespace
 
namespace Mantid::CurveFitting::Algorithms {
 
using namespace API;
using namespace Kernel;
 
//----------------------------------------------------------------------------------------------
 
template <> const std::string ConvolutionFit<QENSFitSequential>::name() const { return "ConvolutionFitSequential"; }
 
template <> const std::string ConvolutionFit<QENSFitSimultaneous>::name() const { return "ConvolutionFitSimultaneous"; }
 
template <typename Base> const std::string ConvolutionFit<Base>::name() const { return "ConvolutionFit"; }
 
template <typename Base> int ConvolutionFit<Base>::version() const { return 1; }
 
template <typename Base> const std::string ConvolutionFit<Base>::category() const { return "Workflow\\MIDAS"; }
 
template <> const std::string ConvolutionFit<QENSFitSequential>::summary() const {
  return "Performs a sequential fit for a convolution workspace";
}
 
template <> const std::string ConvolutionFit<QENSFitSimultaneous>::summary() const {
  return "Performs a simultaneous fit across convolution workspaces";
}
 
template <typename Base> const std::string ConvolutionFit<Base>::summary() const {
  return "Performs a convolution fit";
}
 
template <> const std::vector<std::string> ConvolutionFit<QENSFitSequential>::seeAlso() const {
  return {"QENSFitSequential"};
}
 
template <> const std::vector<std::string> ConvolutionFit<QENSFitSimultaneous>::seeAlso() const {
  return {"QENSFitSimultaneous"};
}
 
template <typename Base> const std::vector<std::string> ConvolutionFit<Base>::seeAlso() const { return {}; }
 
template <typename Base> std::map<std::string, std::string> ConvolutionFit<Base>::validateInputs() {
  auto errors = Base::validateInputs();
  IFunction_sptr function = Base::getProperty("Function");
  if (!containsFunction(function, "Convolution") || !containsFunction(function, "Resolution"))
    errors["Function"] = "Function provided does not contain convolution with "
                         "a resolution function.";
  return errors;
}
 
template <typename Base> bool ConvolutionFit<Base>::isFitParameter(const std::string &name) const {
  bool isBackgroundParameter = name.rfind("A0") != std::string::npos || name.rfind("A1") != std::string::npos;
  return name.rfind("Centre") == std::string::npos && !isBackgroundParameter;
}
 
template <typename Base> bool ConvolutionFit<Base>::throwIfElasticQConversionFails() const { return true; }
 
template <typename Base>
ITableWorkspace_sptr ConvolutionFit<Base>::processParameterTable(ITableWorkspace_sptr parameterTable) {
  IFunction_sptr function = Base::getProperty("Function");
  m_deltaUsed = containsFunction(function, "DeltaFunction");
  if (m_deltaUsed)
    addEISFToTable(parameterTable);
  return parameterTable;
}
 
template <typename Base> std::map<std::string, std::string> ConvolutionFit<Base>::getAdditionalLogStrings() const {
  IFunction_sptr function = Base::getProperty("Function");
  auto logs = Base::getAdditionalLogStrings();
  logs["delta_function"] = m_deltaUsed ? "true" : "false";
  logs["background"] = extractBackgroundType(function);
  return logs;
}
 
template <typename Base> std::map<std::string, std::string> ConvolutionFit<Base>::getAdditionalLogNumbers() const {
  auto logs = Base::getAdditionalLogNumbers();
  IFunction_sptr function = Base::getProperty("Function");
  logs["lorentzians"] = boost::lexical_cast<std::string>(numberOfFunctions(function, "Lorentzian"));
  return logs;
}
 
template <typename Base> std::vector<std::string> ConvolutionFit<Base>::getFitParameterNames() const {
  auto names = Base::getFitParameterNames();
  if (m_deltaUsed)
    names.emplace_back("EISF");
  return names;
}
 
// Register the algorithms into the AlgorithmFactory
template class ConvolutionFit<QENSFitSequential>;
template class ConvolutionFit<QENSFitSimultaneous>;
 
using ConvolutionFitSequential = ConvolutionFit<QENSFitSequential>;
using ConvolutionFitSimultaneous = ConvolutionFit<QENSFitSimultaneous>;
 
DECLARE_ALGORITHM(ConvolutionFitSequential)
DECLARE_ALGORITHM(ConvolutionFitSimultaneous)
 
} // namespace Mantid::CurveFitting::Algorithms