From 6518198db61c145f5fd883a1e05e90d09b7b429e Mon Sep 17 00:00:00 2001 From: Wojciech Jarosz Date: Mon, 18 Dec 2023 00:38:46 +1100 Subject: [PATCH] removing SOA; padding Jittered + LP sampler --- CMakeLists.txt | 2 - include/sampler/GrayCode.h | 15 +--- include/sampler/Jittered.h | 47 +++++++----- include/sampler/LP.h | 2 +- include/sampler/MultiJittered.h | 12 ++- include/sampler/SOA.h | 78 -------------------- src/app.cpp | 4 +- src/sampler/GrayCode.cpp | 30 +++++++- src/sampler/Jittered.cpp | 39 +++++----- src/sampler/LP.cpp | 17 +++-- src/sampler/MultiJittered.cpp | 3 +- src/sampler/SOA.cpp | 126 -------------------------------- 12 files changed, 101 insertions(+), 274 deletions(-) delete mode 100644 include/sampler/SOA.h delete mode 100644 src/sampler/SOA.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index dee4c7d..cb33753 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -272,7 +272,6 @@ add_library( include/sampler/Random.h include/sampler/RandomPermutation.h include/sampler/Sampler.h - include/sampler/SOA.h include/sampler/Sobol.h include/sampler/Sudoku.h include/sampler/xi.h @@ -294,7 +293,6 @@ add_library( src/sampler/OABush.cpp src/sampler/OACMJND.cpp src/sampler/Random.cpp - src/sampler/SOA.cpp src/sampler/Sobol.cpp src/sampler/Sudoku.cpp src/sampler/XiSequence.cpp diff --git a/include/sampler/GrayCode.h b/include/sampler/GrayCode.h index 6764570..cf849a0 100644 --- a/include/sampler/GrayCode.h +++ b/include/sampler/GrayCode.h @@ -3,7 +3,6 @@ */ #pragma once -#include #include #include #include @@ -18,7 +17,7 @@ class GrayCode : public TSamplerDim<2> { public: - GrayCode(unsigned n = 3); + GrayCode(unsigned n = 2); void sample(float[], unsigned i) override; @@ -31,17 +30,7 @@ class GrayCode : public TSamplerDim<2> { return N; } - int setNumSamples(unsigned num) override - { - int log2N = std::round(std::log2(num)); - if (log2N & 1) - log2N++; // Only allow even powers of 2. - N = 1 << log2N; // Make N a power of 4, if not. - log2n = log2N / 2; - n = 1 << log2n; - regenerate(); - return N; - } + int setNumSamples(unsigned num) override; bool randomized() const override { diff --git a/include/sampler/Jittered.h b/include/sampler/Jittered.h index ea7a0ec..c8fddf2 100644 --- a/include/sampler/Jittered.h +++ b/include/sampler/Jittered.h @@ -7,35 +7,41 @@ #include /// Encapsulate a 2D stratified or "jittered" point set. -class Jittered : public TSamplerDim<2> +class Jittered : public TSamplerMinMaxDim<1, 1024> { public: Jittered(unsigned resX, unsigned resY, float jitter = 1.0f); - void reset() override; void sample(float[], unsigned i) override; - std::string name() const override + unsigned dimensions() const override { - return "Jittered"; + return m_numDimensions; + } + void setDimensions(unsigned n) override + { + m_numDimensions = n; } int numSamples() const override { - return m_resX * m_resY; + return m_numSamples; } int setNumSamples(unsigned n) override { + if (n == m_numSamples) + return m_numSamples; int sqrtVal = (n == 0) ? 1 : (int)(std::sqrt((float)n) + 0.5f); - m_resX = m_resY = sqrtVal; - reset(); - return m_resX * m_resY; + setNumSamples(sqrtVal, sqrtVal); + return m_numSamples; } void setNumSamples(unsigned x, unsigned y) { - m_resX = x; - m_resY = y; - reset(); + m_resX = x == 0 ? 1 : x; + m_resY = y == 0 ? 1 : y; + m_numSamples = m_resX * m_resY; + m_xScale = 1.0f / m_resX; + m_yScale = 1.0f / m_resY; } bool randomized() const override @@ -44,8 +50,9 @@ class Jittered : public TSamplerDim<2> } void setRandomized(bool r = true) override { - if ((m_randomize = r)) - m_seed++; + m_randomize = r; + m_rand.seed(m_permutation); + m_permutation = r ? m_rand.nextUInt() : 0; } /// Gets/sets how much the points are jittered @@ -58,12 +65,18 @@ class Jittered : public TSamplerDim<2> return m_maxJit = j; } + std::string name() const override + { + return "Jittered"; + } + private: - unsigned m_resX, m_resY; + unsigned m_resX, m_resY, m_numSamples, m_numDimensions; float m_maxJit; - - float m_xScale, m_yScale; bool m_randomize = true; pcg32 m_rand; - unsigned m_seed = 13; + unsigned m_seed = 13; + unsigned m_permutation = 13; + + float m_xScale, m_yScale; }; diff --git a/include/sampler/LP.h b/include/sampler/LP.h index 1f6a79f..f289279 100644 --- a/include/sampler/LP.h +++ b/include/sampler/LP.h @@ -24,7 +24,7 @@ in P. L'Ecuyer and A. Owen (eds.), Monte Carlo and Quasi-Monte Carlo Methods 2008, Springer-Verlag, 2009. */ -class LarcherPillichshammerGK : public TSamplerMinMaxDim<1, 3> +class LarcherPillichshammerGK : public TSamplerMinMaxDim<1, 1024> { public: LarcherPillichshammerGK(unsigned dimensions = 2, unsigned numSamples = 64, bool randomize = false); diff --git a/include/sampler/MultiJittered.h b/include/sampler/MultiJittered.h index d96ab8d..e3125b3 100644 --- a/include/sampler/MultiJittered.h +++ b/include/sampler/MultiJittered.h @@ -20,7 +20,7 @@ class MultiJittered : public TSamplerDim<2> { public: - MultiJittered(unsigned, unsigned, bool randomize = true, float jitter = 0.0f); + MultiJittered(unsigned x, unsigned y, bool randomize = true, float jitter = 0.0f); ~MultiJittered() override; void clear(); @@ -225,9 +225,7 @@ class CorrelatedMultiJitteredInPlace : public TSamplerMinMaxDim<1, 1024> } float setJitter(float j = 1.0f) override { - m_maxJit = j; - reset(); - return m_maxJit; + return m_maxJit = j; } std::string name() const override @@ -238,9 +236,9 @@ class CorrelatedMultiJitteredInPlace : public TSamplerMinMaxDim<1, 1024> protected: unsigned m_resX, m_resY, m_numSamples, m_numDimensions; float m_maxJit; - bool m_randomize; + bool m_randomize = true; pcg32 m_rand; - unsigned m_seed = 13; - unsigned m_permutation; + unsigned m_seed = 13; + unsigned m_permutation = 13; unsigned m_decorrelate; }; \ No newline at end of file diff --git a/include/sampler/SOA.h b/include/sampler/SOA.h deleted file mode 100644 index 8148dec..0000000 --- a/include/sampler/SOA.h +++ /dev/null @@ -1,78 +0,0 @@ -/** \file SOA.h - \author Afnan Enayet - \brief Strong orthogonal array point sets - - This module generates strong orthogonal arrays from regular orthogonal - arrays. - */ - -#pragma once - -#include - -#include -#include - -/** - A strong orthogonal array sampler class - - The strong orthogonal array uses the construction techniques used in the - regular orthogonal array classes (found in SamplerOrthogonalArray.[ch]), - and converts them to ordered orthogonal arrays, then converting them into - strong orthogonal arrays. - - Currently, SOAs up to strength 5 are supported. - */ -class StrongOrthogonalArray : public TSamplerMinMaxDim<3, 5> -{ -public: - StrongOrthogonalArray(unsigned, unsigned strength, bool randomize = false, - float jitter = 0.0f); - ~StrongOrthogonalArray() override; - void clear(); - - void reset() override; - void sample(float point[], unsigned i) override; - - unsigned dimensions() const override {return m_numDimensions;} - void setDimensions(unsigned d) override {m_numDimensions = d; reset();} - - bool randomized() const override {return m_randomize;} - void setRandomized(bool r) override {m_randomize = r; reset();} - - std::string name() const override; - - int numSamples() const override {return m_numSamples;} - - /** - Set the number of samples to generate points for - - This method will find the smallest prime greater than or equal to `n`, - and then set the number of samples equal to n raised to the desired - strength. - - \param n The base for the number of samples desired, where the number - of samples is \f$ n^{strength} \f$ - */ - int setNumSamples(unsigned n) override; - void setNumSamples(unsigned x, unsigned); - -private: - /** - Construct an OA using the Bush construction - - This constructs an orthogonal array for the Bush construction technique - as described by Art Owen, using the relevant parameters set for the SOA. - This OA will not be randomized or converted back to the canonical point - set. - - \param numDimensions The number of dimensions required for the new - orthogonal array - \returns An orthogonal array - */ - std::vector> bushOA(unsigned numDimensions); - unsigned m_resX, m_numSamples, m_t, m_numDimensions; - bool m_randomize; - std::vector> m_pts; - float m_maxJit; -}; diff --git a/src/app.cpp b/src/app.cpp index a595785..d93a86d 100644 --- a/src/app.cpp +++ b/src/app.cpp @@ -117,7 +117,7 @@ SampleViewer::SampleViewer() m_samplers.emplace_back(new Hammersley(m_num_dimensions, 1)); m_samplers.emplace_back(new Hammersley(m_num_dimensions, 1)); m_samplers.emplace_back(new LarcherPillichshammerGK(3, 1, false)); - m_samplers.emplace_back(new GrayCode(4)); + m_samplers.emplace_back(new GrayCode(1)); m_samplers.emplace_back(new XiSequence(1)); m_samplers.emplace_back(new CSVFile()); @@ -602,7 +602,7 @@ void SampleViewer::draw_editor() if (big_header(ICON_FA_SLIDERS_H " Sampler settings")) // ========================================================= { - if (ImGui::BeginCombo("##Sampler combo", m_samplers[m_sampler]->name().c_str())) + if (ImGui::BeginCombo("##Sampler combo", m_samplers[m_sampler]->name().c_str(), ImGuiComboFlags_HeightLargest)) { for (int n = 0; n < (int)m_samplers.size(); n++) { diff --git a/src/sampler/GrayCode.cpp b/src/sampler/GrayCode.cpp index e289b72..889556e 100644 --- a/src/sampler/GrayCode.cpp +++ b/src/sampler/GrayCode.cpp @@ -4,6 +4,7 @@ #include #include +#include GrayCode::GrayCode(unsigned n) { @@ -17,9 +18,11 @@ void GrayCode::regenerate() m_samples.clear(); m_samples.resize(N); - double res = 1.0 / N; + double offset = m_randomize ? 0. : 0.5; + double res = 1.0 / N; + // Sample in first stratum of each row and column for (unsigned V = 0; V < n; V++) - m_samples[V].x = m_samples[V * n].y = res * (V * n); // Sample in first stratum of each row and column + m_samples[V].x = m_samples[V * n].y = res * (V * n + offset); uint32_t u = 0; // Iterate through strata up the column @@ -39,8 +42,29 @@ void GrayCode::regenerate() // Iterate through columns for (uint32_t V = 0; V < n; V++) - m_samples[U * n + V].x = m_samples[V * n + U].y = res * (V * n + u); + m_samples[U * n + V].x = m_samples[V * n + U].y = res * (V * n + u + offset); } + + if (m_randomize) + { + // perform additional xor scrambling + auto s1 = m_randomize ? m_rand.nextUInt() : 0; + auto s2 = m_randomize ? m_rand.nextUInt() : 0; + for (unsigned i = 0; i < N; ++i) + m_samples[i] = Point{randomDigitScramble(m_samples[i].x, s1), randomDigitScramble(m_samples[i].y, s2)}; + } +} + +int GrayCode::setNumSamples(unsigned num) +{ + int log2N = std::round(std::log2(num)); + if (log2N & 1) + log2N++; // Only allow even powers of 2. + N = 1 << log2N; // Make N a power of 4, if not. + log2n = log2N / 2; + n = 1 << log2n; + regenerate(); + return N; } void GrayCode::sample(float r[], unsigned i) diff --git a/src/sampler/Jittered.cpp b/src/sampler/Jittered.cpp index 937dbd6..201b499 100644 --- a/src/sampler/Jittered.cpp +++ b/src/sampler/Jittered.cpp @@ -3,29 +3,34 @@ */ #include +#include // for permute -Jittered::Jittered(unsigned x, unsigned y, float jitter) : m_resX(x), m_resY(y), m_maxJit(jitter) +Jittered::Jittered(unsigned x, unsigned y, float jitter) : m_maxJit(jitter) { - reset(); -} - -void Jittered::reset() -{ - if (m_resX == 0) - m_resX = 1; - if (m_resY == 0) - m_resY = 1; - m_xScale = 1.0f / m_resX; - m_yScale = 1.0f / m_resY; - // m_rand.seed(m_seed); + setNumSamples(x, y); } void Jittered::sample(float r[], unsigned i) { + i %= m_numSamples; + if (i == 0) m_rand.seed(m_seed); - float jx = 0.5f + int(m_randomize) * m_maxJit * (m_rand.nextFloat() - 0.5f); - float jy = 0.5f + int(m_randomize) * m_maxJit * (m_rand.nextFloat() - 0.5f); - r[0] = ((i % m_resX) + jx) * m_xScale; - r[1] = ((i / m_resY) + jy) * m_yScale; + + for (unsigned d = 0; d < dimensions(); d += 2) + { + int s = permute(i, m_numSamples, m_permutation * 0x51633e2d * (d + 1)); + + // horizontal and vertical indices of the stratum in the jittered grid + int x = s % m_resX; + int y = s / m_resX; + + // jitter in the d and d+1 dimensions + float jx = 0.5f + m_randomize * m_maxJit * (m_rand.nextFloat() - 0.5f); + float jy = 0.5f + m_randomize * m_maxJit * (m_rand.nextFloat() - 0.5f); + + r[d] = (x + jx) * m_xScale; + if (d + 1 < dimensions()) + r[d + 1] = (y + jy) * m_yScale; + } } diff --git a/src/sampler/LP.cpp b/src/sampler/LP.cpp index 1c196c6..e190980 100644 --- a/src/sampler/LP.cpp +++ b/src/sampler/LP.cpp @@ -17,12 +17,17 @@ LarcherPillichshammerGK::~LarcherPillichshammerGK() void LarcherPillichshammerGK::sample(float r[], unsigned i) { - if (m_numDimensions > 0) - r[0] = randomDigitScramble(i * m_inv, m_scramble1); + for (unsigned d = 0; d < dimensions(); d += 3) + { + int s = permute(i, m_numSamples, d); + unsigned ds = 0x68bc21eb * (d + 1); + if (d < dimensions()) + r[d] = randomDigitScramble(s * m_inv, m_scramble1 * ds); - if (m_numDimensions > 1) - r[1] = LarcherPillichshammerRI(i, m_scramble2); + if (d + 1 < dimensions()) + r[d + 1] = LarcherPillichshammerRI(s, m_scramble2 * ds); - if (m_numDimensions > 2) - r[2] = GruenschlossKellerRI(i, m_scramble3); + if (d + 2 < dimensions()) + r[d + 2] = GruenschlossKellerRI(s, m_scramble3 * ds); + } } diff --git a/src/sampler/MultiJittered.cpp b/src/sampler/MultiJittered.cpp index 151ef1c..647d59e 100644 --- a/src/sampler/MultiJittered.cpp +++ b/src/sampler/MultiJittered.cpp @@ -199,8 +199,7 @@ CorrelatedMultiJitteredInPlace::CorrelatedMultiJitteredInPlace(unsigned x, unsig void CorrelatedMultiJitteredInPlace::sample(float r[], unsigned i) { - if (i >= m_numSamples) - i = 0; + i %= m_numSamples; if (i == 0) m_rand.seed(m_seed); diff --git a/src/sampler/SOA.cpp b/src/sampler/SOA.cpp deleted file mode 100644 index 0df83f6..0000000 --- a/src/sampler/SOA.cpp +++ /dev/null @@ -1,126 +0,0 @@ -/** - \file SOA.cpp - \author Afnan Enayet - */ - -#include -#include -#include - -StrongOrthogonalArray::StrongOrthogonalArray(unsigned, unsigned strength, - bool randomize, float jitter) -{ - m_t = strength; - m_randomize = randomize; - m_maxJit = jitter; -} - -StrongOrthogonalArray::~StrongOrthogonalArray() { } - -void StrongOrthogonalArray::clear() -{ - m_pts.clear(); -} - -void StrongOrthogonalArray::reset() -{ - // set the number of dimensions to have a max of p + 1 - m_numDimensions = std::min(m_numDimensions, m_resX + 1); - - // initialize m_pts for the number of samples and dimensions - m_pts.resize(m_numSamples); - - for (auto &point : m_pts) { - point = std::vector(); - point.resize(m_numDimensions); - } - - // first, we need to determine the new m' value - // this equation is given by He and Tang (2012) - // we derive it by solving for m, instead of m' - // our m' is the number of dimensions, and m is - // the dimensionality of the original orthogonal array - unsigned m_prime = m_numDimensions; - - // m gives us the dimensionality of the OA we need to generate - unsigned m; - - if (m_t % 2 == 0) - m = unsigned(std::floor(m_t * m_prime / 2)); - else - m = unsigned(std::floor((m_prime * m_t - 2) / 2)); - - auto oa = bushOA(m); - - // construct the column selection vector that will generate the ordered - // orthogonal array - // this column will yield which column of the original OA corresponds to - // the same column in the ordered orthogonal array - std::vector ooaCols(m_prime * m_t); - - // TODO set up the column selection vector for each strength, using He and Tang's method - switch(m_t) { - case 2: - // determine - break; - case 3: - break; - case 4: - break; - case 5: - break; - } -} - -void StrongOrthogonalArray::sample(float point[], unsigned i) -{ - if (i > m_numSamples) - i = 0; - - for (unsigned d = 0; d < m_numDimensions; d++) { - point[d] = m_pts[i][d]; - } -} - -std::string StrongOrthogonalArray::name() const -{ - return "Strong Orthogonal Array"; -} - -int StrongOrthogonalArray::setNumSamples(unsigned n) -{ - // We want to set the number of samples closest to some power of a prime - int rtVal = (n == 0) ? 1 : (int) (std::pow((float) n, 1.f/m_t) + 0.5f); - setNumSamples(rtVal, rtVal); - return m_numSamples; -} - -void StrongOrthogonalArray::setNumSamples(unsigned x, unsigned) -{ - m_resX = primeGE(x); - m_numSamples = pow(m_resX, m_t); - reset(); -} - -std::vector> StrongOrthogonalArray::bushOA(unsigned numDimensions) { - std::vector> oa; - oa.resize(m_numSamples); - - for (unsigned i = 0; i < m_numSamples; i++) { - oa[i].resize(numDimensions); - auto coeffs = iToPolyCoeffs(i, m_resX, m_t); - - // The Bush construction technique uses a particular construction for - // columns 0...p, and another procedure for p + 1, so that will be - // handled after the loop, if necessary - auto upperDim = std::min(numDimensions, m_resX); - - for (unsigned j = 0; j < upperDim; j++) { - oa[i][j] = polyEval(coeffs, j) % m_resX; - } - - if (numDimensions >= m_resX + 1) - oa[i][m_resX] = i % m_resX; - } - return oa; -}