From ccad14b318de40cd9d8207d17d957e273df82aaa Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Sat, 23 Sep 2023 22:15:48 +0100 Subject: [PATCH] leave out really trivial manipulations, add example, update to-dos --- README.md | 4 +-- examples/11_algebra/example | Bin 0 -> 22416 bytes examples/11_algebra/example.c | 26 +++++++++++++++ examples/11_algebra/makefile | 53 ++++++++++++++++++++++++++++++ squiggle.c | 59 ++++++++++++---------------------- squiggle.h | 10 ++---- 6 files changed, 104 insertions(+), 48 deletions(-) create mode 100755 examples/11_algebra/example create mode 100644 examples/11_algebra/example.c create mode 100644 examples/11_algebra/makefile diff --git a/README.md b/README.md index 7b46a40..fd2af7a 100644 --- a/README.md +++ b/README.md @@ -345,8 +345,8 @@ It emits one warning about something I already took care of, so by default I've - [x] Add prototypes - [x] Use named structs - [x] Add to header file - - [ ] Test results - [ ] Provide example - - [ ] Add conversion between 90%ci and parameters. + - [ ] Test results + - [ ] Add conversion between 90% ci and parameters. - [ ] Move to own file? Or signpost in file? - [ ] Disambiguate sample_laplace--successes vs failures || successes vs total trials as two distinct and differently named functions diff --git a/examples/11_algebra/example b/examples/11_algebra/example new file mode 100755 index 0000000000000000000000000000000000000000..1ac404bf2002375fdd55c03fde8845f94d1b6abb GIT binary patch literal 22416 zcmeHPe|%Kcm4C@(AV4y2K+yPuGTNcWA50PigeuDnB=Dj^0|X5i7=~mf$(m&1WQGLm zYBU6xLfo|LXSF}vjjp?;+aHVV{=jzIPQnj^S{tlxvDPL0P-h5e92HTJx8HN$y^}XD zGt#yF^sn4cCgWzWZk0LnUQPGcz(IQB2ZbNz^L5M8i0F$H1Mu0^pS9 zOS$+zL;8|54(%+4)B2qnK(AH~XC}jJjyC}kUmjIbz^gR0B%sicAn~OuZ4i<%d}NnDwg@+`il=KRjoB{*%Ak`S4ADIDYLp`~5-&jY%|A zh+nvv+Ek}Yk!ZB}Ka84&BY4RBawEZXQ_=+_6 z>NNQNH27}WH+d^O8ymcm+g;n}Z*sQ;J%OOxt&yfA^tLri z4gOlG<-S1B$Np2(CpG(9rIz|88m^iT=WY(vHwAs}>biR+Z(Ds3gFt@uy@YW4JoODy zqo<+4UrqQX9qoNB{s0kIw|U&Y`X*09{r#jMN#+VRdN7hN;Ppz&%1Re6anE+lcFY~A z-#k)Xkf_c{)NXRjA;V;{YsQ2t&Wz7E)*6rhnG*e%@IMPMgMGAX)(LFrih8URYn(L1+hNwulFpz_h4`BRiI&RL`ubeD`f2GToF>ndkQcn1wY|AuvVF(>Ld ztHX0u7(*S$({jhB#PH2JygqNW>F|{2L}}OI*?P=7w(IcvJm0Ou^QD$?JfXu28&OA( z4o|sElwCS}jz9wL*5TpMMA@UmTM}3#c^zKX;rn&?$vXUi4xcIy=t#$^V-^^*z?cQb zEHGw)znBF+kiYb)(s3$V*bmv1b!)+7QHJPbJ#45s1gWk!aHfARc=-(WZqz z9=k8mrUgG9^Ca4|;fcpqB-*smiN{=tHf;prv7$tqHURP1)rmGO{PEai(Z+1w{41Jm zvBc(!_CiB@uAx24(4J{%PdBuu8QQsqw%O4B495X_eLga@PaE2A8QO0e+OO%`${x(0 zijyM2r^rnM0%dhl^RQ?RFQVqKXl6Toa{d8A^+%Unj_BE8F8-1v{fgp2>6*K1x+H-h zhR>iTW{3Gk^izERwE4W=kGh@LBYVjC(P{7nuop-Ir#{K+hk-xC>k-tC@p^yMrm4^i zgiygtody-Y&lyOCe-_NBKgH`@g>PO^)uITy5_}GM#uhS=#Nxf^bZPj|dL?{J={PZ5 zUQyO%-b5pJ<*IoE7VZM4I_;5LCF#R%G{29E(q(>{z=E^EX&7of)p4m!PzU81?Id_F zhoIlg=;sCqTz5eU4=cUzE>e2Wnv{$K%ApIvsSq%i3&>UneJm~nKV6UQ3#nTQ-LP8e zSU82MO8Bi{jK~mPyDe~baH*D9ha>#KI8TPt@`AG z-qs045x{==@yK4^h;if@PIkTQ7_Wr;l-^^xN>QH@9#DFtlSOUNtP~A|$A=?&eU;G+ z;}Z;E9ZBAU=%g~C zk0|=$awVg>jHr?(#O8ymB2TYlnlVSFd`RhyL8BLxFltAk5o(vo+qc6L82KQL`2q|& z?R!@u@5pOVpIm0oMnK5h>D)xGd0+H}bE-;Hcm%|ZTV5NH^AF38{tx68Z+sw6uCiyf zZi2XCdp6N+1YPtkFd^>yNz@1O+G2Yax_S}Jd3&{ZfodJc$PNdULosE-A=2p}>~Jtv zO-@b(d|G?O+@Omt$2dbBI+Ibnz*ZQW`vXCc#$MywuF=ob{H-g<8db2ynqFwq2bX`2 zzlL5)R@--A#V$PVcX^T`#Eys#q4+FtY*K*LHkG8uMK3;bKtqS{bE^|xE4?8;{K1a-#Uyo z%6_ZjxFnVEImDq&yOg6d*-WgssxgmOa0L#qP?sl9h3D7yM$PC~c9;+fQ|;&lcONG4 z0>vzx-jB-wQM}&#FP9;A7Y(q;ri(kbix4R?o)(q5teYwYT_SIP(iqlu2;^=lL#R;- z`qrqpUiPu(7xmJ#ez{8C_7M138nQy5`3=xU*Ml+1Iai?z7dL_%)2Ym|Bere##A-n% z*1xoQ>5hxU_ezTIUz~w|S`G2Nid4}2Q;22Dyk{CVQ;)P`mkK-Q-UdaXwwWND^oFJ( z5^}yzb{r+~lXXH^wlx3qfZ?*#!k zFhL62>A6fI+>j3pmeK-xOD;lG7v6)Ejr#{oQrmA}ofZ2HtRYxuZ-za< zWWJsaiFku9mA%p1u{#5rHLcTg(aG`z1rtSGlF?FZe8{>l8WvV3ghZ=9Ah{#jB@{Dm zgoa*%PtU`PBijSEA6~75KSS<={){>d5;lWuCHx+6S$48~ zG@zlU(lu{ghq+RQgKP5p^jj`kk4EHGnsaby3K==}=xp$70+}hS=J#>V3Ncys=w%}# za@la!dU9W{UbfhucAj~#yyHSfusj@B!e{Bqg=>0`Whoh-c0+NvFT#9*x^kABqAchO zy&WFz7`)K&b_O=mPtdT>v;7JFwIejy^}m;P#qULij_3@;_s_IB8Q1Y?#)H>%oXp(A z?vFb9GLbi!Br?(gjG93Moz|<&>N#C0PzrM~h2Qq@=xx!EmTfFOEhHZLvdLJaE z6L@MAMQC@$=0*Agh3s_wyPap`M^9m{&_&*I8kKMyT_<3}LnmSDecXmD?k{4;1qb9U zw1487L8JXN!_6cye;A)XaH&(zB&0HZeB3`ix$0$~yoR={ulwX_43hPZB;0ff$+k{? zo)Q28{}oXF9^<@;tbQB=AYrU{rF8*nw9%Bvb@_i7dXY4VZNjBF;@1)iDVHWhLBD|@Cq~iJ~`?>CR-BlWX*|oYf{F$r56@GQ~Dkb~@3Fv%9 ze)J@w0g)kZc^4Hj5g%U^ull4lJCP@9p^J1{2NzSh5dqd zgUw=?`2~`4&OuN=%#7?lx#BZP>3ALi{b}e#X?Pqq7Htn;Y9smk2Y9MdbOzKf4&}i9 zO2>f=*MbipycU9rKf?WT)<*~rTEBWX9WM)i>{2Vt($M5?$O`@h^gu09D;Ay=DJEL| zE}I(LQsmvU&Q;}F=UVS_i#SiDyHR=APwP_RekZlvHfh4YsNG@*u|R0|UF2T!YwR9| zdhZPuQhQ@1Xo~#gswfQ~T%|mGDjSP`jj|(aWWb3&Pf%ql;Q^H50x1Ls2*YjC+A}W0znG;Yd4OiK%m1~t9n#3CA zCT)mNf%5P^8ni4tc#o@UFLQEnxVJ3)0*-US5FUL1@)$b$P1K^(@lTiSd@Ns){^v!E zJK@mjOTh{_?F@somg%NdY@c5g5%(i>lIcIHX;)Hd#7 z=||ZCBPQG^cRmXmuji(V-Uzk?|FC)A^NVQ(<>b&6Na;wV12Cx3M*Z!M2keEAggX9L*sYyV>k>8!xfaC8 zengCH4&tD_oRJYRy0Hp&UVvTDscbeD`6l9t%W$cQ8vd{~$2EEFdUZ_Jh!$Uh;Cssmc z?aqKsM{H;jyL8J?a8BN)-H#L;k++^8i3fQy!{o!Dti@EE};DgzQ5(DR%=>4pP(Y04ys-zaLfA%C#&9jmk_z z*N()`!J#d_?vxP^$5HMgKkI`2F$(&J_ffd;xOc6?79D5L>>duYYh_8c9dww3J}PS+YeBF zgw<@yj)}U!ML^@c85dc&U+EphZCWp@MX?5JwLzRDqBgPZHu#dAZ6X=1F?;ZMB)X4w z$H)SRUVo0gG>jTj15T=d%S#w3bOW-24c4KY!2FCMI*#)5(5Nz@pP!gkMuwh=I-xha z9{;@hK8q9n=&Lx7`E%>r{WjP4E|-=;ekIM#k#INh(JENm8r_Zl((nN#ytk}t7j@CT zBmADyu{YbLzKZK=$=R|6L!nbYL#QHzKK^CT)(lF2XA!W`OTY{qSstH*NW5#eI=@Q} z#pReo%}4h)VSt(Sk%fP?pJXwiL~UmPi0}%0z=F z_M&k1Nc{Z?@_&SRT~xm7?y|0h$8Fe^VB7mYwC(+fw!P<9D_yf*1$tH4b#%X&9ewCH zf}r;^Q*cJn8~nO23&+dC?-qyOaj8>}z%$rEz8N}+d1UMMbh;e}*^8c35+s>KsXe&_ za;%J5V9WxaYk>?2%Oy%ho!8b9YP9)%wmbX*e4#Ucrte1EOy6}jPg9K%Ro`O3Ol17c z0e?-XI!H3g{IyM^#f-*6(kLXhccNIGyvAWRhr^4UH`?l(Y%^PI`7>*-lN^o${I%Q{ zs;{kW@H(ooLvF%XV>fvkY}NjzAiic=-&AW0VvhJ?Y^^tFLr=(SYw`zeK7XjGh6(e9 znyQ2K{w4_C`b0GR^yp#m_1b)C+Y zH^fx+F>0Su)qQ|g+{~Q-YzEB2cQytA^8jtQ`@Ipc2hd4)oFY^KTH~s^3vdAN1hsLZ zotGs^!#K|04%mZ-96tms$67xCXvL#0E57T|444l%0O$mC;vwR8`VI&R&R1EC?JaB@$3f3j_MqID&JJqVZvLiOD%b~$&*^M+oeTUeD(U9 z?AH*Q_^a@-!N&UUyAnQ#Pag1tfTGY?JMkF=ekR#?rNwHRl3}5(0Tq(*Lwp9|o2viyBAQe2ylFsAP4ICR!vDh*XPPAk%o5W<~L5FAH z?|X?Dh3Ic2z2IxQ&yW;LVdfFgmri1R3UGG#-`HmZDeYJl#KaTCS7NbtWZrHmG~JXj z(PDL3vdbpj2=sRJ?W1^sUj>cRRG+EQES;43xK?xFj*)D&2lO8zo=Q+B`mLGEEVg@@ zK-Pa4{dU~ntVKkf#xKqGw<~Xv$9N;R4RP`3f zkWAAvEVebwPYRml=z9S1eIuZ>)MEYH%*B?xM@=ObTUSK*xoted!ybjYcOHbB%%dV{3 zEW5J;fG)JxPPF6|GgcU=5ajebs1tPka6K6BR+GzO>&RMW*8nO~n@5Rt%EgVZP1MV>b0|j4QS{GuxS^VFX;Wuym;sb}f~D-Y;&pF5?PgJI4y`tEk{>CrO2~N#*=^O|&(kQqK8l zou)!x*`|W87qODh6%_b7UJ-gnI313cSrOyoYr(7tI@~I-`C&MQy1?5}J~5obvVNf# zwoj}b;r+N3VQn!FF8tGNyK?F9#_5akFLXntKXIqZ;`Ai{AiM7$(9liAt`=bXhJZEyrLUcG7)j1 z9;4Mw_@%k{J-)uzgxz=dY&}NMPfPr+Ro`opE=~NdH5q@IlE`t#CTFU2)y+ck2uvl1-Mb_Ty`eV6>i5h94bQ=a`()ez)Q}&!klnA+0ZyVZbRsM)*1@*0rk{UQemj7pB4A z$>j|4g_NCdfTap}+h{S;y;w+Wv??3ot&PbU&9H;vt&;6}t#3Cc`3HgL3*^-t|6|~- zkZ+9V|6uf^)8LPg9K?@NpCLvMzwqnjUOx7xz}s|s zi0V`XSSr7Nh2gEzFyGh;{}yt*ljBRcpj$b<`zsoPodu#=2K+T*FdeKD^jyvg-q*_T zIznFE#_4x+KW-Ex-ZoETGkca&9SF7rLwM#|EseZsvOSJQh(TJTDb8o#@?!N1Yd z;I0Y!11)Y(s7jhn%u zYq>!gOd{1i?Qy%8u5>LgahKdtOfNl2rmd0O#cS?xEiYY?+{50H0uG=VDix`_83HSCv$_D_o1qN=PbuaY=Pch-r? zW4dcSjg6kqeSr&mX-X>I06mU^r9Zf79HyX(lO~pEhC+C zghoNXsArQqd^(z&g}JmB(a?KdbiJ^L=JS>yoZ}|d zwHLymk^K4^cTJ1mU56)-^lF>ZJ3&=b-{cOpcx$+X>QEz|S_+LC>)V2%fR{Hli*59` zNo@$qx_V!59;Qw2yD{j}Q0v_o@ZdE!?D3@=bJ1o+p`CC##UFpsO)@)oZGua}^F{Z9 zLN^{f`s!=EP1Rn+qxLXT!VuNomKL`+5by^QYTYzj@^Bq68E(=RKEp<@rwMv&1eX@E z?4>zw-Aj2!)S=ibFQxbQjNmBx=&e1;7`Gd*4dT}ep8(wH!z0_zd1nyWE?EqB-3QCTa&5c3NM!;Y|gLR_DEZr4VDV;eAh=nF+8C3pqPfMNTsM*|vAv72aXg%~u7!QY& z8gBFjybT^=;I-z4pyXh&;y~3=>qnC+E#7L$5%jj9Nl6dXes;+^ymdUy)zv_nXlXLE ze5Y}UI)?Hz)^ zdsq?Y3j&Jq#W|W$e-F6m9!&6yb2b47$+#$_BUJ^zIPaUJsQqUX6f5X5p&aD#Az+fEzE%4L( zzIfV6oNtNuM;hrxTfq0hE9>~h`H8X+NEAV*wMN+l2(P6wcSS)q@}Kh!OB*C@6BdADis aJBIiX;|kizsHEcm=^Yxy6$S=FQ~E!drvohj literal 0 HcmV?d00001 diff --git a/examples/11_algebra/example.c b/examples/11_algebra/example.c new file mode 100644 index 0000000..b26d0ad --- /dev/null +++ b/examples/11_algebra/example.c @@ -0,0 +1,26 @@ +#include "../../squiggle.h" +#include +#include +#include +#include + +int main() +{ + // set randomness seed + uint64_t* seed = malloc(sizeof(uint64_t)); + *seed = 1000; // xorshift can't start with 0 + + normal_params n1 = { .mean = 1.0, .std = 3.0 }; + normal_params n2 = { .mean = 2.0, .std = 4.0 }; + normal_params sn = algebra_sum_normals(n1, n2); + printf("The sum of Normal(%f, %f) and Normal(%f, %f) is Normal(%f, %f)\n", + n1.mean, n1.std, n2.mean, n2.std, sn.mean, sn.std); + + lognormal_params ln1 = { .logmean = 1.0, .logstd = 3.0 }; + lognormal_params ln2 = { .logmean = 2.0, .logstd = 4.0 }; + lognormal_params sln = algebra_product_lognormals(ln1, ln2); + printf("The product of Lognormal(%f, %f) and Lognormal(%f, %f) is Lognormal(%f, %f)\n", + ln1.logmean, ln1.logstd, ln2.logmean, ln2.logstd, sln.logmean, sln.logstd); + + free(seed); +} diff --git a/examples/11_algebra/makefile b/examples/11_algebra/makefile new file mode 100644 index 0000000..ef385f7 --- /dev/null +++ b/examples/11_algebra/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 1e03bb2..3bde4da 100644 --- a/squiggle.c +++ b/squiggle.c @@ -78,28 +78,29 @@ double sample_lognormal(double logmean, double logstd, uint64_t* seed) return exp(sample_normal(logmean, logstd, seed)); } -inline double sample_normal_from_95_confidence_interval(double low, double high, uint64_t* seed){ +inline double sample_normal_from_95_confidence_interval(double low, double high, uint64_t* seed) +{ // Explanation of key idea: - // 1. We know that the 90% confidence interval of the unit normal is + // 1. We know that the 90% confidence interval of the unit normal is // [-1.6448536269514722, 1.6448536269514722] // see e.g.: https://stackoverflow.com/questions/20626994/how-to-calculate-the-inverse-of-the-normal-cumulative-distribution-function-in-p // 2. So if we take a unit normal and multiply it by - // L / 1.6448536269514722, its new 90% confidence interval will be + // L / 1.6448536269514722, its new 90% confidence interval will be // [-L, L], i.e., length 2 * L - // 3. Instead, if we want to get a confidence interval of length L, - // we should multiply the unit normal by + // 3. Instead, if we want to get a confidence interval of length L, + // we should multiply the unit normal by // L / (2 * 1.6448536269514722) // Meaning that its standard deviation should be multiplied by that amount // see: https://en.wikipedia.org/wiki/Normal_distribution?lang=en#Operations_on_a_single_normal_variable // 4. So we have learnt that Normal(0, L / (2 * 1.6448536269514722)) // has a 90% confidence interval of length L - // 5. If we want a 90% confidence interval from high to low, - // we can set mean = (high + low)/2; the midpoint, and L = high-low, + // 5. If we want a 90% confidence interval from high to low, + // we can set mean = (high + low)/2; the midpoint, and L = high-low, // Normal([high + low]/2, [high - low]/(2 * 1.6448536269514722)) const double NORMAL95CONFIDENCE = 1.6448536269514722; double mean = (high + low) / 2.0; - double std = (high - low) / (2.0 * NORMAL95CONFIDENCE ); - return sample_normal(mean, std); + double std = (high - low) / (2.0 * NORMAL95CONFIDENCE); + return sample_normal(mean, std, seed); } double sample_to(double low, double high, uint64_t* seed) @@ -111,7 +112,7 @@ double sample_to(double low, double high, uint64_t* seed) // Then see code for sample_normal_from_95_confidence_interval double loglow = logf(low); double loghigh = logf(high); - return exp(sample_normal_from_95_confidence_interval(loglow, loghigh)); + return exp(sample_normal_from_95_confidence_interval(loglow, loghigh, seed)); } double sample_gamma(double alpha, uint64_t* seed) @@ -165,9 +166,10 @@ double sample_beta(double a, double b, uint64_t* seed) return gamma_a / (gamma_a + gamma_b); } -double sample_laplace(double successes, double failures, uint64_t* seed){ - // see - return sample_beta(successes + 1, failures + 1, seed); +double sample_laplace(double successes, double failures, uint64_t* seed) +{ + // see + return sample_beta(successes + 1, failures + 1, seed); } // Array helpers @@ -470,10 +472,9 @@ struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* se return result; } -// # Small algebra system +// # Small algebra manipulations -// Do algebra over lognormals and normals -// here I discover named structs, +// here I discover named structs, // which mean that I don't have to be typing // struct blah all the time. typedef struct normal_params_t { @@ -481,11 +482,6 @@ typedef struct normal_params_t { double std; } normal_params; -typedef struct lognormal_params_t { - double logmean; - double logstd; -} lognormal_params; - normal_params algebra_sum_normals(normal_params a, normal_params b) { normal_params result = { @@ -494,16 +490,11 @@ normal_params algebra_sum_normals(normal_params a, normal_params b) }; return result; } -normal_params algebra_shift_normal(normal_params a, double shift) -{ - normal_params result = { - .mean = a.mean + shift, - .std = a.std, - }; - return result; -} -// Also add stretching +typedef struct lognormal_params_t { + double logmean; + double logstd; +} lognormal_params; lognormal_params algebra_product_lognormals(lognormal_params a, lognormal_params b) { @@ -513,11 +504,3 @@ lognormal_params algebra_product_lognormals(lognormal_params a, lognormal_params }; return result; } -lognormal_params algebra_scale_lognormal(lognormal_params a, double k) -{ - lognormal_params result = { - .logmean = a.logmean + k, - .logstd = a.logstd, - }; - return result; -} diff --git a/squiggle.h b/squiggle.h index 75a9fca..50ba120 100644 --- a/squiggle.h +++ b/squiggle.h @@ -58,24 +58,18 @@ struct c_i { }; struct c_i get_90_confidence_interval(double (*sampler)(uint64_t*), uint64_t* seed); -// small algebra system +// small algebra manipulations typedef struct normal_params_t { double mean; double std; } normal_params; +normal_params algebra_sum_normals(normal_params a, normal_params b); typedef struct lognormal_params_t { double logmean; double logstd; } lognormal_params; - - -normal_params algebra_sum_normals(normal_params a, normal_params b); -normal_params algebra_shift_normal(normal_params a, double shift); - - lognormal_params algebra_product_lognormals(lognormal_params a, lognormal_params b); -lognormal_params algebra_scale_lognormal(lognormal_params a, double k); #endif