commit c487bdfaf5a3f43a56dd62d26033770914073407
parent 607554f22b67531c8b0e75cf3611257b9bfbc18f
Author: NunoSempere <nuno.sempere@protonmail.com>
Date: Sun, 16 Jul 2023 16:30:38 +0200
formatting pass.
Diffstat:
1 file changed, 149 insertions(+), 144 deletions(-)
diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c
@@ -1,9 +1,9 @@
+#include <float.h> // FLT_MAX, FLT_MIN
#include <limits.h> // INT_MAX
+#include <math.h> // erf, sqrt
#include <stdint.h>
-#include <stdlib.h>
-#include <float.h> // FLT_MAX, FLT_MIN
#include <stdio.h>
-#include <math.h> // erf, sqrt
+#include <stdlib.h>
#include <time.h>
#define EXIT_ON_ERROR 0
@@ -12,7 +12,7 @@
struct box {
int empty;
float content;
- char * error_msg;
+ char* error_msg;
};
// Example cdf
@@ -50,9 +50,9 @@ struct box inverse_cdf(float cdf(float), float p)
{
// given a cdf: [-Inf, Inf] => [0,1]
// returns a box with either
- // x such that cdf(x) = p
- // or an error
- // if EXIT_ON_ERROR is set to 1, it exits instead of providing an error
+ // x such that cdf(x) = p
+ // or an error
+ // if EXIT_ON_ERROR is set to 1, it exits instead of providing an error
struct box result;
float low = -1.0;
@@ -75,17 +75,17 @@ struct box inverse_cdf(float cdf(float), float p)
}
}
- if (!interval_found) {
- if(EXIT_ON_ERROR){
- printf("Interval containing the target value not found, in function inverse_cdf, in %s (%d)", __FILE__, __LINE__);
- exit(1);
- }else{
- char error_msg[200];
- snprintf(error_msg, 200, "Interval containing the target value not found in function inverse_cdf, in %s (%d)", __FILE__, __LINE__);
- result.empty = 1;
- result.error_msg = error_msg;
- return result;
- }
+ if (!interval_found) {
+ if (EXIT_ON_ERROR) {
+ printf("Interval containing the target value not found, in function inverse_cdf, in %s (%d)", __FILE__, __LINE__);
+ exit(1);
+ } else {
+ char error_msg[200];
+ snprintf(error_msg, 200, "Interval containing the target value not found in function inverse_cdf, in %s (%d)", __FILE__, __LINE__);
+ result.empty = 1;
+ result.error_msg = error_msg;
+ return result;
+ }
} else {
int convergence_condition = 0;
@@ -93,9 +93,9 @@ struct box inverse_cdf(float cdf(float), float p)
while (!convergence_condition && (count < (INT_MAX / 2))) {
float mid = (high + low) / 2;
int mid_not_new = (mid == low) || (mid == high);
- // float width = high - low;
+ // float width = high - low;
if (mid_not_new) {
- // if ((width < 1e-8) || mid_not_new){
+ // if ((width < 1e-8) || mid_not_new){
convergence_condition = 1;
} else {
float mid_sign = cdf(mid) - p;
@@ -114,17 +114,17 @@ struct box inverse_cdf(float cdf(float), float p)
result.content = low;
result.empty = 0;
} else {
- if(EXIT_ON_ERROR){
- printf("Search process did not converge, in function inverse_cdf, in %s (%d)", __FILE__, __LINE__);
- exit(1);
- }else{
- char error_msg[200];
- snprintf(error_msg, 200, "Search process did not converge, in function inverse_cdf, in %s (%d)", __FILE__, __LINE__);
- result.empty = 1;
- result.error_msg = error_msg;
- return result;
- }
- }
+ if (EXIT_ON_ERROR) {
+ printf("Search process did not converge, in function inverse_cdf, in %s (%d)", __FILE__, __LINE__);
+ exit(1);
+ } else {
+ char error_msg[200];
+ snprintf(error_msg, 200, "Search process did not converge, in function inverse_cdf, in %s (%d)", __FILE__, __LINE__);
+ result.empty = 1;
+ result.error_msg = error_msg;
+ return result;
+ }
+ }
return result;
}
@@ -170,7 +170,7 @@ float sampler_normal_0_1(uint32_t* seed)
return z;
}
-// to do: add beta.
+// to do: add beta.
// for the cdf, use this incomplete beta function implementation, based on continuous fractions:
// <https://codeplea.com/incomplete-beta-function-c>
// <https://github.com/codeplea/incbeta>
@@ -178,12 +178,13 @@ float sampler_normal_0_1(uint32_t* seed)
#define STOP 1.0e-8
#define TINY 1.0e-30
-struct box incbeta(float a, float b, float x) {
- // Descended from <https://github.com/codeplea/incbeta/blob/master/incbeta.c>,
- // but modified to return a box struct and floats instead of doubles.
- // [x] to do: add attribution in README
- // Original code under this license:
- /*
+struct box incbeta(float a, float b, float x)
+{
+ // Descended from <https://github.com/codeplea/incbeta/blob/master/incbeta.c>,
+ // but modified to return a box struct and floats instead of doubles.
+ // [x] to do: add attribution in README
+ // Original code under this license:
+ /*
* zlib License
*
* Regularized Incomplete Beta Function
@@ -207,127 +208,131 @@ struct box incbeta(float a, float b, float x) {
* misrepresented as being the original software.
* 3. This notice may not be removed or altered from any source distribution.
*/
- struct box result;
-
- if (x < 0.0 || x > 1.0){
- if(EXIT_ON_ERROR){
- printf("x = %f, x out of bounds [0, 1], in function incbeta, in %s (%d)", __FILE__, __LINE__);
- exit(1);
- }else{
- char error_msg[200];
- snprintf(error_msg, 200, "x = %f, x out of bounds [0, 1], in function incbeta, in %s (%d)", x, __FILE__, __LINE__);
- result.empty = 1;
- result.error_msg = error_msg;
- return result;
- }
- }
+ struct box result;
+
+ if (x < 0.0 || x > 1.0) {
+ if (EXIT_ON_ERROR) {
+ printf("x = %f, x out of bounds [0, 1], in function incbeta, in %s (%d)", __FILE__, __LINE__);
+ exit(1);
+ } else {
+ char error_msg[200];
+ snprintf(error_msg, 200, "x = %f, x out of bounds [0, 1], in function incbeta, in %s (%d)", x, __FILE__, __LINE__);
+ result.empty = 1;
+ result.error_msg = error_msg;
+ return result;
+ }
+ }
/*The continued fraction converges nicely for x < (a+1)/(a+b+2)*/
- if (x > (a+1.0)/(a+b+2.0)) {
- struct box symmetric_incbeta = incbeta(b,a,1.0-x);
- if(symmetric_incbeta.empty){
- return symmetric_incbeta; // propagate error
- }else{
- result.empty = 0;
- result.content = 1-symmetric_incbeta.content;
- return result;
- }
+ if (x > (a + 1.0) / (a + b + 2.0)) {
+ struct box symmetric_incbeta = incbeta(b, a, 1.0 - x);
+ if (symmetric_incbeta.empty) {
+ return symmetric_incbeta; // propagate error
+ } else {
+ result.empty = 0;
+ result.content = 1 - symmetric_incbeta.content;
+ return result;
+ }
}
/*Find the first part before the continued fraction.*/
- const float lbeta_ab = lgamma(a)+lgamma(b)-lgamma(a+b);
- const float front = exp(log(x)*a+log(1.0-x)*b-lbeta_ab) / a;
+ const float lbeta_ab = lgamma(a) + lgamma(b) - lgamma(a + b);
+ const float front = exp(log(x) * a + log(1.0 - x) * b - lbeta_ab) / a;
/*Use Lentz's algorithm to evaluate the continued fraction.*/
float f = 1.0, c = 1.0, d = 0.0;
int i, m;
for (i = 0; i <= 200; ++i) {
- m = i/2;
+ m = i / 2;
float numerator;
if (i == 0) {
numerator = 1.0; /*First numerator is 1.0.*/
} else if (i % 2 == 0) {
- numerator = (m*(b-m)*x)/((a+2.0*m-1.0)*(a+2.0*m)); /*Even term.*/
+ numerator = (m * (b - m) * x) / ((a + 2.0 * m - 1.0) * (a + 2.0 * m)); /*Even term.*/
} else {
- numerator = -((a+m)*(a+b+m)*x)/((a+2.0*m)*(a+2.0*m+1)); /*Odd term.*/
+ numerator = -((a + m) * (a + b + m) * x) / ((a + 2.0 * m) * (a + 2.0 * m + 1)); /*Odd term.*/
}
/*Do an iteration of Lentz's algorithm.*/
d = 1.0 + numerator * d;
- if (fabs(d) < TINY) d = TINY;
+ if (fabs(d) < TINY)
+ d = TINY;
d = 1.0 / d;
c = 1.0 + numerator / c;
- if (fabs(c) < TINY) c = TINY;
+ if (fabs(c) < TINY)
+ c = TINY;
- const float cd = c*d;
+ const float cd = c * d;
f *= cd;
/*Check for stop.*/
- if (fabs(1.0-cd) < STOP) {
- result.content = front * (f-1.0);
- result.empty = 0;
- return result;
+ if (fabs(1.0 - cd) < STOP) {
+ result.content = front * (f - 1.0);
+ result.empty = 0;
+ return result;
}
}
- if(EXIT_ON_ERROR){
- printf("More loops needed, did not converge, in function incbeta, in %s (%d)", __FILE__, __LINE__);
- exit(1);
- }else{
- char error_msg[200];
- snprintf(error_msg, 200, "More loops needed, did not converge, in function incbeta, in %s (%d)", __FILE__, __LINE__);
- result.empty = 1;
- result.error_msg = error_msg;
- return result;
- }
+ if (EXIT_ON_ERROR) {
+ printf("More loops needed, did not converge, in function incbeta, in %s (%d)", __FILE__, __LINE__);
+ exit(1);
+ } else {
+ char error_msg[200];
+ snprintf(error_msg, 200, "More loops needed, did not converge, in function incbeta, in %s (%d)", __FILE__, __LINE__);
+ result.empty = 1;
+ result.error_msg = error_msg;
+ return result;
+ }
}
-struct box cdf_beta(float x){
- if(x < 0){
- struct box result = { .empty = 0, .content = 0};
- return result;
- } else if(x > 1){
- struct box result = { .empty = 0, .content = 1};
- return result;
- } else {
- float successes = 1, failures = (2023-1945);
- return incbeta(successes, failures, x);
- }
+struct box cdf_beta(float x)
+{
+ if (x < 0) {
+ struct box result = { .empty = 0, .content = 0 };
+ return result;
+ } else if (x > 1) {
+ struct box result = { .empty = 0, .content = 1 };
+ return result;
+ } else {
+ float successes = 1, failures = (2023 - 1945);
+ return incbeta(successes, failures, x);
+ }
}
-float cdf_dangerous_beta(float x){
- // So the thing is, this works
- // But it will propagate through the code
- // So it doesn't feel like a great architectural choice;
- // I prefer my choice of setting a variable which will determine whether to exit on failure or not.
- // Ok, so the proper thing to do would be to refactor inverse_cdf
- // but, I could also use a GOTO? <https://stackoverflow.com/questions/245742/examples-of-good-gotos-in-c-or-c>
- // Ok, alternatives are:
- // - Refactor inverse_cdf to take a box, take the small complexity + penalty. Add a helper
- // - Duplicate the code, have a refactored inverse_cdf as well as a normal cdf
- // - Do something hacky
- // a. dangerous beta, which exits
- // b. clever & hacky go-to statements
- // i. They actually look fun to implement
- // ii. But they would be hard for others to use.
- if(x < 0){
- return 0;
- } else if(x > 1){
- return 1;
- } else {
- float successes = 100, failures = 100;
- struct box result = incbeta(successes, failures, x);
- if(result.empty){
- printf("%s\n", result.error_msg);
- exit(1);
- return 1;
- }else{
- return result.content;
- }
- }
+float cdf_dangerous_beta(float x)
+{
+ // So the thing is, this works
+ // But it will propagate through the code
+ // So it doesn't feel like a great architectural choice;
+ // I prefer my choice of setting a variable which will determine whether to exit on failure or not.
+ // Ok, so the proper thing to do would be to refactor inverse_cdf
+ // but, I could also use a GOTO? <https://stackoverflow.com/questions/245742/examples-of-good-gotos-in-c-or-c>
+ // Ok, alternatives are:
+ // - Refactor inverse_cdf to take a box, take the small complexity + penalty. Add a helper
+ // - Duplicate the code, have a refactored inverse_cdf as well as a normal cdf
+ // - Do something hacky
+ // a. dangerous beta, which exits
+ // b. clever & hacky go-to statements
+ // i. They actually look fun to implement
+ // ii. But they would be hard for others to use.
+ if (x < 0) {
+ return 0;
+ } else if (x > 1) {
+ return 1;
+ } else {
+ float successes = 100, failures = 100;
+ struct box result = incbeta(successes, failures, x);
+ if (result.empty) {
+ printf("%s\n", result.error_msg);
+ exit(1);
+ return 1;
+ } else {
+ return result.content;
+ }
+ }
}
int main()
@@ -353,7 +358,7 @@ int main()
printf("Inverse of %s at %f is: %f\n", name_2, 0.5, result_2.content);
}
- // Get the inverse of a normal(0,1) cdf distribution
+ // Get the inverse of a normal(0,1) cdf distribution
struct box result_3 = inverse_cdf(cdf_normal_0_1, 0.5);
char* name_3 = "cdf_normal_0_1";
if (result_3.empty) {
@@ -363,15 +368,15 @@ int main()
printf("Inverse of %s at %f is: %f\n", name_3, 0.5, result_3.content);
}
- // Use the sampler on a normal(0,1)
+ // Use the sampler on a normal(0,1)
// set randomness seed
uint32_t* seed = malloc(sizeof(uint32_t));
*seed = 1000; // xorshift can't start with 0
int n = 100;
printf("\n\nGetting some samples from %s:\n", name_3);
- clock_t begin = clock();
- for (int i = 0; i < n; i++) {
+ clock_t begin = clock();
+ for (int i = 0; i < n; i++) {
struct box sample = sampler(cdf_normal_0_1, seed);
if (sample.empty) {
printf("Error in sampler function");
@@ -379,23 +384,23 @@ int main()
printf("%f\n", sample.content);
}
}
- clock_t end = clock();
- float time_spent = (float)(end - begin) / CLOCKS_PER_SEC;
- printf("Time spent: %f", time_spent);
+ clock_t end = clock();
+ float time_spent = (float)(end - begin) / CLOCKS_PER_SEC;
+ printf("Time spent: %f", time_spent);
- // Get some normal samples using the previous method.
- clock_t begin_2 = clock();
+ // 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);
+ float normal_sample = sampler_normal_0_1(seed);
+ printf("%f\n", normal_sample);
}
- clock_t end_2 = clock();
- float time_spent_2 = (float)(end_2 - begin_2) / CLOCKS_PER_SEC;
- printf("Time spent: %f", time_spent_2);
-
- // Get some beta samples
- clock_t begin_3 = clock();
+ clock_t end_2 = clock();
+ float time_spent_2 = (float)(end_2 - begin_2) / CLOCKS_PER_SEC;
+ printf("Time spent: %f", time_spent_2);
+
+ // Get some beta samples
+ clock_t begin_3 = clock();
printf("\n\nGetting some samples from box sampler_dangerous_beta\n");
for (int i = 0; i < n; i++) {
struct box sample = sampler(cdf_dangerous_beta, seed);
@@ -405,8 +410,8 @@ int main()
printf("%f\n", sample.content);
}
}
- clock_t end_3 = clock();
- float time_spent_3 = (float)(end_3 - begin_3) / CLOCKS_PER_SEC;
- printf("Time spent: %f", time_spent_3);
+ clock_t end_3 = clock();
+ float time_spent_3 = (float)(end_3 - begin_3) / CLOCKS_PER_SEC;
+ printf("Time spent: %f", time_spent_3);
return 0;
}
