diff --git a/README.md b/README.md index 658ac0a..4185179 100644 --- a/README.md +++ b/README.md @@ -128,7 +128,7 @@ In squiggle.c, this ambiguity doesn't exist, at the cost of much greater overhea int main(){ // set randomness seed - uint32_t* seed = malloc(sizeof(uint32_t)); + uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with a seed of 0 float a = sample_to(1, 10, seed); @@ -153,7 +153,7 @@ vs #include #include -float draw_xyz(uint32_t* seed){ +float draw_xyz(uint64_t* seed){ // function could also be placed inside main with gcc nested functions extension. return sample_to(1, 20, seed); } @@ -161,7 +161,7 @@ float draw_xyz(uint32_t* seed){ int main(){ // set randomness seed - uint32_t* seed = malloc(sizeof(uint32_t)); + uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with a seed of 0 float a = draw_xyz(seed); diff --git a/examples/01_one_sample/example b/examples/01_one_sample/example index ce80991..5038416 100755 Binary files a/examples/01_one_sample/example and b/examples/01_one_sample/example differ diff --git a/examples/01_one_sample/example.c b/examples/01_one_sample/example.c index 4badc5e..f1687ae 100644 --- a/examples/01_one_sample/example.c +++ b/examples/01_one_sample/example.c @@ -4,29 +4,29 @@ #include // Estimate functions -float sample_0(uint32_t* seed) +float sample_0(uint64_t* seed) { return 0; } -float sample_1(uint32_t* seed) +float sample_1(uint64_t* seed) { return 1; } -float sample_few(uint32_t* seed) +float sample_few(uint64_t* seed) { return sample_to(1, 3, seed); } -float sample_many(uint32_t* seed) +float sample_many(uint64_t* seed) { return sample_to(2, 10, seed); } int main(){ // set randomness seed - uint32_t* seed = malloc(sizeof(uint32_t)); + uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with 0 float p_a = 0.8; @@ -35,7 +35,7 @@ int main(){ int n_dists = 4; float weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 }; - float (*samplers[])(uint32_t*) = { sample_0, sample_1, sample_few, sample_many }; + float (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many }; float result_one = sample_mixture(samplers, weights, n_dists, seed); printf("result_one: %f\n", result_one); diff --git a/examples/02_many_samples/example b/examples/02_many_samples/example index 6bc060e..4b781cf 100755 Binary files a/examples/02_many_samples/example and b/examples/02_many_samples/example differ diff --git a/examples/02_many_samples/example.c b/examples/02_many_samples/example.c index c77a17c..7f3faa0 100644 --- a/examples/02_many_samples/example.c +++ b/examples/02_many_samples/example.c @@ -4,29 +4,29 @@ #include "../../squiggle.h" // Estimate functions -float sample_0(uint32_t* seed) +float sample_0(uint64_t* seed) { return 0; } -float sample_1(uint32_t* seed) +float sample_1(uint64_t* seed) { return 1; } -float sample_few(uint32_t* seed) +float sample_few(uint64_t* seed) { return sample_to(1, 3, seed); } -float sample_many(uint32_t* seed) +float sample_many(uint64_t* seed) { return sample_to(2, 10, seed); } int main(){ // set randomness seed - uint32_t* seed = malloc(sizeof(uint32_t)); + uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with 0 float p_a = 0.8; @@ -35,7 +35,7 @@ int main(){ int n_dists = 4; float weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 }; - float (*samplers[])(uint32_t*) = { sample_0, sample_1, sample_few, sample_many }; + float (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many }; int n_samples = 1000000; float* result_many = (float *) malloc(n_samples * sizeof(float)); diff --git a/examples/03_gcc_nested_function/example b/examples/03_gcc_nested_function/example index ef1f7ac..7c4bfc3 100755 Binary files a/examples/03_gcc_nested_function/example and b/examples/03_gcc_nested_function/example differ diff --git a/examples/03_gcc_nested_function/example.c b/examples/03_gcc_nested_function/example.c index 7d1a0b5..62c8e53 100644 --- a/examples/03_gcc_nested_function/example.c +++ b/examples/03_gcc_nested_function/example.c @@ -5,7 +5,7 @@ int main(){ // set randomness seed - uint32_t* seed = malloc(sizeof(uint32_t)); + uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with 0 float p_a = 0.8; @@ -14,12 +14,12 @@ int main(){ int n_dists = 4; - float sample_0(uint32_t* seed){ return 0; } - float sample_1(uint32_t* seed) { return 1; } - float sample_few(uint32_t* seed){ return sample_to(1, 3, seed); } - float sample_many(uint32_t* seed){ return sample_to(2, 10, seed); } + float sample_0(uint64_t* seed){ return 0; } + float sample_1(uint64_t* seed) { return 1; } + float sample_few(uint64_t* seed){ return sample_to(1, 3, seed); } + float sample_many(uint64_t* seed){ return sample_to(2, 10, seed); } - float (*samplers[])(uint32_t*) = { sample_0, sample_1, sample_few, sample_many }; + float (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many }; float weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 }; int n_samples = 1000000; diff --git a/examples/04_sample_from_cdf_simple/example b/examples/04_sample_from_cdf_simple/example index 97dc565..84ed2c2 100755 Binary files a/examples/04_sample_from_cdf_simple/example and b/examples/04_sample_from_cdf_simple/example differ diff --git a/examples/04_sample_from_cdf_simple/example.c b/examples/04_sample_from_cdf_simple/example.c index 8f20a4b..0d21ec5 100644 --- a/examples/04_sample_from_cdf_simple/example.c +++ b/examples/04_sample_from_cdf_simple/example.c @@ -49,7 +49,7 @@ void test_inverse_cdf_float(char* cdf_name, float cdf_float(float)) } } -void test_and_time_sampler_float(char* cdf_name, float cdf_float(float), uint32_t* seed) +void test_and_time_sampler_float(char* cdf_name, float cdf_float(float), uint64_t* seed) { printf("\nGetting some samples from %s:\n", cdf_name); clock_t begin = clock(); @@ -75,7 +75,7 @@ int main() // Testing samplers // set randomness seed - uint32_t* seed = malloc(sizeof(uint32_t)); + uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with 0 // Test float sampler diff --git a/examples/05_sample_from_cdf_beta/example b/examples/05_sample_from_cdf_beta/example index 7f3c14e..853558f 100755 Binary files a/examples/05_sample_from_cdf_beta/example and b/examples/05_sample_from_cdf_beta/example differ diff --git a/examples/05_sample_from_cdf_beta/example.c b/examples/05_sample_from_cdf_beta/example.c index ce023ca..1328ba2 100644 --- a/examples/05_sample_from_cdf_beta/example.c +++ b/examples/05_sample_from_cdf_beta/example.c @@ -131,7 +131,7 @@ void test_inverse_cdf_box(char* cdf_name, struct box cdf_box(float)) } } -void test_and_time_sampler_box(char* cdf_name, struct box cdf_box(float), uint32_t* seed) +void test_and_time_sampler_box(char* cdf_name, struct box cdf_box(float), uint64_t* seed) { printf("\nGetting some samples from %s:\n", cdf_name); clock_t begin = clock(); @@ -154,7 +154,7 @@ int main() test_inverse_cdf_box("cdf_beta", cdf_beta); // Test box sampler - uint32_t* seed = malloc(sizeof(uint32_t)); + uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with 0 test_and_time_sampler_box("cdf_beta", cdf_beta, seed); // Ok, this is slower than python!! diff --git a/examples/06_gamma_beta/example b/examples/06_gamma_beta/example index a7e08b3..b3743d3 100755 Binary files a/examples/06_gamma_beta/example and b/examples/06_gamma_beta/example differ diff --git a/examples/06_gamma_beta/example.c b/examples/06_gamma_beta/example.c index 9d5335e..8c11342 100644 --- a/examples/06_gamma_beta/example.c +++ b/examples/06_gamma_beta/example.c @@ -8,7 +8,7 @@ int main() { // set randomness seed - uint32_t* seed = malloc(sizeof(uint32_t)); + uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with 0 int n = 1000 * 1000; diff --git a/squiggle.c b/squiggle.c index 9b8d0df..9d649d3 100644 --- a/squiggle.c +++ b/squiggle.c @@ -14,14 +14,14 @@ const float PI = 3.14159265358979323846; // M_PI in gcc gnu99 // Pseudo Random number generator -uint32_t xorshift32(uint32_t* seed) +uint64_t xorshift32(uint32_t* seed) { // Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" - // See + // See // https://en.wikipedia.org/wiki/Xorshift // Also some drama: , - uint32_t x = *seed; + uint64_t x = *seed; x ^= x << 13; x ^= x >> 17; x ^= x << 5; @@ -31,7 +31,7 @@ uint32_t xorshift32(uint32_t* seed) uint64_t xorshift64(uint64_t* seed) { // Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" - // See + // See // https://en.wikipedia.org/wiki/Xorshift // Also some drama: , @@ -44,13 +44,13 @@ uint64_t xorshift64(uint64_t* seed) // Distribution & sampling functions // Unit distributions -float sample_unit_uniform(uint32_t* seed) +float sample_unit_uniform(uint64_t* seed) { // samples uniform from [0,1] interval. - return ((float)xorshift32(seed)) / ((float)UINT32_MAX); + return ((float)xorshift64(seed)) / ((float)UINT64_MAX); } -float sample_unit_normal(uint32_t* seed) +float sample_unit_normal(uint64_t* seed) { // See: float u1 = sample_unit_uniform(seed); @@ -60,22 +60,22 @@ float sample_unit_normal(uint32_t* seed) } // Composite distributions -float sample_uniform(float start, float end, uint32_t* seed) +float sample_uniform(float start, float end, uint64_t* seed) { return sample_unit_uniform(seed) * (end - start) + start; } -float sample_normal(float mean, float sigma, uint32_t* seed) +float sample_normal(float mean, float sigma, uint64_t* seed) { return (mean + sigma * sample_unit_normal(seed)); } -float sample_lognormal(float logmean, float logsigma, uint32_t* seed) +float sample_lognormal(float logmean, float logsigma, uint64_t* seed) { return expf(sample_normal(logmean, logsigma, seed)); } -float sample_to(float low, float high, uint32_t* seed) +float sample_to(float low, float high, uint64_t* seed) { // Given a (positive) 90% confidence interval, // returns a sample from a lognormal @@ -88,7 +88,7 @@ float sample_to(float low, float high, uint32_t* seed) return sample_lognormal(logmean, logsigma, seed); } -float sample_gamma(float alpha, uint32_t* seed) +float sample_gamma(float alpha, uint64_t* seed) { // A Simple Method for Generating Gamma Variables, Marsaglia and Wan Tsang, 2001 @@ -125,7 +125,7 @@ float sample_gamma(float alpha, uint32_t* seed) } } -float sample_beta(float a, float b, uint32_t* seed) +float sample_beta(float a, float b, uint64_t* seed) { float gamma_a = sample_gamma(a, seed); float gamma_b = sample_gamma(b, seed); @@ -168,7 +168,7 @@ float array_std(float* array, int length) } // Mixture function -float sample_mixture(float (*samplers[])(uint32_t*), float* weights, int n_dists, uint32_t* seed) +float sample_mixture(float (*samplers[])(uint64_t*), float* weights, int n_dists, uint64_t* seed) { // You can see a simpler version of this function in the git history // or in C-02-better-algorithm-one-thread/ @@ -364,13 +364,13 @@ struct box inverse_cdf_box(struct box cdf_box(float), float p) } // Sampler based on inverse cdf and randomness function -struct box sampler_cdf_box(struct box cdf(float), uint32_t* seed) +struct box sampler_cdf_box(struct box cdf(float), uint64_t* seed) { float p = sample_unit_uniform(seed); struct box result = inverse_cdf_box(cdf, p); return result; } -struct box sampler_cdf_float(float cdf(float), uint32_t* seed) +struct box sampler_cdf_float(float cdf(float), uint64_t* seed) { float p = sample_unit_uniform(seed); struct box result = inverse_cdf_float(cdf, p); @@ -378,7 +378,7 @@ struct box sampler_cdf_float(float cdf(float), uint32_t* seed) } /* Could also define other variations, e.g., -float sampler_danger(struct box cdf(float), uint32_t* seed) +float sampler_danger(struct box cdf(float), uint64_t* seed) { float p = sample_unit_uniform(seed); struct box result = inverse_cdf_box(cdf, p); diff --git a/squiggle.h b/squiggle.h index a94158f..aa3399b 100644 --- a/squiggle.h +++ b/squiggle.h @@ -1,24 +1,24 @@ #ifndef SQUIGGLEC #define SQUIGGLEC -// uint32_t header +// uint64_t header #include // Pseudo Random number generator -uint32_t xorshift32(uint32_t* seed); +uint64_t xorshift64(uint64_t* seed); // Basic distribution sampling functions -float sample_unit_uniform(uint32_t* seed); -float sample_unit_normal(uint32_t* seed); +float sample_unit_uniform(uint64_t* seed); +float sample_unit_normal(uint64_t* seed); // Composite distribution sampling functions -float sample_uniform(float start, float end, uint32_t* seed); -float sample_normal(float mean, float sigma, uint32_t* seed); -float sample_lognormal(float logmean, float logsigma, uint32_t* seed); -float sample_to(float low, float high, uint32_t* seed); +float sample_uniform(float start, float end, uint64_t* seed); +float sample_normal(float mean, float sigma, uint64_t* seed); +float sample_lognormal(float logmean, float logsigma, uint64_t* seed); +float sample_to(float low, float high, uint64_t* seed); -float sample_gamma(float alpha, uint32_t* seed); -float sample_beta(float a, float b, uint32_t* seed); +float sample_gamma(float alpha, uint64_t* seed); +float sample_beta(float a, float b, uint64_t* seed); // Array helpers float array_sum(float* array, int length); @@ -27,7 +27,7 @@ float array_mean(float* array, int length); float array_std(float* array, int length); // Mixture function -float sample_mixture(float (*samplers[])(uint32_t*), float* weights, int n_dists, uint32_t* seed); +float sample_mixture(float (*samplers[])(uint64_t*), float* weights, int n_dists, uint64_t* seed); // Box struct box { @@ -47,7 +47,7 @@ struct box inverse_cdf_float(float cdf(float), float p); struct box inverse_cdf_box(struct box cdf_box(float), float p); // Samplers from cdf -struct box sampler_cdf_float(float cdf(float), uint32_t* seed); -struct box sampler_cdf_box(struct box cdf(float), uint32_t* seed); +struct box sampler_cdf_float(float cdf(float), uint64_t* seed); +struct box sampler_cdf_box(struct box cdf(float), uint64_t* seed); #endif diff --git a/test/test b/test/test index 932cb62..708e1ff 100755 Binary files a/test/test and b/test/test differ diff --git a/test/test.c b/test/test.c index 7c5943d..7e472f6 100644 --- a/test/test.c +++ b/test/test.c @@ -4,9 +4,9 @@ #include #include -#define N 10 * 1000 * 1000 +#define N 1000 * 1000 -void test_unit_uniform(uint32_t* seed){ +void test_unit_uniform(uint64_t* seed){ float* unit_uniform_array = malloc(sizeof(float) * N); for(int i=0; i