Lau1DCubicSpline.cc

Go to the documentation of this file.

// Copyright University of Warwick 2004 - 2015.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

// Authors:
// Thomas Latham
// John Back
// Paul Harrison

#include <cmath>
#include <cstdlib>
#include <iostream>

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

#include "Lau1DCubicSpline.hh"

Lau1DCubicSpline::Lau1DCubicSpline(const std::vector<Double_t>& xs, const std::vector<Double_t>& ys,
                                   LauSplineType type, 
                                   LauSplineBoundaryType leftBound, 
                                   LauSplineBoundaryType rightBound,
                                   Double_t dydx0, Double_t dydxn) :
        nKnots_(xs.size()),
        x_(xs),
        y_(ys),
        type_(type),
        leftBound_(leftBound),
        rightBound_(rightBound),
        dydx0_(dydx0),
        dydxn_(dydxn)
{
        init();
}

Lau1DCubicSpline::~Lau1DCubicSpline()
{
}

Double_t Lau1DCubicSpline::evaluate(Double_t x) const
{
        // do not attempt to extrapolate the spline
        if( x<x_[0] || x>x_[nKnots_-1] ) {
                std::cout << "WARNING in Lau1DCubicSpline::evaluate : function is only defined between " << x_[0] << " and " << x_[nKnots_-1] << std::endl;
                std::cout << "                                        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
        Int_t cell(0);

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

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

        if(type_ == Lau1DCubicSpline::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
        Double_t t = (x - xLow) / (xHigh - xLow);
        Double_t a =     dydx_[cell]   * (xHigh - xLow) - (yHigh - yLow);
        Double_t b = -1.*dydx_[cell+1] * (xHigh - xLow) + (yHigh - yLow);

        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::updateType(LauSplineType type)
{
        if(type_ != type) {
                type_ = type;
                this->calcDerivatives();
        }
}

void Lau1DCubicSpline::updateBoundaryConditions(LauSplineBoundaryType leftBound,
                LauSplineBoundaryType rightBound,
                Double_t dydx0, Double_t dydxn)
{
        Bool_t updateDerivatives(kFALSE);

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

        if(dydx0_ != dydx0) {
                dydx0_ = dydx0;
                if(leftBound_ == Lau1DCubicSpline::Clamped) updateDerivatives = kTRUE;
        }

        if(dydxn_ != dydxn) {
                dydxn_ = dydxn;
                if(rightBound_ == Lau1DCubicSpline::Clamped) updateDerivatives = kTRUE;
        }

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

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

        dydx_.insert(dydx_.begin(),nKnots_,0.);

        a_.insert(a_.begin(),nKnots_,0.);
        b_.insert(b_.begin(),nKnots_,0.);
        c_.insert(c_.begin(),nKnots_,0.);
        d_.insert(d_.begin(),nKnots_,0.);

        this->calcDerivatives();
}

void Lau1DCubicSpline::calcDerivatives()
{
        switch ( type_ ) {
                case Lau1DCubicSpline::StandardSpline :
                        this->calcDerivativesStandard();
                        break;
                case Lau1DCubicSpline::AkimaSpline :
                        this->calcDerivativesAkima();
                        break;
                case Lau1DCubicSpline::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
        if(leftBound_ == Lau1DCubicSpline::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]));

        } else if(leftBound_ == Lau1DCubicSpline::NotAKnot) {
                // define the width, h, and the 'slope', delta, of the first cell
                Double_t h1(x_[1]-x_[0]), h2(x_[2]-x_[0]);
                Double_t delta1((y_[1]-y_[0])/h1), 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.*h2*h2 + 3.*h1*h2)/(h1+h2) + delta2*5.*h1*h1/(h1+h2);

        } else { //Clamped
                b_[0] = 1.;
                c_[0] = 0.;
                d_[0] = dydx0_;
        }

        // set right boundary condition
        if(rightBound_ == Lau1DCubicSpline::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]));

        } else if(rightBound_ == Lau1DCubicSpline::NotAKnot) {
                // define the width, h, and the 'slope', delta, of the last cell
                Double_t hnm1(x_[nKnots_-1]-x_[nKnots_-2]), hnm2(x_[nKnots_-2]-x_[nKnots_-3]);
                Double_t deltanm1((y_[nKnots_-1]-y_[nKnots_-2])/hnm1), 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] = deltanm2*hnm1*hnm1/(hnm2+hnm1) + deltanm1*(2.*hnm2*hnm2 + 3.*hnm2*hnm1)/(hnm2+hnm1);

        } else { //Clamped
                a_[nKnots_-1] = 0.;
                b_[nKnots_-1] = 1.;
                d_[nKnots_-1] = dydxn_;
        }

        // the remaining equations ensure that f_i-1''(x_i) = f''_i(x_i) for all internal knots
        for(UInt_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(UInt_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_t i=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.), an(0.), anp1(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(UInt_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(UInt_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(UInt_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] );
        }
}