commit b0f711e9a6e11278330fe9670506649208e6a8fa parent 3f9027a530edb48f361e4de3ce781794acdd5d19 Author: NunoSempere <nuno.sempere@protonmail.com> Date: Sun, 19 Nov 2023 14:47:19 +0000 reorg: move all examples using squiggle_more to one makefile Diffstat:
45 files changed, 726 insertions(+), 666 deletions(-)
diff --git a/examples/more/00_example_template/example b/examples/more/00_example_template/example Binary files differ. diff --git a/examples/more/00_example_template/example.c b/examples/more/00_example_template/example.c @@ -0,0 +1,16 @@ +#include "../../../squiggle.h" +#include "../../../squiggle_more.h" +#include <stdint.h> +#include <stdio.h> +#include <stdlib.h> + +int main() +{ + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 + + // ... + + free(seed); +} diff --git a/examples/more/01_sample_from_cdf/example b/examples/more/01_sample_from_cdf/example Binary files differ. diff --git a/examples/more/01_sample_from_cdf/example.c b/examples/more/01_sample_from_cdf/example.c @@ -0,0 +1,103 @@ +#include "../../../squiggle.h" +#include "../../../squiggle_more.h" +#include <math.h> +#include <stdint.h> +#include <stdio.h> +#include <stdlib.h> +#include <time.h> + +#define NUM_SAMPLES 1000000 + +// Example cdf +double cdf_uniform_0_1(double x) +{ + if (x < 0) { + return 0; + } else if (x > 1) { + return 1; + } else { + return x; + } +} + +double cdf_squared_0_1(double x) +{ + if (x < 0) { + return 0; + } else if (x > 1) { + return 1; + } else { + return x * x; + } +} + +double cdf_normal_0_1(double x) +{ + double mean = 0; + double std = 1; + return 0.5 * (1 + erf((x - mean) / (std * sqrt(2)))); // erf from math.h +} + +// Some testers +void test_inverse_cdf_double(char* cdf_name, double cdf_double(double)) +{ + struct box result = inverse_cdf_double(cdf_double, 0.5); + if (result.empty) { + printf("Inverse for %s not calculated\n", cdf_name); + exit(1); + } else { + printf("Inverse of %s at %f is: %f\n", cdf_name, 0.5, result.content); + } +} + +void test_and_time_sampler_double(char* cdf_name, double cdf_double(double), uint64_t* seed) +{ + printf("\nGetting some samples from %s:\n", cdf_name); + clock_t begin = clock(); + for (int i = 0; i < NUM_SAMPLES; i++) { + struct box sample = sampler_cdf_double(cdf_double, seed); + if (sample.empty) { + printf("Error in sampler function for %s", cdf_name); + } else { + // printf("%f\n", sample.content); + } + } + clock_t end = clock(); + double time_spent = (double)(end - begin) / CLOCKS_PER_SEC; + printf("Time spent: %f\n", time_spent); +} + +int main() +{ + // Test inverse cdf double + test_inverse_cdf_double("cdf_uniform_0_1", cdf_uniform_0_1); + test_inverse_cdf_double("cdf_squared_0_1", cdf_squared_0_1); + test_inverse_cdf_double("cdf_normal_0_1", cdf_normal_0_1); + + // Testing samplers + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 + + // Test double sampler + test_and_time_sampler_double("cdf_uniform_0_1", cdf_uniform_0_1, seed); + test_and_time_sampler_double("cdf_squared_0_1", cdf_squared_0_1, seed); + test_and_time_sampler_double("cdf_normal_0_1", cdf_normal_0_1, seed); + + // Get some normal samples using a previous approach + printf("\nGetting some samples from sample_unit_normal\n"); + + clock_t begin_2 = clock(); + double* normal_samples = malloc(NUM_SAMPLES * sizeof(double)); + for (int i = 0; i < NUM_SAMPLES; i++) { + normal_samples[i] = sample_unit_normal(seed); + // printf("%f\n", normal_sample); + } + + clock_t end_2 = clock(); + double time_spent_2 = (double)(end_2 - begin_2) / CLOCKS_PER_SEC; + printf("Time spent: %f\n", time_spent_2); + + free(seed); + return 0; +} diff --git a/examples/more/02_sample_from_cdf_beta/example b/examples/more/02_sample_from_cdf_beta/example Binary files differ. diff --git a/examples/more/02_sample_from_cdf_beta/example.c b/examples/more/02_sample_from_cdf_beta/example.c @@ -0,0 +1,169 @@ +#include "../../../squiggle.h" +#include "../../../squiggle_more.h" +#include <math.h> +#include <stdint.h> +#include <stdio.h> +#include <stdlib.h> +#include <time.h> + +#define NUM_SAMPLES 10000 +#define STOP_BETA 1.0e-8 +#define TINY_BETA 1.0e-30 + +// Incomplete beta function +struct box incbeta(double a, double b, double x) +{ + // Descended from <https://github.com/codeplea/incbeta/blob/master/incbeta.c>, + // <https://codeplea.com/incomplete-beta-function-c> + // but modified to return a box struct and doubles instead of doubles. + // [ ] to do: add attribution in README + // Original code under this license: + /* + * zlib License + * + * Regularized Incomplete Beta Function + * + * Copyright (c) 2016, 2017 Lewis Van Winkle + * http://CodePlea.com + * + * This software is provided 'as-is', without any express or implied + * warranty. In no event will the authors be held liable for any damages + * arising from the use of this software. + * + * Permission is granted to anyone to use this software for any purpose, + * including commercial applications, and to alter it and redistribute it + * freely, subject to the following restrictions: + * + * 1. The origin of this software must not be misrepresented; you must not + * claim that you wrote the original software. If you use this software + * in a product, an acknowledgement in the product documentation would be + * appreciated but is not required. + * 2. Altered source versions must be plainly marked as such, and must not be + * misrepresented as being the original software. + * 3. This notice may not be removed or altered from any source distribution. + */ + if (x < 0.0 || x > 1.0) { + return PROCESS_ERROR("x out of bounds [0, 1], in function incbeta"); + } + + /*The continued fraction converges nicely for x < (a+1)/(a+b+2)*/ + if (x > (a + 1.0) / (a + b + 2.0)) { + struct box symmetric_incbeta = incbeta(b, a, 1.0 - x); + if (symmetric_incbeta.empty) { + return symmetric_incbeta; // propagate error + } else { + struct box result = { + .empty = 0, + .content = 1 - symmetric_incbeta.content + }; + return result; + } + } + + /*Find the first part before the continued fraction.*/ + const double lbeta_ab = lgamma(a) + lgamma(b) - lgamma(a + b); + const double front = exp(log(x) * a + log(1.0 - x) * b - lbeta_ab) / a; + + /*Use Lentz's algorithm to evaluate the continued fraction.*/ + double f = 1.0, c = 1.0, d = 0.0; + + int i, m; + for (i = 0; i <= 200; ++i) { + m = i / 2; + + double numerator; + if (i == 0) { + numerator = 1.0; /*First numerator is 1.0.*/ + } else if (i % 2 == 0) { + numerator = (m * (b - m) * x) / ((a + 2.0 * m - 1.0) * (a + 2.0 * m)); /*Even term.*/ + } else { + numerator = -((a + m) * (a + b + m) * x) / ((a + 2.0 * m) * (a + 2.0 * m + 1)); /*Odd term.*/ + } + + /*Do an iteration of Lentz's algorithm.*/ + d = 1.0 + numerator * d; + if (fabs(d) < TINY_BETA) + d = TINY_BETA; + d = 1.0 / d; + + c = 1.0 + numerator / c; + if (fabs(c) < TINY_BETA) + c = TINY_BETA; + + const double cd = c * d; + f *= cd; + + /*Check for stop.*/ + if (fabs(1.0 - cd) < STOP_BETA) { + struct box result = { + .empty = 0, + .content = front * (f - 1.0) + }; + return result; + } + } + + return PROCESS_ERROR("More loops needed, did not converge, in function incbeta"); +} + +struct box cdf_beta(double x) +{ + if (x < 0) { + struct box result = { .empty = 0, .content = 0 }; + return result; + } else if (x > 1) { + struct box result = { .empty = 0, .content = 1 }; + return result; + } else { + double successes = 1, failures = (2023 - 1945); + return incbeta(successes, failures, x); + } +} + +// Some testers +void test_inverse_cdf_box(char* cdf_name, struct box cdf_box(double)) +{ + struct box result = inverse_cdf_box(cdf_box, 0.5); + if (result.empty) { + printf("Inverse for %s not calculated\n", cdf_name); + exit(1); + } else { + printf("Inverse of %s at %f is: %f\n", cdf_name, 0.5, result.content); + } +} + +void test_and_time_sampler_box(char* cdf_name, struct box cdf_box(double), uint64_t* seed) +{ + printf("\nGetting some samples from %s:\n", cdf_name); + clock_t begin = clock(); + for (int i = 0; i < NUM_SAMPLES; i++) { + struct box sample = sampler_cdf_box(cdf_box, seed); + if (sample.empty) { + printf("Error in sampler function for %s", cdf_name); + } else { + // printf("%f\n", sample.content); + } + } + clock_t end = clock(); + double time_spent = (double)(end - begin) / CLOCKS_PER_SEC; + printf("Time spent: %f\n", time_spent); +} + +int main() +{ + // Test inverse cdf box + test_inverse_cdf_box("cdf_beta", cdf_beta); + + // Test box sampler + 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!! + // Partly this is because I am using a more general algorithm, + // which applies to any cdf + // But I am also using absurdly precise convergence conditions. + // This could be optimized. + + free(seed); + return 0; +} diff --git a/examples/more/07_ci_beta/example b/examples/more/03_ci_beta/example Binary files differ. diff --git a/examples/more/03_ci_beta/example.c b/examples/more/03_ci_beta/example.c @@ -0,0 +1,23 @@ +#include "../../../squiggle.h" +#include "../../../squiggle_more.h" +#include <stdint.h> +#include <stdio.h> +#include <stdlib.h> + +// Estimate functions +double beta_1_2_sampler(uint64_t* seed) +{ + return sample_beta(1, 2.0, seed); +} + +int main() +{ + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 + + ci beta_1_2_ci_90 = get_90_confidence_interval(beta_1_2_sampler, seed); + printf("90%% confidence interval of beta(1,2) is [%f, %f]\n", beta_1_2_ci_90.low, beta_1_2_ci_90.high); + + free(seed); +} diff --git a/examples/more/04_nuclear_war/example b/examples/more/04_nuclear_war/example Binary files differ. diff --git a/examples/more/04_nuclear_war/example.c b/examples/more/04_nuclear_war/example.c @@ -0,0 +1,69 @@ +#include "../../../squiggle.h" +#include "../../../squiggle_more.h" +#include <math.h> +#include <stdint.h> +#include <stdio.h> +#include <stdlib.h> + +double probability_of_dying_nuno(uint64_t* seed) +{ + double first_year_russian_nuclear_weapons = 1953; + double current_year = 2022; + double laplace_probability_nuclear_exchange_year = sample_beta(1, current_year - first_year_russian_nuclear_weapons + 1, seed); + double laplace_probability_nuclear_exchange_month = 1 - pow(1 - laplace_probability_nuclear_exchange_year, (1.0 / 12.0)); + + double london_hit_conditional_on_russia_nuclear_weapon_usage = sample_beta(7.67, 69.65, seed); + // I.e., a beta distribution with a range of 0.05 to 0.16 into: https://nunosempere.com/blog/2023/03/15/fit-beta/ + // 0.05 were my estimate and Samotsvety's estimate in March 2022, respectively: + // https://forum.effectivealtruism.org/posts/KRFXjCqqfGQAYirm5/samotsvety-nuclear-risk-forecasts-march-2022#Nu_o_Sempere + double informed_actor_not_able_to_escape = sample_beta(3.26212166586967, 3.26228162008564, seed); + // 0.2 to 0.8, i.e., 20% to 80%, again using the previous tool + double proportion_which_die_if_bomb_drops_in_london = sample_beta(10.00, 2.45, seed); // 60% to 95% + + double probability_of_dying = laplace_probability_nuclear_exchange_month * london_hit_conditional_on_russia_nuclear_weapon_usage * informed_actor_not_able_to_escape * proportion_which_die_if_bomb_drops_in_london; + return probability_of_dying; +} + +double probability_of_dying_eli(uint64_t* seed) +{ + double russia_nato_nuclear_exchange_in_next_month = sample_beta(1.30, 1182.99, seed); // .0001 to .003 + double london_hit_conditional = sample_beta(3.47, 8.97, seed); // 0.1 to 0.5 + double informed_actors_not_able_to_escape = sample_beta(2.73, 5.67, seed); // .1 to .6 + double proportion_which_die_if_bomb_drops_in_london = sample_beta(3.00, 1.46, seed); // 0.3 to 0.95; + + double probability_of_dying = russia_nato_nuclear_exchange_in_next_month * london_hit_conditional * informed_actors_not_able_to_escape * proportion_which_die_if_bomb_drops_in_london; + return probability_of_dying; +} + +double mixture(uint64_t* seed) +{ + double (*samplers[])(uint64_t*) = { probability_of_dying_nuno, probability_of_dying_eli }; + double weights[] = { 0.5, 0.5 }; + return sample_mixture(samplers, weights, 2, seed); +} + +int main() +{ + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 + + int n = 1000 * 1000; + + double* mixture_result = malloc(sizeof(double) * n); + for (int i = 0; i < n; i++) { + mixture_result[i] = mixture(seed); + } + + printf("mixture_result: [ "); + for (int i = 0; i < 9; i++) { + printf("%.6f, ", mixture_result[i]); + } + printf("... ]\n"); + + ci ci_90 = get_90_confidence_interval(mixture, seed); + printf("mean: %f\n", array_mean(mixture_result, n)); + printf("90%% confidence interval: [%f, %f]\n", ci_90.low, ci_90.high); + + free(seed); +} diff --git a/examples/more/08_nuclear_war/scratchpad/example b/examples/more/04_nuclear_war/scratchpad/example Binary files differ. diff --git a/examples/more/08_nuclear_war/scratchpad/example.c b/examples/more/04_nuclear_war/scratchpad/example.c diff --git a/examples/more/08_nuclear_war/scratchpad/makefile b/examples/more/04_nuclear_war/scratchpad/makefile diff --git a/examples/more/04_sample_from_cdf_simple/example b/examples/more/04_sample_from_cdf_simple/example Binary files differ. diff --git a/examples/more/04_sample_from_cdf_simple/example.c b/examples/more/04_sample_from_cdf_simple/example.c @@ -1,103 +0,0 @@ -#include "../../../squiggle.h" -#include "../../squiggle_more.h" -#include <math.h> -#include <stdint.h> -#include <stdio.h> -#include <stdlib.h> -#include <time.h> - -#define NUM_SAMPLES 1000000 - -// Example cdf -double cdf_uniform_0_1(double x) -{ - if (x < 0) { - return 0; - } else if (x > 1) { - return 1; - } else { - return x; - } -} - -double cdf_squared_0_1(double x) -{ - if (x < 0) { - return 0; - } else if (x > 1) { - return 1; - } else { - return x * x; - } -} - -double cdf_normal_0_1(double x) -{ - double mean = 0; - double std = 1; - return 0.5 * (1 + erf((x - mean) / (std * sqrt(2)))); // erf from math.h -} - -// Some testers -void test_inverse_cdf_double(char* cdf_name, double cdf_double(double)) -{ - struct box result = inverse_cdf_double(cdf_double, 0.5); - if (result.empty) { - printf("Inverse for %s not calculated\n", cdf_name); - exit(1); - } else { - printf("Inverse of %s at %f is: %f\n", cdf_name, 0.5, result.content); - } -} - -void test_and_time_sampler_double(char* cdf_name, double cdf_double(double), uint64_t* seed) -{ - printf("\nGetting some samples from %s:\n", cdf_name); - clock_t begin = clock(); - for (int i = 0; i < NUM_SAMPLES; i++) { - struct box sample = sampler_cdf_double(cdf_double, seed); - if (sample.empty) { - printf("Error in sampler function for %s", cdf_name); - } else { - // printf("%f\n", sample.content); - } - } - clock_t end = clock(); - double time_spent = (double)(end - begin) / CLOCKS_PER_SEC; - printf("Time spent: %f\n", time_spent); -} - -int main() -{ - // Test inverse cdf double - test_inverse_cdf_double("cdf_uniform_0_1", cdf_uniform_0_1); - test_inverse_cdf_double("cdf_squared_0_1", cdf_squared_0_1); - test_inverse_cdf_double("cdf_normal_0_1", cdf_normal_0_1); - - // Testing samplers - // set randomness seed - uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with 0 - - // Test double sampler - test_and_time_sampler_double("cdf_uniform_0_1", cdf_uniform_0_1, seed); - test_and_time_sampler_double("cdf_squared_0_1", cdf_squared_0_1, seed); - test_and_time_sampler_double("cdf_normal_0_1", cdf_normal_0_1, seed); - - // Get some normal samples using a previous approach - printf("\nGetting some samples from sample_unit_normal\n"); - - clock_t begin_2 = clock(); - double* normal_samples = malloc(NUM_SAMPLES * sizeof(double)); - for (int i = 0; i < NUM_SAMPLES; i++) { - normal_samples[i] = sample_unit_normal(seed); - // printf("%f\n", normal_sample); - } - - clock_t end_2 = clock(); - double time_spent_2 = (double)(end_2 - begin_2) / CLOCKS_PER_SEC; - printf("Time spent: %f\n", time_spent_2); - - free(seed); - return 0; -} diff --git a/examples/more/05_burn_10kg_fat/example b/examples/more/05_burn_10kg_fat/example Binary files differ. diff --git a/examples/more/05_burn_10kg_fat/example.c b/examples/more/05_burn_10kg_fat/example.c @@ -0,0 +1,49 @@ +#include "../../../squiggle.h" +#include "../../../squiggle_more.h" +#include <math.h> +#include <stdint.h> +#include <stdio.h> +#include <stdlib.h> + +double sample_minutes_per_day_jumping_rope_needed_to_burn_10kg(uint64_t* seed) +{ + double kcal_jumping_rope_minute = sample_to(15, 20, seed); + double kcal_jumping_rope_hour = kcal_jumping_rope_minute * 60; + + double kcal_in_kg_of_fat = 7700; + double num_kg_of_fat_to_lose = 10; + + double hours_jumping_rope_needed = kcal_in_kg_of_fat * num_kg_of_fat_to_lose / kcal_jumping_rope_hour; + + double days_until_end_of_year = 152; // as of 2023-08-01 + double hours_per_day = hours_jumping_rope_needed / days_until_end_of_year; + double minutes_per_day = hours_per_day * 60; + return minutes_per_day; +} + +int main() +{ + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 + + int n = 1000 * 1000; + + double* result = malloc(sizeof(double) * n); + for (int i = 0; i < n; i++) { + result[i] = sample_minutes_per_day_jumping_rope_needed_to_burn_10kg(seed); + } + + printf("## How many minutes per day do I have to jump rope to lose 10kg of fat by the end of the year?\n"); + printf("Mean: %f\n", array_mean(result, n)); + printf("A few samples: [ "); + for (int i = 0; i < 9; i++) { + printf("%.6f, ", result[i]); + } + printf("... ]\n"); + + ci ci_90 = get_90_confidence_interval(sample_minutes_per_day_jumping_rope_needed_to_burn_10kg, seed); + printf("90%% confidence interval: [%f, %f]\n", ci_90.low, ci_90.high); + + free(seed); +} diff --git a/examples/more/05_sample_from_cdf_beta/example b/examples/more/05_sample_from_cdf_beta/example Binary files differ. diff --git a/examples/more/05_sample_from_cdf_beta/example.c b/examples/more/05_sample_from_cdf_beta/example.c @@ -1,169 +0,0 @@ -#include "../../../squiggle.h" -#include "../../squiggle_more.h" -#include <math.h> -#include <stdint.h> -#include <stdio.h> -#include <stdlib.h> -#include <time.h> - -#define NUM_SAMPLES 10000 -#define STOP_BETA 1.0e-8 -#define TINY_BETA 1.0e-30 - -// Incomplete beta function -struct box incbeta(double a, double b, double x) -{ - // Descended from <https://github.com/codeplea/incbeta/blob/master/incbeta.c>, - // <https://codeplea.com/incomplete-beta-function-c> - // but modified to return a box struct and doubles instead of doubles. - // [ ] to do: add attribution in README - // Original code under this license: - /* - * zlib License - * - * Regularized Incomplete Beta Function - * - * Copyright (c) 2016, 2017 Lewis Van Winkle - * http://CodePlea.com - * - * This software is provided 'as-is', without any express or implied - * warranty. In no event will the authors be held liable for any damages - * arising from the use of this software. - * - * Permission is granted to anyone to use this software for any purpose, - * including commercial applications, and to alter it and redistribute it - * freely, subject to the following restrictions: - * - * 1. The origin of this software must not be misrepresented; you must not - * claim that you wrote the original software. If you use this software - * in a product, an acknowledgement in the product documentation would be - * appreciated but is not required. - * 2. Altered source versions must be plainly marked as such, and must not be - * misrepresented as being the original software. - * 3. This notice may not be removed or altered from any source distribution. - */ - if (x < 0.0 || x > 1.0) { - return PROCESS_ERROR("x out of bounds [0, 1], in function incbeta"); - } - - /*The continued fraction converges nicely for x < (a+1)/(a+b+2)*/ - if (x > (a + 1.0) / (a + b + 2.0)) { - struct box symmetric_incbeta = incbeta(b, a, 1.0 - x); - if (symmetric_incbeta.empty) { - return symmetric_incbeta; // propagate error - } else { - struct box result = { - .empty = 0, - .content = 1 - symmetric_incbeta.content - }; - return result; - } - } - - /*Find the first part before the continued fraction.*/ - const double lbeta_ab = lgamma(a) + lgamma(b) - lgamma(a + b); - const double front = exp(log(x) * a + log(1.0 - x) * b - lbeta_ab) / a; - - /*Use Lentz's algorithm to evaluate the continued fraction.*/ - double f = 1.0, c = 1.0, d = 0.0; - - int i, m; - for (i = 0; i <= 200; ++i) { - m = i / 2; - - double numerator; - if (i == 0) { - numerator = 1.0; /*First numerator is 1.0.*/ - } else if (i % 2 == 0) { - numerator = (m * (b - m) * x) / ((a + 2.0 * m - 1.0) * (a + 2.0 * m)); /*Even term.*/ - } else { - numerator = -((a + m) * (a + b + m) * x) / ((a + 2.0 * m) * (a + 2.0 * m + 1)); /*Odd term.*/ - } - - /*Do an iteration of Lentz's algorithm.*/ - d = 1.0 + numerator * d; - if (fabs(d) < TINY_BETA) - d = TINY_BETA; - d = 1.0 / d; - - c = 1.0 + numerator / c; - if (fabs(c) < TINY_BETA) - c = TINY_BETA; - - const double cd = c * d; - f *= cd; - - /*Check for stop.*/ - if (fabs(1.0 - cd) < STOP_BETA) { - struct box result = { - .empty = 0, - .content = front * (f - 1.0) - }; - return result; - } - } - - return PROCESS_ERROR("More loops needed, did not converge, in function incbeta"); -} - -struct box cdf_beta(double x) -{ - if (x < 0) { - struct box result = { .empty = 0, .content = 0 }; - return result; - } else if (x > 1) { - struct box result = { .empty = 0, .content = 1 }; - return result; - } else { - double successes = 1, failures = (2023 - 1945); - return incbeta(successes, failures, x); - } -} - -// Some testers -void test_inverse_cdf_box(char* cdf_name, struct box cdf_box(double)) -{ - struct box result = inverse_cdf_box(cdf_box, 0.5); - if (result.empty) { - printf("Inverse for %s not calculated\n", cdf_name); - exit(1); - } else { - printf("Inverse of %s at %f is: %f\n", cdf_name, 0.5, result.content); - } -} - -void test_and_time_sampler_box(char* cdf_name, struct box cdf_box(double), uint64_t* seed) -{ - printf("\nGetting some samples from %s:\n", cdf_name); - clock_t begin = clock(); - for (int i = 0; i < NUM_SAMPLES; i++) { - struct box sample = sampler_cdf_box(cdf_box, seed); - if (sample.empty) { - printf("Error in sampler function for %s", cdf_name); - } else { - // printf("%f\n", sample.content); - } - } - clock_t end = clock(); - double time_spent = (double)(end - begin) / CLOCKS_PER_SEC; - printf("Time spent: %f\n", time_spent); -} - -int main() -{ - // Test inverse cdf box - test_inverse_cdf_box("cdf_beta", cdf_beta); - - // Test box sampler - 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!! - // Partly this is because I am using a more general algorithm, - // which applies to any cdf - // But I am also using absurdly precise convergence conditions. - // This could be optimized. - - free(seed); - return 0; -} diff --git a/examples/more/06_nuclear_recovery/example b/examples/more/06_nuclear_recovery/example Binary files differ. diff --git a/examples/more/06_nuclear_recovery/example.c b/examples/more/06_nuclear_recovery/example.c @@ -0,0 +1,86 @@ +#include "../../../squiggle.h" +#include "../../../squiggle_more.h" +#include <math.h> +#include <stdint.h> +#include <stdio.h> +#include <stdlib.h> + +double yearly_probability_nuclear_collapse(double year, uint64_t* seed) +{ + double successes = 0; + double failures = (year - 1960); + return sample_laplace(successes, failures, seed); + // ^ can change to (successes + 1)/(trials + 2) + // to get a probability, + // rather than sampling from a distribution over probabilities. +} +double yearly_probability_nuclear_collapse_2023(uint64_t* seed) +{ + return yearly_probability_nuclear_collapse(2023, seed); +} + +double yearly_probability_nuclear_collapse_after_recovery(double year, double rebuilding_period_length_years, uint64_t* seed) +{ + // assumption: nuclear + double successes = 1.0; + double failures = (year - rebuilding_period_length_years - 1960 - 1); + return sample_laplace(successes, failures, seed); +} +double yearly_probability_nuclear_collapse_after_recovery_example(uint64_t* seed) +{ + double year = 2070; + double rebuilding_period_length_years = 30; + // So, there was a nuclear collapse in 2040, + // then a recovery period of 30 years + // and it's now 2070 + return yearly_probability_nuclear_collapse_after_recovery(year, rebuilding_period_length_years, seed); +} + +double yearly_probability_nuclear_collapse_after_recovery_antiinductive(uint64_t* seed) +{ + return yearly_probability_nuclear_collapse(2023, seed) / 2; +} + +int main() +{ + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 + + int num_samples = 1000000; + + // Before a first nuclear collapse + printf("## Before the first nuclear collapse\n"); + ci ci_90_2023 = get_90_confidence_interval(yearly_probability_nuclear_collapse_2023, seed); + printf("90%% confidence interval: [%f, %f]\n", ci_90_2023.low, ci_90_2023.high); + + double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * num_samples); + for (int i = 0; i < num_samples; i++) { + yearly_probability_nuclear_collapse_2023_samples[i] = yearly_probability_nuclear_collapse_2023(seed); + } + printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_2023_samples, num_samples)); + + // After the first nuclear collapse + printf("\n## After the first nuclear collapse\n"); + ci ci_90_2070 = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_example, seed); + printf("90%% confidence interval: [%f, %f]\n", ci_90_2070.low, ci_90_2070.high); + + double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * num_samples); + for (int i = 0; i < num_samples; i++) { + yearly_probability_nuclear_collapse_after_recovery_samples[i] = yearly_probability_nuclear_collapse_after_recovery_example(seed); + } + printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_after_recovery_samples, num_samples)); + + // After the first nuclear collapse (antiinductive) + printf("\n## After the first nuclear collapse (antiinductive)\n"); + ci ci_90_antiinductive = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_antiinductive, seed); + printf("90%% confidence interval: [%f, %f]\n", ci_90_antiinductive.low, ci_90_antiinductive.high); + + double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * num_samples); + for (int i = 0; i < num_samples; i++) { + yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples[i] = yearly_probability_nuclear_collapse_after_recovery_antiinductive(seed); + } + printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples, num_samples)); + + free(seed); +} diff --git a/examples/more/07_algebra/example b/examples/more/07_algebra/example Binary files differ. diff --git a/examples/more/07_algebra/example.c b/examples/more/07_algebra/example.c @@ -0,0 +1,27 @@ +#include "../../../squiggle.h" +#include "../../../squiggle_more.h" +#include <math.h> +#include <stdint.h> +#include <stdio.h> +#include <stdlib.h> + +int main() +{ + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 + + normal_params n1 = { .mean = 1.0, .std = 3.0 }; + normal_params n2 = { .mean = 2.0, .std = 4.0 }; + normal_params sn = algebra_sum_normals(n1, n2); + printf("The sum of Normal(%f, %f) and Normal(%f, %f) is Normal(%f, %f)\n", + n1.mean, n1.std, n2.mean, n2.std, sn.mean, sn.std); + + lognormal_params ln1 = { .logmean = 1.0, .logstd = 3.0 }; + lognormal_params ln2 = { .logmean = 2.0, .logstd = 4.0 }; + lognormal_params sln = algebra_product_lognormals(ln1, ln2); + printf("The product of Lognormal(%f, %f) and Lognormal(%f, %f) is Lognormal(%f, %f)\n", + ln1.logmean, ln1.logstd, ln2.logmean, ln2.logstd, sln.logmean, sln.logstd); + + free(seed); +} diff --git a/examples/more/07_ci_beta/example.c b/examples/more/07_ci_beta/example.c @@ -1,23 +0,0 @@ -#include "../../../squiggle.h" -#include "../../squiggle_more.h" -#include <stdint.h> -#include <stdio.h> -#include <stdlib.h> - -// Estimate functions -double beta_1_2_sampler(uint64_t* seed) -{ - return sample_beta(1, 2.0, seed); -} - -int main() -{ - // set randomness seed - uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with 0 - - ci beta_1_2_ci_90 = get_90_confidence_interval(beta_1_2_sampler, seed); - printf("90%% confidence interval of beta(1,2) is [%f, %f]\n", beta_1_2_ci_90.low, beta_1_2_ci_90.high); - - free(seed); -} diff --git a/examples/more/08_algebra_and_conversion/example b/examples/more/08_algebra_and_conversion/example Binary files differ. diff --git a/examples/more/08_algebra_and_conversion/example.c b/examples/more/08_algebra_and_conversion/example.c @@ -0,0 +1,34 @@ +#include "../../../squiggle.h" +#include "../../../squiggle_more.h" +#include <math.h> +#include <stdint.h> +#include <stdio.h> +#include <stdlib.h> + +int main() +{ + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 + + // Convert to 90% confidence interval form and back + lognormal_params ln1 = { .logmean = 1.0, .logstd = 3.0 }; + ci ln1_ci = convert_lognormal_params_to_ci(ln1); + printf("The 90%% confidence interval of Lognormal(%f, %f) is [%f, %f]\n", + ln1.logmean, ln1.logstd, + ln1_ci.low, ln1_ci.high); + lognormal_params ln1_params2 = convert_ci_to_lognormal_params(ln1_ci); + printf("The lognormal which has 90%% confidence interval [%f, %f] is Lognormal(%f, %f)\n", + ln1_ci.low, ln1_ci.high, + ln1_params2.logmean, ln1_params2.logstd); + + lognormal_params ln2 = convert_ci_to_lognormal_params((ci) { .low = 1, .high = 10 }); + lognormal_params ln3 = convert_ci_to_lognormal_params((ci) { .low = 5, .high = 50 }); + + lognormal_params sln = algebra_product_lognormals(ln2, ln3); + ci sln_ci = convert_lognormal_params_to_ci(sln); + + printf("Result of some lognormal products: to(%f, %f)\n", sln_ci.low, sln_ci.high); + + free(seed); +} diff --git a/examples/more/08_nuclear_war/example b/examples/more/08_nuclear_war/example Binary files differ. diff --git a/examples/more/08_nuclear_war/example.c b/examples/more/08_nuclear_war/example.c @@ -1,69 +0,0 @@ -#include "../../../squiggle.h" -#include "../../squiggle_more.h" -#include <math.h> -#include <stdint.h> -#include <stdio.h> -#include <stdlib.h> - -double probability_of_dying_nuno(uint64_t* seed) -{ - double first_year_russian_nuclear_weapons = 1953; - double current_year = 2022; - double laplace_probability_nuclear_exchange_year = sample_beta(1, current_year - first_year_russian_nuclear_weapons + 1, seed); - double laplace_probability_nuclear_exchange_month = 1 - pow(1 - laplace_probability_nuclear_exchange_year, (1.0 / 12.0)); - - double london_hit_conditional_on_russia_nuclear_weapon_usage = sample_beta(7.67, 69.65, seed); - // I.e., a beta distribution with a range of 0.05 to 0.16 into: https://nunosempere.com/blog/2023/03/15/fit-beta/ - // 0.05 were my estimate and Samotsvety's estimate in March 2022, respectively: - // https://forum.effectivealtruism.org/posts/KRFXjCqqfGQAYirm5/samotsvety-nuclear-risk-forecasts-march-2022#Nu_o_Sempere - double informed_actor_not_able_to_escape = sample_beta(3.26212166586967, 3.26228162008564, seed); - // 0.2 to 0.8, i.e., 20% to 80%, again using the previous tool - double proportion_which_die_if_bomb_drops_in_london = sample_beta(10.00, 2.45, seed); // 60% to 95% - - double probability_of_dying = laplace_probability_nuclear_exchange_month * london_hit_conditional_on_russia_nuclear_weapon_usage * informed_actor_not_able_to_escape * proportion_which_die_if_bomb_drops_in_london; - return probability_of_dying; -} - -double probability_of_dying_eli(uint64_t* seed) -{ - double russia_nato_nuclear_exchange_in_next_month = sample_beta(1.30, 1182.99, seed); // .0001 to .003 - double london_hit_conditional = sample_beta(3.47, 8.97, seed); // 0.1 to 0.5 - double informed_actors_not_able_to_escape = sample_beta(2.73, 5.67, seed); // .1 to .6 - double proportion_which_die_if_bomb_drops_in_london = sample_beta(3.00, 1.46, seed); // 0.3 to 0.95; - - double probability_of_dying = russia_nato_nuclear_exchange_in_next_month * london_hit_conditional * informed_actors_not_able_to_escape * proportion_which_die_if_bomb_drops_in_london; - return probability_of_dying; -} - -double mixture(uint64_t* seed) -{ - double (*samplers[])(uint64_t*) = { probability_of_dying_nuno, probability_of_dying_eli }; - double weights[] = { 0.5, 0.5 }; - return sample_mixture(samplers, weights, 2, seed); -} - -int main() -{ - // set randomness seed - uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with 0 - - int n = 1000 * 1000; - - double* mixture_result = malloc(sizeof(double) * n); - for (int i = 0; i < n; i++) { - mixture_result[i] = mixture(seed); - } - - printf("mixture_result: [ "); - for (int i = 0; i < 9; i++) { - printf("%.6f, ", mixture_result[i]); - } - printf("... ]\n"); - - ci ci_90 = get_90_confidence_interval(mixture, seed); - printf("mean: %f\n", array_mean(mixture_result, n)); - printf("90%% confidence interval: [%f, %f]\n", ci_90.low, ci_90.high); - - free(seed); -} diff --git a/examples/more/09_burn_10kg_fat/example b/examples/more/09_burn_10kg_fat/example Binary files differ. diff --git a/examples/more/09_burn_10kg_fat/example.c b/examples/more/09_burn_10kg_fat/example.c @@ -1,49 +0,0 @@ -#include "../../../squiggle.h" -#include "../../squiggle_more.h" -#include <math.h> -#include <stdint.h> -#include <stdio.h> -#include <stdlib.h> - -double sample_minutes_per_day_jumping_rope_needed_to_burn_10kg(uint64_t* seed) -{ - double kcal_jumping_rope_minute = sample_to(15, 20, seed); - double kcal_jumping_rope_hour = kcal_jumping_rope_minute * 60; - - double kcal_in_kg_of_fat = 7700; - double num_kg_of_fat_to_lose = 10; - - double hours_jumping_rope_needed = kcal_in_kg_of_fat * num_kg_of_fat_to_lose / kcal_jumping_rope_hour; - - double days_until_end_of_year = 152; // as of 2023-08-01 - double hours_per_day = hours_jumping_rope_needed / days_until_end_of_year; - double minutes_per_day = hours_per_day * 60; - return minutes_per_day; -} - -int main() -{ - // set randomness seed - uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with 0 - - int n = 1000 * 1000; - - double* result = malloc(sizeof(double) * n); - for (int i = 0; i < n; i++) { - result[i] = sample_minutes_per_day_jumping_rope_needed_to_burn_10kg(seed); - } - - printf("## How many minutes per day do I have to jump rope to lose 10kg of fat by the end of the year?\n"); - printf("Mean: %f\n", array_mean(result, n)); - printf("A few samples: [ "); - for (int i = 0; i < 9; i++) { - printf("%.6f, ", result[i]); - } - printf("... ]\n"); - - ci ci_90 = get_90_confidence_interval(sample_minutes_per_day_jumping_rope_needed_to_burn_10kg, seed); - printf("90%% confidence interval: [%f, %f]\n", ci_90.low, ci_90.high); - - free(seed); -} diff --git a/examples/more/09_ergonomic_algebra/example b/examples/more/09_ergonomic_algebra/example Binary files differ. diff --git a/examples/more/09_ergonomic_algebra/example.c b/examples/more/09_ergonomic_algebra/example.c @@ -0,0 +1,27 @@ +#include "../../../squiggle.h" +#include "../../../squiggle_more.h" +#include <math.h> +#include <stdint.h> +#include <stdio.h> +#include <stdlib.h> + +#define ln lognormal_params +#define to(...) convert_ci_to_lognormal_params((ci)__VA_ARGS__) +#define from(...) convert_lognormal_params_to_ci((ln)__VA_ARGS__) +#define times(a, b) algebra_product_lognormals(a, b) + +int main() +{ + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 + + ln a = to({ .low = 1, .high = 10 }); + ln b = to({ .low = 5, .high = 500 }); + ln c = times(a, b); + + printf("Result: to(%f, %f)\n", from(c).low, from(c).high); + printf("One sample from it is: %f\n", sample_lognormal(c.logmean, c.logstd, seed)); + + free(seed); +} diff --git a/examples/more/10_nuclear_recovery/example b/examples/more/10_nuclear_recovery/example Binary files differ. diff --git a/examples/more/10_nuclear_recovery/example.c b/examples/more/10_nuclear_recovery/example.c @@ -1,86 +0,0 @@ -#include "../../../squiggle.h" -#include "../../squiggle_more.h" -#include <math.h> -#include <stdint.h> -#include <stdio.h> -#include <stdlib.h> - -double yearly_probability_nuclear_collapse(double year, uint64_t* seed) -{ - double successes = 0; - double failures = (year - 1960); - return sample_laplace(successes, failures, seed); - // ^ can change to (successes + 1)/(trials + 2) - // to get a probability, - // rather than sampling from a distribution over probabilities. -} -double yearly_probability_nuclear_collapse_2023(uint64_t* seed) -{ - return yearly_probability_nuclear_collapse(2023, seed); -} - -double yearly_probability_nuclear_collapse_after_recovery(double year, double rebuilding_period_length_years, uint64_t* seed) -{ - // assumption: nuclear - double successes = 1.0; - double failures = (year - rebuilding_period_length_years - 1960 - 1); - return sample_laplace(successes, failures, seed); -} -double yearly_probability_nuclear_collapse_after_recovery_example(uint64_t* seed) -{ - double year = 2070; - double rebuilding_period_length_years = 30; - // So, there was a nuclear collapse in 2040, - // then a recovery period of 30 years - // and it's now 2070 - return yearly_probability_nuclear_collapse_after_recovery(year, rebuilding_period_length_years, seed); -} - -double yearly_probability_nuclear_collapse_after_recovery_antiinductive(uint64_t* seed) -{ - return yearly_probability_nuclear_collapse(2023, seed) / 2; -} - -int main() -{ - // set randomness seed - uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with 0 - - int num_samples = 1000000; - - // Before a first nuclear collapse - printf("## Before the first nuclear collapse\n"); - ci ci_90_2023 = get_90_confidence_interval(yearly_probability_nuclear_collapse_2023, seed); - printf("90%% confidence interval: [%f, %f]\n", ci_90_2023.low, ci_90_2023.high); - - double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * num_samples); - for (int i = 0; i < num_samples; i++) { - yearly_probability_nuclear_collapse_2023_samples[i] = yearly_probability_nuclear_collapse_2023(seed); - } - printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_2023_samples, num_samples)); - - // After the first nuclear collapse - printf("\n## After the first nuclear collapse\n"); - ci ci_90_2070 = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_example, seed); - printf("90%% confidence interval: [%f, %f]\n", ci_90_2070.low, ci_90_2070.high); - - double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * num_samples); - for (int i = 0; i < num_samples; i++) { - yearly_probability_nuclear_collapse_after_recovery_samples[i] = yearly_probability_nuclear_collapse_after_recovery_example(seed); - } - printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_after_recovery_samples, num_samples)); - - // After the first nuclear collapse (antiinductive) - printf("\n## After the first nuclear collapse (antiinductive)\n"); - ci ci_90_antiinductive = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_antiinductive, seed); - printf("90%% confidence interval: [%f, %f]\n", ci_90_antiinductive.low, ci_90_antiinductive.high); - - double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * num_samples); - for (int i = 0; i < num_samples; i++) { - yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples[i] = yearly_probability_nuclear_collapse_after_recovery_antiinductive(seed); - } - printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples, num_samples)); - - free(seed); -} diff --git a/examples/more/10_twitter_thread_example/example b/examples/more/10_twitter_thread_example/example Binary files differ. diff --git a/examples/more/10_twitter_thread_example/example.c b/examples/more/10_twitter_thread_example/example.c @@ -0,0 +1,48 @@ +#include "../../../squiggle.h" +#include "../../../squiggle_more.h" +#include <stdint.h> +#include <stdio.h> +#include <stdlib.h> + +double sample_0(uint64_t* seed) +{ + return 0; +} + +double sample_1(uint64_t* seed) +{ + return 1; +} + +double sample_normal_mean_1_std_2(uint64_t* seed) +{ + return sample_normal(1, 2, seed); +} + +double sample_1_to_3(uint64_t* seed) +{ + return sample_to(1, 3, seed); +} + +int main() +{ + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 + + int n_dists = 4; + double weights[] = { 1, 2, 3, 4 }; + double (*samplers[])(uint64_t*) = { + sample_0, + sample_1, + sample_normal_mean_1_std_2, + sample_1_to_3 + }; + + int n_samples = 10; + for (int i = 0; i < n_samples; i++) { + printf("Sample #%d: %f\n", i, sample_mixture(samplers, weights, n_dists, seed)); + } + + free(seed); +} diff --git a/examples/more/11_algebra/example b/examples/more/11_algebra/example Binary files differ. diff --git a/examples/more/11_algebra/example.c b/examples/more/11_algebra/example.c @@ -1,27 +0,0 @@ -#include "../../../squiggle.h" -#include "../../squiggle_more.h" -#include <math.h> -#include <stdint.h> -#include <stdio.h> -#include <stdlib.h> - -int main() -{ - // set randomness seed - uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with 0 - - normal_params n1 = { .mean = 1.0, .std = 3.0 }; - normal_params n2 = { .mean = 2.0, .std = 4.0 }; - normal_params sn = algebra_sum_normals(n1, n2); - printf("The sum of Normal(%f, %f) and Normal(%f, %f) is Normal(%f, %f)\n", - n1.mean, n1.std, n2.mean, n2.std, sn.mean, sn.std); - - lognormal_params ln1 = { .logmean = 1.0, .logstd = 3.0 }; - lognormal_params ln2 = { .logmean = 2.0, .logstd = 4.0 }; - lognormal_params sln = algebra_product_lognormals(ln1, ln2); - printf("The product of Lognormal(%f, %f) and Lognormal(%f, %f) is Lognormal(%f, %f)\n", - ln1.logmean, ln1.logstd, ln2.logmean, ln2.logstd, sln.logmean, sln.logstd); - - free(seed); -} diff --git a/examples/more/12_algebra_and_conversion/example b/examples/more/12_algebra_and_conversion/example Binary files differ. diff --git a/examples/more/12_algebra_and_conversion/example.c b/examples/more/12_algebra_and_conversion/example.c @@ -1,34 +0,0 @@ -#include "../../../squiggle.h" -#include "../../squiggle_more.h" -#include <math.h> -#include <stdint.h> -#include <stdio.h> -#include <stdlib.h> - -int main() -{ - // set randomness seed - uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with 0 - - // Convert to 90% confidence interval form and back - lognormal_params ln1 = { .logmean = 1.0, .logstd = 3.0 }; - ci ln1_ci = convert_lognormal_params_to_ci(ln1); - printf("The 90%% confidence interval of Lognormal(%f, %f) is [%f, %f]\n", - ln1.logmean, ln1.logstd, - ln1_ci.low, ln1_ci.high); - lognormal_params ln1_params2 = convert_ci_to_lognormal_params(ln1_ci); - printf("The lognormal which has 90%% confidence interval [%f, %f] is Lognormal(%f, %f)\n", - ln1_ci.low, ln1_ci.high, - ln1_params2.logmean, ln1_params2.logstd); - - lognormal_params ln2 = convert_ci_to_lognormal_params((ci) { .low = 1, .high = 10 }); - lognormal_params ln3 = convert_ci_to_lognormal_params((ci) { .low = 5, .high = 50 }); - - lognormal_params sln = algebra_product_lognormals(ln2, ln3); - ci sln_ci = convert_lognormal_params_to_ci(sln); - - printf("Result of some lognormal products: to(%f, %f)\n", sln_ci.low, sln_ci.high); - - free(seed); -} diff --git a/examples/more/13_ergonomic_algebra/example b/examples/more/13_ergonomic_algebra/example Binary files differ. diff --git a/examples/more/13_ergonomic_algebra/example.c b/examples/more/13_ergonomic_algebra/example.c @@ -1,27 +0,0 @@ -#include "../../../squiggle.h" -#include "../../squiggle_more.h" -#include <math.h> -#include <stdint.h> -#include <stdio.h> -#include <stdlib.h> - -#define ln lognormal_params -#define to(...) convert_ci_to_lognormal_params((ci)__VA_ARGS__) -#define from(...) convert_lognormal_params_to_ci((ln)__VA_ARGS__) -#define times(a, b) algebra_product_lognormals(a, b) - -int main() -{ - // set randomness seed - uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with 0 - - ln a = to({ .low = 1, .high = 10 }); - ln b = to({ .low = 5, .high = 500 }); - ln c = times(a, b); - - printf("Result: to(%f, %f)\n", from(c).low, from(c).high); - printf("One sample from it is: %f\n", sample_lognormal(c.logmean, c.logstd, seed)); - - free(seed); -} diff --git a/examples/more/14_twitter_thread_example/example b/examples/more/14_twitter_thread_example/example Binary files differ. diff --git a/examples/more/14_twitter_thread_example/example.c b/examples/more/14_twitter_thread_example/example.c @@ -1,48 +0,0 @@ -#include "../../../squiggle.h" -#include "../../squiggle_more.h" -#include <stdint.h> -#include <stdio.h> -#include <stdlib.h> - -double sample_0(uint64_t* seed) -{ - return 0; -} - -double sample_1(uint64_t* seed) -{ - return 1; -} - -double sample_normal_mean_1_std_2(uint64_t* seed) -{ - return sample_normal(1, 2, seed); -} - -double sample_1_to_3(uint64_t* seed) -{ - return sample_to(1, 3, seed); -} - -int main() -{ - // set randomness seed - uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with 0 - - int n_dists = 4; - double weights[] = { 1, 2, 3, 4 }; - double (*samplers[])(uint64_t*) = { - sample_0, - sample_1, - sample_normal_mean_1_std_2, - sample_1_to_3 - }; - - int n_samples = 10; - for (int i = 0; i < n_samples; i++) { - printf("Sample #%d: %f\n", i, sample_mixture(samplers, weights, n_dists, seed)); - } - - free(seed); -} diff --git a/examples/more/makefile b/examples/more/makefile @@ -1,58 +1,102 @@ # Interface: -# make -# make build -# make format -# make run +# make all +# make format-all +# make run-all +# make one DIR=06_nuclear_recovery +# make format-one DIR=06_nuclear_recovery +# make run-one DIR=06_nuclear_recovery +# make time-linux-one DIR=06_nuclear_recovery +# make profile-one DIR=06_nuclear_recovery # Compiler -CC=gcc # required for nested functions +CC=gcc # CC=tcc # <= faster compilation # Main file -SRC=example.c ../../squiggle.c ../../squiggle_more.c -OUTPUT=./example +SRC=example.c +OUTPUT=example ## Dependencies -OPENMP=-fopenmp +SQUIGGLE=../../squiggle.c +SQUIGGLE_MORE=../../squiggle_more.c MATH=-lm -DEPENDENCIES=$(MATH) $(OPENMP) -# OPENMP=-fopenmp +OPENMP=-fopenmp +DEPS=$(SQUIGGLE) $(SQUIGGLE_MORE) $(MATH) $(OPENMP) ## Flags DEBUG= #'-g' -STANDARD=-std=c99 ## gnu99 allows for nested functions. -EXTENSIONS= #-fnested-functions +STANDARD=-std=c99 WARNINGS=-Wall -OPTIMIZED=-O3#-Ofast -CFLAGS=$(DEBUG) $(STANDARD) $(EXTENSIONS) $(WARNINGS) $(OPTIMIZED) +OPTIMIZED=-O3 #-Ofast ## Formatter STYLE_BLUEPRINT=webkit FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT) -## make build -build: $(SRC) - # gcc -std=gnu99 example.c -lm -o example - $(CC) $(CFLAGS) $(SRC) $(DEPENDENCIES) -o $(OUTPUT) +## make all +all: + $(CC) $(OPTIMIZED) $(DEBUG) 00_example_template/$(SRC) $(DEPS) -o 00_example_template/$(OUTPUT) + $(CC) $(OPTIMIZED) $(DEBUG) 01_sample_from_cdf/$(SRC) $(DEPS) -o 01_sample_from_cdf/$(OUTPUT) + $(CC) $(OPTIMIZED) $(DEBUG) 02_sample_from_cdf_beta/$(SRC) $(DEPS) -o 02_sample_from_cdf_beta/$(OUTPUT) + $(CC) $(OPTIMIZED) $(DEBUG) 03_ci_beta/$(SRC) $(DEPS) -o 03_ci_beta/$(OUTPUT) + $(CC) $(OPTIMIZED) $(DEBUG) 04_nuclear_war/$(SRC) $(DEPS) -o 04_nuclear_war/$(OUTPUT) + $(CC) $(OPTIMIZED) $(DEBUG) 05_burn_10kg_fat/$(SRC) $(DEPS) -o 05_burn_10kg_fat/$(OUTPUT) + $(CC) $(OPTIMIZED) $(DEBUG) 06_nuclear_recovery/$(SRC) $(DEPS) -o 06_nuclear_recovery/$(OUTPUT) + $(CC) $(OPTIMIZED) $(DEBUG) 07_algebra/$(SRC) $(DEPS) -o 07_algebra/$(OUTPUT) + $(CC) $(OPTIMIZED) $(DEBUG) 08_algebra_and_conversion/$(SRC) $(DEPS) -o 08_algebra_and_conversion/$(OUTPUT) + $(CC) $(OPTIMIZED) $(DEBUG) 09_ergonomic_algebra/$(SRC) $(DEPS) -o 09_ergonomic_algebra/$(OUTPUT) + $(CC) $(OPTIMIZED) $(DEBUG) 10_twitter_thread_example/$(SRC) $(DEPS) -o 10_twitter_thread_example/$(OUTPUT) -format: $(SRC) - $(FORMATTER) $(SRC) +format-all: + $(FORMATTER) 00_example_template/$(SRC) + $(FORMATTER) 01_sample_from_cdf/$(SRC) + $(FORMATTER) 02_sample_from_cdf_beta/$(SRC) + $(FORMATTER) 03_ci_beta/$(SRC) + $(FORMATTER) 04_nuclear_war/$(SRC) + $(FORMATTER) 05_burn_10kg_fat/$(SRC) + $(FORMATTER) 06_nuclear_recovery/$(SRC) + $(FORMATTER) 07_algebra/$(SRC) + $(FORMATTER) 08_algebra_and_conversion/$(SRC) + $(FORMATTER) 09_ergonomic_algebra/$(SRC) + $(FORMATTER) 10_twitter_thread_example/$(SRC) -run: $(SRC) $(OUTPUT) - ./$(OUTPUT) && echo +run-all: + 00_example_template/$(OUTPUT) + 01_sample_from_cdf/$(OUTPUT) + 02_sample_from_cdf_beta/$(OUTPUT) + 03_ci_beta/$(OUTPUT) + 04_nuclear_war/$(OUTPUT) + 05_burn_10kg_fat/$(OUTPUT) + 06_nuclear_recovery/$(OUTPUT) + 07_algebra/$(OUTPUT) + 08_algebra_and_conversion/$(OUTPUT) + 09_ergonomic_algebra/$(OUTPUT) + 10_twitter_thread_example/$(OUTPUT) -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 +## make one DIR=06_nuclear_recovery +one: $(DIR)/$(SRC) + $(CC) $(OPTIMIZED) $(DEBUG) $(DIR)/$(SRC) $(DEPS) -o $(DIR)/$(OUTPUT) -## Profiling +## make format-one DIR=06_nuclear_recovery +format-one: $(DIR)/$(SRC) + $(FORMATTER) $(DIR)/$(SRC) + +## make run-one DIR=06_nuclear_recovery +run-one: $(DIR)/$(OUTPUT) + $(DIR)/$(OUTPUT) && echo + +## make time-linux-one DIR=06_nuclear_recovery +time-linux-one: $(DIR)/$(OUTPUT) + @echo "Requires /bin/time, found on GNU/Linux systems" && echo + @echo "Running 100x and taking avg time $(DIR)/$(OUTPUT)" + @t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do $(DIR)/$(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 -profile-linux: +## e.g., make profile-linux-one DIR=06_nuclear_recovery +profile-linux-one: 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) + $(CC) $(OPTIMIZED) $(DEBUG) $(DIR)/$(SRC) $(DEPS) -o $(DIR)/$(OUTPUT) + # $(CC) $(SRC) $(DEPS) -o $(OUTPUT) + sudo perf record $(DIR)/$(OUTPUT) sudo perf report rm perf.data