diff --git a/QuantLib.vcxproj b/QuantLib.vcxproj index a4efc7bcc40..635e70a9716 100644 --- a/QuantLib.vcxproj +++ b/QuantLib.vcxproj @@ -2237,6 +2237,7 @@ + @@ -2271,6 +2272,7 @@ + diff --git a/QuantLib.vcxproj.filters b/QuantLib.vcxproj.filters index 6f8cacb3663..c2b7264935d 100644 --- a/QuantLib.vcxproj.filters +++ b/QuantLib.vcxproj.filters @@ -4841,6 +4841,9 @@ math + + math + math @@ -7162,6 +7165,9 @@ pricingengines\asian + + math\interpolations + math\interpolations diff --git a/ql/CMakeLists.txt b/ql/CMakeLists.txt index 34e7332439f..41354a4770b 100644 --- a/ql/CMakeLists.txt +++ b/ql/CMakeLists.txt @@ -340,6 +340,7 @@ set(QL_SOURCES legacy/libormarketmodels/lmlinexpvolmodel.cpp legacy/libormarketmodels/lmvolmodel.cpp math/abcdmathfunction.cpp + math/array.cpp math/bernsteinpolynomial.cpp math/beta.cpp math/bspline.cpp @@ -375,6 +376,7 @@ set(QL_SOURCES math/integrals/kronrodintegral.cpp math/integrals/segmentintegral.cpp math/interpolations/chebyshevinterpolation.cpp + math/interpolations/cubicinterpolation.cpp math/matrix.cpp math/matrixutilities/basisincompleteordered.cpp math/matrixutilities/bicgstab.cpp diff --git a/ql/cashflows/dividend.hpp b/ql/cashflows/dividend.hpp index ab6cef8bbca..e27a93c5e98 100644 --- a/ql/cashflows/dividend.hpp +++ b/ql/cashflows/dividend.hpp @@ -26,6 +26,7 @@ #define quantlib_dividend_hpp #include +#include #include #include diff --git a/ql/cashflows/simplecashflow.cpp b/ql/cashflows/simplecashflow.cpp index 4d8e902d956..8626b725bef 100644 --- a/ql/cashflows/simplecashflow.cpp +++ b/ql/cashflows/simplecashflow.cpp @@ -18,6 +18,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include namespace QuantLib { diff --git a/ql/currency.hpp b/ql/currency.hpp index 2d9296bca21..e96c452fc7d 100644 --- a/ql/currency.hpp +++ b/ql/currency.hpp @@ -26,6 +26,7 @@ #define quantlib_currency_hpp #include +#include #include #include #include diff --git a/ql/errors.cpp b/ql/errors.cpp index aa6b0081b49..d61913e0c5c 100644 --- a/ql/errors.cpp +++ b/ql/errors.cpp @@ -93,10 +93,11 @@ namespace QuantLib { Error::Error(const std::string& file, long line, const std::string& function, const std::string& message) { - message_ = ext::make_shared( - format(file, line, function, message)); + message_ = format(file, line, function, message); } - const char* Error::what() const noexcept { return message_->c_str(); } + const char* Error::what() const noexcept { + return message_.c_str(); + } } diff --git a/ql/errors.hpp b/ql/errors.hpp index 96966e4fde8..8279981bb0f 100644 --- a/ql/errors.hpp +++ b/ql/errors.hpp @@ -25,9 +25,6 @@ #ifndef quantlib_errors_hpp #define quantlib_errors_hpp -#include -#include -#include #include #include #include @@ -49,94 +46,58 @@ namespace QuantLib { const char* what() const noexcept override; private: - ext::shared_ptr message_; + std::string message_; }; } -#define QL_MULTILINE_FAILURE_BEGIN do { - /* Disable warning C4127 (conditional expression is constant) when wrapping macros with the do { ... } while(false) construct on MSVC */ -#if defined(BOOST_MSVC) - #define QL_MULTILINE_FAILURE_END \ - __pragma(warning(push)) \ - __pragma(warning(disable:4127)) \ - } while(false) \ - __pragma(warning(pop)) -#else - #define QL_MULTILINE_FAILURE_END } while(false) -#endif - - -#define QL_MULTILINE_ASSERTION_BEGIN do { -/* Disable warning C4127 (conditional expression is constant) when - wrapping macros with the do { ... } while(false) construct on MSVC +/*! \def QL_ASSERT + \brief throw an error if the given condition is not verified */ #if defined(BOOST_MSVC) - #define QL_MULTILINE_ASSERTION_END \ - __pragma(warning(push)) \ - __pragma(warning(disable:4127)) \ - } while(false) \ - __pragma(warning(pop)) + +# define QL_ASSERT(condition, message) \ + do { \ + if (!(condition)) { \ + std::ostringstream _ql_msg_stream; \ + _ql_msg_stream << message; \ + throw QuantLib::Error(__FILE__, __LINE__, BOOST_CURRENT_FUNCTION, \ + _ql_msg_stream.str()); \ + } \ + __pragma(warning(push)) __pragma(warning(disable : 4127)) \ + } while (false) __pragma(warning(pop)) + #else - #define QL_MULTILINE_ASSERTION_END } while(false) -#endif +# define QL_ASSERT(condition, message) \ + do { \ + if (!(condition)) { \ + std::ostringstream _ql_msg_stream; \ + _ql_msg_stream << message; \ + throw QuantLib::Error(__FILE__, __LINE__, BOOST_CURRENT_FUNCTION, \ + _ql_msg_stream.str()); \ + } \ + } while (false) + +#endif /*! \def QL_FAIL \brief throw an error (possibly with file and line information) */ -#define QL_FAIL(message) \ -QL_MULTILINE_FAILURE_BEGIN \ - std::ostringstream _ql_msg_stream; \ - _ql_msg_stream << message; \ - throw QuantLib::Error(__FILE__,__LINE__, \ - BOOST_CURRENT_FUNCTION,_ql_msg_stream.str()); \ -QL_MULTILINE_FAILURE_END - - -/*! \def QL_ASSERT - \brief throw an error if the given condition is not verified -*/ -#define QL_ASSERT(condition,message) \ -QL_MULTILINE_ASSERTION_BEGIN \ -if (!(condition)) { \ - std::ostringstream _ql_msg_stream; \ - _ql_msg_stream << message; \ - throw QuantLib::Error(__FILE__,__LINE__, \ - BOOST_CURRENT_FUNCTION,_ql_msg_stream.str()); \ -} \ -QL_MULTILINE_ASSERTION_END +#define QL_FAIL(message) QL_ASSERT(false, message) /*! \def QL_REQUIRE \brief throw an error if the given pre-condition is not verified */ -#define QL_REQUIRE(condition,message) \ -QL_MULTILINE_ASSERTION_BEGIN \ -if (!(condition)) { \ - std::ostringstream _ql_msg_stream; \ - _ql_msg_stream << message; \ - throw QuantLib::Error(__FILE__,__LINE__, \ - BOOST_CURRENT_FUNCTION,_ql_msg_stream.str()); \ -} \ -QL_MULTILINE_ASSERTION_END +#define QL_REQUIRE(condition, message) QL_ASSERT(condition, message) /*! \def QL_ENSURE \brief throw an error if the given post-condition is not verified */ -#define QL_ENSURE(condition,message) \ -QL_MULTILINE_ASSERTION_BEGIN \ -if (!(condition)) { \ - std::ostringstream _ql_msg_stream; \ - _ql_msg_stream << message; \ - throw QuantLib::Error(__FILE__,__LINE__, \ - BOOST_CURRENT_FUNCTION,_ql_msg_stream.str()); \ -} \ -QL_MULTILINE_ASSERTION_END - +#define QL_ENSURE(condition, message) QL_ASSERT(condition, message) #endif - diff --git a/ql/event.cpp b/ql/event.cpp index bc764b4df86..bba4710d7e4 100644 --- a/ql/event.cpp +++ b/ql/event.cpp @@ -18,6 +18,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include #include diff --git a/ql/experimental/commodities/paymentterm.cpp b/ql/experimental/commodities/paymentterm.cpp index e3ae2cdb80e..d88b77d09ba 100644 --- a/ql/experimental/commodities/paymentterm.cpp +++ b/ql/experimental/commodities/paymentterm.cpp @@ -19,6 +19,8 @@ #include +#include + namespace QuantLib { std::map > diff --git a/ql/experimental/math/isotropicrandomwalk.hpp b/ql/experimental/math/isotropicrandomwalk.hpp index 72d6fdcdb56..dbda3af1a4e 100644 --- a/ql/experimental/math/isotropicrandomwalk.hpp +++ b/ql/experimental/math/isotropicrandomwalk.hpp @@ -24,6 +24,7 @@ #ifndef quantlib_isotropic_random_walk_hpp #define quantlib_isotropic_random_walk_hpp +#include #include #include #include diff --git a/ql/experimental/math/multidimintegrator.hpp b/ql/experimental/math/multidimintegrator.hpp index 336999eb209..2515b0ba7dc 100644 --- a/ql/experimental/math/multidimintegrator.hpp +++ b/ql/experimental/math/multidimintegrator.hpp @@ -20,6 +20,7 @@ #ifndef quantlib_math_multidimintegrator_hpp #define quantlib_math_multidimintegrator_hpp +#include #include #include #include diff --git a/ql/experimental/volatility/noarbsabr.hpp b/ql/experimental/volatility/noarbsabr.hpp index e55071d084c..15761e178ae 100644 --- a/ql/experimental/volatility/noarbsabr.hpp +++ b/ql/experimental/volatility/noarbsabr.hpp @@ -51,6 +51,7 @@ #define quantlib_noarb_sabr #include +#include #include #include diff --git a/ql/handle.hpp b/ql/handle.hpp index 63934035c4d..c8a2cc14081 100644 --- a/ql/handle.hpp +++ b/ql/handle.hpp @@ -25,6 +25,7 @@ #ifndef quantlib_handle_hpp #define quantlib_handle_hpp +#include #include namespace QuantLib { diff --git a/ql/instrument.cpp b/ql/instrument.cpp index 68c747df63e..427b86f7cbc 100644 --- a/ql/instrument.cpp +++ b/ql/instrument.cpp @@ -18,6 +18,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include @@ -46,4 +47,74 @@ namespace QuantLib { QL_FAIL("Instrument::setupArguments() not implemented"); } + void Instrument::calculate() const { + if (!calculated_) { + if (isExpired()) { + setupExpired(); + calculated_ = true; + } else { + LazyObject::calculate(); + } + } + } + + void Instrument::setupExpired() const { + NPV_ = errorEstimate_ = 0.0; + valuationDate_ = Date(); + additionalResults_.clear(); + } + + void Instrument::performCalculations() const { + QL_REQUIRE(engine_, "null pricing engine"); + engine_->reset(); + setupArguments(engine_->getArguments()); + engine_->getArguments()->validate(); + engine_->calculate(); + fetchResults(engine_->getResults()); + } + + void Instrument::fetchResults( + const PricingEngine::results* r) const { + const auto* results = dynamic_cast(r); + QL_ENSURE(results != nullptr, "no results returned from pricing engine"); + + NPV_ = results->value; + errorEstimate_ = results->errorEstimate; + valuationDate_ = results->valuationDate; + + additionalResults_ = results->additionalResults; + } + + Real Instrument::NPV() const { + calculate(); + QL_REQUIRE(NPV_ != Null(), "NPV not provided"); + return NPV_; + } + + Real Instrument::errorEstimate() const { + calculate(); + QL_REQUIRE(errorEstimate_ != Null(), + "error estimate not provided"); + return errorEstimate_; + } + + const Date& Instrument::valuationDate() const { + calculate(); + QL_REQUIRE(valuationDate_ != Date(), + "valuation date not provided"); + return valuationDate_; + } + + const std::map& + Instrument::additionalResults() const { + calculate(); + return additionalResults_; + } + + void Instrument::results::reset() { + value = errorEstimate = Null(); + valuationDate = Date(); + additionalResults.clear(); + } + } diff --git a/ql/instrument.hpp b/ql/instrument.hpp index 833665bb0eb..80e88a837ca 100644 --- a/ql/instrument.hpp +++ b/ql/instrument.hpp @@ -25,6 +25,7 @@ #ifndef quantlib_instrument_hpp #define quantlib_instrument_hpp +#include #include #include #include @@ -112,11 +113,7 @@ namespace QuantLib { class Instrument::results : public virtual PricingEngine::results { public: - void reset() override { - value = errorEstimate = Null(); - valuationDate = Date(); - additionalResults.clear(); - } + void reset() override; Real value; Real errorEstimate; Date valuationDate; @@ -124,82 +121,14 @@ namespace QuantLib { }; - // inline definitions - - inline void Instrument::calculate() const { - if (!calculated_) { - if (isExpired()) { - setupExpired(); - calculated_ = true; - } else { - LazyObject::calculate(); - } - } - } - - inline void Instrument::setupExpired() const { - NPV_ = errorEstimate_ = 0.0; - valuationDate_ = Date(); - additionalResults_.clear(); - } - - inline void Instrument::performCalculations() const { - QL_REQUIRE(engine_, "null pricing engine"); - engine_->reset(); - setupArguments(engine_->getArguments()); - engine_->getArguments()->validate(); - engine_->calculate(); - fetchResults(engine_->getResults()); - } - - inline void Instrument::fetchResults( - const PricingEngine::results* r) const { - const auto* results = dynamic_cast(r); - QL_ENSURE(results != nullptr, "no results returned from pricing engine"); - - NPV_ = results->value; - errorEstimate_ = results->errorEstimate; - valuationDate_ = results->valuationDate; - - additionalResults_ = results->additionalResults; - } - - inline Real Instrument::NPV() const { - calculate(); - QL_REQUIRE(NPV_ != Null(), "NPV not provided"); - return NPV_; - } - - inline Real Instrument::errorEstimate() const { - calculate(); - QL_REQUIRE(errorEstimate_ != Null(), - "error estimate not provided"); - return errorEstimate_; - } - - inline const Date& Instrument::valuationDate() const { - calculate(); - QL_REQUIRE(valuationDate_ != Date(), - "valuation date not provided"); - return valuationDate_; - } - template - inline T Instrument::result(const std::string& tag) const { + T Instrument::result(const std::string& tag) const { calculate(); std::map::const_iterator value = additionalResults_.find(tag); - QL_REQUIRE(value != additionalResults_.end(), - tag << " not provided"); + QL_REQUIRE(value != additionalResults_.end(), tag << " not provided"); return ext::any_cast(value->second); } - inline const std::map& - Instrument::additionalResults() const { - calculate(); - return additionalResults_; - } - } - #endif diff --git a/ql/instruments/bond.hpp b/ql/instruments/bond.hpp index 8e0b0da6e8f..5c91e091bab 100644 --- a/ql/instruments/bond.hpp +++ b/ql/instruments/bond.hpp @@ -29,6 +29,7 @@ #ifndef quantlib_bond_hpp #define quantlib_bond_hpp +#include #include #include diff --git a/ql/instruments/compositeinstrument.cpp b/ql/instruments/compositeinstrument.cpp index 0ee7ddc2167..b78d8c346c0 100644 --- a/ql/instruments/compositeinstrument.cpp +++ b/ql/instruments/compositeinstrument.cpp @@ -17,6 +17,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include namespace QuantLib { diff --git a/ql/instruments/swap.hpp b/ql/instruments/swap.hpp index 5e39da6c04c..cbd460ec828 100644 --- a/ql/instruments/swap.hpp +++ b/ql/instruments/swap.hpp @@ -26,6 +26,8 @@ #ifndef quantlib_swap_hpp #define quantlib_swap_hpp +#include +#include #include #include #include diff --git a/ql/math/Makefile.am b/ql/math/Makefile.am index f7c8778f873..b481ee3cbca 100644 --- a/ql/math/Makefile.am +++ b/ql/math/Makefile.am @@ -40,6 +40,7 @@ this_include_HEADERS = \ cpp_files = \ abcdmathfunction.cpp \ + array.cpp \ bernsteinpolynomial.cpp \ beta.cpp \ bspline.cpp \ diff --git a/ql/math/array.cpp b/ql/math/array.cpp new file mode 100644 index 00000000000..6c2fe75a2cc --- /dev/null +++ b/ql/math/array.cpp @@ -0,0 +1,618 @@ +/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ + +/* + Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl + Copyright (C) 2003, 2004, 2005, 2006, 2009 StatPro Italia srl + Copyright (C) 2004 Ferdinando Ametrano + + 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 + . The license is also available online at + . + + 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. +*/ + +#include +#include +#include + +namespace QuantLib { + + Array::Array(Size size) + : data_(size != 0U ? new Real[size] : (Real*)nullptr), n_(size) {} + + Array::Array(Size size, Real value) + : data_(size != 0U ? new Real[size] : (Real*)nullptr), n_(size) { + std::fill(begin(),end(),value); + } + + Array::Array(Size size, Real value, Real increment) + : data_(size != 0U ? new Real[size] : (Real*)nullptr), n_(size) { + for (iterator i=begin(); i!=end(); ++i, value+=increment) + *i = value; + } + + Array::Array(const Array& from) + : data_(from.n_ != 0U ? new Real[from.n_] : (Real*)nullptr), n_(from.n_) { + if (data_) + std::copy(from.begin(),from.end(),begin()); + } + + Array::Array(Array&& from) noexcept + : data_((Real*)nullptr), n_(0) { + swap(from); + } + + Array::Array(std::initializer_list init) { + detail::_fill_array_(*this, data_, n_, init.begin(), init.end(), + std::false_type()); + } + + Array& Array::operator=(const Array& from) { + // strong guarantee + Array temp(from); + swap(temp); + return *this; + } + + Array& Array::operator=(Array&& from) noexcept { + swap(from); + return *this; + } + + bool Array::operator==(const Array& to) const { + return (n_ == to.n_) && std::equal(begin(), end(), to.begin()); + } + + bool Array::operator!=(const Array& to) const { + return !(this->operator==(to)); + } + + const Array& Array::operator+=(const Array& v) { + QL_REQUIRE(n_ == v.n_, + "arrays with different sizes (" << n_ << ", " + << v.n_ << ") cannot be added"); + std::transform(begin(),end(),v.begin(),begin(), std::plus<>()); + return *this; + } + + + const Array& Array::operator+=(Real x) { + std::transform(begin(), end(), begin(), [=](Real y) -> Real { return y + x; }); + return *this; + } + + const Array& Array::operator-=(const Array& v) { + QL_REQUIRE(n_ == v.n_, + "arrays with different sizes (" << n_ << ", " + << v.n_ << ") cannot be subtracted"); + std::transform(begin(), end(), v.begin(), begin(), std::minus<>()); + return *this; + } + + const Array& Array::operator-=(Real x) { + std::transform(begin(),end(),begin(), [=](Real y) -> Real { return y - x; }); + return *this; + } + + const Array& Array::operator*=(const Array& v) { + QL_REQUIRE(n_ == v.n_, + "arrays with different sizes (" << n_ << ", " + << v.n_ << ") cannot be multiplied"); + std::transform(begin(), end(), v.begin(), begin(), std::multiplies<>()); + return *this; + } + + const Array& Array::operator*=(Real x) { + std::transform(begin(), end(), begin(), [=](Real y) -> Real { return y * x; }); + return *this; + } + + const Array& Array::operator/=(const Array& v) { + QL_REQUIRE(n_ == v.n_, + "arrays with different sizes (" << n_ << ", " + << v.n_ << ") cannot be divided"); + std::transform(begin(), end(), v.begin(), begin(), std::divides<>()); + return *this; + } + + const Array& Array::operator/=(Real x) { + #if defined(QL_EXTRA_SAFETY_CHECKS) + QL_REQUIRE(x != 0.0, "division by zero"); + #endif + std::transform(begin(), end(), begin(), [=](Real y) -> Real { return y / x; }); + return *this; + } + + Real Array::operator[](Size i) const { + #if defined(QL_EXTRA_SAFETY_CHECKS) + QL_REQUIRE(i0, "null Array: array access out of range"); + #endif + return data_.get()[0]; + } + + Real Array::back() const { + #if defined(QL_EXTRA_SAFETY_CHECKS) + QL_REQUIRE(n_>0, "null Array: array access out of range"); + #endif + return data_.get()[n_-1]; + } + + Real& Array::operator[](Size i) { + #if defined(QL_EXTRA_SAFETY_CHECKS) + QL_REQUIRE(i0, "null Array: array access out of range"); + #endif + return data_.get()[0]; + } + + Real& Array::back() { + #if defined(QL_EXTRA_SAFETY_CHECKS) + QL_REQUIRE(n_>0, "null Array: array access out of range"); + #endif + return data_.get()[n_-1]; + } + + Size Array::size() const { + return n_; + } + + bool Array::empty() const { + return n_ == 0; + } + + Array::const_iterator Array::begin() const { + return data_.get(); + } + + Array::iterator Array::begin() { + return data_.get(); + } + + Array::const_iterator Array::end() const { + return data_.get()+n_; + } + + Array::iterator Array::end() { + return data_.get()+n_; + } + + Array::const_reverse_iterator Array::rbegin() const { + return const_reverse_iterator(end()); + } + + Array::reverse_iterator Array::rbegin() { + return reverse_iterator(end()); + } + + Array::const_reverse_iterator Array::rend() const { + return const_reverse_iterator(begin()); + } + + Array::reverse_iterator Array::rend() { + return reverse_iterator(begin()); + } + + void Array::resize(Size n) { + if (n > n_) { + Array swp(n); + std::copy(begin(), end(), swp.begin()); + swap(swp); + } + else if (n < n_) { + n_ = n; + } + } + + void Array::swap(Array& from) noexcept { + data_.swap(from.data_); + std::swap(n_, from.n_); + } + + // dot product and norm + + Real DotProduct(const Array& v1, const Array& v2) { + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be multiplied"); + return std::inner_product(v1.begin(),v1.end(),v2.begin(),Real(0.0)); + } + + Real Norm2(const Array& v) { + return std::sqrt(DotProduct(v, v)); + } + + // overloaded operators + + // unary + + Array operator+(const Array& v) { + Array result = v; + return result; + } + + Array operator+(Array&& v) { + return std::move(v); + } + + Array operator-(const Array& v) { + Array result(v.size()); + std::transform(v.begin(), v.end(), result.begin(), std::negate<>()); + return result; + } + + Array operator-(Array&& v) { + Array result = std::move(v); + std::transform(result.begin(), result.end(), result.begin(), std::negate<>()); + return result; + } + + // binary operators + + Array operator+(const Array& v1, const Array& v2) { + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be added"); + Array result(v1.size()); + std::transform(v1.begin(),v1.end(),v2.begin(),result.begin(), std::plus<>()); + return result; + } + + Array operator+(const Array& v1, Array&& v2) { + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be added"); + Array result = std::move(v2); + std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::plus<>()); + return result; + } + + Array operator+(Array&& v1, const Array& v2) { + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be added"); + Array result = std::move(v1); + std::transform(result.begin(), result.end(), v2.begin(), result.begin(), std::plus<>()); + return result; + } + + Array operator+(Array&& v1, Array&& v2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved) + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be added"); + Array result = std::move(v2); + std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::plus<>()); + return result; + } + + Array operator+(const Array& v1, Real a) { + Array result(v1.size()); + std::transform(v1.begin(), v1.end(), result.begin(), [=](Real y) -> Real { return y + a; }); + return result; + } + + Array operator+(Array&& v1, Real a) { + Array result = std::move(v1); + std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return y + a; }); + return result; + } + + Array operator+(Real a, const Array& v2) { + Array result(v2.size()); + std::transform(v2.begin(),v2.end(),result.begin(), [=](Real y) -> Real { return a + y; }); + return result; + } + + Array operator+(Real a, Array&& v2) { + Array result = std::move(v2); + std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return a + y; }); + return result; + } + + Array operator-(const Array& v1, const Array& v2) { + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be subtracted"); + Array result(v1.size()); + std::transform(v1.begin(), v1.end(), v2.begin(), result.begin(), std::minus<>()); + return result; + } + + Array operator-(const Array& v1, Array&& v2) { + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be subtracted"); + Array result = std::move(v2); + std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::minus<>()); + return result; + } + + Array operator-(Array&& v1, const Array& v2) { + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be subtracted"); + Array result = std::move(v1); + std::transform(result.begin(), result.end(), v2.begin(), result.begin(), std::minus<>()); + return result; + } + + Array operator-(Array&& v1, Array&& v2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved) + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be subtracted"); + Array result = std::move(v2); + std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::minus<>()); + return result; + } + + Array operator-(const Array& v1, Real a) { + Array result(v1.size()); + std::transform(v1.begin(),v1.end(),result.begin(), [=](Real y) -> Real { return y - a; }); + return result; + } + + Array operator-(Array&& v1, Real a) { + Array result = std::move(v1); + std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return y - a; }); + return result; + } + + Array operator-(Real a, const Array& v2) { + Array result(v2.size()); + std::transform(v2.begin(),v2.end(),result.begin(), [=](Real y) -> Real { return a - y; }); + return result; + } + + Array operator-(Real a, Array&& v2) { + Array result = std::move(v2); + std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return a - y; }); + return result; + } + + Array operator*(const Array& v1, const Array& v2) { + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be multiplied"); + Array result(v1.size()); + std::transform(v1.begin(), v1.end(), v2.begin(), result.begin(), std::multiplies<>()); + return result; + } + + Array operator*(const Array& v1, Array&& v2) { + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be multiplied"); + Array result = std::move(v2); + std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::multiplies<>()); + return result; + } + + Array operator*(Array&& v1, const Array& v2) { + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be multiplied"); + Array result = std::move(v1); + std::transform(result.begin(), result.end(), v2.begin(), result.begin(), std::multiplies<>()); + return result; + } + + Array operator*(Array&& v1, Array&& v2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved) + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be multiplied"); + Array result = std::move(v2); + std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::multiplies<>()); + return result; + } + + Array operator*(const Array& v1, Real a) { + Array result(v1.size()); + std::transform(v1.begin(),v1.end(),result.begin(), [=](Real y) -> Real { return y * a; }); + return result; + } + + Array operator*(Array&& v1, Real a) { + Array result = std::move(v1); + std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return y * a; }); + return result; + } + + Array operator*(Real a, const Array& v2) { + Array result(v2.size()); + std::transform(v2.begin(),v2.end(),result.begin(), [=](Real y) -> Real { return a * y; }); + return result; + } + + Array operator*(Real a, Array&& v2) { + Array result = std::move(v2); + std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return a * y; }); + return result; + } + + Array operator/(const Array& v1, const Array& v2) { + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be divided"); + Array result(v1.size()); + std::transform(v1.begin(), v1.end(), v2.begin(), result.begin(), std::divides<>()); + return result; + } + + Array operator/(const Array& v1, Array&& v2) { + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be divided"); + Array result = std::move(v2); + std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::divides<>()); + return result; + } + + Array operator/(Array&& v1, const Array& v2) { + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be divided"); + Array result = std::move(v1); + std::transform(result.begin(), result.end(), v2.begin(), result.begin(), std::divides<>()); + return result; + } + + Array operator/(Array&& v1, Array&& v2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved) + QL_REQUIRE(v1.size() == v2.size(), + "arrays with different sizes (" << v1.size() << ", " + << v2.size() << ") cannot be divided"); + Array result = std::move(v2); + std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::divides<>()); + return result; + } + + Array operator/(const Array& v1, Real a) { + Array result(v1.size()); + std::transform(v1.begin(),v1.end(),result.begin(), [=](Real y) -> Real { return y / a; }); + return result; + } + + Array operator/(Array&& v1, Real a) { + Array result = std::move(v1); + std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return y / a; }); + return result; + } + + Array operator/(Real a, const Array& v2) { + Array result(v2.size()); + std::transform(v2.begin(),v2.end(),result.begin(), [=](Real y) -> Real { return a / y; }); + return result; + } + + Array operator/(Real a, Array&& v2) { + Array result = std::move(v2); + std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return a / y; }); + return result; + } + + // functions + + Array Abs(const Array& v) { + Array result(v.size()); + std::transform(v.begin(), v.end(), result.begin(), + [](Real x) -> Real { return std::fabs(x); }); + return result; + } + + Array Abs(Array&& v) { + Array result = std::move(v); + std::transform(result.begin(), result.end(), result.begin(), + [](Real x) -> Real { return std::fabs(x); }); + return result; + } + + Array Sqrt(const Array& v) { + Array result(v.size()); + std::transform(v.begin(),v.end(),result.begin(), + [](Real x) -> Real { return std::sqrt(x); }); + return result; + } + + Array Sqrt(Array&& v) { + Array result = std::move(v); + std::transform(result.begin(), result.end(), result.begin(), + [](Real x) -> Real { return std::sqrt(x); }); + return result; + } + + Array Log(const Array& v) { + Array result(v.size()); + std::transform(v.begin(),v.end(),result.begin(), + [](Real x) -> Real { return std::log(x); }); + return result; + } + + Array Log(Array&& v) { + Array result = std::move(v); + std::transform(result.begin(), result.end(), result.begin(), + [](Real x) -> Real { return std::log(x); }); + return result; + } + + Array Exp(const Array& v) { + Array result(v.size()); + std::transform(v.begin(), v.end(), result.begin(), + [](Real x) -> Real { return std::exp(x); }); + return result; + } + + Array Exp(Array&& v) { + Array result = std::move(v); + std::transform(result.begin(), result.end(), result.begin(), + [](Real x) -> Real { return std::exp(x); }); + return result; + } + + Array Pow(const Array& v, Real alpha) { + Array result(v.size()); + std::transform(v.begin(), v.end(), result.begin(), + [=](Real x) -> Real { return std::pow(x, alpha); }); + return result; + } + + Array Pow(Array&& v, Real alpha) { + Array result = std::move(v); + std::transform(result.begin(), result.end(), result.begin(), + [=](Real x) -> Real { return std::pow(x, alpha); }); + return result; + } + + void swap(Array& v, Array& w) noexcept { + v.swap(w); + } + + std::ostream& operator<<(std::ostream& out, const Array& a) { + std::streamsize width = out.width(); + out << "[ "; + if (!a.empty()) { + for (Size n=0; n -#include #include #include #include @@ -35,7 +34,6 @@ #include #include #include -#include #include #include @@ -278,38 +276,10 @@ namespace QuantLib { /*! \relates Array */ std::ostream& operator<<(std::ostream&, const Array&); - - // inline definitions - - inline Array::Array(Size size) - : data_(size != 0U ? new Real[size] : (Real*)nullptr), n_(size) {} - - inline Array::Array(Size size, Real value) - : data_(size != 0U ? new Real[size] : (Real*)nullptr), n_(size) { - std::fill(begin(),end(),value); - } - - inline Array::Array(Size size, Real value, Real increment) - : data_(size != 0U ? new Real[size] : (Real*)nullptr), n_(size) { - for (iterator i=begin(); i!=end(); ++i, value+=increment) - *i = value; - } - - inline Array::Array(const Array& from) - : data_(from.n_ != 0U ? new Real[from.n_] : (Real*)nullptr), n_(from.n_) { - if (data_) - std::copy(from.begin(),from.end(),begin()); - } - - inline Array::Array(Array&& from) noexcept - : data_((Real*)nullptr), n_(0) { - swap(from); - } - namespace detail { template - inline void _fill_array_(Array& a, + void _fill_array_(Array& a, std::unique_ptr& data_, Size& n_, I begin, I end, @@ -326,7 +296,7 @@ namespace QuantLib { } template - inline void _fill_array_(Array& a, + void _fill_array_(Array& a, std::unique_ptr& data_, Size& n_, I begin, I end, @@ -343,579 +313,15 @@ namespace QuantLib { } - inline Array::Array(std::initializer_list init) { - detail::_fill_array_(*this, data_, n_, init.begin(), init.end(), - std::false_type()); - } - + /*! \relates Array */ template - inline Array::Array(ForwardIterator begin, ForwardIterator end) { + Array::Array(ForwardIterator begin, ForwardIterator end) { // Unfortunately, calls such as Array(3, 4) match this constructor. // We have to detect integral types and dispatch. detail::_fill_array_(*this, data_, n_, begin, end, std::is_integral()); } - inline Array& Array::operator=(const Array& from) { - // strong guarantee - Array temp(from); - swap(temp); - return *this; - } - - inline Array& Array::operator=(Array&& from) noexcept { - swap(from); - return *this; - } - - inline bool Array::operator==(const Array& to) const { - return (n_ == to.n_) && std::equal(begin(), end(), to.begin()); - } - - inline bool Array::operator!=(const Array& to) const { - return !(this->operator==(to)); - } - - inline const Array& Array::operator+=(const Array& v) { - QL_REQUIRE(n_ == v.n_, - "arrays with different sizes (" << n_ << ", " - << v.n_ << ") cannot be added"); - std::transform(begin(),end(),v.begin(),begin(), std::plus<>()); - return *this; - } - - - inline const Array& Array::operator+=(Real x) { - std::transform(begin(), end(), begin(), [=](Real y) -> Real { return y + x; }); - return *this; - } - - inline const Array& Array::operator-=(const Array& v) { - QL_REQUIRE(n_ == v.n_, - "arrays with different sizes (" << n_ << ", " - << v.n_ << ") cannot be subtracted"); - std::transform(begin(), end(), v.begin(), begin(), std::minus<>()); - return *this; - } - - inline const Array& Array::operator-=(Real x) { - std::transform(begin(),end(),begin(), [=](Real y) -> Real { return y - x; }); - return *this; - } - - inline const Array& Array::operator*=(const Array& v) { - QL_REQUIRE(n_ == v.n_, - "arrays with different sizes (" << n_ << ", " - << v.n_ << ") cannot be multiplied"); - std::transform(begin(), end(), v.begin(), begin(), std::multiplies<>()); - return *this; - } - - inline const Array& Array::operator*=(Real x) { - std::transform(begin(), end(), begin(), [=](Real y) -> Real { return y * x; }); - return *this; - } - - inline const Array& Array::operator/=(const Array& v) { - QL_REQUIRE(n_ == v.n_, - "arrays with different sizes (" << n_ << ", " - << v.n_ << ") cannot be divided"); - std::transform(begin(), end(), v.begin(), begin(), std::divides<>()); - return *this; - } - - inline const Array& Array::operator/=(Real x) { - #if defined(QL_EXTRA_SAFETY_CHECKS) - QL_REQUIRE(x != 0.0, "division by zero"); - #endif - std::transform(begin(), end(), begin(), [=](Real y) -> Real { return y / x; }); - return *this; - } - - inline Real Array::operator[](Size i) const { - #if defined(QL_EXTRA_SAFETY_CHECKS) - QL_REQUIRE(i0, "null Array: array access out of range"); - #endif - return data_.get()[0]; - } - - inline Real Array::back() const { - #if defined(QL_EXTRA_SAFETY_CHECKS) - QL_REQUIRE(n_>0, "null Array: array access out of range"); - #endif - return data_.get()[n_-1]; - } - - inline Real& Array::operator[](Size i) { - #if defined(QL_EXTRA_SAFETY_CHECKS) - QL_REQUIRE(i0, "null Array: array access out of range"); - #endif - return data_.get()[0]; - } - - inline Real& Array::back() { - #if defined(QL_EXTRA_SAFETY_CHECKS) - QL_REQUIRE(n_>0, "null Array: array access out of range"); - #endif - return data_.get()[n_-1]; - } - - inline Size Array::size() const { - return n_; - } - - inline bool Array::empty() const { - return n_ == 0; - } - - inline Array::const_iterator Array::begin() const { - return data_.get(); - } - - inline Array::iterator Array::begin() { - return data_.get(); - } - - inline Array::const_iterator Array::end() const { - return data_.get()+n_; - } - - inline Array::iterator Array::end() { - return data_.get()+n_; - } - - inline Array::const_reverse_iterator Array::rbegin() const { - return const_reverse_iterator(end()); - } - - inline Array::reverse_iterator Array::rbegin() { - return reverse_iterator(end()); - } - - inline Array::const_reverse_iterator Array::rend() const { - return const_reverse_iterator(begin()); - } - - inline Array::reverse_iterator Array::rend() { - return reverse_iterator(begin()); - } - - inline void Array::resize(Size n) { - if (n > n_) { - Array swp(n); - std::copy(begin(), end(), swp.begin()); - swap(swp); - } - else if (n < n_) { - n_ = n; - } - } - - inline void Array::swap(Array& from) noexcept { - data_.swap(from.data_); - std::swap(n_, from.n_); - } - - // dot product and norm - - inline Real DotProduct(const Array& v1, const Array& v2) { - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be multiplied"); - return std::inner_product(v1.begin(),v1.end(),v2.begin(),Real(0.0)); - } - - inline Real Norm2(const Array& v) { - return std::sqrt(DotProduct(v, v)); - } - - // overloaded operators - - // unary - - inline Array operator+(const Array& v) { - Array result = v; - return result; - } - - inline Array operator+(Array&& v) { - return std::move(v); - } - - inline Array operator-(const Array& v) { - Array result(v.size()); - std::transform(v.begin(), v.end(), result.begin(), std::negate<>()); - return result; - } - - inline Array operator-(Array&& v) { - Array result = std::move(v); - std::transform(result.begin(), result.end(), result.begin(), std::negate<>()); - return result; - } - - // binary operators - - inline Array operator+(const Array& v1, const Array& v2) { - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be added"); - Array result(v1.size()); - std::transform(v1.begin(),v1.end(),v2.begin(),result.begin(), std::plus<>()); - return result; - } - - inline Array operator+(const Array& v1, Array&& v2) { - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be added"); - Array result = std::move(v2); - std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::plus<>()); - return result; - } - - inline Array operator+(Array&& v1, const Array& v2) { - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be added"); - Array result = std::move(v1); - std::transform(result.begin(), result.end(), v2.begin(), result.begin(), std::plus<>()); - return result; - } - - inline Array operator+(Array&& v1, Array&& v2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved) - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be added"); - Array result = std::move(v2); - std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::plus<>()); - return result; - } - - inline Array operator+(const Array& v1, Real a) { - Array result(v1.size()); - std::transform(v1.begin(), v1.end(), result.begin(), [=](Real y) -> Real { return y + a; }); - return result; - } - - inline Array operator+(Array&& v1, Real a) { - Array result = std::move(v1); - std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return y + a; }); - return result; - } - - inline Array operator+(Real a, const Array& v2) { - Array result(v2.size()); - std::transform(v2.begin(),v2.end(),result.begin(), [=](Real y) -> Real { return a + y; }); - return result; - } - - inline Array operator+(Real a, Array&& v2) { - Array result = std::move(v2); - std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return a + y; }); - return result; - } - - inline Array operator-(const Array& v1, const Array& v2) { - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be subtracted"); - Array result(v1.size()); - std::transform(v1.begin(), v1.end(), v2.begin(), result.begin(), std::minus<>()); - return result; - } - - inline Array operator-(const Array& v1, Array&& v2) { - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be subtracted"); - Array result = std::move(v2); - std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::minus<>()); - return result; - } - - inline Array operator-(Array&& v1, const Array& v2) { - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be subtracted"); - Array result = std::move(v1); - std::transform(result.begin(), result.end(), v2.begin(), result.begin(), std::minus<>()); - return result; - } - - inline Array operator-(Array&& v1, Array&& v2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved) - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be subtracted"); - Array result = std::move(v2); - std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::minus<>()); - return result; - } - - inline Array operator-(const Array& v1, Real a) { - Array result(v1.size()); - std::transform(v1.begin(),v1.end(),result.begin(), [=](Real y) -> Real { return y - a; }); - return result; - } - - inline Array operator-(Array&& v1, Real a) { - Array result = std::move(v1); - std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return y - a; }); - return result; - } - - inline Array operator-(Real a, const Array& v2) { - Array result(v2.size()); - std::transform(v2.begin(),v2.end(),result.begin(), [=](Real y) -> Real { return a - y; }); - return result; - } - - inline Array operator-(Real a, Array&& v2) { - Array result = std::move(v2); - std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return a - y; }); - return result; - } - - inline Array operator*(const Array& v1, const Array& v2) { - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be multiplied"); - Array result(v1.size()); - std::transform(v1.begin(), v1.end(), v2.begin(), result.begin(), std::multiplies<>()); - return result; - } - - inline Array operator*(const Array& v1, Array&& v2) { - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be multiplied"); - Array result = std::move(v2); - std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::multiplies<>()); - return result; - } - - inline Array operator*(Array&& v1, const Array& v2) { - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be multiplied"); - Array result = std::move(v1); - std::transform(result.begin(), result.end(), v2.begin(), result.begin(), std::multiplies<>()); - return result; - } - - inline Array operator*(Array&& v1, Array&& v2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved) - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be multiplied"); - Array result = std::move(v2); - std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::multiplies<>()); - return result; - } - - inline Array operator*(const Array& v1, Real a) { - Array result(v1.size()); - std::transform(v1.begin(),v1.end(),result.begin(), [=](Real y) -> Real { return y * a; }); - return result; - } - - inline Array operator*(Array&& v1, Real a) { - Array result = std::move(v1); - std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return y * a; }); - return result; - } - - inline Array operator*(Real a, const Array& v2) { - Array result(v2.size()); - std::transform(v2.begin(),v2.end(),result.begin(), [=](Real y) -> Real { return a * y; }); - return result; - } - - inline Array operator*(Real a, Array&& v2) { - Array result = std::move(v2); - std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return a * y; }); - return result; - } - - inline Array operator/(const Array& v1, const Array& v2) { - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be divided"); - Array result(v1.size()); - std::transform(v1.begin(), v1.end(), v2.begin(), result.begin(), std::divides<>()); - return result; - } - - inline Array operator/(const Array& v1, Array&& v2) { - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be divided"); - Array result = std::move(v2); - std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::divides<>()); - return result; - } - - inline Array operator/(Array&& v1, const Array& v2) { - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be divided"); - Array result = std::move(v1); - std::transform(result.begin(), result.end(), v2.begin(), result.begin(), std::divides<>()); - return result; - } - - inline Array operator/(Array&& v1, Array&& v2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved) - QL_REQUIRE(v1.size() == v2.size(), - "arrays with different sizes (" << v1.size() << ", " - << v2.size() << ") cannot be divided"); - Array result = std::move(v2); - std::transform(v1.begin(), v1.end(), result.begin(), result.begin(), std::divides<>()); - return result; - } - - inline Array operator/(const Array& v1, Real a) { - Array result(v1.size()); - std::transform(v1.begin(),v1.end(),result.begin(), [=](Real y) -> Real { return y / a; }); - return result; - } - - inline Array operator/(Array&& v1, Real a) { - Array result = std::move(v1); - std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return y / a; }); - return result; - } - - inline Array operator/(Real a, const Array& v2) { - Array result(v2.size()); - std::transform(v2.begin(),v2.end(),result.begin(), [=](Real y) -> Real { return a / y; }); - return result; - } - - inline Array operator/(Real a, Array&& v2) { - Array result = std::move(v2); - std::transform(result.begin(), result.end(), result.begin(), [=](Real y) -> Real { return a / y; }); - return result; - } - - // functions - - inline Array Abs(const Array& v) { - Array result(v.size()); - std::transform(v.begin(), v.end(), result.begin(), - [](Real x) -> Real { return std::fabs(x); }); - return result; - } - - inline Array Abs(Array&& v) { - Array result = std::move(v); - std::transform(result.begin(), result.end(), result.begin(), - [](Real x) -> Real { return std::fabs(x); }); - return result; - } - - inline Array Sqrt(const Array& v) { - Array result(v.size()); - std::transform(v.begin(),v.end(),result.begin(), - [](Real x) -> Real { return std::sqrt(x); }); - return result; - } - - inline Array Sqrt(Array&& v) { - Array result = std::move(v); - std::transform(result.begin(), result.end(), result.begin(), - [](Real x) -> Real { return std::sqrt(x); }); - return result; - } - - inline Array Log(const Array& v) { - Array result(v.size()); - std::transform(v.begin(),v.end(),result.begin(), - [](Real x) -> Real { return std::log(x); }); - return result; - } - - inline Array Log(Array&& v) { - Array result = std::move(v); - std::transform(result.begin(), result.end(), result.begin(), - [](Real x) -> Real { return std::log(x); }); - return result; - } - - inline Array Exp(const Array& v) { - Array result(v.size()); - std::transform(v.begin(), v.end(), result.begin(), - [](Real x) -> Real { return std::exp(x); }); - return result; - } - - inline Array Exp(Array&& v) { - Array result = std::move(v); - std::transform(result.begin(), result.end(), result.begin(), - [](Real x) -> Real { return std::exp(x); }); - return result; - } - - inline Array Pow(const Array& v, Real alpha) { - Array result(v.size()); - std::transform(v.begin(), v.end(), result.begin(), - [=](Real x) -> Real { return std::pow(x, alpha); }); - return result; - } - - inline Array Pow(Array&& v, Real alpha) { - Array result = std::move(v); - std::transform(result.begin(), result.end(), result.begin(), - [=](Real x) -> Real { return std::pow(x, alpha); }); - return result; - } - - inline void swap(Array& v, Array& w) noexcept { - v.swap(w); - } - - inline std::ostream& operator<<(std::ostream& out, const Array& a) { - std::streamsize width = out.width(); - out << "[ "; - if (!a.empty()) { - for (Size n=0; n #include #include #include diff --git a/ql/math/integrals/gaussianquadratures.hpp b/ql/math/integrals/gaussianquadratures.hpp index 9eb9e70f38b..4545c400429 100644 --- a/ql/math/integrals/gaussianquadratures.hpp +++ b/ql/math/integrals/gaussianquadratures.hpp @@ -25,9 +25,11 @@ #ifndef quantlib_gaussian_quadratures_hpp #define quantlib_gaussian_quadratures_hpp +#include #include #include #include +#include namespace QuantLib { class GaussianOrthogonalPolynomial; diff --git a/ql/math/interpolations/Makefile.am b/ql/math/interpolations/Makefile.am index b56c0fb21ed..e4f2d01e980 100644 --- a/ql/math/interpolations/Makefile.am +++ b/ql/math/interpolations/Makefile.am @@ -26,7 +26,8 @@ this_include_HEADERS = \ xabrinterpolation.hpp cpp_files = \ - chebyshevinterpolation.cpp + chebyshevinterpolation.cpp \ + cubicinterpolation.cpp if UNITY_BUILD diff --git a/ql/math/interpolations/cubicinterpolation.cpp b/ql/math/interpolations/cubicinterpolation.cpp new file mode 100644 index 00000000000..89019c49342 --- /dev/null +++ b/ql/math/interpolations/cubicinterpolation.cpp @@ -0,0 +1,334 @@ +/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ + +/* + Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl + Copyright (C) 2001, 2002, 2003 Nicolas Di Césaré + Copyright (C) 2004, 2008, 2009, 2011 Ferdinando Ametrano + Copyright (C) 2009 Sylvain Bertrand + Copyright (C) 2013 Peter Caspers + Copyright (C) 2016 Nicholas Bertocchi + + 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 + . The license is also available online at + . + + 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. +*/ + +#include + +namespace QuantLib { + + namespace detail { + + void cubicInterpolationSplineOM1(Array& tmp_, std::vector& dx_, Array& Y_, Size n_) { + Matrix T_(n_ - 2, n_, 0.0); + for (Size i = 0; i < n_ - 2; ++i) { + T_[i][i] = dx_[i] / 6.0; + T_[i][i + 1] = (dx_[i + 1] + dx_[i]) / 3.0; + T_[i][i + 2] = dx_[i + 1] / 6.0; + } + Matrix S_(n_ - 2, n_, 0.0); + for (Size i = 0; i < n_ - 2; ++i) { + S_[i][i] = 1.0 / dx_[i]; + S_[i][i + 1] = -(1.0 / dx_[i + 1] + 1.0 / dx_[i]); + S_[i][i + 2] = 1.0 / dx_[i + 1]; + } + Matrix Up_(n_, 2, 0.0); + Up_[0][0] = 1; + Up_[n_ - 1][1] = 1; + Matrix Us_(n_, n_ - 2, 0.0); + for (Size i = 0; i < n_ - 2; ++i) + Us_[i + 1][i] = 1; + Matrix Z_ = Us_ * inverse(T_ * Us_); + Matrix I_(n_, n_, 0.0); + for (Size i = 0; i < n_; ++i) + I_[i][i] = 1; + Matrix V_ = (I_ - Z_ * T_) * Up_; + Matrix W_ = Z_ * S_; + Matrix Q_(n_, n_, 0.0); + Q_[0][0] = 1.0 / (n_ - 1) * dx_[0] * dx_[0] * dx_[0]; + Q_[0][1] = 7.0 / 8 * 1.0 / (n_ - 1) * dx_[0] * dx_[0] * dx_[0]; + for (Size i = 1; i < n_ - 1; ++i) { + Q_[i][i - 1] = 7.0 / 8 * 1.0 / (n_ - 1) * dx_[i - 1] * dx_[i - 1] * dx_[i - 1]; + Q_[i][i] = 1.0 / (n_ - 1) * dx_[i] * dx_[i] * dx_[i] + + 1.0 / (n_ - 1) * dx_[i - 1] * dx_[i - 1] * dx_[i - 1]; + Q_[i][i + 1] = 7.0 / 8 * 1.0 / (n_ - 1) * dx_[i] * dx_[i] * dx_[i]; + } + Q_[n_ - 1][n_ - 2] = 7.0 / 8 * 1.0 / (n_ - 1) * dx_[n_ - 2] * dx_[n_ - 2] * dx_[n_ - 2]; + Q_[n_ - 1][n_ - 1] = 1.0 / (n_ - 1) * dx_[n_ - 2] * dx_[n_ - 2] * dx_[n_ - 2]; + Matrix J_ = (I_ - V_ * inverse(transpose(V_) * Q_ * V_) * transpose(V_) * Q_) * W_; + Array D_ = J_ * Y_; + for (Size i = 0; i < n_ - 1; ++i) + tmp_[i] = (Y_[i + 1] - Y_[i]) / dx_[i] - (2.0 * D_[i] + D_[i + 1]) * dx_[i] / 6.0; + tmp_[n_ - 1] = tmp_[n_ - 2] + D_[n_ - 2] * dx_[n_ - 2] + + (D_[n_ - 1] - D_[n_ - 2]) * dx_[n_ - 2] / 2.0; + } + + void cubicInterpolationSplineOM2(Array& tmp_, std::vector& dx_, Array& Y_, Size n_) { + Matrix T_(n_ - 2, n_, 0.0); + for (Size i = 0; i < n_ - 2; ++i) { + T_[i][i] = dx_[i] / 6.0; + T_[i][i + 1] = (dx_[i] + dx_[i + 1]) / 3.0; + T_[i][i + 2] = dx_[i + 1] / 6.0; + } + Matrix S_(n_ - 2, n_, 0.0); + for (Size i = 0; i < n_ - 2; ++i) { + S_[i][i] = 1.0 / dx_[i]; + S_[i][i + 1] = -(1.0 / dx_[i + 1] + 1.0 / dx_[i]); + S_[i][i + 2] = 1.0 / dx_[i + 1]; + } + Matrix Up_(n_, 2, 0.0); + Up_[0][0] = 1; + Up_[n_ - 1][1] = 1; + Matrix Us_(n_, n_ - 2, 0.0); + for (Size i = 0; i < n_ - 2; ++i) + Us_[i + 1][i] = 1; + Matrix Z_ = Us_ * inverse(T_ * Us_); + Matrix I_(n_, n_, 0.0); + for (Size i = 0; i < n_; ++i) + I_[i][i] = 1; + Matrix V_ = (I_ - Z_ * T_) * Up_; + Matrix W_ = Z_ * S_; + Matrix Q_(n_, n_, 0.0); + Q_[0][0] = 1.0 / (n_ - 1) * dx_[0]; + Q_[0][1] = 1.0 / 2 * 1.0 / (n_ - 1) * dx_[0]; + for (Size i = 1; i < n_ - 1; ++i) { + Q_[i][i - 1] = 1.0 / 2 * 1.0 / (n_ - 1) * dx_[i - 1]; + Q_[i][i] = 1.0 / (n_ - 1) * dx_[i] + 1.0 / (n_ - 1) * dx_[i - 1]; + Q_[i][i + 1] = 1.0 / 2 * 1.0 / (n_ - 1) * dx_[i]; + } + Q_[n_ - 1][n_ - 2] = 1.0 / 2 * 1.0 / (n_ - 1) * dx_[n_ - 2]; + Q_[n_ - 1][n_ - 1] = 1.0 / (n_ - 1) * dx_[n_ - 2]; + Matrix J_ = (I_ - V_ * inverse(transpose(V_) * Q_ * V_) * transpose(V_) * Q_) * W_; + Array D_ = J_ * Y_; + for (Size i = 0; i < n_ - 1; ++i) + tmp_[i] = (Y_[i + 1] - Y_[i]) / dx_[i] - (2.0 * D_[i] + D_[i + 1]) * dx_[i] / 6.0; + tmp_[n_ - 1] = tmp_[n_ - 2] + D_[n_ - 2] * dx_[n_ - 2] + + (D_[n_ - 1] - D_[n_ - 2]) * dx_[n_ - 2] / 2.0; + } + + void cubicInterpolationSplineParabolic(Array& tmp_, + std::vector& dx_, + std::vector& S_, + Size n_) { + // intermediate points + for (Size i = 1; i < n_ - 1; ++i) + tmp_[i] = (dx_[i - 1] * S_[i] + dx_[i] * S_[i - 1]) / (dx_[i] + dx_[i - 1]); + // end points + tmp_[0] = ((2.0 * dx_[0] + dx_[1]) * S_[0] - dx_[0] * S_[1]) / (dx_[0] + dx_[1]); + tmp_[n_ - 1] = + ((2.0 * dx_[n_ - 2] + dx_[n_ - 3]) * S_[n_ - 2] - dx_[n_ - 2] * S_[n_ - 3]) / + (dx_[n_ - 2] + dx_[n_ - 3]); + } + + void cubicInterpolationSplineFritschButland(Array& tmp_, + std::vector& dx_, + std::vector& S_, + Size n_) { + // intermediate points + for (Size i = 1; i < n_ - 1; ++i) { + Real Smin = std::min(S_[i - 1], S_[i]); + Real Smax = std::max(S_[i - 1], S_[i]); + if (Smax + 2.0 * Smin == 0) { + if (Smin * Smax < 0) + tmp_[i] = QL_MIN_REAL; + else if (Smin * Smax == 0) + tmp_[i] = 0; + else + tmp_[i] = QL_MAX_REAL; + } else + tmp_[i] = 3.0 * Smin * Smax / (Smax + 2.0 * Smin); + } + // end points + tmp_[0] = ((2.0 * dx_[0] + dx_[1]) * S_[0] - dx_[0] * S_[1]) / (dx_[0] + dx_[1]); + tmp_[n_ - 1] = + ((2.0 * dx_[n_ - 2] + dx_[n_ - 3]) * S_[n_ - 2] - dx_[n_ - 2] * S_[n_ - 3]) / + (dx_[n_ - 2] + dx_[n_ - 3]); + } + + void cubicInterpolationSplineAkima(Array& tmp_, + std::vector& dx_, + std::vector& S_, + Size n_) { + tmp_[0] = + (std::abs(S_[1] - S_[0]) * 2 * S_[0] * S_[1] + + std::abs(2 * S_[0] * S_[1] - 4 * S_[0] * S_[0] * S_[1]) * S_[0]) / + (std::abs(S_[1] - S_[0]) + std::abs(2 * S_[0] * S_[1] - 4 * S_[0] * S_[0] * S_[1])); + tmp_[1] = + (std::abs(S_[2] - S_[1]) * S_[0] + std::abs(S_[0] - 2 * S_[0] * S_[1]) * S_[1]) / + (std::abs(S_[2] - S_[1]) + std::abs(S_[0] - 2 * S_[0] * S_[1])); + for (Size i = 2; i < n_ - 2; ++i) { + if ((S_[i - 2] == S_[i - 1]) && (S_[i] != S_[i + 1])) + tmp_[i] = S_[i - 1]; + else if ((S_[i - 2] != S_[i - 1]) && (S_[i] == S_[i + 1])) + tmp_[i] = S_[i]; + else if (S_[i] == S_[i - 1]) + tmp_[i] = S_[i]; + else if ((S_[i - 2] == S_[i - 1]) && (S_[i - 1] != S_[i]) && (S_[i] == S_[i + 1])) + tmp_[i] = (S_[i - 1] + S_[i]) / 2.0; + else + tmp_[i] = (std::abs(S_[i + 1] - S_[i]) * S_[i - 1] + + std::abs(S_[i - 1] - S_[i - 2]) * S_[i]) / + (std::abs(S_[i + 1] - S_[i]) + std::abs(S_[i - 1] - S_[i - 2])); + } + tmp_[n_ - 2] = (std::abs(2 * S_[n_ - 2] * S_[n_ - 3] - S_[n_ - 2]) * S_[n_ - 3] + + std::abs(S_[n_ - 3] - S_[n_ - 4]) * S_[n_ - 2]) / + (std::abs(2 * S_[n_ - 2] * S_[n_ - 3] - S_[n_ - 2]) + + std::abs(S_[n_ - 3] - S_[n_ - 4])); + tmp_[n_ - 1] = + (std::abs(4 * S_[n_ - 2] * S_[n_ - 2] * S_[n_ - 3] - 2 * S_[n_ - 2] * S_[n_ - 3]) * + S_[n_ - 2] + + std::abs(S_[n_ - 2] - S_[n_ - 3]) * 2 * S_[n_ - 2] * S_[n_ - 3]) / + (std::abs(4 * S_[n_ - 2] * S_[n_ - 2] * S_[n_ - 3] - 2 * S_[n_ - 2] * S_[n_ - 3]) + + std::abs(S_[n_ - 2] - S_[n_ - 3])); + } + + void cubicInterpolationSplineKruger(Array& tmp_, + std::vector& dx_, + std::vector& S_, + Size n_) { + // intermediate points + for (Size i = 1; i < n_ - 1; ++i) { + if (S_[i - 1] * S_[i] < 0.0) + // slope changes sign at point + tmp_[i] = 0.0; + else + // slope will be between the slopes of the adjacent + // straight lines and should approach zero if the + // slope of either line approaches zero + tmp_[i] = 2.0 / (1.0 / S_[i - 1] + 1.0 / S_[i]); + } + // end points + tmp_[0] = (3.0 * S_[0] - tmp_[1]) / 2.0; + tmp_[n_ - 1] = (3.0 * S_[n_ - 2] - tmp_[n_ - 2]) / 2.0; + } + + void cubicInterpolationSplineHarmonic(Array& tmp_, + std::vector& dx_, + std::vector& S_, + Size n_) { + // intermediate points + for (Size i = 1; i < n_ - 1; ++i) { + Real w1 = 2 * dx_[i] + dx_[i - 1]; + Real w2 = dx_[i] + 2 * dx_[i - 1]; + if (S_[i - 1] * S_[i] <= 0.0) + // slope changes sign at point + tmp_[i] = 0.0; + else + // weighted harmonic mean of S_[i] and S_[i-1] if they + // have the same sign; otherwise 0 + tmp_[i] = (w1 + w2) / (w1 / S_[i - 1] + w2 / S_[i]); + } + // end points [0] + tmp_[0] = ((2 * dx_[0] + dx_[1]) * S_[0] - dx_[0] * S_[1]) / (dx_[1] + dx_[0]); + if (tmp_[0] * S_[0] < 0.0) { + tmp_[0] = 0; + } else if (S_[0] * S_[1] < 0) { + if (std::fabs(tmp_[0]) > std::fabs(3 * S_[0])) { + tmp_[0] = 3 * S_[0]; + } + } + // end points [n-1] + tmp_[n_ - 1] = + ((2 * dx_[n_ - 2] + dx_[n_ - 3]) * S_[n_ - 2] - dx_[n_ - 2] * S_[n_ - 3]) / + (dx_[n_ - 3] + dx_[n_ - 2]); + if (tmp_[n_ - 1] * S_[n_ - 2] < 0.0) { + tmp_[n_ - 1] = 0; + } else if (S_[n_ - 2] * S_[n_ - 3] < 0) { + if (std::fabs(tmp_[n_ - 1]) > std::fabs(3 * S_[n_ - 2])) { + tmp_[n_ - 1] = 3 * S_[n_ - 2]; + } + } + } + + void cubicInterpolationSplineMonotonicFilter(Array& tmp_, + std::vector& dx_, + std::vector& S_, + Size n_, + std::vector& monotonicityAdjustments_) { + Real correction; + Real pm, pu, pd, M; + for (Size i = 0; i < n_; ++i) { + if (i == 0) { + if (tmp_[i] * S_[0] > 0.0) { + correction = tmp_[i] / std::fabs(tmp_[i]) * + std::min(std::fabs(tmp_[i]), std::fabs(3.0 * S_[0])); + } else { + correction = 0.0; + } + if (correction != tmp_[i]) { + tmp_[i] = correction; + monotonicityAdjustments_[i] = true; + } + } else if (i == n_ - 1) { + if (tmp_[i] * S_[n_ - 2] > 0.0) { + correction = + tmp_[i] / std::fabs(tmp_[i]) * + std::min(std::fabs(tmp_[i]), std::fabs(3.0 * S_[n_ - 2])); + } else { + correction = 0.0; + } + if (correction != tmp_[i]) { + tmp_[i] = correction; + monotonicityAdjustments_[i] = true; + } + } else { + pm = (S_[i - 1] * dx_[i] + S_[i] * dx_[i - 1]) / (dx_[i - 1] + dx_[i]); + M = 3.0 * + std::min(std::min(std::fabs(S_[i - 1]), std::fabs(S_[i])), std::fabs(pm)); + if (i > 1) { + if ((S_[i - 1] - S_[i - 2]) * (S_[i] - S_[i - 1]) > 0.0) { + pd = (S_[i - 1] * (2.0 * dx_[i - 1] + dx_[i - 2]) - + S_[i - 2] * dx_[i - 1]) / + (dx_[i - 2] + dx_[i - 1]); + if (pm * pd > 0.0 && pm * (S_[i - 1] - S_[i - 2]) > 0.0) { + M = std::max(M, 1.5 * std::min(std::fabs(pm), std::fabs(pd))); + } + } + } + if (i < n_ - 2) { + if ((S_[i] - S_[i - 1]) * (S_[i + 1] - S_[i]) > 0.0) { + pu = (S_[i] * (2.0 * dx_[i] + dx_[i + 1]) - S_[i + 1] * dx_[i]) / + (dx_[i] + dx_[i + 1]); + if (pm * pu > 0.0 && -pm * (S_[i] - S_[i - 1]) > 0.0) { + M = std::max(M, 1.5 * std::min(std::fabs(pm), std::fabs(pu))); + } + } + } + if (tmp_[i] * pm > 0.0) { + correction = tmp_[i] / std::fabs(tmp_[i]) * std::min(std::fabs(tmp_[i]), M); + } else { + correction = 0.0; + } + if (correction != tmp_[i]) { + tmp_[i] = correction; + monotonicityAdjustments_[i] = true; + } + } + } + } + + Real cubicInterpolatingPolynomialDerivative( + Real a, Real b, Real c, Real d, Real u, Real v, Real w, Real z, Real x) { + return (-((((a - c) * (b - c) * (c - x) * z - (a - d) * (b - d) * (d - x) * w) * + (a - x + b - x) + + ((a - c) * (b - c) * z - (a - d) * (b - d) * w) * (a - x) * (b - x)) * + (a - b) + + ((a - c) * (a - d) * v - (b - c) * (b - d) * u) * (c - d) * (c - x) * + (d - x) + + ((a - c) * (a - d) * (a - x) * v - (b - c) * (b - d) * (b - x) * u) * + (c - x + d - x) * (c - d))) / + ((a - b) * (a - c) * (a - d) * (b - c) * (b - d) * (c - d)); + } + + } // namespace detail + +} // namespace QuantLib diff --git a/ql/math/interpolations/cubicinterpolation.hpp b/ql/math/interpolations/cubicinterpolation.hpp index 0bb40efbfd4..45386abd188 100644 --- a/ql/math/interpolations/cubicinterpolation.hpp +++ b/ql/math/interpolations/cubicinterpolation.hpp @@ -29,9 +29,10 @@ #ifndef quantlib_cubic_interpolation_hpp #define quantlib_cubic_interpolation_hpp -#include #include +#include #include + #include namespace QuantLib { @@ -359,6 +360,51 @@ namespace QuantLib { namespace detail { + void cubicInterpolationSplineOM1(Array& tmp_, + std::vector& dx_, + Array& Y, + Size n_); + + + void cubicInterpolationSplineOM2(Array& tmp_, + std::vector& dx_, + Array& Y, + Size n_); + + void cubicInterpolationSplineParabolic(Array& tmp_, + std::vector& dx_, + std::vector& S_, + Size n_); + + void cubicInterpolationSplineFritschButland(Array& tmp_, + std::vector& dx_, + std::vector& S_, + Size n_); + + void cubicInterpolationSplineAkima(Array& tmp_, + std::vector& dx_, + std::vector& S_, + Size n_); + + void cubicInterpolationSplineKruger(Array& tmp_, + std::vector& dx_, + std::vector& S_, + Size n_); + + void cubicInterpolationSplineHarmonic(Array& tmp_, + std::vector& dx_, + std::vector& S_, + Size n_); + + void cubicInterpolationSplineMonotonicFilter(Array& tmp_, + std::vector& dx_, + std::vector& S_, + Size n_, + std::vector& monotonicityAdjustments_); + + Real cubicInterpolatingPolynomialDerivative( + Real a, Real b, Real c, Real d, Real u, Real v, Real w, Real z, Real x); + template class CubicInterpolationImpl final : public CoefficientHolder, public Interpolation::templateImpl { @@ -470,92 +516,16 @@ namespace QuantLib { // solve the system L_.solveFor(tmp_, tmp_); } else if (da_==CubicInterpolation::SplineOM1) { - Matrix T_(n_-2, n_, 0.0); - for (Size i=0; iyBegin_[i]; - Array D_ = J_*Y_; - for (Size i=0; iyBegin_[i]; - Array D_ = J_*Y_; - for (Size i=0; istd::fabs(3*S_[0])) { - tmp_[0] = 3*S_[0]; - } - } - // end points [n-1] - tmp_[n_-1] = ((2*dx_[n_-2]+dx_[n_-3])*S_[n_-2]-dx_[n_-2]*S_[n_-3])/(dx_[n_-3]+dx_[n_-2]); - if (tmp_[n_-1]*S_[n_-2]<0.0) { - tmp_[n_-1] = 0; - } - else if (S_[n_-2]*S_[n_-3]<0) { - if (std::fabs(tmp_[n_-1])>std::fabs(3*S_[n_-2])) { - tmp_[n_-1] = 3*S_[n_-2]; - } - } + cubicInterpolationSplineHarmonic(tmp_, dx_, S_, n_); break; default: QL_FAIL("unknown scheme"); @@ -670,72 +559,8 @@ namespace QuantLib { monotonicityAdjustments_.end(), false); // Hyman monotonicity constrained filter if (monotonic_) { - Real correction; - Real pm, pu, pd, M; - for (Size i=0; i0.0) { - correction = tmp_[i]/std::fabs(tmp_[i]) * - std::min(std::fabs(tmp_[i]), - std::fabs(3.0*S_[0])); - } else { - correction = 0.0; - } - if (correction!=tmp_[i]) { - tmp_[i] = correction; - monotonicityAdjustments_[i] = true; - } - } else if (i==n_-1) { - if (tmp_[i]*S_[n_-2]>0.0) { - correction = tmp_[i]/std::fabs(tmp_[i]) * - std::min(std::fabs(tmp_[i]), - std::fabs(3.0*S_[n_-2])); - } else { - correction = 0.0; - } - if (correction!=tmp_[i]) { - tmp_[i] = correction; - monotonicityAdjustments_[i] = true; - } - } else { - pm=(S_[i-1]*dx_[i]+S_[i]*dx_[i-1])/ - (dx_[i-1]+dx_[i]); - M = 3.0 * std::min(std::min(std::fabs(S_[i-1]), - std::fabs(S_[i])), - std::fabs(pm)); - if (i>1) { - if ((S_[i-1]-S_[i-2])*(S_[i]-S_[i-1])>0.0) { - pd=(S_[i-1]*(2.0*dx_[i-1]+dx_[i-2]) - -S_[i-2]*dx_[i-1])/ - (dx_[i-2]+dx_[i-1]); - if (pm*pd>0.0 && pm*(S_[i-1]-S_[i-2])>0.0) { - M = std::max(M, 1.5*std::min( - std::fabs(pm),std::fabs(pd))); - } - } - } - if (i0.0) { - pu=(S_[i]*(2.0*dx_[i]+dx_[i+1])-S_[i+1]*dx_[i])/ - (dx_[i]+dx_[i+1]); - if (pm*pu>0.0 && -pm*(S_[i]-S_[i-1])>0.0) { - M = std::max(M, 1.5*std::min( - std::fabs(pm),std::fabs(pu))); - } - } - } - if (tmp_[i]*pm>0.0) { - correction = tmp_[i]/std::fabs(tmp_[i]) * - std::min(std::fabs(tmp_[i]), M); - } else { - correction = 0.0; - } - if (correction!=tmp_[i]) { - tmp_[i] = correction; - monotonicityAdjustments_[i] = true; - } - } - } + cubicInterpolationSplineMonotonicFilter(tmp_, dx_, S_, n_, + monotonicityAdjustments_); } @@ -786,17 +611,6 @@ namespace QuantLib { mutable Array tmp_; mutable std::vector dx_, S_; mutable TridiagonalOperator L_; - - inline Real cubicInterpolatingPolynomialDerivative( - Real a, Real b, Real c, Real d, - Real u, Real v, Real w, Real z, Real x) const { - return (-((((a-c)*(b-c)*(c-x)*z-(a-d)*(b-d)*(d-x)*w)*(a-x+b-x) - +((a-c)*(b-c)*z-(a-d)*(b-d)*w)*(a-x)*(b-x))*(a-b)+ - ((a-c)*(a-d)*v-(b-c)*(b-d)*u)*(c-d)*(c-x)*(d-x) - +((a-c)*(a-d)*(a-x)*v-(b-c)*(b-d)*(b-x)*u) - *(c-x+d-x)*(c-d)))/ - ((a-b)*(a-c)*(a-d)*(b-c)*(b-d)*(c-d)); - } }; } diff --git a/ql/math/matrix.cpp b/ql/math/matrix.cpp index 81f03b5a21b..23030d8402d 100644 --- a/ql/math/matrix.cpp +++ b/ql/math/matrix.cpp @@ -22,6 +22,8 @@ */ #include +#include + #if defined(QL_PATCH_MSVC) #pragma warning(push) #pragma warning(disable:4180) @@ -100,4 +102,506 @@ namespace QuantLib { return retVal; } + // definitions + + Matrix::Matrix() : data_((Real*)nullptr) {} + + Matrix::Matrix(Size rows, Size columns) + : data_(rows * columns > 0 ? new Real[rows * columns] : (Real*)nullptr), rows_(rows), + columns_(columns) {} + + Matrix::Matrix(Size rows, Size columns, Real value) + : data_(rows * columns > 0 ? new Real[rows * columns] : (Real*)nullptr), rows_(rows), + columns_(columns) { + std::fill(begin(),end(),value); + } + + Matrix::Matrix(const Matrix& from) + : data_(!from.empty() ? new Real[from.rows_ * from.columns_] : (Real*)nullptr), + rows_(from.rows_), columns_(from.columns_) { + #if defined(QL_PATCH_MSVC) && defined(QL_DEBUG) + if (!from.empty()) + #endif + std::copy(from.begin(),from.end(),begin()); + } + + Matrix::Matrix(Matrix&& from) noexcept + : data_((Real*)nullptr) { + swap(from); + } + + Matrix::Matrix(std::initializer_list> data) + : data_(data.size() == 0 || data.begin()->size() == 0 ? + (Real*)nullptr : new Real[data.size() * data.begin()->size()]), + rows_(data.size()), columns_(data.size() == 0 ? 0 : data.begin()->size()) { + Size i=0; + for (const auto& row : data) { + #if defined(QL_EXTRA_SAFETY_CHECKS) + QL_REQUIRE(row.size() == columns_, + "a matrix needs the same number of elements for each row"); + #endif + std::copy(row.begin(), row.end(), row_begin(i)); + ++i; + } + } + + Matrix& Matrix::operator=(const Matrix& from) { + // strong guarantee + Matrix temp(from); + swap(temp); + return *this; + } + + Matrix& Matrix::operator=(Matrix&& from) noexcept { + swap(from); + return *this; + } + + bool Matrix::operator==(const Matrix& to) const { + return rows_ == to.rows_ && columns_ == to.columns_ && + std::equal(begin(), end(), to.begin()); + } + + bool Matrix::operator!=(const Matrix& to) const { + return !this->operator==(to); + } + + void Matrix::swap(Matrix& from) noexcept { + data_.swap(from.data_); + std::swap(rows_, from.rows_); + std::swap(columns_, from.columns_); + } + + const Matrix& Matrix::operator+=(const Matrix& m) { + QL_REQUIRE(rows_ == m.rows_ && columns_ == m.columns_, + "matrices with different sizes (" << + m.rows_ << "x" << m.columns_ << ", " << + rows_ << "x" << columns_ << ") cannot be " + "added"); + std::transform(begin(), end(), m.begin(), begin(), std::plus<>()); + return *this; + } + + const Matrix& Matrix::operator-=(const Matrix& m) { + QL_REQUIRE(rows_ == m.rows_ && columns_ == m.columns_, + "matrices with different sizes (" << + m.rows_ << "x" << m.columns_ << ", " << + rows_ << "x" << columns_ << ") cannot be " + "subtracted"); + std::transform(begin(), end(), m.begin(), begin(), std::minus<>()); + return *this; + } + + const Matrix& Matrix::operator*=(Real x) { + std::transform(begin(), end(), begin(), [=](Real y) -> Real { return y * x; }); + return *this; + } + + const Matrix& Matrix::operator/=(Real x) { + std::transform(begin(),end(),begin(), [=](Real y) -> Real { return y / x; }); + return *this; + } + + Matrix::const_iterator Matrix::begin() const { + return data_.get(); + } + + Matrix::iterator Matrix::begin() { + return data_.get(); + } + + Matrix::const_iterator Matrix::end() const { + return data_.get()+rows_*columns_; + } + + Matrix::iterator Matrix::end() { + return data_.get()+rows_*columns_; + } + + Matrix::const_reverse_iterator Matrix::rbegin() const { + return const_reverse_iterator(end()); + } + + Matrix::reverse_iterator Matrix::rbegin() { + return reverse_iterator(end()); + } + + Matrix::const_reverse_iterator Matrix::rend() const { + return const_reverse_iterator(begin()); + } + + Matrix::reverse_iterator Matrix::rend() { + return reverse_iterator(begin()); + } + + Matrix::const_row_iterator + Matrix::row_begin(Size i) const { + #if defined(QL_EXTRA_SAFETY_CHECKS) + QL_REQUIRE(i(rows(), columns()); + Array tmp(arraySize); + for(Size i = 0; i < arraySize; i++) + tmp[i] = (*this)[i][i]; + return tmp; + } + + Real &Matrix::operator()(Size i, Size j) const { + return data_[i*columns()+j]; + } + + Size Matrix::rows() const { + return rows_; + } + + Size Matrix::columns() const { + return columns_; + } + + Size Matrix::size1() const { + return rows(); + } + + Size Matrix::size2() const { + return columns(); + } + + bool Matrix::empty() const { + return rows_ == 0 || columns_ == 0; + } + + Matrix operator+(const Matrix& m1, const Matrix& m2) { + QL_REQUIRE(m1.rows() == m2.rows() && + m1.columns() == m2.columns(), + "matrices with different sizes (" << + m1.rows() << "x" << m1.columns() << ", " << + m2.rows() << "x" << m2.columns() << ") cannot be " + "added"); + Matrix temp(m1.rows(),m1.columns()); + std::transform(m1.begin(), m1.end(), m2.begin(), temp.begin(), std::plus<>()); + return temp; + } + + Matrix operator+(const Matrix& m1, Matrix&& m2) { + QL_REQUIRE(m1.rows() == m2.rows() && + m1.columns() == m2.columns(), + "matrices with different sizes (" << + m1.rows() << "x" << m1.columns() << ", " << + m2.rows() << "x" << m2.columns() << ") cannot be " + "added"); + std::transform(m1.begin(), m1.end(), m2.begin(), m2.begin(), std::plus<>()); + return std::move(m2); + } + + Matrix operator+(Matrix&& m1, const Matrix& m2) { + QL_REQUIRE(m1.rows() == m2.rows() && + m1.columns() == m2.columns(), + "matrices with different sizes (" << + m1.rows() << "x" << m1.columns() << ", " << + m2.rows() << "x" << m2.columns() << ") cannot be " + "added"); + std::transform(m1.begin(), m1.end(), m2.begin(), m1.begin(), std::plus<>()); + return std::move(m1); + } + + Matrix operator+(Matrix&& m1, Matrix&& m2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved) + QL_REQUIRE(m1.rows() == m2.rows() && + m1.columns() == m2.columns(), + "matrices with different sizes (" << + m1.rows() << "x" << m1.columns() << ", " << + m2.rows() << "x" << m2.columns() << ") cannot be " + "added"); + std::transform(m1.begin(), m1.end(), m2.begin(), m1.begin(), std::plus<>()); + return std::move(m1); + } + + Matrix operator-(const Matrix& m1) { + Matrix temp(m1.rows(), m1.columns()); + std::transform(m1.begin(), m1.end(), temp.begin(), std::negate<>()); + return temp; + } + + Matrix operator-(Matrix&& m1) { + std::transform(m1.begin(), m1.end(), m1.begin(), std::negate<>()); + return std::move(m1); + } + + Matrix operator-(const Matrix& m1, const Matrix& m2) { + QL_REQUIRE(m1.rows() == m2.rows() && + m1.columns() == m2.columns(), + "matrices with different sizes (" << + m1.rows() << "x" << m1.columns() << ", " << + m2.rows() << "x" << m2.columns() << ") cannot be " + "subtracted"); + Matrix temp(m1.rows(),m1.columns()); + std::transform(m1.begin(), m1.end(), m2.begin(), temp.begin(), std::minus<>()); + return temp; + } + + Matrix operator-(const Matrix& m1, Matrix&& m2) { + QL_REQUIRE(m1.rows() == m2.rows() && + m1.columns() == m2.columns(), + "matrices with different sizes (" << + m1.rows() << "x" << m1.columns() << ", " << + m2.rows() << "x" << m2.columns() << ") cannot be " + "subtracted"); + std::transform(m1.begin(), m1.end(), m2.begin(), m2.begin(), std::minus<>()); + return std::move(m2); + } + + Matrix operator-(Matrix&& m1, const Matrix& m2) { + QL_REQUIRE(m1.rows() == m2.rows() && + m1.columns() == m2.columns(), + "matrices with different sizes (" << + m1.rows() << "x" << m1.columns() << ", " << + m2.rows() << "x" << m2.columns() << ") cannot be " + "subtracted"); + std::transform(m1.begin(), m1.end(), m2.begin(), m1.begin(), std::minus<>()); + return std::move(m1); + } + + Matrix operator-(Matrix&& m1, Matrix&& m2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved) + QL_REQUIRE(m1.rows() == m2.rows() && + m1.columns() == m2.columns(), + "matrices with different sizes (" << + m1.rows() << "x" << m1.columns() << ", " << + m2.rows() << "x" << m2.columns() << ") cannot be " + "subtracted"); + std::transform(m1.begin(), m1.end(), m2.begin(), m1.begin(), std::minus<>()); + return std::move(m1); + } + + Matrix operator*(const Matrix& m, Real x) { + Matrix temp(m.rows(),m.columns()); + std::transform(m.begin(), m.end(), temp.begin(), [=](Real y) -> Real { return y * x; }); + return temp; + } + + Matrix operator*(Matrix&& m, Real x) { + std::transform(m.begin(), m.end(), m.begin(), [=](Real y) -> Real { return y * x; }); + return std::move(m); + } + + Matrix operator*(Real x, const Matrix& m) { + Matrix temp(m.rows(),m.columns()); + std::transform(m.begin(), m.end(), temp.begin(), [=](Real y) -> Real { return x * y; }); + return temp; + } + + Matrix operator*(Real x, Matrix&& m) { + std::transform(m.begin(), m.end(), m.begin(), [=](Real y) -> Real { return x * y; }); + return std::move(m); + } + + Matrix operator/(const Matrix& m, Real x) { + Matrix temp(m.rows(),m.columns()); + std::transform(m.begin(), m.end(), temp.begin(), [=](Real y) -> Real { return y / x; }); + return temp; + } + + Matrix operator/(Matrix&& m, Real x) { + std::transform(m.begin(), m.end(), m.begin(), [=](Real y) -> Real { return y / x; }); + return std::move(m); + } + + Array operator*(const Array& v, const Matrix& m) { + QL_REQUIRE(v.size() == m.rows(), + "vectors and matrices with different sizes (" + << v.size() << ", " << m.rows() << "x" << m.columns() << + ") cannot be multiplied"); + Array result(m.columns()); + for (Size i=0; i 0 ? new Real[rows * columns] : (Real*)nullptr), rows_(rows), - columns_(columns) {} - - inline Matrix::Matrix(Size rows, Size columns, Real value) - : data_(rows * columns > 0 ? new Real[rows * columns] : (Real*)nullptr), rows_(rows), - columns_(columns) { - std::fill(begin(),end(),value); - } - - template - inline Matrix::Matrix(Size rows, Size columns, Iterator begin, Iterator end) - : data_(rows * columns > 0 ? new Real[rows * columns] : (Real*)nullptr), rows_(rows), - columns_(columns) { - std::copy(begin, end, this->begin()); - } - - inline Matrix::Matrix(const Matrix& from) - : data_(!from.empty() ? new Real[from.rows_ * from.columns_] : (Real*)nullptr), - rows_(from.rows_), columns_(from.columns_) { - #if defined(QL_PATCH_MSVC) && defined(QL_DEBUG) - if (!from.empty()) - #endif - std::copy(from.begin(),from.end(),begin()); - } - - inline Matrix::Matrix(Matrix&& from) noexcept - : data_((Real*)nullptr) { - swap(from); - } - - inline Matrix::Matrix(std::initializer_list> data) - : data_(data.size() == 0 || data.begin()->size() == 0 ? - (Real*)nullptr : new Real[data.size() * data.begin()->size()]), - rows_(data.size()), columns_(data.size() == 0 ? 0 : data.begin()->size()) { - Size i=0; - for (const auto& row : data) { - #if defined(QL_EXTRA_SAFETY_CHECKS) - QL_REQUIRE(row.size() == columns_, - "a matrix needs the same number of elements for each row"); - #endif - std::copy(row.begin(), row.end(), row_begin(i)); - ++i; - } - } - - inline Matrix& Matrix::operator=(const Matrix& from) { - // strong guarantee - Matrix temp(from); - swap(temp); - return *this; - } - - inline Matrix& Matrix::operator=(Matrix&& from) noexcept { - swap(from); - return *this; - } - - inline bool Matrix::operator==(const Matrix& to) const { - return rows_ == to.rows_ && columns_ == to.columns_ && - std::equal(begin(), end(), to.begin()); - } - - inline bool Matrix::operator!=(const Matrix& to) const { - return !this->operator==(to); - } - - inline void Matrix::swap(Matrix& from) noexcept { - data_.swap(from.data_); - std::swap(rows_, from.rows_); - std::swap(columns_, from.columns_); - } - - inline const Matrix& Matrix::operator+=(const Matrix& m) { - QL_REQUIRE(rows_ == m.rows_ && columns_ == m.columns_, - "matrices with different sizes (" << - m.rows_ << "x" << m.columns_ << ", " << - rows_ << "x" << columns_ << ") cannot be " - "added"); - std::transform(begin(), end(), m.begin(), begin(), std::plus<>()); - return *this; - } - - inline const Matrix& Matrix::operator-=(const Matrix& m) { - QL_REQUIRE(rows_ == m.rows_ && columns_ == m.columns_, - "matrices with different sizes (" << - m.rows_ << "x" << m.columns_ << ", " << - rows_ << "x" << columns_ << ") cannot be " - "subtracted"); - std::transform(begin(), end(), m.begin(), begin(), std::minus<>()); - return *this; - } - - inline const Matrix& Matrix::operator*=(Real x) { - std::transform(begin(), end(), begin(), [=](Real y) -> Real { return y * x; }); - return *this; - } - - inline const Matrix& Matrix::operator/=(Real x) { - std::transform(begin(),end(),begin(), [=](Real y) -> Real { return y / x; }); - return *this; - } - - inline Matrix::const_iterator Matrix::begin() const { - return data_.get(); - } - - inline Matrix::iterator Matrix::begin() { - return data_.get(); - } - - inline Matrix::const_iterator Matrix::end() const { - return data_.get()+rows_*columns_; - } - - inline Matrix::iterator Matrix::end() { - return data_.get()+rows_*columns_; - } - - inline Matrix::const_reverse_iterator Matrix::rbegin() const { - return const_reverse_iterator(end()); - } - - inline Matrix::reverse_iterator Matrix::rbegin() { - return reverse_iterator(end()); - } - - inline Matrix::const_reverse_iterator Matrix::rend() const { - return const_reverse_iterator(begin()); - } - - inline Matrix::reverse_iterator Matrix::rend() { - return reverse_iterator(begin()); - } - - inline Matrix::const_row_iterator - Matrix::row_begin(Size i) const { - #if defined(QL_EXTRA_SAFETY_CHECKS) - QL_REQUIRE(i(rows(), columns()); - Array tmp(arraySize); - for(Size i = 0; i < arraySize; i++) - tmp[i] = (*this)[i][i]; - return tmp; - } - - inline Real &Matrix::operator()(Size i, Size j) const { - return data_[i*columns()+j]; - } - - inline Size Matrix::rows() const { - return rows_; - } - - inline Size Matrix::columns() const { - return columns_; - } - - inline Size Matrix::size1() const { - return rows(); - } - - inline Size Matrix::size2() const { - return columns(); - } - - inline bool Matrix::empty() const { - return rows_ == 0 || columns_ == 0; - } - - inline Matrix operator+(const Matrix& m1, const Matrix& m2) { - QL_REQUIRE(m1.rows() == m2.rows() && - m1.columns() == m2.columns(), - "matrices with different sizes (" << - m1.rows() << "x" << m1.columns() << ", " << - m2.rows() << "x" << m2.columns() << ") cannot be " - "added"); - Matrix temp(m1.rows(),m1.columns()); - std::transform(m1.begin(), m1.end(), m2.begin(), temp.begin(), std::plus<>()); - return temp; - } - - inline Matrix operator+(const Matrix& m1, Matrix&& m2) { - QL_REQUIRE(m1.rows() == m2.rows() && - m1.columns() == m2.columns(), - "matrices with different sizes (" << - m1.rows() << "x" << m1.columns() << ", " << - m2.rows() << "x" << m2.columns() << ") cannot be " - "added"); - std::transform(m1.begin(), m1.end(), m2.begin(), m2.begin(), std::plus<>()); - return std::move(m2); - } - - inline Matrix operator+(Matrix&& m1, const Matrix& m2) { - QL_REQUIRE(m1.rows() == m2.rows() && - m1.columns() == m2.columns(), - "matrices with different sizes (" << - m1.rows() << "x" << m1.columns() << ", " << - m2.rows() << "x" << m2.columns() << ") cannot be " - "added"); - std::transform(m1.begin(), m1.end(), m2.begin(), m1.begin(), std::plus<>()); - return std::move(m1); - } - - inline Matrix operator+(Matrix&& m1, Matrix&& m2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved) - QL_REQUIRE(m1.rows() == m2.rows() && - m1.columns() == m2.columns(), - "matrices with different sizes (" << - m1.rows() << "x" << m1.columns() << ", " << - m2.rows() << "x" << m2.columns() << ") cannot be " - "added"); - std::transform(m1.begin(), m1.end(), m2.begin(), m1.begin(), std::plus<>()); - return std::move(m1); - } - - inline Matrix operator-(const Matrix& m1) { - Matrix temp(m1.rows(), m1.columns()); - std::transform(m1.begin(), m1.end(), temp.begin(), std::negate<>()); - return temp; - } - - inline Matrix operator-(Matrix&& m1) { - std::transform(m1.begin(), m1.end(), m1.begin(), std::negate<>()); - return std::move(m1); - } - - inline Matrix operator-(const Matrix& m1, const Matrix& m2) { - QL_REQUIRE(m1.rows() == m2.rows() && - m1.columns() == m2.columns(), - "matrices with different sizes (" << - m1.rows() << "x" << m1.columns() << ", " << - m2.rows() << "x" << m2.columns() << ") cannot be " - "subtracted"); - Matrix temp(m1.rows(),m1.columns()); - std::transform(m1.begin(), m1.end(), m2.begin(), temp.begin(), std::minus<>()); - return temp; - } - - inline Matrix operator-(const Matrix& m1, Matrix&& m2) { - QL_REQUIRE(m1.rows() == m2.rows() && - m1.columns() == m2.columns(), - "matrices with different sizes (" << - m1.rows() << "x" << m1.columns() << ", " << - m2.rows() << "x" << m2.columns() << ") cannot be " - "subtracted"); - std::transform(m1.begin(), m1.end(), m2.begin(), m2.begin(), std::minus<>()); - return std::move(m2); - } - - inline Matrix operator-(Matrix&& m1, const Matrix& m2) { - QL_REQUIRE(m1.rows() == m2.rows() && - m1.columns() == m2.columns(), - "matrices with different sizes (" << - m1.rows() << "x" << m1.columns() << ", " << - m2.rows() << "x" << m2.columns() << ") cannot be " - "subtracted"); - std::transform(m1.begin(), m1.end(), m2.begin(), m1.begin(), std::minus<>()); - return std::move(m1); - } - - inline Matrix operator-(Matrix&& m1, Matrix&& m2) { // NOLINT(cppcoreguidelines-rvalue-reference-param-not-moved) - QL_REQUIRE(m1.rows() == m2.rows() && - m1.columns() == m2.columns(), - "matrices with different sizes (" << - m1.rows() << "x" << m1.columns() << ", " << - m2.rows() << "x" << m2.columns() << ") cannot be " - "subtracted"); - std::transform(m1.begin(), m1.end(), m2.begin(), m1.begin(), std::minus<>()); - return std::move(m1); - } - - inline Matrix operator*(const Matrix& m, Real x) { - Matrix temp(m.rows(),m.columns()); - std::transform(m.begin(), m.end(), temp.begin(), [=](Real y) -> Real { return y * x; }); - return temp; - } - - inline Matrix operator*(Matrix&& m, Real x) { - std::transform(m.begin(), m.end(), m.begin(), [=](Real y) -> Real { return y * x; }); - return std::move(m); - } - - inline Matrix operator*(Real x, const Matrix& m) { - Matrix temp(m.rows(),m.columns()); - std::transform(m.begin(), m.end(), temp.begin(), [=](Real y) -> Real { return x * y; }); - return temp; - } - - inline Matrix operator*(Real x, Matrix&& m) { - std::transform(m.begin(), m.end(), m.begin(), [=](Real y) -> Real { return x * y; }); - return std::move(m); - } - - inline Matrix operator/(const Matrix& m, Real x) { - Matrix temp(m.rows(),m.columns()); - std::transform(m.begin(), m.end(), temp.begin(), [=](Real y) -> Real { return y / x; }); - return temp; - } - - inline Matrix operator/(Matrix&& m, Real x) { - std::transform(m.begin(), m.end(), m.begin(), [=](Real y) -> Real { return y / x; }); - return std::move(m); - } - - inline Array operator*(const Array& v, const Matrix& m) { - QL_REQUIRE(v.size() == m.rows(), - "vectors and matrices with different sizes (" - << v.size() << ", " << m.rows() << "x" << m.columns() << - ") cannot be multiplied"); - Array result(m.columns()); - for (Size i=0; i - inline Matrix outerProduct(Iterator1 v1begin, Iterator1 v1end, Iterator2 v2begin, Iterator2 v2end) { + Matrix outerProduct(Iterator1 v1begin, Iterator1 v1end, Iterator2 v2begin, Iterator2 v2end) { Size size1 = std::distance(v1begin, v1end); QL_REQUIRE(size1>0, "null first vector"); @@ -729,19 +236,12 @@ namespace QuantLib { return result; } - inline void swap(Matrix& m1, Matrix& m2) noexcept { - m1.swap(m2); - } - - inline std::ostream& operator<<(std::ostream& out, const Matrix& m) { - std::streamsize width = out.width(); - for (Size i=0; i + Matrix::Matrix(Size rows, Size columns, Iterator begin, Iterator end) + : data_(rows * columns > 0 ? new Real[rows * columns] : (Real*)nullptr), rows_(rows), + columns_(columns) { + std::copy(begin, end, this->begin()); } } diff --git a/ql/math/matrixutilities/bicgstab.cpp b/ql/math/matrixutilities/bicgstab.cpp index 41bf4b40354..1cd97409e24 100644 --- a/ql/math/matrixutilities/bicgstab.cpp +++ b/ql/math/matrixutilities/bicgstab.cpp @@ -23,6 +23,7 @@ */ +#include #include #include diff --git a/ql/math/matrixutilities/sparsematrix.hpp b/ql/math/matrixutilities/sparsematrix.hpp index d228332b865..be819af56ab 100644 --- a/ql/math/matrixutilities/sparsematrix.hpp +++ b/ql/math/matrixutilities/sparsematrix.hpp @@ -25,6 +25,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. #define quantlib_sparse_matrix_hpp #include +#include #include #if defined(QL_PATCH_MSVC) diff --git a/ql/math/optimization/constraint.hpp b/ql/math/optimization/constraint.hpp index ddbaef85da3..93987f39246 100644 --- a/ql/math/optimization/constraint.hpp +++ b/ql/math/optimization/constraint.hpp @@ -25,7 +25,9 @@ #ifndef quantlib_optimization_constraint_h #define quantlib_optimization_constraint_h +#include #include +#include #include #include diff --git a/ql/math/optimization/projection.cpp b/ql/math/optimization/projection.cpp index d2799f8421f..5a3f5005f9f 100644 --- a/ql/math/optimization/projection.cpp +++ b/ql/math/optimization/projection.cpp @@ -19,6 +19,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include diff --git a/ql/math/statistics/histogram.cpp b/ql/math/statistics/histogram.cpp index 0797415aeb6..cafc81d269d 100644 --- a/ql/math/statistics/histogram.cpp +++ b/ql/math/statistics/histogram.cpp @@ -17,6 +17,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include #include diff --git a/ql/math/statistics/incrementalstatistics.cpp b/ql/math/statistics/incrementalstatistics.cpp index 6495f6cbffe..700ff87603e 100644 --- a/ql/math/statistics/incrementalstatistics.cpp +++ b/ql/math/statistics/incrementalstatistics.cpp @@ -19,29 +19,70 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include -#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include namespace QuantLib { + class IncrementalStatistics::accumulator_set + : public boost::accumulators::accumulator_set< + Real, + boost::accumulators::stats, + Real> {}; + + class IncrementalStatistics::downside_accumulator_set + : public boost::accumulators::accumulator_set< + Real, + boost::accumulators::stats, + boost::accumulators::tag::sum_of_weights>, + Real> {}; + IncrementalStatistics::IncrementalStatistics() { + acc_ = nullptr; + downsideAcc_ = nullptr; reset(); } + IncrementalStatistics::~IncrementalStatistics() { + if (acc_) + delete acc_; + if (downsideAcc_) + delete downsideAcc_; + } + Size IncrementalStatistics::samples() const { return boost::accumulators::extract_result< - boost::accumulators::tag::count>(acc_); + boost::accumulators::tag::count>(*acc_); } Real IncrementalStatistics::weightSum() const { return boost::accumulators::extract_result< - boost::accumulators::tag::sum_of_weights>(acc_); + boost::accumulators::tag::sum_of_weights>(*acc_); } Real IncrementalStatistics::mean() const { QL_REQUIRE(weightSum() > 0.0, "sampleWeight_= 0, unsufficient"); return boost::accumulators::extract_result< - boost::accumulators::tag::weighted_mean>(acc_); + boost::accumulators::tag::weighted_mean>(*acc_); } Real IncrementalStatistics::variance() const { @@ -50,7 +91,7 @@ namespace QuantLib { Real n = static_cast(samples()); return n / (n - 1.0) * boost::accumulators::extract_result< - boost::accumulators::tag::weighted_variance>(acc_); + boost::accumulators::tag::weighted_variance>(*acc_); } Real IncrementalStatistics::standardDeviation() const { @@ -68,20 +109,20 @@ namespace QuantLib { Real r2 = (n - 1.0) / (n - 2.0); return std::sqrt(r1 * r2) * boost::accumulators::extract_result< - boost::accumulators::tag::weighted_skewness>(acc_); + boost::accumulators::tag::weighted_skewness>(*acc_); } Real IncrementalStatistics::kurtosis() const { QL_REQUIRE(samples() > 3, "sample number <= 3, unsufficient"); boost::accumulators::extract_result< - boost::accumulators::tag::weighted_kurtosis>(acc_); + boost::accumulators::tag::weighted_kurtosis>(*acc_); Real n = static_cast(samples()); Real r1 = (n - 1.0) / (n - 2.0); Real r2 = (n + 1.0) / (n - 3.0); Real r3 = (n - 1.0) / (n - 3.0); return ((3.0 + boost::accumulators::extract_result< - boost::accumulators::tag::weighted_kurtosis>(acc_)) * + boost::accumulators::tag::weighted_kurtosis>(*acc_)) * r2 - 3.0 * r3) * r1; @@ -90,23 +131,23 @@ namespace QuantLib { Real IncrementalStatistics::min() const { QL_REQUIRE(samples() > 0, "empty sample set"); return boost::accumulators::extract_result< - boost::accumulators::tag::min>(acc_); + boost::accumulators::tag::min>(*acc_); } Real IncrementalStatistics::max() const { QL_REQUIRE(samples() > 0, "empty sample set"); return boost::accumulators::extract_result< - boost::accumulators::tag::max>(acc_); + boost::accumulators::tag::max>(*acc_); } Size IncrementalStatistics::downsideSamples() const { return boost::accumulators::extract_result< - boost::accumulators::tag::count>(downsideAcc_); + boost::accumulators::tag::count>(*downsideAcc_); } Real IncrementalStatistics::downsideWeightSum() const { return boost::accumulators::extract_result< - boost::accumulators::tag::sum_of_weights>(downsideAcc_); + boost::accumulators::tag::sum_of_weights>(*downsideAcc_); } Real IncrementalStatistics::downsideVariance() const { @@ -116,7 +157,7 @@ namespace QuantLib { Real r1 = n / (n - 1.0); return r1 * boost::accumulators::extract_result< - boost::accumulators::tag::moment<2> >(downsideAcc_); + boost::accumulators::tag::moment<2> >(*downsideAcc_); } Real IncrementalStatistics::downsideDeviation() const { @@ -126,14 +167,18 @@ namespace QuantLib { void IncrementalStatistics::add(Real value, Real valueWeight) { QL_REQUIRE(valueWeight >= 0.0, "negative weight (" << valueWeight << ") not allowed"); - acc_(value, boost::accumulators::weight = valueWeight); + (*acc_)(value, boost::accumulators::weight = valueWeight); if(value < 0.0) - downsideAcc_(value, boost::accumulators::weight = valueWeight); + (*downsideAcc_)(value, boost::accumulators::weight = valueWeight); } void IncrementalStatistics::reset() { - acc_ = accumulator_set(); - downsideAcc_ = downside_accumulator_set(); + if (acc_) + delete acc_; + if (downsideAcc_) + delete downsideAcc_; + acc_ = new accumulator_set; + downsideAcc_ = new downside_accumulator_set; } } diff --git a/ql/math/statistics/incrementalstatistics.hpp b/ql/math/statistics/incrementalstatistics.hpp index d633dbdc3b0..3355bc5dd8b 100644 --- a/ql/math/statistics/incrementalstatistics.hpp +++ b/ql/math/statistics/incrementalstatistics.hpp @@ -29,18 +29,6 @@ #define quantlib_incremental_statistics_hpp #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include namespace QuantLib { @@ -54,6 +42,7 @@ namespace QuantLib { public: typedef Real value_type; IncrementalStatistics(); + ~IncrementalStatistics(); //! \name Inspectors //@{ //! number of samples collected @@ -149,25 +138,10 @@ namespace QuantLib { void reset(); //@} private: - typedef boost::accumulators::accumulator_set< - Real, - boost::accumulators::stats< - boost::accumulators::tag::count, boost::accumulators::tag::min, - boost::accumulators::tag::max, - boost::accumulators::tag::weighted_mean, - boost::accumulators::tag::weighted_variance, - boost::accumulators::tag::weighted_skewness, - boost::accumulators::tag::weighted_kurtosis, - boost::accumulators::tag::sum_of_weights>, - Real> accumulator_set; - accumulator_set acc_; - typedef boost::accumulators::accumulator_set< - Real, boost::accumulators::stats< - boost::accumulators::tag::count, - boost::accumulators::tag::weighted_moment<2>, - boost::accumulators::tag::sum_of_weights>, - Real> downside_accumulator_set; - downside_accumulator_set downsideAcc_; + class accumulator_set; + class downside_accumulator_set; + accumulator_set* acc_; + downside_accumulator_set* downsideAcc_; }; } diff --git a/ql/methods/finitedifferences/meshers/fdmmeshercomposite.cpp b/ql/methods/finitedifferences/meshers/fdmmeshercomposite.cpp index ece9d144505..c3da6840cc9 100644 --- a/ql/methods/finitedifferences/meshers/fdmmeshercomposite.cpp +++ b/ql/methods/finitedifferences/meshers/fdmmeshercomposite.cpp @@ -19,6 +19,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include diff --git a/ql/methods/finitedifferences/meshers/uniformgridmesher.cpp b/ql/methods/finitedifferences/meshers/uniformgridmesher.cpp index e5a47e988d2..ae16a63b724 100644 --- a/ql/methods/finitedifferences/meshers/uniformgridmesher.cpp +++ b/ql/methods/finitedifferences/meshers/uniformgridmesher.cpp @@ -19,6 +19,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include diff --git a/ql/methods/finitedifferences/operators/fdmsquarerootfwdop.hpp b/ql/methods/finitedifferences/operators/fdmsquarerootfwdop.hpp index db0fadfbc1b..9d7e8d730ce 100644 --- a/ql/methods/finitedifferences/operators/fdmsquarerootfwdop.hpp +++ b/ql/methods/finitedifferences/operators/fdmsquarerootfwdop.hpp @@ -26,6 +26,7 @@ #define quantlib_fdm_square_root_fwd_op_hpp #include +#include namespace QuantLib { class FdmMesher; diff --git a/ql/methods/finitedifferences/operators/ninepointlinearop.hpp b/ql/methods/finitedifferences/operators/ninepointlinearop.hpp index 8da0ae1bbac..a806209a0d5 100644 --- a/ql/methods/finitedifferences/operators/ninepointlinearop.hpp +++ b/ql/methods/finitedifferences/operators/ninepointlinearop.hpp @@ -28,6 +28,7 @@ #include #include +#include #include namespace QuantLib { diff --git a/ql/methods/finitedifferences/operators/numericaldifferentiation.cpp b/ql/methods/finitedifferences/operators/numericaldifferentiation.cpp index 44b8a25b523..536b0587f0f 100644 --- a/ql/methods/finitedifferences/operators/numericaldifferentiation.cpp +++ b/ql/methods/finitedifferences/operators/numericaldifferentiation.cpp @@ -19,6 +19,7 @@ /*! \file numericaldifferentiation.cpp */ +#include #include #pragma push_macro("BOOST_DISABLE_ASSERTS") diff --git a/ql/methods/finitedifferences/operators/triplebandlinearop.hpp b/ql/methods/finitedifferences/operators/triplebandlinearop.hpp index 25def57fcb6..3e99a8e1aa4 100644 --- a/ql/methods/finitedifferences/operators/triplebandlinearop.hpp +++ b/ql/methods/finitedifferences/operators/triplebandlinearop.hpp @@ -27,7 +27,9 @@ #ifndef quantlib_triple_band_linear_op_hpp #define quantlib_triple_band_linear_op_hpp +#include #include +#include #include namespace QuantLib { diff --git a/ql/methods/finitedifferences/tridiagonaloperator.hpp b/ql/methods/finitedifferences/tridiagonaloperator.hpp index 36d44267d1b..8b0e7d13847 100644 --- a/ql/methods/finitedifferences/tridiagonaloperator.hpp +++ b/ql/methods/finitedifferences/tridiagonaloperator.hpp @@ -26,6 +26,7 @@ #ifndef quantlib_tridiagonal_operator_hpp #define quantlib_tridiagonal_operator_hpp +#include #include #include #include diff --git a/ql/models/marketmodels/accountingengine.hpp b/ql/models/marketmodels/accountingengine.hpp index d02dd1d77ba..dd7d770050e 100644 --- a/ql/models/marketmodels/accountingengine.hpp +++ b/ql/models/marketmodels/accountingengine.hpp @@ -27,6 +27,7 @@ #include #include +#include #include #include #include diff --git a/ql/models/marketmodels/callability/swapratetrigger.cpp b/ql/models/marketmodels/callability/swapratetrigger.cpp index 102c7b112fd..9356e3c660d 100644 --- a/ql/models/marketmodels/callability/swapratetrigger.cpp +++ b/ql/models/marketmodels/callability/swapratetrigger.cpp @@ -19,6 +19,7 @@ #include #include +#include #include namespace QuantLib { diff --git a/ql/models/marketmodels/callability/upperboundengine.hpp b/ql/models/marketmodels/callability/upperboundengine.hpp index 7dae8ba3182..daf66fed944 100644 --- a/ql/models/marketmodels/callability/upperboundengine.hpp +++ b/ql/models/marketmodels/callability/upperboundengine.hpp @@ -26,6 +26,7 @@ #include #include #include +#include #include #include diff --git a/ql/models/marketmodels/correlations/cotswapfromfwdcorrelation.hpp b/ql/models/marketmodels/correlations/cotswapfromfwdcorrelation.hpp index 1f7adbf8f88..baa67150b05 100644 --- a/ql/models/marketmodels/correlations/cotswapfromfwdcorrelation.hpp +++ b/ql/models/marketmodels/correlations/cotswapfromfwdcorrelation.hpp @@ -22,6 +22,7 @@ #define quantlib_fra_time_dep_corr_struct_hpp #include +#include #include namespace QuantLib { diff --git a/ql/models/marketmodels/curvestate.cpp b/ql/models/marketmodels/curvestate.cpp index e311e659158..e39929fa635 100644 --- a/ql/models/marketmodels/curvestate.cpp +++ b/ql/models/marketmodels/curvestate.cpp @@ -18,6 +18,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include diff --git a/ql/models/marketmodels/curvestates/cmswapcurvestate.cpp b/ql/models/marketmodels/curvestates/cmswapcurvestate.cpp index 32c309a2089..5eb0a0f3107 100644 --- a/ql/models/marketmodels/curvestates/cmswapcurvestate.cpp +++ b/ql/models/marketmodels/curvestates/cmswapcurvestate.cpp @@ -19,6 +19,7 @@ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include namespace QuantLib { diff --git a/ql/models/marketmodels/curvestates/coterminalswapcurvestate.cpp b/ql/models/marketmodels/curvestates/coterminalswapcurvestate.cpp index 907b55fec39..57afae46151 100644 --- a/ql/models/marketmodels/curvestates/coterminalswapcurvestate.cpp +++ b/ql/models/marketmodels/curvestates/coterminalswapcurvestate.cpp @@ -19,6 +19,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include namespace QuantLib { diff --git a/ql/models/marketmodels/curvestates/lmmcurvestate.cpp b/ql/models/marketmodels/curvestates/lmmcurvestate.cpp index b85cad843df..1b3c4f52fe9 100644 --- a/ql/models/marketmodels/curvestates/lmmcurvestate.cpp +++ b/ql/models/marketmodels/curvestates/lmmcurvestate.cpp @@ -21,6 +21,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include namespace QuantLib { diff --git a/ql/models/marketmodels/evolvers/lognormalcmswapratepc.hpp b/ql/models/marketmodels/evolvers/lognormalcmswapratepc.hpp index debf3811173..226adf480df 100644 --- a/ql/models/marketmodels/evolvers/lognormalcmswapratepc.hpp +++ b/ql/models/marketmodels/evolvers/lognormalcmswapratepc.hpp @@ -25,6 +25,7 @@ #include #include #include +#include namespace QuantLib { diff --git a/ql/models/marketmodels/evolvers/lognormalcotswapratepc.hpp b/ql/models/marketmodels/evolvers/lognormalcotswapratepc.hpp index e3804b99d50..8e1f77dc4cb 100644 --- a/ql/models/marketmodels/evolvers/lognormalcotswapratepc.hpp +++ b/ql/models/marketmodels/evolvers/lognormalcotswapratepc.hpp @@ -24,6 +24,7 @@ #include #include #include +#include namespace QuantLib { diff --git a/ql/models/marketmodels/evolvers/lognormalfwdrateballand.hpp b/ql/models/marketmodels/evolvers/lognormalfwdrateballand.hpp index 53e70646695..8de911c52f7 100644 --- a/ql/models/marketmodels/evolvers/lognormalfwdrateballand.hpp +++ b/ql/models/marketmodels/evolvers/lognormalfwdrateballand.hpp @@ -22,6 +22,7 @@ #include #include +#include namespace QuantLib { diff --git a/ql/models/marketmodels/evolvers/lognormalfwdrateeuler.hpp b/ql/models/marketmodels/evolvers/lognormalfwdrateeuler.hpp index d3fdd3df39c..3bf4d65a299 100644 --- a/ql/models/marketmodels/evolvers/lognormalfwdrateeuler.hpp +++ b/ql/models/marketmodels/evolvers/lognormalfwdrateeuler.hpp @@ -23,6 +23,7 @@ #include #include +#include namespace QuantLib { diff --git a/ql/models/marketmodels/evolvers/lognormalfwdrateeulerconstrained.hpp b/ql/models/marketmodels/evolvers/lognormalfwdrateeulerconstrained.hpp index b70b15c01ca..71aa155cfb5 100644 --- a/ql/models/marketmodels/evolvers/lognormalfwdrateeulerconstrained.hpp +++ b/ql/models/marketmodels/evolvers/lognormalfwdrateeulerconstrained.hpp @@ -25,6 +25,7 @@ #include #include #include +#include #include namespace QuantLib { diff --git a/ql/models/marketmodels/evolvers/lognormalfwdrateiballand.hpp b/ql/models/marketmodels/evolvers/lognormalfwdrateiballand.hpp index 154eade4e9a..4dbca455e0f 100644 --- a/ql/models/marketmodels/evolvers/lognormalfwdrateiballand.hpp +++ b/ql/models/marketmodels/evolvers/lognormalfwdrateiballand.hpp @@ -22,6 +22,7 @@ #include #include +#include namespace QuantLib { diff --git a/ql/models/marketmodels/evolvers/lognormalfwdrateipc.hpp b/ql/models/marketmodels/evolvers/lognormalfwdrateipc.hpp index 4b266da245c..c339bf1b354 100644 --- a/ql/models/marketmodels/evolvers/lognormalfwdrateipc.hpp +++ b/ql/models/marketmodels/evolvers/lognormalfwdrateipc.hpp @@ -23,6 +23,7 @@ #include #include +#include namespace QuantLib { diff --git a/ql/models/marketmodels/evolvers/lognormalfwdratepc.hpp b/ql/models/marketmodels/evolvers/lognormalfwdratepc.hpp index 253f5cf2e17..971849fca12 100644 --- a/ql/models/marketmodels/evolvers/lognormalfwdratepc.hpp +++ b/ql/models/marketmodels/evolvers/lognormalfwdratepc.hpp @@ -24,6 +24,7 @@ #include #include #include +#include namespace QuantLib { diff --git a/ql/models/marketmodels/evolvers/normalfwdratepc.hpp b/ql/models/marketmodels/evolvers/normalfwdratepc.hpp index 2c4b8505675..71d04975527 100644 --- a/ql/models/marketmodels/evolvers/normalfwdratepc.hpp +++ b/ql/models/marketmodels/evolvers/normalfwdratepc.hpp @@ -24,6 +24,7 @@ #include #include #include +#include namespace QuantLib { diff --git a/ql/models/marketmodels/evolvers/svddfwdratepc.hpp b/ql/models/marketmodels/evolvers/svddfwdratepc.hpp index 84f9e70a531..cd6bc7d7c52 100644 --- a/ql/models/marketmodels/evolvers/svddfwdratepc.hpp +++ b/ql/models/marketmodels/evolvers/svddfwdratepc.hpp @@ -23,6 +23,7 @@ #include #include #include +#include #include namespace QuantLib { diff --git a/ql/models/marketmodels/historicalratesanalysis.hpp b/ql/models/marketmodels/historicalratesanalysis.hpp index f946711bbd6..826a89c7c50 100644 --- a/ql/models/marketmodels/historicalratesanalysis.hpp +++ b/ql/models/marketmodels/historicalratesanalysis.hpp @@ -25,6 +25,7 @@ #define quantlib_historical_rates_analysis_hpp #include +#include #include namespace QuantLib { diff --git a/ql/models/marketmodels/models/ctsmmcapletcalibration.cpp b/ql/models/marketmodels/models/ctsmmcapletcalibration.cpp index 875dd34efb0..ea1d01afe2b 100644 --- a/ql/models/marketmodels/models/ctsmmcapletcalibration.cpp +++ b/ql/models/marketmodels/models/ctsmmcapletcalibration.cpp @@ -27,6 +27,7 @@ #include #include #include +#include #include namespace QuantLib { diff --git a/ql/models/marketmodels/pathwiseaccountingengine.hpp b/ql/models/marketmodels/pathwiseaccountingengine.hpp index 0cebc948a09..5d0f2ecc32b 100644 --- a/ql/models/marketmodels/pathwiseaccountingengine.hpp +++ b/ql/models/marketmodels/pathwiseaccountingengine.hpp @@ -22,6 +22,7 @@ #ifndef quantlib_pathwise_accounting_engine_hpp #define quantlib_pathwise_accounting_engine_hpp +#include #include #include #include diff --git a/ql/models/marketmodels/products/multistep/multistepinversefloater.cpp b/ql/models/marketmodels/products/multistep/multistepinversefloater.cpp index e940293b2ee..ca091dfa4b0 100644 --- a/ql/models/marketmodels/products/multistep/multistepinversefloater.cpp +++ b/ql/models/marketmodels/products/multistep/multistepinversefloater.cpp @@ -17,6 +17,7 @@ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include #include diff --git a/ql/models/marketmodels/products/multistep/multisteptarn.cpp b/ql/models/marketmodels/products/multistep/multisteptarn.cpp index b6727b4689c..f1bdb415fe2 100644 --- a/ql/models/marketmodels/products/multistep/multisteptarn.cpp +++ b/ql/models/marketmodels/products/multistep/multisteptarn.cpp @@ -16,6 +16,8 @@ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the license for more details. */ + +#include #include #include #include diff --git a/ql/models/marketmodels/products/pathwise/pathwiseproductcaplet.cpp b/ql/models/marketmodels/products/pathwise/pathwiseproductcaplet.cpp index 002e91740fd..8a4135d1fe3 100644 --- a/ql/models/marketmodels/products/pathwise/pathwiseproductcaplet.cpp +++ b/ql/models/marketmodels/products/pathwise/pathwiseproductcaplet.cpp @@ -17,6 +17,7 @@ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include #include diff --git a/ql/models/marketmodels/products/pathwise/pathwiseproductinversefloater.cpp b/ql/models/marketmodels/products/pathwise/pathwiseproductinversefloater.cpp index 32b7d85a1b2..19d7fd1e693 100644 --- a/ql/models/marketmodels/products/pathwise/pathwiseproductinversefloater.cpp +++ b/ql/models/marketmodels/products/pathwise/pathwiseproductinversefloater.cpp @@ -17,6 +17,7 @@ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include #include diff --git a/ql/models/marketmodels/products/pathwise/pathwiseproductswap.cpp b/ql/models/marketmodels/products/pathwise/pathwiseproductswap.cpp index b6f9d54f1bb..eb23ed420d0 100644 --- a/ql/models/marketmodels/products/pathwise/pathwiseproductswap.cpp +++ b/ql/models/marketmodels/products/pathwise/pathwiseproductswap.cpp @@ -17,6 +17,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include #include diff --git a/ql/models/marketmodels/products/pathwise/pathwiseproductswaption.cpp b/ql/models/marketmodels/products/pathwise/pathwiseproductswaption.cpp index 7a4b8e383a2..00a6fdebdcf 100644 --- a/ql/models/marketmodels/products/pathwise/pathwiseproductswaption.cpp +++ b/ql/models/marketmodels/products/pathwise/pathwiseproductswaption.cpp @@ -17,6 +17,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include #include diff --git a/ql/models/marketmodels/proxygreekengine.hpp b/ql/models/marketmodels/proxygreekengine.hpp index def327c367a..01b7730f3a1 100644 --- a/ql/models/marketmodels/proxygreekengine.hpp +++ b/ql/models/marketmodels/proxygreekengine.hpp @@ -26,6 +26,7 @@ #include #include +#include #include namespace QuantLib { diff --git a/ql/option.hpp b/ql/option.hpp index 99cb367eeb6..a0883001b5a 100644 --- a/ql/option.hpp +++ b/ql/option.hpp @@ -24,6 +24,7 @@ #ifndef quantlib_option_hpp #define quantlib_option_hpp +#include #include #include diff --git a/ql/patterns/lazyobject.hpp b/ql/patterns/lazyobject.hpp index f244cb6fd1e..a7eb82ccf58 100644 --- a/ql/patterns/lazyobject.hpp +++ b/ql/patterns/lazyobject.hpp @@ -24,6 +24,7 @@ #ifndef quantlib_lazy_object_h #define quantlib_lazy_object_h +#include #include #include diff --git a/ql/patterns/observable.cpp b/ql/patterns/observable.cpp index 39a7de18100..7872bf6f764 100644 --- a/ql/patterns/observable.cpp +++ b/ql/patterns/observable.cpp @@ -20,11 +20,17 @@ FOR A PARTICULAR PURPOSE. See the license for more details. #include +#include #ifndef QL_ENABLE_THREAD_SAFE_OBSERVER_PATTERN namespace QuantLib { + void ObservableSettings::disableUpdates(bool deferred) { + updatesEnabled_ = false; + updatesDeferred_ = deferred; + } + void ObservableSettings::enableUpdates() { updatesEnabled_ = true; updatesDeferred_ = false; @@ -52,6 +58,8 @@ namespace QuantLib { } } + bool ObservableSettings::updatesEnabled() const { return updatesEnabled_; } + bool ObservableSettings::updatesDeferred() const { return updatesDeferred_; } void Observable::notifyObservers() { if (!ObservableSettings::instance().updatesEnabled()) { @@ -83,6 +91,103 @@ namespace QuantLib { } } + Observable::Observable() = default; + + void ObservableSettings::registerDeferredObservers(const Observable::set_type& observers) { + if (updatesDeferred()) { + deferredObservers_.insert(observers.begin(), observers.end()); + } + } + + void ObservableSettings::unregisterDeferredObserver(Observer* o) { + deferredObservers_.erase(o); + } + + Observable::Observable(const Observable&) { + // the observer set is not copied; no observer asked to + // register with this object + } + + /*! \warning notification is sent before the copy constructor has + a chance of actually change the data + members. Therefore, observers whose update() method + tries to use their observables will not see the + updated values. It is suggested that the update() + method just raise a flag in order to trigger + a later recalculation. + */ + Observable& Observable::operator=(const Observable& o) { + // as above, the observer set is not copied. Moreover, + // observers of this object must be notified of the change + if (&o != this) + notifyObservers(); + return *this; + } + + std::pair + Observable::registerObserver(Observer* o) { + return observers_.insert(o); + } + + Size Observable::unregisterObserver(Observer* o) { + if (ObservableSettings::instance().updatesDeferred()) + ObservableSettings::instance().unregisterDeferredObserver(o); + + return observers_.erase(o); + } + + + Observer::Observer(const Observer& o) + : observables_(o.observables_) { + for (const auto& observable : observables_) + observable->registerObserver(this); + } + + Observer& Observer::operator=(const Observer& o) { + for (const auto& observable : observables_) + observable->unregisterObserver(this); + observables_ = o.observables_; + for (const auto& observable : observables_) + observable->registerObserver(this); + return *this; + } + + Observer::~Observer() { + for (const auto& observable : observables_) + observable->unregisterObserver(this); + } + + std::pair + Observer::registerWith(const ext::shared_ptr& h) { + if (h != nullptr) { + h->registerObserver(this); + return observables_.insert(h); + } + return std::make_pair(observables_.end(), false); + } + + void Observer::registerWithObservables(const ext::shared_ptr& o) { + if (o != nullptr) { + for (const auto& observable : o->observables_) + registerWith(observable); + } + } + + Size Observer::unregisterWith(const ext::shared_ptr& h) { + if (h != nullptr) + h->unregisterObserver(this); + return observables_.erase(h); + } + + void Observer::unregisterWithAll() { + for (const auto& observable : observables_) + observable->unregisterObserver(this); + observables_.clear(); + } + + void Observer::deepUpdate() { + update(); + } } #else @@ -138,6 +243,31 @@ namespace QuantLib { } + void Observer::Proxy::update() const { + std::lock_guard lock(mutex_); + if (active_) { + // c++17 is required if used with std::shared_ptr + const ext::weak_ptr o = observer_->weak_from_this(); + + // check for empty weak reference + // https://stackoverflow.com/questions/45507041/how-to-check-if-weak-ptr-is-empty-non-assigned + const ext::weak_ptr empty; + if (o.owner_before(empty) || empty.owner_before(o)) { + const ext::shared_ptr obs(o.lock()); + if (obs) + obs->update(); + } else { + observer_->update(); + } + } + } + + void Observer::Proxy::deactivate() { + std::lock_guard lock(mutex_); + active_ = false; + } + + void Observable::registerObserver(const ext::shared_ptr& observerProxy) { { std::lock_guard lock(mutex_); @@ -197,6 +327,149 @@ namespace QuantLib { // register with this object } + void ObservableSettings::registerDeferredObservers(const Observable::set_type& observers) { + deferredObservers_.insert(observers.begin(), observers.end()); + } + + void ObservableSettings::unregisterDeferredObserver( + const ext::shared_ptr& o) { + deferredObservers_.erase(o); + } + + void ObservableSettings::enableUpdates() { + std::lock_guard lock(mutex_); + + // if there are outstanding deferred updates, do the notification + updatesType_ = UpdatesEnabled; + + if (deferredObservers_.size()) { + bool successful = true; + std::string errMsg; + + for (auto i=deferredObservers_.begin(); + i!=deferredObservers_.end(); ++i) { + try { + const ext::shared_ptr proxy = i->lock(); + if (proxy) + proxy->update(); + } catch (std::exception& e) { + successful = false; + errMsg = e.what(); + } catch (...) { + successful = false; + } + } + + deferredObservers_.clear(); + + QL_ENSURE(successful, + "could not notify one or more observers: " << errMsg); + } + } + + /*! \warning notification is sent before the copy constructor has + a chance of actually change the data + members. Therefore, observers whose update() method + tries to use their observables will not see the + updated values. It is suggested that the update() + method just raise a flag in order to trigger + a later recalculation. + */ + Observable& Observable::operator=(const Observable& o) { + // as above, the observer set is not copied. Moreover, + // observers of this object must be notified of the change + if (&o != this) + notifyObservers(); + return *this; + } + + Observer::Observer(const Observer& o) { + proxy_.reset(new Proxy(this)); + + { + std::lock_guard lock(o.mutex_); + observables_ = o.observables_; + } + + for (const auto& observable : observables_) + observable->registerObserver(proxy_); + } + + Observer& Observer::operator=(const Observer& o) { + std::lock_guard lock(mutex_); + if (!proxy_) { + proxy_.reset(new Proxy(this)); + } + + for (const auto& observable : observables_) + observable->unregisterObserver(proxy_, true); + + { + std::lock_guard lock(o.mutex_); + observables_ = o.observables_; + } + for (const auto& observable : observables_) + observable->registerObserver(proxy_); + + return *this; + } + + Observer::~Observer() { + std::lock_guard lock(mutex_); + if (proxy_) + proxy_->deactivate(); + + for (const auto& observable : observables_) + observable->unregisterObserver(proxy_, false); + } + + std::pair + Observer::registerWith(const ext::shared_ptr& h) { + std::lock_guard lock(mutex_); + if (!proxy_) { + proxy_.reset(new Proxy(this)); + } + + if (h) { + h->registerObserver(proxy_); + return observables_.insert(h); + } + return std::make_pair(observables_.end(), false); + } + + void + Observer::registerWithObservables(const ext::shared_ptr& o) { + if (o) { + std::lock_guard lock(o->mutex_); + + for (const auto& observable : o->observables_) + registerWith(observable); + } + } + + Size Observer::unregisterWith(const ext::shared_ptr& h) { + std::lock_guard lock(mutex_); + + if (h && proxy_) { + h->unregisterObserver(proxy_, true); + } + + return observables_.erase(h); + } + + void Observer::unregisterWithAll() { + std::lock_guard lock(mutex_); + + for (const auto& observable : observables_) + observable->unregisterObserver(proxy_, true); + + observables_.clear(); + } + + void Observer::deepUpdate() { + update(); + } + } #endif diff --git a/ql/patterns/observable.hpp b/ql/patterns/observable.hpp index 35ead03ee70..46d680cb8ce 100644 --- a/ql/patterns/observable.hpp +++ b/ql/patterns/observable.hpp @@ -29,7 +29,6 @@ FOR A PARTICULAR PURPOSE. See the license for more details. #ifndef quantlib_observable_hpp #define quantlib_observable_hpp -#include #include #include #include @@ -88,14 +87,11 @@ namespace QuantLib { friend class Singleton; friend class Observable; public: - void disableUpdates(bool deferred=false) { - updatesEnabled_ = false; - updatesDeferred_ = deferred; - } + void disableUpdates(bool deferred=false); void enableUpdates(); - bool updatesEnabled() const { return updatesEnabled_; } - bool updatesDeferred() const { return updatesDeferred_; } + bool updatesEnabled() const; + bool updatesDeferred() const; private: ObservableSettings() = default; @@ -155,109 +151,6 @@ namespace QuantLib { set_type observables_; }; - - // inline definitions - - inline Observable::Observable() = default; - - inline void ObservableSettings::registerDeferredObservers(const Observable::set_type& observers) { - if (updatesDeferred()) { - deferredObservers_.insert(observers.begin(), observers.end()); - } - } - - inline void ObservableSettings::unregisterDeferredObserver(Observer* o) { - deferredObservers_.erase(o); - } - - inline Observable::Observable(const Observable&) { - // the observer set is not copied; no observer asked to - // register with this object - } - - /*! \warning notification is sent before the copy constructor has - a chance of actually change the data - members. Therefore, observers whose update() method - tries to use their observables will not see the - updated values. It is suggested that the update() - method just raise a flag in order to trigger - a later recalculation. - */ - inline Observable& Observable::operator=(const Observable& o) { - // as above, the observer set is not copied. Moreover, - // observers of this object must be notified of the change - if (&o != this) - notifyObservers(); - return *this; - } - - inline std::pair - Observable::registerObserver(Observer* o) { - return observers_.insert(o); - } - - inline Size Observable::unregisterObserver(Observer* o) { - if (ObservableSettings::instance().updatesDeferred()) - ObservableSettings::instance().unregisterDeferredObserver(o); - - return observers_.erase(o); - } - - - inline Observer::Observer(const Observer& o) - : observables_(o.observables_) { - for (const auto& observable : observables_) - observable->registerObserver(this); - } - - inline Observer& Observer::operator=(const Observer& o) { - for (const auto& observable : observables_) - observable->unregisterObserver(this); - observables_ = o.observables_; - for (const auto& observable : observables_) - observable->registerObserver(this); - return *this; - } - - inline Observer::~Observer() { - for (const auto& observable : observables_) - observable->unregisterObserver(this); - } - - inline std::pair - Observer::registerWith(const ext::shared_ptr& h) { - if (h != nullptr) { - h->registerObserver(this); - return observables_.insert(h); - } - return std::make_pair(observables_.end(), false); - } - - inline void - Observer::registerWithObservables(const ext::shared_ptr &o) { - if (o != nullptr) { - for (const auto& observable : o->observables_) - registerWith(observable); - } - } - - inline - Size Observer::unregisterWith(const ext::shared_ptr& h) { - if (h != nullptr) - h->unregisterObserver(this); - return observables_.erase(h); - } - - inline void Observer::unregisterWithAll() { - for (const auto& observable : observables_) - observable->unregisterObserver(this); - observables_.clear(); - } - - inline void Observer::deepUpdate() { - update(); - } - } #else @@ -322,31 +215,9 @@ namespace QuantLib { observer_(observer) { } - void update() const { - std::lock_guard lock(mutex_); - if (active_) { - // c++17 is required if used with std::shared_ptr - const ext::weak_ptr o - = observer_->weak_from_this(); - - //check for empty weak reference - //https://stackoverflow.com/questions/45507041/how-to-check-if-weak-ptr-is-empty-non-assigned - const ext::weak_ptr empty; - if (o.owner_before(empty) || empty.owner_before(o)) { - const ext::shared_ptr obs(o.lock()); - if (obs) - obs->update(); - } - else { - observer_->update(); - } - } - } + void update() const; - void deactivate() { - std::lock_guard lock(mutex_); - active_ = false; - } + void deactivate(); private: bool active_; @@ -424,153 +295,6 @@ namespace QuantLib { std::atomic updatesType_; }; - - // inline definitions - - inline void ObservableSettings::registerDeferredObservers(const Observable::set_type& observers) { - deferredObservers_.insert(observers.begin(), observers.end()); - } - - inline void ObservableSettings::unregisterDeferredObserver( - const ext::shared_ptr& o) { - deferredObservers_.erase(o); - } - - inline void ObservableSettings::enableUpdates() { - std::lock_guard lock(mutex_); - - // if there are outstanding deferred updates, do the notification - updatesType_ = UpdatesEnabled; - - if (deferredObservers_.size()) { - bool successful = true; - std::string errMsg; - - for (auto i=deferredObservers_.begin(); - i!=deferredObservers_.end(); ++i) { - try { - const ext::shared_ptr proxy = i->lock(); - if (proxy) - proxy->update(); - } catch (std::exception& e) { - successful = false; - errMsg = e.what(); - } catch (...) { - successful = false; - } - } - - deferredObservers_.clear(); - - QL_ENSURE(successful, - "could not notify one or more observers: " << errMsg); - } - } - - - /*! \warning notification is sent before the copy constructor has - a chance of actually change the data - members. Therefore, observers whose update() method - tries to use their observables will not see the - updated values. It is suggested that the update() - method just raise a flag in order to trigger - a later recalculation. - */ - inline Observable& Observable::operator=(const Observable& o) { - // as above, the observer set is not copied. Moreover, - // observers of this object must be notified of the change - if (&o != this) - notifyObservers(); - return *this; - } - - inline Observer::Observer(const Observer& o) { - proxy_.reset(new Proxy(this)); - - { - std::lock_guard lock(o.mutex_); - observables_ = o.observables_; - } - - for (const auto& observable : observables_) - observable->registerObserver(proxy_); - } - - inline Observer& Observer::operator=(const Observer& o) { - std::lock_guard lock(mutex_); - if (!proxy_) { - proxy_.reset(new Proxy(this)); - } - - for (const auto& observable : observables_) - observable->unregisterObserver(proxy_, true); - - { - std::lock_guard lock(o.mutex_); - observables_ = o.observables_; - } - for (const auto& observable : observables_) - observable->registerObserver(proxy_); - - return *this; - } - - inline Observer::~Observer() { - std::lock_guard lock(mutex_); - if (proxy_) - proxy_->deactivate(); - - for (const auto& observable : observables_) - observable->unregisterObserver(proxy_, false); - } - - inline std::pair - Observer::registerWith(const ext::shared_ptr& h) { - std::lock_guard lock(mutex_); - if (!proxy_) { - proxy_.reset(new Proxy(this)); - } - - if (h) { - h->registerObserver(proxy_); - return observables_.insert(h); - } - return std::make_pair(observables_.end(), false); - } - - inline void - Observer::registerWithObservables(const ext::shared_ptr& o) { - if (o) { - std::lock_guard lock(o->mutex_); - - for (const auto& observable : o->observables_) - registerWith(observable); - } - } - - inline - Size Observer::unregisterWith(const ext::shared_ptr& h) { - std::lock_guard lock(mutex_); - - if (h && proxy_) { - h->unregisterObserver(proxy_, true); - } - - return observables_.erase(h); - } - - inline void Observer::unregisterWithAll() { - std::lock_guard lock(mutex_); - - for (const auto& observable : observables_) - observable->unregisterObserver(proxy_, true); - - observables_.clear(); - } - - inline void Observer::deepUpdate() { - update(); - } } #endif #endif diff --git a/ql/processes/gsrprocesscore.cpp b/ql/processes/gsrprocesscore.cpp index 271c2b6cc60..9bddfc458c2 100644 --- a/ql/processes/gsrprocesscore.cpp +++ b/ql/processes/gsrprocesscore.cpp @@ -17,6 +17,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include diff --git a/ql/termstructures/volatility/abcd.hpp b/ql/termstructures/volatility/abcd.hpp index 9910676f997..8e7f32b66a0 100644 --- a/ql/termstructures/volatility/abcd.hpp +++ b/ql/termstructures/volatility/abcd.hpp @@ -25,6 +25,7 @@ #include #include #include +#include namespace QuantLib { diff --git a/ql/time/asx.cpp b/ql/time/asx.cpp index df046d8eb2c..c5ddffdb6cf 100644 --- a/ql/time/asx.cpp +++ b/ql/time/asx.cpp @@ -21,6 +21,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include #include diff --git a/ql/time/calendar.cpp b/ql/time/calendar.cpp index 9773284fe1a..3ff1000b1a5 100644 --- a/ql/time/calendar.cpp +++ b/ql/time/calendar.cpp @@ -295,4 +295,76 @@ namespace QuantLib { } return result; } + + bool Calendar::empty() const { + return !impl_; + } + + std::string Calendar::name() const { + QL_REQUIRE(impl_, "no calendar implementation provided"); + return impl_->name(); + } + + const std::set& Calendar::addedHolidays() const { + QL_REQUIRE(impl_, "no calendar implementation provided"); + + return impl_->addedHolidays; + } + + const std::set& Calendar::removedHolidays() const { + QL_REQUIRE(impl_, "no calendar implementation provided"); + + return impl_->removedHolidays; + } + + bool Calendar::isBusinessDay(const Date& d) const { + QL_REQUIRE(impl_, "no calendar implementation provided"); + +#ifdef QL_HIGH_RESOLUTION_DATE + const Date _d(d.dayOfMonth(), d.month(), d.year()); +#else + const Date& _d = d; +#endif + + if (!impl_->addedHolidays.empty() && + impl_->addedHolidays.find(_d) != impl_->addedHolidays.end()) + return false; + + if (!impl_->removedHolidays.empty() && + impl_->removedHolidays.find(_d) != impl_->removedHolidays.end()) + return true; + + return impl_->isBusinessDay(_d); + } + + bool Calendar::isEndOfMonth(const Date& d) const { + return (d.month() != adjust(d+1).month()); + } + + Date Calendar::endOfMonth(const Date& d) const { + return adjust(Date::endOfMonth(d), Preceding); + } + + bool Calendar::isHoliday(const Date& d) const { + return !isBusinessDay(d); + } + + bool Calendar::isWeekend(Weekday w) const { + QL_REQUIRE(impl_, "no calendar implementation provided"); + return impl_->isWeekend(w); + } + + bool operator==(const Calendar& c1, const Calendar& c2) { + return (c1.empty() && c2.empty()) + || (!c1.empty() && !c2.empty() && c1.name() == c2.name()); + } + + bool operator!=(const Calendar& c1, const Calendar& c2) { + return !(c1 == c2); + } + + std::ostream& operator<<(std::ostream& out, const Calendar &c) { + return out << c.name(); + } + } diff --git a/ql/time/calendar.hpp b/ql/time/calendar.hpp index 35b603260f2..a6a91b614f5 100644 --- a/ql/time/calendar.hpp +++ b/ql/time/calendar.hpp @@ -28,7 +28,6 @@ #ifndef quantlib_calendar_hpp #define quantlib_calendar_hpp -#include #include #include #include @@ -196,80 +195,6 @@ namespace QuantLib { /*! \relates Calendar */ std::ostream& operator<<(std::ostream&, const Calendar&); - - // inline definitions - - inline bool Calendar::empty() const { - return !impl_; - } - - inline std::string Calendar::name() const { - QL_REQUIRE(impl_, "no calendar implementation provided"); - return impl_->name(); - } - - inline const std::set& Calendar::addedHolidays() const { - QL_REQUIRE(impl_, "no calendar implementation provided"); - - return impl_->addedHolidays; - } - - inline const std::set& Calendar::removedHolidays() const { - QL_REQUIRE(impl_, "no calendar implementation provided"); - - return impl_->removedHolidays; - } - - inline bool Calendar::isBusinessDay(const Date& d) const { - QL_REQUIRE(impl_, "no calendar implementation provided"); - -#ifdef QL_HIGH_RESOLUTION_DATE - const Date _d(d.dayOfMonth(), d.month(), d.year()); -#else - const Date& _d = d; -#endif - - if (!impl_->addedHolidays.empty() && - impl_->addedHolidays.find(_d) != impl_->addedHolidays.end()) - return false; - - if (!impl_->removedHolidays.empty() && - impl_->removedHolidays.find(_d) != impl_->removedHolidays.end()) - return true; - - return impl_->isBusinessDay(_d); - } - - inline bool Calendar::isEndOfMonth(const Date& d) const { - return (d.month() != adjust(d+1).month()); - } - - inline Date Calendar::endOfMonth(const Date& d) const { - return adjust(Date::endOfMonth(d), Preceding); - } - - inline bool Calendar::isHoliday(const Date& d) const { - return !isBusinessDay(d); - } - - inline bool Calendar::isWeekend(Weekday w) const { - QL_REQUIRE(impl_, "no calendar implementation provided"); - return impl_->isWeekend(w); - } - - inline bool operator==(const Calendar& c1, const Calendar& c2) { - return (c1.empty() && c2.empty()) - || (!c1.empty() && !c2.empty() && c1.name() == c2.name()); - } - - inline bool operator!=(const Calendar& c1, const Calendar& c2) { - return !(c1 == c2); - } - - inline std::ostream& operator<<(std::ostream& out, const Calendar &c) { - return out << c.name(); - } - } #endif diff --git a/ql/time/calendars/australia.cpp b/ql/time/calendars/australia.cpp index b5f970c4165..0213b26dc46 100644 --- a/ql/time/calendars/australia.cpp +++ b/ql/time/calendars/australia.cpp @@ -17,6 +17,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include namespace QuantLib { @@ -107,4 +108,4 @@ namespace QuantLib { return true; } -} \ No newline at end of file +} diff --git a/ql/time/date.cpp b/ql/time/date.cpp index 165922d53a3..31701dbc855 100644 --- a/ql/time/date.cpp +++ b/ql/time/date.cpp @@ -1003,4 +1003,83 @@ namespace QuantLib { } #endif } + +#ifndef QL_HIGH_RESOLUTION_DATE + + Weekday Date::weekday() const { + Integer w = serialNumber_ % 7; + return Weekday(w == 0 ? 7 : w); + } + + Day Date::dayOfMonth() const { + return dayOfYear() - monthOffset(month(),isLeap(year())); + } + + Day Date::dayOfYear() const { + return serialNumber_ - yearOffset(year()); + } + + Date::serial_type Date::serialNumber() const { + return serialNumber_; + } + + Date Date::operator+(Date::serial_type days) const { + return Date(serialNumber_+days); + } + + Date Date::operator-(Date::serial_type days) const { + return Date(serialNumber_-days); + } + + Date Date::operator+(const Period& p) const { + return advance(*this,p.length(),p.units()); + } + + Date Date::operator-(const Period& p) const { + return advance(*this,-p.length(),p.units()); + } + + Date Date::endOfMonth(const Date& d) { + Month m = d.month(); + Year y = d.year(); + return {monthLength(m, isLeap(y)), m, y}; + } + + bool Date::isEndOfMonth(const Date& d) { + return (d.dayOfMonth() == monthLength(d.month(), isLeap(d.year()))); + } + + Date::serial_type operator-(const Date& d1, const Date& d2) { + return d1.serialNumber()-d2.serialNumber(); + } + + Time daysBetween(const Date& d1, const Date& d2) { + return Time(d2-d1); + } + + bool operator==(const Date& d1, const Date& d2) { + return (d1.serialNumber() == d2.serialNumber()); + } + + bool operator!=(const Date& d1, const Date& d2) { + return (d1.serialNumber() != d2.serialNumber()); + } + + bool operator<(const Date& d1, const Date& d2) { + return (d1.serialNumber() < d2.serialNumber()); + } + + bool operator<=(const Date& d1, const Date& d2) { + return (d1.serialNumber() <= d2.serialNumber()); + } + + bool operator>(const Date& d1, const Date& d2) { + return (d1.serialNumber() > d2.serialNumber()); + } + + bool operator>=(const Date& d1, const Date& d2) { + return (d1.serialNumber() >= d2.serialNumber()); + } +#endif + } diff --git a/ql/time/date.hpp b/ql/time/date.hpp index 3e0379cfeac..a0b0db2e556 100644 --- a/ql/time/date.hpp +++ b/ql/time/date.hpp @@ -388,85 +388,6 @@ namespace QuantLib { }; #endif - -#ifndef QL_HIGH_RESOLUTION_DATE - // inline definitions - - inline Weekday Date::weekday() const { - Integer w = serialNumber_ % 7; - return Weekday(w == 0 ? 7 : w); - } - - inline Day Date::dayOfMonth() const { - return dayOfYear() - monthOffset(month(),isLeap(year())); - } - - inline Day Date::dayOfYear() const { - return serialNumber_ - yearOffset(year()); - } - - inline Date::serial_type Date::serialNumber() const { - return serialNumber_; - } - - inline Date Date::operator+(Date::serial_type days) const { - return Date(serialNumber_+days); - } - - inline Date Date::operator-(Date::serial_type days) const { - return Date(serialNumber_-days); - } - - inline Date Date::operator+(const Period& p) const { - return advance(*this,p.length(),p.units()); - } - - inline Date Date::operator-(const Period& p) const { - return advance(*this,-p.length(),p.units()); - } - - inline Date Date::endOfMonth(const Date& d) { - Month m = d.month(); - Year y = d.year(); - return {monthLength(m, isLeap(y)), m, y}; - } - - inline bool Date::isEndOfMonth(const Date& d) { - return (d.dayOfMonth() == monthLength(d.month(), isLeap(d.year()))); - } - - inline Date::serial_type operator-(const Date& d1, const Date& d2) { - return d1.serialNumber()-d2.serialNumber(); - } - - inline Time daysBetween(const Date& d1, const Date& d2) { - return Time(d2-d1); - } - - inline bool operator==(const Date& d1, const Date& d2) { - return (d1.serialNumber() == d2.serialNumber()); - } - - inline bool operator!=(const Date& d1, const Date& d2) { - return (d1.serialNumber() != d2.serialNumber()); - } - - inline bool operator<(const Date& d1, const Date& d2) { - return (d1.serialNumber() < d2.serialNumber()); - } - - inline bool operator<=(const Date& d1, const Date& d2) { - return (d1.serialNumber() <= d2.serialNumber()); - } - - inline bool operator>(const Date& d1, const Date& d2) { - return (d1.serialNumber() > d2.serialNumber()); - } - - inline bool operator>=(const Date& d1, const Date& d2) { - return (d1.serialNumber() >= d2.serialNumber()); - } -#endif } namespace std { diff --git a/ql/time/daycounter.hpp b/ql/time/daycounter.hpp index 8c24cd8a32e..3164bb0b6d1 100644 --- a/ql/time/daycounter.hpp +++ b/ql/time/daycounter.hpp @@ -26,6 +26,7 @@ #define quantlib_day_counter_hpp #include +#include #include #include diff --git a/ql/time/ecb.cpp b/ql/time/ecb.cpp index c15496ccb5b..ddaad4e169b 100644 --- a/ql/time/ecb.cpp +++ b/ql/time/ecb.cpp @@ -18,6 +18,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include #include diff --git a/ql/time/imm.cpp b/ql/time/imm.cpp index 58a3afab573..76b1da0d4b8 100644 --- a/ql/time/imm.cpp +++ b/ql/time/imm.cpp @@ -20,6 +20,7 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include #include #include #include