forked from personal/squiggle.c
refactor: struct box => box. Through typedef.
This commit is contained in:
parent
e08ce4334e
commit
edb58bdc1d
Binary file not shown.
Binary file not shown.
|
@ -40,7 +40,7 @@ double cdf_normal_0_1(double x)
|
|||
// Some testers
|
||||
void test_inverse_cdf_double(char* cdf_name, double cdf_double(double))
|
||||
{
|
||||
struct box result = inverse_cdf_double(cdf_double, 0.5);
|
||||
box result = inverse_cdf_double(cdf_double, 0.5);
|
||||
if (result.empty) {
|
||||
printf("Inverse for %s not calculated\n", cdf_name);
|
||||
exit(1);
|
||||
|
@ -54,7 +54,7 @@ void test_and_time_sampler_double(char* cdf_name, double cdf_double(double), uin
|
|||
printf("\nGetting some samples from %s:\n", cdf_name);
|
||||
clock_t begin = clock();
|
||||
for (int i = 0; i < NUM_SAMPLES; i++) {
|
||||
struct box sample = sampler_cdf_double(cdf_double, seed);
|
||||
box sample = sampler_cdf_double(cdf_double, seed);
|
||||
if (sample.empty) {
|
||||
printf("Error in sampler function for %s", cdf_name);
|
||||
} else {
|
||||
|
|
Binary file not shown.
|
@ -10,7 +10,7 @@
|
|||
#define TINY_BETA 1.0e-30
|
||||
|
||||
// Incomplete beta function
|
||||
struct box incbeta(double a, double b, double x)
|
||||
box incbeta(double a, double b, double x)
|
||||
{
|
||||
// Descended from <https://github.com/codeplea/incbeta/blob/master/incbeta.c>,
|
||||
// <https://codeplea.com/incomplete-beta-function-c>
|
||||
|
@ -47,11 +47,11 @@ struct box incbeta(double a, double b, double x)
|
|||
|
||||
/*The continued fraction converges nicely for x < (a+1)/(a+b+2)*/
|
||||
if (x > (a + 1.0) / (a + b + 2.0)) {
|
||||
struct box symmetric_incbeta = incbeta(b, a, 1.0 - x);
|
||||
box symmetric_incbeta = incbeta(b, a, 1.0 - x);
|
||||
if (symmetric_incbeta.empty) {
|
||||
return symmetric_incbeta; // propagate error
|
||||
} else {
|
||||
struct box result = {
|
||||
box result = {
|
||||
.empty = 0,
|
||||
.content = 1 - symmetric_incbeta.content
|
||||
};
|
||||
|
@ -94,7 +94,7 @@ struct box incbeta(double a, double b, double x)
|
|||
|
||||
/*Check for stop.*/
|
||||
if (fabs(1.0 - cd) < STOP_BETA) {
|
||||
struct box result = {
|
||||
box result = {
|
||||
.empty = 0,
|
||||
.content = front * (f - 1.0)
|
||||
};
|
||||
|
@ -105,13 +105,13 @@ struct box incbeta(double a, double b, double x)
|
|||
return PROCESS_ERROR("More loops needed, did not converge, in function incbeta");
|
||||
}
|
||||
|
||||
struct box cdf_beta(double x)
|
||||
box cdf_beta(double x)
|
||||
{
|
||||
if (x < 0) {
|
||||
struct box result = { .empty = 0, .content = 0 };
|
||||
box result = { .empty = 0, .content = 0 };
|
||||
return result;
|
||||
} else if (x > 1) {
|
||||
struct box result = { .empty = 0, .content = 1 };
|
||||
box result = { .empty = 0, .content = 1 };
|
||||
return result;
|
||||
} else {
|
||||
double successes = 1, failures = (2023 - 1945);
|
||||
|
@ -120,9 +120,9 @@ struct box cdf_beta(double x)
|
|||
}
|
||||
|
||||
// Some testers
|
||||
void test_inverse_cdf_box(char* cdf_name, struct box cdf_box(double))
|
||||
void test_inverse_cdf_box(char* cdf_name, box cdf_box(double))
|
||||
{
|
||||
struct box result = inverse_cdf_box(cdf_box, 0.5);
|
||||
box result = inverse_cdf_box(cdf_box, 0.5);
|
||||
if (result.empty) {
|
||||
printf("Inverse for %s not calculated\n", cdf_name);
|
||||
exit(1);
|
||||
|
@ -131,12 +131,12 @@ void test_inverse_cdf_box(char* cdf_name, struct box cdf_box(double))
|
|||
}
|
||||
}
|
||||
|
||||
void test_and_time_sampler_box(char* cdf_name, struct box cdf_box(double), uint64_t* seed)
|
||||
void test_and_time_sampler_box(char* cdf_name, box cdf_box(double), uint64_t* seed)
|
||||
{
|
||||
printf("\nGetting some samples from %s:\n", cdf_name);
|
||||
clock_t begin = clock();
|
||||
for (int i = 0; i < NUM_SAMPLES; i++) {
|
||||
struct box sample = sampler_cdf_box(cdf_box, seed);
|
||||
box sample = sampler_cdf_box(cdf_box, seed);
|
||||
if (sample.empty) {
|
||||
printf("Error in sampler function for %s", cdf_name);
|
||||
} else {
|
||||
|
|
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
|
@ -11,14 +11,14 @@
|
|||
void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_threads, int n_samples)
|
||||
{
|
||||
|
||||
// Division terminology:
|
||||
// Terms of the division:
|
||||
// a = b * quotient + reminder
|
||||
// a = (a/b)*b + (a%b)
|
||||
// a = b * (a/b) + (a%b)
|
||||
// dividend: a
|
||||
// divisor: b
|
||||
// quotient = a/b
|
||||
// reminder = a%b
|
||||
// "divisor's multiple" := (a/b)*b
|
||||
// "divisor's multiple" := b*(a/b)
|
||||
|
||||
// now, we have n_samples and n_threads
|
||||
// to make our life easy, each thread will have a number of samples of: a/b (quotient)
|
||||
|
@ -26,11 +26,9 @@ void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_
|
|||
// to possibly do by Jorge: improve so that the remainder is included in the threads
|
||||
|
||||
int quotient = n_samples / n_threads;
|
||||
/* int remainder = n_samples % n_threads; // not used, comment to avoid lint warning */
|
||||
int divisor_multiple = quotient * n_threads;
|
||||
|
||||
uint64_t** seeds = malloc(n_threads * sizeof(uint64_t*));
|
||||
// printf("UINT64_MAX: %lu\n", UINT64_MAX);
|
||||
srand(1);
|
||||
for (uint64_t i = 0; i < n_threads; i++) {
|
||||
seeds[i] = malloc(sizeof(uint64_t));
|
||||
|
@ -53,7 +51,6 @@ void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_
|
|||
for (i = 0; i < n_threads; i++) {
|
||||
int lower_bound_inclusive = i * quotient;
|
||||
int upper_bound_not_inclusive = ((i + 1) * quotient); // note the < in the for loop below,
|
||||
// printf("Lower bound: %d, upper bound: %d\n", lower_bound, upper_bound);
|
||||
for (int j = lower_bound_inclusive; j < upper_bound_not_inclusive; j++) {
|
||||
results[j] = sampler(seeds[i]);
|
||||
}
|
||||
|
@ -71,9 +68,6 @@ void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_
|
|||
}
|
||||
|
||||
/* Get confidence intervals, given a sampler */
|
||||
// Not in core yet because I'm not sure how much I like the struct
|
||||
// and the built-in 100k samples
|
||||
// to do: add n to function parameters and document
|
||||
|
||||
typedef struct ci_t {
|
||||
double low;
|
||||
|
@ -89,9 +83,7 @@ static void swp(int i, int j, double xs[])
|
|||
|
||||
static int partition(int low, int high, double xs[], int length)
|
||||
{
|
||||
// To understand this function:
|
||||
// - see the note after gt variable definition
|
||||
// - go to commit 578bfa27 and the scratchpad/ folder in it, which has printfs sprinkled throughout
|
||||
// Note: the scratchpad/ folder in commit 578bfa27 has printfs sprinkled throughout
|
||||
int pivot = low + floor((high - low) / 2);
|
||||
double pivot_value = xs[pivot];
|
||||
swp(pivot, high, xs);
|
||||
|
@ -145,9 +137,6 @@ ci array_get_90_ci(double xs[], int 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);
|
||||
}*/
|
||||
sampler_parallel(sampler, xs, 16, n);
|
||||
ci result = array_get_ci(interval, xs, n);
|
||||
free(xs);
|
||||
|
@ -159,9 +148,6 @@ ci sampler_get_90_ci(double (*sampler)(uint64_t*), int n, uint64_t* seed)
|
|||
}
|
||||
|
||||
/* Algebra manipulations */
|
||||
// here I discover named structs,
|
||||
// which mean that I don't have to be typing
|
||||
// struct blah all the time.
|
||||
|
||||
#define NORMAL90CONFIDENCE 1.6448536269514727
|
||||
|
||||
|
@ -221,13 +207,13 @@ ci convert_lognormal_params_to_ci(lognormal_params y)
|
|||
#define EXIT_ON_ERROR 0
|
||||
#define PROCESS_ERROR(error_msg) process_error(error_msg, EXIT_ON_ERROR, __FILE__, __LINE__)
|
||||
|
||||
struct box {
|
||||
typedef struct box_t {
|
||||
int empty;
|
||||
double content;
|
||||
char* error_msg;
|
||||
};
|
||||
} box;
|
||||
|
||||
struct box process_error(const char* error_msg, int should_exit, char* file, int line)
|
||||
box process_error(const char* error_msg, int should_exit, char* file, int line)
|
||||
{
|
||||
if (should_exit) {
|
||||
printf("@, in %s (%d)", file, line);
|
||||
|
@ -235,7 +221,7 @@ struct box process_error(const char* error_msg, int should_exit, char* file, int
|
|||
} else {
|
||||
char error_msg[MAX_ERROR_LENGTH];
|
||||
snprintf(error_msg, MAX_ERROR_LENGTH, "@, in %s (%d)", file, line); // NOLINT: We are being carefull here by considering MAX_ERROR_LENGTH explicitly.
|
||||
struct box error = { .empty = 1, .error_msg = error_msg };
|
||||
box error = { .empty = 1, .error_msg = error_msg };
|
||||
return error;
|
||||
}
|
||||
}
|
||||
|
@ -244,7 +230,7 @@ struct box process_error(const char* error_msg, int should_exit, char* file, int
|
|||
// Version #1:
|
||||
// - input: (cdf: double => double, p)
|
||||
// - output: Box(number|error)
|
||||
struct box inverse_cdf_double(double cdf(double), double p)
|
||||
box inverse_cdf_double(double cdf(double), double p)
|
||||
{
|
||||
// given a cdf: [-Inf, Inf] => [0,1]
|
||||
// returns a box with either
|
||||
|
@ -299,7 +285,7 @@ struct box inverse_cdf_double(double cdf(double), double p)
|
|||
}
|
||||
|
||||
if (convergence_condition) {
|
||||
struct box result = { .empty = 0, .content = low };
|
||||
box result = { .empty = 0, .content = low };
|
||||
return result;
|
||||
} else {
|
||||
return PROCESS_ERROR("Search process did not converge, in function inverse_cdf");
|
||||
|
@ -310,7 +296,7 @@ struct box inverse_cdf_double(double cdf(double), double p)
|
|||
// Version #2:
|
||||
// - input: (cdf: double => Box(number|error), p)
|
||||
// - output: Box(number|error)
|
||||
struct box inverse_cdf_box(struct box cdf_box(double), double p)
|
||||
box inverse_cdf_box(box cdf_box(double), double p)
|
||||
{
|
||||
// given a cdf: [-Inf, Inf] => Box([0,1])
|
||||
// returns a box with either
|
||||
|
@ -326,12 +312,12 @@ struct box inverse_cdf_box(struct box cdf_box(double), double p)
|
|||
while ((!interval_found) && (low > -FLT_MAX / 4) && (high < FLT_MAX / 4)) {
|
||||
// ^ Using FLT_MIN and FLT_MAX is overkill
|
||||
// but it's also the *correct* thing to do.
|
||||
struct box cdf_low = cdf_box(low);
|
||||
box cdf_low = cdf_box(low);
|
||||
if (cdf_low.empty) {
|
||||
return PROCESS_ERROR(cdf_low.error_msg);
|
||||
}
|
||||
|
||||
struct box cdf_high = cdf_box(high);
|
||||
box cdf_high = cdf_box(high);
|
||||
if (cdf_high.empty) {
|
||||
return PROCESS_ERROR(cdf_low.error_msg);
|
||||
}
|
||||
|
@ -361,7 +347,7 @@ struct box inverse_cdf_box(struct box cdf_box(double), double p)
|
|||
// if ((width < 1e-8) || mid_not_new){
|
||||
convergence_condition = 1;
|
||||
} else {
|
||||
struct box cdf_mid = cdf_box(mid);
|
||||
box cdf_mid = cdf_box(mid);
|
||||
if (cdf_mid.empty) {
|
||||
return PROCESS_ERROR(cdf_mid.error_msg);
|
||||
}
|
||||
|
@ -378,7 +364,7 @@ struct box inverse_cdf_box(struct box cdf_box(double), double p)
|
|||
}
|
||||
|
||||
if (convergence_condition) {
|
||||
struct box result = { .empty = 0, .content = low };
|
||||
box result = { .empty = 0, .content = low };
|
||||
return result;
|
||||
} else {
|
||||
return PROCESS_ERROR("Search process did not converge, in function inverse_cdf");
|
||||
|
@ -389,22 +375,22 @@ struct box inverse_cdf_box(struct box cdf_box(double), double p)
|
|||
/* Sample from an arbitrary cdf */
|
||||
// Before: invert an arbitrary cdf at a point
|
||||
// Now: from an arbitrary cdf, get a sample
|
||||
struct box sampler_cdf_box(struct box cdf(double), uint64_t* seed)
|
||||
box sampler_cdf_box(box cdf(double), uint64_t* seed)
|
||||
{
|
||||
double p = sample_unit_uniform(seed);
|
||||
struct box result = inverse_cdf_box(cdf, p);
|
||||
box result = inverse_cdf_box(cdf, p);
|
||||
return result;
|
||||
}
|
||||
struct box sampler_cdf_double(double cdf(double), uint64_t* seed)
|
||||
box sampler_cdf_double(double cdf(double), uint64_t* seed)
|
||||
{
|
||||
double p = sample_unit_uniform(seed);
|
||||
struct box result = inverse_cdf_double(cdf, p);
|
||||
box result = inverse_cdf_double(cdf, p);
|
||||
return result;
|
||||
}
|
||||
double sampler_cdf_danger(struct box cdf(double), uint64_t* seed)
|
||||
double sampler_cdf_danger(box cdf(double), uint64_t* seed)
|
||||
{
|
||||
double p = sample_unit_uniform(seed);
|
||||
struct box result = inverse_cdf_box(cdf, p);
|
||||
box result = inverse_cdf_box(cdf, p);
|
||||
if (result.empty) {
|
||||
exit(1);
|
||||
} else {
|
||||
|
@ -413,7 +399,6 @@ double sampler_cdf_danger(struct box cdf(double), uint64_t* seed)
|
|||
}
|
||||
|
||||
/* array print: potentially useful for debugging */
|
||||
|
||||
void array_print(double xs[], int n)
|
||||
{
|
||||
printf("[");
|
||||
|
|
|
@ -32,23 +32,23 @@ lognormal_params convert_ci_to_lognormal_params(ci x);
|
|||
ci convert_lognormal_params_to_ci(lognormal_params y);
|
||||
|
||||
/* Error handling */
|
||||
struct box {
|
||||
typedef struct box_t {
|
||||
int empty;
|
||||
double content;
|
||||
char* error_msg;
|
||||
};
|
||||
} box;
|
||||
#define MAX_ERROR_LENGTH 500
|
||||
#define EXIT_ON_ERROR 0
|
||||
#define PROCESS_ERROR(error_msg) process_error(error_msg, EXIT_ON_ERROR, __FILE__, __LINE__)
|
||||
struct box process_error(const char* error_msg, int should_exit, char* file, int line);
|
||||
box process_error(const char* error_msg, int should_exit, char* file, int line);
|
||||
void array_print(double* array, int length);
|
||||
|
||||
/* Inverse cdf */
|
||||
struct box inverse_cdf_double(double cdf(double), double p);
|
||||
struct box inverse_cdf_box(struct box cdf_box(double), double p);
|
||||
box inverse_cdf_double(double cdf(double), double p);
|
||||
box inverse_cdf_box(box cdf_box(double), double p);
|
||||
|
||||
/* Samplers from cdf */
|
||||
struct box sampler_cdf_double(double cdf(double), uint64_t* seed);
|
||||
struct box sampler_cdf_box(struct box cdf(double), uint64_t* seed);
|
||||
box sampler_cdf_double(double cdf(double), uint64_t* seed);
|
||||
box sampler_cdf_box(box cdf(double), uint64_t* seed);
|
||||
|
||||
#endif
|
||||
|
|
Loading…
Reference in New Issue
Block a user