diff --git a/README.md b/README.md index eba70a8..51f4bc0 100644 --- a/README.md +++ b/README.md @@ -195,10 +195,6 @@ The bailey: - I think the core interface is not likely to change much, though I've recently changed the interface for parallelism and for getting confidence intervals. - I am using this code for a few important consulting projects, and I trust myself to operate it correctly. -## Licensing - -This project is released under the MIT license, a permissive open-source license. You can see it in the LICENSE.txt file. - ## squiggle extra functions ### Division between core functions and squiggle_more expansions @@ -250,32 +246,7 @@ That said, I found adding parallelism to be an interesting an engaging task. Mos `squiggle_more.c` has some functions to do some simple algebra manipulations: sums of normals and products of lognormals. You can see some example usage [here](examples/more/07_algebra/example.c) and [here](examples/more/08_algebra_and_conversion/example.c). -#### Extra: cdf auxiliary functions +## Licensing -I provide some auxiliary functions that take a cdf, and return a sample from the distribution produced by that cdf. This might make it easier to program models, at the cost of a 20x to 60x slowdown in the parts of the code that use it. - -I don't have any immediate plans or reason to delete these functions, but I probably will, since they are too slow/inefficient for my taste. - -#### Extra: Error propagation vs exiting on error - -The process of taking a cdf and returning a sample might fail, e.g., it's a Newton method which might fail to converge because of cdf artifacts. The cdf itself might also fail, e.g., if a distribution only accepts a range of parameters, but is fed parameters outside that range. - -This library provides two approaches: - -1. Print the line and function in which the error occurred, then exit on error -2. In situations where there might be an error, return a struct containing either the correct value or an error message: - -```C -struct box { - int empty; - double content; - char* error_msg; -}; -``` - -The first approach produces terser programs but might not scale. The second approach seems like it could lead to more robust programs, but is more verbose. - -Behaviour on error can be toggled by the `EXIT_ON_ERROR` variable. This library also provides a convenient macro, `PROCESS_ERROR`, to make error handling in either case much terser—see the usage in example 4 in the examples/ folder. - -Overall, I'd describe the error handling capabilities of this library as pretty rudimentary. For example, this program might fail in surprising ways if you ask for a lognormal with negative standard deviation, because I haven't added error checking for that case yet. +This project is released under the MIT license, a permissive open-source license. You can see it in the LICENSE.txt file. diff --git a/ROADMAP.md b/ROADMAP.md index fd9d8c3..cfa3d3b 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -2,15 +2,17 @@ ## To do -- [ ] Big refactor +- [x] Big refactor - [ ] Come up with a better headline example; fermi paradox paper is too complicated - - [ ] Make README.md less messy + - [x] Make README.md less messy - [ ] Give examples of new functions + - [ ] Reference commit with cdf functions, even though deleted - [ ] Post on suckless subreddit - [ ] Drive in a few more real-life applications - [ ] US election modelling? - [ ] Look into using size_t instead of int for sample numbers - [ ] Reorganize code a little bit to reduce usage of gcc's nested functions +- [ ] Rename examples ## Done diff --git a/examples/more/00_example_template/example b/examples/more/00_example_template/example index d7d88bf..c7304b6 100755 Binary files a/examples/more/00_example_template/example and b/examples/more/00_example_template/example differ diff --git a/examples/more/01_sample_from_cdf/example b/examples/more/01_sample_from_cdf/example deleted file mode 100755 index 852790a..0000000 Binary files a/examples/more/01_sample_from_cdf/example and /dev/null differ diff --git a/examples/more/01_sample_from_cdf/example.c b/examples/more/01_sample_from_cdf/example.c deleted file mode 100644 index f56dd4e..0000000 --- a/examples/more/01_sample_from_cdf/example.c +++ /dev/null @@ -1,102 +0,0 @@ -#include "../../../squiggle.h" -#include "../../../squiggle_more.h" -#include -#include -#include -#include - -#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)) -{ - 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++) { - 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 deleted file mode 100755 index 4156892..0000000 Binary files a/examples/more/02_sample_from_cdf_beta/example and /dev/null differ diff --git a/examples/more/02_sample_from_cdf_beta/example.c b/examples/more/02_sample_from_cdf_beta/example.c deleted file mode 100644 index 1449e8d..0000000 --- a/examples/more/02_sample_from_cdf_beta/example.c +++ /dev/null @@ -1,168 +0,0 @@ -#include "../../../squiggle.h" -#include "../../../squiggle_more.h" -#include -#include -#include -#include - -#define NUM_SAMPLES 10000 -#define STOP_BETA 1.0e-8 -#define TINY_BETA 1.0e-30 - -// Incomplete beta function -box incbeta(double a, double b, double x) -{ - // Descended from , - // - // 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)) { - box symmetric_incbeta = incbeta(b, a, 1.0 - x); - if (symmetric_incbeta.empty) { - return symmetric_incbeta; // propagate error - } else { - 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) { - box result = { - .empty = 0, - .content = front * (f - 1.0) - }; - return result; - } - } - - return PROCESS_ERROR("More loops needed, did not converge, in function incbeta"); -} - -box cdf_beta(double x) -{ - if (x < 0) { - box result = { .empty = 0, .content = 0 }; - return result; - } else if (x > 1) { - 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, box cdf_box(double)) -{ - 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, 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++) { - 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/03_ci_beta/example b/examples/more/03_ci_beta/example index 497dcc9..7022d14 100755 Binary files a/examples/more/03_ci_beta/example and b/examples/more/03_ci_beta/example differ diff --git a/examples/more/03_ci_beta/example.c b/examples/more/03_ci_beta/example.c index 78019f2..fd9581c 100644 --- a/examples/more/03_ci_beta/example.c +++ b/examples/more/03_ci_beta/example.c @@ -4,9 +4,9 @@ #include // Estimate functions -double beta_1_2_sampler(uint64_t* seed) +double sample_beta_3_2(uint64_t* seed) { - return sample_beta(1, 2.0, seed); + return sample_beta(3.0, 2.0, seed); } int main() @@ -15,9 +15,16 @@ int main() uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with 0 - ci beta_1_2_ci_90 = sampler_get_90_ci(beta_1_2_sampler, 1000000, 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); - printf("You can check this in \n"); + int n_samples = 1 * MILLION; + double* xs = malloc(sizeof(double) * (size_t)n_samples); + for (int i = 0; i < n_samples; i++) { + xs[i] = sample_beta_3_2(seed); + } + + printf("\n# Stats\n"); + array_print_stats(xs, n_samples); + printf("\n# Histogram\n"); + array_print_histogram(xs, n_samples, 23); free(seed); } diff --git a/examples/more/04_nuclear_war/example b/examples/more/04_nuclear_war/example index 3c552e8..a1994d0 100755 Binary files a/examples/more/04_nuclear_war/example and b/examples/more/04_nuclear_war/example differ diff --git a/examples/more/04_nuclear_war/example.c b/examples/more/04_nuclear_war/example.c index 9cf0678..065815c 100644 --- a/examples/more/04_nuclear_war/example.c +++ b/examples/more/04_nuclear_war/example.c @@ -34,7 +34,7 @@ double probability_of_dying_eli(uint64_t* seed) return probability_of_dying; } -double mixture(uint64_t* seed) +double sample_nuclear_model(uint64_t* seed) { double (*samplers[])(uint64_t*) = { probability_of_dying_nuno, probability_of_dying_eli }; double weights[] = { 0.5, 0.5 }; @@ -47,23 +47,17 @@ int main() 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) * (size_t)n); + int n = 1 * MILLION; + double* xs = malloc(sizeof(double) * (size_t)n); for (int i = 0; i < n; i++) { - mixture_result[i] = mixture(seed); + xs[i] = sample_nuclear_model(seed); } - printf("mixture_result: [ "); - for (int i = 0; i < 9; i++) { - printf("%.6f, ", mixture_result[i]); - } - printf("... ]\n"); + printf("\n# Stats\n"); + array_print_stats(xs, n); + printf("\n# Histogram\n"); + array_print_90_ci_histogram(xs, n, 20); - ci ci_90 = sampler_get_90_ci(mixture, 1000000, seed); - printf("mean: %f\n", array_mean(mixture_result, n)); - printf("90%% confidence interval: [%f, %f]\n", ci_90.low, ci_90.high); - - free(mixture_result); + free(xs); free(seed); } diff --git a/examples/more/05_burn_10kg_fat/example b/examples/more/05_burn_10kg_fat/example index 02fb690..4a7c6cb 100755 Binary files a/examples/more/05_burn_10kg_fat/example and b/examples/more/05_burn_10kg_fat/example differ diff --git a/examples/more/05_burn_10kg_fat/example.c b/examples/more/05_burn_10kg_fat/example.c index a2b7d82..b690ec5 100644 --- a/examples/more/05_burn_10kg_fat/example.c +++ b/examples/more/05_burn_10kg_fat/example.c @@ -27,22 +27,17 @@ int main() *seed = 1000; // xorshift can't start with 0 int n = 1000 * 1000; - - double* result = malloc(sizeof(double) * (size_t)n); + double* xs = malloc(sizeof(double) * (size_t)n); for (int i = 0; i < n; i++) { - result[i] = sample_minutes_per_day_jumping_rope_needed_to_burn_10kg(seed); + xs[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 = sampler_get_90_ci(sample_minutes_per_day_jumping_rope_needed_to_burn_10kg, 1000000, seed); - printf("90%% confidence interval: [%f, %f]\n", ci_90.low, ci_90.high); + printf("\n# Stats\n"); + array_print_stats(xs, n); + printf("\n# Histogram\n"); + array_print_histogram(xs, n, 23); free(seed); } diff --git a/examples/more/06_nuclear_recovery/example b/examples/more/06_nuclear_recovery/example index 682fb84..ede6020 100755 Binary files a/examples/more/06_nuclear_recovery/example and b/examples/more/06_nuclear_recovery/example differ diff --git a/examples/more/06_nuclear_recovery/example.c b/examples/more/06_nuclear_recovery/example.c index f3cbeaf..9d3e3ae 100644 --- a/examples/more/06_nuclear_recovery/example.c +++ b/examples/more/06_nuclear_recovery/example.c @@ -46,41 +46,37 @@ int main() uint64_t* seed = malloc(sizeof(uint64_t)); *seed = 1000; // xorshift can't start with 0 - int num_samples = 1000000; + int n_samples = 1000000; // Before a first nuclear collapse printf("## Before the first nuclear collapse\n"); - ci ci_90_2023 = sampler_get_90_ci(yearly_probability_nuclear_collapse_2023, 1000000, 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) * (size_t)num_samples); - for (int i = 0; i < num_samples; i++) { + double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * (size_t)n_samples); + for (int i = 0; i < n_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)); + ci ci_90_2023 = array_get_90_ci(yearly_probability_nuclear_collapse_2023_samples, n_samples); + printf("90%% confidence interval: [%f, %f]\n", ci_90_2023.low, ci_90_2023.high); // After the first nuclear collapse printf("\n## After the first nuclear collapse\n"); - ci ci_90_2070 = sampler_get_90_ci(yearly_probability_nuclear_collapse_after_recovery_example, 1000000, 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) * (size_t)num_samples); - for (int i = 0; i < num_samples; i++) { + double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * (size_t)n_samples); + for (int i = 0; i < n_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)); + ci ci_90_2070 = array_get_90_ci(yearly_probability_nuclear_collapse_after_recovery_samples, 1000000); + printf("90%% confidence interval: [%f, %f]\n", ci_90_2070.low, ci_90_2070.high); // After the first nuclear collapse (antiinductive) printf("\n## After the first nuclear collapse (antiinductive)\n"); - ci ci_90_antiinductive = sampler_get_90_ci(yearly_probability_nuclear_collapse_after_recovery_antiinductive, 1000000, 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) * (size_t)num_samples); - for (int i = 0; i < num_samples; i++) { + double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * (size_t)n_samples); + for (int i = 0; i < n_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)); + ci ci_90_antiinductive = array_get_90_ci(yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples, 1000000); + printf("90%% confidence interval: [%f, %f]\n", ci_90_antiinductive.low, ci_90_antiinductive.high); + // free seeds free(yearly_probability_nuclear_collapse_2023_samples); free(yearly_probability_nuclear_collapse_after_recovery_samples); free(yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples); diff --git a/examples/more/07_algebra/example b/examples/more/07_algebra/example index 2b171ce..826e264 100755 Binary files a/examples/more/07_algebra/example and b/examples/more/07_algebra/example differ diff --git a/examples/more/08_algebra_and_conversion/example b/examples/more/08_algebra_and_conversion/example index d00c728..dc0e757 100755 Binary files a/examples/more/08_algebra_and_conversion/example and b/examples/more/08_algebra_and_conversion/example differ diff --git a/examples/more/09_ergonomic_algebra/example b/examples/more/09_ergonomic_algebra/example index fca55dc..6109162 100755 Binary files a/examples/more/09_ergonomic_algebra/example and b/examples/more/09_ergonomic_algebra/example differ diff --git a/examples/more/10_twitter_thread_example/example b/examples/more/10_twitter_thread_example/example index 8fd40c2..6db41c4 100755 Binary files a/examples/more/10_twitter_thread_example/example and b/examples/more/10_twitter_thread_example/example differ diff --git a/examples/more/11_billion_lognormals_paralell/example b/examples/more/11_billion_lognormals_paralell/example index 8065041..6584d76 100755 Binary files a/examples/more/11_billion_lognormals_paralell/example and b/examples/more/11_billion_lognormals_paralell/example differ diff --git a/examples/more/12_time_to_botec_parallel/example b/examples/more/12_time_to_botec_parallel/example index 09b5a24..6faf293 100755 Binary files a/examples/more/12_time_to_botec_parallel/example and b/examples/more/12_time_to_botec_parallel/example differ diff --git a/examples/more/13_parallelize_min/example b/examples/more/13_parallelize_min/example index 155af9d..b63cdc0 100755 Binary files a/examples/more/13_parallelize_min/example and b/examples/more/13_parallelize_min/example differ diff --git a/examples/more/14_check_confidence_interval/example b/examples/more/14_check_confidence_interval/example index 3f880bd..0d584dc 100755 Binary files a/examples/more/14_check_confidence_interval/example and b/examples/more/14_check_confidence_interval/example differ diff --git a/examples/more/makefile b/examples/more/makefile index 14d7ac2..0bcbde1 100644 --- a/examples/more/makefile +++ b/examples/more/makefile @@ -39,8 +39,6 @@ FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT) ## make all all: $(CC) $(OPTIMIZED) $(DEBUG) $(WARN) 00_example_template/$(SRC) $(DEPS) -o 00_example_template/$(OUTPUT) - $(CC) $(OPTIMIZED) $(DEBUG) $(WARN) 01_sample_from_cdf/$(SRC) $(DEPS) -o 01_sample_from_cdf/$(OUTPUT) - $(CC) $(OPTIMIZED) $(DEBUG) $(WARN) 02_sample_from_cdf_beta/$(SRC) $(DEPS) -o 02_sample_from_cdf_beta/$(OUTPUT) $(CC) $(OPTIMIZED) $(DEBUG) $(WARN) 03_ci_beta/$(SRC) $(DEPS) -o 03_ci_beta/$(OUTPUT) $(CC) $(OPTIMIZED) $(DEBUG) $(WARN) 04_nuclear_war/$(SRC) $(DEPS) -o 04_nuclear_war/$(OUTPUT) $(CC) $(OPTIMIZED) $(DEBUG) $(WARN) 05_burn_10kg_fat/$(SRC) $(DEPS) -o 05_burn_10kg_fat/$(OUTPUT) diff --git a/squiggle_more.c b/squiggle_more.c index 06781af..1168845 100644 --- a/squiggle_more.c +++ b/squiggle_more.c @@ -212,15 +212,15 @@ void array_print_stats(double xs[], int n){ double median = array_get_median(xs, n); double mean = array_mean(xs, n); double std = array_std(xs, n); - printf("Mean: %lf\n" - " Std: %lf\n" - " 5%%: %lf\n" - " 10%%: %lf\n" - " 25%%: %lf\n" - " 50%%: %lf\n" - " 75%%: %lf\n" - " 90%%: %lf\n" - " 95%%: %lf\n", + printf("Avg: %lf\n" + "Std: %lf\n" + " 5%%: %lf\n" + "10%%: %lf\n" + "25%%: %lf\n" + "50%%: %lf\n" + "75%%: %lf\n" + "90%%: %lf\n" + "95%%: %lf\n", mean, std, ci_90.low, ci_80.low, ci_50.low, median, ci_50.high, ci_80.high, ci_90.high); } @@ -293,7 +293,7 @@ void array_print_histogram(double* xs, int n_samples, int n_bins) { decimalPlaces = -magnitude; decimalPlaces = decimalPlaces > 10 ? 10 : decimalPlaces; } - printf(" [%*.*f, %*.*f", 4+decimalPlaces, decimalPlaces, bin_start, 4+decimalPlaces, decimalPlaces, bin_end); + printf("[%*.*f, %*.*f", 4+decimalPlaces, decimalPlaces, bin_start, 4+decimalPlaces, decimalPlaces, bin_end); char interval_delimiter = ')'; if(i == (n_bins-1)){ interval_delimiter = ']'; // last bucket is inclusive @@ -377,12 +377,12 @@ void array_print_90_ci_histogram(double* xs, int n_samples, int n_bins){ decimalPlaces = -magnitude; decimalPlaces = decimalPlaces > 10 ? 10 : decimalPlaces; } - printf( " (-∞, %*.*f): %d\n", 4+decimalPlaces, decimalPlaces, min_value, below_min); + printf( "(-∞, %*.*f): %d\n", 4+decimalPlaces, decimalPlaces, min_value, below_min); for (int i = 0; i < n_bins; i++) { double bin_start = min_value + i * bin_width; double bin_end = bin_start + bin_width; - printf(" [%*.*f, %*.*f", 4+decimalPlaces, decimalPlaces, bin_start, 4+decimalPlaces, decimalPlaces, bin_end); + printf("[%*.*f, %*.*f", 4+decimalPlaces, decimalPlaces, bin_start, 4+decimalPlaces, decimalPlaces, bin_end); char interval_delimiter = ')'; if(i == (n_bins-1)){ interval_delimiter = ']'; // last bucket is inclusive @@ -395,7 +395,7 @@ void array_print_90_ci_histogram(double* xs, int n_samples, int n_bins){ } printf(" %d\n", bins[i]); } - printf( " (%*.*f, +∞): %d\n", 4+decimalPlaces, decimalPlaces, max_value, above_max); + printf( "(%*.*f, +∞): %d\n", 4+decimalPlaces, decimalPlaces, max_value, above_max); // Free the allocated memory for bins free(bins); diff --git a/squiggle_more.h b/squiggle_more.h index fa9ef43..54afcde 100644 --- a/squiggle_more.h +++ b/squiggle_more.h @@ -17,10 +17,6 @@ void array_print_stats(double xs[], int n); void array_print_histogram(double* xs, int n_samples, int n_bins); void array_print_90_ci_histogram(double* xs, int n, int n_bins); -// Deprecated: get confidence intervals directly from samplers -ci sampler_get_ci(ci interval, double (*sampler)(uint64_t*), int n, uint64_t* seed); -ci sampler_get_90_ci(double (*sampler)(uint64_t*), int n, uint64_t* seed); - /* Algebra manipulations */ typedef struct normal_params_t { @@ -38,24 +34,9 @@ lognormal_params algebra_product_lognormals(lognormal_params a, lognormal_params lognormal_params convert_ci_to_lognormal_params(ci x); ci convert_lognormal_params_to_ci(lognormal_params y); -/* Error handling */ -typedef struct box_t { - int empty; - double content; - char* error_msg; -} box; -#define MAX_ERROR_LENGTH 500 -#define EXIT_ON_ERROR 0 -#define PROCESS_ERROR(error_msg) process_error(error_msg, EXIT_ON_ERROR, __FILE__, __LINE__) -box process_error(const char* error_msg, int should_exit, char* file, int line); -void array_print(double* array, int length); +/* Utilities */ -/* Inverse cdf */ -box inverse_cdf_double(double cdf(double), double p); -box inverse_cdf_box(box cdf_box(double), double p); - -/* Samplers from cdf */ -box sampler_cdf_double(double cdf(double), uint64_t* seed); -box sampler_cdf_box(box cdf(double), uint64_t* seed); +#define THOUSAND 1000 +#define MILLION 1000000 #endif