-
Notifications
You must be signed in to change notification settings - Fork 32
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
random generators compatibility with parallel processing (#513)
* new random generator * merged from master * send process ident * Fortran ordering in r_in * sync with patpass * restore test_patpass * python common and thread-specific RNGs * C common and thread-specific RNGs * merged from ;aster * Initializers in atrandom.c * removed the id field * taken from master * "inc" initializer must be odd
- Loading branch information
Showing
13 changed files
with
521 additions
and
237 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,24 +1,111 @@ | ||
| /* | ||
| * PCG Random Number Generation for C. | ||
| * | ||
| * Copyright 2014 Melissa O'Neill <[email protected]> | ||
| * | ||
| * Licensed under the Apache License, Version 2.0 (the "License"); | ||
| * you may not use this file except in compliance with the License. | ||
| * You may obtain a copy of the License at | ||
| * | ||
| * http://www.apache.org/licenses/LICENSE-2.0 | ||
| * | ||
| * Unless required by applicable law or agreed to in writing, software | ||
| * distributed under the License is distributed on an "AS IS" BASIS, | ||
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. | ||
| * See the License for the specific language governing permissions and | ||
| * limitations under the License. | ||
| * | ||
| * For additional information about the PCG random number generation scheme, | ||
| * including its license and other licensing options, visit | ||
| * | ||
| * http://www.pcg-random.org | ||
| */ | ||
| #include <math.h> | ||
| #include <stdint.h> | ||
|
|
||
| void init_seed(int seed){ | ||
| static int initseed=1; | ||
| if(initseed) | ||
| { | ||
| srand(seed); | ||
| initseed=0; | ||
| } | ||
| #ifndef TWOPI | ||
| #define TWOPI 6.28318530717959 | ||
| #endif /*TWOPI*/ | ||
|
|
||
| /* initial RNG definitions */ | ||
| #define AT_RNG_STATE 0x853c49e6748fea9bULL | ||
| #define AT_RNG_INC 0xda3e39cb94b95bdbULL | ||
|
|
||
| #define COMMON_PCG32_INITIALIZER { AT_RNG_STATE, AT_RNG_INC } | ||
| #define THREAD_PCG32_INITIALIZER { AT_RNG_STATE, 1ULL } | ||
|
|
||
| struct pcg_state_setseq_64 { // Internals are *Private*. | ||
| uint64_t state; // RNG state. All values are possible. | ||
| uint64_t inc; // Controls which RNG sequence (stream) is | ||
| // selected. Must *always* be odd. | ||
| }; | ||
| typedef struct pcg_state_setseq_64 pcg32_random_t; | ||
|
|
||
| // If you *must* statically initialize it, here's one. | ||
|
|
||
| #define PCG32_INITIALIZER { AT_RNG_STATE, AT_RNG_INC } | ||
|
|
||
| static pcg32_random_t pcg32_global = PCG32_INITIALIZER; | ||
|
|
||
| // pcg32_random() | ||
| // pcg32_random_r(rng) | ||
| // Generate a uniformly distributed 32-bit random number | ||
|
|
||
| static uint32_t pcg32_random_r(pcg32_random_t* rng) | ||
| { | ||
| if (!rng) rng = &pcg32_global; | ||
| uint64_t oldstate; | ||
| #pragma omp atomic read | ||
| oldstate = rng->state; | ||
| #pragma omp atomic write | ||
| rng->state = oldstate * 6364136223846793005ULL + rng->inc; | ||
| uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u; | ||
| uint32_t rot = oldstate >> 59u; | ||
| return (xorshifted >> rot) | (xorshifted << ((-rot) & 31)); | ||
| } | ||
|
|
||
| double atdrand(void) /* uniform distribution, (0..1] */ | ||
| static uint32_t pcg32_random() | ||
| { | ||
| return pcg32_random_r(&pcg32_global); | ||
| } | ||
|
|
||
| // pcg32_srandom(initstate, initseq) | ||
| // pcg32_srandom_r(rng, initstate, initseq): | ||
| // Seed the rng. Specified in two parts, state initializer and a | ||
| // sequence selection constant (a.k.a. stream id) | ||
|
|
||
| static void pcg32_srandom_r(pcg32_random_t* rng, uint64_t initstate, uint64_t initseq) | ||
| { | ||
| return (rand()+1.0)/(RAND_MAX+1.0); | ||
| if (!rng) rng = &pcg32_global; | ||
| rng->state = 0U; | ||
| rng->inc = (initseq << 1u) | 1u; | ||
| pcg32_random_r(rng); | ||
| rng->state += initstate; | ||
| pcg32_random_r(rng); | ||
| } | ||
|
|
||
| double atrandn() /* gaussian distribution, sigma=1, mean=0 */ | ||
| static void pcg32_srandom(uint64_t seed, uint64_t seq) | ||
| { | ||
| pcg32_srandom_r(&pcg32_global, seed, seq); | ||
| } | ||
|
|
||
| /* Functions for uniform, normal and Poisson distributions with | ||
| external state variable */ | ||
|
|
||
| static double atrandd_r(pcg32_random_t* rng) | ||
| /* Uniform [0, 1] distribution */ | ||
| { | ||
| return ldexp(pcg32_random_r(rng), -32); | ||
| } | ||
|
|
||
| static double atrandn_r(pcg32_random_t* rng) | ||
| /* gaussian distribution, sigma=1, mean=0 */ | ||
| { | ||
| /* Marsaglia polar method: https://en.wikipedia.org/wiki/Marsaglia_polar_method */ | ||
|
|
||
| static bool hasSpare = false; | ||
| static double spare; | ||
| static double u, v, s; | ||
| double u, v, s; | ||
|
|
||
| if (hasSpare) { | ||
| hasSpare = false; | ||
|
|
@@ -27,46 +114,52 @@ double atrandn() /* gaussian distribution, sigma=1, mean=0 */ | |
|
|
||
| hasSpare = true; | ||
| do { | ||
| u = (rand() / (double) RAND_MAX) * 2.0 - 1.0; | ||
| v = (rand() / (double) RAND_MAX) * 2.0 - 1.0; | ||
| u = 2.0 * atrandd_r(rng) - 1.0; | ||
| v = 2.0 * atrandd_r(rng) - 1.0; | ||
| s = u * u + v * v; | ||
| } | ||
| while( (s >= 1.0) || (s == 0.0) ); | ||
| while ((s >= 1.0) || (s == 0.0)); | ||
| s = sqrt(-2.0 * log(s) / s); | ||
| spare = v * s; | ||
| return u * s; | ||
| } | ||
|
|
||
| int atrandp(double lamb) /* poisson distribution */ | ||
| static int atrandp_r(pcg32_random_t* rng, double lamb) | ||
| /* poisson distribution */ | ||
| { | ||
|
|
||
| int pk,k; | ||
| double l,p,u,u1,u2,twopi; | ||
|
|
||
| l = -lamb; | ||
| k = 0; | ||
| p = 0.0; | ||
| twopi = 4.0*asin(1.0); | ||
| pk=0; | ||
|
|
||
| if(lamb<11){ | ||
|
|
||
| do{ | ||
| k = k + 1; | ||
| u = atdrand(); | ||
| p = p + log(u); | ||
| }while(p>l); | ||
|
|
||
| int pk; | ||
|
|
||
| if (lamb<11) { | ||
| double l = -lamb; | ||
| int k = 0; | ||
| double p = 0; | ||
| do { | ||
| k += 1; | ||
| p += log(atrandd_r(rng)); | ||
| } while (p>l); | ||
| pk = k-1; | ||
|
|
||
| }else{ | ||
| u1 = atdrand(); | ||
| u2 = atdrand(); | ||
| if(u1==0.0){ | ||
| u1 = 1e-18; | ||
| }; | ||
| pk = (int)floor((sqrt(-2.0*log(u1))*cos(twopi*u2))*sqrt(lamb)+lamb); | ||
| } | ||
|
|
||
| return pk; | ||
| else { /* Gaussian approximation */ | ||
| pk = (int)floor(atrandn_r(rng)*sqrt(lamb)+lamb); | ||
| } | ||
|
|
||
| return pk; | ||
| } | ||
|
|
||
| /* Functions for uniform, normal and Poisson distributions with | ||
| internal state variable */ | ||
|
|
||
| static inline double atrandd(void) | ||
| { | ||
| return atrandd_r(&pcg32_global); | ||
| } | ||
|
|
||
| static inline double atrandn(void) | ||
| { | ||
| return atrandn_r(&pcg32_global); | ||
| } | ||
|
|
||
| static inline int atrandp(double lamb) | ||
| { | ||
| return atrandp_r(&pcg32_global, lamb); | ||
| } | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -7,7 +7,8 @@ at.tracking.atpass | |
|
|
||
| .. autosummary:: | ||
|
|
||
| isopenmp | ||
| ismpi | ||
| atpass | ||
| elempass | ||
| elempass | ||
| reset_rng | ||
| common_rng | ||
| thread_rng | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.