diff --git a/examples/core/06_dissolving_fermi_paradox/example b/examples/core/06_dissolving_fermi_paradox/example index f6ac909..45bb810 100755 Binary files a/examples/core/06_dissolving_fermi_paradox/example and b/examples/core/06_dissolving_fermi_paradox/example differ diff --git a/examples/core/06_dissolving_fermi_paradox/example.c b/examples/core/06_dissolving_fermi_paradox/example.c index 0e9a5e6..8252569 100644 --- a/examples/core/06_dissolving_fermi_paradox/example.c +++ b/examples/core/06_dissolving_fermi_paradox/example.c @@ -4,6 +4,8 @@ #include #include +#define VERBOSE 0 + double sample_loguniform(double a, double b, uint64_t* seed) { return exp(sample_uniform(log(a), log(b), seed)); @@ -18,6 +20,7 @@ int main() uint64_t* seed = malloc(sizeof(uint64_t)); *seed = UINT64_MAX / 64; // xorshift can't start with a seed of 0 + // Do this naïvely, without worrying that much about numerical precision double sample_fermi_naive(uint64_t * seed) { double rate_of_star_formation = sample_loguniform(1, 100, seed); @@ -31,14 +34,14 @@ int main() 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); + if(VERBOSE) printf(" rate_of_star_formation = %lf\n", rate_of_star_formation); + if(VERBOSE) printf(" fraction_of_stars_with_planets = %lf\n", fraction_of_stars_with_planets); + if(VERBOSE) printf(" number_of_habitable_planets_per_star_system = %lf\n", number_of_habitable_planets_per_star_system); + if(VERBOSE) printf(" rate_of_life_formation_in_habitable_planets = %.16lf\n", rate_of_life_formation_in_habitable_planets); + if(VERBOSE) printf(" fraction_of_habitable_planets_in_which_any_life_appears = %lf\n", fraction_of_habitable_planets_in_which_any_life_appears); + if(VERBOSE) printf(" fraction_of_planets_with_life_in_which_intelligent_life_appears = %lf\n", fraction_of_planets_with_life_in_which_intelligent_life_appears); + if(VERBOSE) printf(" fraction_of_intelligent_planets_which_are_detectable_as_such = %lf\n", fraction_of_intelligent_planets_which_are_detectable_as_such); + if(VERBOSE) printf(" longevity_of_detectable_civilizations = %lf\n", longevity_of_detectable_civilizations); // Expected number of civilizations in the Milky way; // see footnote 3 (p. 5) @@ -47,7 +50,7 @@ int main() return n; } - double sample_fermi_paradox_naive(uint64_t * seed) + double sample_are_we_alone_naive(uint64_t * seed) { double n = sample_fermi_naive(seed); return ((n > 1) ? 1 : 0); @@ -56,13 +59,14 @@ int main() 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); + double result = sample_are_we_alone_naive(seed); + if(VERBOSE) 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 + // Taking care of numerical precision double sample_fermi_logspace(uint64_t * seed) { double log_rate_of_star_formation = sample_uniform(log(1), log(100), seed); @@ -72,64 +76,63 @@ int main() 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); + if(VERBOSE) printf(" log_rate_of_star_formation = %lf\n", log_rate_of_star_formation); + if(VERBOSE) printf(" log_fraction_of_stars_with_planets = %lf\n", log_fraction_of_stars_with_planets); + if(VERBOSE) printf(" log_number_of_habitable_planets_per_star_system = %lf\n", log_number_of_habitable_planets_per_star_system); + if(VERBOSE) 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); + if(VERBOSE) printf(" log_fraction_of_intelligent_planets_which_are_detectable_as_such = %lf\n", log_fraction_of_intelligent_planets_which_are_detectable_as_such); + if(VERBOSE) 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); + if(VERBOSE) printf("first part of calculation: %lf\n", log_n1); - /* Consider fraction_of_habitable_planets_in_which_any_life_appears separately. - Imprecisely, we could do: + /* + Consider: + a = underlying normal + b = rate_of_life_formation_in_habitable_planets = exp(underlying normal) = exp(a) + c = 1 - exp(-b) = fraction_of_habitable_planets_in_which_any_life_appears + d = log(c) - 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)) + Now, is there some way we can get d more efficiently/precisely? - However, we lose all precision here. + 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 is often the case), this is close to b. - 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) + But now, if b ~ 0 + c ~ b + and d = log(c) ~ log(b) = log(exp(a)) = a - 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 + Now, we could play around with estimating errors, + and indeed if we want b^2/2 = exp(a)^2/2 < 10^(-n), i.e., to have n decimal digits of precision, + we could compute this as e.g., a < (nlog(10) + log(2))/2 + so for example if we want ten digits of precision, that's a < -11 + + But more empirically, the two numbers do become really close around 11 or so, and at 38 that calculation results in a -inf (so probably an overflow.) + So we should be using that formula for somewhere between -38 << a < -11 + I chose -16 for the sake of it after playing with: + */ 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); + if(VERBOSE) 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) { + if (log_rate_of_life_formation_in_habitable_planets < -16) { 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); + if(VERBOSE) 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 sample_are_we_alone_logspace(uint64_t * seed) { double n = sample_fermi_logspace(seed); return ((n > 0) ? 1 : 0); @@ -137,11 +140,21 @@ int main() double logspace_fermi_proportion = 0; for (int i = 0; i < n; i++) { - double result = sample_fermi_paradox_logspace(seed); - // printf("result: %lf\n", result); + double result = sample_are_we_alone_logspace(seed); + if(VERBOSE) 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); + + /* + double invert(double x){ + return log(1-exp(-exp(-x))); + } + for(int i=0; i<64; i++){ + double j = i; + printf("for %lf, log(1-exp(-exp(-x))) is calculated as... %lf\n", j, invert(j)); + } + */ } diff --git a/examples/core/06_dissolving_fermi_paradox/makefile b/examples/core/06_dissolving_fermi_paradox/makefile deleted file mode 100644 index d515383..0000000 --- a/examples/core/06_dissolving_fermi_paradox/makefile +++ /dev/null @@ -1,56 +0,0 @@ -# Interface: -# make -# make build -# make format -# make run - -# Compiler -CC=gcc -# CC=tcc # <= faster compilation - -# Main file -SRC=scratchpad.c ../squiggle.c ../squiggle_more.c -OUTPUT=scratchpad - -## Dependencies -MATH=-lm - -## Flags -DEBUG= #'-g' -STANDARD=-std=c99 -WARNINGS=-Wall -OPTIMIZED=-O3 #-Ofast -# OPENMP=-fopenmp - -## Formatter -STYLE_BLUEPRINT=webkit -FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT) - -## make build -build: $(SRC) - $(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(MATH) -o $(OUTPUT) - -format: $(SRC) - $(FORMATTER) $(SRC) - -run: $(SRC) $(OUTPUT) - ./$(OUTPUT) - -verify: $(SRC) $(OUTPUT) - ./$(OUTPUT) | grep "NOT passed" -A 2 --group-separator='' || true - -time-linux: - @echo "Requires /bin/time, found on GNU/Linux systems" && echo - - @echo "Running 100x and taking avg time $(OUTPUT)" - @t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 1 thread: |" | sed 's|$$|ms|' && echo - -## Profiling - -profile-linux: - echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar" - echo "Must be run as sudo" - $(CC) $(SRC) $(MATH) -o $(OUTPUT) - sudo perf record ./$(OUTPUT) - sudo perf report - rm perf.data