Sobol is a quasi-random number generator which produces multi-dimensional low discrepancy sequences. For s-dimensional space, it generates sequences in unit cube \(\ \ \ I^s = [0,1]^s\). Grey code implementation for \(s^{th}\) dimension has form:
where \(v_i\) is direction number, \(c_n\) is the right most 0th position in binary representation of \(n\), \(\ x_0 = 0,\ x_1 = v_{c_0},\ i.e \ x_1 = v_1\). This is used to recursively generate sobol sequence in each dimension.
For specific dimension, direction numbers \(v_i\) are computed using set of initialization numbers \(m_i\) and primitive polynomials. Primitive Polynomials are unique for each dimension.
Primitive polynomials P of degree d are in form of :
where \(a_i\) is either 0 or 1.
For d-degree polynomial, initialization numbers {\(m_1, m_2 \ .... \ m_d\)} are chosen such that \(m_i \ is \ odd \ \) and \(\ m_i < 2^i\), then next set of initialization numbers are computed using following recursive relation:
Directions numbers are computed using equation:
[Sobol_Link] provide implementation with set of polynomials, initialization numbers for 40 dimensions, which are used to generate sobol sequences upto length \(2^{30}\). [Sobol]
Refer code snippet provided below to use Sobol base generator.
/* Generate sobol sequences of dimesions 40 and length 100 */
#define SOBOL_GENID 7
#define SOBOL_STATE_LEN 51
void generate_sobol_sequence (int dimension,
int sequence_length,
double *outbuf, int outbuf_len)
{
rng_int_t genid, subid, ldimension, lstate, info;
rng_int_t dimension[ldimension], state[lstate];
double a, b;
/* Check */
if (outbuf_len != dimension * sequence_length) {
printf("Out buffer length is not sufficient !!!\n");
return;
}
/* Use sobol Generator as the base generator */
genid = SOBOL_GENID;
subid = 0;
ldimension = 1;
lstate = SOBOL_STATE_LEN;
dimension[0] = dimension;
/* Initialize sobol as base generator */
drandinitialize(genid, subid, dimension,
&ldimension, state, &lstate, &info);
/* Generate a sequence from a uniform U(0,1) distribution */
a = 0.0;
b = 1.0;
dranduniform((dimension * sequence_length),
a, b, state, outbuf, &info);
/*
* Output sequence is expected to be in following form
* outbuf[0], outbuf[1], outbuf[2] ... outbuf[dim-1]
* first member of the sequence;
* outbuf[dim], outbuf[dim+1], outbuf[dim+2] ... outbuf[2*dim-1]
* second member of the sequence continues like this till given
* length of sequence is generated.
*/
/* Print the sequence */
printf("Sobol Sequence: \n");
for (int i=0; i<sequence_length; i++) {
for (int j=0; j<dimension; j++) {
printf("%10.4f", outbuf[(i*dim)+ j] );
}
if ((i+1)%4 == 0)
printf("\n");
}
}