diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c index 7b2f4ae..20e5330 100644 --- a/scratchpad/scratchpad.c +++ b/scratchpad/scratchpad.c @@ -1,9 +1,9 @@ -#include // INT_MAX -#include -#include #include // FLT_MAX, FLT_MIN -#include +#include // INT_MAX #include // erf, sqrt +#include +#include +#include #include #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: // // @@ -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 , - // 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 , + // 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; + 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; - } - } + 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? - // 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? + // 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; }