/src/quantlib/ql/methods/finitedifferences/solvers/fdmndimsolver.hpp

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

/*
 Copyright (C) 2011 Klaus Spanderen

 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
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <https://www.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 fdmndimsolver.hpp
*/

#ifndef quantlib_fdm_n_dim_solver_hpp
#define quantlib_fdm_n_dim_solver_hpp

#include <ql/patterns/lazyobject.hpp>
#include <ql/math/interpolations/multicubicspline.hpp>
#include <ql/methods/finitedifferences/finitedifferencemodel.hpp>
#include <ql/methods/finitedifferences/meshers/fdmmesher.hpp>
#include <ql/methods/finitedifferences/solvers/fdmsolverdesc.hpp>
#include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
#include <ql/methods/finitedifferences/solvers/fdmbackwardsolver.hpp>
#include <ql/methods/finitedifferences/stepconditions/fdmsnapshotcondition.hpp>
#include <ql/methods/finitedifferences/utilities/fdminnervaluecalculator.hpp>
#include <ql/methods/finitedifferences/stepconditions/fdmstepconditioncomposite.hpp>

#include <numeric>

namespace QuantLib {

    template <Size N>
    class FdmNdimSolver : public LazyObject {
      public:
        FdmNdimSolver(const FdmSolverDesc& solverDesc,
                      const FdmSchemeDesc& schemeDesc,
                      ext::shared_ptr<FdmLinearOpComposite> op);

        void performCalculations() const override;

        Real interpolateAt(const std::vector<Real>& x) const;
        Real thetaAt(const std::vector<Real>& x) const;

        // template meta programming
        typedef typename MultiCubicSpline<N>::data_table data_table;
        void static setValue(data_table& f,
                             const std::vector<Size>& x, Real value);

      private:
        const FdmSolverDesc solverDesc_;
        const FdmSchemeDesc schemeDesc_;
        const ext::shared_ptr<FdmLinearOpComposite> op_;

        const ext::shared_ptr<FdmSnapshotCondition> thetaCondition_;
        const ext::shared_ptr<FdmStepConditionComposite> conditions_;

        std::vector<std::vector<Real> > x_;
        std::vector<Real> initialValues_;
        const std::vector<bool> extrapolation_;

        mutable ext::shared_ptr<data_table> f_;
        mutable ext::shared_ptr<MultiCubicSpline<N> > interp_;
    };


    template <Size N>
    inline FdmNdimSolver<N>::FdmNdimSolver(const FdmSolverDesc& solverDesc,
                                           const FdmSchemeDesc& schemeDesc,
                                           ext::shared_ptr<FdmLinearOpComposite> op)
    : solverDesc_(solverDesc), schemeDesc_(schemeDesc), op_(std::move(op)),
      thetaCondition_(new FdmSnapshotCondition(
          0.99 * std::min(1.0 / 365.0,
                          solverDesc.condition->stoppingTimes().empty() ?
                              solverDesc.maturity :
                              solverDesc.condition->stoppingTimes().front()))),
      conditions_(FdmStepConditionComposite::joinConditions(thetaCondition_, solverDesc.condition)),
      x_(solverDesc.mesher->layout()->dim().size()),
      initialValues_(solverDesc.mesher->layout()->size()),
      extrapolation_(std::vector<bool>(N, false)) {

        QL_REQUIRE(solverDesc.mesher->layout()->dim().size() == N, "solver dim " << N
                    << "does not fit to layout dim " << solverDesc.mesher->layout()->size());

        for (Size i=0; i < N; ++i) {
            x_[i].reserve(solverDesc.mesher->layout()->dim()[i]);
        }

        for (const auto& iter : *solverDesc.mesher->layout()) {
            initialValues_[iter.index()] = solverDesc_.calculator
                                ->avgInnerValue(iter, solverDesc.maturity);

            const std::vector<Size>& c = iter.coordinates();
            for (Size i=0; i < N; ++i) {
                if ((std::accumulate(c.begin(), c.end(), 0UL) - c[i]) == 0U) {
                    x_[i].push_back(solverDesc.mesher->location(iter, i));
                }
            }
        }

        f_ = ext::shared_ptr<data_table>(new data_table(x_));
    }


    template <Size N> inline
    void FdmNdimSolver<N>::performCalculations() const {
        Array rhs(initialValues_.size());
        std::copy(initialValues_.begin(), initialValues_.end(), rhs.begin());

        FdmBackwardSolver(op_, solverDesc_.bcSet, conditions_, schemeDesc_)
                 .rollback(rhs, solverDesc_.maturity, 0.0,
                           solverDesc_.timeSteps, solverDesc_.dampingSteps);

        for (const auto& iter : *solverDesc_.mesher->layout()) {
            setValue(*f_, iter.coordinates(), rhs[iter.index()]);
        }

        interp_ = ext::shared_ptr<MultiCubicSpline<N> >(
            new MultiCubicSpline<N>(x_, *f_, extrapolation_));
    }


    template <Size N> inline
    Real FdmNdimSolver<N>::thetaAt(const std::vector<Real>& x) const {
        if (conditions_->stoppingTimes().front() == 0.0)
            return Null<Real>();

        calculate();
        const Array& rhs = thetaCondition_->getValues();

        data_table f(x_);

        for (const auto& iter : *solverDesc_.mesher->layout()) {
            setValue(f, iter.coordinates(), rhs[iter.index()]);
        }

        return (MultiCubicSpline<N>(x_, f)(x)
                        - interpolateAt(x)) / thetaCondition_->getTime();
    }

    template <Size N> inline
    Real FdmNdimSolver<N>::interpolateAt(const std::vector<Real>& x) const {
        calculate();

        return (*interp_)(x);
    }

    template <Size N> inline
    void FdmNdimSolver<N>::setValue(data_table& f,
                                    const std::vector<Size>& x, Real value) {
        FdmNdimSolver<N-1>::setValue(f[x[x.size()-N]], x, value);
    }

    template <> inline
    void FdmNdimSolver<1>::setValue(data_table& f,
                                    const std::vector<Size>& x, Real value) {
        f[x.back()] = value;
    }
}

#endif

Line	Count	Source
1		/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3		/*
4		Copyright (C) 2011 Klaus Spanderen
5
6		This file is part of QuantLib, a free-software/open-source library
7		for financial quantitative analysts and developers - http://quantlib.org/
8
9		QuantLib is free software: you can redistribute it and/or modify it
10		under the terms of the QuantLib license. You should have received a
11		copy of the license along with this program; if not, please email
12		<quantlib-dev@lists.sf.net>. The license is also available online at
13		<https://www.quantlib.org/license.shtml>.
14
15		This program is distributed in the hope that it will be useful, but WITHOUT
16		ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17		FOR A PARTICULAR PURPOSE. See the license for more details.
18		*/
19
20		/*! \file fdmndimsolver.hpp
21		*/
22
23		#ifndef quantlib_fdm_n_dim_solver_hpp
24		#define quantlib_fdm_n_dim_solver_hpp
25
26		#include <ql/patterns/lazyobject.hpp>
27		#include <ql/math/interpolations/multicubicspline.hpp>
28		#include <ql/methods/finitedifferences/finitedifferencemodel.hpp>
29		#include <ql/methods/finitedifferences/meshers/fdmmesher.hpp>
30		#include <ql/methods/finitedifferences/solvers/fdmsolverdesc.hpp>
31		#include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
32		#include <ql/methods/finitedifferences/solvers/fdmbackwardsolver.hpp>
33		#include <ql/methods/finitedifferences/stepconditions/fdmsnapshotcondition.hpp>
34		#include <ql/methods/finitedifferences/utilities/fdminnervaluecalculator.hpp>
35		#include <ql/methods/finitedifferences/stepconditions/fdmstepconditioncomposite.hpp>
36
37		#include <numeric>
38
39		namespace QuantLib {
40
41		template <Size N>
42		class FdmNdimSolver : public LazyObject {
43		public:
44		FdmNdimSolver(const FdmSolverDesc& solverDesc,
45		const FdmSchemeDesc& schemeDesc,
46		ext::shared_ptr<FdmLinearOpComposite> op);
47
48		void performCalculations() const override;
49
50		Real interpolateAt(const std::vector<Real>& x) const;
51		Real thetaAt(const std::vector<Real>& x) const;
52
53		// template meta programming
54		typedef typename MultiCubicSpline<N>::data_table data_table;
55		void static setValue(data_table& f,
56		const std::vector<Size>& x, Real value);
57
58		private:
59		const FdmSolverDesc solverDesc_;
60		const FdmSchemeDesc schemeDesc_;
61		const ext::shared_ptr<FdmLinearOpComposite> op_;
62
63		const ext::shared_ptr<FdmSnapshotCondition> thetaCondition_;
64		const ext::shared_ptr<FdmStepConditionComposite> conditions_;
65
66		std::vector<std::vector<Real> > x_;
67		std::vector<Real> initialValues_;
68		const std::vector<bool> extrapolation_;
69
70		mutable ext::shared_ptr<data_table> f_;
71		mutable ext::shared_ptr<MultiCubicSpline<N> > interp_;
72		};
73
74
75		template <Size N>
76		inline FdmNdimSolver<N>::FdmNdimSolver(const FdmSolverDesc& solverDesc,
77		const FdmSchemeDesc& schemeDesc,
78		ext::shared_ptr<FdmLinearOpComposite> op)
79	0	: solverDesc_(solverDesc), schemeDesc_(schemeDesc), op_(std::move(op)),
80	0	thetaCondition_(new FdmSnapshotCondition(
81	0	0.99 * std::min(1.0 / 365.0,
82	0	solverDesc.condition->stoppingTimes().empty() ?
83	0	solverDesc.maturity :
84	0	solverDesc.condition->stoppingTimes().front()))),
85	0	conditions_(FdmStepConditionComposite::joinConditions(thetaCondition_, solverDesc.condition)),
86	0	x_(solverDesc.mesher->layout()->dim().size()),
87	0	initialValues_(solverDesc.mesher->layout()->size()),
88	0	extrapolation_(std::vector<bool>(N, false)) {
89
90	0	QL_REQUIRE(solverDesc.mesher->layout()->dim().size() == N, "solver dim " << N
91	0	<< "does not fit to layout dim " << solverDesc.mesher->layout()->size());
92
93	0	for (Size i=0; i < N; ++i) {
94	0	x_[i].reserve(solverDesc.mesher->layout()->dim()[i]);
95	0	}
96
97	0	for (const auto& iter : *solverDesc.mesher->layout()) {
98	0	initialValues_[iter.index()] = solverDesc_.calculator
99	0	->avgInnerValue(iter, solverDesc.maturity);
100
101	0	const std::vector<Size>& c = iter.coordinates();
102	0	for (Size i=0; i < N; ++i) {
103	0	if ((std::accumulate(c.begin(), c.end(), 0UL) - c[i]) == 0U) {
104	0	x_[i].push_back(solverDesc.mesher->location(iter, i));
105	0	}
106	0	}
107	0	}
108
109	0	f_ = ext::shared_ptr<data_table>(new data_table(x_));
110	0	} Unexecuted instantiation: QuantLib::FdmNdimSolver<3ul>::FdmNdimSolver(QuantLib::FdmSolverDesc const&, QuantLib::FdmSchemeDesc const&, boost::shared_ptr<QuantLib::FdmLinearOpComposite>) Unexecuted instantiation: QuantLib::FdmNdimSolver<4ul>::FdmNdimSolver(QuantLib::FdmSolverDesc const&, QuantLib::FdmSchemeDesc const&, boost::shared_ptr<QuantLib::FdmLinearOpComposite>) Unexecuted instantiation: QuantLib::FdmNdimSolver<1ul>::FdmNdimSolver(QuantLib::FdmSolverDesc const&, QuantLib::FdmSchemeDesc const&, boost::shared_ptr<QuantLib::FdmLinearOpComposite>) Unexecuted instantiation: QuantLib::FdmNdimSolver<2ul>::FdmNdimSolver(QuantLib::FdmSolverDesc const&, QuantLib::FdmSchemeDesc const&, boost::shared_ptr<QuantLib::FdmLinearOpComposite>)
111
112
113		template <Size N> inline
114	0	void FdmNdimSolver<N>::performCalculations() const {
115	0	Array rhs(initialValues_.size());
116	0	std::copy(initialValues_.begin(), initialValues_.end(), rhs.begin());
117
118	0	FdmBackwardSolver(op_, solverDesc_.bcSet, conditions_, schemeDesc_)
119	0	.rollback(rhs, solverDesc_.maturity, 0.0,
120	0	solverDesc_.timeSteps, solverDesc_.dampingSteps);
121
122	0	for (const auto& iter : *solverDesc_.mesher->layout()) {
123	0	setValue(*f_, iter.coordinates(), rhs[iter.index()]);
124	0	}
125
126	0	interp_ = ext::shared_ptr<MultiCubicSpline<N> >(
127	0	new MultiCubicSpline<N>(x_, *f_, extrapolation_));
128	0	} Unexecuted instantiation: QuantLib::FdmNdimSolver<3ul>::performCalculations() const Unexecuted instantiation: QuantLib::FdmNdimSolver<4ul>::performCalculations() const Unexecuted instantiation: QuantLib::FdmNdimSolver<1ul>::performCalculations() const Unexecuted instantiation: QuantLib::FdmNdimSolver<2ul>::performCalculations() const
129
130
131		template <Size N> inline
132		Real FdmNdimSolver<N>::thetaAt(const std::vector<Real>& x) const {
133		if (conditions_->stoppingTimes().front() == 0.0)
134		return Null<Real>();
135
136		calculate();
137		const Array& rhs = thetaCondition_->getValues();
138
139		data_table f(x_);
140
141		for (const auto& iter : *solverDesc_.mesher->layout()) {
142		setValue(f, iter.coordinates(), rhs[iter.index()]);
143		}
144
145		return (MultiCubicSpline<N>(x_, f)(x)
146		- interpolateAt(x)) / thetaCondition_->getTime();
147		}
148
149		template <Size N> inline
150	0	Real FdmNdimSolver<N>::interpolateAt(const std::vector<Real>& x) const {
151	0	calculate();
152
153	0	return (*interp_)(x);
154	0	} Unexecuted instantiation: QuantLib::FdmNdimSolver<3ul>::interpolateAt(std::__1::vector<double, std::__1::allocator<double> > const&) const Unexecuted instantiation: QuantLib::FdmNdimSolver<4ul>::interpolateAt(std::__1::vector<double, std::__1::allocator<double> > const&) const Unexecuted instantiation: QuantLib::FdmNdimSolver<1ul>::interpolateAt(std::__1::vector<double, std::__1::allocator<double> > const&) const Unexecuted instantiation: QuantLib::FdmNdimSolver<2ul>::interpolateAt(std::__1::vector<double, std::__1::allocator<double> > const&) const
155
156		template <Size N> inline
157		void FdmNdimSolver<N>::setValue(data_table& f,
158	0	const std::vector<Size>& x, Real value) {
159	0	FdmNdimSolver<N-1>::setValue(f[x[x.size()-N]], x, value);
160	0	} Unexecuted instantiation: QuantLib::FdmNdimSolver<3ul>::setValue(QuantLib::detail::DataTable<QuantLib::detail::DataTable<QuantLib::detail::DataTable<double> > >&, std::__1::vector<unsigned long, std::__1::allocator<unsigned long> > const&, double) Unexecuted instantiation: QuantLib::FdmNdimSolver<2ul>::setValue(QuantLib::detail::DataTable<QuantLib::detail::DataTable<double> >&, std::__1::vector<unsigned long, std::__1::allocator<unsigned long> > const&, double) Unexecuted instantiation: QuantLib::FdmNdimSolver<4ul>::setValue(QuantLib::detail::DataTable<QuantLib::detail::DataTable<QuantLib::detail::DataTable<QuantLib::detail::DataTable<double> > > >&, std::__1::vector<unsigned long, std::__1::allocator<unsigned long> > const&, double)
161
162		template <> inline
163		void FdmNdimSolver<1>::setValue(data_table& f,
164	0	const std::vector<Size>& x, Real value) {
165	0	f[x.back()] = value;
166	0	}
167		}
168
169		#endif

Coverage Report

Created: 2026-06-08 06:47