squiggle.c

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

commit 6387c0df703c3039a72c51d1e2dfca3d68202b9d
parent 61851a321a1d1f82a62eaf8ee10bf7cf52e060d3
Author: NunoSempere <nuno.sempere@protonmail.com>
Date:   Sat, 18 Nov 2023 21:10:21 +0000

fix: reorder headers to fix compilation error

Diffstat:
Mexamples/01_one_sample/example.c | 13+++++++------
Mexamples/03_gcc_nested_function/example.c | 50+++++++++++++++++++++++++-------------------------
Mexamples/04_sample_from_cdf_simple/example | 0
Mexamples/04_sample_from_cdf_simple/example.c | 2+-
Mexamples/05_sample_from_cdf_beta/example | 0
Mexamples/06_gamma_beta/example.c | 33++++++++++++++++-----------------
Mexamples/07_ci_beta/example | 0
Mexamples/07_ci_beta/example.c | 13+++++++------
Mexamples/08_nuclear_war/example | 0
Mexamples/09_burn_10kg_fat/example | 0
Mexamples/09_burn_10kg_fat/example.c | 2--
Mexamples/10_nuclear_recovery/example | 0
Mexamples/10_nuclear_recovery/example.c | 66+++++++++++++++++++++++++++++++++++-------------------------------
Mexamples/11_algebra/example | 0
Mexamples/12_algebra_and_conversion/example | 0
Mexamples/12_algebra_and_conversion/example.c | 14+++++++-------
Mexamples/13_ergonomic_algebra/example | 0
Mexamples/13_ergonomic_algebra/example.c | 10+++++-----
Mexamples/14_twitter_thread_example/example | 0
Mexamples/14_twitter_thread_example/example.c | 22+++++++++++++---------
Aexamples/15_plotting-scratchpad/makefile | 3+++
Mexamples/16_100_lognormal_samples/example.c | 21+++++++++++----------
Mextra.c | 3+--
23 files changed, 131 insertions(+), 121 deletions(-)

diff --git a/examples/01_one_sample/example.c b/examples/01_one_sample/example.c @@ -1,7 +1,7 @@ #include "../../squiggle.h" #include <stdint.h> -#include <stdlib.h> #include <stdio.h> +#include <stdlib.h> // Estimate functions double sample_0(uint64_t* seed) @@ -24,10 +24,11 @@ double sample_many(uint64_t* seed) return sample_to(2, 10, seed); } -int main(){ +int main() +{ // set randomness seed - uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with 0 + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 double p_a = 0.8; double p_b = 0.5; @@ -38,6 +39,6 @@ int main(){ double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many }; double result_one = sample_mixture(samplers, weights, n_dists, seed); - printf("result_one: %f\n", result_one); - free(seed); + printf("result_one: %f\n", result_one); + free(seed); } diff --git a/examples/03_gcc_nested_function/example.c b/examples/03_gcc_nested_function/example.c @@ -1,38 +1,38 @@ +#include "../../squiggle.h" #include <stdint.h> -#include <stdlib.h> #include <stdio.h> -#include "../../squiggle.h" +#include <stdlib.h> -int main(){ +int main() +{ // set randomness seed - uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with 0 + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 double p_a = 0.8; double p_b = 0.5; double p_c = p_a * p_b; int n_dists = 4; - - double sample_0(uint64_t* seed){ return 0; } - double sample_1(uint64_t* 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); } - + + double sample_0(uint64_t * seed) { return 0; } + double sample_1(uint64_t * 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); } + double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many }; - double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 }; + double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 }; - int n_samples = 1000000; - double* result_many = (double *) malloc(n_samples * sizeof(double)); - for(int i=0; i<n_samples; i++){ - result_many[i] = sample_mixture(samplers, weights, n_dists, seed); - } - - printf("result_many: ["); - for(int i=0; i<100; i++){ - printf("%.2f, ", result_many[i]); - } - printf("]\n"); - free(seed); -} + int n_samples = 1000000; + double* result_many = (double*)malloc(n_samples * sizeof(double)); + for (int i = 0; i < n_samples; i++) { + result_many[i] = sample_mixture(samplers, weights, n_dists, seed); + } + printf("result_many: ["); + for (int i = 0; i < 100; i++) { + printf("%.2f, ", result_many[i]); + } + printf("]\n"); + free(seed); +} diff --git a/examples/04_sample_from_cdf_simple/example b/examples/04_sample_from_cdf_simple/example Binary files differ. diff --git a/examples/04_sample_from_cdf_simple/example.c b/examples/04_sample_from_cdf_simple/example.c @@ -88,7 +88,7 @@ int main() printf("\nGetting some samples from sample_unit_normal\n"); clock_t begin_2 = clock(); - double* normal_samples = malloc(NUM_SAMPLES * sizeof(double)); + 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); diff --git a/examples/05_sample_from_cdf_beta/example b/examples/05_sample_from_cdf_beta/example Binary files differ. diff --git a/examples/06_gamma_beta/example.c b/examples/06_gamma_beta/example.c @@ -11,8 +11,8 @@ int main() uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with 0 - int n = 1000 * 1000; - /* + int n = 1000 * 1000; + /* for (int i = 0; i < n; i++) { double gamma_0 = sample_gamma(0.0, seed); // printf("sample_gamma(0.0): %f\n", gamma_0); @@ -20,26 +20,25 @@ int main() printf("\n"); */ - double* gamma_1_array = malloc(sizeof(double) * n); + double* gamma_1_array = malloc(sizeof(double) * n); for (int i = 0; i < n; i++) { double gamma_1 = sample_gamma(1.0, seed); // printf("sample_gamma(1.0): %f\n", gamma_1); - gamma_1_array[i] = gamma_1; + gamma_1_array[i] = gamma_1; } - printf("gamma(1) summary statistics = mean: %f, std: %f\n", array_mean(gamma_1_array, n), array_std(gamma_1_array, n)); - free(gamma_1_array); - printf("\n"); - - double* beta_1_2_array = malloc(sizeof(double) * n); - for (int i = 0; i < n; i++) { + printf("gamma(1) summary statistics = mean: %f, std: %f\n", array_mean(gamma_1_array, n), array_std(gamma_1_array, n)); + free(gamma_1_array); + printf("\n"); + + double* beta_1_2_array = malloc(sizeof(double) * n); + for (int i = 0; i < n; i++) { double beta_1_2 = sample_beta(1, 2.0, seed); // printf("sample_beta(1.0, 2.0): %f\n", beta_1_2); - beta_1_2_array[i] = beta_1_2; + beta_1_2_array[i] = beta_1_2; } - printf("beta(1,2) summary statistics: mean: %f, std: %f\n", array_mean(beta_1_2_array, n), array_std(beta_1_2_array, n)); - free(beta_1_2_array); - printf("\n"); - - free(seed); -} + printf("beta(1,2) summary statistics: mean: %f, std: %f\n", array_mean(beta_1_2_array, n), array_std(beta_1_2_array, n)); + free(beta_1_2_array); + printf("\n"); + free(seed); +} diff --git a/examples/07_ci_beta/example b/examples/07_ci_beta/example Binary files differ. diff --git a/examples/07_ci_beta/example.c b/examples/07_ci_beta/example.c @@ -5,8 +5,9 @@ #include <stdlib.h> // Estimate functions -double beta_1_2_sampler(uint64_t* seed){ - return sample_beta(1, 2.0, seed); +double beta_1_2_sampler(uint64_t* seed) +{ + return sample_beta(1, 2.0, seed); } int main() @@ -15,8 +16,8 @@ int main() 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); + 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/08_nuclear_war/example b/examples/08_nuclear_war/example Binary files differ. diff --git a/examples/09_burn_10kg_fat/example b/examples/09_burn_10kg_fat/example Binary files differ. diff --git a/examples/09_burn_10kg_fat/example.c b/examples/09_burn_10kg_fat/example.c @@ -45,7 +45,5 @@ int main() 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/10_nuclear_recovery/example b/examples/10_nuclear_recovery/example Binary files differ. diff --git a/examples/10_nuclear_recovery/example.c b/examples/10_nuclear_recovery/example.c @@ -7,34 +7,38 @@ 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 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_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(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_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; +double yearly_probability_nuclear_collapse_after_recovery_antiinductive(uint64_t* seed) +{ + return yearly_probability_nuclear_collapse(2023, seed) / 2; } int main() @@ -45,8 +49,8 @@ int main() int num_samples = 1000000; - // Before a first nuclear collapse - printf("## Before the first nuclear collapse\n"); + // 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); @@ -54,10 +58,10 @@ int main() 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)); + 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"); + // 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); @@ -65,10 +69,10 @@ int main() 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)); + 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"); + // 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); @@ -76,7 +80,7 @@ int main() 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)); + printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples, num_samples)); free(seed); } diff --git a/examples/11_algebra/example b/examples/11_algebra/example Binary files differ. diff --git a/examples/12_algebra_and_conversion/example b/examples/12_algebra_and_conversion/example Binary files differ. diff --git a/examples/12_algebra_and_conversion/example.c b/examples/12_algebra_and_conversion/example.c @@ -10,20 +10,20 @@ 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); + ln1.logmean, ln1.logstd, + ln1_ci.low, ln1_ci.high); lognormal_params ln1_ci_paramas = 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.logmean, ln1.logstd); + ln1_ci.low, ln1_ci.high, + ln1.logmean, ln1.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 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); diff --git a/examples/13_ergonomic_algebra/example b/examples/13_ergonomic_algebra/example Binary files differ. diff --git a/examples/13_ergonomic_algebra/example.c b/examples/13_ergonomic_algebra/example.c @@ -6,9 +6,9 @@ #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) +#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() { @@ -16,8 +16,8 @@ int main() 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 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); diff --git a/examples/14_twitter_thread_example/example b/examples/14_twitter_thread_example/example Binary files differ. diff --git a/examples/14_twitter_thread_example/example.c b/examples/14_twitter_thread_example/example.c @@ -4,19 +4,23 @@ #include <stdio.h> #include <stdlib.h> -double sample_0(uint64_t* seed){ +double sample_0(uint64_t* seed) +{ return 0; } -double sample_1(uint64_t* seed){ +double sample_1(uint64_t* seed) +{ return 1; } -double sample_normal_mean_1_std_2(uint64_t* seed){ +double sample_normal_mean_1_std_2(uint64_t* seed) +{ return sample_normal(1, 2, seed); } -double sample_1_to_3(uint64_t* seed){ +double sample_1_to_3(uint64_t* seed) +{ return sample_to(1, 3, seed); } @@ -28,11 +32,11 @@ int main() 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 + double (*samplers[])(uint64_t*) = { + sample_0, + sample_1, + sample_normal_mean_1_std_2, + sample_1_to_3 }; int n_samples = 10; diff --git a/examples/15_plotting-scratchpad/makefile b/examples/15_plotting-scratchpad/makefile @@ -0,0 +1,3 @@ +build: + +format: diff --git a/examples/16_100_lognormal_samples/example.c b/examples/16_100_lognormal_samples/example.c @@ -1,17 +1,18 @@ #include "../../squiggle.h" #include <stdint.h> -#include <stdlib.h> #include <stdio.h> +#include <stdlib.h> // Estimate functions -int main(){ +int main() +{ // set randomness seed - uint64_t* seed = malloc(sizeof(uint64_t)); - *seed = 1000; // xorshift can't start with 0 - - for(int i=0; i<100; i++){ - double sample = sample_lognormal(0,10, seed); - printf("%f\n", sample); - } - free(seed); + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 + + for (int i = 0; i < 100; i++) { + double sample = sample_lognormal(0, 10, seed); + printf("%f\n", sample); + } + free(seed); } diff --git a/extra.c b/extra.c @@ -1,3 +1,4 @@ +#include "squiggle.h" #include <float.h> #include <limits.h> #include <math.h> @@ -6,7 +7,6 @@ #include <stdlib.h> #include <sys/types.h> #include <time.h> -#include "squiggle.h" // math constants #define PI 3.14159265358979323846 // M_PI in gcc gnu99 @@ -58,7 +58,6 @@ ci get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed) return result; } - // ## Sample from an arbitrary cdf struct box { int empty;