Page MenuHome
Paste P1092

(An Untitled Masterwork)
ActivePublic

Authored by Brecht Van Lommel (brecht) on Sep 3 2019, 6:56 PM.
/*
* Sobol sequence direction vectors.
*
* This file contains code to create direction vectors for generating sobol
* sequences in high dimensions. It is adapted from code on this webpage:
*
* http://web.maths.unsw.edu.au/~fkuo/sobol/
*
* From these papers:
*
* S. Joe and F. Y. Kuo, Remark on Algorithm 659: Implementing Sobol's quasirandom
* sequence generator, ACM Trans. Math. Softw. 29, 49-57 (2003)
*
* S. Joe and F. Y. Kuo, Constructing Sobol sequences with better two-dimensional
* projections, SIAM J. Sci. Comput. 30, 2635-2654 (2008)
*/
/* Copyright (c) 2008, Frances Y. Kuo and Stephen Joe
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* * Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
*
* * Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* * Neither the names of the copyright holders nor the names of the
* University of New South Wales and the University of Waikato
* and its contributors may be used to endorse or promote products derived
* from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY
* EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY
* DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
CCL_NAMESPACE_BEGIN
#define SOBOL_MAX_NUMBER 32
#define SOBOL_BITS 32
#define SOBOL_MAX_DIMENSIONS 10
/* Skip initial numbers that for some dimensions have clear patterns that
* don't cover the entire sample space. Ideally we would have a better
* progressive pattern that doesn't suffer from this problem, because even
* with this offset some dimensions are quite poor.
*/
#define SOBOL_SKIP 64
typedef struct SobolDirectionNumbers {
uint d, s, a;
uint m[SOBOL_MAX_NUMBER];
} SobolDirectionNumbers;
/* Note: this file is skipped by clang-format. */
/* Keep simple alignment. */
/* clang-format off */
static const SobolDirectionNumbers SOBOL_NUMBERS[SOBOL_MAX_DIMENSIONS - 1] = {
{2, 1, 0, {1}},
{3, 2, 1, {1, 3}},
{4, 3, 1, {1, 3, 1}},
{5, 3, 2, {1, 1, 1}},
{6, 4, 1, {1, 1, 3, 3}},
{7, 4, 4, {1, 3, 5, 13}},
{8, 5, 2, {1, 1, 5, 5, 17}},
{9, 5, 4, {1, 1, 5, 5, 5}},
{10, 5, 7, {1, 1, 7, 11, 19}},
{11, 5, 11, {1, 1, 5, 1, 1}},
};
/* clang-format on */
void sobol_generate_direction_vectors(uint vectors[][SOBOL_BITS], int dimensions)
{
assert(dimensions <= SOBOL_MAX_DIMENSIONS);
const uint L = SOBOL_BITS;
/* first dimension is exception */
uint *v = vectors[0];
for (uint i = 0; i < L; i++)
v[i] = 1 << (31 - i); // all m's = 1
for (int dim = 1; dim < dimensions; dim++) {
const SobolDirectionNumbers *numbers = &SOBOL_NUMBERS[dim - 1];
const uint s = numbers->s;
const uint a = numbers->a;
const uint *m = numbers->m;
v = vectors[dim];
if (L <= s) {
for (uint i = 0; i < L; i++)
v[i] = m[i] << (31 - i);
}
else {
for (uint i = 0; i < s; i++)
v[i] = m[i] << (31 - i);
for (uint i = s; i < L; i++) {
v[i] = v[i - s] ^ (v[i - s] >> s);
for (uint k = 1; k < s; k++)
v[i] ^= (((a >> (s - 1 - k)) & 1) * v[i - k]);
}
}
}
}
float sobol_sample(uint vectors[][SOBOL_BITS], int index, int dimension)
{
uint result = 0;
uint i = index + SOBOL_SKIP;
for (uint j = 0; i; i >>= 1, j++) {
if (i & 1) {
result ^= vectors[dimension][j];
}
}
return result * (1.0f / (float)0xFFFFFFFF);
}
CCL_NAMESPACE_END

Event Timeline