RooCubicSplineKnot Class Reference

#include <RooCubicSplineKnot.h>

Classes
class	BoundaryConditions
class	S2_edge
class	S2_jk
class	S_edge
class	S_jk

Public Member Functions
	RooCubicSplineKnot (const double *array, int nEntries)
template<typename Iter >
	RooCubicSplineKnot (Iter begin, Iter end)
	RooCubicSplineKnot (const std::vector< double > &knots)
double	u (int i) const
int	size () const
double	evaluate (double _u, const RooArgList &b) const
double	analyticalIntegral (const RooArgList &b) const
void	computeCoefficients (std::vector< double > &y, BoundaryConditions bc=BoundaryConditions()) const
void	smooth (std::vector< double > &y, const std::vector< double > &dy, double lambda) const
const std::vector< double > &	knots () const
RooCubicSplineKnot::S_jk	S_jk_sum (int i, const RooArgList &b) const
RooCubicSplineKnot::S2_jk	S2_jk_sum (int i, const RooArgList &b1, const RooArgList &b2) const
RooCubicSplineKnot::S_edge	S_jk_edge (bool left, const RooArgList &b) const
RooCubicSplineKnot::S2_edge	S2_jk_edge (bool left, const RooArgList &b1, const RooArgList &b2) const
double	expIntegral (const TH1 *hist, double gamma, TVectorD &coefficients, TMatrixD &covarianceMatrix) const

Private Member Functions
int	index (double u_) const
double	A (double u_, int i) const
double	B (double u_, int i) const
double	C (double u_, int i) const
double	D (double u_, int i) const
double	ma (int i, bool bc[]) const
double	mc (int i, bool bc[]) const
double	mb (int i, bool bc[]) const
double	P (int i) const
double	Q (int i) const
double	R (int i) const
double	S (int i) const
void	fillPQRS () const
double	d (double u_, int j) const
double	d (double u_, int i, int j, int k) const
double	h (int i, int j) const
double	h (int i) const
double	r (int i) const
double	f (int i) const
double	p (int i) const

Static Private Member Functions
static double	sqr (double x)
static double	cub (double x)
static double	qua (double x)

Private Attributes
std::vector< double >	_u
std::vector< double >	_PQRS
std::vector< double >	_IABCD
std::vector< RooCubicSplineKnot::S_jk >	_S_jk
std::vector< RooCubicSplineKnot::S2_jk >	_S2_jk

Detailed Description

Definition at line 12 of file RooCubicSplineKnot.h.

Constructor & Destructor Documentation

RooCubicSplineKnot() [1/3]

RooCubicSplineKnot::RooCubicSplineKnot ( const double *  array, int  nEntries  )

inline

Definition at line 26 of file RooCubicSplineKnot.h.

: _u( array, array+nEntries) {}

RooCubicSplineKnot() [2/3]

template<typename Iter >

RooCubicSplineKnot::RooCubicSplineKnot ( Iter  begin, Iter  end  )

inline

Definition at line 27 of file RooCubicSplineKnot.h.

: _u(begin,end) {}

RooCubicSplineKnot() [3/3]

RooCubicSplineKnot::RooCubicSplineKnot ( const std::vector< double > &  knots )

inline

Definition at line 28 of file RooCubicSplineKnot.h.

        : _u(knots) {}

Member Function Documentation

A()

double RooCubicSplineKnot::A ( double  u_, int  i  ) const

inlineprivate

Definition at line 210 of file RooCubicSplineKnot.h.

{ return -cub(d(u_,i+1))/P(i); }

analyticalIntegral()

double RooCubicSplineKnot::analyticalIntegral

(

const RooArgList &

b

)

const

Definition at line 147 of file RooCubicSplineKnot.cpp.

                                                                       {
    using RooCubicSplineKnot_aux::push_back;
    using RooCubicSplineKnot_aux::get;
    if (_IABCD.empty()) {
        // the integrals of A,B,C,D from u(i) to u(i+1) only depend on the knot vector...
        // so we create them 'on demand' and cache the result
        _IABCD.reserve(4*size());
        for (int j=0;j<size();++j) {
            push_back(_IABCD,   qua(h(j,j+1))/(4*P(j))
                            , - cub(h(j,j+1))*(3*u(j)-4*u(j-2)+u(j+1))/(12*P(j))
                              - sqr(h(j,j+1))*(3*sqr(u(j))-2*u(j-1)*u(j+1)+sqr(u(j+1))+u(j)*(-4*u(j-1)+2*u(j+1)-4*u(j+2)) +6*u(j-1)*u(j+2)-2*u(j+1)*u(j+2) )/(12*Q(j))
                              + sqr(h(j,j+1))*(3*sqr(u(j+1))+sqr(u(j  ))+2*u(j)*u(j+1)-8*u(j+1)*u(j+2)-4*u(j  )*u(j+2)+6*sqr(u(j+2)))/(12*R(j))
                            ,   sqr(h(j,j+1))*(3*sqr(u(j  ))+sqr(u(j+1))+2*u(j+1)*u(j)-8*u(j  )*u(j-1)-4*u(j-1)*u(j+1)+6*sqr(u(j-1)))/(12*Q(j))
                              - sqr(h(j,j+1))*(3*sqr(u(j+1))+sqr(u(j))-4*u(j-1)*u(j+1)+6*u(j-1)*u(j+2)-4*u(j+1)*u(j+2)-2*u(j)*(u(j-1)-u(j+1)+u(j+2)))/(12*R(j))
                              + cub(h(j,j+1))*(3*u(j+1)-4*u(j+3)+u(j))/(12*S(j))
                            ,   qua(h(j,j+1))/(4*S(j)) );
        }
    }
    // FIXME: this assumes the integration range goes from first knot to last knot...
    assert(b.getSize()-2==size());
    double norm(0);
    for (int i=0; i < size()-1; ++i) for (int k=0;k<4;++k) {
        norm += get(b,i,k)*RooCubicSplineKnot_aux::get(_IABCD,i,k) ;
    }
    return norm;
}

B()

double RooCubicSplineKnot::B ( double  u_, int  i  ) const

inlineprivate

Definition at line 211 of file RooCubicSplineKnot.h.

{ return  sqr(d(u_,i+1))*d(u_,i-2)/P(i) + d(u_,i-1)*d(u_,i+2)*d(u_,i+1)/Q(i) + d(u_,i  )*sqr(d(u_,i+2))/R(i); }

C()

double RooCubicSplineKnot::C ( double  u_, int  i  ) const

inlineprivate

Definition at line 212 of file RooCubicSplineKnot.h.

{ return -sqr(d(u_,i-1))*d(u_,i+1)/Q(i) - d(u_,i  )*d(u_,i+2)*d(u_,i-1)/R(i) - d(u_,i+3)*sqr(d(u_,i  ))/S(i); }

computeCoefficients()

void RooCubicSplineKnot::computeCoefficients	(	std::vector< double > &	y,
		BoundaryConditions	bc = BoundaryConditions()
	)		const

Definition at line 84 of file RooCubicSplineKnot.cpp.

{
 // see http://en.wikipedia.org/wiki/Spline_interpolation
 // for the derivation of the linear system...
 // see http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
 // for the O(N) algorithm to solve the relevant linear system...
    int n = size();
    assert(int(y.size())==size());

    // Weights of first and last spline are fully determined by values at knots
    double bf = y.front();
    double bb = y.back();

    // Set boundary conditions for linear system
    y.front() = bc.value[0] - bf * ( bc.secondDerivative[0] ? ( double(6) / sqr(h(0)) )  : ( -double(3) / h(0)   ) );
    y.back()  = bc.value[1] - bb * ( bc.secondDerivative[1] ? ( double(6) / sqr(h(n-2))) : (  double(3) / h(n-2) ) );

    // Solve linear system
    std::vector<double> c; c.reserve(n);
    c.push_back( mc(0, bc.secondDerivative) / mb(0, bc.secondDerivative) );
    y[0] /=  mb(0, bc.secondDerivative);
    // TODO: rewrite the special cases i==1 and i==n-1 (move latter out of loop)
    //       so we can avoid 'if' statements in the loop body...
    for (int i = 1; i < n; ++i) {
        double m = double(1) / ( mb(i, bc.secondDerivative ) - ma(i, bc.secondDerivative) * c.back() );
        c.push_back( mc(i,bc.secondDerivative) * m );
        y[i] -=  ma(i, bc.secondDerivative) * y[i - 1];
        y[i] *=  m;
    }
    for (int i = n-1 ; i-- > 0; ) y[i] -= c[i] * y[i + 1];
    y.push_back(bb); y.insert(y.begin(),bf); // ouch... expensive!
}

cub()

static double RooCubicSplineKnot::cub ( double  x )

inlinestaticprivate

Definition at line 252 of file RooCubicSplineKnot.h.

{ return x*sqr(x); }

D()

double RooCubicSplineKnot::D ( double  u_, int  i  ) const

inlineprivate

Definition at line 213 of file RooCubicSplineKnot.h.

{ return  cub(d(u_,i  ))/S(i); }

d() [1/2]

double RooCubicSplineKnot::d ( double  u_, int  j  ) const

inlineprivate

Definition at line 254 of file RooCubicSplineKnot.h.

{ return u_-u(j); }

d() [2/2]

double RooCubicSplineKnot::d ( double  u_, int  i, int  j, int  k  ) const

inlineprivate

Definition at line 255 of file RooCubicSplineKnot.h.

{ return d(u_,i)*d(u_,j)*d(u_,k); }

evaluate()

double RooCubicSplineKnot::evaluate	(	double	_u,
		const RooArgList &	b
	)		const

Definition at line 133 of file RooCubicSplineKnot.cpp.

                                                                       {
    using RooCubicSplineKnot_aux::get;
    int i = index(x); // location in knot vector: index of last knot less than x...
    assert(-1<=i && i<=size());
    // we're a 'natural' spline. If beyond first/last knot, Extrapolate using the derivative at the first (last) knot
    if (i==-1    ) { ++i; return get(b,i,0) - d(x,i)*r(i  )*(get(b,i,0)-get(b,i,1)); }
    if (i==size()) { --i; return get(b,i,2) + d(x,i)*r(i-1)*(get(b,i,2)-get(b,i,1)); }
    assert( u(i) <= x && x<= u(i+1) );
    return get(b,i,0)*A(x,i) // TODO: substitute A,B,C,D 'in situ'
         + get(b,i,1)*B(x,i)
         + get(b,i,2)*C(x,i)
         + get(b,i,3)*D(x,i);
}

expIntegral()

double RooCubicSplineKnot::expIntegral	(	const TH1 *	hist,
		double	gamma,
		TVectorD &	coefficients,
		TMatrixD &	covarianceMatrix
	)		const

Definition at line 367 of file RooCubicSplineKnot.cpp.

                                                                                                                               {
    assert(hist!=0);
    int nknots = size();
    int nsplines = nknots+2;
    int nbins = hist->GetNbinsX();
    TMatrixD matrix(nsplines,nbins); for (int i=0;i<nsplines;++i) for (int j=0;j<nbins;++j) { matrix(i,j)=0; }
    TVectorD Y(nbins), DY(nbins);

    for (int i=0;i<nbins ;++i) {
        Y(i)  = hist->GetBinContent(1+i);
        DY(i) = hist->GetBinError(1+i);


        Y(i) /= DY(i)*DY(i);

        double lo = hist->GetBinLowEdge(1+i);
        double hi = lo + hist->GetBinWidth(1+i);

        double Norm =   (exp(-gamma*lo)-exp(-gamma*hi))/gamma;

        // go for flat parent distribution;
        // Norm = hi-lo;

#if 0
        // compute contributions beyond the first and last knot...
        if (lo < _u.front()) {
                // b-spline 0:  1 - d(x,0)*r(0)  // d(x,0) = x - u(0)
                // b-spline 1:    + d(x,0)*r(0)
                matrix(0, 0) += eIn( lo,hi, 1., -r(0), u(0) )/Norm;
                matrix(1, 0) += eIn( lo,hi, 0., +r(0), u(0) )/Norm;

        }
        if (hi > _u.back()) {
                //                                           c + alpha  *    (x - x0)
                matrix(nknots  , nbins-1) += eIn( lo, hi, 0., -r(nknots-2), u(nknots-1) )/Norm;
                matrix(nknots+1, nbins-1) += eIn( lo, hi, 1., +r(nknots-2), u(nknots-1) )/Norm;

        }
#else
        //  assume nothing beyond spline interval
        assert( lo >= _u.front() && hi <= _u.back() );
#endif

        // this bit is by construction fully contained between the first and last knot
        for (int j=0;j<nknots-1;++j) {
            double l = std::max(lo,u(j));
            double h = std::min(hi,u(j+1));
            if (l>=h) continue;
            // in the knot interval [ u(j),  u(j+1) ], the splines j..j+3 contribute...
            matrix(j  ,i) += (-eI(l,h, u(j+1),u(j+1),u(j+1), gamma)/P(j) ) / Norm;
            matrix(j+1,i) += ( eI(l,h, u(j-2),u(j+1),u(j+1), gamma)/P(j)
                              +eI(l,h, u(j-1),u(j+1),u(j+2), gamma)/Q(j)
                              +eI(l,h, u(j  ),u(j+2),u(j+2), gamma)/R(j) ) / Norm;
            matrix(j+2,i) += (-eI(l,h, u(j-1),u(j-1),u(j+1), gamma)/Q(j)
                              -eI(l,h, u(j-1),u(j  ),u(j+2), gamma)/R(j)
                              -eI(l,h, u(j  ),u(j  ),u(j+3), gamma)/S(j) ) / Norm;
            matrix(j+3,i) += ( eI(l,h, u(j  ),u(j  ),u(j  ), gamma)/S(j) ) / Norm;
        }

    }

    TMatrixDSym D(nsplines);
    for (int i=0;i<nsplines;++i) for (int j=0;j<=i;++j) {
        D(i,j) = 0;
        for (int k=0;k<nbins;++k) D(i,j) += matrix(i,k)*matrix(j,k)/(DY(k)*DY(k));
        if (i!=j) D(j,i) = D(i,j);
    }
    TDecompBK solver(D);

    Y *= matrix;
    TVectorD YY = Y;

    bool ok;
    coefficients.ResizeTo(nsplines);
    coefficients = solver.Solve(Y,ok);
    if (!ok) {
        std::cout << "WARNING: bad linear system solution... " << std::endl;
        return -1;
    }

    D.Invert();
    covarianceMatrix.ResizeTo(D);
    covarianceMatrix = D;

    double chisq = 0;
    for (int i=0;i<nbins;++i) {
        double y = hist->GetBinContent(1+i);
        double dy = hist->GetBinError(1+i);
        double f = 0; for (int j=0;j<nsplines;++j) f+=matrix(j,i)*coefficients(j);
        double c = sqr( (y-f)/dy );
        // double lo = hist->GetBinLowEdge(1+i);
        // double hi = lo + hist->GetBinWidth(1+i);
        // std::cout << " bin " << i  << " [  " << lo << " , " << hi << " ]  y = " << y << " dy = " << dy << " f = " << f << " chi^2 = " << c << std::endl;
        chisq+=c;
    }

    return chisq;
    // silence compiler warning about fI being unused
    (void) fI;
}

f()

double RooCubicSplineKnot::f ( int  i ) const

inlineprivate

Definition at line 259 of file RooCubicSplineKnot.h.

{ return -r(i-1)-r(i); }

fillPQRS()

void RooCubicSplineKnot::fillPQRS ( ) const

private

Definition at line 120 of file RooCubicSplineKnot.cpp.

                                        {
    assert(_PQRS.empty());
    // P,Q,R,S only depend on the knot vector, so build at construction, and cache them...
    _PQRS.reserve(4*size());
    for (int i=0;i<size();++i) {
        _PQRS.push_back( h(i+1,i-2)*h(i+1,i-1)*h(i+1,i) );
        _PQRS.push_back( h(i+1,i-1)*h(i+2,i-1)*h(i+1,i) );
        _PQRS.push_back( h(i+2,i  )*h(i+2,i-1)*h(i+1,i) );
        _PQRS.push_back( h(i+2,i  )*h(i+3,i  )*h(i+1,i) );
    }

}

h() [1/2]

double RooCubicSplineKnot::h ( int  i, int  j  ) const

inlineprivate

Definition at line 256 of file RooCubicSplineKnot.h.

{ return u(i)-u(j); }

h() [2/2]

double RooCubicSplineKnot::h ( int  i ) const

inlineprivate

Definition at line 257 of file RooCubicSplineKnot.h.

{ return h(i+1,i); }

index()

int RooCubicSplineKnot::index ( double  u_ ) const

private

Definition at line 175 of file RooCubicSplineKnot.cpp.

{
    if (u>_u.back()) return size();
    std::vector<double>::const_iterator i = --std::upper_bound(_u.begin(),_u.end()-1,u);
    return std::distance(_u.begin(),i);
}

knots()

const std::vector<double>& RooCubicSplineKnot::knots ( ) const

inline

Definition at line 43 of file RooCubicSplineKnot.h.

{ return _u; }

ma()

double RooCubicSplineKnot::ma ( int  i, bool  bc[]  ) const

inlineprivate

Definition at line 216 of file RooCubicSplineKnot.h.

                                      {  // subdiagonal

        return (i!=size()-1)  ?  A(u(i),i)
             : (bc[1] )       ? double(6)/(h(i,i-2)*h(i-1))
                              : double(0) ;
    }

mb()

double RooCubicSplineKnot::mb ( int  i, bool  bc[]  ) const

inlineprivate

Definition at line 229 of file RooCubicSplineKnot.h.

                                      {  // diagonal
        if (i!=0 && i!=size()-1) return B(u(i),i);

        if (i==0) {
            return bc[0] ? -(double(6)/h(i)+double(6)/h(i+2,i))/h(i)
                         :   double(3)/h(i);
        } else  {
            return bc[1] ? -(double(6)/h(i-1)+double(6)/h(i,i-2))/h(i-1)
                         : - double(3)/h(i-1);
        }
    }

mc()

double RooCubicSplineKnot::mc ( int  i, bool  bc[]  ) const

inlineprivate

Definition at line 223 of file RooCubicSplineKnot.h.

                                      {  // superdiagonal
          return (i!=0) ? C(u(i),i)
               :  bc[0] ? double(6)/(h(2,0)*h(0))
                        : double(0);
    }

P()

double RooCubicSplineKnot::P ( int  i ) const

inlineprivate

Definition at line 244 of file RooCubicSplineKnot.h.

{ if (_PQRS.empty()) fillPQRS(); assert(4*i  <int(_PQRS.size())); return  _PQRS[4*i  ]; }

p()

double RooCubicSplineKnot::p ( int  i ) const

inlineprivate

Definition at line 260 of file RooCubicSplineKnot.h.

{ return 2*h(i-1)+h(i); }

Q()

double RooCubicSplineKnot::Q ( int  i ) const

inlineprivate

Definition at line 245 of file RooCubicSplineKnot.h.

{ if (_PQRS.empty()) fillPQRS(); assert(4*i+1<int(_PQRS.size())); return  _PQRS[4*i+1]; }

qua()

static double RooCubicSplineKnot::qua ( double  x )

inlinestaticprivate

Definition at line 253 of file RooCubicSplineKnot.h.

{ return sqr(sqr(x)); }

R()

double RooCubicSplineKnot::R ( int  i ) const

inlineprivate

Definition at line 246 of file RooCubicSplineKnot.h.

{ if (_PQRS.empty()) fillPQRS(); assert(4*i+2<int(_PQRS.size())); return  _PQRS[4*i+2]; }

r()

double RooCubicSplineKnot::r ( int  i ) const

inlineprivate

Definition at line 258 of file RooCubicSplineKnot.h.

{ return double(3)/h(i); }

S()

double RooCubicSplineKnot::S ( int  i ) const

inlineprivate

Definition at line 247 of file RooCubicSplineKnot.h.

{ if (_PQRS.empty()) fillPQRS(); assert(4*i+3<int(_PQRS.size())); return  _PQRS[4*i+3]; }

S2_jk_edge()

RooCubicSplineKnot::S2_edge RooCubicSplineKnot::S2_jk_edge	(	bool	left,
		const RooArgList &	b1,
		const RooArgList &	b2
	)		const

Definition at line 297 of file RooCubicSplineKnot.cpp.

                                                                                                                    {
       using RooCubicSplineKnot_aux::get;
       if (left) {
         double alpha1 = -r(0)*(get(b1,0,0)-get(b1,0,1));
         double beta1 = evaluate(u(0),b1)-alpha1*u(0);
         double alpha2 = -r(0)*(get(b2,0,0)-get(b2,0,1));
         double beta2 = evaluate(u(0),b2)-alpha2*u(0);
         return S2_edge(alpha1,beta1,alpha2,beta2);
       } else {
         int i = size()-1;
         double alpha1 = r(i-1)*(get(b1,i,2)-get(b1,i,1));
         double beta1 = evaluate(u(i),b1)-alpha1*u(i);
         double alpha2 = r(i-1)*(get(b2,i,2)-get(b2,i,1));
         double beta2 = evaluate(u(i),b2)-alpha2*u(i);
         return S2_edge(alpha1,beta1,alpha2,beta2);
       }
}

S2_jk_sum()

RooCubicSplineKnot::S2_jk RooCubicSplineKnot::S2_jk_sum	(	int	i,
		const RooArgList &	b1,
		const RooArgList &	b2
	)		const

Definition at line 207 of file RooCubicSplineKnot.cpp.

{
    if (_S2_jk.empty()) {
        _S2_jk.reserve(size()*10);
        for(int j=0;j<size();++j) {
            // This 'table' should be compatible with the definitions of A,B,C, and D...
            _S2_jk.push_back(  RooCubicSplineKnot::S2_jk(u(j+1),u(j+1),u(j+1),u(j+1),u(j+1),u(j+1)) / (P(j)*P(j)) ); // (A * A)

            _S2_jk.push_back(  RooCubicSplineKnot::S2_jk(u(j-2),u(j+1),u(j+1),u(j-2),u(j+1),u(j+1)) / (P(j)*P(j))
                              +RooCubicSplineKnot::S2_jk(u(j-1),u(j+1),u(j+2),u(j-1),u(j+1),u(j+2)) / (Q(j)*Q(j))
                              +RooCubicSplineKnot::S2_jk(u(j  ),u(j+2),u(j+2),u(j  ),u(j+2),u(j+2)) / (R(j)*R(j))
                              +RooCubicSplineKnot::S2_jk(u(j-2),u(j+1),u(j+1),u(j-1),u(j+1),u(j+2)) / (P(j)*Q(j)) * 2.
                              +RooCubicSplineKnot::S2_jk(u(j-2),u(j+1),u(j+1),u(j  ),u(j+2),u(j+2)) / (P(j)*R(j)) * 2.
                              +RooCubicSplineKnot::S2_jk(u(j-1),u(j+1),u(j+2),u(j  ),u(j+2),u(j+2)) / (Q(j)*R(j)) * 2. ); // (B * B)

            _S2_jk.push_back(  RooCubicSplineKnot::S2_jk(u(j-1),u(j-1),u(j+1),u(j-1),u(j-1),u(j+1)) / (Q(j)*Q(j))
                              +RooCubicSplineKnot::S2_jk(u(j-1),u(j  ),u(j+2),u(j-1),u(j  ),u(j+2)) / (R(j)*R(j))
                              +RooCubicSplineKnot::S2_jk(u(j  ),u(j  ),u(j+3),u(j  ),u(j  ),u(j+3)) / (S(j)*S(j))
                              +RooCubicSplineKnot::S2_jk(u(j-1),u(j-1),u(j+1),u(j-1),u(j  ),u(j+2)) / (Q(j)*R(j)) * 2.
                              +RooCubicSplineKnot::S2_jk(u(j-1),u(j-1),u(j+1),u(j  ),u(j  ),u(j+3)) / (Q(j)*S(j)) * 2.
                              +RooCubicSplineKnot::S2_jk(u(j-1),u(j  ),u(j+2),u(j  ),u(j  ),u(j+3)) / (R(j)*S(j)) * 2. ); // (C * C)

            _S2_jk.push_back(  RooCubicSplineKnot::S2_jk(u(j  ),u(j  ),u(j  ),u(j  ),u(j  ),u(j  )) / (S(j)*S(j)) ); // (D * D)


            _S2_jk.push_back( -RooCubicSplineKnot::S2_jk(u(j+1),u(j+1),u(j+1),u(j-2),u(j+1),u(j+1)) / (P(j)*P(j))
                              -RooCubicSplineKnot::S2_jk(u(j+1),u(j+1),u(j+1),u(j-1),u(j+1),u(j+2)) / (P(j)*Q(j))
                              -RooCubicSplineKnot::S2_jk(u(j+1),u(j+1),u(j+1),u(j  ),u(j+2),u(j+2)) / (P(j)*R(j)) ); // (A * B)

            _S2_jk.push_back(  RooCubicSplineKnot::S2_jk(u(j+1),u(j+1),u(j+1),u(j-1),u(j-1),u(j+1)) / (P(j)*Q(j))
                              +RooCubicSplineKnot::S2_jk(u(j+1),u(j+1),u(j+1),u(j-1),u(j  ),u(j+2)) / (P(j)*R(j))
                              +RooCubicSplineKnot::S2_jk(u(j+1),u(j+1),u(j+1),u(j  ),u(j  ),u(j+3)) / (P(j)*S(j)) ); // (A * C)

            _S2_jk.push_back( -RooCubicSplineKnot::S2_jk(u(j+1),u(j+1),u(j+1),u(j  ),u(j  ),u(j  )) / (P(j)*S(j)) ); // (A * D)


            _S2_jk.push_back( -RooCubicSplineKnot::S2_jk(u(j-2),u(j+1),u(j+1),u(j-1),u(j-1),u(j+1)) / (P(j)*Q(j))
                              -RooCubicSplineKnot::S2_jk(u(j-2),u(j+1),u(j+1),u(j-1),u(j  ),u(j+2)) / (P(j)*R(j))
                              -RooCubicSplineKnot::S2_jk(u(j-2),u(j+1),u(j+1),u(j  ),u(j  ),u(j+3)) / (P(j)*S(j))
                              -RooCubicSplineKnot::S2_jk(u(j-1),u(j+1),u(j+2),u(j-1),u(j-1),u(j+1)) / (Q(j)*Q(j))
                              -RooCubicSplineKnot::S2_jk(u(j-1),u(j+1),u(j+2),u(j-1),u(j  ),u(j+2)) / (Q(j)*R(j))
                              -RooCubicSplineKnot::S2_jk(u(j-1),u(j+1),u(j+2),u(j  ),u(j  ),u(j+3)) / (Q(j)*S(j))
                              -RooCubicSplineKnot::S2_jk(u(j  ),u(j+2),u(j+2),u(j-1),u(j-1),u(j+1)) / (R(j)*Q(j))
                              -RooCubicSplineKnot::S2_jk(u(j  ),u(j+2),u(j+2),u(j-1),u(j  ),u(j+2)) / (R(j)*R(j))
                              -RooCubicSplineKnot::S2_jk(u(j  ),u(j+2),u(j+2),u(j  ),u(j  ),u(j+3)) / (R(j)*S(j)) ); // (B * C)

            _S2_jk.push_back(  RooCubicSplineKnot::S2_jk(u(j-2),u(j+1),u(j+1),u(j  ),u(j  ),u(j  )) / (P(j)*S(j))
                              +RooCubicSplineKnot::S2_jk(u(j-1),u(j+1),u(j+2),u(j  ),u(j  ),u(j  )) / (Q(j)*S(j))
                              +RooCubicSplineKnot::S2_jk(u(j  ),u(j+2),u(j+2),u(j  ),u(j  ),u(j  )) / (R(j)*S(j)) ); // (B * D)

            _S2_jk.push_back( -RooCubicSplineKnot::S2_jk(u(j-1),u(j-1),u(j+1),u(j  ),u(j  ),u(j  )) / (Q(j)*S(j))
                              -RooCubicSplineKnot::S2_jk(u(j-1),u(j  ),u(j+2),u(j  ),u(j  ),u(j  )) / (R(j)*S(j))
                              -RooCubicSplineKnot::S2_jk(u(j  ),u(j  ),u(j+3),u(j  ),u(j  ),u(j  )) / (S(j)*S(j)) ); // (C * D)
        }
    }
    using RooCubicSplineKnot_aux::get;
    using RooCubicSplineKnot_aux::get2;
    return get2(_S2_jk,i,0)*get(b1,i,0)*get(b2,i,0)
         + get2(_S2_jk,i,1)*get(b1,i,1)*get(b2,i,1)
         + get2(_S2_jk,i,2)*get(b1,i,2)*get(b2,i,2)
         + get2(_S2_jk,i,3)*get(b1,i,3)*get(b2,i,3)
         + get2(_S2_jk,i,4)*(get(b1,i,0)*get(b2,i,1)+get(b1,i,1)*get(b2,i,0))
         + get2(_S2_jk,i,5)*(get(b1,i,0)*get(b2,i,2)+get(b1,i,2)*get(b2,i,0))
         + get2(_S2_jk,i,6)*(get(b1,i,0)*get(b2,i,3)+get(b1,i,3)*get(b2,i,0))
         + get2(_S2_jk,i,7)*(get(b1,i,1)*get(b2,i,2)+get(b1,i,2)*get(b2,i,1))
         + get2(_S2_jk,i,8)*(get(b1,i,1)*get(b2,i,3)+get(b1,i,3)*get(b2,i,1))
         + get2(_S2_jk,i,9)*(get(b1,i,2)*get(b2,i,3)+get(b1,i,3)*get(b2,i,2));
}

S_jk_edge()

RooCubicSplineKnot::S_edge RooCubicSplineKnot::S_jk_edge	(	bool	left,
		const RooArgList &	b
	)		const

Definition at line 281 of file RooCubicSplineKnot.cpp.

                                                                                           {
       using RooCubicSplineKnot_aux::get;
       if (left) {
         // efficiency = return evaluate(u(0),b) - d(x,0)*r(0)*(get(b,0,0)-get(b,0,1));
         double alpha = -r(0)*(get(b,0,0)-get(b,0,1));
         double beta = evaluate(u(0),b)-alpha*u(0);
         return S_edge( alpha, beta );
       } else {
         int i = size()-1;
         double alpha = r(i-1)*(get(b,i,2)-get(b,i,1));
         double beta = evaluate(u(i),b)-alpha*u(i);
         return S_edge(alpha,beta);
       }
}

S_jk_sum()

RooCubicSplineKnot::S_jk RooCubicSplineKnot::S_jk_sum	(	int	i,
		const RooArgList &	b
	)		const

Definition at line 183 of file RooCubicSplineKnot.cpp.

{
    if (_S_jk.empty()) {
        _S_jk.reserve(size()*4);
        for(int j=0;j<size();++j) {
            // This 'table' should be compatible with the definitions of A,B,C, and D...
            _S_jk.push_back( -RooCubicSplineKnot::S_jk(u(j+1),u(j+1),u(j+1))/P(j) ); // A
            _S_jk.push_back(  RooCubicSplineKnot::S_jk(u(j-2),u(j+1),u(j+1))/P(j)
                             +RooCubicSplineKnot::S_jk(u(j-1),u(j+1),u(j+2))/Q(j)
                             +RooCubicSplineKnot::S_jk(u(j  ),u(j+2),u(j+2))/R(j) ); // B
            _S_jk.push_back( -RooCubicSplineKnot::S_jk(u(j-1),u(j-1),u(j+1))/Q(j)
                             -RooCubicSplineKnot::S_jk(u(j-1),u(j  ),u(j+2))/R(j)
                             -RooCubicSplineKnot::S_jk(u(j  ),u(j  ),u(j+3))/S(j) ); // C
            _S_jk.push_back(  RooCubicSplineKnot::S_jk(u(j  ),u(j  ),u(j  ))/S(j) ); // D
        }
    }
    using RooCubicSplineKnot_aux::get;
    return get(_S_jk,i,0)*get(b,i,0)
         + get(_S_jk,i,1)*get(b,i,1)
         + get(_S_jk,i,2)*get(b,i,2)
         + get(_S_jk,i,3)*get(b,i,3);
}

size()

int RooCubicSplineKnot::size ( ) const

inline

Definition at line 36 of file RooCubicSplineKnot.h.

{ return _u.size(); }

smooth()

void RooCubicSplineKnot::smooth	(	std::vector< double > &	y,
		const std::vector< double > &	dy,
		double	lambda
	)		const

Definition at line 26 of file RooCubicSplineKnot.cpp.

                                                                                                      {
// implementation as described in D.S.G. Pollock, Smoothing Splines"
// see http://r.789695.n4.nabble.com/file/n905996/SPLINES.PDF
    using namespace std;
    int n = y.size();
    vector<double> uu(n-2),vv(n-3),ww(n-4), q(n-2);
    assert( dy.size()==y.size());
    assert( size()==n);
    // lambda = 0 : no smoothing ; lambda -> 1: straight line (the ultimate smooth curve)
    if (lambda<0|| !(lambda<1)) {
        throw std::string("RooCubicSplineKnot::smooth: smoothing parameter must be in range [0,1)");
    }
    double mu =2*lambda/(3*(1-lambda));
    for (int i=0;i<n-2;++i) {
        q[i]  = r(i)*(y[i]-y[i+1])-r(i+1)*(y[i+1]-y[i+2]);
        uu[i] = p(i+1)+mu*(sqr(r(i)*dy[i])+sqr(f(i+1)*dy[i+1])+sqr(r(i+1)*dy[i+2]) );
        if (i>n-4) continue;
        vv[i] = h(i+1)+mu*(f(i+1)*r(i+1)*sqr(dy[i+1]) + f(i+1)*r(i+2)*sqr(dy[i+2]));
        if (i>n-5) continue ;
        ww[i] =        mu*r(i+1)*r(i+2)*sqr(dy[i+2]);
    }

    // Solve A b = q  with  A(i,i)=uu(i), A(i+1,i)=A(i,i+1) = vv(i), A(i+2,i),A(i,i+2)=ww(i)
    // start with factorization: A -> L D LT
    vv[0] /= uu[0]; ww[0] /= uu[0];
    uu[1] -= uu[0]*sqr(vv[0]);
    vv[1] -= uu[0]*vv[0]*ww[0];
    vv[1] /= uu[1]; ww[1] /= uu[1];
    for (int i=2;i<n-4;++i) {
        uu[i] -= uu[i-1]*sqr(vv[i-1])+uu[i-2]*sqr(ww[i-2]);
        vv[i] -= uu[i-1]*vv[i-1]*ww[i-1];  vv[i] /= uu[i];
        ww[i] /= uu[i];
    }
    uu[n-4] -= uu[n-5]*sqr(vv[n-5])+uu[n-6]*sqr(ww[n-6]);
    vv[n-4] -= uu[n-5]*vv[n-5]*ww[n-5];  vv[n-4] /= uu[n-4];
    uu[n-3] -= uu[n-4]*sqr(vv[n-4])+uu[n-5]*sqr(ww[n-5]);

    // forward substitution (i.e. solve L ( D LT  b )  = q  for D LT b -> q
    q[1]                         -= vv[0  ]*q[0  ];
    for (int i=2;i<n-2;++i) q[i] -= vv[i-1]*q[i-1]+ww[i-2]*q[i-2];
    // rescale  (i.e. solve D (LT b ) = q  for LT b
    for (int i=0;i<n-2;++i) { q[i] /= uu[i] ; }
    // backward substitution (i.e. solve LT b = q  for b -> q )
    q[n-4]                          -= vv[n-4]*q[n-3];
    for (int i=n-5;i>=0;--i) q[i  ] -= vv[i  ]*q[i+1]+ww[i]*q[i+2];
    // solve for y...
    int i=0;
    y[i]                 -=mu*sqr(dy[i])*(r(i)*q[i]                          ); ++i;
    y[i]                 -=mu*sqr(dy[i])*(r(i)*q[i]+f(i)*q[i-1]              ); ++i;
    for (;i<n-2;++i) y[i]-=mu*sqr(dy[i])*(r(i)*q[i]+f(i)*q[i-1]+r(i-1)*q[i-2]);
    y[i]                 -=mu*sqr(dy[i])*(          f(i)*q[i-1]+r(i-1)*q[i-2]); ++i;
    y[i]                 -=mu*sqr(dy[i])*(                      r(i-1)*q[i-2]);

}

sqr()

static double RooCubicSplineKnot::sqr ( double  x )

inlinestaticprivate

Definition at line 251 of file RooCubicSplineKnot.h.

{ return x*x; }

u()

double RooCubicSplineKnot::u ( int  i ) const

inline

Definition at line 31 of file RooCubicSplineKnot.h.

                          {
        assert(size());
        assert(i>-3);
        assert(i<size()+3);
        return _u[std::min(std::max(0,i),size()-1)]; }

Member Data Documentation

_IABCD

std::vector<double> RooCubicSplineKnot::_IABCD

mutableprivate

Definition at line 264 of file RooCubicSplineKnot.h.

_PQRS

std::vector<double> RooCubicSplineKnot::_PQRS

mutableprivate

Definition at line 263 of file RooCubicSplineKnot.h.

_S2_jk

std::vector<RooCubicSplineKnot::S2_jk> RooCubicSplineKnot::_S2_jk

mutableprivate

Definition at line 266 of file RooCubicSplineKnot.h.

_S_jk

std::vector<RooCubicSplineKnot::S_jk> RooCubicSplineKnot::_S_jk

mutableprivate

Definition at line 265 of file RooCubicSplineKnot.h.

_u

std::vector<double> RooCubicSplineKnot::_u

private

Definition at line 262 of file RooCubicSplineKnot.h.

---

The documentation for this class was generated from the following files:

• Mint/RooCubicSplineKnot.h
• src/Mojito/TimeDependent/RooCubicSplineKnot.cpp