From a4389e605fb477030dc3a0e52f2943b42359866a Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Sat, 23 Sep 2023 23:08:51 +0100 Subject: [PATCH] move to ci named struct. --- README.md | 3 ++ examples/08_nuclear_war/example | Bin 22672 -> 26760 bytes examples/08_nuclear_war/example.c | 4 +- examples/10_nuclear_recovery/example | Bin 22920 -> 22912 bytes examples/10_nuclear_recovery/example.c | 12 ++--- examples/12_algebra_and_conversion/example | Bin 0 -> 22496 bytes examples/12_algebra_and_conversion/example.c | 33 ++++++++++++ examples/12_algebra_and_conversion/makefile | 53 +++++++++++++++++++ squiggle.c | 14 ++--- squiggle.h | 10 ++-- 10 files changed, 109 insertions(+), 20 deletions(-) create mode 100755 examples/12_algebra_and_conversion/example create mode 100644 examples/12_algebra_and_conversion/example.c create mode 100644 examples/12_algebra_and_conversion/makefile diff --git a/README.md b/README.md index 607f1f1..d541ca7 100644 --- a/README.md +++ b/README.md @@ -348,6 +348,9 @@ It emits one warning about something I already took care of, so by default I've - [x] Provide example algebra - [x] Add conversion between 90% ci and parameters. - [ ] Use that conversion in conjuction with small algebra. + - [ ] Consider ergonomics of using ci instead of c_i + - use named struct instead + - demonstrate and document feeding a struct directly to a function; my_function((struct c_i){.low = 1, .high = 2}); - [ ] Test results - [x] Move to own file? Or signpost in file? => signposted in file. - [ ] Disambiguate sample_laplace--successes vs failures || successes vs total trials as two distinct and differently named functions diff --git a/examples/08_nuclear_war/example b/examples/08_nuclear_war/example index 8df8418cbb9c1e2bb2e6352217827f24f1636b91..b8009cff5d590bff641d3ac3fd7d5d837d83fe4b 100755 GIT binary patch delta 4907 zcmZu#4OEm>8vecy#|Dv^;m;Xj~^3D1QlgmjmQ$nSkm$|l~%CD zX1==}E6>q7dv;y7Iv8h*(hSqxL%U^YS+tq59?cPJHLP7q?e25GAN!iT=Y03M@AJI( zz4yEKy_av;vPbMZB5v_W1O8{SPXuaRf>Lx8Y^p5+#_^_?S14YkzJ{_{wRJT18xN^u zzN(ij|2Ymxsy(bCJt0O%X|A&7Dlezq;)hW_)v@my(r=_j-S?-k^n=tAdsoXJRlhd!8Mf`RY*b#<0k70WGFBTcIabjp_~I+ce(nS9RTAw}%>4p1)d{+EnUc z13$sQ2Q@yxJ46_Si-rzS27b)IYjaUlzsBDHi|>hEnhm_Rq>9vL;IV%J-CtO55KPxt zL3r4}^U3O^)4=n=^|Hml>#M}IZ3aGEcjn{oGzb%~34*Z8z(*SR9s_SN@O=i}YT)|~ zyq3I*I>7ljJbd=i3Z+(qjlIu@uN&`NAPDkf{p}JPls6qRA5=tElnr{6rROZ59NX2v zt%6{6mI=ZoYswE`(?jg!5IZi!jtsF)A@)_IzA^rHA$H*F5aLXT{Y8jRi5SzEX`35(esc`}hh zspm(=Slp>9LHKqH%k~p9Wh``|K{??s+_z%o1{q(irMtgr0=CKvS@?dExnONJM1meyezUha-`HOdM7E~ zI_~}sbXHbriAl<$)s8o$86HZG-!8q7Nqgh7rTa6v_63$D@GU#$bwp5&< z^NtLE z4Ur?T_$!R0(e~s_F_jJ_PnKHK=xlO|bpf>CqnG%KEiqN5QEZA^nvzC0r4+(%ZOU;e zB9*F~iQ$uZh$SXxDz!QbL>Im2Op#uhOb48|rMS2hRqe(q)*lE~T~$_JF>T3fT~y$j zBMrN#!?icbsd!Fw2MdvA`r<1Iw9wf*fAYRCsVb>AB7Ko4CXhd)0rlG%H-*P5vyXc= zVQc0i()oDWlu6R!cyf3)=9!^&MC09x#-r6_&AViMq8&%W*O!B<9LJ^HAqWLWtGmwP z5L7wnpl3t4T?rSt&_U%{Z%A!+`Z24o{L4vtpwnm}(5Ur^mk>zqZg98Bo4ztP--m!! zr&oB3|BZO_TOPjaWI=y*ccd9Ts=G~>;^=brI|YqU&L!gnDYwo~{a16|^h^FxI-k=l ziE&gr;#SVsCJtn9)-xy!~tZ%KSr%9B4^{4#xA|fw=x%PT+6m)2d&P{cqjD= zMr;q%$Qs2nDp2OeA#FeazQMs4^yvJ=G;m410qXzBM+t6J=X4q(PBMZc)Y9Q5>qEnG z7sgzA6wN{Tz$up;KB&CO!KPe&aK5(uZ-B3^y1R7QR6+19uITzC*;iZ9^__2#uj@*! z+#{O#3QSNVH~Y%+!=n`!XIXdiD5IRVeN5x zgyzGHHKX19A!A#B>w&wtJ;+!R4vXU`V>PhNA2HSs+yxxwejhU?Kpz9P0X@eU+YGD+ zZUc4#U*Y+@J&_!nayVBhDARdf3!V;h0Z!0o^R;67l< zFk^$j^}y4>G2j^IvB0tT=x9F8m>=jlgMQpT3xA+I!k8D|9%Caj8C#9UuJep_0-L{L zYyda_90q#6MnK>;paUNzj{h*059|c4r*nm#cq~aCPl~b|lt1}yG^gPid!Mpz%nW?M znDI-1e^0gJ832C~s2BXr5}v~2J;K-yPL^A2k4y6``BIm~=Chc&{xY6^=mS5~UsrUT zUg3Q_ZP*wdUyv`@p5_1X-eX$d^u{Lq&wZ zgAIW@bc_y9OB6q$k!fk-5Lt@y#LwxLA}8t&*-2dcZve zZiB*L!+?jtjUA`AiW0>W^hr^bI83S2XQEy){YLQ=ZJpjKo*}6?4Rr?UFR8M4R^Us< z>UfM+i))p|EbyIs8JRO2!tDChGyf`FS` zDf%xEeDNqn{}xWrZ0~Yh^7schYN18;h)i5mB#o0zi_m30w`5M(QVUnQz5?W^D@e{`` z=s(AdLKt>d=HJyJC5$hNKZJxCCL{$7{#T&YZ^ z(K&hgZ#Kmb6+dF)btZI7arot%@lz*UXN4|<6|eIGgsV#uCaReXAKF=Flf3!#T3O44 zV@P6`@^X}2em^$vrSe=!&Znd0jb7yf5(Fi`&}}+gqj4IBb77z;c?qo5_!=#9lv4

gSa9ZeAC?r}gSd<*@`|rnYb&fwx^}wTHpGkWDdx3ur}2@q*c} zbi7tw;k+TCa14|Zf!F8rzBpO6^>gje`81WEg8sZ@YP>K-wR=<$id7z=@D0O@ez(o0 zXXe=}YPD|t{A{8q7cQ_x;~&{M*z4BAlUiu9e1>UA)w>F3d z$(p|lCzls~A@h<-z6!f|Gq2;F+F1GcG~=aR1!2CTHX;JwKbv2pHMpBztgxFNfU4}X z_&X)yA4#mV$M46x_3p0@inRaMyL){LJ`{T` z3r~ef=_xe2D0hN>q2aV-*A|C+wJCXN4J@B(_4#ObZLw{(YUx*3J$+M~AFlm}rer>4 zFD}+C?|fSGlO@lm-LSL^NnSuBu>9IUHI!X<&9c?dnz~tH4ehR*W#)GeN)3Gjhb>E0 TE0GH7i^Kn-#Y&((w>ke0OeJs= delta 4691 zcmZu#dvH|M9X|IaE4v}tU9w3wo5yB>guEchA`HP8$w~qjNlbu55=BUOZl~cPjSwu; zML{8Irj*|}rOH&Ge~2B20k&uagebI=wrYb_hj9>sm9h{jF@OdUPQP;>cv(C%_k8E~ zeZTX#=lssOXLrjY+sVrYDE;kO6RS0nr`?T_NiKp7q1ti z71g}Y#K)QVh|Y&Zh2bV)z*J#`i61iY`dC!8S?7nr66<25po!O~Rh3qmc&uT_@MrBN z!E%RXY=emxgEi6?6E8Y9(pD31%#zTyn|PbyEc)MR65{U&jO{V;2_}A@iMN~hE)(xC z@!cj~4`5a85qvTpG5SQ6(z4##@T<1f_hn6COnJV$P3FDI#uL`#s>ml6fF9&Ig;S`| zc`Nh~V~(s!#;!WNKY-1NvNNOX`zN&KzkrMDAvrA#`L=tllsVv&?VIj&OV^Id7k*SVIDe<2-{*_q%_ zp*-(4=REnVGRCfLRXWRU4+_yD3z$6AmaM#FIa)@kNf}mnD4i2ODWeCIGNdwUO7iA| zVow(~C>zUHi=xtTzAavPseBELNXL(yM;4>)zgd=-(at1)C>@LCa9TD2BUiA`U?d1q z>9{iJ7*UyR73vzHp48Q+OsZ3;L7|3qb&*L;73$+ceMM92MRQH{+0LboydH>^*;JS8 z)5nmI3D2d9V?}plrskgP0t2@M4d?Zu+6R%3Z*nfk`%FElqr8=do}BtPhF}=)ppNp2 zPzQ0uh`N_dYN`v2<(!Z|)#S__X(XroT1uneq}-6tduW|ITR!HY*WD}RS`Rr=3*`b2 z9n0`IV%M%gWp$RCEO8!MocfL&NTp2A4*5@RI_xQsLvG=?EiGPdc9T7=UtW|#H`5B_ zDJcf;OQA>8`&^G=rz!qkW73-vCMHuLqgMVfiMC}F%5NmWzd+uYMEx1Bx?PC+=QiNj z8g$G{7=$HWVbDc8GIOOwI*~b2UhAR@nO?_SXp!fxiWzIN%yvn?5n0TbGeiJ zzNvD=No#zEGjOQDGd?>q9${uozPdEC9O0UgM`GkT@!{b)^^z1v&12>ve{akKHk@@B zeexgUw&gx6U$xQ3T$1P8C^c_=u@l;wM5q((C@pmqUv+G5!`29N4Z~56?NWsU)c;9! z#|3PHnN~WUx6X#c3KbHXtW=f%j@&+sZs&JZeS?EpYv?jcXlT%}c{dtTI_G7tQZ`<; z25ZrvBddZP65oh-uZiY6&icEnI}@y^QQc{|YN4A2$NXhbE@ol}si)3Q^(zE#xtE_$ z{e?lc=sDH-x5Tg`|LP-NeDXj=qhw7TS6P2^q^*J)1C1dmJ*I2| zhLnT137+9Rh;r~O=Uaj8z&*n5<2(Z~l=>Ctb+E1f;k+BT2iPb4zUG`k9|Eod=AGw! z6R;7u9k>Pfme7ISLchTIY2fw?Jj4f3I2GYM9v}PEOPuEcyZSk=7WNg+*8_vV9l##o z5n$;6=e@vo;Mc$*;E>=4IZwj3ICzcoW?*AIp2ap z5OLT8>;d)xPXVt1OA*BJSXKq-19ky}bn)Rl52mCTQ^qQXsb9iu6er*rI!y&rb3=p9TIJ?~%7|!ZoNAOghp)*tGgy1q#xFCt4iFMooE?;FdZZWv7D6Spc zr{EUqexl2l!Ifa`bAag5C;vUpZk;|{ZCAugCF01uPT=M>7lAD)+poK)ktE@qqWbBD>uGdYBDZ=`0}W!qJ|FCgrwT=_($E6R73ma6l$)0 zQ>vr5y1ZfjeBb1v2}N{nmW#S&XLw%SCdp%Tp|_wyjJT`8Eji?3#@+vJJ&6sBH#1r} z26sJUh-Nj0r7+{Re>d-D#;=aM`BbKetK&E`i*k&Z2@$?gDTbvnXB{1y{f0+ZMLdZl zOiAsYo|#jODf#o9Y4U_p8k|#X{4P=bkcs^wv3p1Knd)%wn~5#N*u67!6Rd>_BDuYoK51Cs)E8gfBB+)_ zPd(JRw?>R?_4~U!EQ#T1EW) z(yKJ<;dWRTUZZCl-JaL6c#Jq<+lwA~ni9tFm9?HqeKB%FCS6)w6la`t*jA&KOtR_W%t~k> zEPZ0T{Iq+?B&UA18E2eA7nhXS^uH@g%W2foNrqKXP78mslyce)OB`KD$`twvmNCB- zG-}x-W44_YstdEyN!+Q|FUzs&M4?I?`H8@-V%a3yo<}tQ1j-LiiaV(HkS0+`2tL97sbO$n zo+AhD;~Iw-uHc7`9QP4XDmKErt7@qr>Hq+!cOfcrxl=$ybv%p;?#9rb`|sK)Oi}!9k)yy)^a@M|9BgI zd-ng5Bh8QKY#`R%ckH~*~J0{m& zjK~drqq42cV~WQKrnV_0@-mSF5jm@ORIdAOR9-V-uf)ex{S}|c4eCbSao+lBpX}P1 zSoWAa?IGpPB6h;zgwnRX*poq*DFJ1-#{U`5ZlxM^oU$cWTKzq*62E5jbe+i69up;} zRJDTD+R_>b!E(**DgkSHAM_MeTT=2b@1X(sh8#t!0Y+h)}vG z8;Y5RV^HPq25nBBdo`rs@mNpuTuJY>vU66`~~qkBMn%0-rqIApb?kL~d3l z^sQpZ*#^z3ut!zwZwhumJbJ5OTX2ZjHo>}STHIP_Vcd41tyF=o#Ce*Z_VYs7AqF5E zyNzAM{w>%+0n^$NmXDW8m6WCQvnEysuOS2#h+Nm1+x_%=E zo)&2j%9cvFYDK?RSpEZOsVj$*q_3TLv%+VNiG}^n4 N*)RoA{Z**XcqZ}Pm*_lYxWYm^YOVoLqXyLVS#(eYG5yp@68 zTi%L5?_=*CZ|_=#zc247#k+UqJY_ZSyTUl=lJTlxX2Kl;y(Q1LXEn`Gj3G}*x#-CN zYtNdhuIZ>#v`<(^$e*;&rhU^>RXwaMymTQpXr@o$d@Lb7PgIq<5ee!ynWXF(WqTQE zG1}{ldctyNDF5+4xzq$QPucVnd{9uL&!yG=pfWBM$176?&|zVbA1ulTMdO0rlVjiN z%y+tRTu4{;7as#EguNyKa&hUkLHHU6d^RxSxUWUHDPa1JQLy0bj#{&Y?SYU&YB38J zzOz$=5}c)NIS%g>rFwyojdXkVL@gsz>&WRy&IoI_O!J&jx+b!VWbqT_W+D+1(#Fwt zj&eKQX?)-z{fZMLqJJi1m2p~^xV1Udg-ARXqwS_YO&jNPF5*C`C}nYV$f4meuD*=l zQPs=Uu}vDz;%XNzm5Mp1xaxLlIC44kb1b?h=k>TyY7=h8V7cEZ-{T_5(>1E-X&o{0{7!Ci4(Aev{i2k|)fs;1BU~-SxPU16xN4_vv1q3? zSOm};5IuENWk!ZP8s)zy(h!}+j}?8N@3UEOg0<3Fp3~0_;~lwczRE%L1VyQXt5ew;X8y70MwzPTqN+Ut z$s1ZnCRjzdri=s&1;zGyi6DjoL~d@ZKhP@Iyxa6vLtD+;Ev<*-W~?onNsRJ3xw*ck zxusE)6lZ7#Cd%)`2#-_nZiUle;df<0sq#Ag`GqT0fYTVRDlps*YK3oM+P(tAD1XXH IIK12T54*r`TL1t6 diff --git a/examples/10_nuclear_recovery/example.c b/examples/10_nuclear_recovery/example.c index 9b25e97..3891cd7 100644 --- a/examples/10_nuclear_recovery/example.c +++ b/examples/10_nuclear_recovery/example.c @@ -46,8 +46,8 @@ int main() // Before a first nuclear collapse printf("## Before the first nuclear collapse\n"); - struct c_i c_i_90_2023 = get_90_confidence_interval(yearly_probability_nuclear_collapse_2023, seed); - printf("90%% confidence interval: [%f, %f]\n", c_i_90_2023.low, c_i_90_2023.high); + ci ci_90_2023 = get_90_confidence_interval(yearly_probability_nuclear_collapse_2023, seed); + printf("90%% confidence interval: [%f, %f]\n", ci_90_2023.low, ci_90_2023.high); double* yearly_probability_nuclear_collapse_2023_samples = malloc(sizeof(double) * num_samples); for (int i = 0; i < num_samples; i++) { @@ -57,8 +57,8 @@ int main() // After the first nuclear collapse printf("\n## After the first nuclear collapse\n"); - struct c_i c_i_90_2070 = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_example, seed); - printf("90%% confidence interval: [%f, %f]\n", c_i_90_2070.low, c_i_90_2070.high); + ci ci_90_2070 = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_example, seed); + printf("90%% confidence interval: [%f, %f]\n", ci_90_2070.low, ci_90_2070.high); double* yearly_probability_nuclear_collapse_after_recovery_samples = malloc(sizeof(double) * num_samples); for (int i = 0; i < num_samples; i++) { @@ -68,8 +68,8 @@ int main() // After the first nuclear collapse (antiinductive) printf("\n## After the first nuclear collapse (antiinductive)\n"); - struct c_i c_i_90_antiinductive = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_antiinductive, seed); - printf("90%% confidence interval: [%f, %f]\n", c_i_90_antiinductive.low, c_i_90_antiinductive.high); + ci ci_90_antiinductive = get_90_confidence_interval(yearly_probability_nuclear_collapse_after_recovery_antiinductive, seed); + printf("90%% confidence interval: [%f, %f]\n", ci_90_antiinductive.low, ci_90_antiinductive.high); double* yearly_probability_nuclear_collapse_after_recovery_antiinductive_samples = malloc(sizeof(double) * num_samples); for (int i = 0; i < num_samples; i++) { diff --git a/examples/12_algebra_and_conversion/example b/examples/12_algebra_and_conversion/example new file mode 100755 index 0000000000000000000000000000000000000000..7fc289f310c1f0f22a2c7031f00dfb5fb532816a GIT binary patch literal 22496 zcmeHP4|J5(mH(2-K!9YvK|takbfkj@l}UntP-U5c1ioOD06`N73_~*cQn)!o|N)Ot?Y)9oH=?H=quCkY^^wZT7Ktg-~e`V9ezL6IUb`@8SYalRrO`J|b*&gXHo3|@5uyZuWOK9XH+c@X*cIW-kurJ0gl8?~DS19-j1)qe9 z;FD0=lWantxx${Y(1td;Z0uW6V)=>%Uo6ZO^b%4Zbglf*U;b74D-?WPIb2RGypar! zOG4S-a`2H|{?bdjd4=dNkspG>o)U`gm8_cT4RdB?)s)Yusjh2knb9(L&Wt&;9F0E5 zOr`?qMSW0AEqP!Sxg$t8VMgYQiK0&JJ8ym9Bj@vLI$wC}wne%5o@Y1L|EEpLpgxI) zF5;IiraF~rDiRGJ|Klie9Kl20H);tcx}HbFHE|mMm;gTxcv~F$pC_QNN`T*)0AG{< zUzPxWAOZe1;BENCUd2Gf8}APi;2%wZuStO4mHj7v1kJ_cnMds~i2^hNAly)%fbXMV<{cUd8RMto7Bo8~vUJzuV18;}UvX z>XjN_rPBCFgTI1(Q?){=_cbex)pgWendr`4-%wrWuW*-DJ*;?Js{QB$^2;73guB91 zU8B@`YHECCgs(HuKGNuGAmXwXkGrC}&QnwUC@Dyig@Uyn^ioYM&mnKp>GA>DS%1r!{@9M*rw~Nu~Ljy zj-pKSc-xczfax_#84&dVtEVXEP^XLd>j8ch^&Lz@>@&F7Y;jPc7u_5-;nAmvKDh0-0+F({I3wC6E#72K-ot;a`IWyc`p? zoHyXb7JxBS2|O)#Vn_^MZ@?S#R*L~oc}`xf20UAjMZoex^ScYM9cd5xlf!*l+uqb(rbBr4OrVFt)9wOLZ&rxD zidu$3$HL)oC#%zjU?_AbQm2i;P^ddnr;W-`=x32SZ3u=!KaA9ABQO+tGE%1vz)q)r=wp-@hwP8)!s&@GWVE&M~F@v@H5 z9{D>OZLY$`i~9VS`s|qcjF|egnEI5M`lOh8dQ9CMQ@@1c0HZ#i#ney7)IW@=ABm~I zW2|d?Fn(Nv3szQWYF)pqr#dUt%)Nvg?c88&pEAQ_Q5HrYU<`-g}OP# zPNm+LDH619-F8-8Jm+5z#j$vF?%YPy~Dg1EmRKfq0rDL z!xL2P)h}HK&!|?Afy7A0xQRrlnZL(pKjHegsMjxt!!({&gs;uT z`p|blDB-%CfD55XF4rW8VwznK@IGtosp!Q{qv*ruU3~W8^SV9;?fTpJf*lU>5qHck z5XR_c#@ItM0L?+IJ4CjCdDOZ?D7>oQLkhL+wdMrZ^oseEb2>k8Xq^@~uC<*QEG#PM zFjv6RT1WaBQnX`!F<8QrX5OnPpLU{Jf`ZmzUO`~qdH*B~-a=(e-7(>Q_2yRetw1(| zpr6j@XHO+C{jwGq)Ve=jpmm=&X-WIFLzn#%A>hM%xPa7fe+7#($xrjUbv`vKO?Rx) z+UECAQ44(NPt$hH-;auZ@M1WumqW$<=HH;6w3fzY#!_iP?N3GqDqv=xOUuVy-@PQXA7l z6g@*~F{3(q1E`_~#Kwcklws5{L7&4zK7=rWMz3lCl#W3olqRd&x5GvhgVg7%Fz6&7 z`h#{Tk1xQCX%DK~=^VzWd5`|+7vV4sVJV0ix4Js0W*$}@2hOU?-#e?0FSVyMZ-h9^ z3ZlCmbo$?e3BiAWL?MiG?I~#LMjU7K@_2!2ZO6$D2em^XZOkFk=^*TIFjPiPj`%d8 zy+Ur#>6b2q!vk#wlTp0DR_I&!13{3$UZY#48b@mOU=j_frLaaxH#F&i%fG~51FuK& zEBup*K3QV0LSbgk3?6_9`o5yCo{IHVFUvB_S7~i2)XBo=W56FUhrMzzVr`E4-EoTktuS*l!P!-#QcfYc=`ne)JalCQ4Es z^Vi?Per*|%SWjrb@15uNbAvWPs2lC)1eyL5@Lbt`p%+2M<0(@zd*aT!vvCvcC>=K2ZMN(Virz6 zfN6jzUT3bGjNF~m$0D01cls(q#L0M|fK?4yH(m;wpl<&`Y*?SVHsY27gc_xwe;g6l zpP%D?c?HA?{BonZ?c3mEX~;^!{Ugxo77R1WIX6av8^bA1u_Lx^YeJJjCfC1&dFg>G z#P0ShmdDCt6%>WjZ2GOAlI`Ijn2pv=eG)@hr_?ge(o&kYw0dUtf&d&CAb3lryu&5u+B1Q{ zQkn%3>eeo7_?RGI$wi3jY}zg6vdPN#8(3%Megi8U*4evZ4=|Z8;FuA*VD?cb`hAxe znpLe+j-Zj{2?{2Px+tT?+4zukUyE5-Q=zmz;Vj7=)((Q=Z-fS3*KY%DXe(2MHYR!{ zo&@k`_uMwn+w}6 zC;1BlLt5ZG&0M&q`*@0$bfFW9!+k;K3zW6<

dQPt%_QgKhnn+y0bro6)T zC-`bbXtMd=OS|GT!GWU&1JQUOhc+jp+Abt*p4xUYc@MkmYU@cx-e8i*Nc+)i5_NRi zs4}xJ^erZR%m3&4I?wG-@{bk@Lpf~mZhr%|={D$nkd%(#sge_<-4z=b=?@gL)9oL( zpHrVWg|R{xb<1fK0z+sz0UI7V30v6VaXv>N_y-5s`^^Qi|G=XGWC%!@n0D=D|DF29ljv%XlfDVu_ zmcQ9N4<*`YD&)G%_Xl1hO+p(nHHZCrg0<|Wp)#typ~YyYn*uZ-KTv|wa|c<_qloaP z29k>TPxf;?0~?4Iv< z3W?W}hOS`G+(179Gk0FUx@7?*aEtcdUv?3iA&=#Ee?&<%#a<8y>+8-iEhoS*SbGEe z06~H=!$RY&`8!5M(u=rpi@pzYwOvfw{4G~Uk^M5ik7hB<`~pd7UqH|Sm>Jo9e9V!IsA`wBh{Sh^I6;=Rp11KpO0?we3%G&HHroZ4i|E8Sa-;K0|oW z`qjPhgMz^4u5giA85rLQS^gJ557eS?(fso=#pn}0Wvx_Ghnes7g(&I2gLp%>t6=$x^uH)$+zE$HxfZN|lYYw}t!0K`<=W@w z1m*n*%^V6JT5=J-)-dL&wT2z-pavXnix3n#+0dyzy)XP0TZ9qAh-9?$4UiT5>Q*ZT z?c8!&-7*>lFdRp7`UI%(OELmN?*rt8*O3oiQ@ko}tUWEy1R1pvSurJ+?fGI3-Vloq3Wy1s-y-^rP*75hHF?+kXZcQsBIxl`Eqn_+=@5>M$yl490S|F@>J_84N&#UI~^H$ARIno5V;R zeGzxJn^+eVWL%d5HY1=c8#?J#!1O}gK{mrfJ6k(eLT2su27`{+&?0u}mLcm4bsN7Q z$vUcTJ;!uq=cP=`&?<5f*$tHX#z)DED5YDVIggMulz<}!X&t5Zddfppjp&Wy zfn_6>f{=YEC#7d2ZX1rDg9BSCI^#w>Z7q?z$j^qL ze~NF(VidJuLXJwI(AVL?K=Wzw6@+BFd(1JsUxZU>KZ{tTeR&B^BxsN8pjgKKbos$% zGfF+Z7!xf>v;$?C3};LFZ_qCi0G%>XdW-Bu;q2k~`(Mak`%Wr1#iGp9Cj`C*+jZ;Da+OY?f4+Bzq}9k6oBQD;6p{Wm;|Yer zO~Y$xWWV-MLC5?b-Hcrrw)uacZT@Gp&A+%x>zMHr&?^h3)BSO(zTh~5zWb8Ne{;@z zBF_~Bh6(~7=LSA zu)N_as`A?AWlx)CEA!P=RF`|}%DguG>eJh>(NkmdRoDuAm36)b{LnLVTE(5VX%*9L z)s40_yu210DB{NtNI8wIxvILX%2wrROe|CKlAIC3$0{qljZHOv>bcQZYZP7I;45z` z^EZCo=J!cm&_BOYog5A>aNY^^ZPOZUnbXRrD-K5%{u&=?s;;c8@jA+|c>VYp zYIR+u&5tqQhqIMlzYU{q^4jWrep`jFsji#}t7xh#^H=-oAb6v6+@S-FNwgZw_^d`rykz_J_m2e?`D0 zobSSdq#I24&~M?VqM`M_!r?-c>yL)REr6|noq+4XxeKuWSU603#6F~jK0xJzaQFmZ zHXz;Utp`klyZZobC{x_=QrtEz{C+^TL0(Jr3 zh~IQraTl8nSO}O4*at}eKVZe#g&)v46b=Ug6&&|nq&g0ovyoG6c<9ju*ax^5uoY{! z4(Plr3<#nB8uG6>;}h<;w*eJkA%a6gtu51TQZiajf$s;D7s=g;&u;kWTV(rX7HeB_sl{e`Ajx7~WT7npT_o#8e9l9!TL^)R za0l?22>(qd1TMn8k56_CZV# z(1rTyv6hTn#tH*vgS;QS7s2ZkdNAItCYQz5ma@dM-t-tcwn3bWG2J(=20RmC`_xY` zM+8p?>AusHYq4!bCz);LrIs$UsVm8{%baJ~ZLR@C-{$1CV=b8=a!Ku7^w)wpz8hsez9k5%W;i;J#Uat`0!<1YiX@`B7c{SfrV}*3LYd?|L3($X zV9=c@`IgMBWb?LBxfW+~E3-6=fVm4x=UY-~UZx?3Sdn8MB7UF6Rt9ep-VS_m44z`| zUf{ccuOfODs}EuFEM>t(e5Gg~M0*3Z(>z&fNzFykPxA0!ki{H}jr=&Vq4C&o*2N19 z7N@_4oL0yYsf>ELdR?`^RSR6Tz*P(UAF;rCS}1VY1f;D5U9qs9qm?uP*#dyF{GX## z2XDYnxY#BCU+e}^m;1--MP2^i7_DP;$@b4KhkaDKznxI1Xl) z>>t1CWtXJGEdU!Ix??DuhqAgZ28MG}R+oC=R)p1qq8+yytS-kP#Q)Vm(fa?XU7qC^pP64|F+$~^Gz&-){1soJm!H*ce$|e85aLJ-YU$m3yK8y3B}UX!_sGxpDC(MgIS}Q9Au^=VNi{&5<~d zrXQuq_>RVpR^&X0#-}PWKBDns6d701__2zN-)Q_eMaFqF9H8jYW<#G4-`+1b6SZKuXvPrpKvBt&AY&pzBX(>hdBs`4aCCYR0@c7Xl zqSLDE#&4XA7v&?u_s8Jx5O|sYC65b*B!iwJJ+k2>#j#Jbkh5MUP$8&|$r)+q=YWr= zf1+{yl*yT>?7UkDHh_3?*gY7szg%Y|`Rh#1NI$;A2uA9uGyKT7`hv+B>GuaP&y#7D zL`WO2@i5>ugcR$fT-R=4cq1h*&rg6~Cgfy`P+{j3U?~OOHd2i8uUJTI%E&luVRA+? z>|l7SQhz&d+bu}`S>lC%ZxQ%^2i^+#vGM#1Mn6*i4TMh~Ij{F6pg+d&R;5#HEZEr% z#106&bq4Q)-Q%J-7ULy5ciqYHRU9bO1>Pzabap<3w!0<0ofA4yr)w$0^Rn@KwV>$3#@l*He{5$xOu#V-UtC=T{yJd8t}lbk*g5AxjaT80-w4Pd&U zcM88eB--fQmE?2^KY9h8;yTWfpi~L^Lb1V;ex!qEEAmxFp`lY0Hwk(hB|2O@FW_uU z;QPPE@$761#h;4t%J&Cl@n+>YLEjZK-k(Z2LcSbVx4;hyJ8u$l=!}T;jGZrs1-%?J zJA(z!5r!Ac2Jm?Td@39r&u-TPPjM34@2!j;^Mw`;Tn}@g%m&^j%0{p(fqajU;}qkS zetcNqJH}OPl6YG@ zwe{>7N?C)y(cgq;iDk<0+eqC0T6Y<}kE9W=`Y89gD{FikJT>lezptUu?P+RJ%6!D` z_m(^6X3xoqPeiW^sdjrB8a$ia-a3E7CcGrXQ|on?H`Uf|0*j&M2BklW)bNzY?OwdX zb$^~a?}1!;8A>#5x#G?(dBAml{-Wp>_U;s*@J1eYo+cz{xhoX+l7eLmT?OuCix;oV zD{>dP78c}@RQ5uXvc@K+HR{Hfpg5iK%_j7=5_;9iFoxcB(#VYIuJqK_dcO38hilz5aQ>PZRH9#s<1TOXxvM;N<@Caw{AHjjudZ`9HG0d1gtDgE#->`S zQEhdLzp241s@&cid@V`~Lb|HD!aoN_qqpe9=u%VZ-O%8{%X8TC%3%Z1#E66pMXz{$ zPBgm+jEGqkJe73M%XSYxo<%hC$FmB$DDyTpy1fkzzJ`cecg<8hLI(^+n6y!hu-5CT zgB}~e#p9j5ipOnu)z7dx6uEdM&@h4mmtF)Ejd8p2CL%G_g!b@d1s)oI$!mznaM5Cf z4~Naa!HWns#t+#(`mI9I5~S%o^);R{FEp(;v`5F8R8Vm=ZmRWrHURn?IINN-X6Aau zfrr6fM`c};1Lm&Byx3$wZD^{-o5HHg1;vsW5Up1+P*Z<(ulz1&96y-a=22 z@jN-I;znD8x5h&ZqEuhwR~#%*94I;}eW+5Q(OafC{N5H+DPMxx$7Z|3TO~4NRXL=| z8kfOSCFhW3bmgh7hA_^KpimshB(+GrGV6b3NZV{|HL<7Xmpl)UP@Z3jW=1Q^1iyw4 z9TiFbxmdU9lF-I#!(Z~8NMh!oVB{|r{1VzFp@3NSqON?@jQsK(LPA>76M^K_OHOVi}NhHh^lv005C*fwa(;1zt%kvZoTLr(gvD9Cnl@1ihrjlQtyGR%m zgtC9hFUSA1;9n@{fY|lT7~zZ@k3tBB+;9GLB!KM@g6- zq2z_w`Lj>(OaDjCt+F|N7YR39(ti^5I~0gS^2_rr3FUhxC4FrF@5S)Tb1(_}Vi;ok z{}BB2{xFetBIjlDeU!2EvM%8#;8hL$@;pYH4 +#include +#include +#include + +int main() +{ + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 + + // Convert to 90% confidence interval form and back + lognormal_params ln1 = { .logmean = 1.0, .logstd = 3.0 }; + ci ln1_ci = convert_lognormal_params_to_ci(ln1); + printf("The 90%% confidence interval of Lognormal(%f, %f) is [%f, %f]\n", + ln1.logmean, ln1.logstd, + ln1_ci.low, ln1_ci.high); + lognormal_params ln1_ci_paramas = convert_ci_to_lognormal_params(ln1_ci); + printf("The lognormal which has 90%% confidence interval [%f, %f] is Lognormal(%f, %f)\n", + ln1_ci.low, ln1_ci.high, + ln1.logmean, ln1.logstd); + + lognormal_params ln2 = convert_ci_to_lognormal_params((ci){.low = 1, .high = 10}); + lognormal_params ln3 = convert_ci_to_lognormal_params((ci){.low = 5, .high = 50}); + + lognormal_params sln = algebra_product_lognormals(ln2, ln3); + ci sln_ci = convert_lognormal_params_to_ci(sln); + + printf("Result of some lognormal products: to(%f, %f)\n", sln_ci.low, sln_ci.high); + + free(seed); +} diff --git a/examples/12_algebra_and_conversion/makefile b/examples/12_algebra_and_conversion/makefile new file mode 100644 index 0000000..ef385f7 --- /dev/null +++ b/examples/12_algebra_and_conversion/makefile @@ -0,0 +1,53 @@ +# Interface: +# make +# make build +# make format +# make run + +# Compiler +CC=gcc +# CC=tcc # <= faster compilation + +# Main file +SRC=example.c ../../squiggle.c +OUTPUT=example + +## Dependencies +MATH=-lm + +## Flags +DEBUG= #'-g' +STANDARD=-std=c99 +WARNINGS=-Wall +OPTIMIZED=-O3 #-Ofast +# OPENMP=-fopenmp + +## Formatter +STYLE_BLUEPRINT=webkit +FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT) + +## make build +build: $(SRC) + $(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(MATH) -o $(OUTPUT) + +format: $(SRC) + $(FORMATTER) $(SRC) + +run: $(SRC) $(OUTPUT) + OMP_NUM_THREADS=1 ./$(OUTPUT) && echo + +time-linux: + @echo "Requires /bin/time, found on GNU/Linux systems" && echo + + @echo "Running 100x and taking avg time $(OUTPUT)" + @t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 1 thread: |" | sed 's|$$|ms|' && echo + +## Profiling + +profile-linux: + echo "Requires perf, which depends on the kernel version, and might be in linux-tools package or similar" + echo "Must be run as sudo" + $(CC) $(SRC) $(MATH) -o $(OUTPUT) + sudo perf record ./$(OUTPUT) + sudo perf report + rm perf.data diff --git a/squiggle.c b/squiggle.c index 9151ccd..65e808d 100644 --- a/squiggle.c +++ b/squiggle.c @@ -435,10 +435,10 @@ double sampler_danger(struct box cdf(double), uint64_t* seed) // Get confidence intervals, given a sampler -struct c_i { +typedef struct ci_t { float low; float high; -}; +} ci; int compare_doubles(const void* p, const void* q) { // https://wikiless.esmailelbob.xyz/wiki/Qsort?lang=en @@ -454,7 +454,7 @@ int compare_doubles(const void* p, const void* q) return 0; } -struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed) +ci get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed) { int n = 100 * 1000; double* samples_array = malloc(n * sizeof(double)); @@ -463,7 +463,7 @@ struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* se } qsort(samples_array, n, sizeof(double), compare_doubles); - struct c_i result = { + ci result = { .low = samples_array[5000], .high = samples_array[94999], }; @@ -505,7 +505,7 @@ lognormal_params algebra_product_lognormals(lognormal_params a, lognormal_params return result; } -lognormal_params convert_ci_to_lognormal_params(struct c_i x) +lognormal_params convert_ci_to_lognormal_params(ci x) { double loghigh = logf(x.high); double loglow = logf(x.low); @@ -515,12 +515,12 @@ lognormal_params convert_ci_to_lognormal_params(struct c_i x) return result; } -struct c_i convert_lognormal_params_to_ci(lognormal_params y) +ci convert_lognormal_params_to_ci(lognormal_params y) { double h = y.logstd * NORMAL95CONFIDENCE; double loghigh = y.logmean + h; double loglow = y.logmean - h; - struct c_i result = { .low=exp(loglow), .high=exp(loghigh)}; + ci result = { .low=exp(loglow), .high=exp(loghigh)}; return result; } diff --git a/squiggle.h b/squiggle.h index 8a93325..b97664c 100644 --- a/squiggle.h +++ b/squiggle.h @@ -52,11 +52,11 @@ struct box sampler_cdf_double(double cdf(double), uint64_t* seed); struct box sampler_cdf_box(struct box cdf(double), uint64_t* seed); // Get 90% confidence interval -struct c_i { +typedef struct ci_t { float low; float high; -}; -struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed); +} ci; +ci get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed); // small algebra manipulations @@ -73,9 +73,9 @@ typedef struct lognormal_params_t { lognormal_params algebra_product_lognormals(lognormal_params a, lognormal_params b); -lognormal_params convert_ci_to_lognormal_params(struct c_i x); +lognormal_params convert_ci_to_lognormal_params(ci x); -struct c_i convert_lognormal_params_to_ci(lognormal_params y); +ci convert_lognormal_params_to_ci(lognormal_params y); #endif