PawleyFunction.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 "MantidCurveFitting/Functions/PawleyFunction.h"
 
#include "MantidAPI/Axis.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/MatrixWorkspace.h"
 
#include "MantidCurveFitting/Constraints/BoundaryConstraint.h"
 
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/UnitConversion.h"
#include "MantidKernel/UnitFactory.h"
 
#include <boost/algorithm/string.hpp>
#include <memory>
 
namespace Mantid::CurveFitting::Functions {
 
using namespace CurveFitting;
using namespace Constraints;
 
DECLARE_FUNCTION(PawleyParameterFunction)
 
using namespace API;
 
using namespace Geometry;
 
using namespace Kernel;
 
PawleyParameterFunction::PawleyParameterFunction()
    : ParamFunction(), m_latticeSystem(PointGroup::LatticeSystem::Triclinic), m_profileFunctionCenterParameterName() {}
 
void PawleyParameterFunction::setAttribute(const std::string &attName, const Attribute &attValue) {
  if (attName == "LatticeSystem") {
    setLatticeSystem(attValue.asString());
  } else if (attName == "ProfileFunction") {
    setProfileFunction(attValue.asString());
  }
 
  ParamFunction::setAttribute(attName, attValue);
}
 
PointGroup::LatticeSystem PawleyParameterFunction::getLatticeSystem() const { return m_latticeSystem; }
 
UnitCell PawleyParameterFunction::getUnitCellFromParameters() const {
  switch (m_latticeSystem) {
  case PointGroup::LatticeSystem::Cubic: {
    double a = getParameter("a");
    double aErr = getError(0);
    UnitCell uc(a, a, a);
    uc.setError(aErr, aErr, aErr, 0.0, 0.0, 0.0);
    return uc;
  }
  case PointGroup::LatticeSystem::Tetragonal: {
    double a = getParameter("a");
    double aErr = getError(0);
    UnitCell uc(a, a, getParameter("c"));
    uc.setError(aErr, aErr, getError(1), 0.0, 0.0, 0.0);
    return uc;
  }
  case PointGroup::LatticeSystem::Hexagonal: {
    double a = getParameter("a");
    double aErr = getError(0);
    UnitCell uc(a, a, getParameter("c"), 90, 90, 120);
    uc.setError(aErr, aErr, getError(1), 0.0, 0.0, 0.0);
    return uc;
  }
  case PointGroup::LatticeSystem::Rhombohedral: {
    double a = getParameter("a");
    double alpha = getParameter("Alpha");
    double aErr = getError(0);
    double alphaErr = getError(1);
    UnitCell uc(a, a, a, alpha, alpha, alpha);
    uc.setError(aErr, aErr, aErr, alphaErr, alphaErr, alphaErr);
    return uc;
  }
  case PointGroup::LatticeSystem::Orthorhombic: {
    UnitCell uc(getParameter("a"), getParameter("b"), getParameter("c"));
    uc.setError(getError(0), getError(1), getError(2), 0.0, 0.0, 0.0);
    return uc;
  }
  case PointGroup::LatticeSystem::Monoclinic: {
    UnitCell uc(getParameter("a"), getParameter("b"), getParameter("c"), 90, getParameter("Beta"), 90);
    uc.setError(getError(0), getError(1), getError(2), 0.0, getError(3), 0.0);
    return uc;
  }
  case PointGroup::LatticeSystem::Triclinic: {
    UnitCell uc(getParameter("a"), getParameter("b"), getParameter("c"), getParameter("Alpha"), getParameter("Beta"),
                getParameter("Gamma"));
    uc.setError(getError(0), getError(1), getError(2), getError(3), getError(4), getError(5));
    return uc;
  }
  }
 
  return UnitCell();
}
 
void PawleyParameterFunction::setParametersFromUnitCell(const UnitCell &cell) {
  // Parameter "a" exists in all crystal systems.
  setParameter("a", cell.a());
 
  try {
    setParameter("b", cell.b());
  } catch (const std::invalid_argument &) {
    // do nothing.
  }
 
  try {
    setParameter("c", cell.c());
  } catch (const std::invalid_argument &) {
    // do nothing
  }
 
  try {
    setParameter("Alpha", cell.alpha());
  } catch (const std::invalid_argument &) {
    // do nothing.
  }
  try {
    setParameter("Beta", cell.beta());
  } catch (const std::invalid_argument &) {
    // do nothing.
  }
  try {
    setParameter("Gamma", cell.gamma());
  } catch (const std::invalid_argument &) {
    // do nothing.
  }
}
 
void PawleyParameterFunction::function(const FunctionDomain &domain, FunctionValues &values) const {
  UNUSED_ARG(domain);
  UNUSED_ARG(values);
}
 
void PawleyParameterFunction::functionDeriv(const FunctionDomain &domain, Jacobian &jacobian) {
  UNUSED_ARG(domain)
  UNUSED_ARG(jacobian);
}
 
void PawleyParameterFunction::init() {
  declareAttribute("LatticeSystem", IFunction::Attribute("Triclinic"));
  declareAttribute("ProfileFunction", IFunction::Attribute("Gaussian"));
 
  setLatticeSystem("Triclinic");
  setProfileFunction("Gaussian");
}
 
void PawleyParameterFunction::setProfileFunction(const std::string &profileFunction) {
  IPeakFunction_sptr peakFunction =
      std::dynamic_pointer_cast<IPeakFunction>(FunctionFactory::Instance().createFunction(profileFunction));
 
  if (!peakFunction) {
    throw std::invalid_argument("PawleyFunction can only use IPeakFunctions to "
                                "calculate peak profiles.");
  }
 
  setCenterParameterNameFromFunction(peakFunction);
}
 
void PawleyParameterFunction::setLatticeSystem(const std::string &latticeSystem) {
  m_latticeSystem = Geometry::getLatticeSystemFromString(latticeSystem);
 
  createLatticeSystemParameters(m_latticeSystem);
}
 
void PawleyParameterFunction::createLatticeSystemParameters(PointGroup::LatticeSystem latticeSystem) {
 
  clearAllParameters();
  switch (latticeSystem) {
  case PointGroup::LatticeSystem::Cubic:
    declareParameter("a", 1.0);
    addLengthConstraint("a");
    break;
 
  case PointGroup::LatticeSystem::Hexagonal:
  case PointGroup::LatticeSystem::Tetragonal:
    declareParameter("a", 1.0);
    declareParameter("c", 1.0);
    addLengthConstraint("a");
    addLengthConstraint("c");
    break;
 
  case PointGroup::LatticeSystem::Orthorhombic:
    declareParameter("a", 1.0);
    declareParameter("b", 1.0);
    declareParameter("c", 1.0);
    addLengthConstraint("a");
    addLengthConstraint("b");
    addLengthConstraint("c");
    break;
 
  case PointGroup::LatticeSystem::Monoclinic:
    declareParameter("a", 1.0);
    declareParameter("b", 1.0);
    declareParameter("c", 1.0);
    addLengthConstraint("a");
    addLengthConstraint("b");
    addLengthConstraint("c");
 
    declareParameter("Beta", 90.0);
    addAngleConstraint("Beta");
    break;
 
  case PointGroup::LatticeSystem::Rhombohedral:
    declareParameter("a", 1.0);
    declareParameter("Alpha", 90.0);
    addLengthConstraint("a");
    addAngleConstraint("Alpha");
    break;
 
  default:
    // triclinic
    declareParameter("a", 1.0);
    declareParameter("b", 1.0);
    declareParameter("c", 1.0);
    addLengthConstraint("a");
    addLengthConstraint("b");
    addLengthConstraint("c");
 
    declareParameter("Alpha", 90.0);
    declareParameter("Beta", 90.0);
    declareParameter("Gamma", 90.0);
    addAngleConstraint("Alpha");
    addAngleConstraint("Beta");
    addAngleConstraint("Gamma");
    break;
  }
 
  declareParameter("ZeroShift", 0.0);
}
 
void PawleyParameterFunction::addLengthConstraint(const std::string &parameterName) {
  auto cellEdgeConstraint = std::make_unique<BoundaryConstraint>(this, parameterName, 0.0, true);
  cellEdgeConstraint->setPenaltyFactor(1e12);
  addConstraint(std::move(cellEdgeConstraint));
}
 
void PawleyParameterFunction::addAngleConstraint(const std::string &parameterName) {
  auto cellAngleConstraint = std::make_unique<BoundaryConstraint>(this, parameterName, 0.0, 180.0, true);
  cellAngleConstraint->setPenaltyFactor(1e12);
  addConstraint(std::move(cellAngleConstraint));
}
 
void PawleyParameterFunction::setCenterParameterNameFromFunction(const IPeakFunction_sptr &profileFunction) {
  m_profileFunctionCenterParameterName.clear();
  if (profileFunction) {
    m_profileFunctionCenterParameterName = profileFunction->getCentreParameterName();
  }
}
 
DECLARE_FUNCTION(PawleyFunction)
 
 
PawleyFunction::PawleyFunction()
    : IPawleyFunction(), m_compositeFunction(), m_pawleyParameterFunction(), m_peakProfileComposite(), m_hkls(),
      m_dUnit(), m_wsUnit(), m_peakRadius(5) {
  auto peakRadius = Kernel::ConfigService::Instance().getValue<int>("curvefitting.peakRadius");
  m_peakRadius = peakRadius.get_value_or(5);
}
 
void PawleyFunction::setMatrixWorkspace(std::shared_ptr<const MatrixWorkspace> workspace, size_t wi, double startX,
                                        double endX) {
  if (workspace) {
    Axis *xAxis = workspace->getAxis(wi);
    Kernel::Unit_sptr wsUnit = xAxis->unit();
 
    if (std::dynamic_pointer_cast<Units::Empty>(wsUnit) || std::dynamic_pointer_cast<Units::dSpacing>(wsUnit)) {
      m_wsUnit = m_dUnit;
    } else {
      double factor, power;
      if (wsUnit->quickConversion(*m_dUnit, factor, power)) {
        m_wsUnit = wsUnit;
      } else {
        throw std::invalid_argument("Can not use quick conversion for unit.");
      }
    }
  }
 
  m_wrappedFunction->setMatrixWorkspace(workspace, wi, startX, endX);
}
 
void PawleyFunction::setLatticeSystem(const std::string &latticeSystem) {
  m_pawleyParameterFunction->setAttributeValue("LatticeSystem", latticeSystem);
  m_compositeFunction->checkFunction();
}
 
void PawleyFunction::setProfileFunction(const std::string &profileFunction) {
  m_pawleyParameterFunction->setAttributeValue("ProfileFunction", profileFunction);
 
  /* At this point PawleyParameterFunction guarantees that it's an IPeakFunction
   * and all existing profile functions are replaced.
   */
  for (size_t i = 0; i < m_peakProfileComposite->nFunctions(); ++i) {
    IPeakFunction_sptr oldFunction = std::dynamic_pointer_cast<IPeakFunction>(m_peakProfileComposite->getFunction(i));
 
    IPeakFunction_sptr newFunction = std::dynamic_pointer_cast<IPeakFunction>(
        FunctionFactory::Instance().createFunction(m_pawleyParameterFunction->getProfileFunctionName()));
 
    newFunction->setCentre(oldFunction->centre());
    try {
      newFunction->setFwhm(oldFunction->fwhm());
    } catch (...) {
      // do nothing.
    }
    newFunction->setHeight(oldFunction->height());
 
    m_peakProfileComposite->replaceFunction(i, newFunction);
  }
 
  // Update exposed parameters.
  m_compositeFunction->checkFunction();
}
 
void PawleyFunction::setUnitCell(const std::string &unitCellString) {
  m_pawleyParameterFunction->setParametersFromUnitCell(strToUnitCell(unitCellString));
}
 
double PawleyFunction::getTransformedCenter(double d) const {
  if ((m_dUnit && m_wsUnit) && m_dUnit != m_wsUnit) {
    return UnitConversion::run(*m_dUnit, *m_wsUnit, d, 0, DeltaEMode::Elastic);
  }
 
  return d;
}
 
void PawleyFunction::setPeakPositions(const std::string &centreName, double zeroShift, const UnitCell &cell) const {
  for (size_t i = 0; i < m_hkls.size(); ++i) {
    double centre = getTransformedCenter(cell.d(m_hkls[i]));
 
    m_peakProfileComposite->getFunction(i)->setParameter(centreName, centre + zeroShift);
  }
}
 
size_t PawleyFunction::calculateFunctionValues(const API::IPeakFunction_sptr &peak, const API::FunctionDomain1D &domain,
                                               API::FunctionValues &localValues) const {
  size_t domainSize = domain.size();
  const double *domainBegin = domain.getPointerAt(0);
  const double *domainEnd = domain.getPointerAt(domainSize);
 
  double centre = peak->centre();
  double dx = m_peakRadius * peak->fwhm();
 
  auto lb = std::lower_bound(domainBegin, domainEnd, centre - dx);
  auto ub = std::upper_bound(lb, domainEnd, centre + dx);
 
  size_t n = std::distance(lb, ub);
 
  if (n == 0) {
    throw std::invalid_argument("Null-domain");
  }
 
  FunctionDomain1DView localDomain(lb, n);
  localValues.reset(localDomain);
 
  peak->functionLocal(localValues.getPointerToCalculated(0), localDomain.getPointerAt(0), n);
 
  return std::distance(domainBegin, lb);
}
 
void PawleyFunction::function(const FunctionDomain &domain, FunctionValues &values) const {
  values.zeroCalculated();
  try {
    const auto &domain1D = dynamic_cast<const FunctionDomain1D &>(domain);
 
    UnitCell cell = m_pawleyParameterFunction->getUnitCellFromParameters();
    double zeroShift = m_pawleyParameterFunction->getParameter("ZeroShift");
    std::string centreName = m_pawleyParameterFunction->getProfileFunctionCenterParameterName();
 
    setPeakPositions(centreName, zeroShift, cell);
 
    FunctionValues localValues;
 
    for (size_t i = 0; i < m_peakProfileComposite->nFunctions(); ++i) {
      IPeakFunction_sptr peak = std::dynamic_pointer_cast<IPeakFunction>(m_peakProfileComposite->getFunction(i));
 
      try {
        size_t offset = calculateFunctionValues(peak, domain1D, localValues);
        values.addToCalculated(offset, localValues);
      } catch (const std::invalid_argument &) {
        // do nothing
      }
    }
 
    setPeakPositions(centreName, 0.0, cell);
  } catch (const std::bad_cast &) {
    // do nothing
  }
}
 
void PawleyFunction::clearPeaks() {
  m_peakProfileComposite =
      std::dynamic_pointer_cast<CompositeFunction>(FunctionFactory::Instance().createFunction("CompositeFunction"));
  m_compositeFunction->replaceFunction(1, m_peakProfileComposite);
  m_hkls.clear();
}
 
void PawleyFunction::setPeaks(const std::vector<Kernel::V3D> &hkls, double fwhm, double height) {
  clearPeaks();
 
  for (const auto &hkl : hkls) {
    addPeak(hkl, fwhm, height);
  }
}
 
void PawleyFunction::addPeak(const Kernel::V3D &hkl, double fwhm, double height) {
  m_hkls.emplace_back(hkl);
 
  IPeakFunction_sptr peak = std::dynamic_pointer_cast<IPeakFunction>(
      FunctionFactory::Instance().createFunction(m_pawleyParameterFunction->getProfileFunctionName()));
 
  peak->fix(peak->parameterIndex(m_pawleyParameterFunction->getProfileFunctionCenterParameterName()));
 
  try {
    peak->setFwhm(fwhm);
  } catch (...) {
    // do nothing.
  }
 
  peak->setHeight(height);
 
  m_peakProfileComposite->addFunction(peak);
 
  m_compositeFunction->checkFunction();
}
 
size_t PawleyFunction::getPeakCount() const { return m_hkls.size(); }
 
IPeakFunction_sptr PawleyFunction::getPeakFunction(size_t i) const {
  if (i >= m_hkls.size()) {
    throw std::out_of_range("Peak index out of range.");
  }
 
  return std::dynamic_pointer_cast<IPeakFunction>(m_peakProfileComposite->getFunction(i));
}
 
Kernel::V3D PawleyFunction::getPeakHKL(size_t i) const {
  if (i >= m_hkls.size()) {
    throw std::out_of_range("Peak index out of range.");
  }
 
  return m_hkls[i];
}
 
PawleyParameterFunction_sptr PawleyFunction::getPawleyParameterFunction() const { return m_pawleyParameterFunction; }
 
void PawleyFunction::init() {
  setDecoratedFunction("CompositeFunction");
 
  if (!m_compositeFunction) {
    throw std::runtime_error("PawleyFunction could not construct internal CompositeFunction.");
  }
 
  m_dUnit = UnitFactory::Instance().create("dSpacing");
}
 
void PawleyFunction::beforeDecoratedFunctionSet(const API::IFunction_sptr &fn) {
  CompositeFunction_sptr composite = std::dynamic_pointer_cast<CompositeFunction>(fn);
 
  if (!composite) {
    throw std::invalid_argument("PawleyFunction only works with "
                                "CompositeFunction. Selecting another "
                                "decorated function is not possible.");
  }
 
  m_compositeFunction = composite;
 
  if (m_compositeFunction->nFunctions() == 0) {
    m_peakProfileComposite =
        std::dynamic_pointer_cast<CompositeFunction>(FunctionFactory::Instance().createFunction("CompositeFunction"));
 
    m_pawleyParameterFunction = std::dynamic_pointer_cast<PawleyParameterFunction>(
        FunctionFactory::Instance().createFunction("PawleyParameterFunction"));
 
    m_compositeFunction->addFunction(m_pawleyParameterFunction);
    m_compositeFunction->addFunction(m_peakProfileComposite);
  } else {
    m_pawleyParameterFunction = std::dynamic_pointer_cast<PawleyParameterFunction>(m_compositeFunction->getFunction(0));
    m_peakProfileComposite = std::dynamic_pointer_cast<CompositeFunction>(m_compositeFunction->getFunction(1));
  }
}
 
} // namespace Mantid::CurveFitting::Functions