diff --git a/C/samples.c b/C/samples.c index 6b779a6a..c32fa7ab 100644 --- a/C/samples.c +++ b/C/samples.c @@ -18,18 +18,6 @@ void array_print(float* array, int length) printf("\n"); } -void array_fill(float* array, int length, float item) -{ - int i; - #pragma omp private(i) - { - #pragma omp for - for (i = 0; i < length; i++) { - array[i] = item; - } - } -} - float array_sum(float* array, int length) { float output = 0.0; @@ -47,159 +35,12 @@ void array_cumsum(float* array_to_sum, float* array_cumsummed, int length) } } -float rand_float(float to, unsigned int* seed) -{ - return ((float)rand_r(seed) / (float)RAND_MAX) * to; - // See: - // rand() is not thread-safe, as it relies on (shared) hidden state. -} - -float ur_normal(unsigned int* seed) -{ - float u1 = rand_float(1.0, seed); - float u2 = rand_float(1.0, seed); - float z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2); - return z; -} - -inline float random_uniform(float from, float to, unsigned int* seed) -{ - return ((float) rand_r(seed) / (float)RAND_MAX) * (to - from) + from; -} - -inline float random_normal(float mean, float sigma, unsigned int* seed) -{ - return (mean + sigma * ur_normal(seed)); -} - -inline float random_lognormal(float logmean, float logsigma, unsigned int* seed) -{ - return expf(random_normal(logmean, logsigma, seed)); -} - -inline float random_to(float low, float high, unsigned int* seed) -{ - const float NORMAL95CONFIDENCE = 1.6448536269514722; - float loglow = logf(low); - float loghigh = logf(high); - float logmean = (loglow + loghigh) / 2; - float logsigma = (loghigh - loglow) / (2.0 * NORMAL95CONFIDENCE); - return random_lognormal(logmean, logsigma, seed); -} - +// Split array helpers int split_array_get_my_length(int index, int total_length, int n_threads) { return (total_length % n_threads > index ? total_length / n_threads + 1 : total_length / n_threads); } -//Old version, don't use it!! Optimized version is called mixture_f. This one is just for display -/* -void mixture(float* dists[], float* weights, int n_dists, float* results) -{ - float sum_weights = array_sum(weights, n_dists); - float* normalized_weights = malloc(n_dists * sizeof(float)); - for (int i = 0; i < n_dists; i++) { - normalized_weights[i] = weights[i] / sum_weights; - } - - float* cummulative_weights = malloc(n_dists * sizeof(float)); - array_cumsum(normalized_weights, cummulative_weights, n_dists); - - //create var holders - float p1, p2; - int index_found, index_counter, sample_index, i; - - #pragma omp parallel private(i, p1, p2, index_found, index_counter, sample_index) - { - #pragma omp for - for (i = 0; i < N; i++) { - p1 = random_uniform(0, 1); - p2 = random_uniform(0, 1); - - index_found = 0; - index_counter = 0; - - while ((index_found == 0) && (index_counter < n_dists)) { - if (p1 < cummulative_weights[index_counter]) { - index_found = 1; - } else { - index_counter++; - } - } - if (index_found != 0) { - sample_index = (int)(p2 * N); - results[i] = dists[index_counter][sample_index]; - } else - printf("This shouldn't be able to happen.\n"); - } - } - free(normalized_weights); - free(cummulative_weights); -} -*/ -void mixture_f(float (*samplers[])(unsigned int* ), float* weights, int n_dists, float** results, int n_threads) -{ - float sum_weights = array_sum(weights, n_dists); - float* normalized_weights = malloc(n_dists * sizeof(float)); - for (int i = 0; i < n_dists; i++) { - normalized_weights[i] = weights[i] / sum_weights; - } - - float* cummulative_weights = malloc(n_dists * sizeof(float)); - array_cumsum(normalized_weights, cummulative_weights, n_dists); - - //create var holders - float p1; - int sample_index, i, own_length; - unsigned int* seeds[n_threads]; - for(unsigned int i=0; i for why to use rand_r: + // rand() is not thread-safe, as it relies on (shared) hidden state. +} + +float ur_normal(unsigned int* seed) +{ + float u1 = rand_float(1.0, seed); + float u2 = rand_float(1.0, seed); + float z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2); + return z; +} + +inline float random_uniform(float from, float to, unsigned int* seed) +{ + return ((float) rand_r(seed) / (float)RAND_MAX) * (to - from) + from; +} + +inline float random_normal(float mean, float sigma, unsigned int* seed) +{ + return (mean + sigma * ur_normal(seed)); +} + +inline float random_lognormal(float logmean, float logsigma, unsigned int* seed) +{ + return expf(random_normal(logmean, logsigma, seed)); +} + +inline float random_to(float low, float high, unsigned int* seed) +{ + const float NORMAL95CONFIDENCE = 1.6448536269514722; + float loglow = logf(low); + float loghigh = logf(high); + float logmean = (loglow + loghigh) / 2; + float logsigma = (loghigh - loglow) / (2.0 * NORMAL95CONFIDENCE); + return random_lognormal(logmean, logsigma, seed); +} + +// Mixture function +void mixture(float (*samplers[])(unsigned int* ), float* weights, int n_dists, float** results, int n_threads) +{ + // You can see a simpler version of this function in the git history + // or in C-02-better-algorithm-one-thread/ + float sum_weights = array_sum(weights, n_dists); + float* normalized_weights = malloc(n_dists * sizeof(float)); + for (int i = 0; i < n_dists; i++) { + normalized_weights[i] = weights[i] / sum_weights; + } + + float* cummulative_weights = malloc(n_dists * sizeof(float)); + array_cumsum(normalized_weights, cummulative_weights, n_dists); + + //create var holders + float p1; + int sample_index, i, own_length; + unsigned int* seeds[n_threads]; + for(unsigned int i=0; i