diff --git a/examples/04_sample_from_cdf/example b/examples/04_sample_from_cdf/example deleted file mode 100755 index 03e23d9..0000000 Binary files a/examples/04_sample_from_cdf/example and /dev/null differ diff --git a/examples/04_sample_from_cdf_simple/example b/examples/04_sample_from_cdf_simple/example new file mode 100755 index 0000000..37f99c9 Binary files /dev/null and b/examples/04_sample_from_cdf_simple/example differ diff --git a/examples/04_sample_from_cdf_simple/example.c b/examples/04_sample_from_cdf_simple/example.c new file mode 100644 index 0000000..ab6ac61 --- /dev/null +++ b/examples/04_sample_from_cdf_simple/example.c @@ -0,0 +1,102 @@ +#include +#include +#include +#include +#include +#include "../../squiggle.h" + +#define NUM_SAMPLES 1000000 + +// Example cdf +float cdf_uniform_0_1(float x) +{ + if (x < 0) { + return 0; + } else if (x > 1) { + return 1; + } else { + return x; + } +} + +float cdf_squared_0_1(float x) +{ + if (x < 0) { + return 0; + } else if (x > 1) { + return 1; + } else { + return x * x; + } +} + +float cdf_normal_0_1(float x) +{ + float mean = 0; + float std = 1; + return 0.5 * (1 + erf((x - mean) / (std * sqrt(2)))); // erf from math.h +} + +// Some testers +void test_inverse_cdf_float(char* cdf_name, float cdf_float(float)) +{ + struct box result = inverse_cdf_float(cdf_float, 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_float(char* cdf_name, float cdf_float(float), uint32_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_float_cdf(cdf_float, seed); + if (sample.empty) { + printf("Error in sampler function for %s", cdf_name); + } else { + // printf("%f\n", sample.content); + } + } + clock_t end = clock(); + float time_spent = (float)(end - begin) / CLOCKS_PER_SEC; + printf("Time spent: %f\n", time_spent); +} + +int main() +{ + // Test inverse cdf float + test_inverse_cdf_float("cdf_uniform_0_1", cdf_uniform_0_1); + test_inverse_cdf_float("cdf_squared_0_1", cdf_squared_0_1); + test_inverse_cdf_float("cdf_normal_0_1", cdf_normal_0_1); + + // Testing samplers + // set randomness seed + uint32_t* seed = malloc(sizeof(uint32_t)); + *seed = 1000; // xorshift can't start with 0 + + // Test float sampler + test_and_time_sampler_float("cdf_uniform_0_1", cdf_uniform_0_1, seed); + test_and_time_sampler_float("cdf_squared_0_1", cdf_squared_0_1, seed); + test_and_time_sampler_float("cdf_normal_0_1", cdf_normal_0_1, seed); + + // Get some normal samples using a previous approach + printf("\nGetting some samples from unit_normal\n"); + + clock_t begin_2 = clock(); + + for (int i = 0; i < NUM_SAMPLES; i++) { + float normal_sample = unit_normal(seed); + // printf("%f\n", normal_sample); + } + + clock_t end_2 = clock(); + float time_spent_2 = (float)(end_2 - begin_2) / CLOCKS_PER_SEC; + printf("Time spent: %f\n", time_spent_2); + + free(seed); + return 0; +} diff --git a/examples/04_sample_from_cdf/makefile b/examples/04_sample_from_cdf_simple/makefile similarity index 100% rename from examples/04_sample_from_cdf/makefile rename to examples/04_sample_from_cdf_simple/makefile diff --git a/examples/05_sample_from_cdf_beta/example b/examples/05_sample_from_cdf_beta/example new file mode 100755 index 0000000..1ad5fec Binary files /dev/null and b/examples/05_sample_from_cdf_beta/example differ diff --git a/examples/04_sample_from_cdf/example.c b/examples/05_sample_from_cdf_beta/example.c similarity index 68% rename from examples/04_sample_from_cdf/example.c rename to examples/05_sample_from_cdf_beta/example.c index c99986b..87797a6 100644 --- a/examples/04_sample_from_cdf/example.c +++ b/examples/05_sample_from_cdf_beta/example.c @@ -1,4 +1,4 @@ -#include // erf, sqrt +#include #include #include #include @@ -9,36 +9,7 @@ #define STOP_BETA 1.0e-8 #define TINY_BETA 1.0e-30 -// Example cdf -float cdf_uniform_0_1(float x) -{ - if (x < 0) { - return 0; - } else if (x > 1) { - return 1; - } else { - return x; - } -} - -float cdf_squared_0_1(float x) -{ - if (x < 0) { - return 0; - } else if (x > 1) { - return 1; - } else { - return x * x; - } -} - -float cdf_normal_0_1(float x) -{ - float mean = 0; - float std = 1; - return 0.5 * (1 + erf((x - mean) / (std * sqrt(2)))); // erf from math.h -} - +// Incomplete beta function struct box incbeta(float a, float b, float x) { // Descended from , @@ -149,16 +120,6 @@ struct box cdf_beta(float x) } // Some testers -void test_inverse_cdf_float(char* cdf_name, float cdf_float(float)) -{ - struct box result = inverse_cdf_float(cdf_float, 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_inverse_cdf_box(char* cdf_name, struct box cdf_box(float)) { struct box result = inverse_cdf_box(cdf_box, 0.5); @@ -170,23 +131,6 @@ void test_inverse_cdf_box(char* cdf_name, struct box cdf_box(float)) } } -void test_and_time_sampler_float(char* cdf_name, float cdf_float(float), uint32_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_float_cdf(cdf_float, seed); - if (sample.empty) { - printf("Error in sampler function for %s", cdf_name); - } else { - // printf("%f\n", sample.content); - } - } - clock_t end = clock(); - float time_spent = (float)(end - begin) / CLOCKS_PER_SEC; - printf("Time spent: %f\n", time_spent); -} - void test_and_time_sampler_box(char* cdf_name, struct box cdf_box(float), uint32_t* seed) { printf("\nGetting some samples from %s:\n", cdf_name); @@ -206,44 +150,17 @@ void test_and_time_sampler_box(char* cdf_name, struct box cdf_box(float), uint32 int main() { - // Test inverse cdf float - test_inverse_cdf_float("cdf_uniform_0_1", cdf_uniform_0_1); - test_inverse_cdf_float("cdf_squared_0_1", cdf_squared_0_1); - test_inverse_cdf_float("cdf_normal_0_1", cdf_normal_0_1); - // Test inverse cdf box test_inverse_cdf_box("cdf_beta", cdf_beta); - // Testing samplers - // set randomness seed + // Test box sampler uint32_t* seed = malloc(sizeof(uint32_t)); *seed = 1000; // xorshift can't start with 0 - - // Test float sampler - test_and_time_sampler_float("cdf_uniform_0_1", cdf_uniform_0_1, seed); - test_and_time_sampler_float("cdf_squared_0_1", cdf_squared_0_1, seed); - test_and_time_sampler_float("cdf_normal_0_1", cdf_normal_0_1, seed); - - // Get some normal samples using a previous approach - printf("\nGetting some samples from unit_normal\n"); - - clock_t begin_2 = clock(); - - for (int i = 0; i < NUM_SAMPLES; i++) { - float normal_sample = unit_normal(seed); - // printf("%f\n", normal_sample); - } - - clock_t end_2 = clock(); - float time_spent_2 = (float)(end_2 - begin_2) / CLOCKS_PER_SEC; - printf("Time spent: %f\n", time_spent_2); - - // Test box sampler 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 really anal convergence conditions. + // But I am also using absurdly precise convergence conditions. // This could be optimized. free(seed); diff --git a/examples/05_sample_from_cdf_beta/makefile b/examples/05_sample_from_cdf_beta/makefile new file mode 100644 index 0000000..69adb6d --- /dev/null +++ b/examples/05_sample_from_cdf_beta/makefile @@ -0,0 +1,57 @@ +# Interface: +# make +# make build +# make format +# make run + +# Compiler +CC=gcc # required for nested functions +# CC=tcc # <= faster compilation + +# Main file +SRC=example.c ../../squiggle.c +OUTPUT=./example + +## Dependencies +MATH=-lm +DEPENDENCIES=$(MATH) +# OPENMP=-fopenmp + +## Flags +DEBUG= #'-g' +STANDARD=-std=gnu99 ## allows for nested functions. +EXTENSIONS= #-fnested-functions +WARNINGS=-Wall +OPTIMIZED=-O3#-Ofast +CFLAGS=$(DEBUG) $(STANDARD) $(EXTENSIONS) $(WARNINGS) $(OPTIMIZED) + +## 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) + +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