diff --git a/README.md b/README.md index 91fdb5c..7c2a6b8 100644 --- a/README.md +++ b/README.md @@ -2,6 +2,8 @@ squiggle.c is a self-contained C99 library that provides functions for simple Monte Carlo estimation, based on [Squiggle](https://www.squiggle-language.com/). +![](./core.png) + ## Why C? - Because it is fast diff --git a/core.png b/core.png new file mode 100644 index 0000000..3a35db1 Binary files /dev/null and b/core.png differ diff --git a/examples/02_many_samples_time_to_botec/example.c b/examples/02_many_samples_time_to_botec/example.c index 9bb21e7..2ddac68 100644 --- a/examples/02_many_samples_time_to_botec/example.c +++ b/examples/02_many_samples_time_to_botec/example.c @@ -1,7 +1,7 @@ -#include -#include -#include #include "../../squiggle.h" +#include +#include +#include // 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; @@ -37,18 +38,18 @@ int main(){ double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 }; double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many }; - int n_samples = 1000000; - double* result_many = (double *) malloc(n_samples * sizeof(double)); - for(int i=0; i +#include +#include + +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/14_twitter_thread_example/makefile b/examples/14_twitter_thread_example/makefile new file mode 100644 index 0000000..2d251e0 --- /dev/null +++ b/examples/14_twitter_thread_example/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) + ./$(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/scratchpad/core.c b/scratchpad/core.c new file mode 100644 index 0000000..20b59e5 --- /dev/null +++ b/scratchpad/core.c @@ -0,0 +1,27 @@ + +uint64_t xorshift64(uint64_t* seed) +{ + // Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" + // + uint64_t x = *seed; + x ^= x << 13; + x ^= x >> 7; + x ^= x << 17; + return *seed = x; +} + +double sample_unit_uniform(uint64_t* seed) +{ + // samples uniform from [0,1] interval. + return ((double)xorshift64(seed)) / ((double)UINT64_MAX); +} + +double sample_unit_normal(uint64_t* seed) +{ + // // See: + double u1 = sample_unit_uniform(seed); + double u2 = sample_unit_uniform(seed); + double z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2); + return z; +} + diff --git a/squiggle.c b/squiggle.c index 697aeee..6222d95 100644 --- a/squiggle.c +++ b/squiggle.c @@ -26,8 +26,8 @@ uint64_t xorshift32(uint32_t* seed) // See: // // , - // Also some drama: - // , + // Also some drama: + // , // uint64_t x = *seed; x ^= x << 13; @@ -57,7 +57,7 @@ double sample_unit_uniform(uint64_t* seed) double sample_unit_normal(uint64_t* seed) { // // See: - // double u1 = sample_unit_uniform(seed); + double u1 = sample_unit_uniform(seed); double u2 = sample_unit_uniform(seed); double z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2); return z; @@ -109,7 +109,7 @@ double sample_to(double low, double high, uint64_t* seed) // returns a sample from a lognorma with a matching 90% c.i. // Key idea: If we want a lognormal with 90% confidence interval [a, b] // we need but get a normal with 90% confidence interval [log(a), log(b)]. - // Then see code for sample_normal_from_95_confidence_interval + // Then see code for sample_normal_from_90_confidence_interval double loglow = logf(low); double loghigh = logf(high); return exp(sample_normal_from_90_confidence_interval(loglow, loghigh, seed)); @@ -511,7 +511,7 @@ lognormal_params convert_ci_to_lognormal_params(ci x) double loglow = logf(x.low); double logmean = (loghigh + loglow) / 2.0; double logstd = (loghigh - loglow) / (2.0 * NORMAL90CONFIDENCE); - lognormal_params result = { .logmean = logmean, .logstd = logstd}; + lognormal_params result = { .logmean = logmean, .logstd = logstd }; return result; } @@ -520,8 +520,6 @@ ci convert_lognormal_params_to_ci(lognormal_params y) double h = y.logstd * NORMAL90CONFIDENCE; double loghigh = y.logmean + h; double loglow = y.logmean - h; - ci result = { .low=exp(loglow), .high=exp(loghigh)}; + ci result = { .low = exp(loglow), .high = exp(loghigh) }; return result; - } -