BrookRandomNumberGenerator.cpp 41.5 KB
Newer Older
Mark Friedrichs's avatar
Mark Friedrichs committed
1
2
3
4
5
6
7
8
/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
9
10
 * Portions copyright (c) 2009 Stanford University and the Authors.           *
 * Authors: Mark Friedrichs, Mike Houston                                     *
Mark Friedrichs's avatar
Mark Friedrichs committed
11
12
 * Contributors:                                                              *
 *                                                                            *
13
14
15
16
 * This program is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU Lesser General Public License as published   *
 * by the Free Software Foundation, either version 3 of the License, or       *
 * (at your option) any later version.                                        *
Mark Friedrichs's avatar
Mark Friedrichs committed
17
 *                                                                            *
18
19
20
21
 * This program is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
 * GNU Lesser General Public License for more details.                        *
Mark Friedrichs's avatar
Mark Friedrichs committed
22
 *                                                                            *
23
24
 * You should have received a copy of the GNU Lesser General Public License   *
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
Mark Friedrichs's avatar
Mark Friedrichs committed
25
26
27
28
 * -------------------------------------------------------------------------- */

#include <sstream>
#include "BrookRandomNumberGenerator.h"
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
29
#include "../../reference/src/SimTKUtilities/SimTKOpenMMUtilities.h"
30
#include "openmm/OpenMMException.h"
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
31
#include "kernels/kupdatesd.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

using namespace OpenMM;
using namespace std;

/** 
 *
 * Constructor
 * 
 */

BrookRandomNumberGenerator::BrookRandomNumberGenerator( ){

// ---------------------------------------------------------------------------------------

   //static const std::string methodName      = "BrookRandomNumberGenerator::BrookRandomNumberGenerator";

// ---------------------------------------------------------------------------------------

   // fixed for now

   _numberOfRandomNumberStreams     = 2;
Mark Friedrichs's avatar
Mark Friedrichs committed
53
   _randomNumberGeneratorStreams    = NULL;
Mark Friedrichs's avatar
Mark Friedrichs committed
54
55
56
57
58
59
60

   // mark stream dimension variables as unset

   _randomNumberStreamWidth         = -1;
   _randomNumberStreamHeight        = -1;
   _randomNumberStreamSize          = -1;

Mark Friedrichs's avatar
Mark Friedrichs committed
61
62
63
   _rvStreamIndex                   = 0;
   _rvStreamOffset                  = 0;
   _numberOfShuffles                = 0;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
64
65
   //_maxShuffles                     = 0;
   _maxShuffles                     = 100;
Mark Friedrichs's avatar
Mark Friedrichs committed
66
67
68
69

   _loadBuffer                      = NULL;
   _shuffleIndices                  = NULL;

Mark Friedrichs's avatar
Mark Friedrichs committed
70
   for( int ii = 0; ii < LastStreamIndex; ii++ ){
Mark Friedrichs's avatar
Mark Friedrichs committed
71
      _auxiliaryStreams[ii]   = NULL;
Mark Friedrichs's avatar
Mark Friedrichs committed
72
73
   }

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
74
   // set randomNumber seed & generator
Mark Friedrichs's avatar
Mark Friedrichs committed
75

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
76
   _randomNumberSeed      = 1393;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
77
   //_randomNumberGenerator = Mersenne;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
78
   _randomNumberGenerator = Kiss;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
79
   //_randomNumberGenerator = FixedValue;
Mark Friedrichs's avatar
Mark Friedrichs committed
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

   //SimTKOpenMMUtilities::setRandomNumberSeed( randomNumberSeed );
}   
 
/** 
 * Destructor
 * 
 */

BrookRandomNumberGenerator::~BrookRandomNumberGenerator( ){

// ---------------------------------------------------------------------------------------

   //static const std::string methodName      = "BrookRandomNumberGenerator::~BrookRandomNumberGenerator";

// ---------------------------------------------------------------------------------------

   for( int ii = 0; ii < LastStreamIndex; ii++ ){
Mark Friedrichs's avatar
Mark Friedrichs committed
98
      delete _auxiliaryStreams[ii];
Mark Friedrichs's avatar
Mark Friedrichs committed
99
100
   }

Mark Friedrichs's avatar
Mark Friedrichs committed
101
102
103
104
   delete[] _randomNumberGeneratorStreams;

   delete[] _loadBuffer;
   delete[] _shuffleIndices;
Mark Friedrichs's avatar
Mark Friedrichs committed
105
106
107
108
109
110
111
112
113
114
115
116
117
118
}

/** 
 * Get number of random number streams
 * 
 * @return     number of random number streams 
 *
 */

int BrookRandomNumberGenerator::getNumberOfRandomNumberStreams( void ) const {
   return _numberOfRandomNumberStreams;
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
119
 * Get random number stream width
Mark Friedrichs's avatar
Mark Friedrichs committed
120
 * 
Mark Friedrichs's avatar
Mark Friedrichs committed
121
 * @return     nndom number stream width
Mark Friedrichs's avatar
Mark Friedrichs committed
122
123
124
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
125
126
int BrookRandomNumberGenerator::getRandomNumberStreamWidth( void ) const {
   return _randomNumberStreamWidth;
Mark Friedrichs's avatar
Mark Friedrichs committed
127
}
Mark Friedrichs's avatar
Mark Friedrichs committed
128
129
130
131
132
133
134
135
136
137
138
139

/** 
 * Get random number stream height
 * 
 * @return     nndom number stream height
 *
 */

int BrookRandomNumberGenerator::getRandomNumberStreamHeight( void ) const {
   return _randomNumberStreamHeight;
}

Mark Friedrichs's avatar
Mark Friedrichs committed
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
/** 
 * Get random number seed
 *
 * @return random number seed
 */
    
unsigned long int BrookRandomNumberGenerator::getRandomNumberSeed( void ) const {
   return _randomNumberSeed;
}
          
/** 
 * Increment random number seed
 *
 * @param increment    amount to increment random number seed; default = 1
 *
 * @return updated random number seed
 */
     
unsigned long int BrookRandomNumberGenerator::incrementRandomNumberSeed( unsigned long int  increment ){
   _randomNumberSeed += increment;
   return _randomNumberSeed;
}

/** 
 * Set random number seed
 *
 * @param new random number seed; default = 1
 *
 * @return random number seed
 */
    
Mark Friedrichs's avatar
Mark Friedrichs committed
171
172
unsigned long int BrookRandomNumberGenerator::setRandomNumberSeed( unsigned long int seed ){
   _randomNumberSeed  = seed;
Mark Friedrichs's avatar
Mark Friedrichs committed
173
174
175
176
177
178
179
180
181
182
183
184
   return _randomNumberSeed;
}

/** 
 * Generate a random number using algorithm in Gromacs
 * 
 * @param ig seed
 *
 * @return  random number
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
185
BrookOpenMMFloat BrookRandomNumberGenerator::_generateGromacsRandomNumber( unsigned long int* ig ){
Mark Friedrichs's avatar
Mark Friedrichs committed
186
187
188

   // ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
189
   // static const std::string methodName = "\nBrookRandomNumberGenerator::_generateGromacsRandomNumber";
Mark Friedrichs's avatar
Mark Friedrichs committed
190
191
192
193
194
195
196
197
198
199
200
201
202

   int  irand;
 
   int  m     = 100000000;
   float rm   = 100000000.0;  /* same number as m, but real format */
   int  m1    = 10000;
   int  mult  = 31415821;
   
   BrookOpenMMFloat r;
   int  irandh,irandl,multh,multl;
 
   // ---------------------------------------------------------------------------------------
 
Mark Friedrichs's avatar
Mark Friedrichs committed
203
204
   unsigned long int igg = (*ig > 0) ? *ig : -1*(*ig);
   irand = igg % m; 
Mark Friedrichs's avatar
Mark Friedrichs committed
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
   
   /* multiply irand by mult, but take into account that overflow
    * must be discarded, and do not generate an error.
    */

   irandh = irand / m1;
   irandl = irand % m1;
   multh  = mult / m1;
   multl  = mult % m1;
   irand  = ((irandh*multl+irandl*multh) % m1) * m1 + irandl*multl;
   irand  = (irand + 1) % m; 
 
   /* convert irand to a real random number between 0 and 1. */

   r = (BrookOpenMMFloat) (irand / 10); 
   r = r * 10 / rm;
   if ((r <= 0) || (r > 1))
     r = 0.0; 
   *ig = irand;
   
   return r;
}     

inline int MaxInt( unsigned int x, unsigned int y ){ return x > y ? x : y; }

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
231
 * Generate a random number using algorithm in Kiss code
Mark Friedrichs's avatar
Mark Friedrichs committed
232
233
234
235
236
237
238
239
240
 * http://www.helsbreth.org/random/rng_kiss.html
 * 
 * @param randomV1   output random value
 * @param randomV2   output random value
 * @param randomV3   output random value
 * @param state      state
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
241
void BrookRandomNumberGenerator::_generateRandomsKiss( float* randomV1, float* randomV2, float* randomV3, 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
242
                                                       unsigned int state[4] ){
Mark Friedrichs's avatar
Mark Friedrichs committed
243
244
245
    
   // ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
246
   // static const std::string methodName = "\nBrookRandomNumberGenerator::_generateRandomsKiss";
Mark Friedrichs's avatar
Mark Friedrichs committed
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348

   unsigned int carry          = 0;

   // ---------------------------------------------------------------------------------------
 
   state[0]             = state[0] * 69069 + 1;
   state[1]            ^= state[1] << 13;
   state[1]            ^= state[1] >> 17;
   state[1]            ^= state[1] << 5;
   unsigned int k       = (state[2] >> 2) + (state[3] >> 3) + (carry >> 2);
   unsigned int m       = state[3] + state[3] + state[2] + carry;
   state[2]             = state[3];
   state[3]             = m;
   carry                = k >> 30;
   unsigned int z1      = MaxInt(state[0] + state[1] + state[3], 0x00000001);
   float x1             = (float) ( (double) z1 / (double)UINT_MAX );

/*
if( x1 < 0.0f ){
   unsigned int z1 = MaxInt(state[0] + state[1] + state[3], 0x00000001);
   double       z2 = (double) z1/((double) UINT_MAX);
   (void) fprintf( logFile, "x1=%.6e state[%u %u %u] sum %u %.6e den=%u z2=%.6e\n",
                   x1, state[0], state[1], state[3], 
                   (state[0] + state[1] + state[3]), (float) (state[0] + state[1] + state[3]), z1, z2 );
}
*/

   state[0]             = state[0] * 69069 + 1;
   state[1]            ^= state[1] << 13;
   state[1]            ^= state[1] >> 17;
   state[1]            ^= state[1] << 5;
   x1                   = sqrt(-2.0f * log(x1));
   k                    = (state[2] >> 2) + (state[3] >> 3) + (carry >> 2);
   m                    = state[3] + state[3] + state[2] + carry;
   state[2]             = state[3];
   state[3]             = m;
   carry                = k >> 30;
   float x2             = (float)(state[0] + state[1] + state[3]) / (float)UINT_MAX;
   
   state[0]             = state[0] * 69069 + 1;
   state[1]            ^= state[1] << 13;
   state[1]            ^= state[1] >> 17;
   state[1]            ^= state[1] << 5;
   *randomV1            = x1 * cos(2.0f * 3.14159265f * x2);

   k                    = (state[2] >> 2) + (state[3] >> 3) + (carry >> 2);
   m                    = state[3] + state[3] + state[2] + carry;
   state[2]             = state[3];
   state[3]             = m;
   carry                = k >> 30;

   unsigned int z3      = MaxInt(state[0] + state[1] + state[3], 0x00000001);
   float x3             = (float) ( (double) z3 / (double)UINT_MAX );
   //float x3             = (float)MaxInt(state[0] + state[1] + state[3], 0x00000001) / (float)UINT_MAX;

   state[0]             = state[0] * 69069 + 1;
   state[1]            ^= state[1] << 13;
   state[1]            ^= state[1] >> 17;
   state[1]            ^= state[1] << 5;
   x3                   = sqrt(-2.0f * log(x3));
   k                    = (state[2] >> 2) + (state[3] >> 3) + (carry >> 2);
   m                    = state[3] + state[3] + state[2] + carry;
   state[2]             = state[3];
   state[3]             = m;
   carry                = k >> 30;
   float x4             = (float)(state[0] + state[1] + state[3]) / (float)UINT_MAX;
   
   state[0]             = state[0] * 69069 + 1;
   state[1]            ^= state[1] << 13;
   state[1]            ^= state[1] >> 17;
   state[1]            ^= state[1] << 5;
   *randomV2            = x3 * cos(2.0f * 3.14159265f * x4);
   k                    = (state[2] >> 2) + (state[3] >> 3) + (carry >> 2);
   m                    = state[3] + state[3] + state[2] + carry;
   state[2]             = state[3];
   state[3]             = m;
   carry                = k >> 30;

   //float x5             = (float)MaxInt(state[0] + state[1] + state[3], 0x00000001) / (float)UINT_MAX;
   unsigned int z5      = MaxInt(state[0] + state[1] + state[3], 0x00000001);
   float x5             = (float) ( (double) z5 / (double)UINT_MAX );

   state[0]             = state[0] * 69069 + 1;
   state[1]            ^= state[1] << 13;
   state[1]            ^= state[1] >> 17;
   state[1]            ^= state[1] << 5;
   x5                   = sqrt(-2.0f * log(x5));
   k                    = (state[2] >> 2) + (state[3] >> 3) + (carry >> 2);
   m                    = state[3] + state[3] + state[2] + carry;
   state[2]             = state[3];
   state[3]             = m;
   carry                = k >> 30;
   float x6             = (float)(state[0] + state[1] + state[3]) / (float)UINT_MAX;
   *randomV3            = x5 * cos(2.0f * 3.14159265f * x6); 
   
//(void) fprintf( logFile, "rv=%.6e %.6e %.6e\n", *randomV1, *randomV2, *randomV3 );
//exit(0);

   return;
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
349
 * Load random number streams using Kiss algorithm
Mark Friedrichs's avatar
Mark Friedrichs committed
350
351
352
353
354
 * 
 *
 * @return DefaultReturnValue;
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
355
int BrookRandomNumberGenerator::_loadRandomNumberStreamsKiss( void ){
Mark Friedrichs's avatar
Mark Friedrichs committed
356
357
358
359
360

   // ---------------------------------------------------------------------------------------

   static unsigned int state[4];
   static int stateInitialized    = 0;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
361
362
   int printOn                    = 0;
   FILE* log;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
363
   static const int reseed        = 10000;
Mark Friedrichs's avatar
Mark Friedrichs committed
364
 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
365
   static std::string methodName  = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsKiss";
Mark Friedrichs's avatar
Mark Friedrichs committed
366
367
368

   // ---------------------------------------------------------------------------------------
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
369
370
371
372
373
374
   if( printOn && getLog() ){
       log = getLog();
   } else {
      printOn = 0;
   }   

Mark Friedrichs's avatar
Mark Friedrichs committed
375
376
377
378
379
380
381
382
383
   // periodically reset seeds

   if( !stateInitialized || !(stateInitialized % reseed) ){

      state[0] = rand();
      state[1] = rand();
      state[2] = rand();
      state[3] = rand();

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
384
385
      if( printOn ){
         (void) fprintf( log, "%s reset state seeds stateInitialized=%d reseed=%d [%u %u %u %u]\n",
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
386
                         methodName.c_str(), stateInitialized, reseed,  state[0], state[1], state[2], state[3] );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
387
         (void) fflush( log );
Mark Friedrichs's avatar
Mark Friedrichs committed
388
389
390
391
392
393
      }
   }
   stateInitialized++;

   // allocate memory once for download of random nos.

Mark Friedrichs's avatar
Mark Friedrichs committed
394
   float* loadBuffer = _getLoadBuffer();
Mark Friedrichs's avatar
Mark Friedrichs committed
395

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
396
   if( printOn ){
Mark Friedrichs's avatar
Mark Friedrichs committed
397
	   static float count   = 0.0f;
Mark Friedrichs's avatar
Mark Friedrichs committed
398
      float block          = (float) (3*getRandomNumberStreamSize() );
Mark Friedrichs's avatar
Mark Friedrichs committed
399
      count               += 1.0f;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
400
      (void) fprintf( log, "%s: count=%.1f ttl=%.3e no./count=%.1f %d %d\n", methodName.c_str(),
Mark Friedrichs's avatar
Mark Friedrichs committed
401
                      count, block*count, block, getRandomNumberStreamSize(), getNumberOfRandomNumberStreams() );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
402
      (void) fflush( log );
Mark Friedrichs's avatar
Mark Friedrichs committed
403
404
405
406
407
   }
	
   for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){
      for( int ii = 0; ii < 3*getRandomNumberStreamSize(); ii += 3 ){
         float v1,v2,v3;
Mark Friedrichs's avatar
Mark Friedrichs committed
408
409
410
411
         _generateRandomsKiss( &v1, &v2, &v3, state );
         loadBuffer[ii]   = v1;
         loadBuffer[ii+1] = v2;
         loadBuffer[ii+2] = v3;
Mark Friedrichs's avatar
Mark Friedrichs committed
412
   	}
Mark Friedrichs's avatar
Mark Friedrichs committed
413
	   getRandomNumberStream( jj )->loadFromArray( loadBuffer );
Mark Friedrichs's avatar
Mark Friedrichs committed
414
415
   }

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
416
   if( printOn ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
417
      FILE* log = getLog() ? getLog() : stderr;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
418
419
      (void) fprintf( log, "%s: stats\n%s\n", methodName.c_str(), getStatisticsString().c_str() );
      (void) fflush( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
   }

   return DefaultReturnValue;
}

/** 
 * Load random number streams using Mersenne algorithm
 * 
 *
 * @return DefaultReturnValue;
 */

int BrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne( void ){

   // ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
436
   static const std::string methodName  = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne";
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
437
   int printOn                          = 0;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453

   // ---------------------------------------------------------------------------------------
   
   // allocate memory once for download of random nos.

   float* loadBuffer = _getLoadBuffer();

   for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){
      for( int ii = 0; ii < 3*getRandomNumberStreamSize(); ii += 3 ){
         loadBuffer[ii]   = (float) SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
         loadBuffer[ii+1] = (float) SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
         loadBuffer[ii+2] = (float) SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
   	}
	   getRandomNumberStream( jj )->loadFromArray( loadBuffer );
   }

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
454
   if( printOn ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
455
      FILE* log = getLog() ? getLog() : stderr;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
456
457
      (void) fprintf( log, "%s: stats\n%s\n", methodName.c_str(), getStatisticsString().c_str() );
      (void) fflush( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
458
459
   }

Mark Friedrichs's avatar
Mark Friedrichs committed
460
461
462
   return DefaultReturnValue;
}

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
463
464
465
466
467
468
469
470
471
472
473
474
/** 
 * Load random number streams w/ fixed value
 * 
 *
 * @return DefaultReturnValue;
 */

int BrookRandomNumberGenerator::_loadRandomNumberStreamsFixedValue( void ){

   // ---------------------------------------------------------------------------------------

   static const std::string methodName  = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsFixedValue";
Mark Friedrichs's avatar
Mark Friedrichs committed
475
   int printOn                          = 0;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494

   // ---------------------------------------------------------------------------------------
   
   // load fixed value

   float fixedValue  = 0.1f;
   for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){
	   getRandomNumberStream( jj )->fillWithValue( &fixedValue );
   }

   if( printOn ){
      FILE* log = getLog() ? getLog() : stderr;
      (void) fprintf( log, "%s: stats\n%s\n", methodName.c_str(), getStatisticsString().c_str() );
      (void) fflush( log );
   }

   return DefaultReturnValue;
}

Mark Friedrichs's avatar
Mark Friedrichs committed
495
496
497
498
499
500
501
/** 
 * Load random number streams using original gpu algorithm
 * 
 *
 * @return DefaultReturnValue;
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
502
int BrookRandomNumberGenerator::_loadGVStreamsOriginal( void ){
Mark Friedrichs's avatar
Mark Friedrichs committed
503
504
505
506
507

   // ---------------------------------------------------------------------------------------

	unsigned long int jran;

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
508
   // static const std::string methodName = "\nBrookRandomNumberGenerator::_loadGVStreamsOriginal";
Mark Friedrichs's avatar
Mark Friedrichs committed
509
510
511

   // ---------------------------------------------------------------------------------------
	
Mark Friedrichs's avatar
Mark Friedrichs committed
512
   float* loadBuffer = _getLoadBuffer();
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
513
	jran              = getRandomNumberSeed();
Mark Friedrichs's avatar
Mark Friedrichs committed
514
515
	
   for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){
Mark Friedrichs's avatar
Mark Friedrichs committed
516
517
      for( int ii = 0; ii < 3*getRandomNumberStreamSize(); ii++ ){
			loadBuffer[ii] = _generateGromacsRandomNumber( &jran );
Mark Friedrichs's avatar
Mark Friedrichs committed
518
		}
Mark Friedrichs's avatar
Mark Friedrichs committed
519
	   getRandomNumberStream( jj )->loadFromArray( loadBuffer );
Mark Friedrichs's avatar
Mark Friedrichs committed
520
521
522
523
524
525
526
	}
	
   incrementRandomNumberSeed( 1 );

   return DefaultReturnValue;
}

Mark Friedrichs's avatar
Mark Friedrichs committed
527
528
529
530
531
532
533
534
535
536
537
538
539
/** 
 * Get load buffer
 *
 * @return ptr to load buffer
 *
 * @throw OpenMMException if rv stream size is < 1
 *
 **/

float* BrookRandomNumberGenerator::_getLoadBuffer( void ){

   // ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
540
   static const std::string methodName = "\nBrookRandomNumberGenerator::_getLoadBuffer";
Mark Friedrichs's avatar
Mark Friedrichs committed
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571

   // ---------------------------------------------------------------------------------------
	
   if( getRandomNumberStreamSize() < 1 ){
      std::stringstream message;
      message << methodName << " rv stream size=" << getRandomNumberStreamSize() << " is less than 1.";
      throw OpenMMException( message.str() );
      return NULL;
   }

   if( _loadBuffer == NULL ){
   	_loadBuffer = (float*) malloc( sizeof(float)*3*getRandomNumberStreamSize() );
   }

   return _loadBuffer;

}	

/** 
 * Get ptr to shuffle indices
 *
 * @return ptr to shuffle indices
 *
 * @throw OpenMMException if size is < 1
 *
 **/

int* BrookRandomNumberGenerator::_getShuffleIndices( int size ){

   // ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
572
   static const std::string methodName = "\nBrookRandomNumberGenerator::_getShuffleIndices";
Mark Friedrichs's avatar
Mark Friedrichs committed
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590

   // ---------------------------------------------------------------------------------------
	
   if( size < 1 ){
      std::stringstream message;
      message << methodName << " size=" << size << " is less than 1.";
      throw OpenMMException( message.str() );
      return NULL;
   }

   if( _shuffleIndices == NULL ){
   	_shuffleIndices = (int*) malloc( sizeof(int)*size );
   }

   return _shuffleIndices;

}	

Mark Friedrichs's avatar
Mark Friedrichs committed
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
/** 
 * Loads a permutation of indices from 0 to gvSize-1 in
 * sdp->strShuffle. To make sure that the order of the
 * permutation is atleast NGVSHUFFLE, we create the
 * permutation by introducing a random number of p-cycles
 * where p is randomly determined from 2,3,5,7 and 11.
 * The LCM of these numbers is 2310. 
 * Ofcourse the p-cycles are not necessarily disjoint
 * the way it's done here, but there's a good chance 
 * there will enough disjoint cycles to make the 
 * order of the permutation larger than NGVSHUFFLE
 *
 *
 * This function is only called once at startup
 *
 * @return DefaultReturnValue;
 **/

Mark Friedrichs's avatar
Mark Friedrichs committed
609
int BrookRandomNumberGenerator::_loadGVShuffle( void ){
Mark Friedrichs's avatar
Mark Friedrichs committed
610
611
612

   // ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
613
614
615
	const int p[]   = { 2, 3, 5, 7, 11 };
	const int np    = sizeof(p) / sizeof(p[0]);
   const int pmax  = p[np-1];
Mark Friedrichs's avatar
Mark Friedrichs committed
616

Mark Friedrichs's avatar
Mark Friedrichs committed
617
   static const std::string methodName  = "BrookRandomNumberGenerator::loadGVShuffle";
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
618
   int printOn                          = 0;
Mark Friedrichs's avatar
Mark Friedrichs committed
619
620
621

   // ---------------------------------------------------------------------------------------
	
Mark Friedrichs's avatar
Mark Friedrichs committed
622
623
	float* loadBuffer  = _getLoadBuffer();
   int* indices       = _getShuffleIndices( pmax );
Mark Friedrichs's avatar
Mark Friedrichs committed
624
625
	
   int rvSize = getRandomNumberStreamSize();
Mark Friedrichs's avatar
Mark Friedrichs committed
626
627
	for ( int ii = 0; ii < rvSize; ii++ ) {
		loadBuffer[ii] = (float) ii;
Mark Friedrichs's avatar
Mark Friedrichs committed
628
629
630
631
632
	}

	//How to come up with this number here?

   unsigned long int seed = getRandomNumberSeed();
Mark Friedrichs's avatar
Mark Friedrichs committed
633
634
	for( int iter = 0; iter < 1000000; iter++ ) {

Mark Friedrichs's avatar
Mark Friedrichs committed
635
		//for each p
Mark Friedrichs's avatar
Mark Friedrichs committed
636
637
638

		for( int ii = 0; ii < np; ii++ ){

Mark Friedrichs's avatar
Mark Friedrichs committed
639
640
641
			//Come up with p random indices
			//Not checking that they are distinct
			//because that should be fairly rare
Mark Friedrichs's avatar
Mark Friedrichs committed
642
643

			for ( int jj = 0; jj < p[ii]; jj++ ) {
Mark Friedrichs's avatar
Mark Friedrichs committed
644
				//indices[j] = (int) ( gmx_rando( &sdp->seed ) * sdp->gvSize );
Mark Friedrichs's avatar
Mark Friedrichs committed
645
				indices[jj] = (int) ( _generateGromacsRandomNumber( &seed )*rvSize );
Mark Friedrichs's avatar
Mark Friedrichs committed
646
			}
Mark Friedrichs's avatar
Mark Friedrichs committed
647
648
649
650
651
652

			// do a p-cycle

			float tmp = loadBuffer[ indices[0] ];
			for ( int jj = 0; jj < p[ii]-1; jj++ ) {
				loadBuffer[ indices[jj] ] = loadBuffer[ indices[jj+1] ];
Mark Friedrichs's avatar
Mark Friedrichs committed
653
			}
Mark Friedrichs's avatar
Mark Friedrichs committed
654
			loadBuffer[ indices[ p[ii]-1 ] ] = tmp;
Mark Friedrichs's avatar
Mark Friedrichs committed
655
656
		}
	}
Mark Friedrichs's avatar
Mark Friedrichs committed
657
   _getShuffleStream()->loadFromArray( loadBuffer );
Mark Friedrichs's avatar
Mark Friedrichs committed
658
	
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
   if( printOn ){
      FILE* log = getLog() ? getLog() : stderr;
      (void) fprintf( log, "%s: sz=%d pmax=%d np=%d sample indices:\n", methodName.c_str(), rvSize, pmax, np );
      float maxIdx = -1.0f;
      float minIdx =  (float) (rvSize)*2.0f;
      for( int ii = 0; ii < rvSize; ii++ ){
         if( ii < 30 || ii > (rvSize-30) ){
            (void) fprintf( log, "  %d %.8f\n", ii, loadBuffer[ii] );
         }
         if( loadBuffer[ii] < minIdx ){
            minIdx = loadBuffer[ii];
         }
         if( loadBuffer[ii] > maxIdx ){
            maxIdx = loadBuffer[ii];
         }
      }
      (void) fprintf( log, "%s: min-max indices: %.8f %.8f\n", methodName.c_str(), minIdx, maxIdx );
   }
      
Mark Friedrichs's avatar
Mark Friedrichs committed
678
679
680
681
682
683
684
685
686
   return DefaultReturnValue;
}

/** 
 * Shuffle streams
 *
 * @return DefaultReturnValue;
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
687
int BrookRandomNumberGenerator::_shuffleGVStreams( void ){
Mark Friedrichs's avatar
Mark Friedrichs committed
688
689
690

   // ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
691
   // static const std::string methodName = "\nBrookRandomNumberGenerator::_shuffleGVStreams";
Mark Friedrichs's avatar
Mark Friedrichs committed
692
693
694
695

   // ---------------------------------------------------------------------------------------
	
   int numberOfRvStreams = getNumberOfRandomNumberStreams();
Mark Friedrichs's avatar
Mark Friedrichs committed
696

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
697
	for( int ii = 0; ii < (numberOfRvStreams - 1); ii++ ){
Mark Friedrichs's avatar
Mark Friedrichs committed
698
699
700
701
		kpermute_vectors( (float) getRandomNumberStreamWidth(),
                         _getShuffleStream()->getBrookStream(),
                         getRandomNumberStream( ii + 1 )->getBrookStream(), 
                         getRandomNumberStream( ii  )->getBrookStream() );
Mark Friedrichs's avatar
Mark Friedrichs committed
702
703
   }

Mark Friedrichs's avatar
Mark Friedrichs committed
704
705
706
707
	kpermute_vectors( (float) getRandomNumberStreamWidth(),
                      _getShuffleStream()->getBrookStream(),
                      getRandomNumberStream( 0 )->getBrookStream(), 
                      getRandomNumberStream( numberOfRvStreams - 1 )->getBrookStream() );
Mark Friedrichs's avatar
Mark Friedrichs committed
708
709
710
711

   return DefaultReturnValue;
}

Mark Friedrichs's avatar
Mark Friedrichs committed
712
/** 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
713
 * Advances the current position by 2*gpu->particles
Mark Friedrichs's avatar
Mark Friedrichs committed
714
715
716
 * If we run out of rand numbers, we shuffle and 
 * reuse a few times before reloading from the cpu
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
717
 * @param numberOfRandomValuesConsumedPerIteration    number of random values consumed/iteration
Mark Friedrichs's avatar
Mark Friedrichs committed
718
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
719
 * @return DefaultReturnValue;
Mark Friedrichs's avatar
Mark Friedrichs committed
720
721
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
722
int BrookRandomNumberGenerator::advanceGVCursor( int numberOfRandomValuesConsumedPerIteration ){
Mark Friedrichs's avatar
Mark Friedrichs committed
723
724
725

   // ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
726
   static const std::string methodName = "BrookRandomNumberGenerator::advanceGVCursor";
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
727
   int printOn                         = 0;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
728
   FILE* log;
Mark Friedrichs's avatar
Mark Friedrichs committed
729
730

   // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
731
	
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
732
//setLog( stderr );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
733
   if( printOn && getLog() ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
734
      log = getLog();
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
735
736
737
738
   } else {
      printOn = 0;
   }   

Mark Friedrichs's avatar
Mark Friedrichs committed
739
   int rvStreamSize  = getRandomNumberStreamSize();
Mark Friedrichs's avatar
Mark Friedrichs committed
740

Mark Friedrichs's avatar
Mark Friedrichs committed
741
   // use 2 random values per sd 
Mark Friedrichs's avatar
Mark Friedrichs committed
742

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
743
	_rvStreamOffset  += numberOfRandomValuesConsumedPerIteration;
Mark Friedrichs's avatar
Mark Friedrichs committed
744

Mark Friedrichs's avatar
Mark Friedrichs committed
745
	//Check if we've used up this texture
Mark Friedrichs's avatar
Mark Friedrichs committed
746

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
747
   char* action = "none";
Mark Friedrichs's avatar
Mark Friedrichs committed
748
	if ( _rvStreamOffset > rvStreamSize - numberOfRandomValuesConsumedPerIteration ){
Mark Friedrichs's avatar
Mark Friedrichs committed
749

Mark Friedrichs's avatar
Mark Friedrichs committed
750
		// next one if available
Mark Friedrichs's avatar
Mark Friedrichs committed
751

Mark Friedrichs's avatar
Mark Friedrichs committed
752
      _rvStreamOffset = 0;
Mark Friedrichs's avatar
Mark Friedrichs committed
753

Mark Friedrichs's avatar
Mark Friedrichs committed
754
755
		if ( _rvStreamIndex < _numberOfRandomNumberStreams - 1 ){
         _rvStreamIndex++;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
756
         action = "incremented stream index";
Mark Friedrichs's avatar
Mark Friedrichs committed
757
		} else {
Mark Friedrichs's avatar
Mark Friedrichs committed
758

Mark Friedrichs's avatar
Mark Friedrichs committed
759
			//No more textures, need to shuffle
Mark Friedrichs's avatar
Mark Friedrichs committed
760

Mark Friedrichs's avatar
Mark Friedrichs committed
761
			if( _numberOfShuffles < _maxShuffles ){
Mark Friedrichs's avatar
Mark Friedrichs committed
762

Mark Friedrichs's avatar
Mark Friedrichs committed
763
764
				_shuffleGVStreams( );
            _numberOfShuffles++;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
765
            action = "shuffled streams";
Mark Friedrichs's avatar
Mark Friedrichs committed
766

Mark Friedrichs's avatar
Mark Friedrichs committed
767
			} else { //Need to refresh random numbers from cpu
Mark Friedrichs's avatar
Mark Friedrichs committed
768

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
769
            if( _randomNumberGenerator == Mersenne ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
770
               action = "loaded new values to GPU using Mersenne rng";
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
771
772
               _loadRandomNumberStreamsMersenne( );
            } else if( _randomNumberGenerator == Kiss ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
773
               action = "loaded new values to GPU using KISS rng";
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
774
775
               _loadRandomNumberStreamsKiss( );
            } else if( _randomNumberGenerator == Original ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
776
               action = "loaded new values to GPU using original Gromac's rng";
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
777
               _loadGVStreamsOriginal( );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
778
779
780
            } else if( _randomNumberGenerator == FixedValue ){
               action = "loaded new fixed values to GPU";
               _loadRandomNumberStreamsFixedValue( );
Mark Friedrichs's avatar
Mark Friedrichs committed
781
782
783
784
785
            }
            _numberOfShuffles  = 0;
			}
         _rvStreamIndex = 0;
		}
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
786

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
787
      if( printOn ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
788
789
790
         (void) fprintf( log, "%s offset=%d consume/itr=%d StrmSz=%d idx=%d shffle=%d action=%s\n",
                         methodName.c_str(), _rvStreamOffset, numberOfRandomValuesConsumedPerIteration,
                         rvStreamSize, _rvStreamIndex, _numberOfShuffles, action );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
791
         (void) fflush( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
792
      }
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810

/*
      // check rng distribution

      static int count = 0;
      if( count++ < 2 ){

         // accumulate rng -- stats will be in cumulative fields in stat string

         int testIterations = 1000;
         for( int ii = 0; ii < testIterations; ii++ ){
            //_loadRandomNumberStreamsKiss( );
            _loadRandomNumberStreamsMersenne( );
            getStatisticsString();
         }
         (void) fprintf( log, "%s: stats\n%s\n", methodName.c_str(), getStatisticsString().c_str() );
      }
*/
Mark Friedrichs's avatar
Mark Friedrichs committed
811
	}
Mark Friedrichs's avatar
Mark Friedrichs committed
812
813
814
815
816

   return DefaultReturnValue;
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
817
 * Get random number stream size
Mark Friedrichs's avatar
Mark Friedrichs committed
818
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
819
 * @return  random number stream size
Mark Friedrichs's avatar
Mark Friedrichs committed
820
821
822
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
823
824
int BrookRandomNumberGenerator::getRandomNumberStreamSize( void ) const {
   return _randomNumberStreamSize;
Mark Friedrichs's avatar
Mark Friedrichs committed
825
826
827
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
828
 * Get Shuffle stream 
Mark Friedrichs's avatar
Mark Friedrichs committed
829
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
830
 * @return  shuffle stream
Mark Friedrichs's avatar
Mark Friedrichs committed
831
832
833
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
834
835
BrookFloatStreamInternal* BrookRandomNumberGenerator::_getShuffleStream( void ) const {
   return _auxiliaryStreams[ShuffleStream];
Mark Friedrichs's avatar
Mark Friedrichs committed
836
837
838
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
839
 * Get random number stream 
Mark Friedrichs's avatar
Mark Friedrichs committed
840
 *
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
841
 * @param index     random number stream index     
Mark Friedrichs's avatar
Mark Friedrichs committed
842
843
 *
 * @return  random number stream
Mark Friedrichs's avatar
Mark Friedrichs committed
844
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
845
 * @throw OpenMMException if index is invalid or _randomNumberGeneratorStreams not set
Mark Friedrichs's avatar
Mark Friedrichs committed
846
847
848
 *
 */

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
849
BrookFloatStreamInternal* BrookRandomNumberGenerator::getRandomNumberStream( int index ){
Mark Friedrichs's avatar
Mark Friedrichs committed
850

Mark Friedrichs's avatar
Mark Friedrichs committed
851
   // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
852

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
853
   static const std::string methodName = "\nBrookRandomNumberGenerator::getRandomNumberStream";
Mark Friedrichs's avatar
Mark Friedrichs committed
854

Mark Friedrichs's avatar
Mark Friedrichs committed
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
   // ---------------------------------------------------------------------------------------
	
   if( index < 0 || index >= _numberOfRandomNumberStreams ){
      std::stringstream message;
      message << methodName << " index=" << index << " is out of range: [0," << _numberOfRandomNumberStreams;
      throw OpenMMException( message.str() );
      return NULL;
   }

   if( _randomNumberGeneratorStreams == NULL ){
      std::stringstream message;
      message << methodName << " rv streams not initialized; input index=" << index;
      throw OpenMMException( message.str() );
      return NULL;
   }

   return _randomNumberGeneratorStreams[index];
Mark Friedrichs's avatar
Mark Friedrichs committed
872
873
874
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
875
 * Get random value stream index
Mark Friedrichs's avatar
Mark Friedrichs committed
876
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
877
 * @return  random value stream index
Mark Friedrichs's avatar
Mark Friedrichs committed
878
879
880
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
881
882
int BrookRandomNumberGenerator::getRvStreamIndex( void ) const {
   return _rvStreamIndex;
Mark Friedrichs's avatar
Mark Friedrichs committed
883
884
885
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
886
 * Get random value stream offset
Mark Friedrichs's avatar
Mark Friedrichs committed
887
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
888
 * @return  random value stream offset
Mark Friedrichs's avatar
Mark Friedrichs committed
889
890
891
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
892
893
int BrookRandomNumberGenerator::getRvStreamOffset( void ) const {
   return _rvStreamOffset;
Mark Friedrichs's avatar
Mark Friedrichs committed
894
895
896
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
897
 * Get max shuffles
Mark Friedrichs's avatar
Mark Friedrichs committed
898
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
899
 * @return  max shuffles
Mark Friedrichs's avatar
Mark Friedrichs committed
900
901
902
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
903
904
int BrookRandomNumberGenerator::getMaxShuffles( void ) const {
   return _maxShuffles;
Mark Friedrichs's avatar
Mark Friedrichs committed
905
906
907
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
908
 * Get number of shuffles
Mark Friedrichs's avatar
Mark Friedrichs committed
909
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
910
 * @return  number of shuffles
Mark Friedrichs's avatar
Mark Friedrichs committed
911
912
913
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
914
915
int BrookRandomNumberGenerator::_getNumberOfShuffles( void ) const {
   return _numberOfShuffles;
Mark Friedrichs's avatar
Mark Friedrichs committed
916
917
918
919
920
}

/** 
 * Initialize stream dimensions
 * 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
921
 * @param numberOfParticles         number of particles
Mark Friedrichs's avatar
Mark Friedrichs committed
922
923
924
925
926
927
 * @param platform                  platform
 *      
 * @return ErrorReturnValue if error, else DefaultReturnValue
 *
 */

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
928
int BrookRandomNumberGenerator::_initializeStreamSizes( int numberOfParticles, const Platform& platform ){
Mark Friedrichs's avatar
Mark Friedrichs committed
929
930
931
932
933
934
935

// ---------------------------------------------------------------------------------------

   //static const std::string methodName      = "BrookRandomNumberGenerator::_initializeStreamSizes";

// ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
936
   setNumberOfParticles( numberOfParticles );
Mark Friedrichs's avatar
Mark Friedrichs committed
937

Mark Friedrichs's avatar
Mark Friedrichs committed
938
   // get randomNumber stream dimensions
Mark Friedrichs's avatar
Mark Friedrichs committed
939

940
   const BrookPlatform& brookPlatform           = dynamic_cast<const BrookPlatform&> (platform);
Mark Friedrichs's avatar
Mark Friedrichs committed
941
942
943
944
   const BrookStreamFactory& brookStreamFactory = dynamic_cast<const BrookStreamFactory&> (platform.getDefaultStreamFactory() );
   _randomNumberStreamWidth                     = brookStreamFactory.getDefaultRandomNumberStreamWidth();
   _randomNumberStreamSize                      = brookStreamFactory.getDefaultRandomNumberStreamSize();
   _randomNumberStreamHeight                    = (int) ( ((float) _randomNumberStreamSize)/( (float) _randomNumberStreamWidth) + 0.001);
Mark Friedrichs's avatar
Mark Friedrichs committed
945

Mark Friedrichs's avatar
Mark Friedrichs committed
946
   _randomNumberStreamSize                      =  _randomNumberStreamWidth*_randomNumberStreamHeight;
Mark Friedrichs's avatar
Mark Friedrichs committed
947

Mark Friedrichs's avatar
Mark Friedrichs committed
948
   return DefaultReturnValue;
Mark Friedrichs's avatar
Mark Friedrichs committed
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
}

/** 
 * Initialize streams
 * 
 * @param platform                  platform
 *
 * @return ErrorReturnValue if error, else DefaultReturnValue
 *
 */

int BrookRandomNumberGenerator::_initializeStreams( const Platform& platform ){

// ---------------------------------------------------------------------------------------

   //static const std::string methodName      = "BrookRandomNumberGenerator::_initializeStreams";

Mark Friedrichs's avatar
Mark Friedrichs committed
966
   float dangleValue                          = 0.0f;
Mark Friedrichs's avatar
Mark Friedrichs committed
967
968
969

// ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
970
971
   int randomNumberStreamSize                = getRandomNumberStreamSize();
   int randomNumberStreamWidth               = getRandomNumberStreamWidth();
Mark Friedrichs's avatar
Mark Friedrichs committed
972

Mark Friedrichs's avatar
Mark Friedrichs committed
973
974
975
   _auxiliaryStreams[ShuffleStream]         = new BrookFloatStreamInternal( BrookCommon::ShuffleStream,
                                                                            randomNumberStreamSize, randomNumberStreamWidth,
                                                                            BrookStreamInternal::Float, dangleValue );
Mark Friedrichs's avatar
Mark Friedrichs committed
976

Mark Friedrichs's avatar
Mark Friedrichs committed
977
978
979
980
981
   // insure number of random number streams is > 0
   // delete if already allocated and then initialize
 
   if( _numberOfRandomNumberStreams < 1 ){
      _numberOfRandomNumberStreams     = 1;
Mark Friedrichs's avatar
Mark Friedrichs committed
982
983
   }

Mark Friedrichs's avatar
Mark Friedrichs committed
984
985
   if( _randomNumberGeneratorStreams ){
      delete[] _randomNumberGeneratorStreams;
Mark Friedrichs's avatar
Mark Friedrichs committed
986
987
   }

Mark Friedrichs's avatar
Mark Friedrichs committed
988
   _randomNumberGeneratorStreams = new BrookFloatStreamInternal*[_numberOfRandomNumberStreams];
Mark Friedrichs's avatar
Mark Friedrichs committed
989

Mark Friedrichs's avatar
Mark Friedrichs committed
990
   for( int ii = 0; ii < _numberOfRandomNumberStreams; ii++ ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
991
      _randomNumberGeneratorStreams[ii]    = new BrookFloatStreamInternal( BrookCommon::RandomValuesStream,
Mark Friedrichs's avatar
Mark Friedrichs committed
992
993
                                                                           randomNumberStreamSize, randomNumberStreamWidth,
                                                                           BrookStreamInternal::Float3, dangleValue );
Mark Friedrichs's avatar
Mark Friedrichs committed
994
995
   }

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
996
997
998
999
   // create shuffle stream

   _loadGVShuffle();

Mark Friedrichs's avatar
Mark Friedrichs committed
1000
1001
1002
1003
   return DefaultReturnValue;
}

/*  
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1004
 * Setup of streams, ... associated w/ random number generator
Mark Friedrichs's avatar
Mark Friedrichs committed
1005
 *
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1006
 * @param numberOfParticles     number of particles
Mark Friedrichs's avatar
Mark Friedrichs committed
1007
1008
1009
1010
1011
1012
 * @param platform              Brook platform
 *
 * @return nonzero value if error
 *
 * */
    
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1013
int BrookRandomNumberGenerator::setup( int numberOfParticles, const Platform& platform ){
Mark Friedrichs's avatar
Mark Friedrichs committed
1014
1015
1016
    
// ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
1017
   //static const std::string methodName      = "BrookRandomNumberGenerator::setup";
Mark Friedrichs's avatar
Mark Friedrichs committed
1018
1019
1020

// ---------------------------------------------------------------------------------------

1021
   const BrookPlatform& brookPlatform         = dynamic_cast<const BrookPlatform&> (platform);
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1022
1023
   setLog( brookPlatform.getLog() );

Mark Friedrichs's avatar
Mark Friedrichs committed
1024
1025
   // set stream sizes and then create streams

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1026
   _initializeStreamSizes( numberOfParticles, platform );
Mark Friedrichs's avatar
Mark Friedrichs committed
1027
   _initializeStreams( platform );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1028
1029
1030
1031

   if( _randomNumberGenerator == Mersenne ){
      _loadRandomNumberStreamsMersenne( );
   } else if( _randomNumberGenerator == Kiss ){
1032
      _loadRandomNumberStreamsKiss( );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1033
1034
   } else if( _randomNumberGenerator == Original ){
      _loadGVStreamsOriginal( );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1035
1036
   } else if( _randomNumberGenerator == FixedValue ){
      _loadRandomNumberStreamsFixedValue( );
1037
   }
Mark Friedrichs's avatar
Mark Friedrichs committed
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072

   return DefaultReturnValue;
}

/* 
 * Get contents of object
 *
 * @param level   level of dump
 *
 * @return string containing contents
 *
 * */

std::string BrookRandomNumberGenerator::getContentsString( int level ) const {

// ---------------------------------------------------------------------------------------

   static const std::string methodName      = "BrookRandomNumberGenerator::getContentsString";

   static const unsigned int MAX_LINE_CHARS = 256;
   char value[MAX_LINE_CHARS];
   static const char* Set                   = "Set";
   static const char* NotSet                = "Not set";

// ---------------------------------------------------------------------------------------

   std::stringstream message;
   std::string tab   = "   ";

#ifdef WIN32
#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );   
#else
#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );   
#endif

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1073
1074
1075
1076
1077
1078
1079
   
   if( _randomNumberGenerator == Mersenne ){
      (void) LOCAL_SPRINTF( value, "%s", "Mersenne Rng");
   } else if( _randomNumberGenerator == Kiss ){
      (void) LOCAL_SPRINTF( value, "%s", "Kiss Rng");
   } else if( _randomNumberGenerator == Original ){
      (void) LOCAL_SPRINTF( value, "%s", "Original Gromacs Rng");
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1080
1081
   } else if( _randomNumberGenerator == FixedValue ){
      (void) LOCAL_SPRINTF( value, "%s", "Fixed value Rng");
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1082
   }
Mark Friedrichs's avatar
Mark Friedrichs committed
1083
   message << _getLine( tab, "Random number generator:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1084

1085
1086
   (void) LOCAL_SPRINTF( value, "%d", getRandomNumberStreamSize() );
   message << _getLine( tab, "RandomNumber stream size:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1087

Mark Friedrichs's avatar
Mark Friedrichs committed
1088
1089
   (void) LOCAL_SPRINTF( value, "%d", getRandomNumberStreamWidth() );
   message << _getLine( tab, "RandomNumber stream width:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1090

Mark Friedrichs's avatar
Mark Friedrichs committed
1091
1092
   (void) LOCAL_SPRINTF( value, "%d", getRandomNumberStreamHeight() );
   message << _getLine( tab, "RandomNumber stream height:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1093

Mark Friedrichs's avatar
Mark Friedrichs committed
1094
1095
   (void) LOCAL_SPRINTF( value, "%d", getNumberOfRandomNumberStreams() );
   message << _getLine( tab, "Number of RandomNumber streams:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1096

Mark Friedrichs's avatar
Mark Friedrichs committed
1097
1098
   (void) LOCAL_SPRINTF( value, "%ld", getRandomNumberSeed() );
   message << _getLine( tab, "Seed:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1099

Mark Friedrichs's avatar
Mark Friedrichs committed
1100
1101
   (void) LOCAL_SPRINTF( value, "%d", getRvStreamIndex() );
   message << _getLine( tab, "Rv stream index:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1102

Mark Friedrichs's avatar
Mark Friedrichs committed
1103
1104
1105
1106
1107
   (void) LOCAL_SPRINTF( value, "%d", getRvStreamOffset() );
   message << _getLine( tab, "Rv stream offset:", value ); 

   (void) LOCAL_SPRINTF( value, "%d", _getNumberOfShuffles() );
   message << _getLine( tab, "Number of rv stream shuffles:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1108

Mark Friedrichs's avatar
Mark Friedrichs committed
1109
1110
   (void) LOCAL_SPRINTF( value, "%d", getMaxShuffles() );
   message << _getLine( tab, "Max number of rv stream shuffles:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1111

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1112
   message << _getLine( tab, "Log:",                  (getLog()                 ? Set : NotSet) ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1113

Mark Friedrichs's avatar
Mark Friedrichs committed
1114
   message << _getLine( tab, "Shuffle:",              (_getShuffleStream()      ? Set : NotSet) ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1115
 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1116
1117
   // show stats

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1118
1119
1120
1121
1122
1123
1124
1125
1126
   message << getStatisticsString( );

   for( int ii = 0; ii < LastStreamIndex; ii++ ){
      message << std::endl;
      if( _auxiliaryStreams[ii] ){
         message << _auxiliaryStreams[ii]->getContentsString( );
      }
   }

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1127
1128
1129
1130
1131
1132
   for( int ii = 0; ii < _numberOfRandomNumberStreams; ii++ ){
      message << std::endl;
      if( _randomNumberGeneratorStreams[ii] ){
         message <<  _randomNumberGeneratorStreams[ii]->getContentsString();
      }
   }
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
#undef LOCAL_SPRINTF

   return message.str();
}

/* 
 * Get statistics
 *
 * @return string containing contents
 *
 * */

std::string BrookRandomNumberGenerator::getStatisticsString( void ) const {

// ---------------------------------------------------------------------------------------

   static const std::string methodName      = "BrookRandomNumberGenerator::getStatisticsString";

   static const unsigned int MAX_LINE_CHARS = 256;
   char value[MAX_LINE_CHARS];
   static const char* Set                   = "Set";
   static const char* NotSet                = "Not set";
   static double cumulativeStatistics[7]    = { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0e+99, -1.0e+99 };

// ---------------------------------------------------------------------------------------

   std::stringstream message;
   std::string tab   = "   ";

#ifdef WIN32
#define LOCAL_SPRINTF(a,b,c) sprintf_s( (a), MAX_LINE_CHARS, (b), (c) );   
#else
#define LOCAL_SPRINTF(a,b,c) sprintf( (a), (b), (c) );   
#endif

1168
1169
1170
   message << "\n   Statistics:\n";

   double statistics[7];
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1171
   getStatistics( statistics, -1, cumulativeStatistics );
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190

   (void) LOCAL_SPRINTF( value, "%.5e", statistics[4] );
   message << _getLine( tab, "Count:",  value ); 

   (void) LOCAL_SPRINTF( value, "%.5e", statistics[0] );
   message << _getLine( tab, "Average:",  value ); 

   (void) LOCAL_SPRINTF( value, "%.5e", statistics[1] );
   message << _getLine( tab, "StdDev:",  value ); 

   (void) LOCAL_SPRINTF( value, "%.5e", statistics[3] );
   message << _getLine( tab, "Kurtosis:",  value ); 

   (void) LOCAL_SPRINTF( value, "%.5e", statistics[5] );
   message << _getLine( tab, "Min:",  value ); 

   (void) LOCAL_SPRINTF( value, "%.6e", statistics[6] );
   message << _getLine( tab, "Max:",  value ); 

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1191
1192
1193
1194
   if( cumulativeStatistics[4] > 1000.0 ){ 

      for( int ii = 0; ii < 7; ii++ ){
         statistics[ii] = cumulativeStatistics[ii];
Mark Friedrichs's avatar
Mark Friedrichs committed
1195
1196
      }

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
      statistics[0] /= statistics[4];
      statistics[1]  = statistics[1] - statistics[4]*statistics[0]*statistics[0];
      if( statistics[4] > 1.0 ){
         statistics[1] = sqrt( statistics[1]/( statistics[4] - 1.0 ) );
      }
      statistics[3]  = (statistics[3]/(statistics[4]*statistics[1]*statistics[1]) ) - 3.0;
      (void) LOCAL_SPRINTF( value, "%.5e", statistics[4] );
      message << _getLine( tab, "Cumulative Count:",  value ); 
   
      (void) LOCAL_SPRINTF( value, "%.5e", statistics[0] );
      message << _getLine( tab, "Cumulative Average:",  value ); 
   
      (void) LOCAL_SPRINTF( value, "%.5e", statistics[1] );
      message << _getLine( tab, "Cumulative StdDev:",  value ); 
   
      (void) LOCAL_SPRINTF( value, "%.5e", statistics[3] );
      message << _getLine( tab, "Cumulative Kurtosis:",  value ); 
   
      (void) LOCAL_SPRINTF( value, "%.5e", statistics[5] );
      message << _getLine( tab, "Cumulative Min:",  value ); 
   
      (void) LOCAL_SPRINTF( value, "%.6e", statistics[6] );
      message << _getLine( tab, "Cumulative Max:",  value ); 
   }
   
Mark Friedrichs's avatar
Mark Friedrichs committed
1222
1223
1224
1225
#undef LOCAL_SPRINTF

   return message.str();
}
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244

/* 
 * Get statistics
 *
 * @param statistics   array of size 5:
 *                       0: mean
 *                       1: std dev
 *                       2: 3rd moment (not normalized)
 *                       3: kurtosis
 *                       4: count 
 *                       5: min
 *                       6: max
 *
 * @param streamIndex  stream index to analyze
 *
 * @return DefaultReturnValue
 *
 * */

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1245
int BrookRandomNumberGenerator::getStatistics( double statistics[7], int streamIndex, double cumulativeStatistics[7] ) const {
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263

// ---------------------------------------------------------------------------------------

   static const std::string methodName      = "BrookRandomNumberGenerator::";

// ---------------------------------------------------------------------------------------

   statistics[0] = 0.0;
   statistics[1] = 0.0;
   statistics[2] = 0.0;
   statistics[3] = 0.0;
   statistics[4] = 0.0;
   statistics[5] =  1.0e+99;
   statistics[6] = -1.0e+99;

   for( int ii = 0; ii < _numberOfRandomNumberStreams; ii++ ){
      if( streamIndex < 0 || ii == streamIndex ){
         void* dataArrayV       = _randomNumberGeneratorStreams[ii]->getData( 1 );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1264
         int numberOfValues     = _randomNumberGeneratorStreams[ii]->getStreamSize()*_randomNumberGeneratorStreams[ii]->getWidth();
1265
         const float* dataArray = (float*) dataArrayV;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1266
         for( int ii = 0; ii < numberOfValues; ii++ ){
1267

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1268
             statistics[0] += dataArray[ii];
1269

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1270
             double rv2     = dataArray[ii]*dataArray[ii];
1271
             statistics[1] += rv2;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1272
             statistics[2] += rv2*dataArray[ii];
1273
1274
             statistics[3] += rv2*rv2;

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1275
1276
             if( statistics[5] > dataArray[ii] ){
                statistics[5] = dataArray[ii];
1277
             }
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1278
1279
             if( statistics[6] < dataArray[ii] ){
                statistics[6] = dataArray[ii];
1280
1281
             }
         } 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1282
         statistics[4] += (double) numberOfValues;
1283
1284
1285
      }   
   }   

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1286
1287
   // accumulate moments, ... in cumulativeStatistics array

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
   for( int ii = 0; ii < 5; ii++ ){
      cumulativeStatistics[ii] += statistics[ii];
   }
   if( statistics[5] < cumulativeStatistics[5] ){
      cumulativeStatistics[5] = statistics[5];
   }
   if( statistics[6] > cumulativeStatistics[6] ){
      cumulativeStatistics[6] = statistics[6];
   }

1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
   if( statistics[4] > 0.0 ){ 
      statistics[0] /= statistics[4];
      statistics[1]  = statistics[1] - statistics[4]*statistics[0]*statistics[0];
      if( statistics[4] > 1.0 ){
         statistics[1] = sqrt( statistics[1]/( statistics[4] - 1.0 ) );
      }
      statistics[3]  = (statistics[3]/(statistics[4]*statistics[1]*statistics[1]) ) - 3.0;
   }

   return DefaultReturnValue;
}