BSpline.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/Functions/BSpline.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidKernel/ArrayBoundedValidator.h"
#include "MantidKernel/BoundedValidator.h"
 
#include <algorithm>
 
namespace Mantid::CurveFitting::Functions {
 
using namespace CurveFitting;
using namespace Kernel;
using namespace API;
 
DECLARE_FUNCTION(BSpline)
 
 
BSpline::BSpline() {
  const size_t nbreak = 10;
 
  auto orderValidator = BoundedValidator<int>();
  orderValidator.setLower(1);
  declareAttribute("Order", Attribute(3), orderValidator);
 
  auto NBreakValidator = BoundedValidator<int>();
  NBreakValidator.setLower(2);
  declareAttribute("NBreak", Attribute(static_cast<int>(nbreak)), NBreakValidator);
 
  auto startXValidator = BoundedValidator<double>();
  startXValidator.setUpper(1.0);
  startXValidator.setUpperExclusive(true);
  declareAttribute("StartX", Attribute(0.0), startXValidator);
 
  auto endXValidator = BoundedValidator<double>();
  endXValidator.setLower(0.0);
  endXValidator.setLowerExclusive(true);
  declareAttribute("EndX", Attribute(1.0), endXValidator);
 
  auto breakPointsValidator = ArrayBoundedValidator<double>(0.0, 1.0);
  breakPointsValidator.setError(1e-8);
  declareAttribute("BreakPoints", Attribute(std::vector<double>(nbreak)), breakPointsValidator);
 
  declareAttribute("Uniform", Attribute(true));
 
  m_spline = Spline1D();
  resetParameters();
  resetKnots();
}
 
void BSpline::resetValidators() {
  auto attStartX = getAttribute("StartX");
  auto attEndX = getAttribute("EndX");
 
  auto startXValidator = std::dynamic_pointer_cast<BoundedValidator<double>>(attStartX.getValidator());
  startXValidator->setUpper(attEndX.asDouble());
 
  auto endXValidator = std::dynamic_pointer_cast<BoundedValidator<double>>(attEndX.getValidator());
  endXValidator->setLower(attStartX.asDouble());
 
  auto breakPointsValidator =
      std::dynamic_pointer_cast<ArrayBoundedValidator<double>>(getAttribute("BreakPoints").getValidator());
  breakPointsValidator->setLower(attStartX.asDouble());
  breakPointsValidator->setUpper(attEndX.asDouble());
}
 
void BSpline::function1D(double *out, const double *xValues, const size_t nData) const {
  size_t np = nParams();
  EigenVector B(np);
  size_t currentBBase = 0;
  double startX = getAttribute("StartX").asDouble();
  double endX = getAttribute("EndX").asDouble();
 
  for (size_t i = 0; i < nData; ++i) {
    double x = xValues[i];
    if (x < startX || x > endX) {
      out[i] = 0.0;
    } else {
      currentBBase = evaluateBasisFunctions(B, x, currentBBase);
 
      double val = 0.0;
      for (size_t j = 0; j < np; ++j) {
        val += getParameter(j) * B.get(j);
      }
      out[i] = val;
    }
  }
}
 
void BSpline::initialiseSpline(const std::vector<double> &breakPoints) {
  m_spline = Spline1D(EigenVector(m_knots).mutator(), EigenVector(breakPoints).mutator());
}
 
size_t BSpline::evaluateBasisFunctions(EigenVector &B, const double x, size_t currentBBase) const {
  int degree = getDegree();
  auto res = m_spline.BasisFunctions(x, degree, m_spline.knots());
 
  currentBBase = getSpanIndex(x, currentBBase);
  B.zero();
  for (int i = 0; i < res.size(); i++) { // Populate B
    B[currentBBase + i] = res.data()[i];
  }
  return currentBBase;
}
 
size_t BSpline::getSpanIndex(const double x, const size_t currentBBase, const bool clamped) const {
  const size_t clampedKnots = clamped ? static_cast<size_t>(getClampedKnots()) : 1u;
  const auto it = std::find_if(m_knots.cbegin() + currentBBase + clampedKnots, m_knots.cend(),
                               [x](double knot) { return x < knot; });
  if (it != m_knots.cend()) {
    return std::distance(m_knots.cbegin(), it) - clampedKnots;
  }
  return m_knots.size() - clampedKnots * 2;
}
 
void BSpline::derivative1D(double *out, const double *xValues, size_t nData, const size_t order) const {
  double startX = getAttribute("StartX").asDouble();
  double endX = getAttribute("EndX").asDouble();
  size_t jstart = 0;
  size_t splineOrder = static_cast<size_t>(getAttribute("Order").asInt());
 
  for (size_t i = 0; i < nData; ++i) {
    double x = xValues[i];
    if (x < startX || x > endX) {
      out[i] = 0.0;
    } else {
      jstart = getSpanIndex(x, jstart);
      size_t jend = jstart + splineOrder;
 
      auto res = evaluateBasisFnDerivatives(x, order);
      double val = 0.0;
      for (size_t j = jstart; j < jend; ++j) {
        val += getParameter(j) * res(order, j - jstart);
      }
      out[i] = val;
    }
  }
}
 
EigenMatrix BSpline::evaluateBasisFnDerivatives(const double x, const size_t derivOrder) const {
  int degree = getDegree();
  EigenMatrix res(derivOrder + 1u, degree + 1);
  res.mutator() = m_spline.BasisFunctionDerivatives(x, derivOrder, degree, m_spline.knots());
  return res;
}
 
void BSpline::setAttribute(const std::string &attName, const API::IFunction::Attribute &att) {
  bool isUniform = attName == "Uniform" && att.asBool();
  storeAttributeValue(attName, att);
 
  if (attName == "BreakPoints") {
    storeAttributeValue("NBreak", Attribute(static_cast<int>(att.asVector().size())));
    resetParameters();
    resetKnots();
  } else if (isUniform || attName == "StartX" || attName == "EndX") {
    resetValidators();
    resetKnots();
  } else if (attName == "NBreak" || attName == "Order") {
    resetParameters();
    resetKnots();
  }
}
 
std::vector<std::string> BSpline::getAttributeNames() const {
  return {"Uniform", "Order", "NBreak", "EndX", "StartX", "BreakPoints"};
}
 
void BSpline::resetParameters() {
  if (nParams() > 0) {
    clearAllParameters();
  }
  size_t np = getNBSplineCoefficients();
  for (size_t i = 0; i < np; ++i) {
    std::string pname = "A" + std::to_string(i);
    declareParameter(pname);
  }
}
 
void BSpline::resetKnots() {
  bool isUniform = getAttribute("Uniform").asBool();
  std::vector<double> breakPoints;
  if (isUniform) {
    // create uniform knots in the interval [StartX, EndX]
    double startX = getAttribute("StartX").asDouble();
    double endX = getAttribute("EndX").asDouble();
    // calc uniform break points
    breakPoints = calcUniformBreakPoints(startX, endX);
    storeAttributeValue("BreakPoints", Attribute(breakPoints));
    // calc uniform knots
    resetKnotVector(breakPoints);
  } else {
    // set the break points from BreakPoints vector attribute, update other attributes
    breakPoints = getAttribute("BreakPoints").asVector();
    // check that points are in ascending order
    double prev = breakPoints[0];
    for (size_t i = 1; i < breakPoints.size(); ++i) {
      double next = breakPoints[i];
      if (next < prev) {
        throw std::invalid_argument("BreakPoints must be in ascending order.");
      }
      prev = next;
    }
    int nbreaks = getNBreakPoints();
    // if number of break points change do necessary updates
    if (static_cast<size_t>(nbreaks) != breakPoints.size()) {
      storeAttributeValue("NBreak", Attribute(static_cast<int>(breakPoints.size())));
      resetParameters();
    }
    resetKnotVector(breakPoints);
  }
  initialiseSpline(breakPoints);
}
 
std::vector<double> BSpline::calcUniformBreakPoints(const double startX, const double endX) {
  const int nBreak = getNBreakPoints();
  std::vector<double> breakPoints(nBreak);
  const double interval = (endX - startX) / (nBreak - 1.0);
  std::generate(breakPoints.begin(), breakPoints.end(),
                [n = 0, &interval, &startX]() mutable { return n++ * interval + startX; });
  return breakPoints;
}
 
void BSpline::resetKnotVector(const std::vector<double> &breakPoints) {
  const int nKnots = getNKnots();
  const int clampedKnots = getClampedKnots();
  m_knots.clear();
  m_knots.resize(nKnots);
 
  for (int i = 0; i < nKnots; i++) {
    if (i < clampedKnots) {
      m_knots[i] = breakPoints[0];
    } else if (i >= nKnots - clampedKnots) {
      m_knots[i] = breakPoints[breakPoints.size() - 1];
    } else {
      m_knots[i] = breakPoints[i - clampedKnots + 1];
    }
  }
}
 
} // namespace Mantid::CurveFitting::Functions