add example of parallelizing a min

examples/more/13_parallelize_min/example.c

@@ -0,0 +1,67 @@
+#include "../../../squiggle.h"
+#include "../../../squiggle_more.h"
+#include <stdio.h>
+#include <stdlib.h>
+
+int main()
+{
+    /* Question: can we parallelize this?
+    A = normal(5,2)
+    B = min(A)
+    B * 20
+    */
+
+    /* Option 1: parallelize taking from n samples */
+    // Question being asked: what is the distribution of sampling 1000 times and taking the min?
+    double sample_min_of_n(uint64_t* seed, int n){
+        double min = sample_normal(5, 2, seed);
+        for(int i=0; i<(n-2); i++){
+            double sample = sample_normal(5, 2, seed);
+            if(sample < min){
+                min = sample;
+            }
+        }
+        return min;
+    }
+    double sample_min_of_1000(uint64_t* seed) { 
+        return sample_min_of_n(seed, 1000);
+    } 
+
+    int n_samples = 1000000, n_threads = 16;
+    double* results = malloc(n_samples * sizeof(double));
+    parallel_sampler(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));
+    free(results);
+
+    /* Option 2: take the min from n samples cleverly using parallelism */
+    // Question being asked: can we take the min of n samples cleverly?
+    double sample_n_parallel(int n){
+
+        int n_threads = 16;
+        int quotient = n / 16;
+        int remainder = n % 16;
+
+        uint64_t seed = 1000;
+        double result_remainder = sample_min_of_n(&seed, remainder);
+        
+        double sample_min_of_quotient(uint64_t* seed) { 
+            return sample_min_of_n(seed, quotient);
+        }
+        double* results_quotient = malloc(quotient * sizeof(double));
+        parallel_sampler(sample_min_of_quotient, results_quotient, n_threads, quotient);
+
+        double min = results_quotient[0];
+        for(int i=1; i<quotient; i++){
+            if(min > results_quotient[i]){
+                min = results_quotient[i];
+            }
+        }
+        if(min > result_remainder){
+            min = results_remainder;
+        }
+        free(results_quotient);
+        return min;
+    }
+    printf("Minimum of 1M samples of normal(5,2): %f\n", sample_n_parallel(1000000));
+
+}

squiggle_more.c

@@ -25,21 +25,11 @@ typedef struct ci_t {
     float low;
     float high;
 } ci;
-int compare_doubles(const void* p, const void* q)
+typedef struct ci_searcher_t {
-{
+    double num;
-    // https://wikiless.esmailelbob.xyz/wiki/Qsort?lang=en
+    int remaining;
-    double x = *(const double*)p;
+} ci_searcher;
-    double y = *(const double*)q;
 
-    /* Avoid returning x - y, which can cause undefined behaviour
-       because of signed integer overflow. */
-    if (x < y)
-        return -1; // Return -1 if you want ascending, 1 if you want descending order.
-    else if (x > y)
-        return 1; // Return 1 if you want ascending, -1 if you want descending order.
-
-    return 0;
-}
 ci get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed)
 {
     int n = 100 * 1000;
@@ -47,7 +37,16 @@ ci get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed)
     for (int i = 0; i < n; i++) {
         samples_array[i] = sampler(seed);
     }
-    qsort(samples_array, n, sizeof(double), compare_doubles);
+    // 10% confidence interval: n/20, n - n/20
+    ci_searcher low = {.x = samples_array[0], .remaining = n/20) };
+    ci_searcher high = {.x = samples_array[0], .remaining = n-(n/20) };
+
+    // test with finding the lowest
+    for(int j=1; i<n; j++){
+        if(low.x > samples_array[i]){
+            low.x = samples_array[i];
+        }
+    }
 
     ci result = {
         .low = samples_array[5000],