squiggle.c

Self-contained Monte Carlo estimation in C99
Log | Files | Refs | README

commit 1d89eb6231aa02415e2d8e906d6ee929d54e7366
parent 199e76bdfb061bed4ec0ed1705c2ec0f6780fa41
Author: NunoSempere <nuno.sempere@protonmail.com>
Date:   Sat, 20 Jan 2024 14:30:20 +0100

formatting pass, upkeep

Diffstat:
Mexamples/core/03_gcc_nested_function/example.c | 12++++++++++--
Aexamples/core/06_dissolving_fermi_paradox/example | 0
Aexamples/core/06_dissolving_fermi_paradox/example.c | 147+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Dexamples/core/06_dissolving_fermi_paradox/scratchpad.c | 157-------------------------------------------------------------------------------
Mscratchpad/scratchpad | 0
Mscratchpad/scratchpad.c | 57+++++++++++++++++++++++++++++++++------------------------
6 files changed, 190 insertions(+), 183 deletions(-)

diff --git a/examples/core/03_gcc_nested_function/example.c b/examples/core/03_gcc_nested_function/example.c @@ -15,8 +15,16 @@ int main() int n_dists = 4; // These are nested functions. They will not compile without gcc. - double sample_0(uint64_t * seed) { UNUSED(seed); return 0; } - double sample_1(uint64_t * seed) { UNUSED(seed); return 1; } + double sample_0(uint64_t * seed) + { + UNUSED(seed); + return 0; + } + double sample_1(uint64_t * seed) + { + UNUSED(seed); + return 1; + } double sample_few(uint64_t * seed) { return sample_to(1, 3, seed); } double sample_many(uint64_t * seed) { return sample_to(2, 10, seed); } diff --git a/examples/core/06_dissolving_fermi_paradox/example b/examples/core/06_dissolving_fermi_paradox/example Binary files differ. diff --git a/examples/core/06_dissolving_fermi_paradox/example.c b/examples/core/06_dissolving_fermi_paradox/example.c @@ -0,0 +1,147 @@ +#include "../../../squiggle.h" +#include <math.h> +#include <stdint.h> +#include <stdio.h> +#include <stdlib.h> + +double sample_loguniform(double a, double b, uint64_t* seed) +{ + return exp(sample_uniform(log(a), log(b), seed)); +} + +int main() +{ + // Replicate <https://arxiv.org/pdf/1806.02404.pdf>, and in particular the red line in page 11. + // Could also be interesting to just produce and save many samples. + + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = UINT64_MAX / 64; // xorshift can't start with a seed of 0 + + double sample_fermi_naive(uint64_t * seed) + { + double rate_of_star_formation = sample_loguniform(1, 100, seed); + double fraction_of_stars_with_planets = sample_loguniform(0.1, 1, seed); + double number_of_habitable_planets_per_star_system = sample_loguniform(0.1, 1, seed); + double rate_of_life_formation_in_habitable_planets = sample_lognormal(1, 50, seed); + double fraction_of_habitable_planets_in_which_any_life_appears = -expm1(-rate_of_life_formation_in_habitable_planets); + // double fraction_of_habitable_planets_in_which_any_life_appears = 1-exp(-rate_of_life_formation_in_habitable_planets); + // but with more precision + double fraction_of_planets_with_life_in_which_intelligent_life_appears = sample_loguniform(0.001, 1, seed); + double fraction_of_intelligent_planets_which_are_detectable_as_such = sample_loguniform(0.01, 1, seed); + double longevity_of_detectable_civilizations = sample_loguniform(100, 10000000000, seed); + + // printf(" rate_of_star_formation = %lf\n", rate_of_star_formation); + // printf(" fraction_of_stars_with_planets = %lf\n", fraction_of_stars_with_planets); + // printf(" number_of_habitable_planets_per_star_system = %lf\n", number_of_habitable_planets_per_star_system); + // printf(" rate_of_life_formation_in_habitable_planets = %.16lf\n", rate_of_life_formation_in_habitable_planets); + // printf(" fraction_of_habitable_planets_in_which_any_life_appears = %lf\n", fraction_of_habitable_planets_in_which_any_life_appears); + // printf(" fraction_of_planets_with_life_in_which_intelligent_life_appears = %lf\n", fraction_of_planets_with_life_in_which_intelligent_life_appears); + // printf(" fraction_of_intelligent_planets_which_are_detectable_as_such = %lf\n", fraction_of_intelligent_planets_which_are_detectable_as_such); + // printf(" longevity_of_detectable_civilizations = %lf\n", longevity_of_detectable_civilizations); + + // Expected number of civilizations in the Milky way; + // see footnote 3 (p. 5) + double n = rate_of_star_formation * fraction_of_stars_with_planets * number_of_habitable_planets_per_star_system * fraction_of_habitable_planets_in_which_any_life_appears * fraction_of_planets_with_life_in_which_intelligent_life_appears * fraction_of_intelligent_planets_which_are_detectable_as_such * longevity_of_detectable_civilizations; + + return n; + } + + double sample_fermi_paradox_naive(uint64_t * seed) + { + double n = sample_fermi_naive(seed); + return ((n > 1) ? 1 : 0); + } + + double n = 1000000; + double naive_fermi_proportion = 0; + for (int i = 0; i < n; i++) { + double result = sample_fermi_paradox_naive(seed); + // printf("result: %lf\n", result); + naive_fermi_proportion += result; + } + printf("Naïve %% that we are not alone: %lf\n", naive_fermi_proportion / n); + + // Thinking in log space + double sample_fermi_logspace(uint64_t * seed) + { + double log_rate_of_star_formation = sample_uniform(log(1), log(100), seed); + double log_fraction_of_stars_with_planets = sample_uniform(log(0.1), log(1), seed); + double log_number_of_habitable_planets_per_star_system = sample_uniform(log(0.1), log(1), seed); + double log_fraction_of_planets_with_life_in_which_intelligent_life_appears = sample_uniform(log(0.001), log(1), seed); + double log_fraction_of_intelligent_planets_which_are_detectable_as_such = sample_uniform(log(0.01), log(1), seed); + double log_longevity_of_detectable_civilizations = sample_uniform(log(100), log(10000000000), seed); + + // printf(" log_rate_of_star_formation = %lf\n", log_rate_of_star_formation); + // printf(" log_fraction_of_stars_with_planets = %lf\n", log_fraction_of_stars_with_planets); + // printf(" log_number_of_habitable_planets_per_star_system = %lf\n", log_number_of_habitable_planets_per_star_system); + // printf(" log_fraction_of_planets_with_life_in_which_intelligent_life_appears = %lf\n", log_fraction_of_planets_with_life_in_which_intelligent_life_appears); + // printf(" log_fraction_of_intelligent_planets_which_are_detectable_as_such = %lf\n", log_fraction_of_intelligent_planets_which_are_detectable_as_such); + // printf(" log_longevity_of_detectable_civilizations = %lf\n", log_longevity_of_detectable_civilizations); + + double log_n1 = log_rate_of_star_formation + log_fraction_of_stars_with_planets + log_number_of_habitable_planets_per_star_system + log_fraction_of_planets_with_life_in_which_intelligent_life_appears + log_fraction_of_intelligent_planets_which_are_detectable_as_such + log_longevity_of_detectable_civilizations; + // printf("first part of calculation: %lf\n", log_n1); + + /* Consider fraction_of_habitable_planets_in_which_any_life_appears separately. + Imprecisely, we could do: + + double rate_of_life_formation_in_habitable_planets = sample_lognormal(1, 50, seed); + double fraction_of_habitable_planets_in_which_any_life_appears = 1- exp(-rate_of_life_formation_in_habitable_planets); + double log_fraction_of_habitable_planets_in_which_any_life_appears = log(1-fraction_of_habitable_planets_in_which_any_life_appears); + double n = exp(log_n1) * fraction_of_habitable_planets_in_which_any_life_appears; + // or: + double n2 = exp(log_n1 + log(fraction_of_habitable_planets_in_which_any_life_appears)) + + However, we lose all precision here. + + Now, say + a = underlying normal + b = rate_of_life_formation_in_habitable_planets = exp(underlying normal) + c = 1 - exp(-b) = fraction_of_habitable_planets_in_which_any_life_appears + d = log(c) + + Now, is there some way we can d more efficiently/precisely? + Turns out there is! + + Looking at the Taylor expansion for c = 1 - exp(-b), it's b - b^2/2 + b^3/6 - x^b/24, etc. + // https://www.wolframalpha.com/input?i=1-exp%28-x%29 + When b ~ 0 (as is often the case), this is close to b. + + But now, if b ~ 0 + c ~ b + and d = log(c) ~ log(b) = log(exp(a)) = a + */ + double log_rate_of_life_formation_in_habitable_planets = sample_normal(1, 50, seed); + // printf("log_rate_of_life_formation_in_habitable_planets: %lf\n", log_rate_of_life_formation_in_habitable_planets); + + double log_fraction_of_habitable_planets_in_which_any_life_appears; + if (log_rate_of_life_formation_in_habitable_planets < -32) { + log_fraction_of_habitable_planets_in_which_any_life_appears = log_rate_of_life_formation_in_habitable_planets; + } else { + double rate_of_life_formation_in_habitable_planets = exp(log_rate_of_life_formation_in_habitable_planets); + double fraction_of_habitable_planets_in_which_any_life_appears = -expm1(-rate_of_life_formation_in_habitable_planets); + log_fraction_of_habitable_planets_in_which_any_life_appears = log(fraction_of_habitable_planets_in_which_any_life_appears); + } + // printf(" log_fraction_of_habitable_planets_in_which_any_life_appears: %lf\n", log_fraction_of_habitable_planets_in_which_any_life_appears); + + double log_n = log_n1 + log_fraction_of_habitable_planets_in_which_any_life_appears; + + return log_n; + } + + double sample_fermi_paradox_logspace(uint64_t * seed) + { + double n = sample_fermi_logspace(seed); + return ((n > 0) ? 1 : 0); + } + + double logspace_fermi_proportion = 0; + for (int i = 0; i < n; i++) { + double result = sample_fermi_paradox_logspace(seed); + // printf("result: %lf\n", result); + logspace_fermi_proportion += result; + } + printf("Using more accurate logspace computations, %% that we are not alone: %lf\n", logspace_fermi_proportion / n); + + free(seed); +} diff --git a/examples/core/06_dissolving_fermi_paradox/scratchpad.c b/examples/core/06_dissolving_fermi_paradox/scratchpad.c @@ -1,157 +0,0 @@ -#include "../squiggle.h" -// #include "../squiggle_more.h" -#include <math.h> -#include <stdint.h> -#include <stdio.h> -#include <stdlib.h> - -double sample_loguniform(double a, double b, uint64_t* seed){ - return exp(sample_uniform(log(a), log(b), seed)); -} - -int main() -{ - // Replicate <https://arxiv.org/pdf/1806.02404.pdf>, and in particular the red line in page 11. - // Could also be interesting to just produce and save many samples. - - // set randomness seed - uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = UINT64_MAX/64; // xorshift can't start with a seed of 0 - - double sample_fermi_naive(uint64_t* seed){ - double rate_of_star_formation = sample_loguniform(1,100, seed); - double fraction_of_stars_with_planets = sample_loguniform(0.1, 1, seed); - double number_of_habitable_planets_per_star_system = sample_loguniform(0.1, 1, seed); - double rate_of_life_formation_in_habitable_planets = sample_lognormal(1, 50, seed); - double fraction_of_habitable_planets_in_which_any_life_appears = -expm1(-rate_of_life_formation_in_habitable_planets); - // double fraction_of_habitable_planets_in_which_any_life_appears = 1-exp(-rate_of_life_formation_in_habitable_planets); - // but with more precision - double fraction_of_planets_with_life_in_which_intelligent_life_appears = sample_loguniform(0.001, 1, seed); - double fraction_of_intelligent_planets_which_are_detectable_as_such = sample_loguniform(0.01, 1, seed); - double longevity_of_detectable_civilizations = sample_loguniform(100, 10000000000, seed); - - // printf(" rate_of_star_formation = %lf\n", rate_of_star_formation); - // printf(" fraction_of_stars_with_planets = %lf\n", fraction_of_stars_with_planets); - // printf(" number_of_habitable_planets_per_star_system = %lf\n", number_of_habitable_planets_per_star_system); - // printf(" rate_of_life_formation_in_habitable_planets = %.16lf\n", rate_of_life_formation_in_habitable_planets); - // printf(" fraction_of_habitable_planets_in_which_any_life_appears = %lf\n", fraction_of_habitable_planets_in_which_any_life_appears); - // printf(" fraction_of_planets_with_life_in_which_intelligent_life_appears = %lf\n", fraction_of_planets_with_life_in_which_intelligent_life_appears); - // printf(" fraction_of_intelligent_planets_which_are_detectable_as_such = %lf\n", fraction_of_intelligent_planets_which_are_detectable_as_such); - // printf(" longevity_of_detectable_civilizations = %lf\n", longevity_of_detectable_civilizations); - - // Expected number of civilizations in the Milky way; - // see footnote 3 (p. 5) - double n = rate_of_star_formation * - fraction_of_stars_with_planets * - number_of_habitable_planets_per_star_system * - fraction_of_habitable_planets_in_which_any_life_appears * - fraction_of_planets_with_life_in_which_intelligent_life_appears * - fraction_of_intelligent_planets_which_are_detectable_as_such * - longevity_of_detectable_civilizations; - - return n; - } - - double sample_fermi_paradox_naive(uint64_t* seed){ - double n = sample_fermi_naive(seed); - return ((n > 1) ? 1 : 0); - } - - double n = 1000000; - double naive_fermi_proportion = 0; - for(int i=0; i<n; i++){ - double result = sample_fermi_paradox_naive(seed); - // printf("result: %lf\n", result); - naive_fermi_proportion+=result; - } - printf("Naïve %% that we are not alone: %lf\n", naive_fermi_proportion/n); - - - // Thinking in log space - double sample_fermi_logspace(uint64_t* seed){ - double log_rate_of_star_formation = sample_uniform(log(1), log(100), seed); - double log_fraction_of_stars_with_planets = sample_uniform(log(0.1), log(1), seed); - double log_number_of_habitable_planets_per_star_system = sample_uniform(log(0.1), log(1), seed); - double log_fraction_of_planets_with_life_in_which_intelligent_life_appears = sample_uniform(log(0.001), log(1), seed); - double log_fraction_of_intelligent_planets_which_are_detectable_as_such = sample_uniform(log(0.01), log(1), seed); - double log_longevity_of_detectable_civilizations = sample_uniform(log(100), log(10000000000), seed); - - // printf(" log_rate_of_star_formation = %lf\n", log_rate_of_star_formation); - // printf(" log_fraction_of_stars_with_planets = %lf\n", log_fraction_of_stars_with_planets); - // printf(" log_number_of_habitable_planets_per_star_system = %lf\n", log_number_of_habitable_planets_per_star_system); - // printf(" log_fraction_of_planets_with_life_in_which_intelligent_life_appears = %lf\n", log_fraction_of_planets_with_life_in_which_intelligent_life_appears); - // printf(" log_fraction_of_intelligent_planets_which_are_detectable_as_such = %lf\n", log_fraction_of_intelligent_planets_which_are_detectable_as_such); - // printf(" log_longevity_of_detectable_civilizations = %lf\n", log_longevity_of_detectable_civilizations); - - double log_n1 = - log_rate_of_star_formation + - log_fraction_of_stars_with_planets + - log_number_of_habitable_planets_per_star_system + - log_fraction_of_planets_with_life_in_which_intelligent_life_appears + - log_fraction_of_intelligent_planets_which_are_detectable_as_such + - log_longevity_of_detectable_civilizations; - // printf("first part of calculation: %lf\n", log_n1); - - /* Consider fraction_of_habitable_planets_in_which_any_life_appears separately. - Imprecisely, we could do: - - double rate_of_life_formation_in_habitable_planets = sample_lognormal(1, 50, seed); - double fraction_of_habitable_planets_in_which_any_life_appears = 1- exp(-rate_of_life_formation_in_habitable_planets); - double log_fraction_of_habitable_planets_in_which_any_life_appears = log(1-fraction_of_habitable_planets_in_which_any_life_appears); - double n = exp(log_n1) * fraction_of_habitable_planets_in_which_any_life_appears; - // or: - double n2 = exp(log_n1 + log(fraction_of_habitable_planets_in_which_any_life_appears)) - - However, we lose all precision here. - - Now, say - a = underlying normal - b = rate_of_life_formation_in_habitable_planets = exp(underlying normal) - c = 1 - exp(-b) = fraction_of_habitable_planets_in_which_any_life_appears - d = log(c) - - Now, is there some way we can d more efficiently/precisely? - Turns out there is! - - Looking at the Taylor expansion for c = 1 - exp(-b), it's b - b^2/2 + b^3/6 - x^b/24, etc. - // https://www.wolframalpha.com/input?i=1-exp%28-x%29 - When b ~ 0 (as is often the case), this is close to b. - - But now, if b ~ 0 - c ~ b - and d = log(c) ~ log(b) = log(exp(a)) = a - */ - double log_rate_of_life_formation_in_habitable_planets = sample_normal(1, 50, seed); - // printf("log_rate_of_life_formation_in_habitable_planets: %lf\n", log_rate_of_life_formation_in_habitable_planets); - - double log_fraction_of_habitable_planets_in_which_any_life_appears; - if(log_rate_of_life_formation_in_habitable_planets < -32){ - log_fraction_of_habitable_planets_in_which_any_life_appears = log_rate_of_life_formation_in_habitable_planets; - } else{ - double rate_of_life_formation_in_habitable_planets = exp(log_rate_of_life_formation_in_habitable_planets); - double fraction_of_habitable_planets_in_which_any_life_appears = -expm1(-rate_of_life_formation_in_habitable_planets); - log_fraction_of_habitable_planets_in_which_any_life_appears = log(fraction_of_habitable_planets_in_which_any_life_appears); - } - // printf(" log_fraction_of_habitable_planets_in_which_any_life_appears: %lf\n", log_fraction_of_habitable_planets_in_which_any_life_appears); - - double log_n = log_n1 + log_fraction_of_habitable_planets_in_which_any_life_appears; - - return log_n; - } - - double sample_fermi_paradox_logspace(uint64_t* seed){ - double n = sample_fermi_logspace(seed); - return ((n > 0) ? 1 : 0); - } - - double logspace_fermi_proportion = 0; - for(int i=0; i<n; i++){ - double result = sample_fermi_paradox_logspace(seed); - // printf("result: %lf\n", result); - logspace_fermi_proportion+=result; - } - printf("Using more accurate logspace computations, %% that we are not alone: %lf\n", logspace_fermi_proportion/n); - double result2; - - free(seed); -} diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad Binary files differ. diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c @@ -1,5 +1,5 @@ #include "../squiggle.h" -#include "../squiggle_more.h" +// #include "../squiggle_more.h" #include <math.h> #include <stdint.h> #include <stdio.h> @@ -7,16 +7,18 @@ double sample_loguniform(double a, double b, uint64_t* seed){ return exp(sample_uniform(log(a), log(b), seed)); - } int main() { + // Replicate <https://arxiv.org/pdf/1806.02404.pdf>, and in particular the red line in page 11. + // Could also be interesting to just produce and save many samples. + // set randomness seed uint64_t* seed = malloc(sizeof(uint64_t)); *seed = UINT64_MAX/64; // xorshift can't start with a seed of 0 - double fermi_naive(uint64_t* seed){ + double sample_fermi_naive(uint64_t* seed){ double rate_of_star_formation = sample_loguniform(1,100, seed); double fraction_of_stars_with_planets = sample_loguniform(0.1, 1, seed); double number_of_habitable_planets_per_star_system = sample_loguniform(0.1, 1, seed); @@ -50,21 +52,23 @@ int main() return n; } - double fermi_paradox_naive(uint64_t* seed){ - double n = fermi_naive(seed); - return (n > 1 ? 1 : 0); + double sample_fermi_paradox_naive(uint64_t* seed){ + double n = sample_fermi_naive(seed); + return ((n > 1) ? 1 : 0); } - double result; - for(int i=0; i<1000; i++){ - result = fermi_naive(seed); - printf("result from fermi_naive: %lf\n", result); - printf("\n\n"); + double n = 1000000; + double naive_fermi_proportion = 0; + for(int i=0; i<n; i++){ + double result = sample_fermi_paradox_naive(seed); + // printf("result: %lf\n", result); + naive_fermi_proportion+=result; } - printf("result from naïve implementation: %lf\n", result); + printf("Naïve %% that we are not alone: %lf\n", naive_fermi_proportion/n); + // Thinking in log space - double fermi_logspace(uint64_t* seed){ + double sample_fermi_logspace(uint64_t* seed){ double log_rate_of_star_formation = sample_uniform(log(1), log(100), seed); double log_fraction_of_stars_with_planets = sample_uniform(log(0.1), log(1), seed); double log_number_of_habitable_planets_per_star_system = sample_uniform(log(0.1), log(1), seed); @@ -86,7 +90,7 @@ int main() log_fraction_of_planets_with_life_in_which_intelligent_life_appears + log_fraction_of_intelligent_planets_which_are_detectable_as_such + log_longevity_of_detectable_civilizations; - printf("first part of calculation: %lf\n", log_n1); + // printf("first part of calculation: %lf\n", log_n1); /* Consider fraction_of_habitable_planets_in_which_any_life_appears separately. Imprecisely, we could do: @@ -110,14 +114,15 @@ int main() Turns out there is! Looking at the Taylor expansion for c = 1 - exp(-b), it's b - b^2/2 + b^3/6 - x^b/24, etc. - When b ~ 0 (as it is), this is close to b. + // https://www.wolframalpha.com/input?i=1-exp%28-x%29 + When b ~ 0 (as is often the case), this is close to b. But now, if b ~ 0 c ~ b and d = log(c) ~ log(b) = log(exp(a)) = a */ double log_rate_of_life_formation_in_habitable_planets = sample_normal(1, 50, seed); - printf("log_rate_of_life_formation_in_habitable_planets: %lf\n", log_rate_of_life_formation_in_habitable_planets); + // printf("log_rate_of_life_formation_in_habitable_planets: %lf\n", log_rate_of_life_formation_in_habitable_planets); double log_fraction_of_habitable_planets_in_which_any_life_appears; if(log_rate_of_life_formation_in_habitable_planets < -32){ @@ -127,22 +132,26 @@ int main() double fraction_of_habitable_planets_in_which_any_life_appears = -expm1(-rate_of_life_formation_in_habitable_planets); log_fraction_of_habitable_planets_in_which_any_life_appears = log(fraction_of_habitable_planets_in_which_any_life_appears); } - printf(" log_fraction_of_habitable_planets_in_which_any_life_appears: %lf\n", log_fraction_of_habitable_planets_in_which_any_life_appears); + // printf(" log_fraction_of_habitable_planets_in_which_any_life_appears: %lf\n", log_fraction_of_habitable_planets_in_which_any_life_appears); double log_n = log_n1 + log_fraction_of_habitable_planets_in_which_any_life_appears; return log_n; } - double result2; + double sample_fermi_paradox_logspace(uint64_t* seed){ + double n = sample_fermi_logspace(seed); + return ((n > 0) ? 1 : 0); + } - /* - for(int i=0; i<1000; i++){ - result2 = fermi_logspace(seed); - printf("result from logspace implementation: %lf.2\n", result2); - printf("\n\n"); + double logspace_fermi_proportion = 0; + for(int i=0; i<n; i++){ + double result = sample_fermi_paradox_logspace(seed); + // printf("result: %lf\n", result); + logspace_fermi_proportion+=result; } - */ + printf("Using more accurate logspace computations, %% that we are not alone: %lf\n", logspace_fermi_proportion/n); + double result2; free(seed); }