squiggle.c

commit b1a58f9b74946c6ede7e47fd5035b1879d2578a1
parent 015d33adca87d305784402ce89b0b75821878a5b
Author: NunoSempere <nuno.sempere@protonmail.com>
Date:   Wed, 27 Sep 2023 14:10:40 +0100

fix constant name 95=>90

Diffstat:
A
scratchpad/makefile
 | 
56
++++++++++++++++++++++++++++++++++++++++++++++++++++++++
A
scratchpad/scratchpad
 | 
0

A
scratchpad/scratchpad.c
 | 
20
++++++++++++++++++++
M
squiggle.c
 | 
21
++++++++++++---------

4 files changed, 88 insertions(+), 9 deletions(-)
diff --git a/scratchpad/makefile 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
Binary files differ.
diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c
@@ -0,0 +1,20 @@
+#include "../squiggle.h"
+#include <math.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+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
@@ -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://stackoverflow.com/questions/53886131/how-does-xorshift32-works>
-    // https://en.wikipedia.org/wiki/Xorshift
-    // Also some drama: <https://www.pcg-random.org/posts/on-vignas-pcg-critique.html>, <https://prng.di.unimi.it/>
+    // See:
+    //   <https://en.wikipedia.org/wiki/Xorshift>
+    //   <https://stackoverflow.com/questions/53886131/how-does-xorshift32-works>,
+    // Also some drama: 
+    //   <https://www.pcg-random.org/posts/on-vignas-pcg-critique.html>, 
+    //   <https://prng.di.unimi.it/>
     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: <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>
-    double u1 = sample_unit_uniform(seed);
+    // // See: <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>
+    // 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)};