forked from personal/squiggle.c
refactor & recompile for function definitions
This commit is contained in:
parent
023c9f28ac
commit
fb110a35f3
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
|
@ -15,7 +15,7 @@ int main()
|
||||||
uint64_t* seed = malloc(sizeof(uint64_t));
|
uint64_t* seed = malloc(sizeof(uint64_t));
|
||||||
*seed = 1000; // xorshift can't start with 0
|
*seed = 1000; // xorshift can't start with 0
|
||||||
|
|
||||||
ci beta_1_2_ci_90 = get_90_confidence_interval(beta_1_2_sampler, seed);
|
ci beta_1_2_ci_90 = sampler_get_90_ci(beta_1_2_sampler, 1000000, seed);
|
||||||
printf("90%% confidence interval of beta(1,2) is [%f, %f]\n", beta_1_2_ci_90.low, beta_1_2_ci_90.high);
|
printf("90%% confidence interval of beta(1,2) is [%f, %f]\n", beta_1_2_ci_90.low, beta_1_2_ci_90.high);
|
||||||
|
|
||||||
free(seed);
|
free(seed);
|
||||||
|
|
Binary file not shown.
|
@ -60,7 +60,7 @@ int main()
|
||||||
}
|
}
|
||||||
printf("... ]\n");
|
printf("... ]\n");
|
||||||
|
|
||||||
ci ci_90 = get_90_confidence_interval(mixture, seed);
|
ci ci_90 = sampler_get_90_ci(mixture, 1000000, seed);
|
||||||
printf("mean: %f\n", array_mean(mixture_result, n));
|
printf("mean: %f\n", array_mean(mixture_result, n));
|
||||||
printf("90%% confidence interval: [%f, %f]\n", ci_90.low, ci_90.high);
|
printf("90%% confidence interval: [%f, %f]\n", ci_90.low, ci_90.high);
|
||||||
|
|
||||||
|
|
Binary file not shown.
Binary file not shown.
|
@ -41,7 +41,7 @@ int main()
|
||||||
}
|
}
|
||||||
printf("... ]\n");
|
printf("... ]\n");
|
||||||
|
|
||||||
ci ci_90 = get_90_confidence_interval(sample_minutes_per_day_jumping_rope_needed_to_burn_10kg, seed);
|
ci ci_90 = sampler_get_90_ci(sample_minutes_per_day_jumping_rope_needed_to_burn_10kg, 1000000, seed);
|
||||||
printf("90%% confidence interval: [%f, %f]\n", ci_90.low, ci_90.high);
|
printf("90%% confidence interval: [%f, %f]\n", ci_90.low, ci_90.high);
|
||||||
|
|
||||||
free(seed);
|
free(seed);
|
||||||
|
|
Binary file not shown.
|
@ -50,7 +50,7 @@ int main()
|
||||||
|
|
||||||
// Before a first nuclear collapse
|
// Before a first nuclear collapse
|
||||||
printf("## Before the first nuclear collapse\n");
|
printf("## Before the first nuclear collapse\n");
|
||||||
ci ci_90_2023 = get_90_confidence_interval(yearly_probability_nuclear_collapse_2023, seed);
|
ci ci_90_2023 = sampler_get_90_ci(yearly_probability_nuclear_collapse_2023, 1000000, seed);
|
||||||
printf("90%% confidence interval: [%f, %f]\n", ci_90_2023.low, ci_90_2023.high);
|
printf("90%% confidence interval: [%f, %f]\n", ci_90_2023.low, ci_90_2023.high);
|
||||||
|
|
||||||
double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * num_samples);
|
double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * num_samples);
|
||||||
|
@ -61,7 +61,7 @@ int main()
|
||||||
|
|
||||||
// After the first nuclear collapse
|
// After the first nuclear collapse
|
||||||
printf("\n## After the first nuclear collapse\n");
|
printf("\n## After the first nuclear collapse\n");
|
||||||
ci ci_90_2070 = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_example, seed);
|
ci ci_90_2070 = sampler_get_90_ci(yearly_probability_nuclear_collapse_after_recovery_example, 1000000, seed);
|
||||||
printf("90%% confidence interval: [%f, %f]\n", ci_90_2070.low, ci_90_2070.high);
|
printf("90%% confidence interval: [%f, %f]\n", ci_90_2070.low, ci_90_2070.high);
|
||||||
|
|
||||||
double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * num_samples);
|
double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * num_samples);
|
||||||
|
@ -72,7 +72,7 @@ int main()
|
||||||
|
|
||||||
// After the first nuclear collapse (antiinductive)
|
// After the first nuclear collapse (antiinductive)
|
||||||
printf("\n## After the first nuclear collapse (antiinductive)\n");
|
printf("\n## After the first nuclear collapse (antiinductive)\n");
|
||||||
ci ci_90_antiinductive = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_antiinductive, seed);
|
ci ci_90_antiinductive = sampler_get_90_ci(yearly_probability_nuclear_collapse_after_recovery_antiinductive, 1000000, seed);
|
||||||
printf("90%% confidence interval: [%f, %f]\n", ci_90_antiinductive.low, ci_90_antiinductive.high);
|
printf("90%% confidence interval: [%f, %f]\n", ci_90_antiinductive.low, ci_90_antiinductive.high);
|
||||||
|
|
||||||
double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * num_samples);
|
double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * num_samples);
|
||||||
|
|
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
|
@ -9,7 +9,7 @@ int main()
|
||||||
// set randomness seed
|
// set randomness seed
|
||||||
// uint64_t* seed = malloc(sizeof(uint64_t));
|
// uint64_t* seed = malloc(sizeof(uint64_t));
|
||||||
// *seed = 1000; // xorshift can't start with 0
|
// *seed = 1000; // xorshift can't start with 0
|
||||||
// ^ not necessary, because parallel_sampler takes care of the seed.
|
// ^ not necessary, because sampler_parallel takes care of the seed.
|
||||||
|
|
||||||
int n_samples = 1000 * 1000 * 1000;
|
int n_samples = 1000 * 1000 * 1000;
|
||||||
int n_threads = 16;
|
int n_threads = 16;
|
||||||
|
@ -18,12 +18,12 @@ int main()
|
||||||
}
|
}
|
||||||
double* results = malloc(n_samples * sizeof(double));
|
double* results = malloc(n_samples * sizeof(double));
|
||||||
|
|
||||||
parallel_sampler(sampler, results, n_threads, n_samples);
|
sampler_parallel(sampler, results, n_threads, n_samples);
|
||||||
double avg = array_sum(results, n_samples)/n_samples;
|
double avg = array_sum(results, n_samples)/n_samples;
|
||||||
printf("Average of 1B lognormal(0,10): %f", avg);
|
printf("Average of 1B lognormal(0,10): %f", avg);
|
||||||
|
|
||||||
free(results);
|
free(results);
|
||||||
|
|
||||||
// free(seed);
|
// free(seed);
|
||||||
// ^ not necessary, because parallel_sampler takes care of the seed.
|
// ^ not necessary, because sampler_parallel takes care of the seed.
|
||||||
}
|
}
|
||||||
|
|
Binary file not shown.
|
@ -23,7 +23,7 @@ int main()
|
||||||
|
|
||||||
int n_samples = 1000 * 1000, n_threads = 16;
|
int n_samples = 1000 * 1000, n_threads = 16;
|
||||||
double* results = malloc(n_samples * sizeof(double));
|
double* results = malloc(n_samples * sizeof(double));
|
||||||
parallel_sampler(sampler_result, results, n_threads, n_samples);
|
sampler_parallel(sampler_result, results, n_threads, n_samples);
|
||||||
printf("Avg: %f\n", array_sum(results, n_samples)/n_samples);
|
printf("Avg: %f\n", array_sum(results, n_samples)/n_samples);
|
||||||
free(results);
|
free(results);
|
||||||
}
|
}
|
||||||
|
|
Binary file not shown.
|
@ -29,7 +29,7 @@ int main()
|
||||||
|
|
||||||
int n_samples = 1000000, n_threads = 16;
|
int n_samples = 1000000, n_threads = 16;
|
||||||
double* results = malloc(n_samples * sizeof(double));
|
double* results = malloc(n_samples * sizeof(double));
|
||||||
parallel_sampler(sampler_result, results, n_threads, n_samples);
|
sampler_parallel(sampler_result, results, n_threads, n_samples);
|
||||||
printf("Mean of the distribution of (taking the min of 1000 samples of a normal(5,2)): %f\n", array_mean(results, n_samples));
|
printf("Mean of the distribution of (taking the min of 1000 samples of a normal(5,2)): %f\n", array_mean(results, n_samples));
|
||||||
free(results);
|
free(results);
|
||||||
|
|
||||||
|
@ -48,7 +48,7 @@ int main()
|
||||||
return sample_min_of_n(seed, quotient);
|
return sample_min_of_n(seed, quotient);
|
||||||
}
|
}
|
||||||
double* results_quotient = malloc(quotient * sizeof(double));
|
double* results_quotient = malloc(quotient * sizeof(double));
|
||||||
parallel_sampler(sample_min_of_quotient, results_quotient, n_threads, quotient);
|
sampler_parallel(sample_min_of_quotient, results_quotient, n_threads, quotient);
|
||||||
|
|
||||||
double min = results_quotient[0];
|
double min = results_quotient[0];
|
||||||
for(int i=1; i<quotient; i++){
|
for(int i=1; i<quotient; i++){
|
||||||
|
|
Binary file not shown.
|
@ -8,7 +8,7 @@
|
||||||
#include "squiggle.h"
|
#include "squiggle.h"
|
||||||
|
|
||||||
/* Parallel sampler */
|
/* Parallel sampler */
|
||||||
void parallel_sampler(double (*sampler)(uint64_t* seed), double* results, int n_threads, int n_samples){
|
void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_threads, int n_samples){
|
||||||
if((n_samples % n_threads) != 0){
|
if((n_samples % n_threads) != 0){
|
||||||
fprintf(stderr, "Number of samples isn't divisible by number of threads, aborting\n");
|
fprintf(stderr, "Number of samples isn't divisible by number of threads, aborting\n");
|
||||||
exit(1);
|
exit(1);
|
||||||
|
@ -45,8 +45,8 @@ void parallel_sampler(double (*sampler)(uint64_t* seed), double* results, int n_
|
||||||
// to do: add n to function parameters and document
|
// to do: add n to function parameters and document
|
||||||
|
|
||||||
typedef struct ci_t {
|
typedef struct ci_t {
|
||||||
float low;
|
double low;
|
||||||
float high;
|
double high;
|
||||||
} ci;
|
} ci;
|
||||||
|
|
||||||
static void swp(int i, int j, double xs[])
|
static void swp(int i, int j, double xs[])
|
||||||
|
@ -75,16 +75,16 @@ static int partition(int low, int high, double xs[], int length)
|
||||||
return gt;
|
return gt;
|
||||||
}
|
}
|
||||||
|
|
||||||
static double quickselect(int k, double xs[], int length)
|
static double quickselect(int k, double xs[], int n)
|
||||||
{
|
{
|
||||||
// https://en.wikipedia.org/wiki/Quickselect
|
// https://en.wikipedia.org/wiki/Quickselect
|
||||||
int low = 0;
|
int low = 0;
|
||||||
int high = length - 1;
|
int high = n - 1;
|
||||||
for (;;) {
|
for (;;) {
|
||||||
if (low == high) {
|
if (low == high) {
|
||||||
return xs[low];
|
return xs[low];
|
||||||
}
|
}
|
||||||
int pivot = partition(low, high, xs, length);
|
int pivot = partition(low, high, xs, n);
|
||||||
if (pivot == k) {
|
if (pivot == k) {
|
||||||
return xs[pivot];
|
return xs[pivot];
|
||||||
} else if (k < pivot) {
|
} else if (k < pivot) {
|
||||||
|
@ -95,27 +95,34 @@ static double quickselect(int k, double xs[], int length)
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
ci sampler_get_ci(double (*sampler)(uint64_t*), ci interval, int n, uint64_t* seed){
|
ci array_get_ci(ci interval, double* xs, int n){
|
||||||
double* xs = malloc(n * sizeof(double));
|
|
||||||
for (int i = 0; i < n; i++) {
|
|
||||||
xs[i] = sampler(seed);
|
|
||||||
}
|
|
||||||
|
|
||||||
int low_k = floor(interval.low * n);
|
int low_k = floor(interval.low * n);
|
||||||
int high_k = ceil(interval.high * n);
|
int high_k = ceil(interval.high * n);
|
||||||
|
|
||||||
ci result = {
|
ci result = {
|
||||||
.low = quickselect(low_k, xs, n),
|
.low = quickselect(low_k, xs, n),
|
||||||
.high = quickselect(high_k, xs, n),
|
.high = quickselect(high_k, xs, n),
|
||||||
};
|
};
|
||||||
|
return result;
|
||||||
|
}
|
||||||
|
ci array_get_90_ci(double xs[], int n)
|
||||||
|
{
|
||||||
|
return array_get_ci((ci) {.low = 0.05, .high = 0.95}, xs, n);
|
||||||
|
}
|
||||||
|
|
||||||
|
ci sampler_get_ci(ci interval, double (*sampler)(uint64_t*), int n, uint64_t* seed){
|
||||||
|
double* xs = malloc(n * sizeof(double));
|
||||||
|
for (int i = 0; i < n; i++) {
|
||||||
|
xs[i] = sampler(seed);
|
||||||
|
}
|
||||||
|
ci result = array_get_ci(interval, xs, n);
|
||||||
free(xs);
|
free(xs);
|
||||||
return result;
|
return result;
|
||||||
|
|
||||||
}
|
}
|
||||||
|
ci sampler_get_90_ci(double (*sampler)(uint64_t*), int n, uint64_t* seed)
|
||||||
ci get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed)
|
|
||||||
{
|
{
|
||||||
return sampler_get_ci(sampler, (ci) {.low = 0.05, .high = 0.95}, 1000000, seed);
|
return sampler_get_ci((ci) {.low = 0.05, .high = 0.95}, sampler, n, seed);
|
||||||
}
|
}
|
||||||
|
|
||||||
/* Algebra manipulations */
|
/* Algebra manipulations */
|
||||||
|
|
|
@ -2,14 +2,14 @@
|
||||||
#define SQUIGGLE_C_EXTRA
|
#define SQUIGGLE_C_EXTRA
|
||||||
|
|
||||||
/* Parallel sampling */
|
/* Parallel sampling */
|
||||||
void parallel_sampler(double (*sampler)(uint64_t* seed), double* results, int n_threads, int n_samples);
|
void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_threads, int n_samples);
|
||||||
|
|
||||||
/* Get 90% confidence interval */
|
/* Get 90% confidence interval */
|
||||||
typedef struct ci_t {
|
typedef struct ci_t {
|
||||||
float low;
|
float low;
|
||||||
float high;
|
float high;
|
||||||
} ci;
|
} ci;
|
||||||
ci get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed);
|
ci sampler_get_90_ci(double (*sampler)(uint64_t*), int n, uint64_t* seed);
|
||||||
|
|
||||||
/* Algebra manipulations */
|
/* Algebra manipulations */
|
||||||
|
|
||||||
|
|
Loading…
Reference in New Issue
Block a user