diff --git a/README.md b/README.md index 960eaf2a..0bf64330 100644 --- a/README.md +++ b/README.md @@ -25,7 +25,7 @@ The name of this repository is a pun on two meanings of "time to": "how much tim | Language | Time | Lines of code | |-----------------------------|-----------|---------------| | C | 5.6ms | 252 | -| squiggle.c | 12.7ms | 29* | +| squiggle.c | 10.5ms | 29* | | Nim | 40.8ms | 84 | | Lua (LuaJIT) | 69.9ms | 82 | | OCaml (flambda) | 187.9ms | 123 | diff --git a/squiggle.c/samples b/squiggle.c/samples index 409e05f1..d2edfd2a 100755 Binary files a/squiggle.c/samples and b/squiggle.c/samples differ diff --git a/squiggle.c/samples.c b/squiggle.c/samples.c index a385cb07..33e21674 100644 --- a/squiggle.c/samples.c +++ b/squiggle.c/samples.c @@ -9,8 +9,8 @@ int main() double p_b = 0.5; double p_c = p_a * p_b; - double sample_0(uint64_t * seed) { return 0; } - double sample_1(uint64_t * seed) { return 1; } + double sample_0(uint64_t * seed) { UNUSED(seed); return 0; } + double sample_1(uint64_t * seed) { UNUSED(seed); return 1; } double sample_few(uint64_t * seed) { return sample_to(1, 3, seed); } double sample_many(uint64_t * seed) { return sample_to(2, 10, seed); } @@ -22,8 +22,8 @@ int main() return sample_mixture(samplers, weights, n_dists, seed); } - int n_samples = 1000000, n_threads = 16; - double* results = malloc(n_samples * sizeof(double)); + int n_samples = 1000 * 1000, n_threads = 16; + double* results = malloc((size_t)n_samples * sizeof(double)); sampler_parallel(sampler_result, results, n_threads, n_samples); printf("Avg: %f\n", array_sum(results, n_samples) / n_samples); free(results); diff --git a/squiggle.c/squiggle_c/squiggle.c b/squiggle.c/squiggle_c/squiggle.c index 1d0224c3..ec0d4a23 100644 --- a/squiggle.c/squiggle_c/squiggle.c +++ b/squiggle.c/squiggle_c/squiggle.c @@ -3,12 +3,14 @@ #include #include -// math constants +// Defs #define PI 3.14159265358979323846 // M_PI in gcc gnu99 #define NORMAL90CONFIDENCE 1.6448536269514727 +#define UNUSED(x) (void)(x) +// ^ https://stackoverflow.com/questions/3599160/how-can-i-suppress-unused-parameter-warnings-in-c -// Pseudo Random number generator -static uint64_t xorshift32(uint32_t* seed) +// Pseudo Random number generators +static uint64_t xorshift64(uint64_t* seed) { // Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" // See: @@ -19,19 +21,20 @@ static uint64_t xorshift32(uint32_t* seed) // uint64_t x = *seed; x ^= x << 13; - x ^= x >> 17; - x ^= x << 5; - return *seed = x; -} - -static uint64_t xorshift64(uint64_t* seed) -{ - // same as above, but for generating doubles instead of floats - uint64_t x = *seed; - x ^= x << 13; x ^= x >> 7; x ^= x << 17; return *seed = x; + + /* + // if one wanted to generate 32 bit ints, + // from which to generate floats, + // one could do the following: + uint32_t x = *seed; + x ^= x << 13; + x ^= x >> 17; + x ^= x << 5; + return *seed = x; + */ } // Distribution & sampling functions @@ -47,7 +50,7 @@ double sample_unit_normal(uint64_t* seed) // // See: double u1 = sample_unit_uniform(seed); double u2 = sample_unit_uniform(seed); - double z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2); + double z = sqrt(-2.0 * log(u1)) * sin(2 * PI * u2); return z; } @@ -67,7 +70,7 @@ double sample_lognormal(double logmean, double logstd, uint64_t* seed) return exp(sample_normal(logmean, logstd, seed)); } -inline double sample_normal_from_90_confidence_interval(double low, double high, uint64_t* seed) +double sample_normal_from_90_ci(double low, double high, uint64_t* seed) { // Explanation of key idea: // 1. We know that the 90% confidence interval of the unit normal is @@ -98,10 +101,10 @@ double sample_to(double low, double high, uint64_t* seed) // returns a sample from a lognorma with a matching 90% c.i. // Key idea: If we want a lognormal with 90% confidence interval [a, b] // we need but get a normal with 90% confidence interval [log(a), log(b)]. - // Then see code for sample_normal_from_90_confidence_interval - double loglow = logf(low); - double loghigh = logf(high); - return exp(sample_normal_from_90_confidence_interval(loglow, loghigh, seed)); + // Then see code for sample_normal_from_90_ci + double loglow = log(low); + double loghigh = log(high); + return exp(sample_normal_from_90_ci(loglow, loghigh, seed)); } double sample_gamma(double alpha, uint64_t* seed) @@ -201,7 +204,7 @@ double sample_mixture(double (*samplers[])(uint64_t*), double* weights, int n_di { // Sample from samples with frequency proportional to their weights. double sum_weights = array_sum(weights, n_dists); - double* cumsummed_normalized_weights = (double*)malloc(n_dists * sizeof(double)); + double* cumsummed_normalized_weights = (double*)malloc((size_t)n_dists * sizeof(double)); cumsummed_normalized_weights[0] = weights[0] / sum_weights; for (int i = 1; i < n_dists; i++) { cumsummed_normalized_weights[i] = cumsummed_normalized_weights[i - 1] + weights[i] / sum_weights; diff --git a/squiggle.c/squiggle_c/squiggle.h b/squiggle.c/squiggle_c/squiggle.h index 6f180389..100ea68d 100644 --- a/squiggle.c/squiggle_c/squiggle.h +++ b/squiggle.c/squiggle_c/squiggle.h @@ -15,6 +15,7 @@ double sample_unit_normal(uint64_t* seed); double sample_uniform(double start, double end, uint64_t* seed); double sample_normal(double mean, double sigma, uint64_t* seed); double sample_lognormal(double logmean, double logsigma, uint64_t* seed); +double sample_normal_from_90_ci(double low, double high, uint64_t* seed); double sample_to(double low, double high, uint64_t* seed); double sample_gamma(double alpha, uint64_t* seed); @@ -30,4 +31,7 @@ double array_std(double* array, int length); // Mixture function double sample_mixture(double (*samplers[])(uint64_t*), double* weights, int n_dists, uint64_t* seed); +// Macro to mute "unused variable" warning when -Wall -Wextra is enabled. Useful for nested functions +#define UNUSED(x) (void)(x) + #endif diff --git a/squiggle.c/squiggle_c/squiggle_more.c b/squiggle.c/squiggle_c/squiggle_more.c index 720fdccd..766a66f3 100644 --- a/squiggle.c/squiggle_c/squiggle_more.c +++ b/squiggle.c/squiggle_c/squiggle_more.c @@ -6,19 +6,27 @@ #include #include #include +#include // memcpy /* Parallel sampler */ +#define CACHE_LINE_SIZE 64 +typedef struct seed_cache_box_t { + uint64_t* seed; + char padding[CACHE_LINE_SIZE - sizeof(uint64_t*)]; +} seed_cache_box; +// This avoid false sharing. Dealing with this shaves ~2ms. + void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_threads, int n_samples) { - // Division terminology: - // a = b * quotient + reminder - // a = (a/b)*b + (a%b) + // Terms of the division: + // a = b * quotient + reminder + // a = b * (a/b) + (a%b) // dividend: a // divisor: b - // quotient = a / b - // reminder = a % b - // "divisor's multiple" := (a/b)*b + // quotient = a/b + // reminder = a%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,18 +34,17 @@ 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); + // uint64_t** seeds = malloc((size_t)n_threads * sizeof(uint64_t*)); + seed_cache_box* cache_box = (seed_cache_box*) malloc(sizeof(seed_cache_box) * (size_t)n_threads); srand(1); - for (uint64_t i = 0; i < n_threads; i++) { - seeds[i] = malloc(sizeof(uint64_t)); + for (int i = 0; i < n_threads; i++) { + cache_box[i].seed = malloc(sizeof(uint64_t*)); // Constraints: // - xorshift can't start with 0 // - the seeds should be reasonably separated and not correlated - *seeds[i] = (uint64_t)rand() * (UINT64_MAX / RAND_MAX); + *(cache_box[i].seed) = (uint64_t)rand() * (UINT64_MAX / RAND_MAX); // printf("#%ld: %lu\n",i, *seeds[i]); // Other initializations tried: @@ -47,33 +54,29 @@ void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_ } int i; -#pragma omp parallel private(i) +#pragma omp parallel private(i, quotient) { #pragma omp for 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]); + results[j] = sampler(cache_box[i].seed); } } } for (int j = divisor_multiple; j < n_samples; j++) { - results[j] = sampler(seeds[0]); + results[j] = sampler(cache_box[0].seed); // we can just reuse a seed, this isn't problematic because we are not doing multithreading } - for (uint64_t i = 0; i < n_threads; i++) { - free(seeds[i]); + for (int i = 0; i < n_threads; i++) { + free(cache_box[i].seed); } - free(seeds); + free(cache_box); } /* 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,10 +92,12 @@ 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 - int pivot = low + floor((high - low) / 2); + if (low > high || high >= length) { + printf("Invariant violated for function partition in %s (%d)", __FILE__, __LINE__); + exit(1); + } + // Note: the scratchpad/ folder in commit 578bfa27 has printfs sprinkled throughout + int pivot = low + (int)floor((high - low) / 2); double pivot_value = xs[pivot]; swp(pivot, high, xs); int gt = low; /* This pointer will iterate until finding an element which is greater than the pivot. Then it will move elements that are smaller before it--more specifically, it will move elements to its position and then increment. As a result all elements between gt and i will be greater than the pivot. */ @@ -109,15 +114,24 @@ static int partition(int low, int high, double xs[], int length) static double quickselect(int k, double xs[], int n) { // https://en.wikipedia.org/wiki/Quickselect + + double *ys = malloc((size_t)n * sizeof(double)); + memcpy(ys, xs, (size_t)n * sizeof(double)); + // ^: don't rearrange item order in the original array + int low = 0; int high = n - 1; for (;;) { if (low == high) { - return xs[low]; + double result = ys[low]; + free(ys); + return result; } - int pivot = partition(low, high, xs, n); + int pivot = partition(low, high, ys, n); if (pivot == k) { - return xs[pivot]; + double result = ys[pivot]; + free(ys); + return result; } else if (k < pivot) { high = pivot - 1; } else { @@ -129,8 +143,8 @@ static double quickselect(int k, double xs[], int n) ci array_get_ci(ci interval, double* xs, int n) { - int low_k = floor(interval.low * n); - int high_k = ceil(interval.high * n); + int low_k = (int)floor(interval.low * n); + int high_k = (int)ceil(interval.high * n); ci result = { .low = quickselect(low_k, xs, n), .high = quickselect(high_k, xs, n), @@ -144,10 +158,8 @@ 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); - }*/ + UNUSED(seed); // don't want to use it right now, but want to preserve ability to do so (e.g., remove parallelism from internals). Also nicer for consistency. + double* xs = malloc((size_t)n * sizeof(double)); sampler_parallel(sampler, xs, 16, n); ci result = array_get_ci(interval, xs, n); free(xs); @@ -159,9 +171,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 @@ -195,8 +204,8 @@ lognormal_params algebra_product_lognormals(lognormal_params a, lognormal_params lognormal_params convert_ci_to_lognormal_params(ci x) { - double loghigh = logf(x.high); - double loglow = logf(x.low); + double loghigh = log(x.high); + double loglow = log(x.low); double logmean = (loghigh + loglow) / 2.0; double logstd = (loghigh - loglow) / (2.0 * NORMAL90CONFIDENCE); lognormal_params result = { .logmean = logmean, .logstd = logstd }; @@ -221,21 +230,21 @@ 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); + printf("%s, @, in %s (%d)", error_msg, file, line); exit(1); } 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 +253,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 @@ -257,8 +266,9 @@ struct box inverse_cdf_double(double cdf(double), double p) // 1. Make sure that cdf(low) < p < cdf(high) int interval_found = 0; - while ((!interval_found) && (low > -FLT_MAX / 4) && (high < FLT_MAX / 4)) { - // ^ Using FLT_MIN and FLT_MAX is overkill + while ((!interval_found) && (low > -DBL_MAX / 4) && (high < DBL_MAX / 4)) { + // for floats, use FLT_MAX instead + // Note that this approach is overkill // but it's also the *correct* thing to do. int low_condition = (cdf(low) < p); @@ -299,7 +309,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 +320,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 @@ -323,15 +333,16 @@ struct box inverse_cdf_box(struct box cdf_box(double), double p) // 1. Make sure that cdf(low) < p < cdf(high) int interval_found = 0; - while ((!interval_found) && (low > -FLT_MAX / 4) && (high < FLT_MAX / 4)) { - // ^ Using FLT_MIN and FLT_MAX is overkill + while ((!interval_found) && (low > -DBL_MAX / 4) && (high < DBL_MAX / 4)) { + // for floats, use FLT_MAX instead + // Note that this approach 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 +372,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 +389,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 +400,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 +424,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("["); diff --git a/squiggle.c/squiggle_c/squiggle_more.h b/squiggle.c/squiggle_c/squiggle_more.h index 947b164d..2f88349a 100644 --- a/squiggle.c/squiggle_c/squiggle_more.h +++ b/squiggle.c/squiggle_c/squiggle_more.h @@ -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