tweaks before twitter thread

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

examples/02_many_samples_time_to_botec/example.c

@@ -1,7 +1,7 @@
-#include <stdint.h>
-#include <stdlib.h>
-#include <stdio.h>
 #include "../../squiggle.h"
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
 
 // Estimate functions
 double sample_0(uint64_t* seed)
@@ -24,7 +24,8 @@ 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
@@ -38,8 +39,8 @@ int main(){
     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<n_samples; i++){
+    double* result_many = (double*)malloc(n_samples * sizeof(double));
+    for (int i = 0; i < n_samples; i++) {
         result_many[i] = sample_mixture(samplers, weights, n_dists, seed);
     }
     printf("Mean: %f\n", array_mean(result_many, n_samples));

examples/14_twitter_thread_example/example.c

@@ -0,0 +1,43 @@
+#include "../../squiggle.h"
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+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);
+}

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

scratchpad/core.c

@@ -0,0 +1,27 @@
+
+uint64_t xorshift64(uint64_t* seed)
+{
+    // Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs"
+    // <https://en.wikipedia.org/wiki/Xorshift>
+    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: <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;
+}
+

squiggle.c

@@ -57,7 +57,7 @@ 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);
+    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;
-
 }
-