From 88b51331ea7b70f931602945a4c468793cac4f22 Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Sun, 16 Jul 2023 00:59:27 +0200 Subject: [PATCH] finish inverse! --- scratchpad/scratchpad | Bin 16704 -> 16912 bytes scratchpad/scratchpad.c | 107 ++++++++++++++++++++++++++++++---------- 2 files changed, 81 insertions(+), 26 deletions(-) diff --git a/scratchpad/scratchpad b/scratchpad/scratchpad index 37e8b9d8f84a1eee7909f988e0420dbe02dafae7..e5dacaac26b8e9e1f7d06204af388e18fdc34aff 100755 GIT binary patch delta 3513 zcmZ`+eQZG2_2Qp7(^fg{b@WZknxg*~ zv(|`NOP$vGX*+^l?;pCcIsLNRy?fTHtE0}@j_Gud)$ev7n3*0S`*~;w{jK#0>bDgL zGpNhv5lZ=cCmpe>NR1(f!F%kpU_FbTvP_d(63=f5Y;Oz9C-wF1+Y+tGc>R{<7f91f zi6m+7NN%V5ZN*}Fpn|Nfu3EmLzAR7{sHE+|ne?2kEM+AuRU=&_W8n!Ir6Wmi!bymg ztbh%}I;Jz(Le_w1k6BAL5%z@hblPXs*u;r}JHooqwYf41r|9pvvXm+1oQm}$k3 zK?7&DX0{9;eh9bI%XW{D)h}f2b{>;vIgY4=jA^rv(qd_Ps+kpfZ1jSrY3Xrs<{jSl zvC$z@hG94CMPHNmZMq2_sVJU}5aQhmJEy>~bz5%>c%|Gv(mu4@#Om0fzVa96=Pk0OXkFmhW-vAv+%yAWg0Q&NA}C}XpKur ziT+{CaZpr;cIL1`*mBv-o*Tj_yg#e}qB^{DDwI?IF%_LthjwI%{-HDRG*B%XU;s8MHoeU^{*AVbjft)#%>36_k04Yl1><7v zKQUH|l7qSpt$IJ7LWypu7(gA_gvsZ>0vsm8TD4?VjOYY$9Ud zQ}M0SVwA(bWQP7T`lx!!tX6fX2nN26a7aN^GvX`9n5BAFA3!|Z&mBMIGxKwI8%l81 zEQYI5$r-(MUrT95FW|RFafs?yJ37E7k%s?19Pr}}Q{DfWLdcGM=;qx7C#o$ld1Wcq z6IT3!0XFGppTnA9Qe517J~nhk27og_f+f;ODt^(A)N4Q2yMw@JMrw2`LWsUFN>=;l z&yQnK1KJ%v1AF`L!nsi~=AefFgb5Vzx5P8d8J#;3XXc&AcdNTYyZwYjHdJ-} zC9=M%>vm*ar0d@LSogdZz|U0A{hA#lnO_Na`rmP#%04FBQahUyElqMuTU$G6Y1_F- zo*n;|+?@D+Gvl(!i-|^^scPNPv~7En+}f6uw=}eD>1b(4HZ@}U>@k&XaXH!CByVYq z%MD5R%Zcr}H@Dh{!*t&VnwCP!sl9;M;rc!Nt}GJ@t&7>WFB`w(H*ki6UIHBe&48W< zE%*XON5L z!DK`Z(;vy-#Ha6qoS#Ypl7N*ko$)Y-SMg(`&Ugsx0sk4WCk>3ndynzRW)C}SSLRBG zEs@-U1J>2KHP%rfJGWp}t`y02te6r7U&Z<5Qyz!uHT;JErTMhqfg!P6$?|q~t|Nl= zh6f)El7s2^$El@gxnmvJnsJ<}MWy-uQr_Ar+2&V~%Url!#waBGsMi~GExb-b>#}%> zd9^Zaer-gnF`}7U<22o0>B+ueXk|Bjw{nj6HA6IS-Fz=1(g+=?^mwlu+8$%^@J#^i z7)#rpgY}t*i)jAiGr_S`CYPfaIF?9$7Kwu6m)8D`^D6{}-A=>dhpwfhzNxVx*+3Hd zSyqp8ivBz(i=g$}w{|paYig`7sV}7;1&cY>(VB?2ZF`YpH-hDJA6$PQgk?8hb{Bqt z;Qrh1TW~_?|NX+}1>q7cTCBQyYxH*tuXZ_|lD(TGHHPr?Q=vTm~aBQMG0bn;g<{~oo;Jqkrs33hQRknkBO zOwy3I%V7s6VU+sinHrnT#5uwldZ|bX+w3tysbdA4b|fMhJ|M#f2sHS#26&y_C6|tD zJl6hwn;)=3Co9w%o%B5EeXrTe_B zw=TC7AtGN7=SkT%OkR)hhh_M}Js5QJFa`^)HetEv!7F)w=+sNnrf2XA|f@%ZsVLNqgb)v_$jG!9;EqjL zfHA)860UgDLGLFVCBFWx=Xj$1_T=Sy)TLo&1R<6-s zY#)y|uO&p^v3cxgef#FIdwQon_V;#Uyk;%Bx3yRO!z#*lKTL1+P3Nai<(+j-g{)b_Vp54=4{5@AWN`GcP(o_VlCn+!~|jnu>m|V@<0A}(Z3)4%MpoDAaxda z2+Q$9Y^l5OAhN-88_g5AC`7d4bh6x7ip9w^U&w~Zq(+gmTtz{QW={s@YC_ZmSy5>U zbR|R&L64%!87?8A)&P8D`4K15@*}nUh)lXnj_rs~YHaY`( zv=~f>3ry9M{zVN?TxPD;Pv`ONb9xC5*8H>*-#+J3l4S^cwJ1%%D0&J>2TWr>egib0 zOD0(vu4@XFK#uxp7ruQ?AJ(8FGx6jhkKMXVPMR z+@vjaiYBg@T@WE}@@}!0!v78$l&}QO)Txfn2wbgGf{O=V$OW;yA~AnZ+-fn!Y|bbq z5vd=e9Y#>`BF06s;}r7=5&PbXHM`q z;%oFvq?AaDrajN&lRQ2rXN;dgKTBhA{7)=N_Z{rbeOA*~8e+7xdrxoAP!IpGCKAuK zZ{D!6JE}yLdZ=m&!_KBEm}-iq?>43tOuP^~O)qqn+YT^(5iZcwos0(jb5Qc6K*1V o{I^+^Xndrz3_maf1`-wxz<}Q9?P{>+&WCBe(Xp!uZs`^O0xx#1NdN!< 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; }