diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad index 37e8b9d..e5dacaa 100755 Binary files a/scratchpad/scratchpad and b/scratchpad/scratchpad differ diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c index 90e3565..9ed24e9 100644 --- a/scratchpad/scratchpad.c +++ b/scratchpad/scratchpad.c @@ -2,6 +2,10 @@ #include #include #include // FLT_MAX, FLT_MIN +#include // INT_MAX +#define VERBOSE 1 + +// to do: reuse more informative printing from build-your-own-lisp // Errors // https://mccue.dev/pages/7-27-22-c-errors @@ -28,70 +32,121 @@ float cdf_uniform_0_1(float x){ } } +float cdf_squared_0_1(float x){ + if(x < 0){ + return 0; + } else if (x > 1){ + return 1; + } else { + return x*x; + } +} + // Inverse cdf -struct box inverse_cdf(float (*cdf(float)), float p){ +struct box inverse_cdf(float cdf(float), float p){ // given a cdf: [-Inf, Inf] => [0,1] // returns x such that cdf(x) = p // to do: add bounds, add error checking // [x] maybe return a struct or smth. struct box result; - float start = -1.0; - float end = 1.0; + float low = -1.0; + float high = 1.0; - // 1. Make sure that cdf(start) < p < cdf(end) + // 1. Make sure that cdf(low) < p < cdf(high) // [x] to do: do smth with float min and float max? int interval_found = 0; - while((!interval_found) && (start > FLT_MIN/4) && (end < FLT_MAX/4)){ + while((!interval_found) && (low > -FLT_MAX/4) && (high < FLT_MAX/4)){ // ^ Using FLT_MIN and FLT_MAX is overkill // but it's also the *correct* thing to do. - int start_condition = (*cdf(start) < p); - int end_condition = (p < *cdf(end)); - if( start_condition && end_condition ){ + int low_condition = (cdf(low) < p); + int high_condition = (p < cdf(high)); + if( low_condition && high_condition ){ interval_found = 1; - }else if(!start_condition){ - start = start * 2; - }else if (!end_condition){ - end = end * 2 ; + }else if(!low_condition){ + low = low * 2; + }else if (!high_condition){ + high = high * 2 ; } } + if(0){ + printf("FLT_MIN = %f, FLT_MAX = %f, INT_MAX = %d\n", -FLT_MAX, FLT_MAX, INT_MAX); + printf("low: %f, high: %f\n", low, high); + printf("interval_found? %d\n", interval_found); + int while_condition = (!interval_found) && (low > FLT_MIN/4) && (high < FLT_MAX/4); + printf("while condition: %i\n", while_condition); + } if(!interval_found){ result.empty = 1; return result; } else{ + int convergence_condition = 0; - while(!convergence_condition){ - float mid = (end - start)/2; - int mid_not_new = (mid == start) || (mid == end); + int count = 0; + while(!convergence_condition && (count < (INT_MAX/2) )){ + if(VERBOSE){ + printf("while loop\n"); + } + float mid = (high + low)/2; + int mid_not_new = (mid == low) || (mid == high); + if(VERBOSE){ + printf("low: %f, high: %f\n", low, high); + printf("mid: %f\n", mid); + } + if(mid_not_new){ convergence_condition = 1; - }else{ - float mid_sign = *cdf(mid) - p; + } else{ + float mid_sign = cdf(mid) - p; if(mid_sign < 0){ - start = mid; + low = mid; } else if (mid_sign > 0){ - end = mid; - } else { - result.content = mid; - result.empty = 0; - return result; + high = mid; + } else if (mid_sign == 0){ + low = mid; + high = mid; } - } + } + + if(convergence_condition){ + result.content = low; + result.empty = 0; + } else{ + result.empty = 1; + } + return result; } } // sampler based on inverse cdf - // to-do: integrals // main with an example int main(){ - printf("Hello world"); + + // Uniform: + struct box result = inverse_cdf(cdf_uniform_0_1, 0.5); + if(result.empty){ + printf("Inverse not calculated\n"); + exit(1); + }else{ + printf("Inverse of the cdf at %f is: %f\n", 0.5, result.content); + } + + // Squared cdf + struct box result2 = inverse_cdf(cdf_squared_0_1, 0.5); + if(result2.empty){ + printf("Inverse not calculated\n"); + exit(1); + }else{ + printf("Inverse of the cdf at %f is: %f\n", 0.5, result2.content); + } + return 0; }