From 39652b3bfe328e1cd41cdcef523ed21716c2d01a Mon Sep 17 00:00:00 2001 From: NunoSempere Date: Mon, 26 Jun 2023 18:05:37 +0100 Subject: [PATCH] add squiggle_c beginnings. --- C/squiggle_c/example_simple/example | Bin 0 -> 17696 bytes C/squiggle_c/example_simple/example.c | 62 +++++++++++++++ C/squiggle_c/example_simple/makefile | 53 +++++++++++++ C/squiggle_c/squiggle.h | 109 ++++++++++++++++++++++++++ C/squiggle_c/to-do.md | 9 +++ 5 files changed, 233 insertions(+) create mode 100755 C/squiggle_c/example_simple/example create mode 100644 C/squiggle_c/example_simple/example.c create mode 100644 C/squiggle_c/example_simple/makefile create mode 100644 C/squiggle_c/squiggle.h create mode 100644 C/squiggle_c/to-do.md diff --git a/C/squiggle_c/example_simple/example b/C/squiggle_c/example_simple/example new file mode 100755 index 0000000000000000000000000000000000000000..b3db352fd1b79c50dd3ef781efa445004e27efc9 GIT binary patch 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); + + 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 new file mode 100644 index 00000000..5409e01c --- /dev/null +++ b/C/squiggle_c/squiggle.h @@ -0,0 +1,109 @@ +#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 new file mode 100644 index 00000000..3be7ffb1 --- /dev/null +++ b/C/squiggle_c/to-do.md @@ -0,0 +1,9 @@ + +- [ ] 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