Skip to content

Commit

Permalink
Fix datatypes for consistent Real usage, keeping AD compatibility (#2116
Browse files Browse the repository at this point in the history
)
  • Loading branch information
lballabio authored Nov 18, 2024
2 parents e8e8b06 + 565a570 commit a07bf7b
Show file tree
Hide file tree
Showing 8 changed files with 27 additions and 27 deletions.
4 changes: 2 additions & 2 deletions ql/math/matrixutilities/choleskydecomposition.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -71,13 +71,13 @@ namespace QuantLib {

Array x(n);
for (Size i=0; i < n; ++i) {
x[i] = -std::inner_product(L.row_begin(i), L.row_begin(i)+i, x.begin(), -b[i]);
x[i] = -std::inner_product(L.row_begin(i), L.row_begin(i)+i, x.begin(), Real(-b[i]));
x[i] /= L[i][i];
}

for (Integer i=n-1; i >=0; --i) {
x[i] = -std::inner_product(
L.column_begin(i)+i+1, L.column_end(i), x.begin()+i+1, -x[i]);
L.column_begin(i)+i+1, L.column_end(i), x.begin()+i+1, Real(-x[i]));
x[i] /= L[i][i];
}

Expand Down
4 changes: 2 additions & 2 deletions ql/math/randomnumbers/zigguratgaussianrng.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -91,10 +91,10 @@ namespace QuantLib {
// This saves us generating a whole extra random number, while the added
// precision of using 64 bits for double does not buy us much.
std::uint64_t randomU64 = uint64Generator_.nextInt64();
auto u = 2.0 * (Real(randomU64 >> 11) + 0.5) * (1.0 / Real(1ULL << 53)) - 1.0;
Real u = 2.0 * (Real(randomU64 >> 11) + 0.5) * (1.0 / Real(1ULL << 53)) - 1.0;
auto i = (int)(randomU64 & 0xff);

auto x = u * normX(i);
Real x = u * normX(i);

if (std::abs(x) < normX(i + 1)) {
return x;
Expand Down
10 changes: 5 additions & 5 deletions ql/pricingengines/basket/choibasketengine.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ namespace QuantLib {
//q1 = inverse(C)*vStar1;
for (Size i=0; i < n_; ++i)
q1[i] = (vStar1[i] - std::inner_product(
C.row_begin(i), C.row_begin(i) + i, q1.begin(), 0.0))/C[i][i];
C.row_begin(i), C.row_begin(i) + i, q1.begin(), Real(0.0)))/C[i][i];

vStar1 /= Norm2(q1);
}
Expand Down Expand Up @@ -190,7 +190,7 @@ namespace QuantLib {
Array vq(n_);
for (Size i=0; i < n_; ++i)
vq[i] = 0.5*std::accumulate(
v.row_begin(i), v.row_end(i), 0.0,
v.row_begin(i), v.row_end(i), Real(0.0),
[](Real acc, Real x) -> Real { return acc + x*x; }
);

Expand Down Expand Up @@ -230,7 +230,7 @@ namespace QuantLib {
const auto deltaPricer = [&](const Array& z) -> Real {
const Real d = dStore[dStoreCounter++];
const Real vz = std::inner_product(
v.row_begin(k), v.row_end(k), z.begin(), 0.0);
v.row_begin(k), v.row_end(k), z.begin(), Real(0.0));
const Real f = std::exp(-M_SQRT2*vz - vq[k]);

return std::exp(-DotProduct(z, z)) * f * N(d + vStar1[k]);
Expand All @@ -246,15 +246,15 @@ namespace QuantLib {
for (Size k=0; k < n_; ++k) {
const auto fHatPricer = [&](const Array& z) -> Real {
const Real vz = std::inner_product(
v.row_begin(k), v.row_end(k), z.begin(), 0.0);
v.row_begin(k), v.row_end(k), z.begin(), Real(0.0));
const Real f = std::exp(-M_SQRT2*vz - vq[k]);

return std::exp(-DotProduct(z, z)) * f;
};
fHat[k] = ghq(fHatPricer) * normFactor;
}
const Array cv = fwdDelta*fwd*(fHat-1.0);
results_.value -= std::accumulate(cv.begin(), cv.end(), 0.0);
results_.value -= std::accumulate(cv.begin(), cv.end(), Real(0.0));
}
}
}
Expand Down
14 changes: 7 additions & 7 deletions ql/pricingengines/basket/denglizhoubasketengine.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -126,7 +126,7 @@ namespace QuantLib {
[](const Real& value, const auto& p) -> Real {
return value + std::get<0>(p)*std::get<2>(p);
});
const Real F0 = std::accumulate(F.begin(), F.end(), 0.0);
const Real F0 = std::accumulate(F.begin(), F.end(), Real(0.0));
const DiscountFactor dq_S0 = F0/S0*dr0;

Real v_s = 0.0;
Expand Down Expand Up @@ -174,26 +174,26 @@ namespace QuantLib {
results_.value = std::max(0.0, callValue);
else {
const Real fwd = _s[0]*_dq[0] - dr0*strike
- std::inner_product(_s.begin()+1, _s.end(), _dq.begin()+1, 0.0);
- std::inner_product(_s.begin()+1, _s.end(), _dq.begin()+1, Real(0.0));
results_.value = std::max(0.0, callValue - fwd);
}
}

Real DengLiZhouBasketEngine::I(Real u, Real tF2, const Matrix& D, const Matrix& DF, Size i) {
const Real psi = 1.0/
(1.0 + std::inner_product(
D.row_begin(i), D.row_end(i), D.row_begin(i), 0.0));
D.row_begin(i), D.row_end(i), D.row_begin(i), Real(0.0)));
const Real sqrtPsi = std::sqrt(psi);

const Real n_uSqrtPsi = NormalDistribution()(u*sqrtPsi);
const Real J_0 = CumulativeNormalDistribution()(u*sqrtPsi);

const Real vFv = std::inner_product(
DF.row_begin(i), DF.row_end(i), D.row_begin(i), 0.0);
DF.row_begin(i), DF.row_end(i), D.row_begin(i), Real(0.0));
const Real J_1 = psi*sqrtPsi*(psi*u*u - 1.0) * vFv * n_uSqrtPsi;

const Real vFFv = std::inner_product(
DF.row_begin(i), DF.row_end(i), DF.row_begin(i), 0.0);
DF.row_begin(i), DF.row_end(i), DF.row_begin(i), Real(0.0));
const Real J_2 = u*psi*sqrtPsi*n_uSqrtPsi*(
2 * tF2
+ vFv*vFv*(squared(squared(psi*u))
Expand Down Expand Up @@ -235,7 +235,7 @@ namespace QuantLib {
for (Size i=1; i <= N; ++i)
for (Size j=i; j <= N; ++j)
E(i-1, j-1) = E(j-1, i-1) =
a*(((i==j)? squared(nu[j])*std::exp(mu[j])/(nu[0]*(R+K)) : 0.0)
a*(((i==j)? squared(nu[j])*std::exp(mu[j])/(nu[0]*(R+K)) : Real(0.0))
-nu[i]*nu[j]*std::exp(mu[i]+mu[j])/(nu[0]*squared(R + K)) );

const Matrix F = sqSig11*E*sqSig11;
Expand Down Expand Up @@ -267,7 +267,7 @@ namespace QuantLib {
for (Size k=1; k < N+1; ++k)
C[k] = c + trF + nu[k]*sig11d[k-1] + squared(nu[k])
* std::inner_product(sig11.row_begin(k-1), sig11.row_end(k-1),
Esig11.column_begin(k-1), 0.0);
Esig11.column_begin(k-1), Real(0.0));

std::vector<Array> D(N+2);
D[0] = sqSig11*(d + 2*nu[0]*Esig10);
Expand Down
10 changes: 5 additions & 5 deletions ql/pricingengines/basket/singlefactorbsmbasketengine.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -99,12 +99,12 @@ namespace QuantLib {
"non-positive strikes only allowed for spread options");

// linear approximation
const Real denom = std::accumulate(attr.begin(), attr.end(), 0.0);
const Real denom = std::accumulate(attr.begin(), attr.end(), Real(0.0));
const Real xInit = (std::abs(denom) > 1000*QL_EPSILON)
? std::min(10.0, std::max(-10.0,
(K_ - std::accumulate(a_.begin(), a_.end(), 0.0))/denom)
(K_ - std::accumulate(a_.begin(), a_.end(), Real(0.0)))/denom)
)
: 0.0;
: Real(0.0);

switch(strategy) {
case Brent:
Expand Down Expand Up @@ -165,7 +165,7 @@ namespace QuantLib {
[](Real x) -> bool { return close_enough(x, 0.0); }
)) {
results_.value = dr0*payoff->operator()(
std::accumulate(fwdBasket.begin(), fwdBasket.end(), 0.0));
std::accumulate(fwdBasket.begin(), fwdBasket.end(), Real(0.0)));
}
else {
const Real d = -detail::SumExponentialsRootSolver(
Expand All @@ -178,7 +178,7 @@ namespace QuantLib {
results_.value = cp * dr0 *
std::inner_product(
fwdBasket.begin(), fwdBasket.end(), stdDev.begin(),
-strike*N(cp*d),
Real(-strike*N(cp*d)),
std::plus<>(),
[d, cp, &N](Real x, Real y) -> Real { return x*N(cp*(d+y)); }
);
Expand Down
4 changes: 2 additions & 2 deletions test-suite/basketoption.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1790,7 +1790,7 @@ BOOST_AUTO_TEST_CASE(testRootOfSumExponentials) {
const Real offset = (mt.nextReal() < 0.3)? -1.0 : 0.0;
for (Size j=0; j < n; ++j) {
a[j] = mt.nextReal() + offset;
sig[j] = copysign(1.0, a[j])*mt.nextReal();
sig[j] = std::copysign(Real(1.0), a[j])*mt.nextReal();
}
const Real kMin = detail::SumExponentialsRootSolver(a, sig, 0.0)(-10.0);
const Real kMax = detail::SumExponentialsRootSolver(a, sig, 0.0)( 10.0);
Expand Down Expand Up @@ -1868,7 +1868,7 @@ BOOST_AUTO_TEST_CASE(testSingleFactorBsmBasketEngine) {
);

const Real strike = std::inner_product(
t.weights.begin(), t.weights.end(), t.underlyings.begin(), 0.0
t.weights.begin(), t.weights.end(), t.underlyings.begin(), Real(0.0)
);

const ext::shared_ptr<PlainVanillaPayoff> payoff
Expand Down
4 changes: 2 additions & 2 deletions test-suite/callablebonds.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -655,8 +655,8 @@ BOOST_AUTO_TEST_CASE(testSnappingExerciseDate2ClosestCouponDate) {
<< std::setprecision(1) << tolerance);
}

auto cleanPrice = callableBond->cleanPrice() - 2.0;
auto oas = callableBond->OAS(cleanPrice, termStructure, accrualDCC,
Real cleanPrice = callableBond->cleanPrice() - 2.0;
Real oas = callableBond->OAS(cleanPrice, termStructure, accrualDCC,
QuantLib::Continuous, frequency);
if (prevOAS - oas < expectedOasStep) {
BOOST_ERROR("failed to get expected change in OAS at "
Expand Down
4 changes: 2 additions & 2 deletions test-suite/matrices.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -894,7 +894,7 @@ BOOST_AUTO_TEST_CASE(testCholeskySolverFor) {

const Array diff = Abs(rho*x - b);

BOOST_CHECK_SMALL(std::sqrt(DotProduct(diff, diff)), 20*std::sqrt(n)*QL_EPSILON);
QL_CHECK_SMALL(std::sqrt(DotProduct(diff, diff)), 20*std::sqrt(n)*QL_EPSILON);
}
}

Expand All @@ -917,7 +917,7 @@ namespace {
const Array& actual, const Array& expected, Real tol) {
BOOST_REQUIRE(actual.size() == expected.size());
for (auto i = 0U; i < actual.size(); i++) {
BOOST_CHECK_SMALL(actual[i] - expected[i], tol);
QL_CHECK_SMALL(actual[i] - expected[i], tol);
}
}
}
Expand Down

0 comments on commit a07bf7b

Please sign in to comment.