Warning, /include/Geant4/tools/spline is written in an unsupported language. File is not indexed.
0001 // Copyright (C) 2010, Guy Barrand. All rights reserved.
0002 // See the file tools.license for terms.
0003
0004 #ifndef tools_spline
0005 #define tools_spline
0006
0007 // From Federico Carminati code found in root-6.08.06/TSpline.h, TSpline.cxx.
0008
0009 #include "mnmx"
0010 #include <cstddef>
0011 #include <vector>
0012 #include <ostream>
0013 #include <cmath>
0014
0015 namespace tools {
0016 namespace spline {
0017
0018 class base_poly {
0019 public:
0020 base_poly():fX(0),fY(0) {}
0021 base_poly(double x,double y):fX(x),fY(y) {}
0022 virtual ~base_poly(){}
0023 public:
0024 base_poly(base_poly const &a_from):fX(a_from.fX),fY(a_from.fY) {}
0025 base_poly& operator=(base_poly const &a_from) {
0026 if(this==&a_from) return *this;
0027 fX = a_from.fX;
0028 fY = a_from.fY;
0029 return *this;
0030 }
0031 public:
0032 const double& X() const {return fX;}
0033 const double& Y() const {return fY;}
0034 double &X() {return fX;}
0035 double &Y() {return fY;}
0036 protected:
0037 double fX; // abscissa
0038 double fY; // constant term
0039 };
0040
0041 class cubic_poly : public base_poly {
0042 public:
0043 cubic_poly():fB(0), fC(0), fD(0) {}
0044 cubic_poly(double x, double y, double b, double c, double d):base_poly(x,y), fB(b), fC(c), fD(d) {}
0045 public:
0046 cubic_poly(cubic_poly const &a_from)
0047 :base_poly(a_from), fB(a_from.fB), fC(a_from.fC), fD(a_from.fD) {}
0048 cubic_poly& operator=(cubic_poly const &a_from) {
0049 if(this==&a_from) return *this;
0050 base_poly::operator=(a_from);
0051 fB = a_from.fB;
0052 fC = a_from.fC;
0053 fD = a_from.fD;
0054 return *this;
0055 }
0056 public:
0057 double &B() {return fB;}
0058 double &C() {return fC;}
0059 double &D() {return fD;}
0060 double eval(double x) const {double dx=x-fX;return (fY+dx*(fB+dx*(fC+dx*fD)));}
0061 protected:
0062 double fB; // first order expansion coefficient : fB*1! is the first derivative at x
0063 double fC; // second order expansion coefficient : fC*2! is the second derivative at x
0064 double fD; // third order expansion coefficient : fD*3! is the third derivative at x
0065 };
0066
0067 ////////////////////////////////////////////////////////////////////////////////////
0068 ////////////////////////////////////////////////////////////////////////////////////
0069 ////////////////////////////////////////////////////////////////////////////////////
0070 class base_spline {
0071 protected:
0072 base_spline(std::ostream& a_out):m_out(a_out), fDelta(-1), fXmin(0), fXmax(0), fNp(0), fKstep(false) {}
0073 public:
0074 base_spline(std::ostream& a_out,double delta, double xmin, double xmax, size_t np, bool step)
0075 :m_out(a_out),fDelta(delta), fXmin(xmin),fXmax(xmax), fNp(np), fKstep(step)
0076 {}
0077 virtual ~base_spline() {}
0078 protected:
0079 base_spline(const base_spline& a_from)
0080 :m_out(a_from.m_out)
0081 ,fDelta(a_from.fDelta),fXmin(a_from.fXmin),fXmax(a_from.fXmax),fNp(a_from.fNp),fKstep(a_from.fKstep) {}
0082 base_spline& operator=(const base_spline& a_from) {
0083 if(this==&a_from) return *this;
0084 fDelta=a_from.fDelta;
0085 fXmin=a_from.fXmin;
0086 fXmax=a_from.fXmax;
0087 fNp=a_from.fNp;
0088 fKstep=a_from.fKstep;
0089 return *this;
0090 }
0091 protected:
0092 std::ostream& m_out;
0093 double fDelta; // Distance between equidistant knots
0094 double fXmin; // Minimum value of abscissa
0095 double fXmax; // Maximum value of abscissa
0096 size_t fNp; // Number of knots
0097 bool fKstep; // True of equidistant knots
0098 };
0099
0100
0101 //////////////////////////////////////////////////////////////////////////
0102 // //
0103 // cubic //
0104 // //
0105 // Class to create third splines to interpolate knots //
0106 // Arbitrary conditions can be introduced for first and second //
0107 // derivatives at beginning and ending points //
0108 // //
0109 //////////////////////////////////////////////////////////////////////////
0110
0111 class cubic : public base_spline {
0112 protected:
0113 cubic(std::ostream& a_out) : base_spline(a_out) , fPoly(0), fValBeg(0), fValEnd(0), fBegCond(-1), fEndCond(-1) {}
0114 public:
0115 cubic(std::ostream& a_out,size_t a_n,double a_x[], double a_y[], double a_valbeg = 0, double a_valend = 0)
0116 :base_spline(a_out,-1,0,0,a_n,false)
0117 ,fValBeg(a_valbeg), fValEnd(a_valend), fBegCond(0), fEndCond(0)
0118 {
0119 if(!a_n) {
0120 m_out << "tools::spline::cubic : a_np is null." << std::endl;
0121 return;
0122 }
0123 fXmin = a_x[0];
0124 fXmax = a_x[a_n-1];
0125 fPoly.resize(a_n);
0126 for (size_t i=0; i<a_n; ++i) {
0127 fPoly[i].X() = a_x[i];
0128 fPoly[i].Y() = a_y[i];
0129 }
0130 build_coeff(); // Build the spline coefficients
0131 }
0132 public:
0133 cubic(const cubic& a_from)
0134 :base_spline(a_from)
0135 ,fPoly(a_from.fPoly),fValBeg(a_from.fValBeg),fValEnd(a_from.fValEnd),fBegCond(a_from.fBegCond),fEndCond(a_from.fEndCond)
0136 {}
0137 cubic& operator=(const cubic& a_from) {
0138 if(this==&a_from) return *this;
0139 base_spline::operator=(a_from);
0140 fPoly = a_from.fPoly;
0141 fValBeg=a_from.fValBeg;
0142 fValEnd=a_from.fValEnd;
0143 fBegCond=a_from.fBegCond;
0144 fEndCond=a_from.fEndCond;
0145 return *this;
0146 }
0147 public:
0148 double eval(double x) const {
0149 if(!fNp) return 0;
0150 // Eval this spline at x
0151 size_t klow = find_x(x);
0152 if ( (fNp > 1) && (klow >= (fNp-1))) klow = fNp-2; //see: https://savannah.cern.ch/bugs/?71651
0153 return fPoly[klow].eval(x);
0154 }
0155 protected:
0156 template<typename T>
0157 static int TMath_Nint(T x) {
0158 // Round to nearest integer. Rounds half integers to the nearest even integer.
0159 int i;
0160 if (x >= 0) {
0161 i = int(x + 0.5);
0162 if ( i & 1 && x + 0.5 == T(i) ) i--;
0163 } else {
0164 i = int(x - 0.5);
0165 if ( i & 1 && x - 0.5 == T(i) ) i++;
0166 }
0167 return i;
0168 }
0169 static int TMath_FloorNint(double x) { return TMath_Nint(::floor(x)); }
0170
0171 size_t find_x(double x) const {
0172 int klow=0, khig=int(fNp-1);
0173 //
0174 // If out of boundaries, extrapolate
0175 // It may be badly wrong
0176 if(x<=fXmin) klow=0;
0177 else if(x>=fXmax) klow=khig;
0178 else {
0179 if(fKstep) { // Equidistant knots, use histogramming :
0180 klow = TMath_FloorNint((x-fXmin)/fDelta);
0181 // Correction for rounding errors
0182 if (x < fPoly[klow].X())
0183 klow = max_of<int>(klow-1,0);
0184 else if (klow < khig) {
0185 if (x > fPoly[klow+1].X()) ++klow;
0186 }
0187 } else {
0188 int khalf;
0189 //
0190 // Non equidistant knots, binary search
0191 while((khig-klow)>1) {
0192 khalf = (klow+khig)/2;
0193 if(x>fPoly[khalf].X()) klow=khalf;
0194 else khig=khalf;
0195 }
0196 //
0197 // This could be removed, sanity check
0198 if( (x<fPoly[klow].X()) || (fPoly[klow+1].X()<x) ) {
0199 m_out << "tools::spline::cubic::find_x : Binary search failed"
0200 << " x(" << klow << ") = " << fPoly[klow].X() << " < x= " << x
0201 << " < x(" << klow+1 << ") = " << fPoly[klow+1].X() << "."
0202 << "." << std::endl;
0203 }
0204 }
0205 }
0206 return klow;
0207 }
0208
0209 void build_coeff() {
0210 /// subroutine cubspl ( tau, c, n, ibcbeg, ibcend )
0211 /// from * a practical guide to splines * by c. de boor
0212 /// ************************ input ***************************
0213 /// n = number of data points. assumed to be .ge. 2.
0214 /// (tau(i), c(1,i), i=1,...,n) = abscissae and ordinates of the
0215 /// data points. tau is assumed to be strictly increasing.
0216 /// ibcbeg, ibcend = boundary condition indicators, and
0217 /// c(2,1), c(2,n) = boundary condition information. specifically,
0218 /// ibcbeg = 0 means no boundary condition at tau(1) is given.
0219 /// in this case, the not-a-knot condition is used, i.e. the
0220 /// jump in the third derivative across tau(2) is forced to
0221 /// zero, thus the first and the second cubic polynomial pieces
0222 /// are made to coincide.)
0223 /// ibcbeg = 1 means that the slope at tau(1) is made to equal
0224 /// c(2,1), supplied by input.
0225 /// ibcbeg = 2 means that the second derivative at tau(1) is
0226 /// made to equal c(2,1), supplied by input.
0227 /// ibcend = 0, 1, or 2 has analogous meaning concerning the
0228 /// boundary condition at tau(n), with the additional infor-
0229 /// mation taken from c(2,n).
0230 /// *********************** output **************************
0231 /// c(j,i), j=1,...,4; i=1,...,l (= n-1) = the polynomial coefficients
0232 /// of the cubic interpolating spline with interior knots (or
0233 /// joints) tau(2), ..., tau(n-1). precisely, in the interval
0234 /// (tau(i), tau(i+1)), the spline f is given by
0235 /// f(x) = c(1,i)+h*(c(2,i)+h*(c(3,i)+h*c(4,i)/3.)/2.)
0236 /// where h = x - tau(i). the function program *ppvalu* may be
0237 /// used to evaluate f or its derivatives from tau,c, l = n-1,
0238 /// and k=4.
0239
0240 int j, l;
0241 double divdf1,divdf3,dtau,g=0;
0242 // ***** a tridiagonal linear system for the unknown slopes s(i) of
0243 // f at tau(i), i=1,...,n, is generated and then solved by gauss elim-
0244 // ination, with s(i) ending up in c(2,i), all i.
0245 // c(3,.) and c(4,.) are used initially for temporary storage.
0246 l = int(fNp-1);
0247 // compute first differences of x sequence and store in C also,
0248 // compute first divided difference of data and store in D.
0249 {for (size_t m=1; m<fNp ; ++m) {
0250 fPoly[m].C() = fPoly[m].X() - fPoly[m-1].X();
0251 fPoly[m].D() = (fPoly[m].Y() - fPoly[m-1].Y())/fPoly[m].C();
0252 }}
0253 // construct first equation from the boundary condition, of the form
0254 // D[0]*s[0] + C[0]*s[1] = B[0]
0255 if(fBegCond==0) {
0256 if(fNp == 2) {
0257 // no condition at left end and n = 2.
0258 fPoly[0].D() = 1.;
0259 fPoly[0].C() = 1.;
0260 fPoly[0].B() = 2.*fPoly[1].D();
0261 } else {
0262 // not-a-knot condition at left end and n .gt. 2.
0263 fPoly[0].D() = fPoly[2].C();
0264 fPoly[0].C() = fPoly[1].C() + fPoly[2].C();
0265 fPoly[0].B() = ((fPoly[1].C()+2.*fPoly[0].C())*fPoly[1].D()*fPoly[2].C()+
0266 fPoly[1].C()*fPoly[1].C()*fPoly[2].D())/fPoly[0].C();
0267 }
0268 } else if (fBegCond==1) {
0269 // slope prescribed at left end.
0270 fPoly[0].B() = fValBeg;
0271 fPoly[0].D() = 1.;
0272 fPoly[0].C() = 0.;
0273 } else if (fBegCond==2) {
0274 // second derivative prescribed at left end.
0275 fPoly[0].D() = 2.;
0276 fPoly[0].C() = 1.;
0277 fPoly[0].B() = 3.*fPoly[1].D() - fPoly[1].C()/2.*fValBeg;
0278 }
0279 bool forward_gauss_elimination = true;
0280 if(fNp > 2) {
0281 // if there are interior knots, generate the corresp. equations and car-
0282 // ry out the forward pass of gauss elimination, after which the m-th
0283 // equation reads D[m]*s[m] + C[m]*s[m+1] = B[m].
0284 {for (int m=1; m<l; ++m) {
0285 g = -fPoly[m+1].C()/fPoly[m-1].D();
0286 fPoly[m].B() = g*fPoly[m-1].B() + 3.*(fPoly[m].C()*fPoly[m+1].D()+fPoly[m+1].C()*fPoly[m].D());
0287 fPoly[m].D() = g*fPoly[m-1].C() + 2.*(fPoly[m].C() + fPoly[m+1].C());
0288 }}
0289 // construct last equation from the second boundary condition, of the form
0290 // (-g*D[n-2])*s[n-2] + D[n-1]*s[n-1] = B[n-1]
0291 // if slope is prescribed at right end, one can go directly to back-
0292 // substitution, since c array happens to be set up just right for it
0293 // at this point.
0294 if(fEndCond == 0) {
0295 if (fNp > 3 || fBegCond != 0) {
0296 // not-a-knot and n .ge. 3, and either n.gt.3 or also not-a-knot at
0297 // left end point.
0298 g = fPoly[fNp-2].C() + fPoly[fNp-1].C();
0299 fPoly[fNp-1].B() = ((fPoly[fNp-1].C()+2.*g)*fPoly[fNp-1].D()*fPoly[fNp-2].C()
0300 + fPoly[fNp-1].C()*fPoly[fNp-1].C()*(fPoly[fNp-2].Y()-fPoly[fNp-3].Y())/fPoly[fNp-2].C())/g;
0301 g = -g/fPoly[fNp-2].D();
0302 fPoly[fNp-1].D() = fPoly[fNp-2].C();
0303 } else {
0304 // either (n=3 and not-a-knot also at left) or (n=2 and not not-a-
0305 // knot at left end point).
0306 fPoly[fNp-1].B() = 2.*fPoly[fNp-1].D();
0307 fPoly[fNp-1].D() = 1.;
0308 g = -1./fPoly[fNp-2].D();
0309 }
0310 } else if (fEndCond == 1) {
0311 fPoly[fNp-1].B() = fValEnd;
0312 forward_gauss_elimination = false;
0313 } else if (fEndCond == 2) {
0314 // second derivative prescribed at right endpoint.
0315 fPoly[fNp-1].B() = 3.*fPoly[fNp-1].D() + fPoly[fNp-1].C()/2.*fValEnd;
0316 fPoly[fNp-1].D() = 2.;
0317 g = -1./fPoly[fNp-2].D();
0318 }
0319 } else {
0320 if(fEndCond == 0) {
0321 if (fBegCond > 0) {
0322 // either (n=3 and not-a-knot also at left) or (n=2 and not not-a-
0323 // knot at left end point).
0324 fPoly[fNp-1].B() = 2.*fPoly[fNp-1].D();
0325 fPoly[fNp-1].D() = 1.;
0326 g = -1./fPoly[fNp-2].D();
0327 } else {
0328 // not-a-knot at right endpoint and at left endpoint and n = 2.
0329 fPoly[fNp-1].B() = fPoly[fNp-1].D();
0330 forward_gauss_elimination = false;
0331 }
0332 } else if(fEndCond == 1) {
0333 fPoly[fNp-1].B() = fValEnd;
0334 forward_gauss_elimination = false;
0335 } else if(fEndCond == 2) {
0336 // second derivative prescribed at right endpoint.
0337 fPoly[fNp-1].B() = 3.*fPoly[fNp-1].D() + fPoly[fNp-1].C()/2.*fValEnd;
0338 fPoly[fNp-1].D() = 2.;
0339 g = -1./fPoly[fNp-2].D();
0340 }
0341 }
0342 // complete forward pass of gauss elimination.
0343 if(forward_gauss_elimination) {
0344 fPoly[fNp-1].D() = g*fPoly[fNp-2].C() + fPoly[fNp-1].D();
0345 fPoly[fNp-1].B() = (g*fPoly[fNp-2].B() + fPoly[fNp-1].B())/fPoly[fNp-1].D();
0346 }
0347 // carry out back substitution
0348 j = l-1;
0349 do {
0350 fPoly[j].B() = (fPoly[j].B() - fPoly[j].C()*fPoly[j+1].B())/fPoly[j].D();
0351 --j;
0352 } while (j>=0);
0353 // ****** generate cubic coefficients in each interval, i.e., the deriv.s
0354 // at its left endpoint, from value and slope at its endpoints.
0355 for (size_t i=1; i<fNp; ++i) {
0356 dtau = fPoly[i].C();
0357 divdf1 = (fPoly[i].Y() - fPoly[i-1].Y())/dtau;
0358 divdf3 = fPoly[i-1].B() + fPoly[i].B() - 2.*divdf1;
0359 fPoly[i-1].C() = (divdf1 - fPoly[i-1].B() - divdf3)/dtau;
0360 fPoly[i-1].D() = (divdf3/dtau)/dtau;
0361 }
0362 }
0363
0364 protected:
0365 std::vector<cubic_poly> fPoly; //[fNp] Array of polynomial terms
0366 double fValBeg; // Initial value of first or second derivative
0367 double fValEnd; // End value of first or second derivative
0368 int fBegCond; // 0=no beg cond, 1=first derivative, 2=second derivative
0369 int fEndCond; // 0=no end cond, 1=first derivative, 2=second derivative
0370 };
0371
0372 }}
0373
0374 #endif
