From 4e079e9015a4c3ae6214e337314f11daadca849f Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Thu, 17 Aug 2023 14:19:32 +0200 Subject: [PATCH] move squiggle.c to its own repo, fix readme typo. --- C/squiggle_c/examples/01_one_sample/example | Bin 17696 -> 0 bytes C/squiggle_c/examples/01_one_sample/example.c | 50 -------- C/squiggle_c/examples/01_one_sample/makefile | 53 --------- C/squiggle_c/examples/02_many_samples/example | Bin 17696 -> 0 bytes .../examples/02_many_samples/example.c | 59 ---------- .../examples/02_many_samples/makefile | 53 --------- C/squiggle_c/squiggle.h | 109 ------------------ C/squiggle_c/to-do.md | 9 -- README.md | 2 +- python/beta/beta.py | 6 + python/beta/beta_array.py | 9 ++ squigglepy/scratchpad/samples_correlated.py | 10 ++ 12 files changed, 26 insertions(+), 334 deletions(-) delete mode 100755 C/squiggle_c/examples/01_one_sample/example delete mode 100644 C/squiggle_c/examples/01_one_sample/example.c delete mode 100644 C/squiggle_c/examples/01_one_sample/makefile delete mode 100755 C/squiggle_c/examples/02_many_samples/example delete mode 100644 C/squiggle_c/examples/02_many_samples/example.c delete mode 100644 C/squiggle_c/examples/02_many_samples/makefile delete mode 100644 C/squiggle_c/squiggle.h delete mode 100644 C/squiggle_c/to-do.md create mode 100644 python/beta/beta.py create mode 100644 python/beta/beta_array.py create mode 100644 squigglepy/scratchpad/samples_correlated.py diff --git a/C/squiggle_c/examples/01_one_sample/example b/C/squiggle_c/examples/01_one_sample/example deleted file mode 100755 index b3db352fd1b79c50dd3ef781efa445004e27efc9..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 17696 zcmeHPe{kEym0wwM;si%unb^gILUbT{4iLo-A=s28l4B?GjRGODivvj%#g=1VZOh0K zPPkq$bpxousXJGWAHz*UuhT-gx$)i6yY!|^{38UD0^=}11A(i_@Ixj*j1xkd26Ugd z-+l5|-$*Uj;qIUF&V2QDKl|R>-M4S|ySx6b@2>W*%CT4k2dh{vWNNOCcvQxgWy~Rj zN0f^>_?<5<6EgtY@OV_VnI6f}GnnQnJ}W566;mVw{kjedrd&g!B$q8btO{T%96Be- z{_$}b-%ys&iBysP^UX+yv0t+V~U|A9wY2bD=Y6iA-C zn9>v{*CZaj{r4cIV?+;!UTPznEj)va3o~S10L+;IzZ+O4{y-M|ds*b0pwRA>9 zk=ivYTf^<4+ThmKkO%~}wT0UQozY+<8VKmD85vu{+eGIbk!X`ZtVx8ncQlF4mUb$p zK@}6|h_tjvn*t5Z-xeLQXeX3IZ^O3<2s8y-T18v1wKd!TykUDV(A3f%Y;C!V>Oe|V z9om9W+7tTQ}K$a==G`SweBUU_!2aod^881Yr|-;%G3<}&XoSj zzHC8`1$;DSJa}UBT4s?)4=F7E!@q}n=Zcq9`aYSqir0~*K=RF?Jg%b`ll)Baj=~CC zTISFUIHA(I{v4?vkl=Zj=e;rYers0qE8}|=p64N{sd7M zz^QBw`waMO20;f5IL$E}Mhtj9gP@}Z9340vrY$gSfoTg&TVUD((-!zYu)s;}va`OP z5A%Gxb6;O1gm2eCG$%Rg>p74&EK5&bw+`@R^2+rf+I&I$vr79Nx1ccjv^z?vKn zrqZ;qnjF7{)2QvGvbI%1))(o2a8vE4(hr)_y{7aoQ~GaB>6j_qW=e-l=}o3|oiXjb zbw63e`z>#sx3+FQx}(0)x7&7Q2|!=G%H=B?@b^Yt4u5Z}%b9o-&D2VA|IsdlKYrRD zKk=MMz`!-}1HRrlm%{YE-rR41H2Jjl5p6-0%R#Y0fBaaJ);b8D1t<(Wl@OwwQQ~n* z9xx11sxpuGN^BZI_G+j8d4@*~ADlEye zdPpRT{`DFmPW7Rz*AapJUM5m;CR&7{vxnjsIi^d0^pS_lOpr64akrnj5rX?%4rqihkoy3W_LtdgR>_`A`3-|l6rkVwq> zJef>1qKXG|mmzK01X*E}AMZk@st@TSn;P>^YK~l2z91RltwS#C@v67BCVmRe82`h5 z+H=7vi#5j#t>lPSGI~<0e@i4_;TMS$=aNY%N7TD} z4LbD1>62c~;d0Th-(`c}P^{LKM`WWb|D?8g49^I}ouE|C3yAw{NF1SB12_y|Z5|?T zAX1+n9~m5*Q#J^aD3D|6;%>qBHxC3>YL}Q+WQYvH;+0+ki*y)y|V1nhkR_S)8MhC$shiV@Lo3xS< zy;irO0Z=nY;>5OGveoC_07{J>|Ku$E=o-iz|F`t`rS?4!rC+9f z_kd5fuOEWh+t(!7>T|b(%Jvmd`~3JZw{Q1My)8rWQSu@p@S@T9%;6^aev^GWkg>h8 zH{&zOp6}tUG_3fQ(X)$?z%L;nMF)rfl)b-CtY3{k!cG<)ZfKap+8r-l;~XJFgsnx-0~)(g^pFQ3#lj$~{D*$_j0L(xz1o|gu} zz^EU>Y^2`Ty8odAHsVCy$O?He%pvH-kke$ED+1!A~1}N3FqY*8mgN+aBs& zYpHl`=jFbJmwN_sBvNs3=Um^z6upnv1ynSx2ZqZAD0+Ozuzo=6G2<1?=s9cId3nzV zIr9LjTpg$ z@55r`hKaeS;Kf+B)S>&~4016;ix_wXEnEs5-rY0w>CW_8@%GlbzEJDgp&U}-e+D{e z)gu=-wa?&n^&Cw1+t@L0d_cv99gpB6V~CJl=l2d{(huwGvL>EIffW;0!{Ds+san zRxU*8Z)N=-CzJhz+a0Uy#ebz;(3RIMR?NHMs%u@B6HM|O@sz_hQjkCd!lQUff$s(7 zKs@yBD1oi8(3OF(y?9)}urHRO9KTQFaRR#v^%RwMM{iDzz4!rZmEC!-&1)~|$@SYk zIpeeKC0@JJYcB>zr9IDkfhY$11fD|J{WVbGvpaK!X4~_s^}KxOJCY%f>QREc$Iwq0 zwlWZQBd{lcN&6s881?S~VDxua3REWb^#rgX0;tS?*2^URZsh$G{r*kLLm>Vi0UHH& zgaG0n0l%t$Kk|MAe|XcFhx$8_egS#)@QwMPqT250$ysf8T4z{h+Z|qNB$YRTyju9p zlZaCw%m&|}HSZFM2;l)O2KIIM(|k~|%I?^cQ*AHqwN~4m57?^gCHLlf?WH|4R@*(+ zCg?AP*==cHZQ@qDoX?OP8D(xi?kdAw2clDqks;fhVO5ErCQRu^`<0Q@cibk8rL>_y@F==g z(JhL0DB7)PpQ3ve?N@X_QC-`Y|IPF0ysKBPyun#gw>8!tjX5jaOWdUk%VIKBwsUc* zyL5@WY=MsdY2Eb2XK@BY(l?^z8F8zq&xqT^cQWF+DL*#W#)|Js*?-gVnS%X19nTZo z|LOQF!G4pD&lc>T>G%aHKTgLnsx$P96}Kpy8F9N9&4_DaEFX*Lu8RlteARDBK zm+&nnoTfj*dEVE#g`lteEr$A2{U|izv1G8@Zs4Sc-y_6J;UPR0%)a~{kRX1t)&Jim zf1$XM4x&sUQ+p3co}%;m@m1g!%-ZID{4d3CuFoeDKi{uxXgFqv^ZI=UCgvP^3!40V z9{4XLkQ#`B5OejC(_vCLB?o{P2Qsw?Oi=Owh z;AgVn7r~xn2jNM#vrF1J9i!B>zzcFN7W}?p4iDkfGN+E%RL0G~sozW0LWa-U=!`6r z9m0ygPb~~WO3&R{_(zfNz<4b7C`Lv8PU-Qe1&%!TL>_%$CH+v@neF*wg>O*{FnN{; z{wJAV*`L=dbk0oWn*C=$@$-Bx&j6w4S&6I48OZMx|6b)kM^xTNz{#FCxR-(LIl=nX zLQlN<&n$FMm-wgcnbQO9kv_06?(He3$h_;sP(O3kraC>L8BNFa#mx_i+v@;sR-+day z)Xkhgv@Os;cXR0KOk+5(tu?$g*cxbzh9jMUU~Ic+2$RZasL@?sy0k1a6Wyz62?QgN z;Eq73JsR18+cd$pP@pl^*0uv8hExEY(KOaTVAVSBn(9FHO;vOoD4n5E1gbXNIu)#3F(Aa0IuHBp7ZB#M)b$0H7}E zHx{MO_pr=*Q%0puRX45Dc~pbw&uDCiXjp8A zb2Ya#MXz0q`fTm&RCSa0umXlFT~oT~x)m80c2rhfTf+eCqcB_!n-ZZrTxuAwUCDiM z6^&hu?xB(|)r&Qe9(X^+l=t^4Tk>+8=oNTqPtEeY z?_%mio-yz~nlW0x80G6R|59LTUQY6RM(4t3$%q)`c^}8r$&8AcnXvBDm=w#~@+6zvTP;W`0gH zrTwR7kmqxW=7pqf%3t-=RGQ`K%P{g?t|sq$`TsMpVQBnOU^%wqQ^3f6DS0{=KDvs#xO`rhJCwYcpUN_Z0#&|ZPi3+E)hR@e zu)LA7?-*G -#include -#include - -// Estimate functions -float sample_0(uint32_t* seed) -{ - return 0; -} - -float sample_1(uint32_t* seed) -{ - return 1; -} - -float sample_few(uint32_t* seed) -{ - return random_to(1, 3, seed); -} - -float sample_many(uint32_t* seed) -{ - return random_to(2, 10, seed); -} - -int main(){ - // set randomness seed - uint32_t* seed = malloc(sizeof(uint32_t)); - *seed = 1000; // xorshift can't start with 0 - - float p_a = 0.8; - float p_b = 0.5; - float p_c = p_a * p_b; - - int n_dists = 4; - float weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 }; - float (*samplers[])(uint32_t*) = { sample_0, sample_1, sample_few, sample_many }; - - float result_one = mixture(samplers, weights, n_dists, seed); - printf("result_one: %f\n", result_one); -} - -/* -Aggregation mechanisms: -- Quantiles (requires a sort) -- Sum -- Average -- Std -*/ diff --git a/C/squiggle_c/examples/01_one_sample/makefile b/C/squiggle_c/examples/01_one_sample/makefile deleted file mode 100644 index 302a1499..00000000 --- a/C/squiggle_c/examples/01_one_sample/makefile +++ /dev/null @@ -1,53 +0,0 @@ -# Interface: -# make -# make build -# make format -# make run - -# Compiler -CC=gcc -# CC=tcc # <= faster compilation - -# Main file -SRC=example.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/C/squiggle_c/examples/02_many_samples/example b/C/squiggle_c/examples/02_many_samples/example deleted file mode 100755 index b3db352fd1b79c50dd3ef781efa445004e27efc9..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 17696 zcmeHPe{kEym0wwM;si%unb^gILUbT{4iLo-A=s28l4B?GjRGODivvj%#g=1VZOh0K zPPkq$bpxousXJGWAHz*UuhT-gx$)i6yY!|^{38UD0^=}11A(i_@Ixj*j1xkd26Ugd z-+l5|-$*Uj;qIUF&V2QDKl|R>-M4S|ySx6b@2>W*%CT4k2dh{vWNNOCcvQxgWy~Rj zN0f^>_?<5<6EgtY@OV_VnI6f}GnnQnJ}W566;mVw{kjedrd&g!B$q8btO{T%96Be- z{_$}b-%ys&iBysP^UX+yv0t+V~U|A9wY2bD=Y6iA-C zn9>v{*CZaj{r4cIV?+;!UTPznEj)va3o~S10L+;IzZ+O4{y-M|ds*b0pwRA>9 zk=ivYTf^<4+ThmKkO%~}wT0UQozY+<8VKmD85vu{+eGIbk!X`ZtVx8ncQlF4mUb$p zK@}6|h_tjvn*t5Z-xeLQXeX3IZ^O3<2s8y-T18v1wKd!TykUDV(A3f%Y;C!V>Oe|V z9om9W+7tTQ}K$a==G`SweBUU_!2aod^881Yr|-;%G3<}&XoSj zzHC8`1$;DSJa}UBT4s?)4=F7E!@q}n=Zcq9`aYSqir0~*K=RF?Jg%b`ll)Baj=~CC zTISFUIHA(I{v4?vkl=Zj=e;rYers0qE8}|=p64N{sd7M zz^QBw`waMO20;f5IL$E}Mhtj9gP@}Z9340vrY$gSfoTg&TVUD((-!zYu)s;}va`OP z5A%Gxb6;O1gm2eCG$%Rg>p74&EK5&bw+`@R^2+rf+I&I$vr79Nx1ccjv^z?vKn zrqZ;qnjF7{)2QvGvbI%1))(o2a8vE4(hr)_y{7aoQ~GaB>6j_qW=e-l=}o3|oiXjb zbw63e`z>#sx3+FQx}(0)x7&7Q2|!=G%H=B?@b^Yt4u5Z}%b9o-&D2VA|IsdlKYrRD zKk=MMz`!-}1HRrlm%{YE-rR41H2Jjl5p6-0%R#Y0fBaaJ);b8D1t<(Wl@OwwQQ~n* z9xx11sxpuGN^BZI_G+j8d4@*~ADlEye zdPpRT{`DFmPW7Rz*AapJUM5m;CR&7{vxnjsIi^d0^pS_lOpr64akrnj5rX?%4rqihkoy3W_LtdgR>_`A`3-|l6rkVwq> zJef>1qKXG|mmzK01X*E}AMZk@st@TSn;P>^YK~l2z91RltwS#C@v67BCVmRe82`h5 z+H=7vi#5j#t>lPSGI~<0e@i4_;TMS$=aNY%N7TD} z4LbD1>62c~;d0Th-(`c}P^{LKM`WWb|D?8g49^I}ouE|C3yAw{NF1SB12_y|Z5|?T zAX1+n9~m5*Q#J^aD3D|6;%>qBHxC3>YL}Q+WQYvH;+0+ki*y)y|V1nhkR_S)8MhC$shiV@Lo3xS< zy;irO0Z=nY;>5OGveoC_07{J>|Ku$E=o-iz|F`t`rS?4!rC+9f z_kd5fuOEWh+t(!7>T|b(%Jvmd`~3JZw{Q1My)8rWQSu@p@S@T9%;6^aev^GWkg>h8 zH{&zOp6}tUG_3fQ(X)$?z%L;nMF)rfl)b-CtY3{k!cG<)ZfKap+8r-l;~XJFgsnx-0~)(g^pFQ3#lj$~{D*$_j0L(xz1o|gu} zz^EU>Y^2`Ty8odAHsVCy$O?He%pvH-kke$ED+1!A~1}N3FqY*8mgN+aBs& zYpHl`=jFbJmwN_sBvNs3=Um^z6upnv1ynSx2ZqZAD0+Ozuzo=6G2<1?=s9cId3nzV zIr9LjTpg$ z@55r`hKaeS;Kf+B)S>&~4016;ix_wXEnEs5-rY0w>CW_8@%GlbzEJDgp&U}-e+D{e z)gu=-wa?&n^&Cw1+t@L0d_cv99gpB6V~CJl=l2d{(huwGvL>EIffW;0!{Ds+san zRxU*8Z)N=-CzJhz+a0Uy#ebz;(3RIMR?NHMs%u@B6HM|O@sz_hQjkCd!lQUff$s(7 zKs@yBD1oi8(3OF(y?9)}urHRO9KTQFaRR#v^%RwMM{iDzz4!rZmEC!-&1)~|$@SYk zIpeeKC0@JJYcB>zr9IDkfhY$11fD|J{WVbGvpaK!X4~_s^}KxOJCY%f>QREc$Iwq0 zwlWZQBd{lcN&6s881?S~VDxua3REWb^#rgX0;tS?*2^URZsh$G{r*kLLm>Vi0UHH& zgaG0n0l%t$Kk|MAe|XcFhx$8_egS#)@QwMPqT250$ysf8T4z{h+Z|qNB$YRTyju9p zlZaCw%m&|}HSZFM2;l)O2KIIM(|k~|%I?^cQ*AHqwN~4m57?^gCHLlf?WH|4R@*(+ zCg?AP*==cHZQ@qDoX?OP8D(xi?kdAw2clDqks;fhVO5ErCQRu^`<0Q@cibk8rL>_y@F==g z(JhL0DB7)PpQ3ve?N@X_QC-`Y|IPF0ysKBPyun#gw>8!tjX5jaOWdUk%VIKBwsUc* zyL5@WY=MsdY2Eb2XK@BY(l?^z8F8zq&xqT^cQWF+DL*#W#)|Js*?-gVnS%X19nTZo z|LOQF!G4pD&lc>T>G%aHKTgLnsx$P96}Kpy8F9N9&4_DaEFX*Lu8RlteARDBK zm+&nnoTfj*dEVE#g`lteEr$A2{U|izv1G8@Zs4Sc-y_6J;UPR0%)a~{kRX1t)&Jim zf1$XM4x&sUQ+p3co}%;m@m1g!%-ZID{4d3CuFoeDKi{uxXgFqv^ZI=UCgvP^3!40V z9{4XLkQ#`B5OejC(_vCLB?o{P2Qsw?Oi=Owh z;AgVn7r~xn2jNM#vrF1J9i!B>zzcFN7W}?p4iDkfGN+E%RL0G~sozW0LWa-U=!`6r z9m0ygPb~~WO3&R{_(zfNz<4b7C`Lv8PU-Qe1&%!TL>_%$CH+v@neF*wg>O*{FnN{; z{wJAV*`L=dbk0oWn*C=$@$-Bx&j6w4S&6I48OZMx|6b)kM^xTNz{#FCxR-(LIl=nX zLQlN<&n$FMm-wgcnbQO9kv_06?(He3$h_;sP(O3kraC>L8BNFa#mx_i+v@;sR-+day z)Xkhgv@Os;cXR0KOk+5(tu?$g*cxbzh9jMUU~Ic+2$RZasL@?sy0k1a6Wyz62?QgN z;Eq73JsR18+cd$pP@pl^*0uv8hExEY(KOaTVAVSBn(9FHO;vOoD4n5E1gbXNIu)#3F(Aa0IuHBp7ZB#M)b$0H7}E zHx{MO_pr=*Q%0puRX45Dc~pbw&uDCiXjp8A zb2Ya#MXz0q`fTm&RCSa0umXlFT~oT~x)m80c2rhfTf+eCqcB_!n-ZZrTxuAwUCDiM z6^&hu?xB(|)r&Qe9(X^+l=t^4Tk>+8=oNTqPtEeY z?_%mio-yz~nlW0x80G6R|59LTUQY6RM(4t3$%q)`c^}8r$&8AcnXvBDm=w#~@+6zvTP;W`0gH zrTwR7kmqxW=7pqf%3t-=RGQ`K%P{g?t|sq$`TsMpVQBnOU^%wqQ^3f6DS0{=KDvs#xO`rhJCwYcpUN_Z0#&|ZPi3+E)hR@e zu)LA7?-*G -#include -#include -#include "../../squiggle.h" - -// Estimate functions -float sample_0(uint32_t* seed) -{ - return 0; -} - -float sample_1(uint32_t* seed) -{ - return 1; -} - -float sample_few(uint32_t* seed) -{ - return random_to(1, 3, seed); -} - -float sample_many(uint32_t* seed) -{ - return random_to(2, 10, seed); -} - -int main(){ - // set randomness seed - uint32_t* seed = malloc(sizeof(uint32_t)); - *seed = 1000; // xorshift can't start with 0 - - float p_a = 0.8; - float p_b = 0.5; - float p_c = p_a * p_b; - - int n_dists = 4; - float weights[] = { 1 - p_c, p_c / 2, p_c / 4, p_c / 4 }; - float (*samplers[])(uint32_t*) = { sample_0, sample_1, sample_few, sample_many }; - - int n_samples = 1000000; - float* result_many = (float *) malloc(n_samples * sizeof(float)); - for(int i=0; i&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/C/squiggle_c/squiggle.h b/C/squiggle_c/squiggle.h deleted file mode 100644 index 5409e01c..00000000 --- a/C/squiggle_c/squiggle.h +++ /dev/null @@ -1,109 +0,0 @@ -#include -#include -#include - -const float PI = 3.14159265358979323846; - -// Pseudo Random number generator - -uint32_t xorshift32 -(uint32_t* seed) -{ - // Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" - // See - // https://en.wikipedia.org/wiki/Xorshift - // Also some drama: , - - uint32_t x = *seed; - x ^= x << 13; - x ^= x >> 17; - x ^= x << 5; - return *seed = x; -} - -// Distribution & sampling functions - -float rand_0_to_1(uint32_t* seed){ - return ((float) xorshift32(seed)) / ((float) UINT32_MAX); -} - -float rand_float(float max, uint32_t* seed) -{ - return rand_0_to_1(seed) * max; -} - -float ur_normal(uint32_t* seed) -{ - float u1 = rand_0_to_1(seed); - float u2 = rand_0_to_1(seed); - float z = sqrtf(-2.0 * log(u1)) * sin(2 * PI * u2); - return z; -} - -float random_uniform(float from, float to, uint32_t* seed) -{ - return rand_0_to_1(seed) * (to - from) + from; -} - -float random_normal(float mean, float sigma, uint32_t* seed) -{ - return (mean + sigma * ur_normal(seed)); -} - -float random_lognormal(float logmean, float logsigma, uint32_t* seed) -{ - return expf(random_normal(logmean, logsigma, seed)); -} - -float random_to(float low, float high, uint32_t* seed) -{ - const float NORMAL95CONFIDENCE = 1.6448536269514722; - float loglow = logf(low); - float loghigh = logf(high); - float logmean = (loglow + loghigh) / 2; - float logsigma = (loghigh - loglow) / (2.0 * NORMAL95CONFIDENCE); - return random_lognormal(logmean, logsigma, seed); -} - -// Array helpers -float array_sum(float* array, int length) -{ - float output = 0.0; - for (int i = 0; i < length; i++) { - output += array[i]; - } - return output; -} - -void array_cumsum(float* array_to_sum, float* array_cumsummed, int length) -{ - array_cumsummed[0] = array_to_sum[0]; - for (int i = 1; i < length; i++) { - array_cumsummed[i] = array_cumsummed[i - 1] + array_to_sum[i]; - } -} - -// Mixture function -float mixture(float (*samplers[])(uint32_t*), float* weights, int n_dists, uint32_t* seed) -{ - // You can see a simpler version of this function in the git history - // or in C-02-better-algorithm-one-thread/ - float sum_weights = array_sum(weights, n_dists); - float* cumsummed_normalized_weights = (float*) malloc(n_dists * sizeof(float)); - cumsummed_normalized_weights[0] = weights[0]/sum_weights; - for (int i = 1; i < n_dists; i++) { - cumsummed_normalized_weights[i] = cumsummed_normalized_weights[i - 1] + weights[i]/sum_weights; - } - - float p = random_uniform(0, 1, seed); - float result; - for (int k = 0; k < n_dists; k++) { - if (p < cumsummed_normalized_weights[k]) { - result = samplers[k](seed); - break; - } - } - - free(cumsummed_normalized_weights); - return result; -} diff --git a/C/squiggle_c/to-do.md b/C/squiggle_c/to-do.md deleted file mode 100644 index 3be7ffb1..00000000 --- a/C/squiggle_c/to-do.md +++ /dev/null @@ -1,9 +0,0 @@ - -- [ ] Add example for only one sample -- [ ] Add example for many samples -- [ ] Use gcc extension to define functions nested inside main. -- [ ] Use OpenMP for acceleration -- [ ] Chain various mixture functions -- [ ] Have some more complicated & realistic example -- [ ] Add summarization functions, like mean, std, 90% ci (or all c.i.?) -- [ ] Add beta distribution diff --git a/README.md b/README.md index b9806588..9a8ef576 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ result = mixture(dists, weights) # should be 1M samples mean(result) ``` -I don't particularly care about the speed of this particular example, but rather think that the speed in this simple exaxmple would be indicative of the speed when considering 100x or 1000x more complicated models. As of now, it may also be useful for checking the validity of simple estimations. +I don't particularly care about the speed of this particular example, but rather think that the speed in this simple example would be indicative of the speed when considering 100x or 1000x more complicated models. As of now, it may also be useful for checking the validity of simple estimations. The title of this repository is a pun on two meanings of "time to": "how much time does it take to do x", and "let's do x". "BOTEC" stands for "back of the envelope calculation". diff --git a/python/beta/beta.py b/python/beta/beta.py new file mode 100644 index 00000000..58bffe7e --- /dev/null +++ b/python/beta/beta.py @@ -0,0 +1,6 @@ + +import numpy as np + +for i in range(1000* 1000): + x=np.random.beta(1,2) + diff --git a/python/beta/beta_array.py b/python/beta/beta_array.py new file mode 100644 index 00000000..d182931d --- /dev/null +++ b/python/beta/beta_array.py @@ -0,0 +1,9 @@ +import numpy as np + +n = 1000 * 1000 + +def sample_beta_1_2(): + return np.random.beta(1,2) + +a = np.array([sample_beta_1_2() for _ in range(n)]) +print(np.mean(a)) diff --git a/squigglepy/scratchpad/samples_correlated.py b/squigglepy/scratchpad/samples_correlated.py new file mode 100644 index 00000000..2e790c59 --- /dev/null +++ b/squigglepy/scratchpad/samples_correlated.py @@ -0,0 +1,10 @@ +import squigglepy as sq +import numpy as np + +a = sq.to(1, 3) +b = a / 2 +c = b / a + +c_samples = sq.sample(c, 10) + +print(c_samples)