squiggle.c/scratchpad/scratchpad.c

#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>
#include <float.h> // 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;
}