UT_SobolSequence.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * NAME:        UT_SobolSequence.h ( UT Library, C++)
 *
 * COMMENTS:    A sobol low-discrepancy sequence.
 */

#ifndef __UT_SobolSequence__
#define __UT_SobolSequence__

#include "UT_API.h"
#include "UT_Array.h"
#include "UT_Assert.h"
#include "UT_StackAlloc.h"
#include "UT_Vector3.h"
#include <SYS/SYS_Math.h>
#include <SYS/SYS_StaticAssert.h>
#include <SYS/SYS_Types.h>
#include <type_traits>

/// Random Sobol sequence of dimension _DIM_ with _BITS_ bits of integral
/// *unsigned* type _itype_
template <class itype, uint DIM, uint BITS = sizeof(itype)*8>
class UT_SobolSequence
{
public:
                    UT_SobolSequence();

    constexpr uint  dimension() const   { return DIM; }
    constexpr uint  bits() const        { return BITS; }

    /// Advance to the next sample, storing the result into x.  When myN
    /// exceeds the range of the sequence, it will automatically wrap.
    void         advance();

    /// Return the sequence value for dimension dim.
    itype        x(exint dim) const     { return myX[dim]; }
    fpreal32     xf(exint dim) const
                 {
                     SYS_FPRealUnionF   fval;

                     // Move significant digit to line up with mantissa and
                     // bitwise OR it with 1.0 to get [1.0,2.0).
                     // BITS_FROM_23 is moved out to avoid compiler warning.
                     if (BITS < 23)
                        fval.uval = 0x3f800000 | (uint)(x(dim) << BITS_FROM_23);
                     else
                        fval.uval = 0x3f800000 | (uint)(x(dim) >> BITS_FROM_23);
                     return fval.fval - 1;
                 }

    /// Reset the sequence
    void         rewind(itype n = 0);

private:
    /// Returns an odd random value between 0 and val.
    static constexpr uint
    orand(uint val, uint &seed)
    {
        return ((SYSfastRandomInt(seed) % (val >> 1)) << 1) | 1;
    }
    static constexpr itype
    computeMask()
    {
        itype mask = 0;
        for (uint i = 0; i < BITS; i++)
            mask |= itype(1) << i;
        return mask;
    }

private:
    static constexpr uint   LEN_IV = DIM * (BITS + 1);
    static constexpr uint   LEN_X = DIM;
    static constexpr uint   MAXDIM = 6;
    static constexpr itype  MASK = computeMask();
    static constexpr uint   BITS_FROM_23 = (BITS<23) ? (23-BITS) : (BITS-23);

    UT_Array<itype>         myIV;
    UT_Array<itype>         myX;
    itype                   myN;
};

class UT_VectorSequence : UT_SobolSequence<uint64,3>
{
    using                PARENT = UT_SobolSequence<uint64,3>;

public:
                         UT_VectorSequence() = default;

    using                PARENT::advance;
    using                PARENT::rewind;

    UT_Vector3           getVector()
                         {
                             return UT_Vector3(PARENT::xf(0),
                                               PARENT::xf(1),
                                               PARENT::xf(2));
                         }
};


//////////////////////////////////////////////////////////////////////////////
//
// Implementation
//

template <class itype, uint DIM, uint BITS>
UT_SobolSequence<itype,DIM,BITS>::UT_SobolSequence()
    : myIV(LEN_IV, LEN_IV)
    , myX(LEN_X, LEN_X)
    , myN(0)
{
    constexpr uint    mdeg[MAXDIM] = {1,2,3,3,4,4};
    constexpr uint    ip[MAXDIM] = {0,1,1,2,1,4};

    SYS_STATIC_ASSERT(std::is_unsigned<itype>::value);
    SYS_STATIC_ASSERT(DIM >= 1 && DIM <= MAXDIM);
    SYS_STATIC_ASSERT(BITS >= 1 && BITS <= sizeof(itype)*8);

    // Indexing
    itype **iu = (itype **)UTstackAlloc((BITS+1)*sizeof(itype *));
    for (uint i = 0; i < BITS+1; i++)
        iu[i] = &myIV[i*DIM];

    // Fill in the initial odd random values
    uint seed = 1337;
    for (uint j = 0; j < DIM; j++)
        for (uint i = 0; i < mdeg[j]; i++)
            iu[i][j] = orand(1 << (i+1), seed);

    for (uint k = 0; k < DIM; k++)
    {
        // Initial values only require normalization
        for (uint j = 0; j < mdeg[k]; j++)
            iu[j][k] <<= (BITS-j-1);

        // Calculate the remaining values using the recurrence equation
        for (uint j = mdeg[k]; j < BITS+1; j++)
        {
            itype ipp = ip[k];

            // First and last terms
            itype tmp = iu[j-mdeg[k]][k];
            tmp ^= (tmp >> mdeg[k]);

            // Intermediate terms
            for (uint l = mdeg[k]-1; l >= 1; l--)
            {
                if (ipp & 1)
                    tmp ^= iu[j-l][k];
                ipp >>= 1;
            }
            iu[j][k] = tmp;
        }
    }

    rewind();

    UTstackFree(iu);
}

template <class itype, uint DIM, uint BITS>
void
UT_SobolSequence<itype,DIM,BITS>::advance()
{
    if (myN < MASK)
    {
        // Find the first 0-bit in n
        uint i;
        itype n = myN & MASK;
        for (i = 0; i < BITS; i++)
        {
            if (!(n & 1))
                break;
            n >>= 1;
        }

        // XOR the bits into the sequence entry
        n = i*DIM;
        for (i = 0; i < DIM; i++)
        {
            myX[i] ^= myIV[n+i];
        }
        myN++;
    }
    else
    {
        rewind();
    }
}

template <class itype, uint DIM, uint BITS>
void
UT_SobolSequence<itype,DIM,BITS>::rewind(itype n)
{
    n &= MASK;
    myN = n;
    for (uint i = 0; i < DIM; i++)
        myX[i] = 0;

    uint off = 0;
    while (n)
    {
        if ((n+1) & 0x2)
        {
            for (uint i = 0; i < DIM; i++)
                myX[i] ^= myIV[off+i];
        }
        n >>= 1;
        off += DIM;
    }
}

#endif // __UT_SobolSequence__