squiggle.c

commit 88b51331ea7b70f931602945a4c468793cac4f22
parent d4f61427341c22c17f025c340e2d57df17bcf7a2
Author: NunoSempere <nuno.sempere@protonmail.com>
Date:   Sun, 16 Jul 2023 00:59:27 +0200

finish inverse!

Diffstat:
M
scratchpad/scratchpad
 | 
0

M
scratchpad/scratchpad.c
 | 
107
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------

2 files changed, 81 insertions(+), 26 deletions(-)
diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad
Binary files differ.
diff --git a/scratchpad/scratchpad.c b/scratchpad/scratchpad.c
@@ -2,6 +2,10 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <float.h> // FLT_MAX, FLT_MIN
+#include <limits.h> // 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;
 }