diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c index f42803f..90e3565 100644 --- a/scratchpad/scratchpad.c +++ b/scratchpad/scratchpad.c @@ -1,7 +1,97 @@ #include #include #include +#include // FLT_MAX, FLT_MIN +// Errors +// https://mccue.dev/pages/7-27-22-c-errors +// Huh, first time I've really needed a struct +// https://en.wikipedia.org/wiki/Struct_(C_programming_language) +// https://www.cs.yale.edu/homes/aspnes/pinewiki/C(2f)Structs.html +// https://stackoverflow.com/questions/9653072/return-a-struct-from-a-function-in-c +// options: +// - exit +// - pass structs +struct box { + int empty; + float content; +}; + +// Example cdf +float cdf_uniform_0_1(float x){ + if(x < 0){ + return 0; + } else if (x > 1){ + return 1; + } else { + return x; + } +} + +// Inverse cdf +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; + + // 1. Make sure that cdf(start) < p < cdf(end) + // [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)){ + // ^ 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 ){ + interval_found = 1; + }else if(!start_condition){ + start = start * 2; + }else if (!end_condition){ + end = end * 2 ; + } + } + + 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); + if(mid_not_new){ + convergence_condition = 1; + }else{ + float mid_sign = *cdf(mid) - p; + if(mid_sign < 0){ + start = mid; + } else if (mid_sign > 0){ + end = mid; + } else { + result.content = mid; + result.empty = 0; + return result; + } + + } + } + return result; + } + +} + +// sampler based on inverse cdf + +// to-do: integrals + +// main with an example int main(){ printf("Hello world"); + return 0; }