Compare commits

...

5 Commits

10 changed files with 140 additions and 0 deletions

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

View File

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

View File

@ -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

View File

@ -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)
{

View File

@ -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);