From 08eb790a6de911bb6d5287fe7cc71359de776566 Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Sun, 23 Jul 2023 10:09:03 +0200 Subject: [PATCH] add beta and gamma example/tests --- examples/06_gamma_beta/example | Bin 0 -> 22160 bytes examples/06_gamma_beta/example.c | 50 +++++++++++++++++++++++++++++ examples/06_gamma_beta/makefile | 53 +++++++++++++++++++++++++++++++ 3 files changed, 103 insertions(+) create mode 100755 examples/06_gamma_beta/example create mode 100644 examples/06_gamma_beta/example.c create mode 100644 examples/06_gamma_beta/makefile diff --git a/examples/06_gamma_beta/example b/examples/06_gamma_beta/example new file mode 100755 index 0000000000000000000000000000000000000000..4ed2a6e567ad69da239f55f864c61f82a6f23e26 GIT binary patch literal 22160 zcmeHPe|%Kcm4C@hAV@NA(8NXsnbOdX7A6S-B`C{G2;qeWiIS`#M8+XAA^9{vU@{|M zTZ)DNBSh(EyOrH`yV%`Ycem}w^|ReVi)E7V+g1Z&i$&KMvCe=94oVS`x8HMrOx_GL z>UOuE&;F77$;-L-+;h%7_ndQoy!+m~Z_Tf&$jr!4WHKp#rBG@64IJYX6$5vR41iOa zrR3uOMCDdx6v}3Xw38jVvi7#E+FoxGis91TK z__8P~%E{38IZh{`tdF->%Tr6q3#ynjx6)q+ne;c}AUs~oLJ zDBHUSeA2{%G@O|{MXaq$=IpnY)ONk3A_zga;gU+m7G1mPf8}H zGUBm~!f&SHuB=()_WJz3l`WkCpTB0&yw;9(UyWx)t50#eSGIMuyE_9Of57eLq$vqo zJ60;4tNno{_FqA3D&X>ILjr22LPkK5h2*5ht!Y4^0YJVZK@Y@uVD2W>U^eLiJDRb|;c z_cX^e$Bf~8>2UTd@$B???heOvm?jgWhQ^H_ zgEzNe#Tcy&30&yg-^DP?Qr;Ei&8%!zK17)k@vjDy^GjY!OD>H}GlUY3%QMC01#Tkw znTmxPMUv})oPSr*B!Q$Y;!xs?1wNYw29m_f@`=S9PjdpzBQgyd@YhKsV8no@u`JVh z1Ad%D0$wuU#a4jTHVgb{lth;pew6`loWIr@@N#Y>h7AThTbG#>McHD&^ZAZfY&YPk zZJBl$@X|(9x6^Ow9A0E#Icy{4S3anKV-nm7@_I`13uNb zH;5R?j#yyC0wWd}vA~E0{^Az6px$~}>v|_!+nV*tTx>(S!+}hFKqHpNBzMYk6BM^%ojhATy5R3N4%e3&vqW>N*(}pD$eKKCA4L~gV z?Rc3M{#f+kc$pUbSafx~ObdN1>WPq9yS%E&Q?Q=i_Bs z@MF<&vW(uo_E*%~a)tF5qB5<7++TY|*CvF4@70xxO=faVV z@KM*9jNrMw%9uGy#0Toi(0(n{U)8gdFp~a)YR$^ImGNm^{pLs&_%5j9=3S{QeKj~# z6*~88YB5p^7s&#u%L{WNkHCxIt}4vgOYx=k%&427D7z@OE~xckT!$l9uj)G3B29pM z$?K6;$}yFWA(x~Q^QSd;1LTli0YLPJ>2(}271!$zQkBSa5Go5ouhWROM7=KIkSb9x z(ueeV;bW=SapcfeWDmiveFY`qW%s)pT=(}6n>Oy>lCzbeqsz6>kk)nTe05D#Pu5e| zf@?jwMGVZ*AYIR&beE#Mzg_E@Q;Ce$lT}8b^im)nL;32R$YUsv(FfIu8(8=lmxMk=<8xJ>rFm2eN*Oa>o7Ojkizy(KVSF%l_8>c4zL`ZiS|%OPdDbK9 z1X&=-w*OPxW`eh3@hy!o?qDd1%FfjZvG4R?{2hBu7BaCPw1a3-oM8|#ysByzR%Xb=&wG+=ov+Q z+@(%*GXF2pLLb40R^}~p~5avSF-Wxdp)K(Lad4-`q)qVutL*6X9<8C5(YfKws<3g@^#_vgnRWFwvX%cjgytyf49+;wCCqi?V)zlYI4& zdaM)FJ}!)$`a?0bA=iesVwm?itUt``IR&+dQ;If>Zd;S;7$7xb>OF_m%}+rJaa@GV zn7X`I^azFpt)2X{)dXDXxrq!!fbW(Wvq>`W>z6%TT0^!#=OXRkAWYdQeVSn8Uz5WZpWyC zMV3N|x-9>VuvO6|wQB}uWK8=vGAZLgPhy%hwf#_JHe|70VrsK8N%bet`gl_NBuydp zXvW?L^Pw`uG(tTC!_o=gL}MWnTG@pRlMUXJ;A<$fMvhIFSgQ?~ zH;bNM(}~9VSz`r}yA#iM7vvW^cBTq7s)~TgjRFx1_00Acu4pA8P7_NM( z2X22%t%G1pH#X!D1Vw08-E<2UTV^-}3e1ZLEMC@67N40qg+7kYA32f3SZ0`}IoOauMOcmIkY_m>ou*NCKRN~9E+2qv z`wdGHtw-G#)yKAB)$gg$yDnv_n?lHkVp`}Dx)|;~ovHN>nYE0|Z2eKK5EsVORoAQ9 zkE#cb){F?H0(i)F- zg)lY)EcYeUG}~vCgeiNb&#=EhXR$hr8`X6=WBsJA)0w;2ezL1Clikli z2>EFf)XdOI`+{eU-DUF&_4hl5{{G!xq`#9B`+H+ze@Bb{QWRZKZ+`)<>81WsnK~TD zQK2M!=2=7kfx>)o`yZ(vr>Gxy>s=qhZ*Q=Ep!Otu2ftxp3D+66AXcmJADng?ON)Bo zbgh~XaEyEx;lA6CY_k78TYnJNp~I^ECG_VvT2J}o=D6PEv-)FN*FjAE4G56>Fv1>J zT^}^Yz^xz9N)HBq560uF;}6gb)8gSZw?ZLdmy-~CW{|p9f{C*~+DdaAwuXpkIGlEJ z{}_vp1Z*vcllj|qpX;wHLoc{$D?=Z+YFwe?wPFP>=su=C_9&vKr;JDb7_1L495*5M z(L-jmsb$2B_QUN|*otn#L`PKnfO=qnh5*avqZB4HQVHNF4T0a5hXyeQTo?D50woj^ zF14=aZxxIJSLt!}vFp*i2n}B0D&4n!45EWd-O8fIE9hK9Va~#yyPku8X^O&-j_i96 zU3IC}ML^JeikmuQKnMHL#x~ovb(CDN5_aLPo|?j|{0ZvfhaP;VD)c8k z@4}fJsz<&|uleXwd}@k3dY(B6vlQ+B_+680xyubNV}DDght1P<`?tjWXLY;lN{0Hl zj>#W=QoHH?mU`x5mWVIb|B9o3>RogKuV<0=QnUhW@i@Nqva1I7bK1l2;4zTm)wY@^ z+u;{2w6X9Kd|49iz8G9w89If~%>E@1z}$mdzSr@zVfXteSB2jGHNnX0bGlwvuWjXg zUhNJpC^%Jncpr(b3JtND&DGE^@0kA{LL#q#k0B!Gkc&Kof5Y)f^SKOfJA5+14Hm#$ zY9HRk=9|ysG)H656yQG=WgP?7+Dg%b7_u_DJhce3LBCUdV!!^9TEuqkyb|%22SPBb zn~tM{&W&f)jr0IlwZ8#;ek+>#C5mLqH?a_6JinxVFRae*+Za~6A0v9nl^oN02W`vf-gX1TkjVltRZ;9gS~*XuMY%tO ztOv&dKjZtzsYMsm&A&igkbvDWaA-aGV^C3YaT)H(^O3th`XXDjnLO1##HO^mH=~U# zw#4}hgdaZh zHdGlpu694dN#Pkp3*+IWY%U^)b;If1Akumxq+e@cp4NND)}*$FaXqgV?N>K7Q71{W zACZc71P$|F)+Fd7B5CP9?Bs`|=s$?iaY8y$?fxeq&^iyC8X4EXAb6PS@c`})vyLG0 z$JJe4lX)ChWmEm%6!^$ANA4 zY?9NHp8+h7gIwfrE)?C58)7!@5H9HY`hk*T%P8Pb$LBqW=WE7b)MWIp>UdFo^gWt? zy3f#m+%S+Xw7LvNUm*;SSCWgFW%n1M0c$${8#wBsrWew^`0#n3p!;&1u13&JBD(D) zXLlv<#%$C;7759c`m^&3z9_IWo}@E|_)Pi$ctfM2_hSSYcAMh-=QHrzIRDfH<(b_P z@Ap;p%vth9oN!`U|6P1pzf!C9Of3g`N!4w%KhBPvi(q{AeryU%EIBFcR27O zz3tNTUdKrx*16Y$r_t+d!%Ae(nEFd6=LAV6QBrv%jaXpB0wWd}vA~E0Ml3L5f&VWS z$e=q%q?MkwHjlmJHd|*98UH#P-tY|I{nD1kPTL$?o6ple+cvpr3eaB4V7JB73SR(h zQ>LY2oSl|$170RjoX#n>mUi3ZPMdwQ_cq1hD8gUo>R`*tm90KUBbQ&Ik-bnZTxCdt;!bqaDH@^!Yn|?nZBu za#F2y;;qA?lj>nWW$yMz@1^oyUH=v~#ot;F9HnqAOD~gtNRyoOUs2ABa<-Ed6|1vt z|J7%|9@{DKgL5-(G2QtpD)9hHIq_e*ZpD`P(CfNB3-APBHQ=s6U2g`|hIIWB=w`jC z>rbOx2e=1t3*Z2iPwDzQfI9%mJ!{|6_3?mpe}p_Lzpd*HfICj>`ZmA_;7_Q0M%TlD zCjgIA`5oNu10Fi7>t+*vYk|X%YQSB98v&oj@z);0;}Cg@8S(+0fI9%|0JCwLyc3X; zcpX*Ac<5dwW34sgmYmVq+o2=j@o-hq-q7_TA}k-xUVAJ=b~>)`Yr1}(@D^)@rC^~t zc1`vMW$w+hZ@;5(5}}E|0oN9?86iX>d;r%5;0Fbcz+Jeu;nHB!hta-LVX^jPF0d4A zGr26bO=g$H-j!8pahe{;u-NBWY%WVdnZ;UWp;tyIQ5(O(wFUa(wFH(3H;C&yvJc4P zQIYaKu4jPj15`ANE%RsBS+eJkWirUWUik4&szW0BLg2OocY+W^e+cxhv8>LCx*GVG zVkBM%SvC|KP*(t7PXtuvTdZB^kIghH<2sAgMGZ~?dIRe6C=TE+PO}i5^ks4y*Vs%y za+Q|sMYw?^e>?+vGj3<*BTw|3NMW;n~wT_p?V~;!71PdfSb-R97jh~SAg3s zaMaHn#MIMK_#1VK>4pr8t&Z7OL(z`9J&4CjK=RGwnPrxO$4v7rwjT3*i+x*Gxutm1 zD3@hc*XRY7YSSa=;cFKBZj|=JjvRYxzMuRyxg+GtXX!N z>j7OT*Q4VLnVm`Rr@=D?@%}dH#`=v2=rS#UvA%@13ckurpp9)qezC=>S+c9fwt{pJ ze5c;k^{)xNw`GEFld0Td>oS*vc|DkIi!B8j3#Nr*%TYwsp#{WN#n?!q`5tMEg8ScqJtv2{xr#{GN%nOO&bwKkZ{E(JR}O@Inzw#X>=eUnMe9@35f5 z%e5@Y_T?U3(&2oc^$*Q46mH&F`Ml_uF3M6bY_V8>5O9@%8wA`g;7$Rb5il&^fPjMno)=It{)mq^`I-L9-xV&HH*dDhUb`aL z9the>9n&4fQ%i!ZRI+|rv7>mp17`#b|LJ1smzD)73{Ahwl)q0mN}*qI%HN+S;>~Oi zY{XICgy@!WZmcsYqZAq6iTKfq90!T`Y(>ULB7Te_<0=tJ!ZpdBZlWqIAyV>#LL-on1`k%v+|xk zjNv(&=ThPEE+C@Qs)X^Xy)n^zMELVb_%92*oc|?{3z-aqo??1r-N{H{pEW{`Tu-E& zE+*$%JO2>)RQjju*H4+8Jmty^DR~&ClJhH8oqvrVf5YTl>&KS`y2jMA*8GZ z-gd2cMT4F;<=Xh{Nh9Y6z*|9|9JfD8Lw^h5Gp`-LKTkt{l;N$)PO;#!vn1$!Qs4*0 zLViN@C7MRghk|}UEaZKHp3Zzy*=Ir;{9O!hRR$+;0qm>@*$Uvv4$0$s3Gg=2mJzG~ zeq82Fid>&$U5Ajr6VLA{$upCOiQX{3iQ`LUyPzKw8zwp{qVx{}Z=J}??A!`jI!~o| zdmu^v9){<%#{Vx$`GVdjDql}Se?icPX~IEb=YoI-aal2M6eP`Wbbib9`8-i83Of@6 ztP*nM_j-z0EocS(2H{7!kD4IxHsP1ug5C~1wVNED(*?c!CKAOU!&BC?fFZzrYdf{1Dkoz)zBm8DJhdPJY}e{388@4{LbZR3Y;q7p*Z0%U#X?1%89sW+YC%9H=>>&1l&+C{~JhLP< z5q*QB#qII?J?q@Q_JDsKzP#aS^SQlv=CckghLRhUfdo>+lNq0u0@se5^LBuNr1waZrt-VAwes@S8*?>T3qI;axbo^STeuHUE?aNnom;M zS3(*)gG_6bjo%b;I_2j==o29H{g7b{eL$p>8Pm<4YRRvYBvO3RS5)XDCVsJk>aGfl0I!SjH5`Q585PR+-`hQ zN5r zVYAHPBRxLUFGu!))`PMf-U%_Rf|`kTxFG01D=%%BX9O}HY*PA zx^^_dVZhI8=s7B$s3sKLsPp?;J;We#t6Bq!gZbNmtYc*dij?W}H7bsPZ!L;61A)4O zjVFh%SYr|w_m}CDpS(@Il1;QSbCwc!W zp*(L9)r?l=3DPn^M;nrV7S;hu653d4_@6x2keHds82Rf2zl4R7P(ZnkN_1rsKRr95 zB%zI!hW{n=--jwZ_NDC=OY;6&LfO7JPKlp833A1Q4BdlBet91*;UF0oiTX%c*?)Q7 z(*hi+EcxaAx`gta3M4EEak9w1%|vHFl3(7}OK2CAQh!M&A-xDdXG*dx@Bby-Ao!(^ zr2Yb}^dLhvmHhIYK*F#fl|@=Lr(Zb2j29o@D+ffg_tp`PgdkB%s*7 z2qf1vUO?vixQzTZ!6)H)BQf$x;*}pK@ym0JIQ$g;XOj5kd5469pTbZ3S)=~gZ1E&u zt=Jz*F;ah7-iwTpU!I>xD9>vpee(F(FZiYZQH-t0H> z$L~{}oJ+PNWk^0*&H*#wWc%_y#p>h~$@J8gG0`~(;Y7H^t0X_2HyM-UH$oXZMmndH xmo#y}PPcz_?;-V(^AA-^-8l#Uzl9bSBzsc8$#!`y-Q%U^e`XP<__X=^zX5KkdUF5( literal 0 HcmV?d00001 diff --git a/examples/06_gamma_beta/example.c b/examples/06_gamma_beta/example.c new file mode 100644 index 0000000..77c86e3 --- /dev/null +++ b/examples/06_gamma_beta/example.c @@ -0,0 +1,50 @@ +#include "../../squiggle.h" +#include +#include +#include + +// Estimate functions + +int main() +{ + // set randomness seed + uint32_t* seed = malloc(sizeof(uint32_t)); + *seed = 1000; // xorshift can't start with 0 + + int n = 1000 * 1000; + for (int i = 0; i < n; i++) { + float gamma_0 = sample_gamma(0.0, seed); + // printf("sample_gamma(0.0): %f\n", gamma_0); + } + // printf("\n"); + + float* gamma_1_array = malloc(sizeof(float) * n); + for (int i = 0; i < n; i++) { + float gamma_1 = sample_gamma(1.0, seed); + // printf("sample_gamma(1.0): %f\n", gamma_1); + gamma_1_array[i] = gamma_1; + } + printf("gamma(1) summary statistics = mean: %f, std: %f\n", array_mean(gamma_1_array, n), array_std(gamma_1_array, n)); + free(gamma_1_array); + printf("\n"); + + float* beta_1_2_array = malloc(sizeof(float) * n); + for (int i = 0; i < n; i++) { + float beta_1_2 = sample_beta(1, 2.0, seed); + // printf("sample_beta(1.0, 2.0): %f\n", beta_1_2); + beta_1_2_array[i] = beta_1_2; + } + printf("beta(1,2) summary statistics: mean: %f, std: %f\n", array_mean(beta_1_2_array, n), array_std(beta_1_2_array, n)); + free(beta_1_2_array); + printf("\n"); + + free(seed); +} + +/* +Aggregation mechanisms: +- Quantiles (requires a sort) +- Sum +- Average +- Std +*/ diff --git a/examples/06_gamma_beta/makefile b/examples/06_gamma_beta/makefile new file mode 100644 index 0000000..ef385f7 --- /dev/null +++ b/examples/06_gamma_beta/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