forked from personal/squiggle.c
reformat & remake
This commit is contained in:
parent
186b10cddf
commit
ca1f81444e
|
@ -12,10 +12,10 @@ 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_few(uint64_t* seed) { return sample_to(1, 3, seed); }
|
||||
double sample_many(uint64_t* seed) { return sample_to(2, 10, seed); }
|
||||
double sample_0(uint64_t * seed) { return 0; }
|
||||
double sample_1(uint64_t * 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); }
|
||||
|
||||
int n_dists = 4;
|
||||
double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
|
||||
|
|
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.
|
@ -13,13 +13,14 @@ int main()
|
|||
|
||||
int n_samples = 1000 * 1000 * 1000;
|
||||
int n_threads = 16;
|
||||
double sampler(uint64_t* seed){
|
||||
double sampler(uint64_t * seed)
|
||||
{
|
||||
return sample_lognormal(0, 10, seed);
|
||||
}
|
||||
double* results = malloc(n_samples * sizeof(double));
|
||||
|
||||
sampler_parallel(sampler, results, n_threads, n_samples);
|
||||
double avg = array_sum(results, n_samples)/n_samples;
|
||||
double avg = array_sum(results, n_samples) / n_samples;
|
||||
printf("Average of 1B lognormal(0,10): %f", avg);
|
||||
|
||||
free(results);
|
||||
|
|
Binary file not shown.
|
@ -9,21 +9,22 @@ 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_few(uint64_t* seed) { return sample_to(1, 3, seed); }
|
||||
double sample_many(uint64_t* seed) { return sample_to(2, 10, seed); }
|
||||
double sample_0(uint64_t * seed) { return 0; }
|
||||
double sample_1(uint64_t * 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); }
|
||||
|
||||
int n_dists = 4;
|
||||
double weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 };
|
||||
double (*samplers[])(uint64_t*) = { sample_0, sample_1, sample_few, sample_many };
|
||||
double sampler_result(uint64_t* seed) {
|
||||
double sampler_result(uint64_t * seed)
|
||||
{
|
||||
return sample_mixture(samplers, weights, n_dists, seed);
|
||||
}
|
||||
}
|
||||
|
||||
int n_samples = 1000 * 1000, n_threads = 16;
|
||||
double* results = malloc(n_samples * sizeof(double));
|
||||
sampler_parallel(sampler_result, results, n_threads, n_samples);
|
||||
printf("Avg: %f\n", array_sum(results, n_samples)/n_samples);
|
||||
printf("Avg: %f\n", array_sum(results, n_samples) / n_samples);
|
||||
free(results);
|
||||
}
|
||||
|
|
Binary file not shown.
|
@ -13,19 +13,21 @@ int main()
|
|||
|
||||
/* 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 sample_min_of_n(uint64_t * seed, int n)
|
||||
{
|
||||
double min = sample_normal(5, 2, seed);
|
||||
for(int i=0; i<(n-2); i++){
|
||||
for (int i = 0; i < (n - 2); i++) {
|
||||
double sample = sample_normal(5, 2, seed);
|
||||
if(sample < min){
|
||||
if (sample < min) {
|
||||
min = sample;
|
||||
}
|
||||
}
|
||||
return min;
|
||||
}
|
||||
double sample_min_of_1000(uint64_t* seed) {
|
||||
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));
|
||||
|
@ -35,7 +37,8 @@ int main()
|
|||
|
||||
/* 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){
|
||||
double sample_n_parallel(int n)
|
||||
{
|
||||
|
||||
int n_threads = 16;
|
||||
int quotient = n / 16;
|
||||
|
@ -43,25 +46,25 @@ int main()
|
|||
|
||||
uint64_t seed = 1000;
|
||||
double result_remainder = sample_min_of_n(&seed, remainder);
|
||||
|
||||
double sample_min_of_quotient(uint64_t* seed) {
|
||||
|
||||
double sample_min_of_quotient(uint64_t * seed)
|
||||
{
|
||||
return sample_min_of_n(seed, quotient);
|
||||
}
|
||||
double* results_quotient = malloc(quotient * sizeof(double));
|
||||
sampler_parallel(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]){
|
||||
for (int i = 1; i < quotient; i++) {
|
||||
if (min > results_quotient[i]) {
|
||||
min = results_quotient[i];
|
||||
}
|
||||
}
|
||||
if(min > result_remainder){
|
||||
if (min > result_remainder) {
|
||||
min = result_remainder;
|
||||
}
|
||||
free(results_quotient);
|
||||
return min;
|
||||
}
|
||||
printf("Minimum of 1M samples of normal(5,2): %f\n", sample_n_parallel(1000000));
|
||||
|
||||
}
|
||||
|
|
Binary file not shown.
|
@ -8,7 +8,7 @@ int main()
|
|||
// set randomness seed
|
||||
uint64_t* seed = malloc(sizeof(uint64_t));
|
||||
*seed = 1000; // xorshift can't start with a seed of 0
|
||||
|
||||
|
||||
int n = 1000000;
|
||||
double* xs = malloc(sizeof(double) * n);
|
||||
for (int i = 0; i < n; i++) {
|
||||
|
|
4
makefile
4
makefile
|
@ -13,8 +13,8 @@ format-examples:
|
|||
cd examples/more && make format-all
|
||||
|
||||
format: squiggle.c squiggle.h
|
||||
$(FORMATTER) squiggle.c
|
||||
$(FORMATTER) squiggle.h
|
||||
$(FORMATTER) squiggle.c squiggle.h
|
||||
$(FORMATTER) squiggle_more.c squiggle_more.h
|
||||
|
||||
lint:
|
||||
clang-tidy squiggle.c -- -lm
|
||||
|
|
|
@ -1,15 +1,16 @@
|
|||
#include "squiggle.h"
|
||||
#include <float.h>
|
||||
#include <math.h>
|
||||
#include <limits.h>
|
||||
#include <math.h>
|
||||
#include <omp.h>
|
||||
#include <stdint.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include "squiggle.h"
|
||||
|
||||
/* Parallel sampler */
|
||||
void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_threads, int n_samples){
|
||||
if((n_samples % n_threads) != 0){
|
||||
void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_threads, int n_samples)
|
||||
{
|
||||
if ((n_samples % n_threads) != 0) {
|
||||
fprintf(stderr, "Number of samples isn't divisible by number of threads, aborting\n");
|
||||
exit(1);
|
||||
}
|
||||
|
@ -20,12 +21,12 @@ void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_
|
|||
}
|
||||
|
||||
int i;
|
||||
#pragma omp parallel private(i)
|
||||
#pragma omp parallel private(i)
|
||||
{
|
||||
#pragma omp for
|
||||
#pragma omp for
|
||||
for (i = 0; i < n_threads; i++) {
|
||||
int lower_bound = i * (n_samples / n_threads);
|
||||
int upper_bound = ((i+1) * (n_samples / n_threads)) - 1;
|
||||
int upper_bound = ((i + 1) * (n_samples / n_threads)) - 1;
|
||||
// printf("Lower bound: %d, upper bound: %d\n", lower_bound, upper_bound);
|
||||
for (int j = lower_bound; j < upper_bound; j++) {
|
||||
results[j] = sampler(seeds[i]);
|
||||
|
@ -95,7 +96,8 @@ static double quickselect(int k, double xs[], int n)
|
|||
}
|
||||
}
|
||||
|
||||
ci array_get_ci(ci interval, 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);
|
||||
|
@ -107,10 +109,11 @@ ci array_get_ci(ci interval, double* xs, int n){
|
|||
}
|
||||
ci array_get_90_ci(double xs[], int n)
|
||||
{
|
||||
return array_get_ci((ci) {.low = 0.05, .high = 0.95}, xs, n);
|
||||
return array_get_ci((ci) { .low = 0.05, .high = 0.95 }, xs, n);
|
||||
}
|
||||
|
||||
ci sampler_get_ci(ci interval, double (*sampler)(uint64_t*), int n, uint64_t* seed){
|
||||
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);
|
||||
|
@ -118,11 +121,10 @@ ci sampler_get_ci(ci interval, double (*sampler)(uint64_t*), int n, uint64_t* se
|
|||
ci result = array_get_ci(interval, xs, n);
|
||||
free(xs);
|
||||
return result;
|
||||
|
||||
}
|
||||
ci sampler_get_90_ci(double (*sampler)(uint64_t*), int n, uint64_t* seed)
|
||||
{
|
||||
return sampler_get_ci((ci) {.low = 0.05, .high = 0.95}, sampler, n, seed);
|
||||
return sampler_get_ci((ci) { .low = 0.05, .high = 0.95 }, sampler, n, seed);
|
||||
}
|
||||
|
||||
/* Algebra manipulations */
|
||||
|
@ -274,7 +276,7 @@ struct box inverse_cdf_double(double cdf(double), double p)
|
|||
}
|
||||
}
|
||||
|
||||
// Version #2:
|
||||
// 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)
|
||||
|
@ -372,11 +374,11 @@ double sampler_cdf_danger(struct box cdf(double), uint64_t* seed)
|
|||
{
|
||||
double p = sample_unit_uniform(seed);
|
||||
struct box result = inverse_cdf_box(cdf, p);
|
||||
if(result.empty){
|
||||
exit(1);
|
||||
}else{
|
||||
return result.content;
|
||||
}
|
||||
if (result.empty) {
|
||||
exit(1);
|
||||
} else {
|
||||
return result.content;
|
||||
}
|
||||
}
|
||||
|
||||
/* array print: potentially useful for debugging */
|
||||
|
@ -390,4 +392,3 @@ void array_print(double xs[], int n)
|
|||
printf("%f", xs[n - 1]);
|
||||
printf("]\n");
|
||||
}
|
||||
|
||||
|
|
|
@ -1,5 +1,5 @@
|
|||
#ifndef SQUIGGLE_C_EXTRA
|
||||
#define SQUIGGLE_C_EXTRA
|
||||
#define SQUIGGLE_C_EXTRA
|
||||
|
||||
/* Parallel sampling */
|
||||
void sampler_parallel(double (*sampler)(uint64_t* seed), double* results, int n_threads, int n_samples);
|
||||
|
@ -29,7 +29,7 @@ typedef struct lognormal_params_t {
|
|||
lognormal_params algebra_product_lognormals(lognormal_params a, lognormal_params b);
|
||||
|
||||
lognormal_params convert_ci_to_lognormal_params(ci x);
|
||||
ci convert_lognormal_params_to_ci(lognormal_params y);
|
||||
ci convert_lognormal_params_to_ci(lognormal_params y);
|
||||
|
||||
/* Error handling */
|
||||
struct box {
|
||||
|
|
Loading…
Reference in New Issue
Block a user