No OneTemporary
Actions

Size

28 KB

Subscribers

None

View Options

	diff --git a/inc/Lau1DCubicSpline.hh b/inc/Lau1DCubicSpline.hh
	index 83f1bd1..11b89c2 100644
	--- a/inc/Lau1DCubicSpline.hh
	+++ b/inc/Lau1DCubicSpline.hh
	@@ -1,261 +1,276 @@

	/*
	Copyright 2015 University of Warwick

	Licensed under the Apache License, Version 2.0 (the "License");
	you may not use this file except in compliance with the License.
	You may obtain a copy of the License at

	http://www.apache.org/licenses/LICENSE-2.0

	Unless required by applicable law or agreed to in writing, software
	distributed under the License is distributed on an "AS IS" BASIS,
	WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
	See the License for the specific language governing permissions and
	limitations under the License.
	*/

	/*
	Laura++ package authors:
	John Back
	Paul Harrison
	Thomas Latham
	*/

	/*! \file Lau1DCubicSpline.hh
	\brief File containing declaration of Lau1DCubicSpline class.
	*/

	/*! \class Lau1DCubicSpline
	\brief Class for defining a 1D cubic spline based on a set of knots.

	Class for defining a 1D cubic spline based on a set of knots.
	Interpolation between the knots is performed either by one of
	two types of cubic spline (standard or Akima) or by linear interpolation.
	The splines are defined by a piecewise cubic function which, between knots i and i+1, has the form

	f_i(x) = (1 - t)y_i + ty_i+1 + t(1 - t)(a(1 - t) + b*t)

	where t = (x - x_i)/(x_i+1 - x_i),
	a = k_i *(x_i+1 - x_i) - (y_i+1 - y_i),
	b = -k_i+1*(x_i+1 - x_i) + (y_i+1 - y_i)

	and k_i is (by construction) the first derivative at knot i.
	f(x) and f'(x) are continuous at the internal knots by construction.

	For the standard splines, f''(x) is required to be continuous
	at all internal knots placing n-2 constraints on the n parameters, k_i.
	The final two constraints are set by the boundary conditions.
	At each boundary, the function may be:

	(i) Clamped : f'(x) = C at the last knot
	(ii) Natural : f''(x) = 0 at the last knot
	(iii) Not a knot : f'''(x) continuous at the second last knot

	The algorithms used in these splines can be found on:
	http://en.wikipedia.org/wiki/Spline_interpolation#Algorithm_to_find_the_interpolating_cubic_spline
	http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

	For the Akima splines, the values of k_i are determined from the slopes of the four nearest segments (a_i-1, a_i, a_i+1 and a_i+2) as

	k_i = ( \| a_i+2 - a_i+1 \| a_i + \| a_i - a_i-1 \| a_i+1 ) / ( \| a_i+2 - a_i+1 \| + \| a_i - a_i-1 \| )

	and as

	k_i = ( a_i + a_i+1 ) / 2

	in the special case a_i-1 == a_i != a_i+1 == a_i+2.

	Boundary conditions are specified by the relations
	a_2 - a_1 = a_1 - a_0 = a_0 - a_-1 and
	a_n-1 - a_n-2 = a_n - a_n-1 = a_n+1 - a_n

	The algorithms used in these splines can be found in:
	J.ACM vol. 17 no. 4 pp 589-602
	*/

	#ifndef LAU_1DCUBICSPLINE
	#define LAU_1DCUBICSPLINE

	#include <vector>
	#include <array>

	#include "Rtypes.h"
	#include "TF1.h"
	#include "TFitResultPtr.h"

	#include "LauPrint.hh"

	class TH1;
	class LauAbsRValue;
	class LauParameter;

	class Lau1DCubicSpline final {

	public:
	//! Define the allowed interpolation types
	enum class SplineType {
	StandardSpline, /!< standard cubic splines with f''(x) continuous at all internal knots /
	AkimaSpline, /!< Akima cubic splines with f'(x) at each knot defined locally by the positions of only five knots /
	LinearInterpolation /! Linear interpolation between each pair of knots /
	};

	//! Define the allowed boundary condition types
	/*!
	These are only supported by standard splines
	*/
	enum class BoundaryType {
	Clamped, /!< clamped boundary - f'(x) = C /
	Natural, /!< natural boundary - f''(x) = 0 /
	NotAKnot /!< 'not a knot' boundary - f'''(x) continuous at second last knot /
	};

	//! Constructor
	/*!
	\param [in] xs the x-values of the knots
	\param [in] ys the y-values of the knots
	\param [in] type the type of spline (Standard, Akima, Linear)
	\param [in] leftBound the left-hand boundary condition
	\param [in] rightBound the right-hand boundary condition
	\param [in] dydx0 the gradient at the left-hand boundary of a clamped spline
	\param [in] dydxn the gradient at the right-hand boundary of a clamped spline
	*/
	Lau1DCubicSpline(const std::vector<Double_t>& xs, const std::vector<Double_t>& ys,
	const SplineType type = SplineType::StandardSpline,
	const BoundaryType leftBound = BoundaryType::NotAKnot,
	const BoundaryType rightBound = BoundaryType::NotAKnot,
	const Double_t dydx0 = 0.0, const Double_t dydxn = 0.0);

	//! Evaluate the function at given point
	/*!
	\param [in] x the x-coordinate
	\return the value of the spline at x
	*/
	Double_t evaluate(const Double_t x) const;

	//! Update the y-values of the knots
	/*!
	\param [in] ys the y-values of the knots
	*/
	void updateYValues(const std::vector<Double_t>& ys);
	//! Update the y-values of the knots
	/*!
	\param [in] ys the y-values of the knots
	*/
	void updateYValues(const std::vector<LauParameter*>& ys);
	//! Update the y-values of the knots
	/*!
	\param [in] ys the y-values of the knots
	*/
	void updateYValues(const std::vector<LauAbsRValue*>& ys);

	//! Update the type of interpolation to perform
	/*!
	\param [in] type the type of interpolation
	*/
	void updateType(const SplineType type);

	//! Update the boundary conditions for the spline
	/*!
	\param [in] leftBound the left-hand boundary condition
	\param [in] rightBound the right-hand boundary condition
	\param [in] dydx0 the gradient at the left-hand boundary of a clamped spline
	\param [in] dydxn the gradient at the right-hand boundary of a clamped spline
	*/
	void updateBoundaryConditions(const BoundaryType leftBound,
	const BoundaryType rightBound,
	const Double_t dydx0 = 0.0,
	const Double_t dydxn = 0.0);

	//! Get the number of knots
	std::size_t getnKnots() const {return nKnots_;}

	//! Get y values
	const std::vector<Double_t>& getYValues() const {return y_;}

	//! Get x values
	const std::vector<Double_t>& getXValues() const {return x_;}

	//! Get the coefficients of a given spline segment in the form a + bx + cx^2 + dx^3
	/*!
	\param [in] segIndex refers to the index of the knot at the left end of the segment (segIndex = 0 gets the coefficients of the the segment between x_0 and x_1)
	\param [in] normalise if true, the coefficients returned will be normalised by the integral of the complete spline (defaults to false)
	\return coefficients {a, b, c, d}
	*/
	std::array<Double_t,4> getCoefficients(const std::size_t segIndex, const bool normalise = false) const;

	//! Get the integral over all the spline segments
	Double_t integral() const;

	//! Make a TF1 showing the spline with its current knot values
	/*!
	\param [in] normalise whether or not you want the spline normalised
	\return 1D function object
	*/
	TF1* makeTF1(const bool normalise = false) const;

	//! Fit the the normalisation of the spline to a TH1
	/*!
	\param [in] hist the histogram to be fitted
	\param [in] printLevel the level of printout desired from fit
	\return a TF1 fit to the histogram
	*/
	TF1* normaliseToTH1(TH1* hist, const LauOutputLevel printLevel = LauOutputLevel::Standard) const;

	//! Fit the spline to a TH1
	/*!
	\param [in] hist the histogram to be fitted
	\param [in] printLevel the level of printout desired from fit
	\param [in] doWL If true do a weighted likelihood fit, else chi2
	\return a shared-ownership smart pointer to the fit results
	*/
	TFitResultPtr fitToTH1(TH1* hist, const LauOutputLevel printLevel = LauOutputLevel::Standard, const bool doWL = false);

	+ //! Set spline variables to those read from a json file
	+ /*
	+ \param[in] filename to be used
	+ \return the newly constructed spline
	+ */
	+ //Lau1DCubicSpline readFromJson(const TString filename);
	+
	+ //! Write the spline to a json file using nlohmann json
	+ /*!
	+ \param[in] filename to be used
	+ \param[in] name of spline, if empty, no name is used
	+ \param[in] if true and filename already exists, append the spline to the file
	+ */
	+ void writeToJson(const TString filename, const TString splineName = "", const bool append = false) const;
	+
	private:
	//! Initialise the class
	void init();

	//! Calculate the first derivative at each knot
	void calcDerivatives();

	//! Calculate the first derivatives according to the standard method
	void calcDerivativesStandard();
	//! Calculate the first derivatives according to the Akima method
	void calcDerivativesAkima();

	//! The number of knots in the spline
	const std::size_t nKnots_;

	//! The x-value at each knot
	const std::vector<Double_t> x_;
	//! The y-value at each knot
	std::vector<Double_t> y_;
	//! The first derivative at each knot
	std::vector<Double_t> dydx_;

	//! The 'a' coefficients used to determine the derivatives
	std::vector<Double_t> a_;
	//! The 'b' coefficients used to determine the derivatives
	std::vector<Double_t> b_;
	//! The 'c' coefficients used to determine the derivatives
	std::vector<Double_t> c_;
	//! The 'd' coefficients used to determine the derivatives
	std::vector<Double_t> d_;

	//! The type of interpolation to be performed
	SplineType type_;

	//! The left-hand boundary condition on the spline
	BoundaryType leftBound_;
	//! The right-hand boundary condition on the spline
	BoundaryType rightBound_;

	//! The gradient at the left boundary for a clamped spline
	Double_t dydx0_;
	//! The gradient at the right boundary for a clamped spline
	Double_t dydxn_;

	ClassDef(Lau1DCubicSpline, 0) // Class for defining a 1D cubic spline
	};

	#endif
	diff --git a/src/Lau1DCubicSpline.cc b/src/Lau1DCubicSpline.cc
	index 15cd255..38cd91c 100644
	--- a/src/Lau1DCubicSpline.cc
	+++ b/src/Lau1DCubicSpline.cc
	@@ -1,590 +1,620 @@

	/*
	Copyright 2015 University of Warwick

	Licensed under the Apache License, Version 2.0 (the "License");
	you may not use this file except in compliance with the License.
	You may obtain a copy of the License at

	http://www.apache.org/licenses/LICENSE-2.0

	Unless required by applicable law or agreed to in writing, software
	distributed under the License is distributed on an "AS IS" BASIS,
	WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
	See the License for the specific language governing permissions and
	limitations under the License.
	*/

	/*
	Laura++ package authors:
	John Back
	Paul Harrison
	Thomas Latham
	*/

	/*! \file Lau1DCubicSpline.cc
	\brief File containing implementation of Lau1DCubicSpline class.
	*/

	#include <cmath>
	#include <cstdlib>
	#include <iostream>
	+#include <fstream>

	#include <TH2.h>
	#include <TSystem.h>

	#include "Lau1DCubicSpline.hh"
	#include "LauAbsRValue.hh"
	#include "LauParameter.hh"

	+#include <nlohmann/json.hpp>
	+
	ClassImp(Lau1DCubicSpline)


	Lau1DCubicSpline::Lau1DCubicSpline(const std::vector<Double_t>& xs, const std::vector<Double_t>& ys,
	const SplineType type,
	const BoundaryType leftBound,
	const BoundaryType rightBound,
	const Double_t dydx0, const Double_t dydxn) :
	nKnots_{xs.size()},
	x_{xs},
	y_{ys},
	type_{type},
	leftBound_{leftBound},
	rightBound_{rightBound},
	dydx0_{dydx0},
	dydxn_{dydxn}
	{
	init();
	}
	+/*
	+Lau1DCubicSpline Lau1DCubicSpline::readFromJson(const TString filename)
	+{
	+ //return;
	+}
	+*/
	+void Lau1DCubicSpline::writeToJson(const TString filename, const TString splineName, const bool append) const
	+{
	+ using namespace nlohmann;
	+ json j = {
	+ {"nKnots",nKnots_},
	+ {"x", x_},
	+ {"y", y_},
	+ {"type", type_},
	+ {"leftBound", leftBound_},
	+ {"rightBound", rightBound_},
	+ {"dydx0", dydx0_},
	+ {"dydxn", dydxn_}
	+ };
	+ if(splineName != "")
	+ { j = json::array({ splineName, j }); }
	+ std::fstream out;
	+ out.open(filename, append ? std::ios_base::app : std::ios_base::out);
	+ out << j.dump(4);
	+ out.close();
	+ return;
	+}

	Double_t Lau1DCubicSpline::evaluate(const Double_t x) const
	{
	// do not attempt to extrapolate the spline
	if( x<x_[0] \|\| x>x_[nKnots_-1] ) {
	std::cerr << "WARNING in Lau1DCubicSpline::evaluate : function is only defined between " << x_[0] << " and " << x_[nKnots_-1] << std::endl;
	std::cerr << " value at " << x << " returned as 0" << std::endl;
	return 0.;
	}

	// first determine which 'cell' of the spline x is in
	// cell i runs from knot i to knot i+1
	std::size_t cell{0};

	while( x > x_[cell+1] ) {
	++cell;
	}

	// obtain x- and y-values of the neighbouring knots
	const Double_t xLow { x_[cell] };
	const Double_t xHigh { x_[cell+1] };
	const Double_t yLow { y_[cell] };
	const Double_t yHigh { y_[cell+1] };

	if(type_ == SplineType::LinearInterpolation) {
	return yHigh(x-xLow)/(xHigh-xLow) + yLow(xHigh-x)/(xHigh-xLow);
	}

	// obtain t, the normalised x-coordinate within the cell,
	// and the coefficients a and b, which are defined in cell i as:
	//
	// a_i = k_i *(x_i+1 - x_i) - (y_i+1 - y_i),
	// b_i = -k_i+1*(x_i+1 - x_i) + (y_i+1 - y_i)
	//
	// where k_i is (by construction) the first derivative at knot i
	const Double_t t { (x - xLow) / (xHigh - xLow) };
	const Double_t a { dydx_[cell] * (xHigh - xLow) - (yHigh - yLow) };
	const Double_t b { -1.dydx_[cell+1] (xHigh - xLow) + (yHigh - yLow) };

	const Double_t retVal { (1 - t) * yLow + t * yHigh + t * (1 - t) * ( a * (1 - t) + b * t ) };

	return retVal;
	}

	void Lau1DCubicSpline::updateYValues(const std::vector<Double_t>& ys)
	{
	y_ = ys;
	this->calcDerivatives();
	}

	void Lau1DCubicSpline::updateYValues(const std::vector<LauParameter*>& ys)
	{
	for (std::size_t i{0}; i<nKnots_; ++i){
	y_[i] = ys[i]->unblindValue();
	}
	this->calcDerivatives();
	}

	void Lau1DCubicSpline::updateYValues(const std::vector<LauAbsRValue*>& ys)
	{
	for (std::size_t i{0}; i<nKnots_; ++i){
	y_[i] = ys[i]->unblindValue();
	}
	this->calcDerivatives();
	}

	void Lau1DCubicSpline::updateType(const SplineType type)
	{
	if(type_ != type) {
	type_ = type;
	this->calcDerivatives();
	}
	}

	void Lau1DCubicSpline::updateBoundaryConditions(const BoundaryType leftBound,
	const BoundaryType rightBound,
	const Double_t dydx0, const Double_t dydxn)
	{
	bool updateDerivatives{false};

	if(leftBound_ != leftBound \|\| rightBound_ != rightBound) {
	leftBound_ = leftBound;
	rightBound_ = rightBound;
	updateDerivatives = true;
	}

	if(dydx0_ != dydx0) {
	dydx0_ = dydx0;
	if(leftBound_ == BoundaryType::Clamped) updateDerivatives = true;
	}

	if(dydxn_ != dydxn) {
	dydxn_ = dydxn;
	if(rightBound_ == BoundaryType::Clamped) updateDerivatives = true;
	}

	if(updateDerivatives) {
	this->calcDerivatives();
	}
	}

	std::array<Double_t,4> Lau1DCubicSpline::getCoefficients(const std::size_t segIndex, const bool normalise) const
	{
	std::array<Double_t,4> result = {0.,0.,0.,0.};

	if(segIndex >= nKnots_-1)
	{
	std::cerr << "ERROR in Lau1DCubicSpline::getCoefficients requested for too high a knot value" << std::endl;
	return result;
	}

	Double_t xL = x_[segIndex] , xH = x_[segIndex+1];
	Double_t yL = y_[segIndex] , yH = y_[segIndex+1];
	Double_t h = xH-xL; //This number comes up a lot

	switch(type_)
	{
	case SplineType::StandardSpline:
	case SplineType::AkimaSpline:
	{
	Double_t kL = dydx_[segIndex], kH = dydx_[segIndex+1];

	//a and b based on definitions from https://en.wikipedia.org/wiki/Spline_interpolation#Algorithm_to_find_the_interpolating_cubic_spline
	Double_t a = kL*h-(yH-yL);
	Double_t b =-kH*h+(yH-yL);

	Double_t denom = -hhh;//The terms have a common demoninator

	result[0] = -bxLxLxH + axLxHxH + hh(xLyH - xHyL);
	result[1] = -axH(2xL+xH) + bxL(xL+2xH) + hh(yL-yH);
	result[2] = -b(2xL+xH) + a(xL+2xH);
	result[3] = -a+b;

	for(auto& res : result){res /= denom;}
	break;
	}
	/* case SplineType::AkimaSpline: // Double check the Akima description of splines (in evaluate) right now they're the same except for the first derivatives
	{
	//using fomulae from https://asmquantmacro.com/2015/09/01/akima-spline-interpolation-in-excel/
	std::function<Double_t(Int_t)> m = [&](Int_t j) //formula to get the straight line gradient
	{
	if(j < 0){return 2*m(j+1)-m(j+2);}
	if(j >= nKnots_){return 2*m(j-1)-m(j-2);}
	return (y_[j+1]-y_[j]) / (x_[j+1]-x_[j]);
	};

	auto t = [&](Int_t j)
	{
	Double_t res = 0.; //originally res was called 't' but that produced a shadow warning
	Double_t denom = TMath::Abs(m(j+1)-m(j)) + TMath::Abs(m(j-1)-m(j-2));
	if(denom == 0){res = (m(j)-m(j-1))/2;} //Special case for when denom = 0
	else
	{
	res = TMath::Abs(m(j+1)-m(j))m(j-1) + TMath::Abs(m(j-1)-m(j-2))m(j);
	res /= denom;
	}
	return res;
	};

	//These are the p's to get the spline in the form p_k(x-xL)^k
	Double_t pDenom = x_[segIndex+1]/x_[segIndex]; //a denominator used for p[2] and p[3]
	std::array<Double_t,4> p = {y_[segIndex],t(segIndex),0.,0.}; //we'll do the last 2 below

	p[2] = 3m(segIndex)-2t(segIndex)-t(segIndex+1);
	p[2]/= pDenom;

	p[3] = t(segIndex)+t(segIndex+1)-2*m(segIndex);
	p[3]/= pDenom*pDenom;

	//Now finally rearranging the p's into the desired results
	result[0] = p[0]-p[1]xL+p[2]xLxL-p[3]xLxLxL;
	result[1] = p[1]-2p[2]xL+3p[3]xL*xL;
	result[2] = p[2]-3p[3]xL;
	result[3] = p[3];

	break;
	}*/
	case SplineType::LinearInterpolation:
	{
	result[0] = xHyL-xLyH;
	result[1] = yH-yL;

	for(auto& res : result){res /= h;}
	break;
	}
	}

	if (normalise) {
	const Double_t integral { this->integral() };
	for(auto& res : result){res /= integral;}
	}

	return result;
	}

	Double_t Lau1DCubicSpline::integral() const
	{
	Double_t integral{0.0};
	for(std::size_t iKnot{0}; iKnot < nKnots_ -1; ++iKnot)
	{
	const Double_t minAbs { x_[iKnot] };
	const Double_t maxAbs { x_[iKnot+1] };

	const std::array<Double_t,4> coeffs { this->getCoefficients(iKnot, false) };

	auto integralFunc = [&coeffs](Double_t x){return coeffs[0]x + coeffs[1]xx/2.0 + coeffs[2]xxx/3.0 + coeffs[3]xxxx/4.0;};

	integral += integralFunc(maxAbs);
	integral -= integralFunc(minAbs);
	}
	return integral;
	}

	TF1* Lau1DCubicSpline::makeTF1(const bool normalise) const
	{
	auto fitf = [this](Double_t* x, Double_t* par) {//there is only 1 x (the abscissa) and 1 par (a scaling of the entire thing for normalisation)
	return this->evaluate( x[0] ) * par[0];
	};

	//Make the function
	TF1* func = new TF1("FittedFunc",fitf,x_.front(),x_.back(),1);

	func -> SetParNames("Normalisation");
	if(normalise)
	{
	func->SetParameter(0, 1./func->Integral(x_.front(),x_.back()));
	}
	else{func->SetParameter(0,1.);}

	return func;
	}

	TF1* Lau1DCubicSpline::normaliseToTH1(TH1* hist, const LauOutputLevel printLevel) const
	{
	//first define the fit function
	auto fitf = [this](Double_t* x, Double_t* par) {//there is only 1 x (the abscissa) and 1 par (a scaling of the entire thing)
	return this->evaluate( x[0] ) * par[0];
	};
	//Make the function
	TF1* FittedFunc = new TF1("FittedFunc",fitf,x_.front(),x_.back(),1);
	//Set the parameter name
	FittedFunc -> SetParNames("Constant");
	//Set the parameter limits and default value
	FittedFunc -> SetParLimits(0,0.,10.);
	FittedFunc -> SetParameter(0,1.);

	// Options to:
	// - N = do not store the graphics function
	// - 0 = do not plot the fit result
	TString fitOptions {"N 0"};

	switch ( printLevel ) {
	case LauOutputLevel::Verbose :
	fitOptions += " V";
	break;
	case LauOutputLevel::Quiet :
	fitOptions += " Q";
	break;
	case LauOutputLevel::Standard :
	break;
	}

	hist->Fit( "FittedFunc", fitOptions );

	return FittedFunc;
	}

	TFitResultPtr Lau1DCubicSpline::fitToTH1(TH1* hist, const LauOutputLevel printLevel, const bool doWL)
	{
	auto fitf = [this](Double_t* x, Double_t* par) {
	this -> updateYValues( std::vector<Double_t>(par, par + nKnots_) );
	return this -> evaluate( x[0] );
	};

	TF1 FittedFunc("FittedFunc",fitf,x_.front(),x_.back(),nKnots_);

	const Double_t knotMax { hist->GetMaximum() * 1.5 };

	for(std::size_t knot{0}; knot <= nKnots_ ; ++knot)
	{
	FittedFunc.SetParName(knot, Form("Knot%lu",knot));
	FittedFunc.SetParLimits(knot, 0., knotMax);
	FittedFunc.SetParameter(knot, y_[knot]);
	}

	// Options to:
	// - N = do not store the graphics function
	// - 0 = do not plot the fit result
	// - S = save the result of the fit
	// - WL= do a weighted likelihood fit (not chi2)
	TString fitOptions {"N 0 S"};
	if(doWL){fitOptions += " WL";}

	switch ( printLevel ) {
	case LauOutputLevel::Verbose :
	fitOptions += " V";
	break;
	case LauOutputLevel::Quiet :
	fitOptions += " Q";
	break;
	case LauOutputLevel::Standard :
	break;
	}

	return hist->Fit( "FittedFunc", fitOptions );
	}

	void Lau1DCubicSpline::init()
	{
	if( y_.size() != x_.size()) {
	std::cerr << "ERROR in Lau1DCubicSpline::init : The number of y-values given does not match the number of x-values" << std::endl;
	std::cerr << " Found " << y_.size() << ", expected " << x_.size() << std::endl;
	gSystem->Exit(EXIT_FAILURE);
	}

	if( nKnots_ < 3 ) {
	std::cerr << "ERROR in Lau1DCubicSpline::init : The number of knots is too small" << std::endl;
	std::cerr << " Found " << nKnots_ << ", expected at least 3 (to have at least 1 internal knot)" << std::endl;
	gSystem->Exit(EXIT_FAILURE);
	}

	dydx_.assign(nKnots_,0.0);

	a_.assign(nKnots_,0.0);
	b_.assign(nKnots_,0.0);
	c_.assign(nKnots_,0.0);
	d_.assign(nKnots_,0.0);

	this->calcDerivatives();
	}

	void Lau1DCubicSpline::calcDerivatives()
	{
	switch ( type_ ) {
	case SplineType::StandardSpline :
	this->calcDerivativesStandard();
	break;
	case SplineType::AkimaSpline :
	this->calcDerivativesAkima();
	break;
	case SplineType::LinearInterpolation :
	//derivatives not needed for linear interpolation
	break;
	}
	}

	void Lau1DCubicSpline::calcDerivativesStandard()
	{
	// derivatives are determined such that the second derivative is continuous at internal knots

	// derivatives, k_i, are the solutions to a set of linear equations of the form:
	// a_i * k_i-1 + b_i * k+i + c_i * k_i+1 = d_i with a_0 = 0, c_n-1 = 0

	// this is solved using the tridiagonal matrix algorithm as on en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

	// first and last equations give boundary conditions
	// - for natural boundary, require f''(x) = 0 at end knot
	// - for 'not a knot' boundary, require f'''(x) continuous at second knot
	// - for clamped boundary, require predefined value of f'(x) at end knot

	// non-zero values of a_0 and c_n-1 would give cyclic boundary conditions
	a_[0] = 0.;
	c_[nKnots_-1] = 0.;

	// set left boundary condition
	switch ( leftBound_ ) {
	case BoundaryType::Natural :
	{
	b_[0] = 2./(x_[1]-x_[0]);
	c_[0] = 1./(x_[1]-x_[0]);
	d_[0] = 3.(y_[1]-y_[0])/((x_[1]-x_[0])(x_[1]-x_[0]));
	break;
	}

	case BoundaryType::NotAKnot :
	{
	// define the width, h, and the 'slope', delta, of the first cell
	const Double_t h1{x_[1]-x_[0]};
	const Double_t h2{x_[2]-x_[0]};
	const Double_t delta1{(y_[1]-y_[0])/h1};
	const Double_t delta2{(y_[2]-y_[1])/h2};
	// these coefficients can be determined by requiring f'''_0(x_1) = f'''_1(x_1)
	// the requirement f''_0(x_1) = f''_1(x_1) has been used to remove the dependence on k_2
	b_[0] = h2;
	c_[0] = h1+h2;
	d_[0] = delta1(2.h2h2 + 3.h1h2)/(h1+h2) + delta25.h1h1/(h1+h2);
	break;
	}

	case BoundaryType::Clamped :
	{
	b_[0] = 1.;
	c_[0] = 0.;
	d_[0] = dydx0_;
	break;
	}
	}

	// set right boundary condition
	switch ( rightBound_ ) {
	case BoundaryType::Natural :
	{
	a_[nKnots_-1] = 1./(x_[nKnots_-1]-x_[nKnots_-2]);
	b_[nKnots_-1] = 2./(x_[nKnots_-1]-x_[nKnots_-2]);
	d_[nKnots_-1] = 3.(y_[nKnots_-1]-y_[nKnots_-2])/((x_[nKnots_-1]-x_[nKnots_-2])(x_[nKnots_-1]-x_[nKnots_-2]));
	break;
	}
	case BoundaryType::NotAKnot :
	{
	// define the width, h, and the 'slope', delta, of the last cell
	const Double_t hnm1{x_[nKnots_-1]-x_[nKnots_-2]};
	const Double_t hnm2(x_[nKnots_-2]-x_[nKnots_-3]);
	const Double_t deltanm1{(y_[nKnots_-1]-y_[nKnots_-2])/hnm1};
	const Double_t deltanm2{(y_[nKnots_-2]-y_[nKnots_-3])/hnm2};
	// these coefficients can be determined by requiring f'''_n-3(x_n-2) = f'''_n-2(x_n-2)
	// the requirement f''_n-3(x_n-2) = f''_n-2(x_n-2) has been used to remove
	// the dependence on k_n-3
	a_[nKnots_-1] = hnm2 + hnm1;
	b_[nKnots_-1] = hnm1;
	d_[nKnots_-1] = deltanm2hnm1hnm1/(hnm2+hnm1) + deltanm1(2.hnm2hnm2 + 3.hnm2*hnm1)/(hnm2+hnm1);
	break;
	}

	case BoundaryType::Clamped :
	{
	a_[nKnots_-1] = 0.;
	b_[nKnots_-1] = 1.;
	d_[nKnots_-1] = dydxn_;
	break;
	}
	}

	// the remaining equations ensure that f_i-1''(x_i) = f''_i(x_i) for all internal knots
	for(std::size_t i{1}; i<nKnots_-1; ++i) {
	a_[i] = 1./(x_[i]-x_[i-1]);
	b_[i] = 2./(x_[i]-x_[i-1]) + 2./(x_[i+1]-x_[i]);
	c_[i] = 1./(x_[i+1]-x_[i]);
	d_[i] = 3.(y_[i]-y_[i-1])/((x_[i]-x_[i-1])(x_[i]-x_[i-1])) + 3.(y_[i+1]-y_[i])/((x_[i+1]-x_[i])(x_[i+1]-x_[i]));
	}

	// forward sweep - replace c_i and d_i with
	//
	// c'_i = c_i / b_i for i = 0
	// c_i / (b_i - a_i * c_i-1) otherwise
	//
	// and
	//
	// d'_i = d_i / b_i for i = 0
	// (d_i - a_i * d_i-1) / (b_i - a_i * c_i-1) otherwise

	c_[0] /= b_[0];
	d_[0] /= b_[0];

	for(std::size_t i{1}; i<nKnots_-1; ++i) {
	c_[i] = c_[i] / (b_[i] - a_[i]*c_[i-1]);
	d_[i] = (d_[i] - a_[i]d_[i-1]) / (b_[i] - a_[i]c_[i-1]);
	}

	d_[nKnots_-1] = (d_[nKnots_-1] - a_[nKnots_-1]d_[nKnots_-2]) / (b_[nKnots_-1] - a_[nKnots_-1]c_[nKnots_-2]);

	// back substitution - calculate k_i = dy/dx at each knot
	//
	// k_i = d'_i for i = n-1
	// d'_i - c'_i * k_i+1 otherwise

	dydx_[nKnots_-1] = d_[nKnots_-1];

	for(int i {static_cast<int>(nKnots_-2)}; i>=0; --i) {
	dydx_[i] = d_[i] - c_[i]*dydx_[i+1];
	}
	}

	void Lau1DCubicSpline::calcDerivativesAkima() {
	//derivatives are calculated according to the Akima method
	// J.ACM vol. 17 no. 4 pp 589-602

	Double_t am1{0.0}, an{0.0}, anp1{0.0};

	// a[i] is the slope of the segment from point i-1 to point i
	//
	// n.b. segment 0 is before point 0 and segment n is after point n-1
	// internal segments are numbered 1 - n-1
	for(std::size_t i{1}; i<nKnots_; ++i) {
	a_[i] = (y_[i]-y_[i-1])/(x_[i]-x_[i-1]);
	}

	// calculate slopes between additional points on each end according to the Akima method
	// method assumes that the additional points follow a quadratic defined by the last three points
	// this leads to the relations a[2] - a[1] = a[1] - a[0] = a[0] - a[-1] and a[n-1] - a[n-2] = a[n] - a[n-1] = a[n+1] - a[n]
	a_[0] = 2*a_[1] - a_[2];
	am1 = 2*a_[0] - a_[1];

	an = 2*a_[nKnots_-1] - a_[nKnots_-2];
	anp1 = 2*an - a_[nKnots_-1];

	// b[i] is the weight of a[i] towards dydx[i]
	// c[i] is the weight of a[i+1] towards dydx[i]
	// See Appendix A of J.ACM vol. 17 no. 4 pp 589-602 for a justification of these weights
	b_[0] = TMath::Abs(a_[1] - a_[0]);
	b_[1] = TMath::Abs(a_[2] - a_[1]);

	c_[0] = TMath::Abs(a_[0] - am1 );
	c_[1] = TMath::Abs(a_[1] - a_[0]);

	for(std::size_t i{2}; i<nKnots_-2; ++i) {
	b_[i] = TMath::Abs(a_[i+2] - a_[i+1]);
	c_[i] = TMath::Abs(a_[i] - a_[i-1]);
	}

	b_[nKnots_-2] = TMath::Abs(an - a_[nKnots_-1]);
	b_[nKnots_-1] = TMath::Abs(anp1 - an );

	c_[nKnots_-2] = TMath::Abs(a_[nKnots_-2] - a_[nKnots_-3]);
	c_[nKnots_-1] = TMath::Abs(a_[nKnots_-1] - a_[nKnots_-2]);

	// dy/dx calculated as the weighted average of a[i] and a[i+1]:
	// dy/dx_i = ( \| a_i+2 - a_i+1 \| a_i + \| a_i - a_i-1 \| a_i+1 ) / ( \| a_i+2 - a_i+1 \| + \| a_i - a_i-1 \| )
	// in the special case a_i-1 == a_i != a_i+1 == a_i+2 this function is undefined so dy/dx is then defined as (a_i + a_i+1) / 2
	for(std::size_t i{0}; i<nKnots_-2; ++i) {
	if(b_[i]==0 && c_[i]==0) {
	dydx_[i] = ( a_[i] + a_[i+1] ) / 2.;
	} else {
	dydx_[i] = ( b_[i] * a_[i] + c_[i] * a_[i+1] ) / ( b_[i] + c_[i] );
	}
	}

	if(b_[nKnots_-1]==0 && c_[nKnots_-1]==0) {
	dydx_[nKnots_-1] = ( a_[nKnots_-1] + an ) / 2.;
	} else {
	dydx_[nKnots_-1] = ( b_[nKnots_-1] * a_[nKnots_-1] + c_[nKnots_-1] * an ) / ( b_[nKnots_-1] + c_[nKnots_-1] );
	}
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 4:09 PM (1 d, 15 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3794163
Default Alt Text: (28 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions