diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad index 0faebb6..5817d41 100755 Binary files a/scratchpad/scratchpad and b/scratchpad/scratchpad differ diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c index 64932b6..bd763c1 100644 --- a/scratchpad/scratchpad.c +++ b/scratchpad/scratchpad.c @@ -5,21 +5,144 @@ #include #include +double sample_loguniform(double a, double b, uint64_t* seed){ + return exp(sample_uniform(log(a), log(b), seed)); + +} + int main() { // set randomness seed uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with a seed of 0 + *seed = UINT64_MAX/64; // xorshift can't start with a seed of 0 - int n = 1000000; - double* xs = malloc(sizeof(double) * n); - for (int i = 0; i < n; i++) { - xs[i] = sample_to(10, 100, seed); + double 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; } - ci ci_90 = array_get_90_ci(xs, n); - printf("Recovering confidence interval of sample_to(10, 100):\n low: %f, high: %f\n", ci_90.low, ci_90.high); + double fermi_paradox_naive(uint64_t* seed){ + double n = 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"); + } + printf("result from naïve implementation: %lf\n", result); + + // Thinking in log space + double 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. + When b ~ 0 (as it is), 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 result2; + + /* + for(int i=0; i<1000; i++){ + result2 = fermi_logspace(seed); + printf("result from logspace implementation: %lf.2\n", result2); + printf("\n\n"); + } + */ - printf("Size of uint64_t: %ld", sizeof(uint64_t*)); free(seed); }