diff --git a/scratchpad/makefile b/scratchpad/makefile new file mode 100644 index 0000000..a2aed9e --- /dev/null +++ b/scratchpad/makefile @@ -0,0 +1,56 @@ +# Interface: +# make +# make build +# make format +# make run + +# Compiler +CC=gcc +# CC=tcc # <= faster compilation + +# Main file +SRC=scratchpad.c ../squiggle.c +OUTPUT=scratchpad + +## 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) + +verify: $(SRC) $(OUTPUT) + ./$(OUTPUT) | grep "NOT passed" -A 2 --group-separator='' || true + +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/scratchpad b/scratchpad/scratchpad new file mode 100755 index 0000000..2c89a47 Binary files /dev/null and b/scratchpad/scratchpad differ diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c new file mode 100644 index 0000000..2d49d32 --- /dev/null +++ b/scratchpad/scratchpad.c @@ -0,0 +1,20 @@ +#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 a seed of 0 + + for (int i = 0; i < 100; i++) { + double draw = sample_unit_uniform(seed); + printf("%f\n", draw); + + } + + free(seed); +} diff --git a/squiggle.c b/squiggle.c index b332b41..697aeee 100644 --- a/squiggle.c +++ b/squiggle.c @@ -13,7 +13,7 @@ #define PROCESS_ERROR(error_msg) process_error(error_msg, EXIT_ON_ERROR, __FILE__, __LINE__) #define PI 3.14159265358979323846 // M_PI in gcc gnu99 -#define NORMAL95CONFIDENCE 1.6448536269514722 +#define NORMAL90CONFIDENCE 1.6448536269514722 // # Key functionality // Define the minimum number of functions needed to do simple estimation @@ -23,9 +23,12 @@ uint64_t xorshift32(uint32_t* seed) { // Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" - // See - // https://en.wikipedia.org/wiki/Xorshift - // Also some drama: , + // See: + // + // , + // Also some drama: + // , + // uint64_t x = *seed; x ^= x << 13; x ^= x >> 17; @@ -53,8 +56,8 @@ double sample_unit_uniform(uint64_t* seed) double sample_unit_normal(uint64_t* seed) { - // See: - double u1 = sample_unit_uniform(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; @@ -96,7 +99,7 @@ inline double sample_normal_from_90_confidence_interval(double low, double high, // we can set mean = (high + low)/2; the midpoint, and L = high-low, // Normal([high + low]/2, [high - low]/(2 * 1.6448536269514722)) double mean = (high + low) / 2.0; - double std = (high - low) / (2.0 * NORMAL95CONFIDENCE); + double std = (high - low) / (2.0 * NORMAL90CONFIDENCE); return sample_normal(mean, std, seed); } @@ -507,14 +510,14 @@ lognormal_params convert_ci_to_lognormal_params(ci x) double loghigh = logf(x.high); double loglow = logf(x.low); double logmean = (loghigh + loglow) / 2.0; - double logstd = (loghigh - loglow) / (2.0 * NORMAL95CONFIDENCE); + double logstd = (loghigh - loglow) / (2.0 * NORMAL90CONFIDENCE); lognormal_params result = { .logmean = logmean, .logstd = logstd}; return result; } ci convert_lognormal_params_to_ci(lognormal_params y) { - double h = y.logstd * NORMAL95CONFIDENCE; + double h = y.logstd * NORMAL90CONFIDENCE; double loghigh = y.logmean + h; double loglow = y.logmean - h; ci result = { .low=exp(loglow), .high=exp(loghigh)};