diff --git a/.gitignore b/.gitignore index 51b1a3a..34b4884 100644 --- a/.gitignore +++ b/.gitignore @@ -310,10 +310,15 @@ $RECYCLE.BIN/ *.*.*.bak *.*.*.*.bak *.*.*.*.*.bak -*.txt + *.db *.opendb *.*.opendb *.*.*.opendb *.*.*.*.opendb +build +out +.cache +.vscode +build_v* \ No newline at end of file diff --git a/AADExpr.h b/AADExpr.h deleted file mode 100644 index af44a25..0000000 --- a/AADExpr.h +++ /dev/null @@ -1,1135 +0,0 @@ - -/* -Written by Antoine Savine in 2018 - -This code is the strict IP of Antoine Savine - -License to use and alter this code for personal and commercial applications -is freely granted to any person or company who purchased a copy of the book - -Modern Computational Finance: AAD and Parallel Simulations -Antoine Savine -Wiley, 2018 - -As long as this comment is preserved at the top of the file -*/ - -#pragma once - -// Implementation of AAD with expression templates -// (AADET, chapter 15) - -// Defines expressions and the Number type - -#include -#include "AADTape.h" - -// Base CRTP expression class -// Note: overloaded operators catch all expressions and nothing else -template -struct Expression -{ - // CRTP "virtualization" - double value() const - { - return static_cast(this)->value(); - } - - // Just another interface - explicit operator double() const - { - return value(); - } -}; - -// Note that Number is a leaf expression -// Defined in the bottom of the file - -// Binary expression -// LHS : the expression on the left -// RHS : the expression on the right -// OP : the binary operator -template -class BinaryExpression - // CRTP - : public Expression> -{ - const double myValue; - - const LHS lhs; - const RHS rhs; - -public: - - // Constructor out of 2 expressions - // Note: eager evaluation on construction - explicit BinaryExpression( - const Expression& l, - const Expression& r) - : myValue(OP::eval(l.value(), r.value())), - lhs(static_cast(l)), - rhs(static_cast(r)) - {} - - // Value accessors - double value() const { return myValue; } - - // Expression template magic - // Expressions know - // AT COMPILE TIME - // the number of active inputs in their sub-expressions - enum { numNumbers = LHS::numNumbers + RHS::numNumbers }; - - // Push adjoint down the expression - // N : total number of active inputs in the expression - // n : number of active inputs already processed - template - void pushAdjoint( - // Node for the complete expression being processed - Node& exprNode, - // Adjoint cumulated for this binary node, or 1 if top - const double adjoint) - const - { - // Push on the left - if (LHS::numNumbers > 0) - { - lhs.pushAdjoint( - exprNode, - adjoint * OP::leftDerivative(lhs.value(), rhs.value(), value())); - } - - // Push on the right - if (RHS::numNumbers > 0) - { - // Note left push processed LHS::numNumbers numbers - // So the next number to be processed is n + LHS::numNumbers - rhs.pushAdjoint( - exprNode, - adjoint * OP::rightDerivative(lhs.value(), rhs.value(), value())); - } - } -}; - -// "Concrete" binaries, we only need to define operations and derivatives -struct OPMult -{ - static const double eval(const double l, const double r) - { - return l * r; - } - - static const double leftDerivative - (const double l, const double r, const double v) - { - return r; - } - - static const double rightDerivative - (const double l, const double r, const double v) - { - return l; - } -}; - -struct OPAdd -{ - static const double eval(const double l, const double r) - { - return l + r; - } - - static const double leftDerivative - (const double l, const double r, const double v) - { - return 1.0; - } - - static const double rightDerivative - (const double l, const double r, const double v) - { - return 1.0; - } -}; - -struct OPSub -{ - static const double eval(const double l, const double r) - { - return l - r; - } - - static const double leftDerivative - (const double l, const double r, const double v) - { - return 1.0; - } - - static const double rightDerivative - (const double l, const double r, const double v) - { - return -1.0; - } -}; - -struct OPDiv -{ - static const double eval(const double l, const double r) - { - return l / r; - } - - static const double leftDerivative - (const double l, const double r, const double v) - { - return 1.0 / r; - } - - static const double rightDerivative - (const double l, const double r, const double v) - { - return -l / r / r; - } -}; - -struct OPPow -{ - static const double eval(const double l, const double r) - { - return pow(l, r); - } - - static const double leftDerivative - (const double l, const double r, const double v) - { - return r*v / l; - } - - static const double rightDerivative - (const double l, const double r, const double v) - { - return log(l)*v; - } -}; - -struct OPMax -{ - static const double eval(const double l, const double r) - { - return max(l, r); - } - - static const double leftDerivative - (const double l, const double r, const double v) - { - return l > r ? 1.0 : 0.0; - } - - static const double rightDerivative - (const double l, const double r, const double v) - { - return r > l? 1.0 : 0.0; - } -}; - -struct OPMin -{ - static const double eval(const double l, const double r) - { - return min(l, r); - } - - static const double leftDerivative - (const double l, const double r, const double v) - { - return l < r ? 1.0 : 0.0; - } - - static const double rightDerivative - (const double l, const double r, const double v) - { - return r < l ? 1.0 : 0.0; - } -}; - -// Operator overloading for binary expressions -// build the corresponding expressions - -template - BinaryExpression operator*( - const Expression& lhs, const Expression& rhs) -{ - return BinaryExpression(lhs, rhs); -} - -template - BinaryExpression operator+( - const Expression& lhs, const Expression& rhs) -{ - return BinaryExpression(lhs, rhs); -} - -template - BinaryExpression operator-( - const Expression& lhs, const Expression& rhs) -{ - return BinaryExpression(lhs, rhs); -} - -template - BinaryExpression operator/( - const Expression& lhs, const Expression& rhs) -{ - return BinaryExpression(lhs, rhs); -} - -template - BinaryExpression pow( - const Expression& lhs, const Expression& rhs) -{ - return BinaryExpression(lhs, rhs); -} - -template -BinaryExpression max( - const Expression& lhs, const Expression& rhs) -{ - return BinaryExpression(lhs, rhs); -} - -template -BinaryExpression min( - const Expression& lhs, const Expression& rhs) -{ - return BinaryExpression(lhs, rhs); -} - -// Unary expressions : Same logic with one argument - -// The CRTP class -template -class UnaryExpression - // CRTP - : public Expression> -{ - const double myValue; - - const ARG arg; - - // For binary operators with a double on one side - // we store the double - const double dArg = 0.0; - -public: - - // Constructor - // Note : eager evaluation on construction - explicit UnaryExpression( - const Expression& a) - : myValue(OP::eval(a.value(), 0.0)), arg(static_cast(a)) {} - - // Special constructor for binary expressions with a double on one side - explicit UnaryExpression( - const Expression& a, - const double b) - : myValue(OP::eval(a.value(), b)), arg(static_cast(a)), dArg(b) {} - - // Value accessors - double value() const { return myValue; } - - // Expression template magic - enum { numNumbers = ARG::numNumbers }; - - // Push adjoint down the expression - template - void pushAdjoint( - // Node for the complete expression being processed - Node& exprNode, - // Adjoint cumulated on the node, 1 if top - const double adjoint) - const - { - // Push into argument - if (ARG::numNumbers > 0) - { - arg.pushAdjoint( - exprNode, - adjoint * OP::derivative(arg.value(), value(), dArg)); - } - } -}; - -// The unary operators - -struct OPExp -{ - static const double eval(const double r, const double d) - { - return exp(r); - } - - static const double derivative - (const double r, const double v, const double d) - { - return v; - } -}; - -struct OPLog -{ - static const double eval(const double r, const double d) - { - return log(r); - } - - static const double derivative - (const double r, const double v, const double d) - { - return 1.0 / r; - } -}; - -struct OPSqrt -{ - static const double eval(const double r, const double d) - { - return sqrt(r); - } - - static const double derivative - (const double r, const double v, const double d) - { - return 0.5 / v; - } -}; - -struct OPFabs -{ - static const double eval(const double r, const double d) - { - return fabs(r); - } - - static const double derivative - (const double r, const double v, const double d) - { - return r > 0.0 ? 1.0 : -1.0; - } -}; - -struct OPNormalDens -{ - static const double eval(const double r, const double d) - { - return normalDens(r); - } - - static const double derivative - (const double r, const double v, const double d) - { - return - r * v; - } -}; - -struct OPNormalCdf -{ - static const double eval(const double r, const double d) - { - return normalCdf(r); - } - - static const double derivative - (const double r, const double v, const double d) - { - return normalDens(r); - } -}; - -// Binary operators with a double on one side - -// * double or double * -struct OPMultD -{ - static const double eval(const double r, const double d) - { - return r * d; - } - - static const double derivative - (const double r, const double v, const double d) - { - return d; - } -}; - -// + double or double + -struct OPAddD -{ - static const double eval(const double r, const double d) - { - return r + d; - } - - static const double derivative - (const double r, const double v, const double d) - { - return 1.0; - } -}; - -// double - -struct OPSubDL -{ - static const double eval(const double r, const double d) - { - return d - r; - } - - static const double derivative - (const double r, const double v, const double d) - { - return -1.0; - } -}; - -// - double -struct OPSubDR -{ - static const double eval(const double r, const double d) - { - return r - d; - } - - static const double derivative - (const double r, const double v, const double d) - { - return 1.0; - } -}; - -// double / -struct OPDivDL -{ - static const double eval(const double r, const double d) - { - return d / r; - } - - static const double derivative - (const double r, const double v, const double d) - { - return -d / r / r; - } -}; - -// / double -struct OPDivDR -{ - static const double eval(const double r, const double d) - { - return r / d; - } - - static const double derivative - (const double r, const double v, const double d) - { - return 1.0 / d; - } -}; - -// pow (d,) -struct OPPowDL -{ - static const double eval(const double r, const double d) - { - return pow(d, r); - } - - static const double derivative - (const double r, const double v, const double d) - { - return log(d) * v; - } -}; - -// pow (,d) -struct OPPowDR -{ - static const double eval(const double r, const double d) - { - return pow(r, d); - } - - static const double derivative - (const double r, const double v, const double d) - { - return d * v / r; - } -}; - -// max (d,) -struct OPMaxD -{ - static const double eval(const double r, const double d) - { - return max(r, d); - } - - static const double derivative - (const double r, const double v, const double d) - { - return r > d ? 1.0 : 0.0; - } -}; - -// min (d,) -struct OPMinD -{ - static const double eval(const double r, const double d) - { - return min(r, d); - } - - static const double derivative - (const double r, const double v, const double d) - { - return r < d ? 1.0 : 0.0; - } -}; - -// And overloading - -template - UnaryExpression exp(const Expression& arg) -{ - return UnaryExpression(arg); -} - -template - UnaryExpression log(const Expression& arg) -{ - return UnaryExpression(arg); -} - -template - UnaryExpression sqrt(const Expression& arg) -{ - return UnaryExpression(arg); -} - -template -UnaryExpression fabs(const Expression& arg) -{ - return UnaryExpression(arg); -} - -template -UnaryExpression normalDens(const Expression& arg) -{ - return UnaryExpression(arg); -} - -template -UnaryExpression normalCdf(const Expression& arg) -{ - return UnaryExpression(arg); -} - -// Overloading continued, -// binary operators with a double on one side - -template - UnaryExpression operator*( - const double d, const Expression& rhs) -{ - return UnaryExpression(rhs, d); -} - -template - UnaryExpression operator*( - const Expression& lhs, const double d) -{ - return UnaryExpression(lhs, d); -} - -template - UnaryExpression operator+( - const double d, const Expression& rhs) -{ - return UnaryExpression(rhs, d); -} - -template - UnaryExpression operator+( - const Expression& lhs, const double d) -{ - return UnaryExpression(lhs, d); -} - -template - UnaryExpression operator-( - const double d, const Expression& rhs) -{ - return UnaryExpression(rhs, d); -} - -template - UnaryExpression operator-( - const Expression& lhs, const double d) -{ - return UnaryExpression(lhs, d); -} - -template - UnaryExpression operator/( - const double d, const Expression& rhs) -{ - return UnaryExpression(rhs, d); -} - -template - UnaryExpression operator/( - const Expression& lhs, const double d) -{ - return UnaryExpression(lhs, d); -} - -template - UnaryExpression pow( - const double d, const Expression& rhs) -{ - return UnaryExpression(rhs, d); -} - -template - UnaryExpression pow( - const Expression& lhs, const double d) -{ - return UnaryExpression(lhs, d); -} - - -template -UnaryExpression max( - const double d, const Expression& rhs) -{ - return UnaryExpression(rhs, d); -} - -template -UnaryExpression max( - const Expression& lhs, const double d) -{ - return UnaryExpression(lhs, d); -} - -template -UnaryExpression min( - const double d, const Expression& rhs) -{ - return UnaryExpression(rhs, d); -} - -template -UnaryExpression min( - const Expression& lhs, const double d) -{ - return UnaryExpression(lhs, d); -} - -// Comparison, same as traditional - -template - bool operator==(const Expression& lhs, const Expression& rhs) -{ - return lhs.value() == rhs.value(); -} -template - bool operator==(const Expression& lhs, const double& rhs) -{ - return lhs.value() == rhs; -} -template - bool operator==(const double& lhs, const Expression& rhs) -{ - return lhs == rhs.value(); -} - -template - bool operator!=(const Expression& lhs, const Expression& rhs) -{ - return lhs.value() != rhs.value(); -} -template - bool operator!=(const Expression& lhs, const double& rhs) -{ - return lhs.value() != rhs; -} -template - bool operator!=(const double& lhs, const Expression& rhs) -{ - return lhs != rhs.value(); -} - -template - bool operator<(const Expression& lhs, const Expression& rhs) -{ - return lhs.value() < rhs.value(); -} -template - bool operator<(const Expression& lhs, const double& rhs) -{ - return lhs.value() < rhs; -} -template - bool operator<(const double& lhs, const Expression& rhs) -{ - return lhs < rhs.value(); -} - -template - bool operator>(const Expression& lhs, const Expression& rhs) -{ - return lhs.value() > rhs.value(); -} -template - bool operator>(const Expression& lhs, const double& rhs) -{ - return lhs.value() > rhs; -} -template - bool operator>(const double& lhs, const Expression& rhs) -{ - return lhs > rhs.value(); -} - -template - bool operator<=(const Expression& lhs, const Expression& rhs) -{ - return lhs.value() <= rhs.value(); -} -template - bool operator<=(const Expression& lhs, const double& rhs) -{ - return lhs.value() <= rhs; -} -template - bool operator<=(const double& lhs, const Expression& rhs) -{ - return lhs <= rhs.value(); -} - -template - bool operator>=(const Expression& lhs, const Expression& rhs) -{ - return lhs.value() >= rhs.value(); -} -template - bool operator>=(const Expression& lhs, const double& rhs) -{ - return lhs.value() >= rhs; -} -template - bool operator>=(const double& lhs, const Expression& rhs) -{ - return lhs >= rhs.value(); -} - -// Finally, unary +/- operators - -template -UnaryExpression operator- -(const Expression& rhs) -{ - return 0.0 - rhs; -} - -template -Expression operator+ -(const Expression& rhs) -{ - return rhs; -} - -// The Number type, also an expression - -class Number : public Expression -{ - // The value and node for this number, same as traditional - double myValue; - Node* myNode; - - // Node creation on tape - - template - Node* createMultiNode() - { - return tape->recordNode(); - } - - // Flattening: - // This is where, on assignment or construction from an expression, - // that derivatives are pushed through the expression's DAG - template - void fromExpr( - // RHS expression, will be flattened into this Number - const Expression& e) - { - // Build expression node on tape - auto* node = createMultiNode(); - - // Push adjoints through expression with adjoint = 1 on top - static_cast(e).pushAdjoint(*node, 1.0); - - // Set my node - myNode = node; - } - -public: - - // Expression template magic - enum { numNumbers = 1 }; - - // Push adjoint - // Numbers are expression leaves, - // pushAdjoint() receives their adjoint in the expression - // Numbers don't "push" anything, they register their derivatives on tape - template - void pushAdjoint( - // Node for the complete expression - Node& exprNode, - // Adjoint accumulated for this number, in the expression - const double adjoint) - const - { - // adjoint = d (expression) / d (thisNumber) : register on tape - // note n: index of this number on the node on tape - - // Register adjoint - exprNode.pAdjPtrs[n] = Tape::multi? myNode->pAdjoints : &myNode->mAdjoint; - - // Register derivative - exprNode.pDerivatives[n] = adjoint; - } - - // Static access to tape, same as traditional - static thread_local Tape* tape; - - // Constructors - - Number() {} - - explicit Number(const double val) : myValue(val) - { - // Create leaf - myNode = createMultiNode<0>(); - } - - Number& operator=(const double val) - { - myValue = val; - // Create leaf - myNode = createMultiNode<0>(); - return *this; - } - - // No need for copy and assignment - // Default ones do the right thing: - // copy value and pointer to node on tape - - // Construct or assign from expression - - template - Number(const Expression& e) : myValue(e.value()) - { - // Flatten RHS expression - fromExpr(static_cast(e)); - } - - template - Number& operator= - (const Expression& e) - { - myValue = e.value(); - // Flatten RHS expression - fromExpr(static_cast(e)); - return *this; - } - - // Explicit coversion to double - explicit operator double& () { return myValue; } - explicit operator double () const { return myValue; } - - // All the normal accessors and propagators, same as traditional - - // Put on tape - void putOnTape() - { - myNode = createMultiNode<0>(); - } - - // Accessors: value and adjoint - - double& value() - { - return myValue; - } - double value() const - { - return myValue; - } - - // Single dimensional - double& adjoint() - { - return myNode->adjoint(); - } - double adjoint() const - { - return myNode->adjoint(); - } - - // Multi dimensional - double& adjoint(const size_t n) - { - return myNode->adjoint(n); - } - double adjoint(const size_t n) const - { - return myNode->adjoint(n); - } - - // Reset all adjoints on the tape - // note we don't use this method - void resetAdjoints() - { - tape->resetAdjoints(); - } - - // Propagation - - // Propagate adjoints - // from and to both INCLUSIVE - static void propagateAdjoints( - Tape::iterator propagateFrom, - Tape::iterator propagateTo) - { - auto it = propagateFrom; - while (it != propagateTo) - { - it->propagateOne(); - --it; - } - it->propagateOne(); - } - - // Convenient overloads - - // Set the adjoint on this node to 1, - // Then propagate from the node - void propagateAdjoints( - // We start on this number's node - Tape::iterator propagateTo) - { - // Set this adjoint to 1 - adjoint() = 1.0; - // Find node on tape - auto it = tape->find(myNode); - // Reverse and propagate until we hit the stop - while (it != propagateTo) - { - it->propagateOne(); - --it; - } - it->propagateOne(); - } - - // These 2 set the adjoint to 1 on this node - void propagateToStart() - { - propagateAdjoints(tape->begin()); - } - void propagateToMark() - { - propagateAdjoints(tape->markIt()); - } - - // This one only propagates - // Note: propagation starts at mark - 1 - static void propagateMarkToStart() - { - propagateAdjoints(prev(tape->markIt()), tape->begin()); - } - - // Multi-adjoint propagation - - // Multi dimensional case: - // Propagate adjoints from and to both INCLUSIVE - static void propagateAdjointsMulti( - Tape::iterator propagateFrom, - Tape::iterator propagateTo) - { - auto it = propagateFrom; - while (it != propagateTo) - { - it->propagateAll(); - --it; - } - it->propagateAll(); - } - - // Unary operators - - template - Number& operator+=(const Expression& e) - { - *this = *this + e; - return *this; - } - - template - Number& operator*=(const Expression& e) - { - *this = *this * e; - return *this; - } - - template - Number& operator-=(const Expression& e) - { - *this = *this - e; - return *this; - } - - template - Number& operator/=(const Expression& e) - { - *this = *this / e; - return *this; - } - - Number& operator+=(const double& e) - { - *this = *this + e; - return *this; - } - - Number& operator*=(const double& e) - { - *this = *this * e; - return *this; - } - - Number& operator-=(const double& e) - { - *this = *this - e; - return *this; - } - - Number& operator/=(const double& e) - { - *this = *this / e; - return *this; - } -}; - diff --git a/AADNode.h b/AADNode.h deleted file mode 100644 index 4aa6a44..0000000 --- a/AADNode.h +++ /dev/null @@ -1,103 +0,0 @@ - -/* -Written by Antoine Savine in 2018 - -This code is the strict IP of Antoine Savine - -License to use and alter this code for personal and commercial applications -is freely granted to any person or company who purchased a copy of the book - -Modern Computational Finance: AAD and Parallel Simulations -Antoine Savine -Wiley, 2018 - -As long as this comment is preserved at the top of the file -*/ - -#pragma once - -// AAD implementation of chapter 10 -// (With multi-dimensional additions of chapter 14) - -// Implementation of Node = record on tape - -// Unchanged for AADET of chapter 15 - -#include -using namespace std; - -class Node -{ - friend class Tape; - friend class Number; - friend auto setNumResultsForAAD(const bool, const size_t); - friend struct numResultsResetterForAAD; - - // The adjoint(s) - // in single case, self held (chapter 10) - double mAdjoint = 0; - // in multi case, held separately and accessed by pointer (chapter 14) - double* pAdjoints; - - // Data lives in separate memory - - // the n derivatives to arguments, - double* pDerivatives; - - // the n pointers to the adjoints of arguments - double** pAdjPtrs; - - // Number of adjoints (results) to propagate, usually 1 - // See chapter 14 - static size_t numAdj; - - // Number of childs (arguments) - const size_t n; - -public: - - Node(const size_t N = 0) : n(N) {} - - // Access to adjoint(s) - // single - double& adjoint() - { - return mAdjoint; - } - // multi - double& adjoint(const size_t n) { return pAdjoints[n]; } - - // Back-propagate adjoints to arguments adjoints - - // Single case, chapter 10 - void propagateOne() -{ - // Nothing to propagate - if (!n || !mAdjoint) return; - - for (size_t i = 0; i < n; ++i) - { - *(pAdjPtrs[i]) += pDerivatives[i] * mAdjoint; - } - } - - // Multi case, chapter 14 - void propagateAll() -{ - // No adjoint to propagate - if (!n || all_of(pAdjoints, pAdjoints + numAdj, - [](const double& x) { return !x; })) - return; - - for (size_t i = 0; i < n; ++i) - { - double *adjPtrs = pAdjPtrs[i], ders = pDerivatives[i]; - - // Vectorized! - for (size_t j = 0; j < numAdj; ++j) - { - adjPtrs[j] += ders * pAdjoints[j]; - } - } - } -}; \ No newline at end of file diff --git a/AADNumber.h b/AADNumber.h deleted file mode 100644 index a1e6dc5..0000000 --- a/AADNumber.h +++ /dev/null @@ -1,677 +0,0 @@ - -/* -Written by Antoine Savine in 2018 - -This code is the strict IP of Antoine Savine - -License to use and alter this code for personal and commercial applications -is freely granted to any person or company who purchased a copy of the book - -Modern Computational Finance: AAD and Parallel Simulations -Antoine Savine -Wiley, 2018 - -As long as this comment is preserved at the top of the file -*/ - -#pragma once - -// Traditional AAD implementation of chapter 10 -// (With multi-dimensional additions of chapter 14) - -// The custom number type - -#include -#include "AADTape.h" - -class Number -{ - // Value and node are the only data members - double myValue; - Node* myNode; - - // Create node on tape - template - void createNode() - { - myNode = tape->recordNode(); - } - - // Access node (friends only) - // Note const incorectness - Node& node() const - { - -#ifdef _DEBUG - - // Help identify errors when arguments are not on tape - - // Find node on tape - auto it = tape->find(myNode); - - // Not found - if (it == tape->end()) - { - throw runtime_error("Put a breakpoint here"); - } - -#endif - // Const incorrectness - return const_cast(*myNode); - } - - // Convenient access to node data for friends - - double& derivative() { return myNode->pDerivatives[0]; } - double& lDer() { return myNode->pDerivatives[0]; } - double& rDer() { return myNode->pDerivatives[1]; } - - double*& adjPtr() { return myNode->pAdjPtrs[0]; } - double*& leftAdj() { return myNode->pAdjPtrs[0]; } - double*& rightAdj() { return myNode->pAdjPtrs[1]; } - - // Private constructors for operator overloading - - // Unary - Number(Node& arg, const double val) : - myValue(val) - { - createNode<1>(); - - myNode->pAdjPtrs[0] = Tape::multi - ? arg.pAdjoints - : &arg.mAdjoint; - } - - // Binary - Number(Node& lhs, Node& rhs, const double val) : - myValue(val) - { - createNode<2>(); - - if (Tape::multi) - { - myNode->pAdjPtrs[0] = lhs.pAdjoints; - myNode->pAdjPtrs[1] = rhs.pAdjoints; - } - else - { - myNode->pAdjPtrs[0] = &lhs.mAdjoint; - myNode->pAdjPtrs[1] = &rhs.mAdjoint; - } - } - -public: - - // Static access to tape - static thread_local Tape* tape; - - // Public constructors for leaves - - Number() {} - - // Put on tape on construction - explicit Number(const double val) : - myValue(val) - { - createNode<0>(); - } - - // Put on tape on assignment - Number& operator=(const double val) - { - myValue = val; - createNode<0>(); - - return *this; - } - - // Explicitly put existing Number on tape - void putOnTape() - { - createNode<0>(); - } - - // Explicit coversion to double - explicit operator double& () { return myValue; } - explicit operator double() const { return myValue; } - - // Accessors: value and adjoint - - double& value() - { - return myValue; - } - double value() const - { - return myValue; - } - // Single dimensional - double& adjoint() - { - return myNode->adjoint(); - } - double adjoint() const - { - return myNode->adjoint(); - } - // Multi dimensional - double& adjoint(const size_t n) - { - return myNode->adjoint(n); - } - double adjoint(const size_t n) const - { - return myNode->adjoint(n); - } - - // Reset all adjoints on the tape - // note we don't use this method - void resetAdjoints() - { - tape->resetAdjoints(); - } - - // Propagation - - // Propagate adjoints - // from and to both INCLUSIVE - static void propagateAdjoints( - Tape::iterator propagateFrom, - Tape::iterator propagateTo) - { - auto it = propagateFrom; - while (it != propagateTo) - { - it->propagateOne(); - --it; - } - it->propagateOne(); - } - - // Convenient overloads - - // Set the adjoint on this node to 1, - // Then propagate from the node - void propagateAdjoints( - // We start on this number's node - Tape::iterator propagateTo) - { - // Set this adjoint to 1 - adjoint() = 1.0; - // Find node on tape - auto propagateFrom = tape->find(myNode); - propagateAdjoints(propagateFrom, propagateTo); - } - - // These 2 set the adjoint to 1 on this node - void propagateToStart() - { - propagateAdjoints(tape->begin()); - } - void propagateToMark() - { - propagateAdjoints(tape->markIt()); - } - - // This one only propagates - // Note: propagation starts at mark - 1 - static void propagateMarkToStart() - { - propagateAdjoints(prev(tape->markIt()), tape->begin()); - } - - // Multi dimensional case: - // Propagate adjoints from and to both INCLUSIVE - static void propagateAdjointsMulti( - Tape::iterator propagateFrom, - Tape::iterator propagateTo) - { - auto it = propagateFrom; - while (it != propagateTo) - { - it->propagateAll(); - --it; - } - it->propagateAll(); - } - - // Operator overloading - - inline friend Number operator+(const Number& lhs, const Number& rhs) - { - const double e = lhs.value() + rhs.value(); - // Eagerly evaluate and put on tape - Number result(lhs.node(), rhs.node(), e); - // Eagerly compute derivatives - result.lDer() = 1.0; - result.rDer() = 1.0; - - return result; - } - inline friend Number operator+(const Number& lhs, const double& rhs) - { - const double e = lhs.value() + rhs; - // Eagerly evaluate and put on tape - Number result(lhs.node(), e); - // Eagerly compute derivatives - result.derivative() = 1.0; - - return result; - - } - inline friend Number operator+(const double& lhs, const Number& rhs) - { - return rhs + lhs; - } - - inline friend Number operator-(const Number& lhs, const Number& rhs) - { - const double e = lhs.value() - rhs.value(); - // Eagerly evaluate and put on tape - Number result(lhs.node(), rhs.node(), e); - // Eagerly compute derivatives - result.lDer() = 1.0; - result.rDer() = -1.0; - - return result; - } - inline friend Number operator-(const Number& lhs, const double& rhs) - { - const double e = lhs.value() - rhs; - // Eagerly evaluate and put on tape - Number result(lhs.node(), e); - // Eagerly compute derivatives - result.derivative() = 1.0; - - return result; - - } - inline friend Number operator-(const double& lhs, const Number& rhs) - { - const double e = lhs - rhs.value(); - // Eagerly evaluate and put on tape - Number result(rhs.node(), e); - // Eagerly compute derivatives - result.derivative() = -1.0; - - return result; - } - - inline friend Number operator*(const Number& lhs, const Number& rhs) - { - const double e = lhs.value() * rhs.value(); - // Eagerly evaluate and put on tape - Number result(lhs.node(), rhs.node(), e); - // Eagerly compute derivatives - result.lDer() = rhs.value(); - result.rDer() = lhs.value(); - - return result; - } - inline friend Number operator*(const Number& lhs, const double& rhs) - { - const double e = lhs.value() * rhs; - // Eagerly evaluate and put on tape - Number result(lhs.node(), e); - // Eagerly compute derivatives - result.derivative() = rhs; - - return result; - - } - inline friend Number operator*(const double& lhs, const Number& rhs) - { - return rhs * lhs; - } - - inline friend Number operator/(const Number& lhs, const Number& rhs) - { - const double e = lhs.value() / rhs.value(); - // Eagerly evaluate and put on tape - Number result(lhs.node(), rhs.node(), e); - // Eagerly compute derivatives - const double invRhs = 1.0 / rhs.value(); - result.lDer() = invRhs; - result.rDer() = -lhs.value() * invRhs * invRhs; - - return result; - } - inline friend Number operator/(const Number& lhs, const double& rhs) - { - const double e = lhs.value() / rhs; - // Eagerly evaluate and put on tape - Number result(lhs.node(), e); - // Eagerly compute derivatives - result.derivative() = 1.0 / rhs; - - return result; - - } - inline friend Number operator/(const double& lhs, const Number& rhs) - { - const double e = lhs / rhs.value(); - // Eagerly evaluate and put on tape - Number result(rhs.node(), e); - // Eagerly compute derivatives - result.derivative() = -lhs / rhs.value() / rhs.value(); - - return result; - } - - inline friend Number pow(const Number& lhs, const Number& rhs) - { - const double e = pow(lhs.value(), rhs.value()); - // Eagerly evaluate and put on tape - Number result(lhs.node(), rhs.node(), e); - // Eagerly compute derivatives - result.lDer() = rhs.value() * e / lhs.value(); - result.rDer() = log(lhs.value()) * e; - - return result; - } - inline friend Number pow(const Number& lhs, const double& rhs) - { - const double e = pow(lhs.value(), rhs); - // Eagerly evaluate and put on tape - Number result(lhs.node(), e); - // Eagerly compute derivatives - result.derivative() = rhs * e / lhs.value(); - - return result; - } - inline friend Number pow(const double& lhs, const Number& rhs) - { - const double e = pow(lhs, rhs.value()); - // Eagerly evaluate and put on tape - Number result(rhs.node(), e); - // Eagerly compute derivatives - result.derivative() = log(lhs) * e; - - return result; - } - - inline friend Number max(const Number& lhs, const Number& rhs) - { - const bool lmax = lhs.value() > rhs.value(); - // Eagerly evaluate and put on tape - Number result(lhs.node(), rhs.node(), lmax? lhs.value() : rhs.value()); - // Eagerly compute derivatives - if (lmax) - { - result.lDer() = 1.0; - result.rDer() = 0.0; - } - else - { - result.lDer() = 0.0; - result.rDer() = 1.0; - } - - return result; - } - inline friend Number max(const Number& lhs, const double& rhs) - { - const bool lmax = lhs.value() > rhs; - // Eagerly evaluate and put on tape - Number result(lhs.node(), lmax ? lhs.value() : rhs); - // Eagerly compute derivatives - result.derivative() = lmax ? 1.0 : 0.0; - - return result; - } - inline friend Number max(const double& lhs, const Number& rhs) - { - const bool rmax = rhs.value() > lhs; - // Eagerly evaluate and put on tape - Number result(rhs.node(), rmax ? rhs.value() : lhs); - // Eagerly compute derivatives - result.derivative() = rmax ? 1.0 : 0.0; - - return result; - } - - inline friend Number min(const Number& lhs, const Number& rhs) - { - const bool lmin = lhs.value() < rhs.value(); - // Eagerly evaluate and put on tape - Number result(lhs.node(), rhs.node(), lmin ? lhs.value() : rhs.value()); - // Eagerly compute derivatives - if (lmin) - { - result.lDer() = 1.0; - result.rDer() = 0.0; - } - else - { - result.lDer() = 0.0; - result.rDer() = 1.0; - } - - return result; - } - inline friend Number min(const Number& lhs, const double& rhs) - { - const bool lmin = lhs.value() < rhs; - // Eagerly evaluate and put on tape - Number result(lhs.node(), lmin ? lhs.value() : rhs); - // Eagerly compute derivatives - result.derivative() = lmin ? 1.0 : 0.0; - - return result; - } - inline friend Number min(const double& lhs, const Number& rhs) - { - const bool rmin = rhs.value() < lhs; - // Eagerly evaluate and put on tape - Number result(rhs.node(), rmin ? rhs.value() : lhs); - // Eagerly compute derivatives - result.derivative() = rmin ? 1.0 : 0.0; - - return result; - } - - Number& operator+=(const Number& arg) - { - *this = *this + arg; - return *this; - } - Number& operator+=(const double& arg) - { - *this = *this + arg; - return *this; - } - - Number& operator-=(const Number& arg) - { - *this = *this - arg; - return *this; - } - Number& operator-=(const double& arg) - { - *this = *this - arg; - return *this; - } - - Number& operator*=(const Number& arg) - { - *this = *this * arg; - return *this; - } - Number& operator*=(const double& arg) - { - *this = *this * arg; - return *this; - } - - Number& operator/=(const Number& arg) - { - *this = *this / arg; - return *this; - } - Number& operator/=(const double& arg) - { - *this = *this / arg; - return *this; - } - - // Unary +/- - Number operator-() const - { - return 0.0 - *this; - } - Number operator+() const - { - return *this; - } - - // Overloading continued, unary functions - - inline friend Number exp(const Number& arg) - { - const double e = exp(arg.value()); - // Eagerly evaluate and put on tape - Number result(arg.node(), e); - // Eagerly compute derivatives - result.derivative() = e; - - return result; - } - - inline friend Number log(const Number& arg) - { - const double e = log(arg.value()); - // Eagerly evaluate and put on tape - Number result(arg.node(), e); - // Eagerly compute derivatives - result.derivative() = 1.0 / arg.value(); - - return result; - } - - inline friend Number sqrt(const Number& arg) - { - const double e = sqrt(arg.value()); - // Eagerly evaluate and put on tape - Number result(arg.node(), e); - // Eagerly compute derivatives - result.derivative() = 0.5 / e; - - return result; - } - - inline friend Number fabs(const Number& arg) - { - const double e = fabs(arg.value()); - // Eagerly evaluate and put on tape - Number result(arg.node(), e); - // Eagerly compute derivatives - result.derivative() = arg.value() > 0.0 ? 1.0 : -1.0; - - return result; - } - - inline friend Number normalDens(const Number& arg) - { - const double e = normalDens(arg.value()); - // Eagerly evaluate and put on tape - Number result(arg.node(), e); - // Eagerly compute derivatives - result.derivative() = - arg.value() * e; - - return result; - } - - inline friend Number normalCdf(const Number& arg) - { - const double e = normalCdf(arg.value()); - // Eagerly evaluate and put on tape - Number result(arg.node(), e); - // Eagerly compute derivatives - result.derivative() = normalDens(arg.value()); - - return result; - } - - // Finally, comparison - - inline friend bool operator==(const Number& lhs, const Number& rhs) - { - return lhs.value() == rhs.value(); - } - inline friend bool operator==(const Number& lhs, const double& rhs) - { - return lhs.value() == rhs; - } - inline friend bool operator==(const double& lhs, const Number& rhs) - { - return lhs == rhs.value(); - } - - inline friend bool operator!=(const Number& lhs, const Number& rhs) - { - return lhs.value() != rhs.value(); - } - inline friend bool operator!=(const Number& lhs, const double& rhs) - { - return lhs.value() != rhs; - } - inline friend bool operator!=(const double& lhs, const Number& rhs) - { - return lhs != rhs.value(); - } - - inline friend bool operator<(const Number& lhs, const Number& rhs) - { - return lhs.value() < rhs.value(); - } - inline friend bool operator<(const Number& lhs, const double& rhs) - { - return lhs.value() < rhs; - } - inline friend bool operator<(const double& lhs, const Number& rhs) - { - return lhs < rhs.value(); - } - - inline friend bool operator>(const Number& lhs, const Number& rhs) - { - return lhs.value() > rhs.value(); - } - inline friend bool operator>(const Number& lhs, const double& rhs) - { - return lhs.value() > rhs; - } - inline friend bool operator>(const double& lhs, const Number& rhs) - { - return lhs > rhs.value(); - } - - inline friend bool operator<=(const Number& lhs, const Number& rhs) - { - return lhs.value() <= rhs.value(); - } - inline friend bool operator<=(const Number& lhs, const double& rhs) - { - return lhs.value() <= rhs; - } - inline friend bool operator<=(const double& lhs, const Number& rhs) - { - return lhs <= rhs.value(); - } - - inline friend bool operator>=(const Number& lhs, const Number& rhs) - { - return lhs.value() >= rhs.value(); - } - inline friend bool operator>=(const Number& lhs, const double& rhs) - { - return lhs.value() >= rhs; - } - inline friend bool operator>=(const double& lhs, const Number& rhs) - { - return lhs >= rhs.value(); - } -}; - - diff --git a/AADTape.h b/AADTape.h deleted file mode 100644 index 844fdfa..0000000 --- a/AADTape.h +++ /dev/null @@ -1,181 +0,0 @@ - -/* -Written by Antoine Savine in 2018 - -This code is the strict IP of Antoine Savine - -License to use and alter this code for personal and commercial applications -is freely granted to any person or company who purchased a copy of the book - -Modern Computational Finance: AAD and Parallel Simulations -Antoine Savine -Wiley, 2018 - -As long as this comment is preserved at the top of the file -*/ - -#pragma once - -// AAD implementation of chapter 10 -// (With multi-dimensional additions of chapter 14) - -// Implementation of the Tape - -// Unchanged for AADET of chapter 15 - -#include "blocklist.h" -#include "AADNode.h" - -constexpr size_t BLOCKSIZE = 16384; // Number of nodes -constexpr size_t ADJSIZE = 32768; // Number of adjoints -constexpr size_t DATASIZE = 65536; // Data in bytes - -class Tape -{ - // Working with multiple results / adjoints? - static bool multi; - - // Storage for adjoints in multi-dimensional case (chapter 14) - blocklist myAdjointsMulti; - - // Storage for derivatives and child adjoint pointers - blocklist myDers; - blocklist myArgPtrs; - - // Storage for the nodes - blocklist myNodes; - - // Padding so tapes in a vector don't interfere - char myPad[64]; - - friend auto setNumResultsForAAD(const bool, const size_t); - friend struct numResultsResetterForAAD; - friend class Number; - -public: - - // Build note in place and return a pointer - // N : number of childs (arguments) - template - Node* recordNode() - { - // Construct the node in place on tape - Node* node = myNodes.emplace_back(N); - - // Store and zero the adjoint(s) - if (multi) - { - node->pAdjoints = myAdjointsMulti.emplace_back_multi(Node::numAdj); - fill(node->pAdjoints, node->pAdjoints + Node::numAdj, 0.0); - } - - // Store the derivatives and child adjoint pointers unless leaf - if constexpr(N>0) - { - node->pDerivatives = myDers.emplace_back_multi(); - node->pAdjPtrs = myArgPtrs.emplace_back_multi(); - - } - - return node; - } - - // Reset all adjoints to 0 - void resetAdjoints() - { - if (multi) - { - myAdjointsMulti.memset(0); - } - else - { - for (Node& node : myNodes) - { - node.mAdjoint = 0; - } - } - } - - // Clear - void clear() - { - myAdjointsMulti.clear(); - myDers.clear(); - myArgPtrs.clear(); - myNodes.clear(); - } - - // Rewind - void rewind() - { - -#ifdef _DEBUG - - // In debug mode, always wipe - // makes it easier to identify errors - - clear(); - -#else - // In release mode, rewind and reuse - - if (multi) - { - myAdjointsMulti.rewind(); - } - myDers.rewind(); - myArgPtrs.rewind(); - myNodes.rewind(); - -#endif - - } - - // Set mark - void mark() - { - if (multi) - { - myAdjointsMulti.setmark(); - } - myDers.setmark(); - myArgPtrs.setmark(); - myNodes.setmark(); - } - - // Rewind to mark - void rewindToMark() - { - if (multi) - { - myAdjointsMulti.rewind_to_mark(); - } - myDers.rewind_to_mark(); - myArgPtrs.rewind_to_mark(); - myNodes.rewind_to_mark(); - } - - // Iterators - - using iterator = blocklist::iterator; - - auto begin() - { - return myNodes.begin(); - } - - auto end() - { - return myNodes.end(); - } - - auto markIt() - { - return myNodes.mark(); - } - - auto find(Node* node) - { - return myNodes.find(node); - } -}; \ No newline at end of file diff --git a/AutocallPricer.xlsx b/AutocallPricer.xlsx deleted file mode 100644 index e21e212..0000000 Binary files a/AutocallPricer.xlsx and /dev/null differ diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..c840515 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,42 @@ +cmake_minimum_required(VERSION 3.30) + +project(AAD_book VERSION 1.0.0 LANGUAGES C CXX) + +set(CMAKE_CXX_STANDARD 20) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSIONS OFF) +set(CMAKE_EXPORT_COMPILE_COMMANDS ON) + +if(WIN32) + message(STATUS "Windows detected") + message(STATUS "Toolset: ${CMAKE_GENERATOR_TOOLSET}") + set(dlloader_include_dir_platform ${CMAKE_SOURCE_DIR}/dl_loader/windows) + add_compile_definitions(WIN_EXPORT) + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -DNDEBUG /openmp /Oi") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_RELEASE} -D_DEBUG /Od") +endif() + +if(UNIX) + message(STATUS "Unix system detected") + set(dlloader_include_dir_platform ${CMAKE_SOURCE_DIR}/dl_loader/linux) +endif() + +set(dlloader_include_dir ${CMAKE_SOURCE_DIR}/dl_loader) + + +MESSAGE(STATUS "the source directory is: ${CMAKE_SOURCE_DIR}") +MESSAGE(STATUS "the binary directory is: ${CMAKE_BINARY_DIR}") + +set(CMAKE_INSTALL_PREFIX ${CMAKE_SOURCE_DIR}/install) + +set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin) +set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin) +set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin) + +set(CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/") + +if(ENABLE_WARNINGS) + include(Warnings) +endif() + +add_subdirectory(libraries) diff --git a/ConcurrentQueue.h b/ConcurrentQueue.h deleted file mode 100644 index 1edc950..0000000 --- a/ConcurrentQueue.h +++ /dev/null @@ -1,101 +0,0 @@ -#pragma once - -// Concurrent queue of chapter 3, -// Used in the thread pool - -#include -#include -using namespace std; - -template -class ConcurrentQueue -{ - - queue myQueue; - mutable mutex myMutex; - condition_variable myCV; - bool myInterrupt; - -public: - - ConcurrentQueue() : myInterrupt(false) {} - ~ConcurrentQueue() { interrupt(); } - - bool empty() const - { - // Lock - lock_guard lk(myMutex); - // Access underlying queue - return myQueue.empty(); - } // Unlock - - // Pop into argument - bool tryPop(T& t) - { - // Lock - lock_guard lk(myMutex); - if (myQueue.empty()) return false; - // Move from queue - t = move(myQueue.front()); - // Combine front/pop - myQueue.pop(); - - return true; - } // Unlock - - // Pass t byVal or move with push( move( t)) - void push(T t) - { - { - // Lock - lock_guard lk(myMutex); - // Move into queue - myQueue.push(move(t)); - } // Unlock before notification - - // Unlock before notification - myCV.notify_one(); - } - - // Wait if empty - bool pop(T& t) - { - // (Unique) lock - unique_lock lk(myMutex); - - // Wait if empty, release lock until notified - while (!myInterrupt && myQueue.empty()) myCV.wait(lk); - - // Re-acquire lock, resume - - // Check for interruption - if (myInterrupt) return false; - - // Combine front/pop - t = move(myQueue.front()); - myQueue.pop(); - - return true; - - } // Unlock - - void interrupt() - { - { - lock_guard lk(myMutex); - myInterrupt = true; - } - myCV.notify_all(); - } - - void resetInterrupt() - { - myInterrupt = false; - } - - void clear() - { - queue empty; - swap(myQueue, empty); - } -}; \ No newline at end of file diff --git a/NOTES.md b/NOTES.md new file mode 100644 index 0000000..1df812b --- /dev/null +++ b/NOTES.md @@ -0,0 +1,11 @@ +# Some Notes + +Some useful commands, assuming we are in the root of the project: + +For Windows: + +1. mkdir build +2. cd build && cmake .. -G "Visual Studio 17 2022" -A x64 +3. cmake --build . --config Release +4. cmake --build . --config Debug +5. cmake --install . -> this needs a pre setup if you are not admin, you could hit some errors. diff --git a/README.md b/README.md index 0f084d6..b1dc229 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # Companion code for "Modern Computational Finance: AAD and Parallel Simulations" (Antoine Savine, Wiley, 2018) -![Screenshot](back.jpg) +![Screenshot](files/back.jpg) This is the professional implementation in C++ of the book Modern Computational Finance: AAD and Parallel Simulations by Antoine Savine. The code is freely available to anyone. Any person who purchased a copy of the book is authorized to use, modify and distribute the code for any application, as long as the credits remain on the top of the files. diff --git a/analytics.h b/analytics.h deleted file mode 100644 index 3c230f5..0000000 --- a/analytics.h +++ /dev/null @@ -1,123 +0,0 @@ -#pragma once - -#include "gaussians.h" - -// Black and Scholes formula, templated - -// Price -template -inline T blackScholes( - const U spot, - const V strike, - const T vol, - const W mat) -{ - const auto std = vol * sqrt(mat); - if (std <= EPS) return max(T(0.0), T(spot - strike)); - const auto d2 = log(spot / strike) / std - 0.5 * std; - const auto d1 = d2 + std; - return spot * normalCdf(d1) - strike * normalCdf(d2); -} - -// Implied vol, untemplated -inline double blackScholesIvol( - const double spot, - const double strike, - const double prem, - const double mat) -{ - if (prem <= max(0.0, spot - strike) + EPS) return 0.0; - - double p, pu, pl; - double u = 0.5; - while (blackScholes(spot, strike, u, mat) < prem) u *= 2; - double l = 0.05; - while (blackScholes(spot, strike, l, mat) > prem) l /= 2; - pu = blackScholes(spot, strike, u, mat); - blackScholes(spot, strike, l, mat); - - while (u - l > 1.e-12) - { - const double m = 0.5 * (u + l); - p = blackScholes(spot, strike, m, mat); - if (p > prem) - { - u = m; - pu = p; - } - else - { - l = m; - pl = p; - } - } - - return l + (prem - pl) / (pu - pl) * (u - l); -} - -// Merton, templated - -template -inline T merton( - const U spot, - const V strike, - const T vol, - const W mat, - const X intens, - const X meanJmp, - const X stdJmp) -{ - const auto varJmp = stdJmp * stdJmp; - const auto mv2 = meanJmp + 0.5 * varJmp; - const auto comp = intens * (exp(mv2) - 1); - const auto var = vol * vol; - const auto intensT = intens * mat; - - unsigned fact = 1; - X iT = 1.0; - const size_t cut = 10; - T result = 0.0; - for (size_t n = 0; n < cut; ++n) - { - const auto s = spot*exp(n*mv2 - comp*mat); - const auto v = sqrt(var + n * varJmp / mat); - const auto prob = exp(-intensT) * iT / fact; - result += prob * blackScholes(s, strike, v, mat); - fact *= n + 1; - iT *= intensT; - } - - return result; -} - -// Up and out call in Black-Scholes, untemplated - -inline double BlackScholesKO( - const double spot, - const double rate, - const double div, - const double strike, - const double barrier, - const double mat, - const double vol) -{ - const double std = vol * sqrt(mat); - const double fwdFact = exp((rate - div)*mat); - const double fwd = spot * fwdFact; - const double disc = exp(-rate * mat); - const double v = rate - div - 0.5 * vol * vol; - const double d2 = (log(spot / barrier) + v * mat) / std; - const double d2prime = (log(barrier / spot) + v * mat) / std; - - const double bar = blackScholes(fwd, strike, vol, mat) - - blackScholes(fwd, barrier, vol, mat) - - (barrier - strike) * normalCdf(d2) - - pow(barrier / spot, 2 * v / vol / vol) * - ( - blackScholes(fwdFact * barrier* barrier / spot, strike, vol, mat) - - blackScholes(fwdFact * barrier* barrier / spot, barrier, vol, mat) - - (barrier - strike) * normalCdf(d2prime) - ); - - return disc * bar; -} diff --git a/blocklist.h b/blocklist.h deleted file mode 100644 index f18bb69..0000000 --- a/blocklist.h +++ /dev/null @@ -1,339 +0,0 @@ - -/* -Written by Antoine Savine in 2018 - -This code is the strict IP of Antoine Savine - -License to use and alter this code for personal and commercial applications -is freely granted to any person or company who purchased a copy of the book - -Modern Computational Finance: AAD and Parallel Simulations -Antoine Savine -Wiley, 2018 - -As long as this comment is preserved at the top of the file -*/ - -#pragma once -#pragma warning(disable : 4018) - -// Blocklist data structure for AAD memory management -// See chapter 10, unchanged with expression templates of chapter 15 - -#include -#include -#include -using namespace std; - -template -class blocklist -{ - // Container = list of blocks - list> data; - - using list_iter = decltype(data.begin()); - using block_iter = decltype(data.back().begin()); - - // Current block - list_iter cur_block; - - // Last block - list_iter last_block; - - // Next free space in current block - block_iter next_space; - - // Last free space (+1) in current block - block_iter last_space; - - // Mark - list_iter marked_block; - block_iter marked_space; - - // Create new array - void newblock() - { - data.emplace_back(); - cur_block = last_block = prev(data.end()); - next_space = cur_block->begin(); - last_space = cur_block->end(); - } - - // Move on to next array - void nextblock() - { - // This is the last array: create new - if (cur_block == last_block) - { - newblock(); - } - else - { - ++cur_block; - next_space = cur_block->begin(); - last_space = cur_block->end(); - } - } - -public: - - // Create first block on construction - blocklist() - { - newblock(); - } - - // Factory reset - void clear() - { - data.clear(); - newblock(); - } - - // Rewind but keep all blocks - void rewind() - { - cur_block = data.begin(); - next_space = cur_block->begin(); - last_space = cur_block->end(); - } - - // Memset - void memset(unsigned char value = 0) - { - for (auto& arr : data) - { - std::memset(&arr[0], value, block_size * sizeof(T)); - } - } - - // Construct object of type T in place - // in the next free space and return a pointer on it - // Implements perfect forwarding of constructor arguments - template - T* emplace_back(Args&& ...args) - { - // No more space in current array - if (next_space == last_space) - { - nextblock(); - } - // Placement new, construct in memory pointed by next - T* emplaced = new (&*next_space) // memory pointed by next as T* - T(forward(args)...); // perfect forwarding of ctor arguments - - // Advance next - ++next_space; - - // Return - return emplaced; - } - - // Overload for default constructed - T* emplace_back() - { - // No more space in current array - if (next_space == last_space) - { - nextblock(); - } - - // Current space - auto old_next = next_space; - - // Advance next - ++next_space; - - // Return - return &*old_next; - } - - // Stores n default constructed elements - // and returns a pointer on the first - - // Version 1: n known at compile time - template - T* emplace_back_multi() - { - // No more space in current array - if (distance(next_space, last_space) < n) - { - nextblock(); - } - - // Current space - auto old_next = next_space; - - // Advance next - next_space += n; - - // Return - return &*old_next; - } - - // Version 2: n unknown at compile time - T* emplace_back_multi(const size_t n) - { - // No more space in current array - if (distance(next_space, last_space) < n) - { - nextblock(); - } - - // Current space - auto old_next = next_space; - - // Advance next - next_space += n; - - // Return - return &*old_next; - } - - // Marks - - // Set mark - void setmark() - { - if (next_space == last_space) - { - nextblock(); - } - - marked_block = cur_block; - marked_space = next_space; - } - - // Rewind to mark - void rewind_to_mark() - { - cur_block = marked_block; - next_space = marked_space; - last_space = cur_block->end(); - } - - // Iterator - - class iterator - { - // List and block - list_iter cur_block; // current block - block_iter cur_space; // current space - block_iter first_space; // first space in block - block_iter last_space; // last (+1) space in block - - public: - - // iterator traits - using difference_type = ptrdiff_t; - using reference = T&; - using pointer = T*; - using value_type = T; - using iterator_category = bidirectional_iterator_tag; - - // Default constructor - iterator() {} - - // Constructor - iterator(list_iter cb, block_iter cs, block_iter fs, block_iter ls) : - cur_block(cb), cur_space(cs), first_space(fs), last_space(ls) {} - - // Pre-increment (we do not provide post) - iterator& operator++() - { - ++cur_space; - if (cur_space == last_space) - { - ++cur_block; - first_space = cur_block->begin(); - last_space = cur_block->end(); - cur_space = first_space; - } - - return *this; - } - - // Pre-decrement - iterator& operator--() - { - if (cur_space == first_space) - { - --cur_block; - first_space = cur_block->begin(); - last_space = cur_block->end(); - cur_space = last_space; - } - - --cur_space; - - return *this; - } - - // Access to contained elements - T& operator*() - { - return *cur_space; - } - const T& operator*() const - { - return *cur_space; - } - T* operator->() - { - return &*cur_space; - } - const T* operator->() const - { - return &*cur_space; - } - - // Check equality - bool operator ==(const iterator& rhs) const - { - return (cur_block == rhs.cur_block && cur_space == rhs.cur_space); - } - bool operator !=(const iterator& rhs) const - { - return (cur_block != rhs.cur_block|| cur_space != rhs.cur_space); - } - }; - - // Access to iterators - - iterator begin() - { - return iterator(data.begin(), data.begin()->begin(), - data.begin()->begin(), data.begin()->end()); - } - - iterator end() - { - return iterator(cur_block, next_space, - cur_block->begin(), cur_block->end()); - } - - // Iterator on mark - iterator mark() - { - return iterator(marked_block, marked_space, - marked_block->begin(), marked_block->end()); - } - - // Find element, by pointer, searching sequentially from the end - iterator find(const T* const element) - { - // Search from the end - iterator it = end(); - iterator b = begin(); - - while (it != b) - { - --it; - if (&*it == element) return it; - } - - if (&*it == element) return it; - - return end(); - } -}; \ No newline at end of file diff --git a/choldc.h b/choldc.h deleted file mode 100644 index 4fbaa53..0000000 --- a/choldc.h +++ /dev/null @@ -1,50 +0,0 @@ -#pragma once - -// Choleski decomposition, inspired by Numerical Recipes - -#include "matrix.h" - -template -void choldc(const matrix& in, matrix& out) -{ - int n = in.rows(); - T sum; - - fill(out.begin(), out.end(), T(0.0)); - - for(int i=0; i + +namespace dlloader { + +template +class IDLLoader { + public: + virtual ~IDLLoader() = default; + + virtual void DLOpenLib() = 0; + virtual std::shared_ptr DLGetInstance() = 0; + virtual void DLCloseLib() = 0; +}; +} // namespace dlloader + +#endif // IDL_LOADER_HPP \ No newline at end of file diff --git a/dl_loader/dl_loader b/dl_loader/dl_loader new file mode 100644 index 0000000..69be950 --- /dev/null +++ b/dl_loader/dl_loader @@ -0,0 +1,16 @@ +#pragma once + +#ifndef DL_LOADER_INCLUDES +#define DL_LOADER_INCLUDES + +#ifdef WIN32 +#include "./windows/DLLoader.hpp" +#elif defined(__linux__) || defined(__APPLE__) +#include "./linux/DLLoader.hpp" +#else +#endif + +#include "./IDLLoader.hpp" +#include "./library_names.hpp" + +#endif // DL_LOADER_INCLUDES \ No newline at end of file diff --git a/dl_loader/library_names.hpp b/dl_loader/library_names.hpp new file mode 100644 index 0000000..3deb707 --- /dev/null +++ b/dl_loader/library_names.hpp @@ -0,0 +1,19 @@ +#pragma once + +#ifndef DL_LOADER_LIBRARY_NAMES_HPP +#define DL_LOADER_LIBRARY_NAMES_HPP +#include + +#ifdef __APPLE__ +static const std::string LIBRARY_NAME = "toupdate"; +#endif + +#ifdef __linux__ +static const std::string LIBRARY_NAME = "toupdate"; +#endif + +#ifdef WIN32 +static const std::string LIBRARY_NAME = "to_update"; +#endif + +#endif // DL_LOADER_LIBRARY_NAMES_HPP \ No newline at end of file diff --git a/dl_loader/linux/DLLoader.hpp b/dl_loader/linux/DLLoader.hpp new file mode 100644 index 0000000..a24ef95 --- /dev/null +++ b/dl_loader/linux/DLLoader.hpp @@ -0,0 +1,68 @@ +#pragma once +#ifndef IDL_LOADER_UNIX_HPP +#define IDL_LOADER_UNIX_HPP + +#if defined(__linux__) || defined(__APPLE__) +#include + +#include +#include + +#include "../IDLLoader.hpp" + +namespace dlloader { + +template +class DLLoaderUnix : public IDLLoader { + void *handle; + std::string pathToLib; + std::string allocClassSymbol; + std::string deallocClassSymbol; + + public: + DLLoaderUnix(const std::string &pathToLib, + const std::string &allocClassSymbol, + const std::string &deallocClassSymbol) + : handle(nullptr), + pathToLib(pathToLib), + allocClassSymbol(allocClassSymbol), + deallocClassSymbol(deallocClassSymbol) {} + + ~DLLoaderUnix() = default; + + void DLOpenLib() override { + handle = dlopen(pathToLib.c_str(), RTLD_NOW | RTLD_LAZY); + if (!handle) { + // throw std::runtime_error(dlerror()); + std::cerr << "Error opening library: " << dlerror() << std::endl; + } + } + + std::shared_ptr DLGetInstance() override { + using allocClass = T *(*)(); + using deallocClass = void (*)(T *); + + auto allocFunc = + reinterpret_cast(dlsym(handle, allocClassSymbol.c_str())); + + auto deleteFunc = reinterpret_cast( + dlsym(handle, deallocClassSymbol.c_str())); + if (!allocFunc || !deleteFunc) { + DLCloseLib(); + std::cerr << "Error loading symbols: " << dlerror() << std::endl; + } + + return std::shared_ptr(allocFunc(), + [deleteFunc](T *ptr) { deleteFunc(ptr); }); + } + + void DLCloseLib() override { + if (dlclose(handle) != 0) { + std::cerr << "Error closing library: " << dlerror() << std::endl; + } + } +}; +} // namespace dlloader +#endif // defined(__linux__) || defined(__APPLE__) + +#endif // IDL_LOADER_UNIX_HPP \ No newline at end of file diff --git a/dl_loader/windows/DLLoader.hpp b/dl_loader/windows/DLLoader.hpp new file mode 100644 index 0000000..20930c5 --- /dev/null +++ b/dl_loader/windows/DLLoader.hpp @@ -0,0 +1,67 @@ +#pragma once +#ifndef IDL_LOADER_WIN_HPP +#define IDL_LOADER_WIN_HPP + +#ifdef WIN32 + +#include + +#include "../IDLLoader.hpp" + +namespace dlloader { + +template +class DLLoader : public IDLLoader { + HMODULE handle; + std::string pathToLib; + std::string allocClassSymbol; + std::string deallocClassSymbol; + + public: + DLLoader(const std::string &pathToLib, const std::string &allocClassSymbol, + const std::string &deallocClassSymbol) + : handle(nullptr), + pathToLib(pathToLib), + allocClassSymbol(allocClassSymbol), + deallocClassSymbol(deallocClassSymbol) {} + + ~DLLoader() = default; + + void DLOpenLib() override { + handle = LoadLibrary(pathToLib.c_str()); + + if (!handle) { + std::cerr << "Error opening library: " << pathToLib << std::endl; + } + } + + std::shared_ptr DLGetInstance() override { + using allocClass = T *(*)(); + using deallocClass = void (*)(T *); + + auto allocFunc = reinterpret_cast( + GetProcAddress(handle, allocClassSymbol.c_str())); + + auto deleteFunc = reinterpret_cast( + GetProcAddress(handle, deallocClassSymbol.c_str())); + + if (!allocFunc || !deleteFunc) { + DLCloseLib(); + std::cerr << "Error loading symbols: " << pathToLib << std::endl; + } + + return std::shared_ptr(allocFunc(), + [deleteFunc](T *ptr) { deleteFunc(ptr); }); + } + + void DLCloseLib() override { + if (FreeLibrary(handle) == 0) { + std::cerr << "Error closing library: " << pathToLib << std::endl; + } + } +}; +} // namespace dlloader + +#endif // _WIN32 + +#endif // IDL_LOADER_WIN_HPP \ No newline at end of file diff --git a/Intro2AADinMachineLearningAndFinance.pdf b/files/Intro2AADinMachineLearningAndFinance.pdf similarity index 100% rename from Intro2AADinMachineLearningAndFinance.pdf rename to files/Intro2AADinMachineLearningAndFinance.pdf diff --git a/VC_redist.x64.exe b/files/VC_redist.x64.exe similarity index 100% rename from VC_redist.x64.exe rename to files/VC_redist.x64.exe diff --git a/VC_redist.x86.exe b/files/VC_redist.x86.exe similarity index 100% rename from VC_redist.x86.exe rename to files/VC_redist.x86.exe diff --git a/back.jpg b/files/back.jpg similarity index 100% rename from back.jpg rename to files/back.jpg diff --git a/xlComp.sln b/files/xlComp.sln similarity index 100% rename from xlComp.sln rename to files/xlComp.sln diff --git a/xlComp.vcxproj b/files/xlComp.vcxproj similarity index 98% rename from xlComp.vcxproj rename to files/xlComp.vcxproj index fc4d81c..71475a4 100644 --- a/xlComp.vcxproj +++ b/files/xlComp.vcxproj @@ -41,13 +41,13 @@ DynamicLibrary true - v142 + v143 Unicode DynamicLibrary false - v142 + v143 true Unicode @@ -204,6 +204,7 @@ + diff --git a/xlComp.vcxproj.filters b/files/xlComp.vcxproj.filters similarity index 97% rename from xlComp.vcxproj.filters rename to files/xlComp.vcxproj.filters index 33a3801..f15d638 100644 --- a/xlComp.vcxproj.filters +++ b/files/xlComp.vcxproj.filters @@ -133,5 +133,8 @@ Source Files + + Source Files + \ No newline at end of file diff --git a/files/xlComp.xll b/files/xlComp.xll new file mode 100644 index 0000000..56fed1c Binary files /dev/null and b/files/xlComp.xll differ diff --git a/gaussians.h b/gaussians.h deleted file mode 100644 index aee6d2d..0000000 --- a/gaussians.h +++ /dev/null @@ -1,87 +0,0 @@ -#pragma once - -#include -#include -#include -using namespace std; - -#define EPS 1.0e-08 - -// Gaussian functions - -// Normal density -inline double normalDens(const double x) -{ - return x<-10.0 || 10.0 10.0) return 1.0; - if (x < 0.0) return 1.0 - normalCdf(-x); - - static constexpr double p = 0.2316419; - static constexpr double b1 = 0.319381530; - static constexpr double b2 = -0.356563782; - static constexpr double b3 = 1.781477937; - static constexpr double b4 = -1.821255978; - static constexpr double b5 = 1.330274429; - - const auto t = 1.0 / (1.0 + p*x); - - const auto pol = t*(b1 + t*(b2 + t*(b3 + t*(b4 + t*b5)))); - - const auto pdf = normalDens(x); - - return 1.0 - pdf * pol; -} - -// Inverse CDF (for generation of Gaussians out of Uniforms) -// Beasley-Springer-Moro algorithm -// Moro, The full Monte, Risk, 1995 -// See Glasserman, Monte Carlo Methods in Financial Engineering, p 68 -inline double invNormalCdf(const double p) -{ - const bool sup = p > 0.5; - const double up = sup ? 1.0 - p : p; - - static constexpr double a0 = 2.50662823884; - static constexpr double a1 = -18.61500062529; - static constexpr double a2 = 41.39119773534; - static constexpr double a3 = -25.44106049637; - - static constexpr double b0 = -8.47351093090; - static constexpr double b1 = 23.08336743743; - static constexpr double b2 = -21.06224101826; - static constexpr double b3 = 3.13082909833; - - static constexpr double c0 = 0.3374754822726147; - static constexpr double c1 = 0.9761690190917186; - static constexpr double c2 = 0.1607979714918209; - static constexpr double c3 = 0.0276438810333863; - static constexpr double c4 = 0.0038405729373609; - static constexpr double c5 = 0.0003951896511919; - static constexpr double c6 = 0.0000321767881768; - static constexpr double c7 = 0.0000002888167364; - static constexpr double c8 = 0.0000003960315187; - - double x = up - 0.5; - double r; - - if (fabs(x)<0.42) - { - r = x*x; - r = x*(((a3*r + a2)*r + a1)*r + a0) / ((((b3*r + b2)*r + b1)*r + b0)*r + 1.0); - return sup ? -r: r; - } - - r = up; - r = log(-log(r)); - r = c0 + r*(c1 + r*(c2 + r*(c3 + r*(c4 + r*(c5 + r*(c6 + r*(c7 + r*c8))))))); - - return sup? r: -r; -} diff --git a/interp.h b/interp.h deleted file mode 100644 index ae4a978..0000000 --- a/interp.h +++ /dev/null @@ -1,109 +0,0 @@ - -/* -Written by Antoine Savine in 2018 - -This code is the strict IP of Antoine Savine - -License to use and alter this code for personal and commercial applications -is freely granted to any person or company who purchased a copy of the book - -Modern Computational Finance: AAD and Parallel Simulations -Antoine Savine -Wiley, 2018 - -As long as this comment is preserved at the top of the file -*/ - -#pragma once - -#include -using namespace std; - -// Interpolation utility - -// Interpolates the vector y against knots x in value x0 -// Interpolation is linear or smooth, extrapolation is flat -template -inline auto interp( - // sorted on xs - ITX xBegin, - ITX xEnd, - // corresponding ys - ITY yBegin, - ITY yEnd, - // interpolate for point x0 - const T& x0) - ->remove_reference_t -{ - // STL binary search, returns iterator on 1st no less than x0 - // upper_bound guarantees logarithmic search - auto it = upper_bound(xBegin, xEnd, x0); - - // Extrapolation? - if (it == xEnd) return *(yEnd - 1); - if (it == xBegin) return *yBegin; - - // Interpolation - size_t n = distance(xBegin, it) - 1; - auto x1 = xBegin[n]; - auto y1 = yBegin[n]; - auto x2 = xBegin[n + 1]; - auto y2 = yBegin[n + 1]; - - auto t = (x0 - x1) / (x2 - x1); - - if constexpr (smoothStep) - { - return y1 + (y2 - y1) * t * t * (3.0 - 2 * t); - } - else - { - return y1 + (y2 - y1) * t; - } -} - -// 2D -template -inline V interp2D( - // sorted on xs - const vector& x, - // sorted on ys - const vector& y, - // zs in a matrix - const matrix& z, - // interpolate for point (x0,y0) - const W& x0, - const X& y0) -{ - const size_t n = x.size(); - const size_t m = y.size(); - - // STL binary search, returns iterator on 1st no less than x0 - // upper_boung guarantees logarithmic search - auto it = upper_bound(x.begin(), x.end(), x0); - const size_t n2 = distance(x.begin(), it); - - // Extrapolation in x? - if (n2 == n) - return interp(y.begin(), y.end(), z[n2 - 1], z[n2 - 1] + m, y0); - if (n2 == 0) - return interp(y.begin(), y.end(), z[0], z[0] + m, y0); - - // Interpolation in x - const size_t n1 = n2 - 1; - auto x1 = x[n1]; - auto x2 = x[n2]; - auto z1 = interp(y.begin(), y.end(), z[n1], z[n1] + m, y0); - auto z2 = interp(y.begin(), y.end(), z[n2], z[n2] + m, y0); - - // Smooth step - auto t = (x0 - x1) / (x2 - x1); - if constexpr (smoothStep) - { - return z1 + (z2 - z1) * t * t * (3.0 - 2 * t); - } - else - { - return z1 + (z2 - z1) * t;; - } -} diff --git a/libraries/CMakeLists.txt b/libraries/CMakeLists.txt new file mode 100644 index 0000000..7455e62 --- /dev/null +++ b/libraries/CMakeLists.txt @@ -0,0 +1,14 @@ +project(libraries) + +set(common_includes_dir ${CMAKE_CURRENT_SOURCE_DIR}/common/) + +add_subdirectory(aad) +set(aad_dir ${CMAKE_CURRENT_SOURCE_DIR}/aad) + +if(WIN32) + message(STATUS "Windows detected - xll interfaces") + set(excel_sdk_api_dir ${CMAKE_CURRENT_SOURCE_DIR}/excel_sdk_api) + + add_subdirectory(excel_sdk_api) + add_subdirectory(aad_xl_interface) +endif() diff --git a/libraries/aad/CMakeLists.txt b/libraries/aad/CMakeLists.txt new file mode 100644 index 0000000..44f84f7 --- /dev/null +++ b/libraries/aad/CMakeLists.txt @@ -0,0 +1,14 @@ +project(aad) + +set(includes includes/) + +set(sources src/AAD.cpp + src/mcBase.cpp + src/sobol.cpp + src/ThreadPool.cpp +) + +# assuming at the moment that we are only msvc + +add_library(${PROJECT_NAME} ${sources}) +target_include_directories(${PROJECT_NAME} PUBLIC ${includes} ${common_includes_dir}) diff --git a/libraries/aad/aad b/libraries/aad/aad new file mode 100644 index 0000000..ef8f23a --- /dev/null +++ b/libraries/aad/aad @@ -0,0 +1,8 @@ +#ifndef AAD_LIBRARY +#define AAD_LIBRARY + +#include "./includes/AAD.h" // IWYU pragma: keep +#include "./includes/ThreadPool.h" // IWYU pragma: keep +#include "./includes/mcBase.h" // IWYU pragma: keep + +#endif // AAD_LIBRARY \ No newline at end of file diff --git a/AAD.h b/libraries/aad/includes/AAD.h similarity index 56% rename from AAD.h rename to libraries/aad/includes/AAD.h index d50a060..bd7a1c0 100644 --- a/AAD.h +++ b/libraries/aad/includes/AAD.h @@ -4,7 +4,7 @@ Written by Antoine Savine in 2018 This code is the strict IP of Antoine Savine -License to use and alter this code for personal and commercial applications +License to use and alter this code for personal and commercial applications is freely granted to any person or company who purchased a copy of the book Modern Computational Finance: AAD and Parallel Simulations @@ -18,12 +18,9 @@ As long as this comment is preserved at the top of the file #include -// So we can instrument Gaussians like standard math functions -#include "gaussians.h" - // Use traditional AAD of chapter 10 (false) // or expression templated (AADET) of chapter 15 (true) -#define AADET true +#define AADET true #if AADET @@ -39,37 +36,32 @@ As long as this comment is preserved at the top of the file // Set static context for multi-dimensional AAD // RAII: reset dimension 1 on destruction -struct numResultsResetterForAAD -{ - ~numResultsResetterForAAD() - { - Tape::multi = false; - Node::numAdj = 1; - } +struct numResultsResetterForAAD { + ~numResultsResetterForAAD() { + Tape::multi = false; + Node::numAdj = 1; + } }; // Routine: set dimension and get RAII resetter -inline auto setNumResultsForAAD(const bool multi = false, const size_t numResults = 1) -{ - Tape::multi = multi; - Node::numAdj = numResults; - return make_unique(); +inline auto setNumResultsForAAD(const bool multi = false, + const size_t numResults = 1) { + Tape::multi = multi; + Node::numAdj = numResults; + return make_unique(); } // Other utilities // Put collection on tape -template -inline void putOnTape(IT begin, IT end) -{ - for_each(begin, end, [](Number& n) {n.putOnTape(); }); +template inline void putOnTape(IT begin, IT end) { + for_each(begin, end, [](Number &n) { n.putOnTape(); }); } // Convert collection between double and Number or reverse -template -inline void convertCollection(It1 srcBegin, It1 srcEnd, It2 destBegin) -{ - using destType = remove_reference_t; - transform(srcBegin, srcEnd, destBegin, - [](const auto& source) { return destType(source); }); +template +inline void convertCollection(It1 srcBegin, It1 srcEnd, It2 destBegin) { + using destType = remove_reference_t; + transform(srcBegin, srcEnd, destBegin, + [](const auto &source) { return destType(source); }); } diff --git a/libraries/aad/includes/AADExpr.h b/libraries/aad/includes/AADExpr.h new file mode 100644 index 0000000..7f9e5cd --- /dev/null +++ b/libraries/aad/includes/AADExpr.h @@ -0,0 +1,865 @@ + +/* +Written by Antoine Savine in 2018 + +This code is the strict IP of Antoine Savine + +License to use and alter this code for personal and commercial applications +is freely granted to any person or company who purchased a copy of the book + +Modern Computational Finance: AAD and Parallel Simulations +Antoine Savine +Wiley, 2018 + +As long as this comment is preserved at the top of the file +*/ + +#pragma once + +// Implementation of AAD with expression templates +// (AADET, chapter 15) + +// Defines expressions and the Number type + +// So we can instrument Gaussians like standard math functions +#include "AADTape.h" +#include "gaussians.h" +#include + +// Base CRTP expression class +// Note: overloaded operators catch all expressions and nothing else +template struct Expression { + // CRTP "virtualization" + double value() const { return static_cast(this)->value(); } + + // Just another interface + explicit operator double() const { return value(); } +}; + +// Note that Number is a leaf expression +// Defined in the bottom of the file + +// Binary expression +// LHS : the expression on the left +// RHS : the expression on the right +// OP : the binary operator +template +class BinaryExpression + // CRTP + : public Expression> { + const double myValue; + + const LHS lhs; + const RHS rhs; + +public: + // Constructor out of 2 expressions + // Note: eager evaluation on construction + explicit BinaryExpression(const Expression &l, const Expression &r) + : myValue(OP::eval(l.value(), r.value())), + lhs(static_cast(l)), rhs(static_cast(r)) {} + + // Value accessors + double value() const { return myValue; } + + // Expression template magic + // Expressions know + // AT COMPILE TIME + // the number of active inputs in their sub-expressions + enum { numNumbers = LHS::numNumbers + RHS::numNumbers }; + + // Push adjoint down the expression + // N : total number of active inputs in the expression + // n : number of active inputs already processed + template + void pushAdjoint( + // Node for the complete expression being processed + Node &exprNode, + // Adjoint cumulated for this binary node, or 1 if top + const double adjoint) const { + // Push on the left + if (LHS::numNumbers > 0) { + lhs.template pushAdjoint( + exprNode, + adjoint * OP::leftDerivative(lhs.value(), rhs.value(), value())); + } + + // Push on the right + if (RHS::numNumbers > 0) { + // Note left push processed LHS::numNumbers numbers + // So the next number to be processed is n + LHS::numNumbers + rhs.template pushAdjoint( + exprNode, + adjoint * OP::rightDerivative(lhs.value(), rhs.value(), value())); + } + } +}; + +// "Concrete" binaries, we only need to define operations and derivatives +struct OPMult { + static const double eval(const double l, const double r) { return l * r; } + + static const double leftDerivative(const double l, const double r, + const double v) { + return r; + } + + static const double rightDerivative(const double l, const double r, + const double v) { + return l; + } +}; + +struct OPAdd { + static const double eval(const double l, const double r) { return l + r; } + + static const double leftDerivative(const double l, const double r, + const double v) { + return 1.0; + } + + static const double rightDerivative(const double l, const double r, + const double v) { + return 1.0; + } +}; + +struct OPSub { + static const double eval(const double l, const double r) { return l - r; } + + static const double leftDerivative(const double l, const double r, + const double v) { + return 1.0; + } + + static const double rightDerivative(const double l, const double r, + const double v) { + return -1.0; + } +}; + +struct OPDiv { + static const double eval(const double l, const double r) { return l / r; } + + static const double leftDerivative(const double l, const double r, + const double v) { + return 1.0 / r; + } + + static const double rightDerivative(const double l, const double r, + const double v) { + return -l / r / r; + } +}; + +struct OPPow { + static const double eval(const double l, const double r) { return pow(l, r); } + + static const double leftDerivative(const double l, const double r, + const double v) { + return r * v / l; + } + + static const double rightDerivative(const double l, const double r, + const double v) { + return log(l) * v; + } +}; + +struct OPMax { + static const double eval(const double l, const double r) { return max(l, r); } + + static const double leftDerivative(const double l, const double r, + const double v) { + return l > r ? 1.0 : 0.0; + } + + static const double rightDerivative(const double l, const double r, + const double v) { + return r > l ? 1.0 : 0.0; + } +}; + +struct OPMin { + static const double eval(const double l, const double r) { return min(l, r); } + + static const double leftDerivative(const double l, const double r, + const double v) { + return l < r ? 1.0 : 0.0; + } + + static const double rightDerivative(const double l, const double r, + const double v) { + return r < l ? 1.0 : 0.0; + } +}; + +// Operator overloading for binary expressions +// build the corresponding expressions + +template +BinaryExpression operator*(const Expression &lhs, + const Expression &rhs) { + return BinaryExpression(lhs, rhs); +} + +template +BinaryExpression operator+(const Expression &lhs, + const Expression &rhs) { + return BinaryExpression(lhs, rhs); +} + +template +BinaryExpression operator-(const Expression &lhs, + const Expression &rhs) { + return BinaryExpression(lhs, rhs); +} + +template +BinaryExpression operator/(const Expression &lhs, + const Expression &rhs) { + return BinaryExpression(lhs, rhs); +} + +template +BinaryExpression pow(const Expression &lhs, + const Expression &rhs) { + return BinaryExpression(lhs, rhs); +} + +template +BinaryExpression max(const Expression &lhs, + const Expression &rhs) { + return BinaryExpression(lhs, rhs); +} + +template +BinaryExpression min(const Expression &lhs, + const Expression &rhs) { + return BinaryExpression(lhs, rhs); +} + +// Unary expressions : Same logic with one argument + +// The CRTP class +template +class UnaryExpression + // CRTP + : public Expression> { + const double myValue; + + const ARG arg; + + // For binary operators with a double on one side + // we store the double + const double dArg = 0.0; + +public: + // Constructor + // Note : eager evaluation on construction + explicit UnaryExpression(const Expression &a) + : myValue(OP::eval(a.value(), 0.0)), arg(static_cast(a)) {} + + // Special constructor for binary expressions with a double on one side + explicit UnaryExpression(const Expression &a, const double b) + : myValue(OP::eval(a.value(), b)), arg(static_cast(a)), + dArg(b) {} + + // Value accessors + double value() const { return myValue; } + + // Expression template magic + enum { numNumbers = ARG::numNumbers }; + + // Push adjoint down the expression + template + void pushAdjoint( + // Node for the complete expression being processed + Node &exprNode, + // Adjoint cumulated on the node, 1 if top + const double adjoint) const { + // Push into argument + if (ARG::numNumbers > 0) { + arg.template pushAdjoint( + exprNode, adjoint * OP::derivative(arg.value(), value(), dArg)); + } + } +}; + +// The unary operators + +struct OPExp { + static const double eval(const double r, const double d) { return exp(r); } + + static const double derivative(const double r, const double v, + const double d) { + return v; + } +}; + +struct OPLog { + static const double eval(const double r, const double d) { return log(r); } + + static const double derivative(const double r, const double v, + const double d) { + return 1.0 / r; + } +}; + +struct OPSqrt { + static const double eval(const double r, const double d) { return sqrt(r); } + + static const double derivative(const double r, const double v, + const double d) { + return 0.5 / v; + } +}; + +struct OPFabs { + static const double eval(const double r, const double d) { return fabs(r); } + + static const double derivative(const double r, const double v, + const double d) { + return r > 0.0 ? 1.0 : -1.0; + } +}; + +struct OPNormalDens { + static const double eval(const double r, const double d) { + return normalDens(r); + } + + static const double derivative(const double r, const double v, + const double d) { + return -r * v; + } +}; + +struct OPNormalCdf { + static const double eval(const double r, const double d) { + return normalCdf(r); + } + + static const double derivative(const double r, const double v, + const double d) { + return normalDens(r); + } +}; + +// Binary operators with a double on one side + +// * double or double * +struct OPMultD { + static const double eval(const double r, const double d) { return r * d; } + + static const double derivative(const double r, const double v, + const double d) { + return d; + } +}; + +// + double or double + +struct OPAddD { + static const double eval(const double r, const double d) { return r + d; } + + static const double derivative(const double r, const double v, + const double d) { + return 1.0; + } +}; + +// double - +struct OPSubDL { + static const double eval(const double r, const double d) { return d - r; } + + static const double derivative(const double r, const double v, + const double d) { + return -1.0; + } +}; + +// - double +struct OPSubDR { + static const double eval(const double r, const double d) { return r - d; } + + static const double derivative(const double r, const double v, + const double d) { + return 1.0; + } +}; + +// double / +struct OPDivDL { + static const double eval(const double r, const double d) { return d / r; } + + static const double derivative(const double r, const double v, + const double d) { + return -d / r / r; + } +}; + +// / double +struct OPDivDR { + static const double eval(const double r, const double d) { return r / d; } + + static const double derivative(const double r, const double v, + const double d) { + return 1.0 / d; + } +}; + +// pow (d,) +struct OPPowDL { + static const double eval(const double r, const double d) { return pow(d, r); } + + static const double derivative(const double r, const double v, + const double d) { + return log(d) * v; + } +}; + +// pow (,d) +struct OPPowDR { + static const double eval(const double r, const double d) { return pow(r, d); } + + static const double derivative(const double r, const double v, + const double d) { + return d * v / r; + } +}; + +// max (d,) +struct OPMaxD { + static const double eval(const double r, const double d) { return max(r, d); } + + static const double derivative(const double r, const double v, + const double d) { + return r > d ? 1.0 : 0.0; + } +}; + +// min (d,) +struct OPMinD { + static const double eval(const double r, const double d) { return min(r, d); } + + static const double derivative(const double r, const double v, + const double d) { + return r < d ? 1.0 : 0.0; + } +}; + +// And overloading + +template +UnaryExpression exp(const Expression &arg) { + return UnaryExpression(arg); +} + +template +UnaryExpression log(const Expression &arg) { + return UnaryExpression(arg); +} + +template +UnaryExpression sqrt(const Expression &arg) { + return UnaryExpression(arg); +} + +template +UnaryExpression fabs(const Expression &arg) { + return UnaryExpression(arg); +} + +template +UnaryExpression normalDens(const Expression &arg) { + return UnaryExpression(arg); +} + +template +UnaryExpression normalCdf(const Expression &arg) { + return UnaryExpression(arg); +} + +// Overloading continued, +// binary operators with a double on one side + +template +UnaryExpression operator*(const double d, + const Expression &rhs) { + return UnaryExpression(rhs, d); +} + +template +UnaryExpression operator*(const Expression &lhs, + const double d) { + return UnaryExpression(lhs, d); +} + +template +UnaryExpression operator+(const double d, + const Expression &rhs) { + return UnaryExpression(rhs, d); +} + +template +UnaryExpression operator+(const Expression &lhs, + const double d) { + return UnaryExpression(lhs, d); +} + +template +UnaryExpression operator-(const double d, + const Expression &rhs) { + return UnaryExpression(rhs, d); +} + +template +UnaryExpression operator-(const Expression &lhs, + const double d) { + return UnaryExpression(lhs, d); +} + +template +UnaryExpression operator/(const double d, + const Expression &rhs) { + return UnaryExpression(rhs, d); +} + +template +UnaryExpression operator/(const Expression &lhs, + const double d) { + return UnaryExpression(lhs, d); +} + +template +UnaryExpression pow(const double d, const Expression &rhs) { + return UnaryExpression(rhs, d); +} + +template +UnaryExpression pow(const Expression &lhs, const double d) { + return UnaryExpression(lhs, d); +} + +template +UnaryExpression max(const double d, const Expression &rhs) { + return UnaryExpression(rhs, d); +} + +template +UnaryExpression max(const Expression &lhs, const double d) { + return UnaryExpression(lhs, d); +} + +template +UnaryExpression min(const double d, const Expression &rhs) { + return UnaryExpression(rhs, d); +} + +template +UnaryExpression min(const Expression &lhs, const double d) { + return UnaryExpression(lhs, d); +} + +// Comparison, same as traditional + +template +bool operator==(const Expression &lhs, const Expression &rhs) { + return lhs.value() == rhs.value(); +} +template +bool operator==(const Expression &lhs, const double &rhs) { + return lhs.value() == rhs; +} +template +bool operator==(const double &lhs, const Expression &rhs) { + return lhs == rhs.value(); +} + +template +bool operator!=(const Expression &lhs, const Expression &rhs) { + return lhs.value() != rhs.value(); +} +template +bool operator!=(const Expression &lhs, const double &rhs) { + return lhs.value() != rhs; +} +template +bool operator!=(const double &lhs, const Expression &rhs) { + return lhs != rhs.value(); +} + +template +bool operator<(const Expression &lhs, const Expression &rhs) { + return lhs.value() < rhs.value(); +} +template bool operator<(const Expression &lhs, const double &rhs) { + return lhs.value() < rhs; +} +template bool operator<(const double &lhs, const Expression &rhs) { + return lhs < rhs.value(); +} + +template +bool operator>(const Expression &lhs, const Expression &rhs) { + return lhs.value() > rhs.value(); +} +template bool operator>(const Expression &lhs, const double &rhs) { + return lhs.value() > rhs; +} +template bool operator>(const double &lhs, const Expression &rhs) { + return lhs > rhs.value(); +} + +template +bool operator<=(const Expression &lhs, const Expression &rhs) { + return lhs.value() <= rhs.value(); +} +template +bool operator<=(const Expression &lhs, const double &rhs) { + return lhs.value() <= rhs; +} +template +bool operator<=(const double &lhs, const Expression &rhs) { + return lhs <= rhs.value(); +} + +template +bool operator>=(const Expression &lhs, const Expression &rhs) { + return lhs.value() >= rhs.value(); +} +template +bool operator>=(const Expression &lhs, const double &rhs) { + return lhs.value() >= rhs; +} +template +bool operator>=(const double &lhs, const Expression &rhs) { + return lhs >= rhs.value(); +} + +// Finally, unary +/- operators + +template +UnaryExpression operator-(const Expression &rhs) { + return 0.0 - rhs; +} + +template Expression operator+(const Expression &rhs) { + return rhs; +} + +// The Number type, also an expression + +class Number : public Expression { + // The value and node for this number, same as traditional + double myValue; + Node *myNode; + + // Node creation on tape + + template Node *createMultiNode() { return tape->recordNode(); } + + // Flattening: + // This is where, on assignment or construction from an expression, + // that derivatives are pushed through the expression's DAG + template + void fromExpr( + // RHS expression, will be flattened into this Number + const Expression &e) { + // Build expression node on tape + auto *node = createMultiNode(); + + // Push adjoints through expression with adjoint = 1 on top + static_cast(e).template pushAdjoint(*node, + 1.0); + + // Set my node + myNode = node; + } + +public: + // Expression template magic + enum { numNumbers = 1 }; + + // Push adjoint + // Numbers are expression leaves, + // pushAdjoint() receives their adjoint in the expression + // Numbers don't "push" anything, they register their derivatives on tape + template + void pushAdjoint( + // Node for the complete expression + Node &exprNode, + // Adjoint accumulated for this number, in the expression + const double adjoint) const { + // adjoint = d (expression) / d (thisNumber) : register on tape + // note n: index of this number on the node on tape + + // Register adjoint + exprNode.pAdjPtrs[n] = Tape::multi ? myNode->pAdjoints : &myNode->mAdjoint; + + // Register derivative + exprNode.pDerivatives[n] = adjoint; + } + + // Static access to tape, same as traditional + static thread_local Tape *tape; + + // Constructors + + Number() {} + + explicit Number(const double val) : myValue(val) { + // Create leaf + myNode = createMultiNode<0>(); + } + + Number &operator=(const double val) { + myValue = val; + // Create leaf + myNode = createMultiNode<0>(); + return *this; + } + + // No need for copy and assignment + // Default ones do the right thing: + // copy value and pointer to node on tape + + // Construct or assign from expression + + template Number(const Expression &e) : myValue(e.value()) { + // Flatten RHS expression + fromExpr(static_cast(e)); + } + + template Number &operator=(const Expression &e) { + myValue = e.value(); + // Flatten RHS expression + fromExpr(static_cast(e)); + return *this; + } + + // Explicit coversion to double + explicit operator double &() { return myValue; } + explicit operator double() const { return myValue; } + + // All the normal accessors and propagators, same as traditional + + // Put on tape + void putOnTape() { myNode = createMultiNode<0>(); } + + // Accessors: value and adjoint + + double &value() { return myValue; } + double value() const { return myValue; } + + // Single dimensional + double &adjoint() { return myNode->adjoint(); } + double adjoint() const { return myNode->adjoint(); } + + // Multi dimensional + double &adjoint(const size_t n) { return myNode->adjoint(n); } + double adjoint(const size_t n) const { return myNode->adjoint(n); } + + // Reset all adjoints on the tape + // note we don't use this method + void resetAdjoints() { tape->resetAdjoints(); } + + // Propagation + + // Propagate adjoints + // from and to both INCLUSIVE + static void propagateAdjoints(Tape::iterator propagateFrom, + Tape::iterator propagateTo) { + auto it = propagateFrom; + while (it != propagateTo) { + it->propagateOne(); + --it; + } + it->propagateOne(); + } + + // Convenient overloads + + // Set the adjoint on this node to 1, + // Then propagate from the node + void propagateAdjoints( + // We start on this number's node + Tape::iterator propagateTo) { + // Set this adjoint to 1 + adjoint() = 1.0; + // Find node on tape + auto it = tape->find(myNode); + // Reverse and propagate until we hit the stop + while (it != propagateTo) { + it->propagateOne(); + --it; + } + it->propagateOne(); + } + + // These 2 set the adjoint to 1 on this node + void propagateToStart() { propagateAdjoints(tape->begin()); } + void propagateToMark() { propagateAdjoints(tape->markIt()); } + + // This one only propagates + // Note: propagation starts at mark - 1 + static void propagateMarkToStart() { + propagateAdjoints(prev(tape->markIt()), tape->begin()); + } + + // Multi-adjoint propagation + + // Multi dimensional case: + // Propagate adjoints from and to both INCLUSIVE + static void propagateAdjointsMulti(Tape::iterator propagateFrom, + Tape::iterator propagateTo) { + auto it = propagateFrom; + while (it != propagateTo) { + it->propagateAll(); + --it; + } + it->propagateAll(); + } + + // Unary operators + + template Number &operator+=(const Expression &e) { + *this = *this + e; + return *this; + } + + template Number &operator*=(const Expression &e) { + *this = *this * e; + return *this; + } + + template Number &operator-=(const Expression &e) { + *this = *this - e; + return *this; + } + + template Number &operator/=(const Expression &e) { + *this = *this / e; + return *this; + } + + Number &operator+=(const double &e) { + *this = *this + e; + return *this; + } + + Number &operator*=(const double &e) { + *this = *this * e; + return *this; + } + + Number &operator-=(const double &e) { + *this = *this - e; + return *this; + } + + Number &operator/=(const double &e) { + *this = *this / e; + return *this; + } +}; diff --git a/libraries/aad/includes/AADNode.h b/libraries/aad/includes/AADNode.h new file mode 100644 index 0000000..1a4b04c --- /dev/null +++ b/libraries/aad/includes/AADNode.h @@ -0,0 +1,95 @@ + +/* +Written by Antoine Savine in 2018 + +This code is the strict IP of Antoine Savine + +License to use and alter this code for personal and commercial applications +is freely granted to any person or company who purchased a copy of the book + +Modern Computational Finance: AAD and Parallel Simulations +Antoine Savine +Wiley, 2018 + +As long as this comment is preserved at the top of the file +*/ + +#pragma once + +// AAD implementation of chapter 10 +// (With multi-dimensional additions of chapter 14) + +// Implementation of Node = record on tape + +// Unchanged for AADET of chapter 15 + +#include +#include // IWYU pragma: keep +using namespace std; + +class Node { + friend class Tape; + friend class Number; + friend auto setNumResultsForAAD(const bool, const size_t); + friend struct numResultsResetterForAAD; + + // The adjoint(s) + // in single case, self held (chapter 10) + double mAdjoint = 0; + // in multi case, held separately and accessed by pointer (chapter 14) + double *pAdjoints; + + // Data lives in separate memory + + // the n derivatives to arguments, + double *pDerivatives; + + // the n pointers to the adjoints of arguments + double **pAdjPtrs; + + // Number of adjoints (results) to propagate, usually 1 + // See chapter 14 + static size_t numAdj; + + // Number of childs (arguments) + const size_t n; + +public: + Node(const size_t N = 0) : n(N) {} + + // Access to adjoint(s) + // single + double &adjoint() { return mAdjoint; } + // multi + double &adjoint(const size_t n) { return pAdjoints[n]; } + + // Back-propagate adjoints to arguments adjoints + + // Single case, chapter 10 + void propagateOne() { + // Nothing to propagate + if (!n || !mAdjoint) + return; + + for (size_t i = 0; i < n; ++i) { + *(pAdjPtrs[i]) += pDerivatives[i] * mAdjoint; + } + } + + // Multi case, chapter 14 + void propagateAll() { + // No adjoint to propagate + if (!n || all_of(pAdjoints, pAdjoints + numAdj, + [](const double &x) { return !x; })) + return; + + for (size_t i = 0; i < n; ++i) { + double *adjPtrs = pAdjPtrs[i], ders = pDerivatives[i]; + + // Vectorized! + for (size_t j = 0; j < numAdj; ++j) { + adjPtrs[j] += ders * pAdjoints[j]; + } + } + } +}; \ No newline at end of file diff --git a/libraries/aad/includes/AADNumber.h b/libraries/aad/includes/AADNumber.h new file mode 100644 index 0000000..3614b23 --- /dev/null +++ b/libraries/aad/includes/AADNumber.h @@ -0,0 +1,547 @@ + +/* +Written by Antoine Savine in 2018 + +This code is the strict IP of Antoine Savine + +License to use and alter this code for personal and commercial applications +is freely granted to any person or company who purchased a copy of the book + +Modern Computational Finance: AAD and Parallel Simulations +Antoine Savine +Wiley, 2018 + +As long as this comment is preserved at the top of the file +*/ + +#pragma once + +// Traditional AAD implementation of chapter 10 +// (With multi-dimensional additions of chapter 14) + +// The custom number type + +#include "AADTape.h" +#include "gaussians.h" +#include // IWYU pragma: keep + +class Number { + // Value and node are the only data members + double myValue; + Node *myNode; + + // Create node on tape + template void createNode() { myNode = tape->recordNode(); } + + // Access node (friends only) + // Note const incorectness + Node &node() const { + +#ifdef _DEBUG + + // Help identify errors when arguments are not on tape + + // Find node on tape + auto it = tape->find(myNode); + + // Not found + if (it == tape->end()) { + throw runtime_error("Put a breakpoint here"); + } + +#endif + // Const incorrectness + return const_cast(*myNode); + } + + // Convenient access to node data for friends + + double &derivative() { return myNode->pDerivatives[0]; } + double &lDer() { return myNode->pDerivatives[0]; } + double &rDer() { return myNode->pDerivatives[1]; } + + double *&adjPtr() { return myNode->pAdjPtrs[0]; } + double *&leftAdj() { return myNode->pAdjPtrs[0]; } + double *&rightAdj() { return myNode->pAdjPtrs[1]; } + + // Private constructors for operator overloading + + // Unary + Number(Node &arg, const double val) : myValue(val) { + createNode<1>(); + + myNode->pAdjPtrs[0] = Tape::multi ? arg.pAdjoints : &arg.mAdjoint; + } + + // Binary + Number(Node &lhs, Node &rhs, const double val) : myValue(val) { + createNode<2>(); + + if (Tape::multi) { + myNode->pAdjPtrs[0] = lhs.pAdjoints; + myNode->pAdjPtrs[1] = rhs.pAdjoints; + } else { + myNode->pAdjPtrs[0] = &lhs.mAdjoint; + myNode->pAdjPtrs[1] = &rhs.mAdjoint; + } + } + +public: + // Static access to tape + static thread_local Tape *tape; + + // Public constructors for leaves + + Number() {} + + // Put on tape on construction + explicit Number(const double val) : myValue(val) { createNode<0>(); } + + // Put on tape on assignment + Number &operator=(const double val) { + myValue = val; + createNode<0>(); + + return *this; + } + + // Explicitly put existing Number on tape + void putOnTape() { createNode<0>(); } + + // Explicit coversion to double + explicit operator double &() { return myValue; } + explicit operator double() const { return myValue; } + + // Accessors: value and adjoint + + double &value() { return myValue; } + double value() const { return myValue; } + // Single dimensional + double &adjoint() { return myNode->adjoint(); } + double adjoint() const { return myNode->adjoint(); } + // Multi dimensional + double &adjoint(const size_t n) { return myNode->adjoint(n); } + double adjoint(const size_t n) const { return myNode->adjoint(n); } + + // Reset all adjoints on the tape + // note we don't use this method + void resetAdjoints() { tape->resetAdjoints(); } + + // Propagation + + // Propagate adjoints + // from and to both INCLUSIVE + static void propagateAdjoints(Tape::iterator propagateFrom, + Tape::iterator propagateTo) { + auto it = propagateFrom; + while (it != propagateTo) { + it->propagateOne(); + --it; + } + it->propagateOne(); + } + + // Convenient overloads + + // Set the adjoint on this node to 1, + // Then propagate from the node + void propagateAdjoints( + // We start on this number's node + Tape::iterator propagateTo) { + // Set this adjoint to 1 + adjoint() = 1.0; + // Find node on tape + auto propagateFrom = tape->find(myNode); + propagateAdjoints(propagateFrom, propagateTo); + } + + // These 2 set the adjoint to 1 on this node + void propagateToStart() { propagateAdjoints(tape->begin()); } + void propagateToMark() { propagateAdjoints(tape->markIt()); } + + // This one only propagates + // Note: propagation starts at mark - 1 + static void propagateMarkToStart() { + propagateAdjoints(prev(tape->markIt()), tape->begin()); + } + + // Multi dimensional case: + // Propagate adjoints from and to both INCLUSIVE + static void propagateAdjointsMulti(Tape::iterator propagateFrom, + Tape::iterator propagateTo) { + auto it = propagateFrom; + while (it != propagateTo) { + it->propagateAll(); + --it; + } + it->propagateAll(); + } + + // Operator overloading + + inline friend Number operator+(const Number &lhs, const Number &rhs) { + const double e = lhs.value() + rhs.value(); + // Eagerly evaluate and put on tape + Number result(lhs.node(), rhs.node(), e); + // Eagerly compute derivatives + result.lDer() = 1.0; + result.rDer() = 1.0; + + return result; + } + inline friend Number operator+(const Number &lhs, const double &rhs) { + const double e = lhs.value() + rhs; + // Eagerly evaluate and put on tape + Number result(lhs.node(), e); + // Eagerly compute derivatives + result.derivative() = 1.0; + + return result; + } + inline friend Number operator+(const double &lhs, const Number &rhs) { + return rhs + lhs; + } + + inline friend Number operator-(const Number &lhs, const Number &rhs) { + const double e = lhs.value() - rhs.value(); + // Eagerly evaluate and put on tape + Number result(lhs.node(), rhs.node(), e); + // Eagerly compute derivatives + result.lDer() = 1.0; + result.rDer() = -1.0; + + return result; + } + inline friend Number operator-(const Number &lhs, const double &rhs) { + const double e = lhs.value() - rhs; + // Eagerly evaluate and put on tape + Number result(lhs.node(), e); + // Eagerly compute derivatives + result.derivative() = 1.0; + + return result; + } + inline friend Number operator-(const double &lhs, const Number &rhs) { + const double e = lhs - rhs.value(); + // Eagerly evaluate and put on tape + Number result(rhs.node(), e); + // Eagerly compute derivatives + result.derivative() = -1.0; + + return result; + } + + inline friend Number operator*(const Number &lhs, const Number &rhs) { + const double e = lhs.value() * rhs.value(); + // Eagerly evaluate and put on tape + Number result(lhs.node(), rhs.node(), e); + // Eagerly compute derivatives + result.lDer() = rhs.value(); + result.rDer() = lhs.value(); + + return result; + } + inline friend Number operator*(const Number &lhs, const double &rhs) { + const double e = lhs.value() * rhs; + // Eagerly evaluate and put on tape + Number result(lhs.node(), e); + // Eagerly compute derivatives + result.derivative() = rhs; + + return result; + } + inline friend Number operator*(const double &lhs, const Number &rhs) { + return rhs * lhs; + } + + inline friend Number operator/(const Number &lhs, const Number &rhs) { + const double e = lhs.value() / rhs.value(); + // Eagerly evaluate and put on tape + Number result(lhs.node(), rhs.node(), e); + // Eagerly compute derivatives + const double invRhs = 1.0 / rhs.value(); + result.lDer() = invRhs; + result.rDer() = -lhs.value() * invRhs * invRhs; + + return result; + } + inline friend Number operator/(const Number &lhs, const double &rhs) { + const double e = lhs.value() / rhs; + // Eagerly evaluate and put on tape + Number result(lhs.node(), e); + // Eagerly compute derivatives + result.derivative() = 1.0 / rhs; + + return result; + } + inline friend Number operator/(const double &lhs, const Number &rhs) { + const double e = lhs / rhs.value(); + // Eagerly evaluate and put on tape + Number result(rhs.node(), e); + // Eagerly compute derivatives + result.derivative() = -lhs / rhs.value() / rhs.value(); + + return result; + } + + inline friend Number pow(const Number &lhs, const Number &rhs) { + const double e = pow(lhs.value(), rhs.value()); + // Eagerly evaluate and put on tape + Number result(lhs.node(), rhs.node(), e); + // Eagerly compute derivatives + result.lDer() = rhs.value() * e / lhs.value(); + result.rDer() = log(lhs.value()) * e; + + return result; + } + inline friend Number pow(const Number &lhs, const double &rhs) { + const double e = pow(lhs.value(), rhs); + // Eagerly evaluate and put on tape + Number result(lhs.node(), e); + // Eagerly compute derivatives + result.derivative() = rhs * e / lhs.value(); + + return result; + } + inline friend Number pow(const double &lhs, const Number &rhs) { + const double e = pow(lhs, rhs.value()); + // Eagerly evaluate and put on tape + Number result(rhs.node(), e); + // Eagerly compute derivatives + result.derivative() = log(lhs) * e; + + return result; + } + + inline friend Number max(const Number &lhs, const Number &rhs) { + const bool lmax = lhs.value() > rhs.value(); + // Eagerly evaluate and put on tape + Number result(lhs.node(), rhs.node(), lmax ? lhs.value() : rhs.value()); + // Eagerly compute derivatives + if (lmax) { + result.lDer() = 1.0; + result.rDer() = 0.0; + } else { + result.lDer() = 0.0; + result.rDer() = 1.0; + } + + return result; + } + inline friend Number max(const Number &lhs, const double &rhs) { + const bool lmax = lhs.value() > rhs; + // Eagerly evaluate and put on tape + Number result(lhs.node(), lmax ? lhs.value() : rhs); + // Eagerly compute derivatives + result.derivative() = lmax ? 1.0 : 0.0; + + return result; + } + inline friend Number max(const double &lhs, const Number &rhs) { + const bool rmax = rhs.value() > lhs; + // Eagerly evaluate and put on tape + Number result(rhs.node(), rmax ? rhs.value() : lhs); + // Eagerly compute derivatives + result.derivative() = rmax ? 1.0 : 0.0; + + return result; + } + + inline friend Number min(const Number &lhs, const Number &rhs) { + const bool lmin = lhs.value() < rhs.value(); + // Eagerly evaluate and put on tape + Number result(lhs.node(), rhs.node(), lmin ? lhs.value() : rhs.value()); + // Eagerly compute derivatives + if (lmin) { + result.lDer() = 1.0; + result.rDer() = 0.0; + } else { + result.lDer() = 0.0; + result.rDer() = 1.0; + } + + return result; + } + inline friend Number min(const Number &lhs, const double &rhs) { + const bool lmin = lhs.value() < rhs; + // Eagerly evaluate and put on tape + Number result(lhs.node(), lmin ? lhs.value() : rhs); + // Eagerly compute derivatives + result.derivative() = lmin ? 1.0 : 0.0; + + return result; + } + inline friend Number min(const double &lhs, const Number &rhs) { + const bool rmin = rhs.value() < lhs; + // Eagerly evaluate and put on tape + Number result(rhs.node(), rmin ? rhs.value() : lhs); + // Eagerly compute derivatives + result.derivative() = rmin ? 1.0 : 0.0; + + return result; + } + + Number &operator+=(const Number &arg) { + *this = *this + arg; + return *this; + } + Number &operator+=(const double &arg) { + *this = *this + arg; + return *this; + } + + Number &operator-=(const Number &arg) { + *this = *this - arg; + return *this; + } + Number &operator-=(const double &arg) { + *this = *this - arg; + return *this; + } + + Number &operator*=(const Number &arg) { + *this = *this * arg; + return *this; + } + Number &operator*=(const double &arg) { + *this = *this * arg; + return *this; + } + + Number &operator/=(const Number &arg) { + *this = *this / arg; + return *this; + } + Number &operator/=(const double &arg) { + *this = *this / arg; + return *this; + } + + // Unary +/- + Number operator-() const { return 0.0 - *this; } + Number operator+() const { return *this; } + + // Overloading continued, unary functions + + inline friend Number exp(const Number &arg) { + const double e = exp(arg.value()); + // Eagerly evaluate and put on tape + Number result(arg.node(), e); + // Eagerly compute derivatives + result.derivative() = e; + + return result; + } + + inline friend Number log(const Number &arg) { + const double e = log(arg.value()); + // Eagerly evaluate and put on tape + Number result(arg.node(), e); + // Eagerly compute derivatives + result.derivative() = 1.0 / arg.value(); + + return result; + } + + inline friend Number sqrt(const Number &arg) { + const double e = sqrt(arg.value()); + // Eagerly evaluate and put on tape + Number result(arg.node(), e); + // Eagerly compute derivatives + result.derivative() = 0.5 / e; + + return result; + } + + inline friend Number fabs(const Number &arg) { + const double e = fabs(arg.value()); + // Eagerly evaluate and put on tape + Number result(arg.node(), e); + // Eagerly compute derivatives + result.derivative() = arg.value() > 0.0 ? 1.0 : -1.0; + + return result; + } + + inline friend Number normalDens(const Number &arg) { + const double e = normalDens(arg.value()); + // Eagerly evaluate and put on tape + Number result(arg.node(), e); + // Eagerly compute derivatives + result.derivative() = -arg.value() * e; + + return result; + } + + inline friend Number normalCdf(const Number &arg) { + const double e = normalCdf(arg.value()); + // Eagerly evaluate and put on tape + Number result(arg.node(), e); + // Eagerly compute derivatives + result.derivative() = normalDens(arg.value()); + + return result; + } + + // Finally, comparison + + inline friend bool operator==(const Number &lhs, const Number &rhs) { + return lhs.value() == rhs.value(); + } + inline friend bool operator==(const Number &lhs, const double &rhs) { + return lhs.value() == rhs; + } + inline friend bool operator==(const double &lhs, const Number &rhs) { + return lhs == rhs.value(); + } + + inline friend bool operator!=(const Number &lhs, const Number &rhs) { + return lhs.value() != rhs.value(); + } + inline friend bool operator!=(const Number &lhs, const double &rhs) { + return lhs.value() != rhs; + } + inline friend bool operator!=(const double &lhs, const Number &rhs) { + return lhs != rhs.value(); + } + + inline friend bool operator<(const Number &lhs, const Number &rhs) { + return lhs.value() < rhs.value(); + } + inline friend bool operator<(const Number &lhs, const double &rhs) { + return lhs.value() < rhs; + } + inline friend bool operator<(const double &lhs, const Number &rhs) { + return lhs < rhs.value(); + } + + inline friend bool operator>(const Number &lhs, const Number &rhs) { + return lhs.value() > rhs.value(); + } + inline friend bool operator>(const Number &lhs, const double &rhs) { + return lhs.value() > rhs; + } + inline friend bool operator>(const double &lhs, const Number &rhs) { + return lhs > rhs.value(); + } + + inline friend bool operator<=(const Number &lhs, const Number &rhs) { + return lhs.value() <= rhs.value(); + } + inline friend bool operator<=(const Number &lhs, const double &rhs) { + return lhs.value() <= rhs; + } + inline friend bool operator<=(const double &lhs, const Number &rhs) { + return lhs <= rhs.value(); + } + + inline friend bool operator>=(const Number &lhs, const Number &rhs) { + return lhs.value() >= rhs.value(); + } + inline friend bool operator>=(const Number &lhs, const double &rhs) { + return lhs.value() >= rhs; + } + inline friend bool operator>=(const double &lhs, const Number &rhs) { + return lhs >= rhs.value(); + } +}; diff --git a/libraries/aad/includes/AADTape.h b/libraries/aad/includes/AADTape.h new file mode 100644 index 0000000..fe70b41 --- /dev/null +++ b/libraries/aad/includes/AADTape.h @@ -0,0 +1,150 @@ + +/* +Written by Antoine Savine in 2018 + +This code is the strict IP of Antoine Savine + +License to use and alter this code for personal and commercial applications +is freely granted to any person or company who purchased a copy of the book + +Modern Computational Finance: AAD and Parallel Simulations +Antoine Savine +Wiley, 2018 + +As long as this comment is preserved at the top of the file +*/ + +#pragma once + +// AAD implementation of chapter 10 +// (With multi-dimensional additions of chapter 14) + +// Implementation of the Tape + +// Unchanged for AADET of chapter 15 + +#include "AADNode.h" +#include "blocklist.h" + +constexpr size_t BLOCKSIZE = 16384; // Number of nodes +constexpr size_t ADJSIZE = 32768; // Number of adjoints +constexpr size_t DATASIZE = 65536; // Data in bytes + +class Tape { + // Working with multiple results / adjoints? + static bool multi; + + // Storage for adjoints in multi-dimensional case (chapter 14) + blocklist myAdjointsMulti; + + // Storage for derivatives and child adjoint pointers + blocklist myDers; + blocklist myArgPtrs; + + // Storage for the nodes + blocklist myNodes; + + // Padding so tapes in a vector don't interfere + char myPad[64]; + + friend auto setNumResultsForAAD(const bool, const size_t); + friend struct numResultsResetterForAAD; + friend class Number; + +public: + // Build note in place and return a pointer + // N : number of childs (arguments) + template Node *recordNode() { + // Construct the node in place on tape + Node *node = myNodes.emplace_back(N); + + // Store and zero the adjoint(s) + if (multi) { + node->pAdjoints = myAdjointsMulti.emplace_back_multi(Node::numAdj); + fill(node->pAdjoints, node->pAdjoints + Node::numAdj, 0.0); + } + + // Store the derivatives and child adjoint pointers unless leaf + if constexpr (N > 0) { + node->pDerivatives = myDers.emplace_back_multi(); + node->pAdjPtrs = myArgPtrs.emplace_back_multi(); + } + + return node; + } + + // Reset all adjoints to 0 + void resetAdjoints() { + if (multi) { + //default is 0 + myAdjointsMulti.myMemset(); + } else { + for (Node &node : myNodes) { + node.mAdjoint = 0; + } + } + } + + // Clear + void clear() { + myAdjointsMulti.clear(); + myDers.clear(); + myArgPtrs.clear(); + myNodes.clear(); + } + + // Rewind + void rewind() { + +#ifdef _DEBUG + + // In debug mode, always wipe + // makes it easier to identify errors + + clear(); + +#else + // In release mode, rewind and reuse + + if (multi) { + myAdjointsMulti.rewind(); + } + myDers.rewind(); + myArgPtrs.rewind(); + myNodes.rewind(); + +#endif + } + + // Set mark + void mark() { + if (multi) { + myAdjointsMulti.setmark(); + } + myDers.setmark(); + myArgPtrs.setmark(); + myNodes.setmark(); + } + + // Rewind to mark + void rewindToMark() { + if (multi) { + myAdjointsMulti.rewind_to_mark(); + } + myDers.rewind_to_mark(); + myArgPtrs.rewind_to_mark(); + myNodes.rewind_to_mark(); + } + + // Iterators + + using iterator = blocklist::iterator; + + auto begin() { return myNodes.begin(); } + + auto end() { return myNodes.end(); } + + auto markIt() { return myNodes.mark(); } + + auto find(Node *node) { return myNodes.find(node); } +}; \ No newline at end of file diff --git a/libraries/aad/includes/ConcurrentQueue.h b/libraries/aad/includes/ConcurrentQueue.h new file mode 100644 index 0000000..414ad3e --- /dev/null +++ b/libraries/aad/includes/ConcurrentQueue.h @@ -0,0 +1,93 @@ +#pragma once + +// Concurrent queue of chapter 3, +// Used in the thread pool + +#include +#include +#include +using namespace std; + +template class ConcurrentQueue { + + queue myQueue; + mutable mutex myMutex; + condition_variable myCV; + bool myInterrupt; + +public: + ConcurrentQueue() : myInterrupt(false) {} + ~ConcurrentQueue() { interrupt(); } + + bool empty() const { + // Lock + lock_guard lk(myMutex); + // Access underlying queue + return myQueue.empty(); + } // Unlock + + // Pop into argument + bool tryPop(T &t) { + // Lock + lock_guard lk(myMutex); + if (myQueue.empty()) + return false; + // Move from queue + t = std::move(myQueue.front()); + // Combine front/pop + myQueue.pop(); + + return true; + } // Unlock + + // Pass t byVal or move with push( move( t)) + void push(T t) { + { + // Lock + lock_guard lk(myMutex); + // Move into queue + myQueue.push(std::move(t)); + } // Unlock before notification + + // Unlock before notification + myCV.notify_one(); + } + + // Wait if empty + bool pop(T &t) { + // (Unique) lock + unique_lock lk(myMutex); + + // Wait if empty, release lock until notified + while (!myInterrupt && myQueue.empty()) + myCV.wait(lk); + + // Re-acquire lock, resume + + // Check for interruption + if (myInterrupt) + return false; + + // Combine front/pop + t = std::move(myQueue.front()); + myQueue.pop(); + + return true; + + } // Unlock + + void interrupt() { + { + lock_guard lk(myMutex); + myInterrupt = true; + } + myCV.notify_all(); + } + + void resetInterrupt() { myInterrupt = false; } + + void clear() { + queue empty; + swap(myQueue, empty); + } +}; \ No newline at end of file diff --git a/libraries/aad/includes/analytics.h b/libraries/aad/includes/analytics.h new file mode 100644 index 0000000..f832d00 --- /dev/null +++ b/libraries/aad/includes/analytics.h @@ -0,0 +1,100 @@ +#pragma once + +#include "gaussians.h" + +// Black and Scholes formula, templated + +// Price +template +inline T blackScholes(const U spot, const V strike, const T vol, const W mat) { + const auto std = vol * sqrt(mat); + if (std <= EPS) + return max(T(0.0), T(spot - strike)); + const auto d2 = log(spot / strike) / std - 0.5 * std; + const auto d1 = d2 + std; + return spot * normalCdf(d1) - strike * normalCdf(d2); +} + +// Implied vol, untemplated +inline double blackScholesIvol(const double spot, const double strike, + const double prem, const double mat) { + if (prem <= max(0.0, spot - strike) + EPS) + return 0.0; + + double p, pu, pl; + double u = 0.5; + while (blackScholes(spot, strike, u, mat) < prem) + u *= 2; + double l = 0.05; + while (blackScholes(spot, strike, l, mat) > prem) + l /= 2; + pu = blackScholes(spot, strike, u, mat); + blackScholes(spot, strike, l, mat); + + while (u - l > 1.e-12) { + const double m = 0.5 * (u + l); + p = blackScholes(spot, strike, m, mat); + if (p > prem) { + u = m; + pu = p; + } else { + l = m; + pl = p; + } + } + + return l + (prem - pl) / (pu - pl) * (u - l); +} + +// Merton, templated + +template +inline T merton(const U spot, const V strike, const T vol, const W mat, + const X intens, const X meanJmp, const X stdJmp) { + const auto varJmp = stdJmp * stdJmp; + const auto mv2 = meanJmp + 0.5 * varJmp; + const auto comp = intens * (exp(mv2) - 1); + const auto var = vol * vol; + const auto intensT = intens * mat; + + unsigned fact = 1; + X iT = 1.0; + const size_t cut = 10; + T result = 0.0; + for (size_t n = 0; n < cut; ++n) { + const auto s = spot * exp(n * mv2 - comp * mat); + const auto v = sqrt(var + n * varJmp / mat); + const auto prob = exp(-intensT) * iT / fact; + result += prob * blackScholes(s, strike, v, mat); + fact *= n + 1; + iT *= intensT; + } + + return result; +} + +// Up and out call in Black-Scholes, untemplated + +inline double BlackScholesKO(const double spot, const double rate, + const double div, const double strike, + const double barrier, const double mat, + const double vol) { + const double std = vol * sqrt(mat); + const double fwdFact = exp((rate - div) * mat); + const double fwd = spot * fwdFact; + const double disc = exp(-rate * mat); + const double v = rate - div - 0.5 * vol * vol; + const double d2 = (log(spot / barrier) + v * mat) / std; + const double d2prime = (log(barrier / spot) + v * mat) / std; + + const double bar = + blackScholes(fwd, strike, vol, mat) - + blackScholes(fwd, barrier, vol, mat) - + (barrier - strike) * normalCdf(d2) - + pow(barrier / spot, 2 * v / vol / vol) * + (blackScholes(fwdFact * barrier * barrier / spot, strike, vol, mat) - + blackScholes(fwdFact * barrier * barrier / spot, barrier, vol, mat) - + (barrier - strike) * normalCdf(d2prime)); + + return disc * bar; +} diff --git a/libraries/aad/includes/blocklist.h b/libraries/aad/includes/blocklist.h new file mode 100644 index 0000000..109f8fc --- /dev/null +++ b/libraries/aad/includes/blocklist.h @@ -0,0 +1,291 @@ + +/* +Written by Antoine Savine in 2018 + +This code is the strict IP of Antoine Savine + +License to use and alter this code for personal and commercial applications +is freely granted to any person or company who purchased a copy of the book + +Modern Computational Finance: AAD and Parallel Simulations +Antoine Savine +Wiley, 2018 + +As long as this comment is preserved at the top of the file +*/ + +#pragma once +#pragma warning(disable : 4018) + +// Blocklist data structure for AAD memory management +// See chapter 10, unchanged with expression templates of chapter 15 + +#include +#include +#include +#include +using namespace std; + +template class blocklist { + // Container = list of blocks + list> data; + + using list_iter = decltype(data.begin()); + using block_iter = decltype(data.back().begin()); + + // Current block + list_iter cur_block; + + // Last block + list_iter last_block; + + // Next free space in current block + block_iter next_space; + + // Last free space (+1) in current block + block_iter last_space; + + // Mark + list_iter marked_block; + block_iter marked_space; + + // Create new array + void newblock() { + data.emplace_back(); + cur_block = last_block = prev(data.end()); + next_space = cur_block->begin(); + last_space = cur_block->end(); + } + + // Move on to next array + void nextblock() { + // This is the last array: create new + if (cur_block == last_block) { + newblock(); + } else { + ++cur_block; + next_space = cur_block->begin(); + last_space = cur_block->end(); + } + } + +public: + // Create first block on construction + blocklist() { newblock(); } + + // Factory reset + void clear() { + data.clear(); + newblock(); + } + + // Rewind but keep all blocks + void rewind() { + cur_block = data.begin(); + next_space = cur_block->begin(); + last_space = cur_block->end(); + } + + // Memset + void myMemset(unsigned char value = 0) { + for (auto &arr : data) { + + memset(&arr[0], value, block_size * sizeof(T)); + + } + } + + // Construct object of type T in place + // in the next free space and return a pointer on it + // Implements perfect forwarding of constructor arguments + template T *emplace_back(Args &&...args) { + // No more space in current array + if (next_space == last_space) { + nextblock(); + } + // Placement new, construct in memory pointed by next + T *emplaced = new (&*next_space) // memory pointed by next as T* + T(std::forward(args)...); // perfect forwarding of ctor arguments + + // Advance next + ++next_space; + + // Return + return emplaced; + } + + // Overload for default constructed + T *emplace_back() { + // No more space in current array + if (next_space == last_space) { + nextblock(); + } + + // Current space + auto old_next = next_space; + + // Advance next + ++next_space; + + // Return + return &*old_next; + } + + // Stores n default constructed elements + // and returns a pointer on the first + + // Version 1: n known at compile time + template T *emplace_back_multi() { + // No more space in current array + if (distance(next_space, last_space) < n) { + nextblock(); + } + + // Current space + auto old_next = next_space; + + // Advance next + next_space += n; + + // Return + return &*old_next; + } + + // Version 2: n unknown at compile time + T *emplace_back_multi(const size_t n) { + // No more space in current array + if (distance(next_space, last_space) < n) { + nextblock(); + } + + // Current space + auto old_next = next_space; + + // Advance next + next_space += n; + + // Return + return &*old_next; + } + + // Marks + + // Set mark + void setmark() { + if (next_space == last_space) { + nextblock(); + } + + marked_block = cur_block; + marked_space = next_space; + } + + // Rewind to mark + void rewind_to_mark() { + cur_block = marked_block; + next_space = marked_space; + last_space = cur_block->end(); + } + + // Iterator + + class iterator { + // List and block + list_iter cur_block; // current block + block_iter cur_space; // current space + block_iter first_space; // first space in block + block_iter last_space; // last (+1) space in block + + public: + // iterator traits + using difference_type = ptrdiff_t; + using reference = T &; + using pointer = T *; + using value_type = T; + using iterator_category = bidirectional_iterator_tag; + + // Default constructor + iterator() {} + + // Constructor + iterator(list_iter cb, block_iter cs, block_iter fs, block_iter ls) + : cur_block(cb), cur_space(cs), first_space(fs), last_space(ls) {} + + // Pre-increment (we do not provide post) + iterator &operator++() { + ++cur_space; + if (cur_space == last_space) { + ++cur_block; + first_space = cur_block->begin(); + last_space = cur_block->end(); + cur_space = first_space; + } + + return *this; + } + + // Pre-decrement + iterator &operator--() { + if (cur_space == first_space) { + --cur_block; + first_space = cur_block->begin(); + last_space = cur_block->end(); + cur_space = last_space; + } + + --cur_space; + + return *this; + } + + // Access to contained elements + T &operator*() { return *cur_space; } + const T &operator*() const { return *cur_space; } + T *operator->() { return &*cur_space; } + const T *operator->() const { return &*cur_space; } + + // Check equality + bool operator==(const iterator &rhs) const { + return (cur_block == rhs.cur_block && cur_space == rhs.cur_space); + } + bool operator!=(const iterator &rhs) const { + return (cur_block != rhs.cur_block || cur_space != rhs.cur_space); + } + }; + + // Access to iterators + + iterator begin() { + return iterator(data.begin(), data.begin()->begin(), data.begin()->begin(), + data.begin()->end()); + } + + iterator end() { + return iterator(cur_block, next_space, cur_block->begin(), + cur_block->end()); + } + + // Iterator on mark + iterator mark() { + return iterator(marked_block, marked_space, marked_block->begin(), + marked_block->end()); + } + + // Find element, by pointer, searching sequentially from the end + iterator find(const T *const element) { + // Search from the end + iterator it = end(); + iterator b = begin(); + + while (it != b) { + --it; + if (&*it == element) + return it; + } + + if (&*it == element) + return it; + + return end(); + } +}; \ No newline at end of file diff --git a/libraries/aad/includes/choldc.h b/libraries/aad/includes/choldc.h new file mode 100644 index 0000000..131fae8 --- /dev/null +++ b/libraries/aad/includes/choldc.h @@ -0,0 +1,39 @@ +#pragma once + +// Choleski decomposition, inspired by Numerical Recipes + +#include "matrix.h" + +template void choldc(const matrix &in, matrix &out) { + int n = in.rows(); + T sum; + + fill(out.begin(), out.end(), T(0.0)); + + for (int i = 0; i < n; ++i) { + auto *ai = in[i]; + auto *pi = out[i]; + for (int j = 0; j <= i; ++j) { + auto *aj = in[j]; + auto *pj = out[j]; + sum = ai[j]; + for (int k = 0; k < j; ++k) { + sum -= pi[k] * pj[k]; + } + if (i == j) { + if (sum < -1.0e-15) { + throw runtime_error("choldc : matrix not positive definite"); + } + if (sum < 1.0e-15) + sum = 0.0; + pi[i] = sqrt(sum); + } else { + if (fabs(pj[j]) < 1.0e-15) { + pi[j] = 0.0; + } else { + pi[j] = sum / pj[j]; + } + } + } + } +} \ No newline at end of file diff --git a/libraries/aad/includes/gaussians.h b/libraries/aad/includes/gaussians.h new file mode 100644 index 0000000..855f51e --- /dev/null +++ b/libraries/aad/includes/gaussians.h @@ -0,0 +1,91 @@ +#pragma once + +#include // IWYU pragma: keep +#include +#include // IWYU pragma: keep +using namespace std; + +#define EPS 1.0e-08 + +// Gaussian functions + +// Normal density +inline double normalDens(const double x) { + return x < -10.0 || 10.0 < x ? 0.0 : exp(-0.5 * x * x) / 2.506628274631; +} + +// Normal CDF (N in Black-Scholes) +// Zelen and Severo's approximation (1964) +// See +// https://en.wikipedia.org/wiki/Normal_distribution#Numerical_approximations_for_the_normal_CDF +inline double normalCdf(const double x) { + if (x < -10.0) + return 0.0; + if (x > 10.0) + return 1.0; + if (x < 0.0) + return 1.0 - normalCdf(-x); + + static constexpr double p = 0.2316419; + static constexpr double b1 = 0.319381530; + static constexpr double b2 = -0.356563782; + static constexpr double b3 = 1.781477937; + static constexpr double b4 = -1.821255978; + static constexpr double b5 = 1.330274429; + + const auto t = 1.0 / (1.0 + p * x); + + const auto pol = t * (b1 + t * (b2 + t * (b3 + t * (b4 + t * b5)))); + + const auto pdf = normalDens(x); + + return 1.0 - pdf * pol; +} + +// Inverse CDF (for generation of Gaussians out of Uniforms) +// Beasley-Springer-Moro algorithm +// Moro, The full Monte, Risk, 1995 +// See Glasserman, Monte Carlo Methods in Financial Engineering, p 68 +inline double invNormalCdf(const double p) { + const bool sup = p > 0.5; + const double up = sup ? 1.0 - p : p; + + static constexpr double a0 = 2.50662823884; + static constexpr double a1 = -18.61500062529; + static constexpr double a2 = 41.39119773534; + static constexpr double a3 = -25.44106049637; + + static constexpr double b0 = -8.47351093090; + static constexpr double b1 = 23.08336743743; + static constexpr double b2 = -21.06224101826; + static constexpr double b3 = 3.13082909833; + + static constexpr double c0 = 0.3374754822726147; + static constexpr double c1 = 0.9761690190917186; + static constexpr double c2 = 0.1607979714918209; + static constexpr double c3 = 0.0276438810333863; + static constexpr double c4 = 0.0038405729373609; + static constexpr double c5 = 0.0003951896511919; + static constexpr double c6 = 0.0000321767881768; + static constexpr double c7 = 0.0000002888167364; + static constexpr double c8 = 0.0000003960315187; + + double x = up - 0.5; + double r; + + if (fabs(x) < 0.42) { + r = x * x; + r = x * (((a3 * r + a2) * r + a1) * r + a0) / + ((((b3 * r + b2) * r + b1) * r + b0) * r + 1.0); + return sup ? -r : r; + } + + r = up; + r = log(-log(r)); + r = c0 + + r * (c1 + + r * (c2 + + r * (c3 + r * (c4 + r * (c5 + r * (c6 + r * (c7 + r * c8))))))); + + return sup ? r : -r; +} diff --git a/libraries/aad/includes/interp.h b/libraries/aad/includes/interp.h new file mode 100644 index 0000000..178f948 --- /dev/null +++ b/libraries/aad/includes/interp.h @@ -0,0 +1,102 @@ + +/* +Written by Antoine Savine in 2018 + +This code is the strict IP of Antoine Savine + +License to use and alter this code for personal and commercial applications +is freely granted to any person or company who purchased a copy of the book + +Modern Computational Finance: AAD and Parallel Simulations +Antoine Savine +Wiley, 2018 + +As long as this comment is preserved at the top of the file +*/ + +#pragma once + +#include "matrix.h" +#include +#include +using namespace std; + +// Interpolation utility + +// Interpolates the vector y against knots x in value x0 +// Interpolation is linear or smooth, extrapolation is flat +template +inline auto interp( + // sorted on xs + ITX xBegin, ITX xEnd, + // corresponding ys + ITY yBegin, ITY yEnd, + // interpolate for point x0 + const T &x0) -> remove_reference_t { + // STL binary search, returns iterator on 1st no less than x0 + // upper_bound guarantees logarithmic search + auto it = upper_bound(xBegin, xEnd, x0); + + // Extrapolation? + if (it == xEnd) + return *(yEnd - 1); + if (it == xBegin) + return *yBegin; + + // Interpolation + size_t n = distance(xBegin, it) - 1; + auto x1 = xBegin[n]; + auto y1 = yBegin[n]; + auto x2 = xBegin[n + 1]; + auto y2 = yBegin[n + 1]; + + auto t = (x0 - x1) / (x2 - x1); + + if constexpr (smoothStep) { + return y1 + (y2 - y1) * t * t * (3.0 - 2 * t); + } else { + return y1 + (y2 - y1) * t; + } +} + +// 2D +template +inline V interp2D( + // sorted on xs + const vector &x, + // sorted on ys + const vector &y, + // zs in a matrix + const matrix &z, + // interpolate for point (x0,y0) + const W &x0, const X &y0) { + const size_t n = x.size(); + const size_t m = y.size(); + + // STL binary search, returns iterator on 1st no less than x0 + // upper_boung guarantees logarithmic search + auto it = upper_bound(x.begin(), x.end(), x0); + const size_t n2 = distance(x.begin(), it); + + // Extrapolation in x? + if (n2 == n) + return interp(y.begin(), y.end(), z[n2 - 1], z[n2 - 1] + m, y0); + if (n2 == 0) + return interp(y.begin(), y.end(), z[0], z[0] + m, y0); + + // Interpolation in x + const size_t n1 = n2 - 1; + auto x1 = x[n1]; + auto x2 = x[n2]; + auto z1 = interp(y.begin(), y.end(), z[n1], z[n1] + m, y0); + auto z2 = interp(y.begin(), y.end(), z[n2], z[n2] + m, y0); + + // Smooth step + auto t = (x0 - x1) / (x2 - x1); + if constexpr (smoothStep) { + return z1 + (z2 - z1) * t * t * (3.0 - 2 * t); + } else { + return z1 + (z2 - z1) * t; + ; + } +} diff --git a/ivs.h b/libraries/aad/includes/ivs.h similarity index 100% rename from ivs.h rename to libraries/aad/includes/ivs.h diff --git a/libraries/aad/includes/main.h b/libraries/aad/includes/main.h new file mode 100644 index 0000000..3e488a2 --- /dev/null +++ b/libraries/aad/includes/main.h @@ -0,0 +1,638 @@ + +/* +Written by Antoine Savine in 2018 + +This code is the strict IP of Antoine Savine + +License to use and alter this code for personal and commercial applications +is freely granted to any person or company who purchased a copy of the book + +Modern Computational Finance: AAD and Parallel Simulations +Antoine Savine +Wiley, 2018 + +As long as this comment is preserved at the top of the file +*/ + +// Entry points to the library + +#pragma once + +#include "mcBase.h" +#include "mcMdl.h" // IWYU pragma: keep +#include "mcPrd.h" // IWYU pragma: keep +#include "mcPrdMulti.h" // IWYU pragma: keep +#include "mrg32k3a.h" +#include "sobol.h" +#include // IWYU pragma: keep +#include +using namespace std; + +#include "store.h" + +struct NumericalParam { + bool parallel; + bool useSobol; + int numPath; + int seed1 = 12345; + int seed2 = 1234; +}; + +// Price product in model +inline auto value(const Model &model, const Product &product, + // numerical parameters + const NumericalParam &num) { + // Random Number Generator + unique_ptr rng; + if (num.useSobol) + rng = make_unique(); + else + rng = make_unique(num.seed1, num.seed2); + + // Simulate + const auto resultMat = + num.parallel ? mcParallelSimul(product, model, *rng, num.numPath) + : mcSimul(product, model, *rng, num.numPath); + + // We return 2 vectors : the payoff identifiers and their values + struct { + vector identifiers; + vector values; + } results; + + const size_t nPayoffs = product.payoffLabels().size(); + results.identifiers = product.payoffLabels(); + results.values.resize(nPayoffs); + for (size_t i = 0; i < nPayoffs; ++i) { + results.values[i] = + accumulate(resultMat.begin(), resultMat.end(), 0.0, + [i](const double acc, const vector &v) { + return acc + v[i]; + }) / + num.numPath; + } + + return results; +} + +// Overload that picks product and model by name in the store +inline auto value(const string &modelId, const string &productId, + // numerical parameters + const NumericalParam &num) { + // Get model and product + const Model *model = getModel(modelId); + const Product *product = getProduct(productId); + + if (!model || !product) { + throw runtime_error("value() : Could not retrieve model and product"); + } + + return value(*model, *product, num); +} + +// AAD risk, one payoff +inline auto AADriskOne(const string &modelId, const string &productId, + const NumericalParam &num, + const string &riskPayoff = "") { + // Get model and product + const Model *model = getModel(modelId); + const Product *product = getProduct(productId); + + if (!model || !product) { + throw runtime_error("AADrisk() : Could not retrieve model and product"); + } + + // Random Number Generator + unique_ptr rng; + if (num.useSobol) + rng = make_unique(); + else + rng = make_unique(num.seed1, num.seed2); + + // Find the payoff for risk + size_t riskPayoffIdx = 0; + if (!riskPayoff.empty()) { + const vector &allPayoffs = product->payoffLabels(); + auto it = find(allPayoffs.begin(), allPayoffs.end(), riskPayoff); + if (it == allPayoffs.end()) { + throw runtime_error("AADriskOne() : payoff not found"); + } + riskPayoffIdx = distance(allPayoffs.begin(), it); + } + + // Simulate + const auto simulResults = + num.parallel + ? mcParallelSimulAAD(*product, *model, *rng, num.numPath, + [riskPayoffIdx](const vector &v) { + return v[riskPayoffIdx]; + }) + : mcSimulAAD(*product, *model, *rng, num.numPath, + [riskPayoffIdx](const vector &v) { + return v[riskPayoffIdx]; + }); + + // We return: a number and 2 vectors : + // - The payoff identifiers and their values + // - The value of the aggreagte payoff + // - The parameter idenitifiers + // - The sensititivities of the aggregate to parameters + struct { + vector payoffIds; + vector payoffValues; + double riskPayoffValue; + vector paramIds; + vector risks; + } results; + + const size_t nPayoffs = product->payoffLabels().size(); + results.payoffIds = product->payoffLabels(); + results.payoffValues.resize(nPayoffs); + for (size_t i = 0; i < nPayoffs; ++i) { + results.payoffValues[i] = + accumulate(simulResults.payoffs.begin(), simulResults.payoffs.end(), + 0.0, + [i](const double acc, const vector &v) { + return acc + v[i]; + }) / + num.numPath; + } + results.riskPayoffValue = accumulate(simulResults.aggregated.begin(), + simulResults.aggregated.end(), 0.0) / + num.numPath; + results.paramIds = model->parameterLabels(); + results.risks = std::move(simulResults.risks); + + return results; +} + +// AAD risk, aggregate portfolio +inline auto AADriskAggregate(const string &modelId, const string &productId, + const map ¬ionals, + const NumericalParam &num) { + // Get model and product + const Model *model = getModel(modelId); + const Product *product = getProduct(productId); + + if (!model || !product) { + throw runtime_error( + "AADriskAggregate() : Could not retrieve model and product"); + } + + // Random Number Generator + unique_ptr rng; + if (num.useSobol) + rng = make_unique(); + else + rng = make_unique(num.seed1, num.seed2); + + // Vector of notionals + const vector &allPayoffs = product->payoffLabels(); + vector vnots(allPayoffs.size(), 0.0); + for (const auto ¬ional : notionals) { + auto it = find(allPayoffs.begin(), allPayoffs.end(), notional.first); + if (it == allPayoffs.end()) { + throw runtime_error("AADriskAggregate() : payoff not found"); + } + vnots[distance(allPayoffs.begin(), it)] = notional.second; + } + + // Aggregator + auto aggregator = [&vnots](const vector &payoffs) { + // return inner_product(payoffs.begin(), payoffs.end(), vnots.begin(), + // Number(0.0)); -> deprecation in c++20 - custom + // implementation (from cppreference) + auto it_first = payoffs.begin(); + auto it_last = payoffs.end(); + auto it_vnots = vnots.begin(); + auto init = Number(0.0); + + while (it_first != it_last) { + init += (*it_first) * (*it_vnots); + ++it_first; + ++it_vnots; + } + return init; + }; + + // Simulate + const auto simulResults = + num.parallel + ? mcParallelSimulAAD(*product, *model, *rng, num.numPath, aggregator) + : mcSimulAAD(*product, *model, *rng, num.numPath, aggregator); + + // We return: a number and 2 vectors : + // - The payoff identifiers and their values + // - The value of the aggreagte payoff + // - The parameter idenitifiers + // - The sensititivities of the aggregate to parameters + struct { + vector payoffIds; + vector payoffValues; + double riskPayoffValue; + vector paramIds; + vector risks; + } results; + + const size_t nPayoffs = product->payoffLabels().size(); + results.payoffIds = product->payoffLabels(); + results.payoffValues.resize(nPayoffs); + for (size_t i = 0; i < nPayoffs; ++i) { + results.payoffValues[i] = + accumulate(simulResults.payoffs.begin(), simulResults.payoffs.end(), + 0.0, + [i](const double acc, const vector &v) { + return acc + v[i]; + }) / + num.numPath; + } + results.riskPayoffValue = accumulate(simulResults.aggregated.begin(), + simulResults.aggregated.end(), 0.0) / + num.numPath; + results.paramIds = model->parameterLabels(); + results.risks = std::move(simulResults.risks); + + return results; +} + +// Returns a vector of values and a matrix of risks +// with payoffs in columns and parameters in rows +// along with ids of payoffs and parameters + +struct RiskReports { + vector payoffs; + vector params; + vector values; + matrix risks; +}; + +// Itemized AAD risk, one per payoff +inline RiskReports AADriskMulti(const string &modelId, const string &productId, + const NumericalParam &num) { + const Model *model = getModel(modelId); + const Product *product = getProduct(productId); + + if (!model || !product) { + throw runtime_error("AADrisk() : Could not retrieve model and product"); + } + + RiskReports results; + + // Random Number Generator + unique_ptr rng; + if (num.useSobol) + rng = make_unique(); + else + rng = make_unique(num.seed1, num.seed2); + + // Simulate + const auto simulResults = + num.parallel + ? mcParallelSimulAADMulti(*product, *model, *rng, num.numPath) + : mcSimulAADMulti(*product, *model, *rng, num.numPath); + + results.params = model->parameterLabels(); + results.payoffs = product->payoffLabels(); + results.risks = std::move(simulResults.risks); + + // Average values across paths + const size_t nPayoffs = product->payoffLabels().size(); + results.values.resize(nPayoffs); + for (size_t i = 0; i < nPayoffs; ++i) { + results.values[i] = + accumulate(simulResults.payoffs.begin(), simulResults.payoffs.end(), + 0.0, + [i](const double acc, const vector &v) { + return acc + v[i]; + }) / + num.numPath; + } + + return results; +} + +// Bump risk, itemized +// Same result format as AADriskMulti() +inline RiskReports bumpRisk(const string &modelId, const string &productId, + const NumericalParam &num) { + auto *orig = getModel(modelId); + const Product *product = getProduct(productId); + + if (!orig || !product) { + throw runtime_error("bumpRisk() : Could not retrieve model and product"); + } + + RiskReports results; + + // base values + auto baseRes = value(*orig, *product, num); + results.payoffs = baseRes.identifiers; + results.values = baseRes.values; + + // make copy so we don't modify the model in memory + auto model = orig->clone(); + + results.params = model->parameterLabels(); + const vector parameters = model->parameters(); + const size_t n = parameters.size(), m = results.payoffs.size(); + results.risks.resize(n, m); + + // bumps + for (size_t i = 0; i < n; ++i) { + *parameters[i] += 1.e-08; + auto bumpRes = value(*model, *product, num); + *parameters[i] -= 1.e-08; + + for (size_t j = 0; j < m; ++j) { + results.risks[i][j] = 1.0e+08 * (bumpRes.values[j] - baseRes.values[j]); + } + } + + return results; +} + +// Dupire specific + +// Returns a struct with price, delta and vega matrix +inline auto dupireAADRisk( + // model id + const string &modelId, + // product id + const string &productId, const map ¬ionals, + // numerical parameters + const NumericalParam &num) { + // Check that the model is a Dupire + const Model *model = getModel(modelId); + if (!model) { + throw runtime_error("dupireAADRisk() : Model not found"); + } + const Dupire *dupire = dynamic_cast *>(model); + if (!dupire) { + throw runtime_error("dupireAADRisk() : Model not a Dupire"); + } + + // Results + struct { + double value; + double delta; + matrix vega; + } results; + + // Go + auto simulResults = AADriskAggregate(modelId, productId, notionals, num); + + // Find results + + // Value + results.value = simulResults.riskPayoffValue; + + // Delta + + results.delta = simulResults.risks[0]; + + // Vegas + + results.vega.resize(dupire->spots().size(), dupire->times().size()); + copy(next(simulResults.risks.begin()), simulResults.risks.end(), + results.vega.begin()); + + return results; +} + +// Returns spots, times and lVols in a struct +inline auto dupireCalib( + // The local vol grid + // The spots to include + const vector &inclSpots, + // Maximum space between spots + const double maxDs, + // The times to include, note NOT 0 + const vector