BrookRandomNumberGenerator.cpp 41.9 KB
Newer Older
Mark Friedrichs's avatar
Mark Friedrichs committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
/* -------------------------------------------------------------------------- *
 *                                   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.               *
 *                                                                            *
 * Portions copyright (c) 2008 Stanford University and the Authors.           *
 * Authors: Mark Friedrichs                                                   *
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

#include <sstream>
#include "BrookRandomNumberGenerator.h"
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
34
#include "../../reference/src/SimTKUtilities/SimTKOpenMMUtilities.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
35
#include "OpenMMException.h"
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
36
#include "gpu/kupdatesd.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

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
58
   _randomNumberGeneratorStreams    = NULL;
Mark Friedrichs's avatar
Mark Friedrichs committed
59
60
61
62
63
64
65

   // mark stream dimension variables as unset

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

Mark Friedrichs's avatar
Mark Friedrichs committed
66
67
68
   _rvStreamIndex                   = 0;
   _rvStreamOffset                  = 0;
   _numberOfShuffles                = 0;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
69
70
   //_maxShuffles                     = 0;
   _maxShuffles                     = 100;
Mark Friedrichs's avatar
Mark Friedrichs committed
71
72
73
74

   _loadBuffer                      = NULL;
   _shuffleIndices                  = NULL;

Mark Friedrichs's avatar
Mark Friedrichs committed
75
   for( int ii = 0; ii < LastStreamIndex; ii++ ){
Mark Friedrichs's avatar
Mark Friedrichs committed
76
      _auxiliaryStreams[ii]   = NULL;
Mark Friedrichs's avatar
Mark Friedrichs committed
77
78
79
80
   }

   // set randomNumber seed 

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
81
   _randomNumberSeed      = 1393;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
82
   //_randomNumberGenerator = Mersenne;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
83
   _randomNumberGenerator = Kiss;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
84
   //_randomNumberGenerator = FixedValue;
Mark Friedrichs's avatar
Mark Friedrichs committed
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102

   //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
103
      delete _auxiliaryStreams[ii];
Mark Friedrichs's avatar
Mark Friedrichs committed
104
105
   }

Mark Friedrichs's avatar
Mark Friedrichs committed
106
107
108
109
   delete[] _randomNumberGeneratorStreams;

   delete[] _loadBuffer;
   delete[] _shuffleIndices;
Mark Friedrichs's avatar
Mark Friedrichs committed
110
111
112
113
114
115
116
117
118
119
120
121
122
123
}

/** 
 * 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
124
 * Get random number stream width
Mark Friedrichs's avatar
Mark Friedrichs committed
125
 * 
Mark Friedrichs's avatar
Mark Friedrichs committed
126
 * @return     nndom number stream width
Mark Friedrichs's avatar
Mark Friedrichs committed
127
128
129
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
130
131
int BrookRandomNumberGenerator::getRandomNumberStreamWidth( void ) const {
   return _randomNumberStreamWidth;
Mark Friedrichs's avatar
Mark Friedrichs committed
132
}
Mark Friedrichs's avatar
Mark Friedrichs committed
133
134
135
136
137
138
139
140
141
142
143
144

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

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

Mark Friedrichs's avatar
Mark Friedrichs committed
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
171
172
173
174
175
/** 
 * 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
176
177
unsigned long int BrookRandomNumberGenerator::setRandomNumberSeed( unsigned long int seed ){
   _randomNumberSeed  = seed;
Mark Friedrichs's avatar
Mark Friedrichs committed
178
179
180
181
182
183
184
185
186
187
188
189
   return _randomNumberSeed;
}

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

Mark Friedrichs's avatar
Mark Friedrichs committed
190
BrookOpenMMFloat BrookRandomNumberGenerator::_generateGromacsRandomNumber( unsigned long int* ig ){
Mark Friedrichs's avatar
Mark Friedrichs committed
191
192
193

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

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
194
   // static const std::string methodName = "\nBrookRandomNumberGenerator::_generateGromacsRandomNumber";
Mark Friedrichs's avatar
Mark Friedrichs committed
195
196
197
198
199
200
201
202
203
204
205
206
207

   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
208
209
   unsigned long int igg = (*ig > 0) ? *ig : -1*(*ig);
   irand = igg % m; 
Mark Friedrichs's avatar
Mark Friedrichs committed
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
   
   /* 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
236
 * Generate a random number using algorithm in Kiss code
Mark Friedrichs's avatar
Mark Friedrichs committed
237
238
239
240
241
242
243
244
245
 * 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
246
void BrookRandomNumberGenerator::_generateRandomsKiss( float* randomV1, float* randomV2, float* randomV3, 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
247
                                                       unsigned int state[4] ){
Mark Friedrichs's avatar
Mark Friedrichs committed
248
249
250
    
   // ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
251
   // static const std::string methodName = "\nBrookRandomNumberGenerator::_generateRandomsKiss";
Mark Friedrichs's avatar
Mark Friedrichs committed
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
349
350
351
352
353

   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
354
 * Load random number streams using Kiss algorithm
Mark Friedrichs's avatar
Mark Friedrichs committed
355
356
357
358
359
 * 
 *
 * @return DefaultReturnValue;
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
360
int BrookRandomNumberGenerator::_loadRandomNumberStreamsKiss( void ){
Mark Friedrichs's avatar
Mark Friedrichs committed
361
362
363
364
365

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

   static unsigned int state[4];
   static int stateInitialized    = 0;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
366
367
   int printOn                    = 0;
   FILE* log;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
368
   static const int reseed        = 10000;
Mark Friedrichs's avatar
Mark Friedrichs committed
369
 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
370
   static std::string methodName  = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsKiss";
Mark Friedrichs's avatar
Mark Friedrichs committed
371
372
373

   // ---------------------------------------------------------------------------------------
   
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
374
375
376
377
378
379
   if( printOn && getLog() ){
       log = getLog();
   } else {
      printOn = 0;
   }   

Mark Friedrichs's avatar
Mark Friedrichs committed
380
381
382
383
384
385
386
387
388
   // 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
389
390
      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
391
                         methodName.c_str(), stateInitialized, reseed,  state[0], state[1], state[2], state[3] );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
392
         (void) fflush( log );
Mark Friedrichs's avatar
Mark Friedrichs committed
393
394
395
396
397
398
399
400
401
402
403
404
405
406
      }

/*
state[0] = 9578;
state[1] = 29245;
state[2] = 16266;
state[3] = 27587;
*/

   }
   stateInitialized++;

   // allocate memory once for download of random nos.

Mark Friedrichs's avatar
Mark Friedrichs committed
407
   float* loadBuffer = _getLoadBuffer();
Mark Friedrichs's avatar
Mark Friedrichs committed
408

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
409
   if( printOn ){
Mark Friedrichs's avatar
Mark Friedrichs committed
410
	   static float count   = 0.0f;
Mark Friedrichs's avatar
Mark Friedrichs committed
411
      float block          = (float) (3*getRandomNumberStreamSize() );
Mark Friedrichs's avatar
Mark Friedrichs committed
412
      count               += 1.0f;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
413
      (void) fprintf( log, "%s: count=%.1f ttl=%.3e no./count=%.1f %d %d\n", methodName.c_str(),
Mark Friedrichs's avatar
Mark Friedrichs committed
414
                      count, block*count, block, getRandomNumberStreamSize(), getNumberOfRandomNumberStreams() );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
415
      (void) fflush( log );
Mark Friedrichs's avatar
Mark Friedrichs committed
416
417
418
419
420
   }
	
   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
421
422
423
424
         _generateRandomsKiss( &v1, &v2, &v3, state );
         loadBuffer[ii]   = v1;
         loadBuffer[ii+1] = v2;
         loadBuffer[ii+2] = v3;
Mark Friedrichs's avatar
Mark Friedrichs committed
425
   	}
Mark Friedrichs's avatar
Mark Friedrichs committed
426
	   getRandomNumberStream( jj )->loadFromArray( loadBuffer );
Mark Friedrichs's avatar
Mark Friedrichs committed
427
428
   }

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
429
   if( printOn ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
430
      FILE* log = getLog() ? getLog() : stderr;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
431
432
      (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
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
   }

   return DefaultReturnValue;
}


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

int BrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne( void ){

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

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
450
   static const std::string methodName  = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsMersenne";
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
451
   int printOn                          = 0;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467

   // ---------------------------------------------------------------------------------------
   
   // 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
468
   if( printOn ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
469
      FILE* log = getLog() ? getLog() : stderr;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
470
471
      (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
472
473
   }

Mark Friedrichs's avatar
Mark Friedrichs committed
474
475
476
   return DefaultReturnValue;
}

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
/** 
 * Load random number streams w/ fixed value
 * 
 *
 * @return DefaultReturnValue;
 */

int BrookRandomNumberGenerator::_loadRandomNumberStreamsFixedValue( void ){

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

   static const std::string methodName  = "\nBrookRandomNumberGenerator::_loadRandomNumberStreamsFixedValue";
   int printOn                          = 1;

   // ---------------------------------------------------------------------------------------
   
   // 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
509
510
511
512
513
514
515
/** 
 * Load random number streams using original gpu algorithm
 * 
 *
 * @return DefaultReturnValue;
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
516
int BrookRandomNumberGenerator::_loadGVStreamsOriginal( void ){
Mark Friedrichs's avatar
Mark Friedrichs committed
517
518
519
520
521

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

	unsigned long int jran;

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
522
   // static const std::string methodName = "\nBrookRandomNumberGenerator::_loadGVStreamsOriginal";
Mark Friedrichs's avatar
Mark Friedrichs committed
523
524
525

   // ---------------------------------------------------------------------------------------
	
Mark Friedrichs's avatar
Mark Friedrichs committed
526
   float* loadBuffer = _getLoadBuffer();
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
527
	jran              = getRandomNumberSeed();
Mark Friedrichs's avatar
Mark Friedrichs committed
528
529
	
   for( int jj = 0; jj < getNumberOfRandomNumberStreams(); jj++ ){
Mark Friedrichs's avatar
Mark Friedrichs committed
530
531
      for( int ii = 0; ii < 3*getRandomNumberStreamSize(); ii++ ){
			loadBuffer[ii] = _generateGromacsRandomNumber( &jran );
Mark Friedrichs's avatar
Mark Friedrichs committed
532
		}
Mark Friedrichs's avatar
Mark Friedrichs committed
533
	   getRandomNumberStream( jj )->loadFromArray( loadBuffer );
Mark Friedrichs's avatar
Mark Friedrichs committed
534
535
536
537
538
539
540
	}
	
   incrementRandomNumberSeed( 1 );

   return DefaultReturnValue;
}

Mark Friedrichs's avatar
Mark Friedrichs committed
541
542
543
544
545
546
547
548
549
550
551
552
553
/** 
 * 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
554
   static const std::string methodName = "\nBrookRandomNumberGenerator::_getLoadBuffer";
Mark Friedrichs's avatar
Mark Friedrichs committed
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585

   // ---------------------------------------------------------------------------------------
	
   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
586
   static const std::string methodName = "\nBrookRandomNumberGenerator::_getShuffleIndices";
Mark Friedrichs's avatar
Mark Friedrichs committed
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604

   // ---------------------------------------------------------------------------------------
	
   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
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
/** 
 * 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
623
int BrookRandomNumberGenerator::_loadGVShuffle( void ){
Mark Friedrichs's avatar
Mark Friedrichs committed
624
625
626
627
628
629
630

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

	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
Mods  
Mark Friedrichs committed
631
632
   static const std::string methodName  = "\nBrookRandomNumberGenerator::loadGVShuffle";
   int printOn                          = 0;
Mark Friedrichs's avatar
Mark Friedrichs committed
633
634
635

   // ---------------------------------------------------------------------------------------
	
Mark Friedrichs's avatar
Mark Friedrichs committed
636
637
	float* loadBuffer  = _getLoadBuffer();
   int* indices       = _getShuffleIndices( pmax );
Mark Friedrichs's avatar
Mark Friedrichs committed
638
639
	
   int rvSize = getRandomNumberStreamSize();
Mark Friedrichs's avatar
Mark Friedrichs committed
640
641
	for ( int ii = 0; ii < rvSize; ii++ ) {
		loadBuffer[ii] = (float) ii;
Mark Friedrichs's avatar
Mark Friedrichs committed
642
643
644
645
646
	}

	//How to come up with this number here?

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

Mark Friedrichs's avatar
Mark Friedrichs committed
649
		//for each p
Mark Friedrichs's avatar
Mark Friedrichs committed
650
651
652

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

Mark Friedrichs's avatar
Mark Friedrichs committed
653
654
655
			//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
656
657

			for ( int jj = 0; jj < p[ii]; jj++ ) {
Mark Friedrichs's avatar
Mark Friedrichs committed
658
				//indices[j] = (int) ( gmx_rando( &sdp->seed ) * sdp->gvSize );
Mark Friedrichs's avatar
Mark Friedrichs committed
659
				indices[jj] = (int) ( _generateGromacsRandomNumber( &seed )*rvSize );
Mark Friedrichs's avatar
Mark Friedrichs committed
660
			}
Mark Friedrichs's avatar
Mark Friedrichs committed
661
662
663
664
665
666

			// 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
667
			}
Mark Friedrichs's avatar
Mark Friedrichs committed
668
			loadBuffer[ indices[ p[ii]-1 ] ] = tmp;
Mark Friedrichs's avatar
Mark Friedrichs committed
669
670
		}
	}
Mark Friedrichs's avatar
Mark Friedrichs committed
671
   _getShuffleStream()->loadFromArray( loadBuffer );
Mark Friedrichs's avatar
Mark Friedrichs committed
672
	
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
   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
692
693
694
695
696
697
698
699
700
   return DefaultReturnValue;
}

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

Mark Friedrichs's avatar
Mark Friedrichs committed
701
int BrookRandomNumberGenerator::_shuffleGVStreams( void ){
Mark Friedrichs's avatar
Mark Friedrichs committed
702
703
704

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

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
705
   // static const std::string methodName = "\nBrookRandomNumberGenerator::_shuffleGVStreams";
Mark Friedrichs's avatar
Mark Friedrichs committed
706
707
708
709

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

Mark Friedrichs's avatar
Mark Friedrichs committed
711
	for( int ii = 0; ii < numberOfRvStreams - 1; ii++ ){
Mark Friedrichs's avatar
Mark Friedrichs committed
712
713
714
715
		kpermute_vectors( (float) getRandomNumberStreamWidth(),
                         _getShuffleStream()->getBrookStream(),
                         getRandomNumberStream( ii + 1 )->getBrookStream(), 
                         getRandomNumberStream( ii  )->getBrookStream() );
Mark Friedrichs's avatar
Mark Friedrichs committed
716
717
   }

Mark Friedrichs's avatar
Mark Friedrichs committed
718
719
720
721
	kpermute_vectors( (float) getRandomNumberStreamWidth(),
                      _getShuffleStream()->getBrookStream(),
                      getRandomNumberStream( 0 )->getBrookStream(), 
                      getRandomNumberStream( numberOfRvStreams - 1 )->getBrookStream() );
Mark Friedrichs's avatar
Mark Friedrichs committed
722
723
724
725

   return DefaultReturnValue;
}

Mark Friedrichs's avatar
Mark Friedrichs committed
726
/** 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
727
 * Advances the current position by 2*gpu->particles
Mark Friedrichs's avatar
Mark Friedrichs committed
728
729
730
 * 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
731
 * @param numberOfRandomValuesConsumedPerIteration    number of random values consumed/iteration
Mark Friedrichs's avatar
Mark Friedrichs committed
732
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
733
 * @return DefaultReturnValue;
Mark Friedrichs's avatar
Mark Friedrichs committed
734
735
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
736
int BrookRandomNumberGenerator::advanceGVCursor( int numberOfRandomValuesConsumedPerIteration ){
Mark Friedrichs's avatar
Mark Friedrichs committed
737
738
739

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

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
740
   static const std::string methodName = "BrookRandomNumberGenerator::advanceGVCursor";
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
741
   int printOn                         = 1;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
742
   FILE* log;
Mark Friedrichs's avatar
Mark Friedrichs committed
743
744

   // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
745
	
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
746
setLog( stderr );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
747
   if( printOn && getLog() ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
748
      log = getLog();
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
749
750
751
752
   } else {
      printOn = 0;
   }   

Mark Friedrichs's avatar
Mark Friedrichs committed
753
   int rvStreamSize  = getRandomNumberStreamSize();
Mark Friedrichs's avatar
Mark Friedrichs committed
754

Mark Friedrichs's avatar
Mark Friedrichs committed
755
   // use 2 random values per sd 
Mark Friedrichs's avatar
Mark Friedrichs committed
756

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
757
	_rvStreamOffset  += numberOfRandomValuesConsumedPerIteration;
Mark Friedrichs's avatar
Mark Friedrichs committed
758

Mark Friedrichs's avatar
Mark Friedrichs committed
759
	//Check if we've used up this texture
Mark Friedrichs's avatar
Mark Friedrichs committed
760

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
761
   char* action = "none";
Mark Friedrichs's avatar
Mark Friedrichs committed
762
	if ( _rvStreamOffset > rvStreamSize - numberOfRandomValuesConsumedPerIteration ){
Mark Friedrichs's avatar
Mark Friedrichs committed
763

Mark Friedrichs's avatar
Mark Friedrichs committed
764
		// next one if available
Mark Friedrichs's avatar
Mark Friedrichs committed
765

Mark Friedrichs's avatar
Mark Friedrichs committed
766
      _rvStreamOffset = 0;
Mark Friedrichs's avatar
Mark Friedrichs committed
767

Mark Friedrichs's avatar
Mark Friedrichs committed
768
769
		if ( _rvStreamIndex < _numberOfRandomNumberStreams - 1 ){
         _rvStreamIndex++;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
770
         action = "incremented stream index";
Mark Friedrichs's avatar
Mark Friedrichs committed
771
		} else {
Mark Friedrichs's avatar
Mark Friedrichs committed
772

Mark Friedrichs's avatar
Mark Friedrichs committed
773
			//No more textures, need to shuffle
Mark Friedrichs's avatar
Mark Friedrichs committed
774

Mark Friedrichs's avatar
Mark Friedrichs committed
775
			if( _numberOfShuffles < _maxShuffles ){
Mark Friedrichs's avatar
Mark Friedrichs committed
776

Mark Friedrichs's avatar
Mark Friedrichs committed
777
778
				_shuffleGVStreams( );
            _numberOfShuffles++;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
779
            action = "shuffled streams";
Mark Friedrichs's avatar
Mark Friedrichs committed
780

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

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
783
            if( _randomNumberGenerator == Mersenne ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
784
               action = "loaded new values to GPU using Mersenne rng";
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
785
786
               _loadRandomNumberStreamsMersenne( );
            } else if( _randomNumberGenerator == Kiss ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
787
               action = "loaded new values to GPU using KISS rng";
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
788
789
               _loadRandomNumberStreamsKiss( );
            } else if( _randomNumberGenerator == Original ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
790
               action = "loaded new values to GPU using original Gromac's rng";
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
791
               _loadGVStreamsOriginal( );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
792
793
794
            } else if( _randomNumberGenerator == FixedValue ){
               action = "loaded new fixed values to GPU";
               _loadRandomNumberStreamsFixedValue( );
Mark Friedrichs's avatar
Mark Friedrichs committed
795
796
797
798
799
            }
            _numberOfShuffles  = 0;
			}
         _rvStreamIndex = 0;
		}
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
800

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
801
      if( printOn ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
802
803
804
         (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
805
         (void) fflush( log );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
806
      }
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824

/*
      // 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
825
	}
Mark Friedrichs's avatar
Mark Friedrichs committed
826
827
828
829
830

   return DefaultReturnValue;
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
831
 * Get random number stream size
Mark Friedrichs's avatar
Mark Friedrichs committed
832
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
833
 * @return  random number stream size
Mark Friedrichs's avatar
Mark Friedrichs committed
834
835
836
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
837
838
int BrookRandomNumberGenerator::getRandomNumberStreamSize( void ) const {
   return _randomNumberStreamSize;
Mark Friedrichs's avatar
Mark Friedrichs committed
839
840
841
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
842
 * Get Shuffle stream 
Mark Friedrichs's avatar
Mark Friedrichs committed
843
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
844
 * @return  shuffle stream
Mark Friedrichs's avatar
Mark Friedrichs committed
845
846
847
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
848
849
BrookFloatStreamInternal* BrookRandomNumberGenerator::_getShuffleStream( void ) const {
   return _auxiliaryStreams[ShuffleStream];
Mark Friedrichs's avatar
Mark Friedrichs committed
850
851
852
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
853
 * Get random number stream 
Mark Friedrichs's avatar
Mark Friedrichs committed
854
 *
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
855
 * @param index     random number stream index     
Mark Friedrichs's avatar
Mark Friedrichs committed
856
857
 *
 * @return  random number stream
Mark Friedrichs's avatar
Mark Friedrichs committed
858
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
859
 * @throw OpenMMException if index is invalid or _randomNumberGeneratorStreams not set
Mark Friedrichs's avatar
Mark Friedrichs committed
860
861
862
 *
 */

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

Mark Friedrichs's avatar
Mark Friedrichs committed
865
   // ---------------------------------------------------------------------------------------
Mark Friedrichs's avatar
Mark Friedrichs committed
866

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

Mark Friedrichs's avatar
Mark Friedrichs committed
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
   // ---------------------------------------------------------------------------------------
	
   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
886
887
888
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
889
 * Get random value stream index
Mark Friedrichs's avatar
Mark Friedrichs committed
890
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
891
 * @return  random value stream index
Mark Friedrichs's avatar
Mark Friedrichs committed
892
893
894
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
895
896
int BrookRandomNumberGenerator::getRvStreamIndex( void ) const {
   return _rvStreamIndex;
Mark Friedrichs's avatar
Mark Friedrichs committed
897
898
899
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
900
 * Get random value stream offset
Mark Friedrichs's avatar
Mark Friedrichs committed
901
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
902
 * @return  random value stream offset
Mark Friedrichs's avatar
Mark Friedrichs committed
903
904
905
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
906
907
int BrookRandomNumberGenerator::getRvStreamOffset( void ) const {
   return _rvStreamOffset;
Mark Friedrichs's avatar
Mark Friedrichs committed
908
909
910
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
911
 * Get max shuffles
Mark Friedrichs's avatar
Mark Friedrichs committed
912
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
913
 * @return  max shuffles
Mark Friedrichs's avatar
Mark Friedrichs committed
914
915
916
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
917
918
int BrookRandomNumberGenerator::getMaxShuffles( void ) const {
   return _maxShuffles;
Mark Friedrichs's avatar
Mark Friedrichs committed
919
920
921
}

/** 
Mark Friedrichs's avatar
Mark Friedrichs committed
922
 * Get number of shuffles
Mark Friedrichs's avatar
Mark Friedrichs committed
923
 *
Mark Friedrichs's avatar
Mark Friedrichs committed
924
 * @return  number of shuffles
Mark Friedrichs's avatar
Mark Friedrichs committed
925
926
927
 *
 */

Mark Friedrichs's avatar
Mark Friedrichs committed
928
929
int BrookRandomNumberGenerator::_getNumberOfShuffles( void ) const {
   return _numberOfShuffles;
Mark Friedrichs's avatar
Mark Friedrichs committed
930
931
932
933
934
}

/** 
 * Initialize stream dimensions
 * 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
935
 * @param numberOfParticles         number of particles
Mark Friedrichs's avatar
Mark Friedrichs committed
936
937
938
939
940
941
 * @param platform                  platform
 *      
 * @return ErrorReturnValue if error, else DefaultReturnValue
 *
 */

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
942
int BrookRandomNumberGenerator::_initializeStreamSizes( int numberOfParticles, const Platform& platform ){
Mark Friedrichs's avatar
Mark Friedrichs committed
943
944
945
946
947
948
949

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

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

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

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
950
   setNumberOfParticles( numberOfParticles );
Mark Friedrichs's avatar
Mark Friedrichs committed
951

Mark Friedrichs's avatar
Mark Friedrichs committed
952
   // get randomNumber stream dimensions
Mark Friedrichs's avatar
Mark Friedrichs committed
953

Mark Friedrichs's avatar
Mark Friedrichs committed
954
955
956
957
958
   const BrookPlatform brookPlatform            = dynamic_cast<const BrookPlatform&> (platform);
   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
959

Mark Friedrichs's avatar
Mark Friedrichs committed
960
   _randomNumberStreamSize                      =  _randomNumberStreamWidth*_randomNumberStreamHeight;
Mark Friedrichs's avatar
Mark Friedrichs committed
961

Mark Friedrichs's avatar
Mark Friedrichs committed
962
   return DefaultReturnValue;
Mark Friedrichs's avatar
Mark Friedrichs committed
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
}

/** 
 * 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
980
   float dangleValue                          = 0.0f;
Mark Friedrichs's avatar
Mark Friedrichs committed
981
982
983

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

Mark Friedrichs's avatar
Mark Friedrichs committed
984
985
   int randomNumberStreamSize                = getRandomNumberStreamSize();
   int randomNumberStreamWidth               = getRandomNumberStreamWidth();
Mark Friedrichs's avatar
Mark Friedrichs committed
986

Mark Friedrichs's avatar
Mark Friedrichs committed
987
988
989
   _auxiliaryStreams[ShuffleStream]         = new BrookFloatStreamInternal( BrookCommon::ShuffleStream,
                                                                            randomNumberStreamSize, randomNumberStreamWidth,
                                                                            BrookStreamInternal::Float, dangleValue );
Mark Friedrichs's avatar
Mark Friedrichs committed
990

Mark Friedrichs's avatar
Mark Friedrichs committed
991
992
993
994
995
   // 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
996
997
   }

Mark Friedrichs's avatar
Mark Friedrichs committed
998
999
   if( _randomNumberGeneratorStreams ){
      delete[] _randomNumberGeneratorStreams;
Mark Friedrichs's avatar
Mark Friedrichs committed
1000
1001
   }

Mark Friedrichs's avatar
Mark Friedrichs committed
1002
   _randomNumberGeneratorStreams = new BrookFloatStreamInternal*[_numberOfRandomNumberStreams];
Mark Friedrichs's avatar
Mark Friedrichs committed
1003

Mark Friedrichs's avatar
Mark Friedrichs committed
1004
   for( int ii = 0; ii < _numberOfRandomNumberStreams; ii++ ){
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1005
      _randomNumberGeneratorStreams[ii]    = new BrookFloatStreamInternal( BrookCommon::RandomValuesStream,
Mark Friedrichs's avatar
Mark Friedrichs committed
1006
1007
                                                                           randomNumberStreamSize, randomNumberStreamWidth,
                                                                           BrookStreamInternal::Float3, dangleValue );
Mark Friedrichs's avatar
Mark Friedrichs committed
1008
1009
   }

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1010
1011
1012
1013
   // create shuffle stream

   _loadGVShuffle();

Mark Friedrichs's avatar
Mark Friedrichs committed
1014
1015
1016
1017
   return DefaultReturnValue;
}

/*  
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1018
 * Setup of streams, ... associated w/ random number generator
Mark Friedrichs's avatar
Mark Friedrichs committed
1019
 *
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1020
 * @param numberOfParticles     number of particles
Mark Friedrichs's avatar
Mark Friedrichs committed
1021
1022
1023
1024
1025
1026
 * @param platform              Brook platform
 *
 * @return nonzero value if error
 *
 * */
    
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1027
int BrookRandomNumberGenerator::setup( int numberOfParticles, const Platform& platform ){
Mark Friedrichs's avatar
Mark Friedrichs committed
1028
1029
1030
    
// ---------------------------------------------------------------------------------------

Mark Friedrichs's avatar
Mark Friedrichs committed
1031
   //static const std::string methodName      = "BrookRandomNumberGenerator::setup";
Mark Friedrichs's avatar
Mark Friedrichs committed
1032
1033
1034

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

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1035
1036
1037
   const BrookPlatform brookPlatform            = dynamic_cast<const BrookPlatform&> (platform);
   setLog( brookPlatform.getLog() );

Mark Friedrichs's avatar
Mark Friedrichs committed
1038
1039
   // set stream sizes and then create streams

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1040
   _initializeStreamSizes( numberOfParticles, platform );
Mark Friedrichs's avatar
Mark Friedrichs committed
1041
   _initializeStreams( platform );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1042
1043
1044
1045

   if( _randomNumberGenerator == Mersenne ){
      _loadRandomNumberStreamsMersenne( );
   } else if( _randomNumberGenerator == Kiss ){
1046
      _loadRandomNumberStreamsKiss( );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1047
1048
   } else if( _randomNumberGenerator == Original ){
      _loadGVStreamsOriginal( );
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1049
1050
   } else if( _randomNumberGenerator == FixedValue ){
      _loadRandomNumberStreamsFixedValue( );
1051
   }
Mark Friedrichs's avatar
Mark Friedrichs committed
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086

   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
1087
1088
1089
1090
1091
1092
1093
   
   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
1094
1095
   } else if( _randomNumberGenerator == FixedValue ){
      (void) LOCAL_SPRINTF( value, "%s", "Fixed value Rng");
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1096
   }
Mark Friedrichs's avatar
Mark Friedrichs committed
1097
   message << _getLine( tab, "Random number generator:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1098

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

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

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

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

Mark Friedrichs's avatar
Mark Friedrichs committed
1111
1112
   (void) LOCAL_SPRINTF( value, "%ld", getRandomNumberSeed() );
   message << _getLine( tab, "Seed:", value ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1113

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

Mark Friedrichs's avatar
Mark Friedrichs committed
1117
1118
1119
1120
1121
   (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
1122

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

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

Mark Friedrichs's avatar
Mark Friedrichs committed
1128
   message << _getLine( tab, "Shuffle:",              (_getShuffleStream()      ? Set : NotSet) ); 
Mark Friedrichs's avatar
Mark Friedrichs committed
1129
 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1130
1131
   // show stats

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1132
1133
1134
1135
1136
1137
1138
1139
1140
   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
1141
1142
1143
1144
1145
1146
   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
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
#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

1182
1183
1184
   message << "\n   Statistics:\n";

   double statistics[7];
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1185
   getStatistics( statistics, -1, cumulativeStatistics );
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204

   (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
1205
1206
1207
1208
   if( cumulativeStatistics[4] > 1000.0 ){ 

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

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
      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
1236
1237
1238
1239
#undef LOCAL_SPRINTF

   return message.str();
}
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258

/* 
 * 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
1259
int BrookRandomNumberGenerator::getStatistics( double statistics[7], int streamIndex, double cumulativeStatistics[7] ) const {
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277

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

   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
1278
         int numberOfValues     = _randomNumberGeneratorStreams[ii]->getStreamSize()*_randomNumberGeneratorStreams[ii]->getWidth();
1279
         const float* dataArray = (float*) dataArrayV;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1280
         for( int ii = 0; ii < numberOfValues; ii++ ){
1281

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

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1284
             double rv2     = dataArray[ii]*dataArray[ii];
1285
             statistics[1] += rv2;
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1286
             statistics[2] += rv2*dataArray[ii];
1287
1288
             statistics[3] += rv2*rv2;

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1289
1290
             if( statistics[5] > dataArray[ii] ){
                statistics[5] = dataArray[ii];
1291
             }
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1292
1293
             if( statistics[6] < dataArray[ii] ){
                statistics[6] = dataArray[ii];
1294
1295
             }
         } 
Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1296
         statistics[4] += (double) numberOfValues;
1297
1298
1299
      }   
   }   

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1300
1301
   // accumulate moments, ... in cumulativeStatistics array

Mark Friedrichs's avatar
Mods  
Mark Friedrichs committed
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
   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];
   }

1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
   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;
}