From 2a39a04c690e8854bde55432e1374dcf7710a80d Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Sat, 18 Nov 2023 20:10:29 +0000 Subject: [PATCH] add desiderata split_array_refactor --- C/alt/04-factor-out-paralellization/why.md | 1 + C/alt/05-refactor-split-array/desiderata.md | 37 +++ C/alt/05-refactor-split-array/makefile | 106 ++++++++ C/alt/05-refactor-split-array/out/samples | Bin 0 -> 22560 bytes C/alt/05-refactor-split-array/samples.c | 260 ++++++++++++++++++++ 5 files changed, 404 insertions(+) create mode 100644 C/alt/04-factor-out-paralellization/why.md create mode 100644 C/alt/05-refactor-split-array/desiderata.md create mode 100644 C/alt/05-refactor-split-array/makefile create mode 100755 C/alt/05-refactor-split-array/out/samples create mode 100644 C/alt/05-refactor-split-array/samples.c diff --git a/C/alt/04-factor-out-paralellization/why.md b/C/alt/04-factor-out-paralellization/why.md new file mode 100644 index 00000000..e477e1b4 --- /dev/null +++ b/C/alt/04-factor-out-paralellization/why.md @@ -0,0 +1 @@ +So that the mixture distribution could be composable with other functions in squiggle.c diff --git a/C/alt/05-refactor-split-array/desiderata.md b/C/alt/05-refactor-split-array/desiderata.md new file mode 100644 index 00000000..45c300e7 --- /dev/null +++ b/C/alt/05-refactor-split-array/desiderata.md @@ -0,0 +1,37 @@ +Instead of + +```C +#pragma omp parallel private(i, sample_index, split_array_length) +{ + #pragma omp for + for (i = 0; i < n_threads; i++) { + split_array_length = split_array_get_length(i, N, n_threads); + for (int j = 0; j < split_array_length; j++) { + results[i][j] = sampler(seeds[i]); + } + } +} + +``` + +Algo como + +```C +#pragma omp parallel private(i, sample_index, bounds) +{ + #pragma omp for + for (i = 0; i < n_threads; i++) { + int bounds[2] = split_array_get_bounds(i, N, n_threads); + for (int j = bound[0]; j < bounds[1] + 1; j++) { + // o j < bounds[1], no se si el +1 va a ser más elegante + // dentro o fuera, aunque algo me dice que dentro + results[j] = sampler(seeds[i]); + } + } +} + +``` + +Por qué? Porque esto hace que la estructura subyacente sea un solo array, +lo cual implica que no *tengo* que utilizar split_array_functions especializadas +si no quiero diff --git a/C/alt/05-refactor-split-array/makefile b/C/alt/05-refactor-split-array/makefile new file mode 100644 index 00000000..149b8827 --- /dev/null +++ b/C/alt/05-refactor-split-array/makefile @@ -0,0 +1,106 @@ +# Interface: +# make +# make build +# make format +# make run + +# Compiler +CC=gcc +# CC=tcc # <= faster compilation + +# Main file +SRC=samples.c +OUTPUT=out/samples + +SRC_ONE_THREAD=./samples-one-thread.c +OUTPUT_ONE_THREAD=out/samples-one-thread + +## Dependencies +# Has no dependencies +MATH=-lm + +## Flags +DEBUG= #'-g' +STANDARD=-std=c99 +WARNINGS=-Wall +OPTIMIZED=-O3 #-O3 actually gives better performance than -Ofast, at least for this version +OPENMP=-fopenmp + +## Formatter +STYLE_BLUEPRINT=webkit +FORMATTER=clang-format -i -style=$(STYLE_BLUEPRINT) + +## make build +build: $(SRC) + $(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(OPENMP) $(MATH) -o $(OUTPUT) + +static: + $(CC) $(OPTIMIZED) $(DEBUG) $(SRC) $(OPENMP) $(MATH) -o $(OUTPUT) + +format: $(SRC) + $(FORMATTER) $(SRC) + +run: $(SRC) $(OUTPUT) + OMP_NUM_THREADS=1 ./$(OUTPUT) && echo + +multi: + OMP_NUM_THREADS=1 ./$(OUTPUT) && echo + OMP_NUM_THREADS=2 ./$(OUTPUT) && echo + OMP_NUM_THREADS=4 ./$(OUTPUT) && echo + OMP_NUM_THREADS=8 ./$(OUTPUT) && echo + OMP_NUM_THREADS=16 ./$(OUTPUT) && echo + +## Timing + +time-linux: + @echo "Requires /bin/time, found on GNU/Linux systems" && echo + + @echo "Running 100x and taking avg time: OMP_NUM_THREADS=1 $(OUTPUT)" + @t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do OMP_NUM_THREADS=1 $(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 + + @echo "Running 100x and taking avg time: OMP_NUM_THREADS=2 $(OUTPUT)" + @t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do OMP_NUM_THREADS=2 $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 2 threads: |" | sed 's|$$|ms|' && echo + + @echo "Running 100x and taking avg time: OMP_NUM_THREADS=4 $(OUTPUT)" + @t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do OMP_NUM_THREADS=4 $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time for 4 threads: |" | sed 's|$$|ms|' && echo + + @echo "Running 100x and taking avg time: OMP_NUM_THREADS=8 $(OUTPUT)" + @t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do OMP_NUM_THREADS=8 $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 8 threads: |" | sed 's|$$|ms|' && echo + + @echo "Running 100x and taking avg time: OMP_NUM_THREADS=16 $(OUTPUT)" + @t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do OMP_NUM_THREADS=16 $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 16 threads: |" | sed 's|$$|ms|' && echo + +time-linux-fastest: + @echo "Running 100x and taking avg time: OMP_NUM_THREADS=16 $(OUTPUT)" + @t=$$(/usr/bin/time -f "%e" -p bash -c 'for i in {1..100}; do OMP_NUM_THREADS=16 $(OUTPUT); done' 2>&1 >/dev/null | grep real | awk '{print $$2}' ); echo "scale=2; 1000 * $$t / 100" | bc | sed "s|^|Time using 16 threads: |" | sed 's|$$|ms|' && echo + +time-linux-simple: + @echo "Requires /bin/time, found on GNU/Linux systems" && echo + OMP_NUM_THREADS=1 /bin/time -f "Time: %es" ./$(OUTPUT) && echo + OMP_NUM_THREADS=2 /bin/time -f "Time: %es" ./$(OUTPUT) && echo + OMP_NUM_THREADS=4 /bin/time -f "Time: %es" ./$(OUTPUT) && echo + OMP_NUM_THREADS=8 /bin/time -f "Time: %es" ./$(OUTPUT) && echo + OMP_NUM_THREADS=16 /bin/time -f "Time: %es" ./$(OUTPUT) && 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) $(OPENMP) $(MATH) -o $(OUTPUT) + # ./$(OUTPUT) + # gprof: + # gprof $(OUTPUT) gmon.out > analysis.txt + # rm gmon.out + # vim analysis.txt + # rm analysis.txt + # perf: + OMP_NUM_THREADS=16 sudo perf record $(OUTPUT) + sudo perf report + rm perf.data + + +## Install +debian-install-dependencies: + sudo apt-get install libomp-dev + diff --git a/C/alt/05-refactor-split-array/out/samples b/C/alt/05-refactor-split-array/out/samples new file mode 100755 index 0000000000000000000000000000000000000000..7157dc0bada9100dee55a618703db13977e02f3a GIT binary patch literal 22560 zcmeHPe{|HvwV(YFSS0Lssfmq0TK3tcja5Px(PkAS3mf>ZT_8$$Mn#vHY)IPt>SkAh z&lhP^DBp%feCM>qK6|w9^}N#7_lmDSeQGtD*JKvej zcQ>21&*S^6nR9mU+&lN)nYnZ4&U|P0I}elvmYPi_CXu0@Jx6-U+0rhmZ$0~*2`H2EF=RUNobU#qT+?Tm=!Nihl9c{aRH?qp(Ix-nKEW>`)j?_1@kq_bsjpnfZ4q|NgjO;* zQVQxQy9IK>)PnMtfgI*rMST;EL%*=6gra&jx1s)?>*vpHsF~AH-xO_|(^hi*oa^To zwKf;c;}uiisSfg~Wh+)QsOch_Fe4v}Jb_dBAKv)X`^RTL`tbT6bUuGNoV@wT*FSnu z>YzG_hZ4z)CUKnd7vn=bdjG%9=|rP65PGgr&t?C<1cG8Jr!l+RThBaK-W^4b4rVimH1WLM#|u*VxiZfBf+4~nw7Dkc^zy0S~yb6kgH{(ww7AfTHi#~tZQy;p_<$bnc%umB-mKh7L3$| zLsd1cDVmN4o1%?u*~;bR!Ir9URYOClA;o}_ln}I5i^d0A!u3s&+F*6vSJ7HDwEC-* z57t)IH?Wpyq`Iyu%o?Gpxfor!MbcXY z-BwOgW9x*(kF3<`6pLi6Bg&wGev!_?zvBjau7N&mpwBeW&l>15#!+#dpwk!#lW?gP z1Kk)m+YEGyvodWo(D{Bq6l^ij^*E~+Y&X!UZkcu%=+bsn_LzY#wjx|)r-4o}R;FDB zItDP6b{pvS6cvk%fvyu?zb*rxlY!69 zz(13L=V#zEGVq)X{5*~WjO}?Z1Aiw2Kc0cVl7YWq#Jv@(;R1yiWQ0>H1rHd8j*Y>E|zW6B^UON%c-c(}?S78w2uVO3bKHxs2 za?(cl#Xg^lZ%4ESS1$54k`8DtSG*P4o+iiHO2LOzn7MxV$0| z%Xt%9A2l|knS+Hks7&S`Tfo@q?P&Ql*g}ovmMyhX})Y+NX0CAHO#&+ zOj%Vs9?4ZZ7v2jbzTtc_8LvST`*IqAo9=>as(zABtVf|}58<`SsV<{_STU#W3f^Y! zzutfJuErRzQninKo$GSrJ5hnQH+~C&Zl5;z@cZXI1IpGMOk3sW`=RZSX}}lDwfQ^; zl-M>5Pv;WEhZ`UEzJCsUv7dnA^?;+&q!hm7RaOj?X|I%NFMGS+F?&1SIH$d!RuA~Q zk6Qz=yg!v~I;w1Y5E9Ef%^hdW%H~ew16t3QWAh7@!ryOK9A>qn2SZ-DUnv|^`UY1c zdyumEz`BxzQK9M?iT=^29V_#^9=Ter9vSs2mAShDvH8EE(ed=RU((TI9$lhT28yj- z&mY^T>Ump7w+Z3{{*E4-C!uU@gyp>G*0v`hrS0`kU_CeG@ng_J3>tOFqw$$Q(oa^LYDXVUwrtpYvrG99DEl z?}5XniY`#I9}J=WovR#b?8}Z?rRjv$kItQ_RoYJuD^IvRJ<%0+VlH^^^4=}|0V9;i zq307X?>n+8tANXIjFvlxXEvNTh^*t8e9dUK5%X7YIl{n z&xraOQ}w;)@95+8m9gcX3(+?+Dx|6&c$y-Z95ipRSGWDM6wA0in%~|E#IB3aAsEN5 zqw@m4c5Zd8sx);!O`#P7pr(cH5_cEE*(ipCnEH|b6NGA&h91I?C|4oGQmzlAUc#L5 zUw?o>>&$a0g*kT?P&nF)7&7A;WR;mr>!$$IJvhVNqxE;kXL2x@<+Wg@)x;%GFTLM!6d(m%E#`1B4c<+Nt=zo<+Z4>yojkPHY$q zSVpDD&h4<9_7qqJk`;RpEM}s3s9Tv;%G_N;34(kT8G1Fo73vHg%hv-~sGxHyI&pwH z0_FKD&`7PP``8S3KfsBZ930CbK%MHMk{)dUee3B?%yd7a4RpVj<9>$9p3y!%s@wrB zw?Bh+1YEYb2UU*13T`N+N=8tu6dqBE_KzyJ90k)FSMK_yyiYc=){B?WLWiYoEwxoC zq=*ZBcc4#D(@o;U_G`{one*R(FhxP)q4%*aoWs}9FUF=M`N#PFH9_VR>)-9DlKOWG zl>ROK`ygNRZzBXJ?%%7Ntup5&K;FM6&Y*u+Wa;0gGwpTW)xKC$(u>7>_HqcU#-?PI z`zU&D)8~oy5?K*jk+D4mne+!NegbM$?TFgZYw~;cMh~Md+tai|)Ms7dizpwnCof4{ zLbTNU+jK_RwiN=grAfYoX~=6Ms&^@;uyN_8_(}$d2p@7{*>shx_IH2?m(0bNG z*?eFJXg%v!3cHD@ditaN+Hl9~KkNChRC%<|G>oeB^@IFv(;;QskH`kYv=aEW3$(rPYR}UK0lN&$ zHU#b5O^4d=^Ln25#oA1M?L4*#&o8e|x{vz1-!f}2z__};fR?0zLHv!o-%BR_TG?68 z(ESOo=jHaR0@c_ws?`uVj;2fgQ*j$-w7yD13dRS9$bAekl*fgmLtO-oog-Y2rxZa~ z!L6>ukYZP)#v1ksKf*o-YHe}rX?TriD0=En1RX@7#ONsPxwt%fSm-{G7IU#(rZvcY zNQ<~^wENW-xUf}hhabruaQy2~$)kWHV;7r)_80cX2T%i#SDjP!xtktt`?aaP`u>*3 zvaYQ2n@2AVqkntOSFqXb75YSwnOKAbJP_cb-MwS;A$xqW6YD1vKXcP&~BH;N4>Uxcv@|Dr@#Jf`y>6M_-4=?12X5PvK3$Y6R~&g%RN^0(K&&4txS@ORg*a>btsW4r8I8 zz%-@oy>_u_E0b9w~VtpNa4S6hoXx35w z6nxC86!xm75wf%5aENt~?5FfSL+c(D-XzJnoYHPvdL}jJ@tNLo>(@YDXWsaJn6vMn~XRwC%n3?SS@vGH;Z(Gncl& zJKhls18sR7av{^>4H*|!&DKNH$L?XrIg|1rJr$z;%BBs-P%lz_C>DK4d`T}k>c)Py zHxa@^iS8f#`45Y_T8(ED+KSQn5H0+wO6ymd;*Q^K5vQ4avv4@xqvnj&`$cax4g-`KhuUs^}G^2fkxZLT5R0oQ2eHoTM@@8 z#*33&YT;y>w7{eVCM_^&fk_KYT42%wlNR`Uw*dW*u|5)NymNNVUCxEh*|p-&d#j?2 zg{>_O^^sszI9#o4SchRZ)KfEdLH&3ICE!zI=TT7WM;wPSG4+a*ktF0v&Dx=h8y4-Rj{ik;pkR5+c zCh32WTLzQKwSYSS+W^~MOD6O1GqZXunREg#c^wZVfGXgR04o9MtaIn_Wb!Ov!CQDt z!0!f?@no_Euq2U8ZUrpBHhU-F4!}16)gkZ$ZUM~2j}toqoq%=N!bb=Pd<@Wr$BW+r z76SGIRsy~ONY8I)0qgKUmWSt=oq&aaTkuHj1FXZt!dk-d*z!#($GO#hz#}*iET#V7 zJaY^E2Lg%yGfv5L->uBl<}h8JJJm*gBszZ5Wlq>&E;0C~;%tU77e4gQ)yh|s$zdYc z9ZT&6H!IWE+cvUAm)-EiYh70ondH~t9>J#o^g%$Gh-WK4kD)&Ldj%z; zJ%-OCpvnI#E&%>Ze0G4g7&c@}?T&5c#rA@jWwG74)#|kuZqD)Ai#w)Rj+-vB7kceZ z6fU+q7Tawq{+FRb6p6AzSdF-JtZJP zHCR&Gywm9+tQ0y9zL8AM7rM5ZZ?zX}wydzPHNSX~y}-*WUOJuOr!7Wj8H4bZ!^kt$ z?lkvYWVe;+Wm{18qc_Kwm4KCetq#6*AIaTpcf`yv`Bn=wY_|IBg&jFd?Zpv$;nH5a z)7yv69z>&gU>$q2eKEqQ5evP~jBZWN<-Lyrha zq-`=yT42%wlNOk?z@!EKt}P(no09KI(GpK76V~bE6#XfIUJJePSRe<&z)o*fmSm+{SY=kamEE6}o|E!(gL-Dcdz8|@)9F*0oCoQ2 z8C>6?^K?39y&NAyfd#YQnU!v5jx6J7VG46*r9%n+k&~G$ za0VGyjAa&LE;ds1yajKX%lGk(r52nLY;49cK*$&{a6fJR}M=*9T6BPD(& z0Z7lE+4L+0o%s1n`pCLE6QZ_!-&c^?A!M?(f1-Z9#`WYemwZqciYC%S#o6?q;&S=p{Q(^`VRgv#hiSay zF3>?AAb88U9%jXgM#i93no8H0n6bD^p6OB zhiJcCZ>OQ1WS?Q7hx~)mO5!({=&`~Nv$coT zMW>K8f-i9XOPO4cW!Ze98{#&OXn0xzI``+xbzGhymJ9l0Vgts{L%|Z{be-4u9Txm+ z#k>fK!f#B#e=+C|#H)ggIQ(6qr%P;P4+}lN7j~2P3dwQ&H^D#b)EW8N9I}rK`gSqy zav%7ls5f(7^a_4CnfQ4b^w2&a+x$8v_>1WW2c#o9V8fu3ofj1A1Zn3{$uHW&&p^Rq z0R!3J;newG)&Zm60kYLQTkt<7HuC%o88WvC{f;Yj1^oOGS-0Tt7xO|sFWdn7O!N5o zu$0rs`~Mvi(3?Sb>Q$zFYoKRq{|BO-PT_Cz9Q%iyA8|#b&Ex1agH_QsRt>EUp-8BvsHFIMcXlSaYNTHL;1z0$gg4-Vk*dZ} zuqN8rxB(&tEC^2d+nCPUTwB}Fj2lu!(}KaJw|bYC1bsQ%mjzZX_6CA0mo8mZRuQc5E)JBDR(>H$b!(KjoL`4hT57!h zL|mwnPBa>v-=OlBZj_;$VGQD5V!=>NRiuiAtD0(p)Z++?hJ#JbVcgl0y8C79a+uOm zc`J*tR=QNCA=I=kQpfc)HwL3k^|jD}E~Q@yGgkHdly<{)G-K)wz9M`pSgbd)I2dUT zy6K*qF=g{Ss7Z`i@ZKu;NT{iXUxj1fL0>a(oH3N*av69UOrc*$ljg?!`kJapDAi=c zWjJF^9&15uw5eKeLgocN!J5|QU>$5tm;3lvqOO|yreL%+RKxTiW~-x(m^8A>jrDDj zXgH)_?UUKql(UfUiaU1FdGtw%`+k}>$bpG8vo_43y87D4we!&9dsG9dZy1tniF-f0}`)$$|mnB6x8B(3lFkyhraLzc5I4kY+A!JW5R~!x#o5;tC{5 zBdbax-IneD_tfeOi>$Jj~b2F5T!Sf6}4_?M8$xSuny~F4)+{|!Y!<* zsW}oVTGtdUY6&;Dgu;;x2GTvzdfYKqUn5A~#r`>ws&$OB*HyLFv7(v{O{hYLk+5EJ zZz$Ya-`tdj1W^_aHB^y+$h9;?SP>uWB4mr!H3OwgYp9wPMM7;rDdK^P@RZz zbv4i?Fe{Gbn6!H?fB%^>>?$VMUm{DGyV@T*sFM>(DGNnD=CNVrSLQ~Rl%vi*Y0V#tt9 zrMx`1k+4fJ%KD|e?Em+K{9?f;&vzu0=bj`_GMVx}1dVJW^<%53CjrGaN)X47|05_d z$~%Regu^MxG_|Q}ew-mM&y`Z}ugLGpkay111tc8&EAs!7A)ne*ip}K<8RsLSdnv6? zXUNO*ED7a)SPEp$pS?m}j(_Uht5}!VA_WD^tpC{zd3io2p*6$dGwVN^AurF-B;1i9 zApJ)+>o_uarHA`A(d4;~yw5U|U*ZzJ3mL^AFVE#l7Lv9xf7P)uT*{wE0r{>7d#Upp zc|WN`7?$QAB`GKE_%qPRL{eU!^OY|mZpL=$AoVTfC7c3bqkQUlyi^yF^+*L$PU5)` zCYr2YJ~ui_b&gDas>_&Wiu#LPU(B0frtc=iwvUBjXR9 lgUGh)68d-B%{rw2kr|V$SJG&@W|J>nsdHSGA&`Nx{{?iVj!FOk literal 0 HcmV?d00001 diff --git a/C/alt/05-refactor-split-array/samples.c b/C/alt/05-refactor-split-array/samples.c new file mode 100644 index 00000000..f07468fb --- /dev/null +++ b/C/alt/05-refactor-split-array/samples.c @@ -0,0 +1,260 @@ +#include +#include +#include +#include +#include + +const float PI = 3.14159265358979323846; + +#define N 1000000 + +//Array helpers + +void array_print(float* array, int length) +{ + for (int i = 0; i < length; i++) { + printf("item[%d] = %f\n", i, array[i]); + } + printf("\n"); +} + +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]; + } +} + +// Split array helpers +int split_array_get_length(int index, int total_length, int n_threads) +{ + return (total_length % n_threads > index ? total_length / n_threads + 1 : total_length / n_threads); +} + +void split_array_allocate(float** meta_array, int length, int divide_into) +{ + int split_array_length; + + for (int i = 0; i < divide_into; i++) { + split_array_length = split_array_get_length(i, length, divide_into); + meta_array[i] = malloc(split_array_length * sizeof(float)); + } +} + +void split_array_free(float** meta_array, int divided_into) +{ + for (int i = 0; i < divided_into; i++) { + free(meta_array[i]); + } + free(meta_array); +} + +float split_array_sum(float** meta_array, int length, int divided_into) +{ + int i; + float output = 0; + +#pragma omp parallel for reduction(+ \ + : output) + for (int i = 0; i < divided_into; i++) { + float own_partial_sum = 0; + int split_array_length = split_array_get_length(i, length, divided_into); + for (int j = 0; j < split_array_length; j++) { + own_partial_sum += meta_array[i][j]; + } + output += own_partial_sum; + } + return output; +} + +// 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); + /* + uint32_t x = *seed; + x ^= x << 13; + x ^= x >> 17; + x ^= x << 5; + return ((float)(*seed = x))/((float) UINT32_MAX); + */ + // previously: + // ((float)rand_r(seed) / (float)RAND_MAX) + // and before that: rand, but it wasn't thread-safe. + // See: for why to use rand_r: + // rand() is not thread-safe, as it relies on (shared) hidden seed. +} + +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); +} + +// Mixture function + +float mixture(float (*samplers[])(uint32_t*), float* weights, int n_dists, uint32_t* seed) +{ + + // You can see a slightly simpler version of this function in the git history + // or in alt/C-02-better-algorithm-one-thread/ + float sum_weights = array_sum(weights, n_dists); + float* cumsummed_normalized_weights = 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; + } + + //create var holders + float p1, result; + int sample_index, i, own_length; + p1 = random_uniform(0, 1, seed); + for (int i = 0; i < n_dists; i++) { + if (p1 < cumsummed_normalized_weights[i]) { + result = samplers[i](seed); + break; + } + } + free(cumsummed_normalized_weights); + return result; +} + +// Parallization function +void paralellize(float (*sampler)(uint32_t* seed), float** results, int n_threads){ + + int sample_index, i, split_array_length; + uint32_t** seeds = malloc(n_threads * sizeof(uint32_t*)); + for (uint32_t i = 0; i < n_threads; i++) { + seeds[i] = malloc(sizeof(uint32_t)); + *seeds[i] = i + 1; // xorshift can't start with 0 + } + + #pragma omp parallel private(i, sample_index, split_array_length) + { + #pragma omp for + for (i = 0; i < n_threads; i++) { + split_array_length = split_array_get_length(i, N, n_threads); + for (int j = 0; j < split_array_length; j++) { + results[i][j] = sampler(seeds[i]); + } + } + } + + for (uint32_t i = 0; i < n_threads; i++) { + free(seeds[i]); + } + free(seeds); +} + +// Functions used for the BOTEC. +// Their type has to be the same, as we will be passing them around. + +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); +} + +float sample_mixture(uint32_t* seed){ + float p_a, p_b, p_c; + + // Initialize variables + p_a = 0.8; + p_b = 0.5; + p_c = p_a * p_b; + + // Generate mixture + 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 }; + + return mixture(samplers, weights, n_dists, seed); +} + +int main() +{ + int n_threads = omp_get_max_threads(); + // printf("Max threads: %d\n", n_threads); + // omp_set_num_threads(n_threads); + float** split_array_results = malloc(n_threads * sizeof(float*)); + split_array_allocate(split_array_results, N, n_threads); + + paralellize(sample_mixture, split_array_results, n_threads); + printf("Sum(split_array_results, N)/N = %f\n", split_array_sum(split_array_results, N, n_threads) / N); + + split_array_free(split_array_results, n_threads); + return 0; +}