multicubicspline.hpp

有很多的函数库

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2003, 2004 Roman Gitlin

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <[e-mail]>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

/*! \file multicubicspline.hpp
    \brief N-dimensional cubic spline interpolation between discrete points
*/

#ifndef quantlib_multi_cubic_spline_hpp
#define quantlib_multi_cubic_spline_hpp

#include <ql/errors.hpp>
#include <ql/types.hpp>
#include <functional>
#include <vector>

namespace QuantLib {

    namespace detail {

        // data structures

        typedef std::vector<std::vector<Real> > SplineGrid;

        // different termination markers are necessary
        // to maintain separation of possibly overlapping types
        struct EmptyArg {};  // arg_type termination marker
        struct EmptyRes {};  // res_type termination marker
        struct EmptyDim {};  // size_t termination marker

        template<class X> struct DataTable {
            DataTable<X>(const std::vector<Size>::const_iterator &i) {
                std::vector<X> temp(*i, X(i + 1));
                data_table_.swap(temp);
            }
            DataTable<X>(const SplineGrid::const_iterator &i) {
                std::vector<X> temp(i->size(), X(i + 1));
                data_table_.swap(temp);
            }
            template<class U> DataTable<X>(const std::vector<U> &v) {
                DataTable temp(v.begin());
                data_table_.swap(temp.data_table_);
            }
            Size size() const {
                return data_table_.size();
            }
            const X &operator[](Size n) const {return data_table_[n];}
            X &operator[](Size n) {return data_table_[n];}
            std::vector<X> data_table_;
        };

        template<> struct DataTable<Real> {
            DataTable<Real>(Size n) : data_table_(n) {}
            DataTable<Real>(const std::vector<Size>::const_iterator& i)
            : data_table_(*i) {}
            DataTable<Real>(const SplineGrid::const_iterator &i)
            : data_table_(i->size()) {}
            template<class U> DataTable<Real>(const std::vector<U> &v) {
                DataTable temp(v.begin());
                data_table_.swap(temp.data_table_);
            }
            Size size() const {
                return data_table_.size();
            }
            Real operator[](Size n) const {return data_table_[n];}
            Real& operator[](Size n) {return data_table_[n];}
            std::vector<Real> data_table_;
        };

        typedef DataTable<Real> base_data_table;

        template<class X, class Y> struct Data {
            Data<X, Y>()
            : first(X()), second(Y()) {}
            Data<X, Y>(const SplineGrid::const_iterator &i)
            : first(*i), second(i + 1) {}
            Data<X, Y>(const SplineGrid &v)
            : first(v[0]), second(v.begin()+1) {}
            void swap(Data<X, Y> &d) {
                first.swap(d.first);
                second.swap(d.second);
            }
            X first;
            Y second;
        };

        template<> struct Data<std::vector<Real>, EmptyArg> {
            Data<std::vector<Real>, EmptyArg>()
            : first(std::vector<Real>()) {}
            Data<std::vector<Real>, EmptyArg>(const SplineGrid::const_iterator &i)
            : first(*i) {}
            Data<std::vector<Real>, EmptyArg>(const SplineGrid &v)
            : first(v[0]) {}
            Data<std::vector<Real>, EmptyArg>(const std::vector<Real> &v)
            : first(v) {}
            void swap(Data<std::vector<Real>, EmptyArg> &d) {
                first.swap(d.first);
            }
            Real operator[](Size n) const {return first[n];}
            Real& operator[](Size n) {return first[n];}
            std::vector<Real> first;
            EmptyArg second;
        };

        typedef Data<std::vector<Real>, EmptyArg> base_data;

        template<class X, class Y> struct Point {
            typedef X data_type;
            Point<X, Y>()
            : first(data_type()), second(Y()) {}
            Point<X, Y>(const std::vector<Real>::const_iterator &i)
            : first(*i), second(i + 1) {}
            Point<X, Y>(const std::vector<Real> &v)
            : first(v[0]), second(v.begin()+1) {}
            Point<X, Y>(const SplineGrid::const_iterator &i)
            : first(i->size()), second(i + 1) {}
            Point<X, Y>(const SplineGrid &grid)
            : first(grid[0].size()), second(grid.begin()+1) {}
            operator data_type() const {
                return first;
            }
            data_type operator[](Size n) const {
                return n ? second[n - 1] : first;
            }
            data_type& operator[](Size n) {
                return n ? second[n - 1] : first;
            }
            data_type first;
            Y second;
        };

        template<> struct Point<Real, EmptyArg> {
            typedef Real data_type;
            Point<Real, EmptyArg>(data_type s)
            : first(s) {}
            Point<Real, EmptyArg>(const std::vector<Real>::const_iterator &i)
            : first(*i) {}
            Point<Real, EmptyArg>(const std::vector<Real> &v)
            : first(v[0]) {}
            operator data_type() const {return first;}
            data_type operator[](Size n) const {
                QL_REQUIRE(n == 0, "operator[] : access violation");
                return first;
            }
            data_type& operator[](Size n) {
                QL_REQUIRE(n == 0, "operator[] : access violation");
                return first;
            }
            data_type first;
            EmptyArg second;
        };

        typedef Point<Real, EmptyArg> base_arg_type;

        template<> struct Point<Real, EmptyRes> {
            typedef Real data_type;
            Point<Real, EmptyRes>()
            : first(data_type()) {}
            Point<Real, EmptyRes>(data_type s)
            : first(s) {}
            operator data_type() const {return first;}
            const data_type &operator[](Size n) const {
                QL_REQUIRE(n == 0, "operator[] : access violation");
                return first;
            }
            data_type &operator[](Size n) {
                QL_REQUIRE(n == 0, "operator[] : access violation");
                return first;
            }
            data_type first;
            EmptyRes second;
        };

        typedef Point<Real, EmptyRes> base_return_type;

        template<> struct Point<Size, EmptyDim> {
            typedef Size data_type;
            Point<Size, EmptyDim>()
            : first(data_type()) {}
            Point<Size, EmptyDim>(data_type s)
            : first(s) {}
            operator data_type() const {return first;}
            data_type operator[](Size n) const {
                QL_REQUIRE(n == 0, "operator[] : access violation");
                return first;
            }
            data_type& operator[](Size n) {
                QL_REQUIRE(n == 0, "operator[] : access violation");
                return first;
            }
            data_type first;
            EmptyDim second;
        };

        typedef Point<Size, EmptyDim> base_dimensions;

        template<> struct Point<base_data_table, EmptyRes> {
            typedef base_data_table data_type;
            Point<base_data_table, EmptyRes>(data_type s)
            : first(s) {}
            Point<base_data_table, EmptyRes>(const SplineGrid::const_iterator &i)
            : first(i->size()) {}
            Point<base_data_table, EmptyRes>(const SplineGrid &grid)
            : first(grid[0].size()) {}
            Real operator[](Size n) const {return first[n];}
            Real& operator[](Size n) {return first[n];}
            data_type first;
            EmptyRes second;
        };

        typedef Point<base_data_table, EmptyRes> base_output_data;


        // cubic spline iplementations

        // no heap memory is allocated
        // in any of the recursive calls
        class base_cubic_spline : public std::unary_function<Real,Real> {
          public:
            typedef base_data data;
            typedef base_data_table data_table;
            typedef base_output_data output_data;
            base_cubic_spline(const data &d, const data &d2,
                              const data_table& y, data_table &y2,
                              output_data &v) {
                Size dim = d.first.size();
                Size j = 1, k = 2, l = 3;
                result_type  &w = ((y2[dim] = y[1]) -= y[0]) /= d[0],
                    &u = ((y2[0] = y[2]) -= y[1]) /= d[1], &t = v[dim];
                y2[1] = -d[1] / d2[0], v[1] = 6.0 * (u - w) / d2[0];
                for(; k < dim; u = w, j = k, k = l, ++l) {
                    w = (y[l]-y[k])/d[k];
                    u = (u-w)*6.0;
                    (y2[k] = d[k]) /= ((t = -y2[j]) *= d[j]) -= d2[j];
                    (v[k] = (u += d[j] * v[j])) /= t;
                }
                y2[0] = y2[dim] = 0.0;
                while (k) {
                    (y2[k-1] *= y2[l-1]) += v[k-1];
                    --k; --l;
                }
            }
        };

        template<class X>
        class n_cubic_spline {
          public:
            typedef Data<base_data, typename X::data> data;
            typedef DataTable<typename X::data_table> data_table;
            typedef Point<base_output_data, typename X::output_data> output_data;
            n_cubic_spline(const data &d, const data &d2,
                           const data_table &y, data_table &y2, output_data &v)
            :  d_(d), d2_(d2), y_(y), y2_(y2), v_(v) {
                for(Size j = 0, dim = y_.size();  j < dim; ++j)
                    X(d_.second, d2_.second, y_[j], y2_[j], v_.second);
            }
            ~n_cubic_spline(){}
          private:
            const data &d_, &d2_;
            const data_table &y_;
            mutable data_table &y2_;
            mutable output_data &v_;
        };

        class base_cubic_splint : public std::unary_function<base_arg_type,Real> {
          public:
            typedef base_data data;
            typedef base_data_table data_table;
            typedef base_dimensions dimensions;