Skip to content

Commit e239f95

Browse files
pcasperspcaspers
committed
Merge branch 'master' into 'master'
merge from master Closes QPR-12184 and QPR-12140 See merge request qs/quantlib-github!1
2 parents 55f433c + db9a28c commit e239f95

13 files changed

+533
-71
lines changed

‎ql/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -416,6 +416,7 @@ set(QL_SOURCES
416416
math/polynomialmathfunction.cpp
417417
math/primenumbers.cpp
418418
math/quadratic.cpp
419+
math/randomnumbers/burley2020sobolrsg.cpp
419420
math/randomnumbers/faurersg.cpp
420421
math/randomnumbers/haltonrsg.cpp
421422
math/randomnumbers/knuthuniformrng.cpp
@@ -1512,6 +1513,7 @@ set(QL_HEADERS
15121513
math/quadratic.hpp
15131514
math/randomnumbers/boxmullergaussianrng.hpp
15141515
math/randomnumbers/centrallimitgaussianrng.hpp
1516+
math/randomnumbers/burley2020sobolrsg.hpp
15151517
math/randomnumbers/faurersg.hpp
15161518
math/randomnumbers/haltonrsg.hpp
15171519
math/randomnumbers/inversecumulativerng.hpp

‎ql/math/randomnumbers/burley2020sobolrsg.cpp

Lines changed: 151 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,151 @@
1+
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2+
3+
/*
4+
Copyright (C) 2023 Peter Caspers
5+
6+
This file is part of QuantLib, a free-software/open-source library
7+
for financial quantitative analysts and developers - http://quantlib.org/
8+
9+
QuantLib is free software: you can redistribute it and/or modify it
10+
under the terms of the QuantLib license. You should have received a
11+
copy of the license along with this program; if not, please email
12+
<quantlib-dev@lists.sf.net>. The license is also available online at
13+
<http://quantlib.org/license.shtml>.
14+
15+
This program is distributed in the hope that it will be useful, but WITHOUT
16+
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17+
FOR A PARTICULAR PURPOSE. See the license for more details.
18+
*/
19+
20+
#include <ql/math/randomnumbers/burley2020sobolrsg.hpp>
21+
#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
22+
23+
namespace QuantLib {
24+
25+
Burley2020SobolRsg::Burley2020SobolRsg(Size dimensionality,
26+
unsigned long seed,
27+
SobolRsg::DirectionIntegers directionIntegers,
28+
unsigned long scrambleSeed)
29+
: dimensionality_(dimensionality), seed_(seed), directionIntegers_(directionIntegers),
30+
integerSequence_(dimensionality), sequence_(std::vector<Real>(dimensionality), 1.0) {
31+
reset();
32+
group4Seeds_.resize((dimensionality_ - 1) / 4 + 1);
33+
MersenneTwisterUniformRng mt(scrambleSeed);
34+
for (auto& s : group4Seeds_) {
35+
s = static_cast<std::uint32_t>(mt.nextInt32());
36+
}
37+
}
38+
39+
void Burley2020SobolRsg::reset() const {
40+
sobolRsg_ = ext::make_shared<SobolRsg>(dimensionality_, seed_, directionIntegers_, false);
41+
nextSequenceCounter_ = 0;
42+
}
43+
44+
const std::vector<std::uint32_t>& Burley2020SobolRsg::skipTo(std::uint32_t n) const {
45+
reset();
46+
for (Size k = 0; k < n + 1; ++k) {
47+
nextInt32Sequence();
48+
}
49+
return integerSequence_;
50+
}
51+
52+
namespace {
53+
54+
// for reverseBits() see http://graphics.stanford.edu/~seander/bithacks.html#BitReverseTable
55+
56+
static const std::uint8_t bitReverseTable[] = {
57+
0u, 128u, 64u, 192u, 32u, 160u, 96u, 224u, 16u, 144u, 80u, 208u, 48u, 176u,
58+
112u, 240u, 8u, 136u, 72u, 200u, 40u, 168u, 104u, 232u, 24u, 152u, 88u, 216u,
59+
56u, 184u, 120u, 248u, 4u, 132u, 68u, 196u, 36u, 164u, 100u, 228u, 20u, 148u,
60+
84u, 212u, 52u, 180u, 116u, 244u, 12u, 140u, 76u, 204u, 44u, 172u, 108u, 236u,
61+
28u, 156u, 92u, 220u, 60u, 188u, 124u, 252u, 2u, 130u, 66u, 194u, 34u, 162u,
62+
98u, 226u, 18u, 146u, 82u, 210u, 50u, 178u, 114u, 242u, 10u, 138u, 74u, 202u,
63+
42u, 170u, 106u, 234u, 26u, 154u, 90u, 218u, 58u, 186u, 122u, 250u, 6u, 134u,
64+
70u, 198u, 38u, 166u, 102u, 230u, 22u, 150u, 86u, 214u, 54u, 182u, 118u, 246u,
65+
14u, 142u, 78u, 206u, 46u, 174u, 110u, 238u, 30u, 158u, 94u, 222u, 62u, 190u,
66+
126u, 254u, 1u, 129u, 65u, 193u, 33u, 161u, 97u, 225u, 17u, 145u, 81u, 209u,
67+
49u, 177u, 113u, 241u, 9u, 137u, 73u, 201u, 41u, 169u, 105u, 233u, 25u, 153u,
68+
89u, 217u, 57u, 185u, 121u, 249u, 5u, 133u, 69u, 197u, 37u, 165u, 101u, 229u,
69+
21u, 149u, 85u, 213u, 53u, 181u, 117u, 245u, 13u, 141u, 77u, 205u, 45u, 173u,
70+
109u, 237u, 29u, 157u, 93u, 221u, 61u, 189u, 125u, 253u, 3u, 131u, 67u, 195u,
71+
35u, 163u, 99u, 227u, 19u, 147u, 83u, 211u, 51u, 179u, 115u, 243u, 11u, 139u,
72+
75u, 203u, 43u, 171u, 107u, 235u, 27u, 155u, 91u, 219u, 59u, 187u, 123u, 251u,
73+
7u, 135u, 71u, 199u, 39u, 167u, 103u, 231u, 23u, 151u, 87u, 215u, 55u, 183u,
74+
119u, 247u, 15u, 143u, 79u, 207u, 47u, 175u, 111u, 239u, 31u, 159u, 95u, 223u,
75+
63u, 191u, 127u, 255u};
76+
77+
inline std::uint32_t reverseBits(std::uint32_t x) {
78+
return (bitReverseTable[x & 0xff] << 24) | (bitReverseTable[(x >> 8) & 0xff] << 16) |
79+
(bitReverseTable[(x >> 16) & 0xff] << 8) | (bitReverseTable[(x >> 24) & 0xff]);
80+
}
81+
82+
inline std::uint32_t laine_karras_permutation(std::uint32_t x, std::uint32_t seed) {
83+
x += seed;
84+
x ^= x * 0x6c50b47cu;
85+
x ^= x * 0xb82f1e52u;
86+
x ^= x * 0xc7afe638u;
87+
x ^= x * 0x8d22f6e6u;
88+
return x;
89+
}
90+
91+
inline std::uint32_t nested_uniform_scramble(std::uint32_t x, std::uint32_t seed) {
92+
x = reverseBits(x);
93+
x = laine_karras_permutation(x, seed);
94+
x = reverseBits(x);
95+
return x;
96+
}
97+
98+
// the results depend a lot on the details of the hash_combine() function that is used
99+
// we use hash_combine() calling hash(), hash_mix() as implemented here:
100+
// https://github.com/boostorg/container_hash/blob/boost-1.83.0/include/boost/container_hash/hash.hpp#L560
101+
// https://github.com/boostorg/container_hash/blob/boost-1.83.0/include/boost/container_hash/hash.hpp#L115
102+
// https://github.com/boostorg/container_hash/blob/boost-1.83.0/include/boost/container_hash/detail/hash_mix.hpp#L67
103+
104+
inline std::uint64_t local_hash_mix(std::uint64_t x) {
105+
const std::uint64_t m = 0xe9846af9b1a615d;
106+
x ^= x >> 32;
107+
x *= m;
108+
x ^= x >> 32;
109+
x *= m;
110+
x ^= x >> 28;
111+
return x;
112+
}
113+
114+
inline std::uint64_t local_hash(const std::uint64_t v) {
115+
std::uint64_t seed = 0;
116+
seed = (v >> 32) + local_hash_mix(seed);
117+
seed = (v & 0xFFFFFFFF) + local_hash_mix(seed);
118+
return seed;
119+
}
120+
121+
inline std::uint64_t local_hash_combine(std::uint64_t x, const uint64_t v) {
122+
return local_hash_mix(x + 0x9e3779b9 + local_hash(v));
123+
}
124+
}
125+
126+
const std::vector<std::uint32_t>& Burley2020SobolRsg::nextInt32Sequence() const {
127+
auto n = nested_uniform_scramble(nextSequenceCounter_, group4Seeds_[0]);
128+
const auto& seq = sobolRsg_->skipTo(n);
129+
std::copy(seq.begin(), seq.end(), integerSequence_.begin());
130+
Size i = 0, group = 0;
131+
do {
132+
std::uint64_t seed = group4Seeds_[group++];
133+
for (Size g = 0; g < 4 && i < dimensionality_; ++g, ++i) {
134+
seed = local_hash_combine(seed, g);
135+
integerSequence_[i] =
136+
nested_uniform_scramble(integerSequence_[i], static_cast<std::uint32_t>(seed));
137+
}
138+
} while (i < dimensionality_);
139+
++nextSequenceCounter_;
140+
return integerSequence_;
141+
}
142+
143+
const SobolRsg::sample_type& Burley2020SobolRsg::nextSequence() const {
144+
const std::vector<std::uint32_t>& v = nextInt32Sequence();
145+
// normalize to get a double in (0,1)
146+
for (Size k = 0; k < dimensionality_; ++k) {
147+
sequence_.value[k] = static_cast<double>(v[k]) / 4294967296.0;
148+
}
149+
return sequence_;
150+
}
151+
}

‎ql/math/randomnumbers/burley2020sobolrsg.hpp

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2+
3+
/*
4+
Copyright (C) 2023 Peter Caspers
5+
6+
This file is part of QuantLib, a free-software/open-source library
7+
for financial quantitative analysts and developers - http://quantlib.org/
8+
9+
QuantLib is free software: you can redistribute it and/or modify it
10+
under the terms of the QuantLib license. You should have received a
11+
copy of the license along with this program; if not, please email
12+
<quantlib-dev@lists.sf.net>. The license is also available online at
13+
<http://quantlib.org/license.shtml>.
14+
15+
This program is distributed in the hope that it will be useful, but WITHOUT
16+
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17+
FOR A PARTICULAR PURPOSE. See the license for more details.
18+
*/
19+
20+
/*! \file burley2020sobolrsg.hpp
21+
\brief scrambled Sobol sequence following Burley, 2020
22+
*/
23+
24+
#ifndef quantlib_burley2020_scrambled_sobolrsg_hpp
25+
#define quantlib_burley2020_scrambled_sobolrsg_hpp
26+
27+
#include <ql/math/randomnumbers/sobolrsg.hpp>
28+
#include <ql/shared_ptr.hpp>
29+
30+
namespace QuantLib {
31+
32+
//! Scrambled sobol sequence according to Burley, 2020
33+
/*! Reference: Brent Burley: Practical Hash-based Owen Scrambling,
34+
Journal of Computer Graphics Techniques, Vol. 9, No. 4, 2020 */
35+
class Burley2020SobolRsg {
36+
public:
37+
typedef Sample<std::vector<Real>> sample_type;
38+
explicit Burley2020SobolRsg(
39+
Size dimensionality,
40+
unsigned long seed = 42,
41+
SobolRsg::DirectionIntegers directionIntegers = SobolRsg::Jaeckel,
42+
unsigned long scrambleSeed = 43);
43+
const std::vector<std::uint32_t>& skipTo(std::uint32_t n) const;
44+
const std::vector<std::uint32_t>& nextInt32Sequence() const;
45+
const SobolRsg::sample_type& nextSequence() const;
46+
const sample_type& lastSequence() const { return sequence_; }
47+
Size dimension() const { return dimensionality_; }
48+
49+
private:
50+
void reset() const;
51+
Size dimensionality_;
52+
unsigned long seed_;
53+
SobolRsg::DirectionIntegers directionIntegers_;
54+
mutable ext::shared_ptr<SobolRsg> sobolRsg_;
55+
mutable std::vector<std::uint32_t> integerSequence_;
56+
mutable sample_type sequence_;
57+
mutable std::uint32_t nextSequenceCounter_;
58+
mutable std::vector<std::uint32_t> group4Seeds_;
59+
};
60+
61+
}
62+
63+
64+
#endif

‎ql/math/randomnumbers/sobolbrownianbridgersg.cpp

Lines changed: 47 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -25,27 +25,31 @@
2525
#include <ql/math/randomnumbers/sobolbrownianbridgersg.hpp>
2626

2727
namespace QuantLib {
28-
SobolBrownianBridgeRsg::SobolBrownianBridgeRsg(
29-
Size factors, Size steps,
30-
SobolBrownianGenerator::Ordering ordering,
31-
unsigned long seed,
32-
SobolRsg::DirectionIntegers directionIntegers)
33-
: factors_(factors), steps_(steps), dim_(factors*steps),
34-
seq_(sample_type::value_type(factors*steps), 1.0),
35-
gen_(factors, steps, ordering, seed, directionIntegers) {
28+
29+
namespace {
30+
void setNextSequence(SobolBrownianGeneratorBase& gen, std::vector<Real>& seq) {
31+
gen.nextPath();
32+
std::vector<Real> output(gen.numberOfFactors());
33+
for (Size i = 0; i < gen.numberOfSteps(); ++i) {
34+
gen.nextStep(output);
35+
std::copy(output.begin(), output.end(), seq.begin() + i * gen.numberOfFactors());
36+
}
37+
}
3638
}
3739

40+
SobolBrownianBridgeRsg::SobolBrownianBridgeRsg(Size factors,
41+
Size steps,
42+
SobolBrownianGenerator::Ordering ordering,
43+
unsigned long seed,
44+
SobolRsg::DirectionIntegers directionIntegers)
45+
: seq_(sample_type::value_type(factors * steps), 1.0),
46+
gen_(factors, steps, ordering, seed, directionIntegers) {}
47+
3848
const SobolBrownianBridgeRsg::sample_type&
3949
SobolBrownianBridgeRsg::nextSequence() const {
40-
gen_.nextPath();
41-
std::vector<Real> output(factors_);
42-
for (Size i=0; i < steps_; ++i) {
43-
gen_.nextStep(output);
44-
std::copy(output.begin(), output.end(),
45-
seq_.value.begin()+i*factors_);
46-
}
47-
50+
setNextSequence(gen_, seq_.value);
4851
return seq_;
52+
4953
}
5054

5155
const SobolBrownianBridgeRsg::sample_type&
@@ -54,6 +58,32 @@ namespace QuantLib {
5458
}
5559

5660
Size SobolBrownianBridgeRsg::dimension() const {
57-
return dim_;
61+
return gen_.numberOfFactors() * gen_.numberOfSteps();
62+
}
63+
64+
Burley2020SobolBrownianBridgeRsg::Burley2020SobolBrownianBridgeRsg(
65+
Size factors,
66+
Size steps,
67+
SobolBrownianGenerator::Ordering ordering,
68+
unsigned long seed,
69+
SobolRsg::DirectionIntegers directionIntegers,
70+
unsigned long scrambleSeed)
71+
: seq_(sample_type::value_type(factors * steps), 1.0),
72+
gen_(factors, steps, ordering, seed, directionIntegers, scrambleSeed) {}
73+
74+
const Burley2020SobolBrownianBridgeRsg::sample_type&
75+
Burley2020SobolBrownianBridgeRsg::nextSequence() const {
76+
setNextSequence(gen_, seq_.value);
77+
return seq_;
78+
}
79+
80+
const Burley2020SobolBrownianBridgeRsg::sample_type&
81+
Burley2020SobolBrownianBridgeRsg::lastSequence() const {
82+
return seq_;
5883
}
84+
85+
Size Burley2020SobolBrownianBridgeRsg::dimension() const {
86+
return gen_.numberOfFactors() * gen_.numberOfSteps();
87+
}
88+
5989
}

‎ql/math/randomnumbers/sobolbrownianbridgersg.hpp

Lines changed: 21 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -45,10 +45,30 @@ namespace QuantLib {
4545
Size dimension() const;
4646

4747
private:
48-
const Size factors_, steps_, dim_;
4948
mutable sample_type seq_;
5049
mutable SobolBrownianGenerator gen_;
5150
};
51+
52+
class Burley2020SobolBrownianBridgeRsg {
53+
public:
54+
typedef Sample<std::vector<Real> > sample_type;
55+
56+
Burley2020SobolBrownianBridgeRsg(
57+
Size factors,
58+
Size steps,
59+
SobolBrownianGenerator::Ordering ordering = SobolBrownianGenerator::Diagonal,
60+
unsigned long seed = 42,
61+
SobolRsg::DirectionIntegers directionIntegers = SobolRsg::JoeKuoD7,
62+
unsigned long scrambleSeed = 43);
63+
64+
const sample_type& nextSequence() const;
65+
const sample_type& lastSequence() const;
66+
Size dimension() const;
67+
68+
private:
69+
mutable sample_type seq_;
70+
mutable Burley2020SobolBrownianGenerator gen_;
71+
};
5272
}
5373

5474
#endif

0 commit comments

Comments
 (0)