update squiggle.c with speedup after avoiding conflicting cache hits

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  |

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_0(uint64_t * seed) { UNUSED(seed); return 0; }
-    double sample_1(uint64_t * seed) { return 1; }
+    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;
+    int n_samples = 1000 * 1000, n_threads = 16;
-    double* results = malloc(n_samples * sizeof(double));
+    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);

squiggle.c/squiggle_c/squiggle.c

@@ -3,12 +3,14 @@
 #include <stdint.h>
 #include <stdlib.h>
 
-// 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
+// Pseudo Random number generators
-static uint64_t xorshift32(uint32_t* seed)
+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)
     //   <https://prng.di.unimi.it/>
     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: <https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>
     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
+    // Then see code for sample_normal_from_90_ci
-    double loglow = logf(low);
+    double loglow = log(low);
-    double loghigh = logf(high);
+    double loghigh = log(high);
-    return exp(sample_normal_from_90_confidence_interval(loglow, loghigh, seed));
+    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;

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

squiggle.c/squiggle_c/squiggle_more.c

@@ -6,19 +6,27 @@
 #include <stdint.h>
 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h> // 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:
+    // Terms of the division:
-    // a =  b * quotient + reminder
+    // a = b * quotient + reminder
-    // a = (a/b)*b + (a%b)
+    // a = b * (a/b)    + (a%b)
     // dividend: a
     // divisor: b
-    // quotient = a / b
+    // quotient = a/b
-    // reminder = 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,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*));
+    // uint64_t** seeds = malloc((size_t)n_threads * sizeof(uint64_t*));
-    // printf("UINT64_MAX: %lu\n", UINT64_MAX);
+    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++) {
+    for (int i = 0; i < n_threads; i++) {
-        seeds[i] = malloc(sizeof(uint64_t));
+        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++) {
+    for (int i = 0; i < n_threads; i++) {
-        free(seeds[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:
+    if (low > high || high >= length) {
-    // - see the note after gt variable definition
+        printf("Invariant violated for function partition in %s (%d)", __FILE__, __LINE__);
-    // - go to commit 578bfa27 and the scratchpad/ folder in it, which has printfs sprinkled throughout
+        exit(1);
-    int pivot = low + floor((high - low) / 2);
+    }
+    // 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 low_k = (int)floor(interval.low * n);
-    int high_k = ceil(interval.high * 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));
+    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.
-    /*for (int i = 0; i < n; i++) {
+    double* xs = malloc((size_t)n * sizeof(double));
-        xs[i] = sampler(seed);
-    }*/
     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 loghigh = log(x.high);
-    double loglow = logf(x.low);
+    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)) {
+    while ((!interval_found) && (low > -DBL_MAX / 4) && (high < DBL_MAX / 4)) {
-        // ^ Using FLT_MIN and FLT_MAX is overkill
+        // 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)) {
+    while ((!interval_found) && (low > -DBL_MAX / 4) && (high < DBL_MAX / 4)) {
-        // ^ Using FLT_MIN and FLT_MAX is overkill
+        // 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("[");

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);
+box inverse_cdf_double(double cdf(double), double p);
-struct box inverse_cdf_box(struct box cdf_box(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);
+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_box(box cdf(double), uint64_t* seed);
 
 #endif