multicubicspline.hpp

有很多的函数库

            typedef base_output_data output_data;
            typedef base_return_type return_type;
            base_cubic_splint(const return_type &a, const return_type &b,
                              const return_type &a2, const return_type &b2,
                              const dimensions &i,
                              const data&, const data&,
                              const data_table &y, data_table &y2,
                              output_data&,
                              output_data&, output_data&,
                              result_type &res) {
                res = a * y[i] + b * y[i + 1] + a2 * y2[i] + b2 * y2[i + 1];
            }
        };

        template<class X>
        class n_cubic_splint : public
        std::unary_function<Point<Real, typename X::argument_type>, Real> {
          public:
            typedef std::unary_function<Point<Real, typename X::argument_type>,
                                        Real> super;
            typedef Data<base_data, typename X::data> data;
            typedef DataTable<typename X::data_table> data_table;
            typedef Point<Size, typename X::dimensions> dimensions;
            typedef Point<base_output_data, typename X::output_data> output_data;
            typedef Point<typename super::result_type,
                          typename X::return_type> return_type;
            n_cubic_splint(const return_type &a, const return_type &b,
                           const return_type &a2, const return_type &b2,
                           const dimensions &i, const data &d, const data &d2,
                           const data_table &y, data_table &y2, output_data &v,
                           output_data &v1, output_data &v2,
                           typename super::result_type& r)
            :  a_(a), b_(b), a2_(a2), b2_(b2), i_(i), d_(d), d2_(d2),
               y_(y), y2_(y2), v_(v), v1_(v1), v2_(v2) {
                for(Size j = 0, dim = y_.size(); j < dim; ++j)
                    X(a_.second, b_.second, a2_.second, b2_.second, i_.second,
                      d_.second, d2_.second, y_[j], y2_[j], v_.second,
                      v1_.second, v2_.second, v1_.first[j]);
                base_cubic_spline(d_.first, d2_.first,
                                  v1_.first.first, v2_.first.first, v_.first);
                base_cubic_splint(a_.first, b_.first, a2_.first, b2_.first,
                                  i_.first, d_.first, d2_.first,
                                  v1_.first.first, v2_.first.first,
                                  v_.first, v1_.first, v2_.first, r);
            }
            ~n_cubic_splint(){}
          private:
            const return_type &a_, &b_, &a2_, &b2_;
            const dimensions &i_;
            const data &d_, &d2_;
            const data_table &y_;
            mutable data_table &y2_;
            mutable output_data &v_, &v1_, &v2_;
        };

        typedef base_cubic_spline               cubic_spline_01;
        typedef n_cubic_spline<cubic_spline_01> cubic_spline_02;
        typedef n_cubic_spline<cubic_spline_02> cubic_spline_03;
        typedef n_cubic_spline<cubic_spline_03> cubic_spline_04;
        typedef n_cubic_spline<cubic_spline_04> cubic_spline_05;
        typedef n_cubic_spline<cubic_spline_05> cubic_spline_06;
        typedef n_cubic_spline<cubic_spline_06> cubic_spline_07;
        typedef n_cubic_spline<cubic_spline_07> cubic_spline_08;
        typedef n_cubic_spline<cubic_spline_08> cubic_spline_09;
        typedef n_cubic_spline<cubic_spline_09> cubic_spline_10;
        typedef n_cubic_spline<cubic_spline_10> cubic_spline_11;
        typedef n_cubic_spline<cubic_spline_11> cubic_spline_12;
        typedef n_cubic_spline<cubic_spline_12> cubic_spline_13;
        typedef n_cubic_spline<cubic_spline_13> cubic_spline_14;
        typedef n_cubic_spline<cubic_spline_14> cubic_spline_15;

        typedef base_cubic_splint               cubic_splint_01;
        typedef n_cubic_splint<cubic_splint_01> cubic_splint_02;
        typedef n_cubic_splint<cubic_splint_02> cubic_splint_03;
        typedef n_cubic_splint<cubic_splint_03> cubic_splint_04;
        typedef n_cubic_splint<cubic_splint_04> cubic_splint_05;
        typedef n_cubic_splint<cubic_splint_05> cubic_splint_06;
        typedef n_cubic_splint<cubic_splint_06> cubic_splint_07;
        typedef n_cubic_splint<cubic_splint_07> cubic_splint_08;
        typedef n_cubic_splint<cubic_splint_08> cubic_splint_09;
        typedef n_cubic_splint<cubic_splint_09> cubic_splint_10;
        typedef n_cubic_splint<cubic_splint_10> cubic_splint_11;
        typedef n_cubic_splint<cubic_splint_11> cubic_splint_12;
        typedef n_cubic_splint<cubic_splint_12> cubic_splint_13;
        typedef n_cubic_splint<cubic_splint_13> cubic_splint_14;
        typedef n_cubic_splint<cubic_splint_14> cubic_splint_15;

        template<Size i> struct Int2Type {
            typedef cubic_spline_01 c_spline;
            typedef cubic_splint_01 c_splint;
        };

        template<> struct Int2Type<2> {
            typedef cubic_spline_02 c_spline;
            typedef cubic_splint_02 c_splint;
        };

        template<> struct Int2Type<3> {
            typedef cubic_spline_03 c_spline;
            typedef cubic_splint_03 c_splint;
        };

        template<> struct Int2Type<4> {
            typedef cubic_spline_04 c_spline;
            typedef cubic_splint_04 c_splint;
        };

        template<> struct Int2Type<5> {
            typedef cubic_spline_05 c_spline;
            typedef cubic_splint_05 c_splint;
        };

        template<> struct Int2Type<6> {
            typedef cubic_splint_06 c_splint;
            typedef cubic_spline_06 c_spline;
        };

        template<> struct Int2Type<7> {
            typedef cubic_spline_07 c_spline;
            typedef cubic_splint_07 c_splint;
        };

        template<> struct Int2Type<8> {
            typedef cubic_spline_08 c_spline;
            typedef cubic_splint_08 c_splint;
        };

        template<> struct Int2Type<9> {
            typedef cubic_spline_09 c_spline;
            typedef cubic_splint_09 c_splint;
        };

        template<> struct Int2Type<10> {
            typedef cubic_spline_10 c_spline;
            typedef cubic_splint_10 c_splint;
        };

        template<> struct Int2Type<11> {
            typedef cubic_splint_11 c_splint;
            typedef cubic_spline_11 c_spline;
        };

        template<> struct Int2Type<12> {
            typedef cubic_spline_12 c_spline;
            typedef cubic_splint_12 c_splint;
        };

        template<> struct Int2Type<13> {
            typedef cubic_spline_13 c_spline;
            typedef cubic_splint_13 c_splint;
        };

        template<> struct Int2Type<14> {
            typedef cubic_spline_14 c_spline;
            typedef cubic_splint_14 c_splint;
        };

        template<> struct Int2Type<15> {
            typedef cubic_spline_15 c_spline;
            typedef cubic_splint_15 c_splint;
        };

    }


    // Multi-cubic spline

    typedef detail::SplineGrid SplineGrid;

    //! N-dimensional cubic spline interpolation between discrete points
    /*! \test interpolated values are checked against the original
              function.

        \todo
        - allow extrapolation as for the other interpolations
        - investigate if and how to implement Hyman filters and
          different boundary conditions

        \bug cannot interpolate at the grid points on the boundary
             surface of the N-dimensional region
    */
    template <Size i> class MultiCubicSpline {
        typedef typename detail::Int2Type<i>::c_spline c_spline;
        typedef typename detail::Int2Type<i>::c_splint c_splint;
      public:
        typedef typename c_splint::argument_type argument_type;
        typedef typename c_splint::result_type result_type;
        typedef typename c_splint::data_table data_table;
        typedef typename c_splint::return_type return_type;
        typedef typename c_splint::output_data output_data;
        typedef typename c_splint::dimensions dimensions;
        typedef typename c_splint::data data;
        MultiCubicSpline(const SplineGrid& grid, const data_table &y,
                         const std::vector<bool>& ae =
                                             std::vector<bool>(20, false))
        : grid_(grid), y_(y), ae_(ae), v_(grid), v1_(grid),
          v2_(grid), y2_(grid) {
            set_shared_increments();
            c_spline(d_, d2_, y_, y2_, v_);
        }
        result_type operator()(const argument_type& x) const {
            set_shared_coefficients(x);
            c_splint(a_, b_, a2_, b2_, i_, d_, d2_, y_, y2_,
                     v_, v1_, v2_, res_);
            return res_;
        }
        void set_shared_increments() const;
        void set_shared_coefficients(const argument_type &x) const;
      private:
        const SplineGrid &grid_;
        const data_table &y_;
        const std::vector<bool> &ae_;
        mutable return_type a_, b_, a2_, b2_;
        mutable output_data v_, v1_, v2_;
        mutable result_type res_;
        mutable dimensions i_;
        mutable data d_, d2_;
        mutable data_table y2_;
    };

    // the data is checked and, in case of insufficient number of points,
    // exception is thrown BEFORE the main body of interpolation begins
    template <Size i>
    void MultiCubicSpline<i>::set_shared_increments() const {
        SplineGrid x(i), y(i);
        Size k = 0, dim = 0;
        for(Size j = 0; j < i; k = 0, ++j) {
            const std::vector<Real> &v = grid_[j];
            if((dim = v.size() - 1) > 2) {
                std::vector<Real> tmp1(dim);
                x[j].swap(tmp1);
                std::vector<Real> tmp2(dim - 1);
                y[j].swap(tmp2);
                for(; k < dim; ++k) {
                    if((x[j][k] = v[k + 1] - v[k]) <= 0.0) break;
                    if(k) y[j][k - 1] = 2.0 * (v[k + 1] - v[k - 1]);
                }
            }
            QL_REQUIRE(dim >= 3,
                       "Dimension " << j
                       << " : not enough points for interpolation");
            QL_REQUIRE(k >= dim,
                       "Dimension " << j << " : invalid data");
        }

        typename c_splint::data tmp1(x), tmp2(y);
        d_.swap(tmp1);
        d2_.swap(tmp2);
    }

    #ifndef __DOXYGEN__
    // the argument value is checked and, in out of boundaries case,
    // exception is thrown BEFORE the main body of interpolation begins
    template <Size i>
    void MultiCubicSpline<i>::set_shared_coefficients(
                 const typename MultiCubicSpline<i>::argument_type &x) const {
        for(Size j = 0; j < i; ++j) {
            Size &k = i_[j], sz = grid_[j].size() - 1;
            const std::vector<Real> &v = grid_[j];
            if(x[j] < v[0] || x[j] >= v[sz]) {
                QL_REQUIRE(ae_[j],
                           "Dimension " << j
                           << ": extrapolation is not allowed.");
                a_[j] = 1.0, a2_[j] = b_[j] = b2_[j] = 0.0;
                k =  x[j] < v[0] ? 0 : sz;
            }
            else {
                k = v[k] <= x[j] && x[j] < v[k + 1] ? k :
                    std::upper_bound(v.begin(),v.end(),x[j])-v.begin()-1;
                Real h = v[k + 1] - v[k];
                a_[j] = (v[k + 1] - x[j]) / h, b_[j] = (x[j] - v[k]) / h;
                a2_[j] = (a_[j] * a_[j] * a_[j] - a_[j]) * h * h / 6.0,
                    b2_[j] = (b_[j] * b_[j] * b_[j] - b_[j]) * h * h / 6.0;
            }
        }
    }
    #endif

}


#endif