savepoint

examples/10_nuclear_recovery/example.c

@@ -0,0 +1,81 @@
+#include "../../squiggle.h"
+#include <math.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+double yearly_probability_nuclear_collapse(double year, uint64_t* seed)
+{
+	double successes = 0;
+	double failures = (year - 1960);
+	return sample_laplace(successes, failures, seed); 
+	// ^ can change to (successes + 1)/(trials + 2) 
+	// to get a probability, 
+	// rather than sampling from a distribution over probabilities.
+}
+double yearly_probability_nuclear_collapse_2023(uint64_t* seed){
+	return yearly_probability_nuclear_collapse(2023, seed);
+}
+
+double yearly_probability_nuclear_collapse_after_recovery(double year, double rebuilding_period_length_years, uint64_t* seed){
+	// assumption: nuclear 
+	double successes = 1.0;
+	double failures = (year - rebuilding_period_length_years - 1960 - 1);
+	return sample_laplace(successes, failures, seed);
+}
+double yearly_probability_nuclear_collapse_after_recovery_example(uint64_t* seed){
+	double year = 2070;
+	double rebuilding_period_length_years =30;
+	// So, there was a nuclear collapse in 2040, 
+	// then a recovery period of 30 years
+	// and it's now 2070
+	return yearly_probability_nuclear_collapse_after_recovery(year, rebuilding_period_length_years, seed);
+}
+
+double yearly_probability_nuclear_collapse_after_recovery_antiinductive(uint64_t* seed){
+	return yearly_probability_nuclear_collapse(2023, seed)/2;
+}
+
+int main()
+{
+    // set randomness seed
+    uint64_t* seed = malloc(sizeof(uint64_t));
+    *seed = 1000; // xorshift can't start with 0
+
+    int num_samples = 1000000;
+
+		// Before a first nuclear collapse
+		printf("## Before the first nuclear collapse\n");
+    struct c_i c_i_90_2023 = get_90_confidence_interval(yearly_probability_nuclear_collapse_2023, seed);
+    printf("90%% confidence interval: [%f, %f]\n", c_i_90_2023.low, c_i_90_2023.high);
+
+    double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * num_samples);
+    for (int i = 0; i < num_samples; i++) {
+        yearly_probability_nuclear_collapse_2023_samples[i] = yearly_probability_nuclear_collapse_2023(seed);
+    }
+		printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_2023_samples, num_samples));
+
+		// After the first nuclear collapse
+		printf("\n## After the first nuclear collapse\n");
+    struct c_i c_i_90_2070 = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_example, seed);
+    printf("90%% confidence interval: [%f, %f]\n", c_i_90_2070.low, c_i_90_2070.high);
+
+    double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * num_samples);
+    for (int i = 0; i < num_samples; i++) {
+        yearly_probability_nuclear_collapse_after_recovery_samples[i] = yearly_probability_nuclear_collapse_after_recovery_example(seed);
+    }
+		printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_after_recovery_samples, num_samples));
+
+		// After the first nuclear collapse (antiinductive)
+		printf("\n## After the first nuclear collapse (antiinductive)\n");
+    struct c_i c_i_90_antiinductive = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_antiinductive, seed);
+    printf("90%% confidence interval: [%f, %f]\n", c_i_90_antiinductive.low, c_i_90_antiinductive.high);
+
+    double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * num_samples);
+    for (int i = 0; i < num_samples; i++) {
+        yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples[i] = yearly_probability_nuclear_collapse_after_recovery_antiinductive(seed);
+    }
+		printf("mean: %f\n", array_mean(yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples, num_samples));
+
+    free(seed);
+}

examples/10_nuclear_recovery/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)
+	OMP_NUM_THREADS=1 ./$(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

squiggle.c

@@ -134,6 +134,11 @@ double sample_beta(double a, double b, uint64_t* seed)
     return gamma_a / (gamma_a + gamma_b);
 }
 
+double sample_laplace(double successes, double failures, uint64_t* seed){
+	// see <https://wikiless.esmailelbob.xyz/wiki/Beta_distribution?lang=en#Rule_of_succession>
+	return sample_beta(successes + 1, failures + 1, seed);
+}
+
 // Array helpers
 double array_sum(double* array, int length)
 {

squiggle.h

@@ -19,6 +19,7 @@ double sample_to(double low, double high, uint64_t* seed);
 
 double sample_gamma(double alpha, uint64_t* seed);
 double sample_beta(double a, double b, uint64_t* seed);
+double sample_laplace(double successes, double failures, uint64_t* seed);
 
 // Array helpers
 double array_sum(double* array, int length);