- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I have tried to initialize Sobol sequence with my own polynomial and direction number as it is described in the vslNote.
First of all there is a lack of infornation in the documentation as we do not know the convention used to enter polynomial and direction( it is said that they must be 32 bit array, indeed there are many possibility: for polynomial the first coeficient is alwys one so many implementation of sobol take just the other coeficient...)
Anyway I didn't succes just to use my initialisation for the first 40 dimention, whatever I use as initial direction I got the same sequence with the following code.
paramToFill[0] = fGeneratorSize;
paramToFill[1] = VSL_USER_QRNG_INITIAL_VALUES;
//here I use a mix of the VSL note and the header definition. Anyway none of the other possibilityies work
paramToFill[2] = VSL_USER_PRIMITIVE_POLYMS | VSL_USER_INIT_DIRECTION_NUMBERS | VSL_USER_DIRECTION_NUMBERS;
if(overRideFirstDim) paramToFill[2] = paramToFill[2] | VSL_QRNG_OVERRIDE_1ST_DIM_INIT;
paramToFill[3+fGeneratorSize] = maxDegree;
for (size_t compteur = 0; compteur<fGeneratorSize; ++compteur)
{
SSSobolDimension& refSobolDim = sobolDimension[compteur];
size_t upperToFill = refSobolDim.degree;
if(tryComplete)
{
upperToFill = maxDegree;
refSobolDim.CompleteInitialisation(maxDegree);
}
paramToFill[3+compteur] = refSobolDim.polynom;
for (size_t idIndex = 0; idIndex<upperToFill; ++idIndex)
{
paramToFill[3+maxDegree*compteur+idIndex] = refSobolDim.initialDirection[idIndex];
}
}
VSLStreamStatePtr * streamGenerator = (VSLStreamStatePtr *) fGeneratorVector;
int status = 0;
for (size_t j=0; j<fThreadCount; ++j )
status = vslNewStreamEx( &streamGenerator[j], VSL_BRNG_SOBOL, fGeneratorSize,paramToFill));
Did I miss something?
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I found some misunderstanding in my code:
now I am using
vslNewStreamEx
( &streamGenerator[j], VSL_BRNG_SOBOL, parameterSize, paramToFill);where parameterSize is the size of paramToFill.
Anyways if I am using less than or exactly 40dimension my initialisation is not taken into account and the sequences generated are the 40 dimension sequence of Bradley&Fox. If I try 41 or more dimention then myt initilisation is not taken into account but now only one dimension is generated and when I am calling
vdRngUniform
( VSL_METHOD_SUNIFORM_STD_ACCURATE, streamGenerator[threadID-1], fGeneratorSize, fLocalBuffer.GetReference(), 0.0, 1.0 );I get the sequence generated by using a congruence modulo fGeneratorSize on the first sequence which is very very bad in case of sobol...
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thank you for your post and feedback to the VSL documentation.
Ideally, in order to understand why things go wrong during Sobol initialization on your side it could be nice to have working short test case which would demonstrate the issue.
To make initialization of VSL Sobol QRNG you need to pass into the library either polynomials and initial direction numbers or table of direction numbers computed on you side. Initialization of the Sobol QRNG is done using vslNewStreamEx function:
errcode = vslNewStreamEx( &stream, brng, nparams, params );
Prior call to the function nparams and params should be properly initialized. Array of unsigned integers params of size nparams is used to pack Sobol initialization values.
First element is maximal dimension of quasi-random vector you are planning to have:
params[0] = dimen;
The next field defines that QRNG parameters are passed into the library. It should be set to VSL_USER_QRNG_INITIAL_VALUES:
params[1] = VSL_USER_QRNG_INITIAL_VALUES;
Value of the next field determines what you are going to pass in the library, polynomials/initial direction numbers or table of Sobol direction numbers.
Let's set it to VSL_USER_PRIMITIVE_POLYMS | VSL_USER_INIT_DIRECTION_NUMBERS if there is a need to pass both polynomials & initial direction numbers:
params[2] = VSL_USER_INIT_DIRECTION_NUMBERS | VSL_USER_PRIMITIVE_POLYMS;
The next dimen positions are used to pack polynomials:
for ( i = 0; i < dimen; i++ ) params[3+i] = uSobolPoly;
The polynomials are packed in the following way: least significant bit contains coefficient at zero degree, the next bit contains coefficient at 1st degree and so forth.
The library needs to know about the maximal dimension of polynomials (over all polynomials used to initialize Sobol):
params[3+dimen] = maxdeg;
The rest positions are used to hold initial direction numbers:
k = 4+dim;
for ( i = 0; i < dimen-1; i++ )
{
for ( j = 0; j < maxdeg; j++ )
{
params[k++] = uSobolMInit
}
}
where 2 dimensional table uSobolMInit of size num_of_poly*maxdeg contains initial direction numbers. If you do not want to override default initialization for 1st dimension num_of_poly equals to dimen-1 (otherwise, it equals to dimen and params[2] = params[2] | VSL_QRNG_OVERRIDE_1ST_DIM_INIT).
Thus, size nparams of array params should be enough to hold at least 3 + num_of_poly + 1 + num_of_poly * maxdeg:
nparams = 3 + num_of_poly + 1 + num_of_poly * maxdeg;
Value of parameter nparams depends on type of initialization and, thus, number of Sobol parameters to be passed in VSL.
Finally, vslNewStreamEx routine is called to initialize Sobol QRNG.
Please, let me know if you have more questions/concerns on this and, also, how this scheme works for you. Also, please fill free to attach a test case with your initialization of VSL Sobol.
Also, what version of MKL do you use?
Andrey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thank for your explanation.
As far as the convention for polynimial, ifI understand correctly x^3+x+1 is 1011 = 1*8+0*4+1*2+1=11.
For the initial direction it is not clear. Let us says that my first two dimension are 0.5 and 0.75 it can be coded in an int as 1 and 3 or (1 << 31) and (2 << 30). What is the convention?
If I understand correctly the flag VSL_QRNG_OVERRIDE_1ST_DIM_INIT is to say if I overwrite the first dimension. But in any case I can fill dimen polynomial and dimen*maxdeg initial direction. If I do not overwrite I just give 1 dimension useless data. I am right?
I am using MKL 10.0.
Thanks for you help.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I finally manage to initialize sequence:
the initial direction convention is 0.5 and 0.75must be coded as 1 and 3.
when using VSL_QRNG_OVERRIDE_1ST_DIM_INIT the initialisation data are exactly as explained in the vsl note. If you do not use then for a sequence of dimention dimen you must use dimen in parameter 0 but for the positions in vector param you must replace dimen by dimen-1 everywhere.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
As far as the convention for polynimial, if I understand correctly x^3+x+1 is 1011 = 1*8+0*4+1*2+1=11.
Yes, this is correct.
For the initial direction it is not clear. Let us says that my first two dimension are 0.5 and 0.75 it can be coded in an int as 1 and 3 or (1 << 31) and (2 << 30). What is the convention?
First variant, that is 1 and 3, should be used for initial direction numbers.
If I understand correctly the flag VSL_QRNG_OVERRIDE_1ST_DIM_INIT is to say if I overwrite the first dimension.
Yes, that's right.
But in any case I can fill dimen polynomial and dimen*maxdeg initial direction. If I do not overwrite I just give 1 dimension useless data. I am right?
Let us consider the example. You want to pass initial parameters for dimen=40 and do not want to override 1st dimension.In this case you need to pass only 40-1=39 sets of polynomials/initial direction numbers.If you want to override 1st dimension you need to pass 40 sets of polynomials/initial direction numbers.Value of maxdeg should be properly defined in each case. In other words, the library expects only "useful" data. Passing 1st dimension related data into the library when flag VSL_QRNG_OVERRIDE_1ST_DIM_INIT is not set would result in incorrect generation of quasi-random sequences.
Please, let me know how generation of quasi-random vectors with your initialization parameters work. Please, feel free to address any questions/concerns on this. Also, you may want to have a look into soboluserdirnums.c file which is available in examplesvslcsource directory of your MKL installation package.
Thanks,
Andrey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Thanks for the explanation. I will look in the .c files.
In fact my initialisation work only when I am asking for less than 40 dimensions simultaneously. If I am asking for let says 41 dimension then all dimension return the first dimension sobol sequences. I do not know if it is an internal limitation or if it is a remaining mistake in my code.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Can you please provide short test case which would demonstrate the issue and/or output of the generator with parameters of the generation? Thanks, Andrey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
As it is difficult to provide my original code (cpp with many object and method) I have modified the code of the example provided by the MKL: I have added one polynom in the list and the corresponding dimention. Then I initialize a 40 dimensions generator and if you look in the array named r after the initialisation you will see that all the 40 sequences are differents. Then I do the same initialisation with dimension 41 and now if you look in array r you will see that only the first and second dimension are differents, all the others are the same as the first one.
Did I make a mistake somewhere?
//the code:
/* Dimension of quasi-random vectors */
#define
DIM 41/* Number of bits used for storing unsigned integer */
#define
NBITS 32/* Number of primitive polynomials */
#define
NPOLY 40/* Maximum degree of primitive polynomials */
#define
MAXDEG 8/* Number of digits to scramble */
#define
S_DIGITS 32/* Table of initial direction numbers */
static
unsigned int uSobolMInit[NPOLY][MAXDEG] = {1,0,0, 0, 0, 0, 0, 0,
1,1,0, 0, 0, 0, 0, 0,
1,3,7, 0, 0, 0, 0, 0,
1,1,5, 0, 0, 0, 0, 0,
1,3,1, 1, 0, 0, 0, 0,
1,1,3, 7, 0, 0, 0, 0,
1,3,3, 9, 9, 0, 0, 0,
1,3,7,13, 3, 0, 0, 0,
1,1,5,11,27, 0, 0, 0,
1,3,5, 1,15, 0, 0, 0,
1,1,7, 3,29, 0, 0, 0,
1,3,7, 7,21, 0, 0, 0,
1,1,1, 9,23,37, 0, 0,
1,3,3, 5,19,33, 0, 0,
1,1,3,13,11, 7, 0, 0,
1,1,7,13,25, 5, 0, 0,
1,3,5,11, 7,11, 0, 0,
1,1,1, 3,13,39, 0, 0,
1,3,1,15,17,63, 13, 0,
1,1,5, 5, 1,27, 33, 0,
1,3,3, 3,25,17,115, 0,
1,1,3,15,29,15, 41, 0,
1,3,1, 7, 3,23, 79, 0,
1,3,7, 9,31,29, 17, 0,
1,1,5,13,11, 3, 29, 0,
1,3,1, 9, 5,21,119, 0,
1,1,3, 1,23,13, 75, 0,
1,3,3,11,27,31, 73, 0,
1,1,7, 7,19,25,105, 0,
1,3,5, 5,21, 9, 7, 0,
1,1,1,15, 5,49, 59, 0,
1,1,1, 1, 1,33, 65, 0,
1,3,5,15,17,19, 21, 0,
1,1,7,11,13,29, 3, 0,
1,3,7, 5, 7,11,113, 0,
1,1,5, 3,15,19, 61, 0,
1,3,1, 1, 9,27, 89, 7,
1,1,3, 7,31,15, 45,23,
1,3,3, 9, 9,25,107,39,
1,3,1,11,11,11,77,249
};
/* Array of primitive polynom ials */
static
unsigned int uSobolPoly[NPOLY] = {3, 7, 11, 13, 19, 25, 37, 59, 47,
61, 55, 41, 67, 97, 91, 109, 103, 115, 131,
193, 137, 145, 143, 241, 157, 185, 167, 229, 171,
213, 191, 253, 203, 211, 239, 247, 285, 369, 299,
301
};
/* Array of degrees of primitive polynomials */
static
int iSobolPolyDeg[NPOLY] = {1, 2, 3, 3, 4, 4, 5, 5, 5,
5, 5, 5, 6, 6, 6, 6, 6, 6, 7,
7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
7, 7, 7, 7, 7, 7, 7, 8, 8, 8,
8
};
/* Generator matrix, table of direction numbers */
static
unsigned int uSobolV[DIM][NBITS];/* Generator matrix, table of modified direction numbers */
static
unsigned int uSobolSV[DIM][NBITS];/* Routine to initialize matrix of Sobol */
static
void InitSobolMatrix( int dimen );/* Routine to initialize Sobol matrix for Owen-like scrambled sequence */
static
void InitSobolMatrixOs( int dimen );/* Routine to generate sequence of low triangular matrices for
Owen scrambling */
static
void GenerateLSM( int dimen, unsigned int lsm[ DIM][S_DIGITS] );/* Dimension of vectors to generate */
/* Size of params array */
#define
NU DIM*NBITS + 4#define
MaxDim 41#define
NN 10/* Size of buffer for quasi-random sequence */
#define
N 1000*MaxDim/* Buffer for quasi-random sequence */
double
r[N];double
rt[N];int
testGeneration()//main(){
VSLStreamStatePtr stream; MKL_INT brng; int dim; int errcode; double a, b; int i, j, k; MKL_INT nparams; int errnum; /***** Initialize *****/ brng = VSL_BRNG_SOBOL; a = 0.0; b = 1.0; /* dim=40 work fine */ dim = 40; size_t{
params[3+i] = uSobolPoly[i]; for ( j = 0; j < maxDegree; j++ ) params[4+sizeToUse+maxDegree*i+j] = uSobolMInit[i][j< /FONT>];}
/***** Initialize ****/ errcode = vslNewStreamEx( &stream, brng, parameterSize, params ); //CheckVslError( errcode ); /***** Call RNG *****/ errcode = vdRngUniform( VSL_METHOD_SUNIFORM_STD_ACCURATE, stream, 1000*dim, r, a, b ); //CheckVslError( errcode ); delete [] params; /***** Printing results *****/ printf("Sobol sequence with user-defined direction numbers: "); printf("-------------------------------------------------- "); printf("Parameters: "); printf(" a=%.4f ",a); printf(" b=%.4f ",b); k = 12; printf("Results: %d component of quasi-random sequence (first %d of %d): ", k, NN, N/dim); for(i=0;i<NN;i++) { printf("%.4f ",r[k+i*dim]);}
printf(" "); /***** Deinitialize *****/ errcode = vslDeleteStream( &stream ); //CheckVslError( errcode ); /* dim=41 does not work */ dim = 41; sizeToUse = dim-1;// must be dim if we overwrite 1st dim maxDegree = iSobolPolyDeg[sizeToUse-1]; parameterSize = 4+sizeToUse*(1+maxDegree); params = new unsigned int[parameterSize]; params[0] = dim; params[1] = VSL_USER_QRNG_INITIAL_VALUES; params[2] = VSL_USER_PRIMITIVE_POLYMS | VSL_USER_INIT_DIRECTION_NUMBERS; params[3+sizeToUse] = maxDegree; k = 3; for ( i = 0; i < sizeToUse; i++ ){
params[3+i] = uSobolPoly[i]; for ( j = 0; j < maxDegree; j++ ) params[4+sizeToUse+maxDegree*i+j] = uSobolMInit[i][j];}
/***** Initialize ****/ errcode = vslNewStreamEx( &stream, brng, parameterSize, params ); //CheckVslError( errcode ); /***** Call RNG *****/ errcode = vdRngUniform( VSL_METHOD_SUNIFORM_STD_ACCURATE, stream, 1000*dim, r, a, b ); //CheckVslError( errcode ); delete [] params; /***** Printing results *****/ printf("Sobol sequence with user-defined direction numbers: "); printf("-------------------------------------------------- "); printf("Parameters: "); printf(" a=%.4f ",a); printf(" b=%.4f ",b); k = 12; printf("Results: %d component of quasi-random sequence (first %d of %d): ", k, NN, N/dim); for(i=0; i<NN;i++) { printf("%.4f ",r[k+i*dim]);}
printf(" "); /***** Deinitialize *****/ errcode = vslDeleteStream( &stream ); //CheckVslError( errcode ); return 0;}
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Just a remark. I am using the debugger to look in the array. If you look the result (the 12th dimension) then you must see that it is not the same when I am usig 40 and when I am using 41.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
OK I will report the problem tomorrow. I have checked that the workarround allow to define more than 40 simultaneous dimensions.
Thanks for your help.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
"I have tried to initialize Sobol sequence .."
Initialization of the Sobol' sequence or developingthe direction numbers is not easy and may take many months-men work. There are known well tested high-dimensional Sobol' sequence generators.
Currently, the Sobol' sequence generator with highest maximum dimension isSobolSeq6144. It wasdeveloped byBRODA ( www.broda.co.uk ) in collaboration with I. M. Sobol'. Not only this generator has very high dimensionality and employs the super fast generation algorithm but generated Sobol' sequences satisfy homogeneity Property A in all dimensions.
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page