1/*!
2 * \file   include/TFEL/Math/CubicSpline.ixx
3 * \author Castelier Etienne
4 * \date   07/08/2007
5 * \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All rights
6 * reserved.
7 * This project is publicly released under either the GNU GPL Licence
8 * or the CECILL-A licence. A copy of thoses licences are delivered
9 * with the sources of TFEL. CEA or EDF may also distribute this
10 * project under specific licensing conditions.
11 */
12
13#ifndef LIB_TFEL_MATH_CUBICSPLINE_IXX
14#define LIB_TFEL_MATH_CUBICSPLINE_IXX 1
15
16#include <limits>
17#include <algorithm>
18
19#include "TFEL/Math/General/Abs.hxx"
20
21namespace tfel {
22
23  namespace math {
24
25    template <typename real, typename value>
26    bool CubicSpline<real, value>::PointComparator::operator()(
27        const typename CubicSpline<real, value>::Point& p,
28        const real& x) const {
29      return p.x < x;
30    }  // end of CubicSpline<real,value>::Point
31
32    template <typename real, typename value>
33    template <typename AIterator, typename OIterator>
34    void CubicSpline<real, value>::setCollocationPoints(AIterator px,
35                                                        AIterator pxe,
36                                                        OIterator py) {
37      using namespace std;
38      AIterator px2;
39      if (px == pxe) {
40        throw(CubicSplineInvalidAbscissaVectorSize());
41      }
42      this->values.clear();
43      Point p;
44      p.x = *px;
45      p.y = *py;
46      this->values.push_back(p);
47      px2 = px;
48      ++px;
49      ++py;
50      while (px != pxe) {
51        px2 = px - 1u;
52        if (*px2 >= *px) {
53          throw(CubicSplineUnorderedAbscissaVector());
54        }
55        p.x = *px;
56        p.y = *py;
57        this->values.push_back(p);
58        px2 = px;
59        ++px;
60        ++py;
61      }
62      this->buildInterpolation();
63    }  // CubicSpline<real,value>::CubicSpline
64
65    template <typename real, typename value>
66    template <typename AContainer, typename OContainer>
67    void CubicSpline<real, value>::setCollocationPoints(const AContainer& x,
68                                                        const OContainer& y) {
69      if (x.size() < 1) {
70        throw(CubicSplineInvalidAbscissaVectorSize());
71      }
72      if (y.size() < 1) {
73        throw(CubicSplineInvalidOrdinateVectorSize());
74      }
75      if (x.size() != y.size()) {
76        throw(CubicSplineInvalidInputs());
77      }
78      this->setCollocationPoints(x.begin(), x.end(), y.begin());
79    }  // CubicSpline<real,value>::CubicSpline
80
81    template <typename real, typename value>
82    void CubicSpline<real, value>::buildInterpolation() {
83      using namespace std;
84      using std::vector;
85      if (this->values.size() == 1) {
86        this->values[0].d = value(real(0));
87      } else {
88        typename vector<real>::size_type i;
89        typename vector<real>::size_type s = this->values.size() - 1u;
90        vector<real> m(2 * s + 1u);    // matrix (main an upper diagonal)
91        real* const md = &m[0];        // main  diagonal
92        real* const mu = md + s + 1u;  // upper diagonal
93        real ho = real(0);
94        value uo = value(real(0));
95        for (i = 0; i != s; ++i) {
96          const real hn = 1. / (this->values[i + 1].x - this->values[i].x);
97          const value un =
98              3. * hn * hn * (this->values[i + 1].y - this->values[i].y);
99          mu[i] = hn;
100          md[i] = 2. * (hn + ho);
101          this->values[i].d = un + uo;
102          uo = un;
103          ho = hn;
104        }
105        md[s] = 2. * ho;
106        this->values[s].d = uo;
107        solveTridiagonalLinearSystem(mu, md);
108      }
109    }
110
111    template <typename real, typename value>
112    void CubicSpline<real, value>::solveTridiagonalLinearSystem(
113        const real* const c, real* const b) {
114      using namespace std;
115      using std::vector;
116      typename vector<typename CubicSpline<real, value>::Point>::size_type i;
117      typename vector<typename CubicSpline<real, value>::Point>::size_type n;
118      const real prec = 100 * numeric_limits<real>::min();
119      n = this->values.size();
120      for (i = 1; i < n; i++) {
121        if (abs(b[i - 1]) < prec) {
122          throw(CubicSplineNullPivot());
123        }
124        real m = c[i - 1] / b[i - 1];
125        b[i] -= m * c[i - 1];
126        this->values[i].d -= m * this->values[i - 1].d;
127      }
128      if (abs(b[n - 1]) < prec) {
129        throw(CubicSplineNullPivot());
130      }
131      this->values[n - 1].d /= b[n - 1];
132      i = n;
133      i -= 2u;
134      for (; i != 0; i--) {
135        if (abs(b[i]) < prec) {
136          throw(CubicSplineNullPivot());
137        }
138        this->values[i].d =
139            (this->values[i].d - c[i] * (this->values[i + 1].d)) / b[i];
140      }
141      if (abs(b[0]) < prec) {
142        throw(CubicSplineNullPivot());
143      }
144      this->values[0].d =
145          (this->values[0].d - c[0] * (this->values[1].d)) / b[0];
146    }
147
148    template <typename real, typename value>
149    void CubicSpline<real, value>::computeLocalCoefficients(
150        value& a2,
151        value& a3,
152        const typename std::vector<
153            typename CubicSpline<real, value>::Point>::const_iterator pa) {
154      using namespace std;
155      const typename vector<Point>::const_iterator pb = pa + 1;
156      const real usL = 1 / (pb->x - pa->x);
157      const value Dy = (pb->y - pa->y) * usL;
158      a2 = (3 * Dy - pb->d - 2 * pa->d) * usL;
159      a3 = (-2 * Dy + pb->d + pa->d) * usL * usL;
160    }  // end of CubicSpline<real,value>::computeLocalCoefficient
161
162    template <typename real, typename value>
163    value CubicSpline<real, value>::computeLocalIntegral(
164        const real xa,
165        const real xb,
166        const typename std::vector<
167            typename CubicSpline<real, value>::Point>::const_iterator pa) {
168      using namespace std;
169      const value py = pa->y;
170      const value d = pa->d;
171      value a2;
172      value a3;
173      CubicSpline<real, value>::computeLocalCoefficients(a2, a3, pa);
174      const real x0 = xa - pa->x;
175      const real x1 = xb - pa->x;
176      return (3 * a3 * (x1 * x1 * x1 * x1 - x0 * x0 * x0 * x0) +
177              4 * a2 * (x1 * x1 * x1 - x0 * x0 * x0) +
178              6 * d * (x1 * x1 - x0 * x0) + 12 * py * (x1 - x0)) /
179             12;
180    }  // end of
181
182    template <typename real, typename value>
183    value CubicSpline<real, value>::computeMeanValue(const real xa,
184                                                     const real xb) const {
185      return this->computeIntegral(xa, xb) / (xb - xa);
186    }  // end of value CubicSpline<real,value>::computeMeanValue
187
188    template <typename real, typename value>
189    value CubicSpline<real, value>::computeIntegral(const real xa,
190                                                    const real xb) const {
191      using namespace std;
192      typename vector<Point>::const_iterator pa;
193      typename vector<Point>::const_iterator pb;
194      if (this->values.empty()) {
195        throw(CubicSplineUninitialised());
196      }
197      if (this->values.size() == 1) {
198        const auto& f = values.back().y;
199        return f * (xb - xa);
200      }
201      if (xb < xa) {
202        return -this->computeIntegral(xb, xa);
203      }
204      pa = lower_bound(this->values.begin(), this->values.end(), xa,
205                       PointComparator());
206      pb = lower_bound(this->values.begin(), this->values.end(), xb,
207                       PointComparator());
208      if (pa == pb) {
209        if (pb == this->values.begin()) {
210          const real xe = this->values.front().x;
211          const auto& ye = this->values.front().y;
212          const auto& df = this->values.front().d;
213          // l'équation de l'extrapoltion est :
214          // y = ye-df*xe + df*x
215          // l'intégrale est alors:
216          // (ye-df*xe)*(xb-xa)+0.5*df*(xb-xa)*(xb-xa)
217          return ye * (xb - xa) +
218                 0.5 * df * ((xb - xe) * (xb - xe) - (xa - xe) * (xa - xe));
219        } else if (pb == this->values.end()) {
220          const real xe = this->values.back().x;
221          const auto& ye = this->values.back().y;
222          const auto& df = this->values.back().d;
223          // l'équation de l'extrapoltion est :
224          // y = ye-df*xe + df*x
225          // l'intégrale est alors:
226          // (ye-df*xe)*(xb-xa)+0.5*df*(xb-xa)*(xb-xa)
227          return ye * (xb - xa) +
228                 0.5 * df * ((xb - xe) * (xb - xe) - (xa - xe) * (xa - xe));
229        } else {
230          pa = pb - 1;
231          // spline=pa->y+(x-pa->x)*(pa->d+(x-pa->x)*(a2+(x-pa->x)*a3));
232          //       =pa->y+(x-pa->x)*(pa->d+(x-pa->x)*(a2+(x-pa->x)*a3));
233          return CubicSpline<real, value>::computeLocalIntegral(xa, xb, pa);
234        }
235      }
236      value s(real(0));
237      if (pa == this->values.begin()) {
238        const real xe = this->values.front().x;
239        const auto& ye = this->values.front().y;
240        const auto& df = this->values.front().d;
241        // l'équation de l'extrapoltion est :
242        // y = ye-df*xe + df*x
243        // l'intégrale est alors:
244        // (ye-df*xe)*(pa->a-xa)+0.5*df*(pa->a-xa)*(pa->a-xa)
245        s += ye * (xe - xa) - 0.5 * df * ((xa - xe) * (xa - xe));
246      } else {
247        s += CubicSpline<real, value>::computeLocalIntegral(xa, pa->x, pa - 1);
248      }
249      if (pb == this->values.end()) {
250        const real xe = this->values.back().x;
251        const auto& ye = this->values.back().y;
252        const auto& df = this->values.back().d;
253        // l'équation de l'extrapoltion est :
254        // y = ye-df*xe + df*x
255        // l'intégrale est alors:
256        // (ye-df*xe)*(xb-xa)+0.5*df*(xb-xa)*(xb-xa)
257        s += ye * (xb - xe) + 0.5 * df * ((xb - xe) * (xb - xe));
258        //	s += (ye-df*xe)*(xb-pb->x)+0.5*df*(xb-pb->x)*(xb-pb->x);
259      } else {
260        s += CubicSpline<real, value>::computeLocalIntegral((pb - 1)->x, xb,
261                                                            pb - 1);
262      }
263      --pb;
264      while (pa != pb) {
265        s += CubicSpline<real, value>::computeLocalIntegral(pa->x, (pa + 1)->x,
266                                                            pa);
267        ++pa;
268      }
269      return s;
270    }  // end of CubicSpline<real,value>::computeIntegral
271
272    // Valeur de la spline et de sa dérivée
273    // Interpolation cubique
274    // x : Abscisse où calculer la valeur de la spline
275    // f : Valeur de la fonction
276    // df : Valeur de la dérivée
277    template <typename real, typename value>
278    void CubicSpline<real, value>::getValues(value& f,
279                                             value& df,
280                                             real x) const {
281      using namespace std;
282      typename vector<Point>::const_iterator in;
283      if (this->values.empty()) {
284        throw(CubicSplineUninitialised());
285      }
286      // Cas d'un seul couple
287      if (this->values.size() == 1) {
288        f = this->values[0].y;
289        df = value(real(0));
290        return;
291      }
292      in = lower_bound(this->values.begin(), this->values.end(), x,
293                       PointComparator());
294      // Extrapolation
295      if (in == this->values.begin()) {
296        x -= in->x;
297        df = in->d;
298        f = in->y + x * df;
299        return;
300      }
301      if (in == this->values.end()) {
302        --in;
303        x -= in->x;
304        df = in->d;
305        f = in->y + x * df;
306        return;
307      }
308      typename vector<Point>::const_iterator ip = in - 1;
309      value a2;
310      value a3;
311      CubicSpline<real, value>::computeLocalCoefficients(a2, a3, ip);
312      x -= ip->x;
313      f = ip->y + x * (ip->d + x * (a2 + x * a3));
314      df = ip->d + x * (2 * a2 + x * 3 * a3);
315    }
316
317    template <typename real, typename value>
318    value CubicSpline<real, value>::operator()(const real x) const {
319      return this->getValue(x);
320    }
321
322    // Valeur de la spline
323    // Interpolation cubique
324    // x : Abscisse où calculer la valeur de la spline
325    template <typename real, typename value>
326    value CubicSpline<real, value>::getValue(real x) const {
327      using namespace std;
328      typename vector<Point>::const_iterator in;
329      if (this->values.empty()) {
330        throw(CubicSplineUninitialised());
331      }
332      // Cas d'un seul couple
333      if (this->values.size() == 1) {
334        return this->values[0].y;
335      }
336      in = lower_bound(this->values.begin(), this->values.end(), x,
337                       PointComparator());
338      // Extrapolation
339      if (in == this->values.begin()) {
340        x -= in->x;
341        return in->y + x * in->d;
342      }
343      if (in == this->values.end()) {
344        --in;
345        x -= in->x;
346        return in->y + x * in->d;
347      }
348      typename vector<Point>::const_iterator ip = in - 1;
349      value a2;
350      value a3;
351      CubicSpline<real, value>::computeLocalCoefficients(a2, a3, ip);
352      x -= ip->x;
353      return ip->y + x * (ip->d + x * (a2 + x * a3));
354    }
355
356    // Valeur de la spline et de sa dérivée
357    // Interpolation cubique
358    // x : Abscisse où calculer la valeur de la spline
359    // f : Valeur de la fonction
360    // df : Valeur de la dérivée
361    // d2f : Valeur de la dérivée seconde
362    template <typename real, typename value>
363    void CubicSpline<real, value>::getValues(value& f,
364                                             value& df,
365                                             value& d2f,
366                                             real x) const {
367      using namespace std;
368      if (this->values.empty()) {
369        throw(CubicSplineUninitialised());
370      }
371      typename vector<Point>::const_iterator in;
372      in = this->values.begin();
373      // Cas d'un seul couple
374      if (this->values.size() == 1) {
375        f = this->values[0].y;
376        df = value(real(0));
377        d2f = value(real(0));
378        return;
379      }
380      in = lower_bound(this->values.begin(), this->values.end(), x,
381                       PointComparator());
382      // Extrapolation
383      if (in == this->values.begin()) {
384        x -= in->x;
385        df = in->d;
386        f = in->y + x * df;
387        d2f = value(real(0));
388        return;
389      }
390      if (in == this->values.end()) {
391        --in;
392        x -= in->x;
393        df = in->d;
394        f = in->y + x * df;
395        d2f = value(real(0));
396        return;
397      }
398      typename vector<Point>::const_iterator ip = in - 1;
399      value a2;
400      value a3;
401      CubicSpline<real, value>::computeLocalCoefficients(a2, a3, ip);
402      x -= ip->x;
403      f = ip->y + x * (ip->d + x * (a2 + x * a3));
404      df = ip->d + x * (2 * a2 + x * 3 * a3);
405      d2f = 2 * a2 + x * 6 * a3;
406    }
407
408  }  // end of namespace math
409
410}  // end of namespace tfel
411
412#endif /* LIB_TFEL_MATH_CUBICSPLINE_IXX */
413