add comparison with previous sampler.

New approach turns out to be ~50x slower.
Though it is much more general.
This commit is contained in:
NunoSempere 2023-07-16 13:02:11 +02:00
parent aca99a3386
commit 2acf129ef2
3 changed files with 36 additions and 13 deletions

View File

@ -10,7 +10,7 @@ CC=gcc # required for nested functions
# Main file # Main file
SRC=scratchpad.c SRC=scratchpad.c
OUTPUT=scratchpad OUTPUT=./scratchpad
## Dependencies ## Dependencies
MATH=-lm MATH=-lm

Binary file not shown.

View File

@ -4,6 +4,7 @@
#include <float.h> // FLT_MAX, FLT_MIN #include <float.h> // FLT_MAX, FLT_MIN
#include <stdio.h> #include <stdio.h>
#include <math.h> // erf, sqrt #include <math.h> // erf, sqrt
#include <time.h>
#define EXIT_ON_ERROR 0 #define EXIT_ON_ERROR 0
@ -41,9 +42,9 @@ float cdf_squared_0_1(float x)
float cdf_normal_0_1(float x) float cdf_normal_0_1(float x)
{ {
float mean = 0; // float mean = 0;
float std = 1; // float std = 1;
return 0.5 * (1 + erf((x - mean) / (std * sqrt(2)))); // erf from math.h return 0.5 * (1 + erf((x - 0) / (1 * sqrt(2)))); // erf from math.h
} }
// Inverse cdf // Inverse cdf
@ -94,13 +95,9 @@ struct box inverse_cdf(float cdf(float), float p)
while (!convergence_condition && (count < (INT_MAX / 2))) { while (!convergence_condition && (count < (INT_MAX / 2))) {
float mid = (high + low) / 2; float mid = (high + low) / 2;
int mid_not_new = (mid == low) || (mid == high); int mid_not_new = (mid == low) || (mid == high);
if (0) { // float width = high - low;
printf("while loop\n");
printf("low: %f, high: %f\n", low, high);
printf("mid: %f\n", mid);
}
if (mid_not_new) { if (mid_not_new) {
// if ((width < 1e-8) || mid_not_new){
convergence_condition = 1; convergence_condition = 1;
} else { } else {
float mid_sign = cdf(mid) - p; float mid_sign = cdf(mid) - p;
@ -166,6 +163,16 @@ struct box sampler(float cdf(float), uint32_t* seed)
return result; return result;
} }
// For comparison, raw sampler
const float PI = 3.14159265358979323846;
float sampler_normal_0_1(uint32_t* seed)
{
float u1 = rand_0_to_1(seed);
float u2 = rand_0_to_1(seed);
float z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2);
return z;
}
// to do: add beta. // to do: add beta.
// for the cdf, use this incomplete beta function implementation, based on continuous fractions: // for the cdf, use this incomplete beta function implementation, based on continuous fractions:
// <https://codeplea.com/incomplete-beta-function-c> // <https://codeplea.com/incomplete-beta-function-c>
@ -208,16 +215,32 @@ int main()
// set randomness seed // set randomness seed
uint32_t* seed = malloc(sizeof(uint32_t)); uint32_t* seed = malloc(sizeof(uint32_t));
*seed = 1000; // xorshift can't start with 0 *seed = 1000; // xorshift can't start with 0
int n = 1000000;
printf("\n\nGetting some samples from %s:\n", name_3); printf("\n\nGetting some samples from %s:\n", name_3);
int n = 100; clock_t begin = clock();
for (int i = 0; i < n; i++) { for (int i = 0; i < n; i++) {
struct box sample = sampler(cdf_normal_0_1, seed); struct box sample = sampler(cdf_normal_0_1, seed);
if (sample.empty) { if (sample.empty) {
printf("Error in sampler function"); printf("Error in sampler function");
} else { } else {
printf("%f\n", sample.content); // printf("%f\n", sample.content);
} }
} }
clock_t end = clock();
double time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
printf("Time spent: %f", time_spent);
// Get some normal samples using the previous method.
clock_t begin_2 = clock();
printf("\n\nGetting some samples from sampler_normal_0_1\n");
for (int i = 0; i < n; i++) {
float normal_sample = sampler_normal_0_1(seed);
// printf("%f\n", normal_sample);
}
clock_t end_2 = clock();
double time_spent_2 = (double)(end_2 - begin_2) / CLOCKS_PER_SEC;
printf("Time spent: %f", time_spent_2);
return 0; return 0;
} }