diff --git a/README.md b/README.md index 607f1f1..d541ca7 100644 --- a/README.md +++ b/README.md @@ -348,6 +348,9 @@ It emits one warning about something I already took care of, so by default I've - [x] Provide example algebra - [x] Add conversion between 90% ci and parameters. - [ ] Use that conversion in conjuction with small algebra. + - [ ] Consider ergonomics of using ci instead of c_i + - use named struct instead + - demonstrate and document feeding a struct directly to a function; my_function((struct c_i){.low = 1, .high = 2}); - [ ] Test results - [x] Move to own file? Or signpost in file? => signposted in file. - [ ] Disambiguate sample_laplace--successes vs failures || successes vs total trials as two distinct and differently named functions diff --git a/examples/08_nuclear_war/example b/examples/08_nuclear_war/example index 8df8418..b8009cf 100755 Binary files a/examples/08_nuclear_war/example and b/examples/08_nuclear_war/example differ diff --git a/examples/08_nuclear_war/example.c b/examples/08_nuclear_war/example.c index 066dc3e..a17b44f 100644 --- a/examples/08_nuclear_war/example.c +++ b/examples/08_nuclear_war/example.c @@ -60,9 +60,9 @@ int main() } printf("... ]\n"); - struct c_i c_i_90 = get_90_confidence_interval(mixture, seed); + 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", c_i_90.low, c_i_90.high); + 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 index ebae9c9..b8b1ff6 100755 Binary files a/examples/10_nuclear_recovery/example and b/examples/10_nuclear_recovery/example differ diff --git a/examples/10_nuclear_recovery/example.c b/examples/10_nuclear_recovery/example.c index 9b25e97..3891cd7 100644 --- a/examples/10_nuclear_recovery/example.c +++ b/examples/10_nuclear_recovery/example.c @@ -46,8 +46,8 @@ int main() // Before a first nuclear collapse printf("## Before the first nuclear collapse\n"); - struct c_i c_i_90_2023 = get_90_confidence_interval(yearly_probability_nuclear_collapse_2023, seed); - printf("90%% confidence interval: [%f, %f]\n", c_i_90_2023.low, c_i_90_2023.high); + 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++) { @@ -57,8 +57,8 @@ int main() // After the first nuclear collapse printf("\n## After the first nuclear collapse\n"); - struct c_i c_i_90_2070 = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_example, seed); - printf("90%% confidence interval: [%f, %f]\n", c_i_90_2070.low, c_i_90_2070.high); + 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++) { @@ -68,8 +68,8 @@ int main() // After the first nuclear collapse (antiinductive) printf("\n## After the first nuclear collapse (antiinductive)\n"); - struct c_i c_i_90_antiinductive = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_antiinductive, seed); - printf("90%% confidence interval: [%f, %f]\n", c_i_90_antiinductive.low, c_i_90_antiinductive.high); + 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++) { diff --git a/examples/12_algebra_and_conversion/example b/examples/12_algebra_and_conversion/example new file mode 100755 index 0000000..7fc289f Binary files /dev/null and b/examples/12_algebra_and_conversion/example differ diff --git a/examples/12_algebra_and_conversion/example.c b/examples/12_algebra_and_conversion/example.c new file mode 100644 index 0000000..a93860b --- /dev/null +++ b/examples/12_algebra_and_conversion/example.c @@ -0,0 +1,33 @@ +#include "../../squiggle.h" +#include +#include +#include +#include + +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_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); + + 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/12_algebra_and_conversion/makefile b/examples/12_algebra_and_conversion/makefile new file mode 100644 index 0000000..ef385f7 --- /dev/null +++ b/examples/12_algebra_and_conversion/makefile @@ -0,0 +1,53 @@ +# Interface: +# make +# make build +# make format +# make run + +# Compiler +CC=gcc +# CC=tcc # <= faster compilation + +# Main file +SRC=example.c ../../squiggle.c +OUTPUT=example + +## Dependencies +MATH=-lm + +## Flags +DEBUG= #'-g' +STANDARD=-std=c99 +WARNINGS=-Wall +OPTIMIZED=-O3 #-Ofast +# OPENMP=-fopenmp + +## Formatter +STYLE_BLUEPRINT=webkit +FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT) + +## make build +build: $(SRC) + $(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(MATH) -o $(OUTPUT) + +format: $(SRC) + $(FORMATTER) $(SRC) + +run: $(SRC) $(OUTPUT) + OMP_NUM_THREADS=1 ./$(OUTPUT) && echo + +time-linux: + @echo "Requires /bin/time, found on GNU/Linux systems" && echo + + @echo "Running 100x and taking avg time $(OUTPUT)" + @t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 1 thread: |" | sed 's|$$|ms|' && echo + +## Profiling + +profile-linux: + echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar" + echo "Must be run as sudo" + $(CC) $(SRC) $(MATH) -o $(OUTPUT) + sudo perf record ./$(OUTPUT) + sudo perf report + rm perf.data diff --git a/squiggle.c b/squiggle.c index 9151ccd..65e808d 100644 --- a/squiggle.c +++ b/squiggle.c @@ -435,10 +435,10 @@ double sampler_danger(struct box cdf(double), uint64_t* seed) // Get confidence intervals, given a sampler -struct c_i { +typedef struct ci_t { float low; float high; -}; +} ci; int compare_doubles(const void* p, const void* q) { // https://wikiless.esmailelbob.xyz/wiki/Qsort?lang=en @@ -454,7 +454,7 @@ int compare_doubles(const void* p, const void* q) return 0; } -struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed) +ci get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed) { int n = 100 * 1000; double* samples_array = malloc(n * sizeof(double)); @@ -463,7 +463,7 @@ struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* se } qsort(samples_array, n, sizeof(double), compare_doubles); - struct c_i result = { + ci result = { .low = samples_array[5000], .high = samples_array[94999], }; @@ -505,7 +505,7 @@ lognormal_params algebra_product_lognormals(lognormal_params a, lognormal_params return result; } -lognormal_params convert_ci_to_lognormal_params(struct c_i x) +lognormal_params convert_ci_to_lognormal_params(ci x) { double loghigh = logf(x.high); double loglow = logf(x.low); @@ -515,12 +515,12 @@ lognormal_params convert_ci_to_lognormal_params(struct c_i x) return result; } -struct c_i convert_lognormal_params_to_ci(lognormal_params y) +ci convert_lognormal_params_to_ci(lognormal_params y) { double h = y.logstd * NORMAL95CONFIDENCE; double loghigh = y.logmean + h; double loglow = y.logmean - h; - struct c_i result = { .low=exp(loglow), .high=exp(loghigh)}; + ci result = { .low=exp(loglow), .high=exp(loghigh)}; return result; } diff --git a/squiggle.h b/squiggle.h index 8a93325..b97664c 100644 --- a/squiggle.h +++ b/squiggle.h @@ -52,11 +52,11 @@ struct box sampler_cdf_double(double cdf(double), uint64_t* seed); struct box sampler_cdf_box(struct box cdf(double), uint64_t* seed); // Get 90% confidence interval -struct c_i { +typedef struct ci_t { float low; float high; -}; -struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed); +} ci; +ci get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed); // small algebra manipulations @@ -73,9 +73,9 @@ typedef struct lognormal_params_t { lognormal_params algebra_product_lognormals(lognormal_params a, lognormal_params b); -lognormal_params convert_ci_to_lognormal_params(struct c_i x); +lognormal_params convert_ci_to_lognormal_params(ci x); -struct c_i convert_lognormal_params_to_ci(lognormal_params y); +ci convert_lognormal_params_to_ci(lognormal_params y); #endif