Compare commits
No commits in common. "2ea32e2a47f3b25d1214f47b7169a2fec28c57c3" and "91ab0d27dcc32d75eab795204a8095f5fed2ba7c" have entirely different histories.
2ea32e2a47
...
91ab0d27dc
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
|
@ -1,81 +0,0 @@
|
|||
#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);
|
||||
}
|
|
@ -1,53 +0,0 @@
|
|||
# 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
|
|
@ -134,11 +134,6 @@ 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)
|
||||
{
|
||||
|
|
|
@ -19,7 +19,6 @@ 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);
|
||||
|
|
Loading…
Reference in New Issue
Block a user