squiggle.c

Self-contained Monte Carlo estimation in C99
Log | Files | Refs | README

example.c (2228B)


      1 #include "../../../squiggle.h"
      2 #include "../../../squiggle_more.h"
      3 #include <stdio.h>
      4 #include <stdlib.h>
      5 
      6 int main()
      7 {
      8     /* Question: can we parallelize this?
      9     A = normal(5,2)
     10     B = min(A)
     11     B * 20
     12     */
     13 
     14     /* Option 1: parallelize taking from n samples */
     15     // Question being asked: what is the distribution of sampling 1000 times and taking the min?
     16     double sample_min_of_n(uint64_t * seed, int n)
     17     {
     18         double min = sample_normal(5, 2, seed);
     19         for (int i = 0; i < (n - 2); i++) {
     20             double sample = sample_normal(5, 2, seed);
     21             if (sample < min) {
     22                 min = sample;
     23             }
     24         }
     25         return min;
     26     }
     27     double sample_min_of_1000(uint64_t * seed)
     28     {
     29         return sample_min_of_n(seed, 1000);
     30     }
     31 
     32     int n_samples = 1000000, n_threads = 16;
     33     double* results = malloc((size_t)n_samples * sizeof(double));
     34     sampler_parallel(sample_min_of_1000, results, n_threads, n_samples);
     35     printf("Mean of the distribution of (taking the min of 1000 samples of a normal(5,2)): %f\n", array_mean(results, n_samples));
     36     free(results);
     37 
     38     /* Option 2: take the min from n samples cleverly using parallelism */
     39     // Question being asked: can we take the min of n samples cleverly?
     40     double sample_n_parallel(int n)
     41     {
     42 
     43         int n_threads = 16;
     44         int quotient = n / 16;
     45         int remainder = n % 16;
     46 
     47         uint64_t seed = 1000;
     48         double result_remainder = sample_min_of_n(&seed, remainder);
     49 
     50         double sample_min_of_quotient(uint64_t * seed)
     51         {
     52             return sample_min_of_n(seed, quotient);
     53         }
     54         double* results_quotient = malloc((size_t)quotient * sizeof(double));
     55         sampler_parallel(sample_min_of_quotient, results_quotient, n_threads, quotient);
     56 
     57         double min = results_quotient[0];
     58         for (int i = 1; i < quotient; i++) {
     59             if (min > results_quotient[i]) {
     60                 min = results_quotient[i];
     61             }
     62         }
     63         if (min > result_remainder) {
     64             min = result_remainder;
     65         }
     66         free(results_quotient);
     67         return min;
     68     }
     69     printf("Minimum of 1M samples of normal(5,2): %f\n", sample_n_parallel(1000000));
     70 }