From 7e46951e9547d17bf6470dece6a0b0d5ff3fc478 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 3 Aug 2012 14:48:34 -0400 Subject: [PATCH] Added material on fission bank algorithm from NSE paper to documentation. --- docs/img/master-slave.png | Bin 0 -> 6863 bytes docs/img/nearest-neighbor-example.png | Bin 0 -> 34317 bytes docs/img/nearest-neighbor.png | Bin 0 -> 4131 bytes docs/source/methods/criticality.rst | 8 +- docs/source/methods/parallelization.rst | 645 ++++++++++++++++++++++++ docs/source/methods/statistics.rst | 2 + 6 files changed, 654 insertions(+), 1 deletion(-) create mode 100644 docs/img/master-slave.png create mode 100644 docs/img/nearest-neighbor-example.png create mode 100644 docs/img/nearest-neighbor.png diff --git a/docs/img/master-slave.png b/docs/img/master-slave.png new file mode 100644 index 0000000000000000000000000000000000000000..01d0cbe977a81681b025b3e8fa1df7afd1a850a9 GIT binary patch literal 6863 zcmdUTRZtv2udSf;$8Y1P$)6I0RYT2^w621$TE1EWzEK;2L0Ial71C z?|r96q~ci_0+Mj`n5WB2#1FF6Wb7hZ z-mpbR@XFwR7Vptysne3Ktzl+1mEa$w&8?u&#v@81foDy5Z7ktD3H&IPxX(nYY;t={S3 z-D1>fgn+tK9Ue9ji*y`LE)n5pwj2Y?9`(RaK4$6p-LQj&Oh@+njP2tw0|*g10LXHW zg8dTs;(b>4MIZ>iSTbl=r0VhoqLoIIklwV4MwL1OC#0AA!zlf%cHXV{j8Z6PSX?2t!5^M{uJ-Bn)JoL~8qunuhSn7y}r9aSrIBL9Iaa!9~Om z1egcv*I*F_(AJ=wplIe2^9C?ZqN@a}I8o^%EcDnrkzfWhiK8<9W}|s?)+3B0p%Z~Q zMMEZ!rXKD_Ll=RiD~0`oVh8K3#4k18BqR&*)=BpgRMyb4T%i)Q`?o?U27wNe!IN~> zm>FTuVWK^#NRup>}DA?Y&*+dIIDsp{{>@Jo#V4L7Dzy-=h0aP=$q_vw0FD zpyHX(3`cPE%C~Zel^-z5?-?Raq|*ibG;74-FSsC3<5ah`14k#7ABzwXmD># zWa#f2kR_>=5-uZlJ)TCiV1x}+zaP`IpQ~PnT?@H9xUm;v46Qj_d}GIOjdYFgj`}h1 zDS&P($w^i~P9Ma9mKx4Bz_-n@73P%K;?*HwNa+}O^IPsO|KCsNSho(3l#fV$*)$Ty zconE!Ap*T0B;Wo@c8UQhdVZqE*NoVc<}_t@W~&YIjPwi^kbRWTOK4EgS4U+DZusdI z?UoLYfTt)if`8Z#_3XG^IJgkpvEGs01(PZ(sYWvqziX1I(5sNLlRx2mtNOXviTRjW zFv(?@Y?yL5TEI|WDwB;oQ6C}Er@34qHM)|! zS-RQsS!Z%djY4M|pJjJF#DC=WL9xH;4}>AW&IrE*QE?Ilb%q#*Cx&{4at2O$bTU+_ z;2tx2kRoX@84hVq0bzlFdc1mjfnH%)!Ct{`fqfxeAx7a!L3p9q5Xb^(zGZPe>^3Z% zEHh*~0vWO#qDtsw{75I4TA#E%LNuI};?0gT0v_=maZD1U-QW}}&n@{|@FP#H)JfAs zHvX4&Icc#@k>_V;z3)Zlm4xMAirf?q88i}xOc+6IX`M>(bscr}E@3WN7XcR?=%4J=^Lv#M*TyK0YUf#4|{EttX9U84snnU;9O=#S5i3dPd4(I{wvuInN}B1i?tb$iYluU!+K&!cYfj zcT7Xft?Zoajdp^zP}OQ(9e9_*jB;CkTdE&N*_!c`@uMIl~QOPme zQQ1YeCZ_f@UntPt!K&-eW5=BwhDMZI)C{q|?vY;q*AHJ+DoraR_2+7AR~7%Tj+Tup zT8}kgH;UT5Z}{AVR$toef;7+R;|q)=WG|ZhTMk%jdM2Do^16D&aU^Dpe8P?9%LRh)IWVXIZDACk135a&tH7 zGylN*&~ejv`5xxGReYGbxv;GY!+|+o{1Nm61*{9MUs%^?YUeOn3%yFZGC7_I$&*5q{3%7!E8KeowS;+J zruqlGu=ak)V>cC=X079O;T~o!}HU|HZKCR-)fo8m>KaLOYCupj9&(^L3{;l1i-+jFt@DTWC&4SB@OBj99u&Vn& za8$gQgX~_ps@Zz{$wYdyk~TQ+A}`d`4_29jpJKgdIvbsYzNfzR_2wfSGuU~N@^<2IW-!5`x@gj zR+Z5jjE54Ay+mbOsaf*A?8!~w`Wz}$+!hOoxxb)Gqf3`Fx+sk@AeXiTNv^!?xhF6w=?}yb#FXsns4M&SJP8o`;1(fdCKPHQVGz8pKYr!*)xG(p?B=-ka(KU-Q{IrdRlX(7 z-x;PpUbr;El4?~YZA5mayK%W=`*-f4)36h9LBD#f<<#@Zmb&*dxlf z{IcorU#wr^^Yw<3AD9|MdC+3!UVL12w|3g28B7vx7gn`e-|g`vdiMNpa}FAe^GGEO zL?go!jrJ6{V&0tkIglQ=pJPWU6Gtpc)vi(d6M=e&s%@RN8L-k%d4esyO$%9w~>tN zsO#*{29H{M|7tz@&mfi}Rym!jK4N0pw-S3$SlmgRZ_XQ0<#DU`#w4twVBgc{t$PbM zE@x+$7)RcZ4i}c$WB>qxPEF~H%*)FQ!vD5mZ_lo7AK{eJo<90QK|w*s#L6qG5Ye!y zdww1D59tn|NGo3w7}3S`SqX?5&Z4#N-r-~y}cS!Q#zrc zqMMrr1qBN3?t(fxoK{x6tgJNW=ca;!4E6OY3JR>r$afk#h7Z=DpI7 z7xp~`RbtD0^@GjNf~Kumh1CX|pZWHg)(r6}(u(V8G<_Nc&o9*nR4vSCNDTOk2z+Ok zi%HhHw>6h0?(QTQ^R^G)uP-?ZVf97#cJU~`o?LXuyK@E0-ew>UA&mu^2knV-OSsac z=L-H%ngq%-T4U1_@S-H3{r8p=?6D1Tij(-N*H zT;Ne>rSaQi!NXsI*^J=PDKR>dldF;qraX+po7f3yewV4kdVWqE&NzkYjs53?Y|5Ly z?yPIhcB@WT-T}j!aiHT+I~4lKg=KW=o$nsL%TU)z;XbxzuL0;nDc^FK}y^*+NK zt0k%68Q7eHZEh}2g@OhS(O!;A23#Ue98|U(c9Y3yOodK9rB}7 zOqAnf53%3atYbS(Dh09Q;vk(;HnBV|qnMNXdo2fWR6$58x15nntKVbLUJ2>(7)&Jz zcNIv&&4B899<+TCk;%9x5c5osVV0xg7#rpHBAp>6`?M7y?V8sy83B&P&M>3ab;SNY z3Y1%eCW=-Is!ZtS2?$yqZUw>H)-)W_vDqZnvwsO*M=6i>Oe5xdaLCdJsoXL?)Xrr4 zV2srong7hEAi*!rV>@&gN+n@3NTF5PsGWvO&r=!PnK>?f^bz)+5_!!GR*QMu|1kUi zuXg^T0BD!2<2L!beP{hOe5>Tqtdr(7!Q#ew{^ExyZMHv1qyH9&k5mr z&MMnr&!|@-A%rVS0p}EPty_5eH7nl+hRSOE?R%Rb$v?Aq?j#&VFW?p$W-EdU%zzWy z&qd-g&pLV|4%Ews_znPPwJOmNg68a2%LJEo5@_-6!zcVm1h;SeZ&51)3qaN*Pq%Rs z&(bfec3bSE3H`hI+RKyR*nxpY0UZg`<=J`g`{kYBWr%X~6s3H~Ir6dY;_cj)Z{+E( zkGMljE?ujHQ0<9w*zSWWtO>iI|2SsiHY-0*dvgP3(R^>T<1v^h>mi-f>xfUtnbPrj zvaU1`ahXYF4f~9|F|S=tD)9ih(HZ+mVY?q}=g`0=5mh9pZ)>`6;Wi~~8=g{=(2gjg z(S12FNU+@+uZSh<@{pP``lGII%XAk8pIjqBiD@aheWdy3ALoj@k1-n6BU#f6MA5x50;8kY4(2g5^n&hH&JM%x9JY(Le`)96lJD~g+G;^_<<*F0#3<4%M-swL*ItbueKDNLmZ>y0TAQLP)AR(>>{1DhK~Pl?OfJF)f}mkm!Jkv~?L zDMpc2g_)7*SZ5pyZHjvRG=jF;7>?TT@gHIJJ%+{qB$n<y&w5!ba7uv^IvA+k7jX z@^phTuA%c{OvM{BUzA*2)!3qdA3x2g(?T53zmmb%L?ZSsu(|1Ba(@2y#A`l8j!J6= zkM9&ZzfmzcGIZ)&)lq6vT~G9F9bEI~@@s}m2-8IBa+)P*x!OE-v1q@kZ|80x#y|=3_1^+#nkvQEdk$X%b4fbi4D2(? zTG3VJv`AmnN&@ywmQGK<`S@N?tb3rn8P z@3$->H`NCH=kC6}EBRM&m+ya2e6I7uC~$t+&s3v#AVO`IUy9uEK2}vG_r#JGh=sO= zsG7=Z$RLjp{CV$t`85Kg`)URsg@wcE+x&h=Vpme5HrVJRg@Nm?khGxH67Sjih3245 zE&R)4f6GqBmTtjqJu5*aiEH~&MSgN!LfLI(C28~BPPLdyPkZ*yD*z> z$`P@N5_MgBD{&oazxR*Cu4anFN{D{{L#TAM@s$-&uQ)W3zQ3a6tf}c@&&{p zeRv`t4fvD^q`?OtA{y|-QuO#>W_3g(PJoa?hu|9KwA=tLmU6#YW#Uh>?EGse(QbxO zVi{5*$}*pR;9XHG)4YqP&O>}dFl;s5bQp)eW|CgCC72_P{$st8*bhjG>ZjI@iDr^U zYuy_^9W*@?NkKX+RsuIeRAHM1YBUBREiXQMoS^cv&x|P1VEZHe#4z{ydtvYyPB+Jm zZ%BDAq5v9J4^AG)0H+G_44f%Hr1>G|JAO-iN-8UytuW`O!*un5aoZwL z^sY_l)ggRsOe(No4vBB{%xumb1?KP~1RP6+6UJ$=%QhzHC1}``O10jW&?s_JsBD$Y z3Co)TllKB@u!o~K&?0IzzDfvi*D8hX1#mS_Vq&Pz3DGnDF8XYOjhlu3;zTI$oz8)R zl&-bJjAa8z6`>M}P!eT!Rg1oS&}|dAEbpHTjY5$WQL8*>xPn?I*?_agpP(&qY+c za*gm@K{cF=eDt1xB9yT}$Rtn`6`5DW8hbW)uS61*dnnOf*$S(oiQMp=bg8BZJ65=C&A_+>%ws1m}pJ2D4M%}}OI0O~W)@NnvyGHl%^nibATts zba;+#lncnI``{$fW1e;Z^OzHLq{_3x?FZM%!$-|TXd}T_J85PTDbC3ZL3lL{f5@{H z&WCiKk-TA6=H&X}xl=*z2}1z?AUa*#7$>RBBLY{6KG`0>-}Fbf$GC8qn%5-;ic6pxv| z(L@2&^(p)8Hc6FL|H^>s=)0@Y7AhDC0oTq34$9rrxPxtMi*JbF{6V_vpc9)b=H3an z@!Z)M&D;RaQ*DR5>|@KqYD8}m6v%dJ?gk>7yEUy*cSZV6m2d!4Q0vHcuCPjIz=;!A zoUrzQxX7{P`BpvkD2{RMLHaCwWKwFgWwNuC_XHhu_}ftFKsc7BCws+l1V4_^n2L>- z#yWaetl$xQr`oEmQ4l$A zA3d0hF_FT^65Acl24lkRFSupX(+J5Pvy?9uJ|Pgy zR$dNYK*%Fr{4!y=N1C8ZikbaBFM7O6JzgY@l0PGQ#Q@~vx0saeT(U2Fgf_hJO07Ar z25rL)3ZYJI9FE2{wb^5%ct@L^8S+|*Y7+-05SLMB2JnUqA*L5d>-u(UR}Wr3`f=fM zRq7u->E@&jo&~m7UHt1EEop*2%3PEr(w`%oi)xTAQj+Zc(Cp1fDzPaZ2A%~Rp*lTI$$wlL~Q29GGl9Gt(p;aOJi$=q1_~rE`=%;z!&Af}#j?L~MQ`a3w zefe(8F3knR#`ndKRox4GtnCf2qjUK*!CYyDcbq;<&9(|WyiTWw^oJgL)q6dlx~`&% zZ2op2JWjVOU~1kwmy-bDFC7i3=+`hiD(HsXKD=^&q)x6Q7TUq%i)+9pk@h zLWq69-=cr9mr+h7W;Bja0MFEFE{xxi2AROGPh|!-867us6E_PXGZ%|j0pJF4alZ!% zya(}o;o=qI5flQwn!G|FP*7o})BopS?`Upg>GfYbT;9@;zB*99Qc{$C0RqA8fI#pTsK~$-n*qfx z5D1OcK}JT?!O98*0?((X*lE6vCXSl82wK2!$^6ham7~Lop=juzZYqbA0X58^x0LBe zl_kV;h@&kI4cB7DXD>EDMUOL>_T)n0Ks7gjMe@(B|F1m1abaGRPZKa2 z^Q8RnwkrS*q|*}kT-sTcYeAAz?+w97vW4As>yL&dO3EJ+w~ap{HnI+}@AWrBq;S*N zK_GWUmF0eon0S)92e{C4>P)a}>NC=VKk2Q?iYa(pB&3?KsOArlSWGnG ziVXHfANG@Q{2wSNMaa8FXfF$?#-#;~Bb?j`4tj$~BkgD-SiuMrI_O0ikI1!oPWwLakPQ2#MzsjT4c?tgk+E zH}d$ZIb2WNoOTGki!zS4=O_o3fcgb;xd%C1$8h;~ar**#BYnBx{XbROwX6r+=!|&zeo> zqzqz8`SkjQoI}}SVS+YPyF!QPBhpM}w$ze(b#ZpWH(b!uj>zD=`|x|E2)g2@Ys~nZpVUh9N>r?6H>#f)>X^J`Z+p(C z(fHytcU%0MFwdu~sI1opB|r4)U)Ez-T3KSOVR|rlxE&K7zc|)AURdMBmqo)uV~Mqo zm5#N-pJQSsBb7;C?Ox&WeU?osNUE(0ubs9tsA`?Tt-%|eb&`so%AEQ-TR9b&9;@cKD4_$81;{Du3Mler5B@Mt_rI6M2(v zlZi+zT`CHciY+hS(H4loWF3NJzH{_x4WJhT}*JB zbro=`v8A>!wUZdM@C+3CkZVK?#t4XaAtvFJFv>PW=a_hTZ}Gt{-Fn86zd_zhZ8W=- zKW&z6Hf0uTR+^$oEn-ZIqNw_lONq<-eHBsNcKdeQc83}P!IXfwn_1{&zT2`x*S7gi z-N5QchF{xj2ET?P^JOHxGQeodT<`=f8hmg0AfU3hKUTO)k}lnwaQBG^1I5 zxvU>7AAA*zaF6+rcwc=_gYX;S3w%D3Ba$uBIZ`ZgB+6$*PGnK!OcV?mOzEbGMXYoz zHPn4HhA5;+{%FoBua|!AOFAtMV6DsdXW4sRbIVZLcDh%69}8s#|@ z{1pV`xnx`AX5_|BD8bQSx8nU`_F_A$DyuJ6rB>BeU#;xcs@Adl#|9S$^#>=Igc!dw zLR7rfsx*vM|EPuT&v5ZoxwTy0>DqPtY1wJp?+>B7YrBV_NX65F1^7OKf7i{NNB)xjxJ7U7s zJ|8C=cPVGpUhO;g5PCzm8eJq_d7*Q7rAKKJSk_-x^2xgOt-k-3!h4#^VB?<@?SAE> zRlkQNdON=dC7&-GNE&jP=Icwm>fV#mM`|x_S?lXbL~+y4)|sU#ru-tpSoNoBw?eNH zs^4KuV8;8>$!tN3#jLr^_1OJ%HHE=fVSbNjS7o=0QcdW(rRPIYrPllUKcuS$=k@w( z&(#oygoli>ZE{qyyYdyq>BM&YK6MR1K77#a{9J|Q^3^d%oG5tr-eG83eBXZ`i>fLq zmDinlZ*;GKKG@%TcHMay2(t5&-K{+lvb+Hk#kta&tAk5d3k~~;BcuzApbou`JjuO> zz_sE(1M@1WnP0fHxkyY+TT}gPRvoWrTAe)Rm~RBm#oMVp?<~699wxDpDVkkHmLU?| zR4txguV)7qGMb+*t}gC$I{7v`O?uazkW8D2M2|+b9lTrKYjyMJ@&k_@ZQ?DYIN?yz zECz0$A02GZe_K@uOQueaOh%(okj#ITdb_bS{ku$MqA7WOD83+Gyf%>Tq4-#T-%K>Z zwpWE9iQx6qJ2Q!#4YkfI&XeZNieCLw{q$+&X-?ijUJ|=JyTdNi3;$iU>zog*#mfUB zwE>l9=eIPzA(-RQqZ%FST}D?8ryZj=1T?A=r<|7@kjtz;9uD=FX`^Ln1;Pa=orodt zOura1bLj`o+|S>PtuI_EfBEty1jJSR2BKMZP!Dh?WJ?tVS8UD41B z_^J5#PkzKn0t5=^Rg!(K1D@aS2mw>g1Z_PuVG%bu_Yw6Dau1KBk0sZJk=DqhzL0KL zoQ&eT^M9^BFUGkamA+syzjSawF!J+2VtK#m^VSKipMQ&q?vq;QpRxOL6@;PuT7k<=7}yhuOq?lwZ)w%6G{Vxfi)7 z5?d$ld{%praHUb{7(+(KuMtt{B5+^O#PIVpB1^QUCcRqMU|9A$=O!83IArFv7AzUGP;gCme z_Fmgvo0+UeHcV7#)-U8xQ&1m>4Rj6H*j}@MW)6~K9ryEOgsDMFOB{O!_e3~J8)0&% z?O?v!P8rk&s)|2ZG8U2t)*1V|6fwJFE+3wS&eZHB2_ z9)sNJun(xNM>X^l(e2w05s*m_BD26sFAt!733k%Dw{TL!f=8}Z()!y0LSG_ir=2-l zzz=5dRa>>pqhR=$nMI1#&4q(ACKf}K!ZjiFg1t6=(Eia>DEU+TN|ds8Kuuhl&Pa>2 zv%-O8JL?$8`9DsXlSh5{88O~4+?cgrd`1UP4UE0byKgd2HSXBi4OODF( zWe-mhVVZwTWZHw$PO$V@6joo(ITymVmx*7-QF&rodb#-l;3@|WfgXy@k&i0Y{svfX zvf}gu2F~pt;iSmnv0#QSvA~UgXHC|h4X%bS%*e1GC6o<}-mG^+wxS(7DP#d+B~gYf zp9Fo#!@vLSZ#}T7gTgF}zASvc_lzDhkI9p%5^hAcnKusI`qrV=ecX@%<)M|tU1T`# zUIN$qg?`!JucSOp?hVmo;B|AGUgVrjBBEBn*jbB==dC-i61U=&5><%X;!i zi;1jk#vYJ7U^BJyD+by!fZXxkNvwr6;6P7((4aEv2iG|MQzb=s@3-i{g7ctK);%1! z{`<-2APTCXiYbIV(+Ctcb}lTy29#xAeU@gzCVcuROlp%J^lkZqj^u}4ul1B+m{xfq zm9;CqpGK&iY#T+@!KcPA>nuYNez9PG$|)P|3UW7THc5DM1NYWu!I!dZ!Xvl(sNa^8 zCfic`585d%B-lC+Ec6A7&|zx(?JAmArklqQZs^k5FNO@_qH595rm`LZ7wiXhOumeE<~W3 z^~RP|ChMgvFFY)zpC_$U)qnHIWrMKJ1%vCkt{~Iu7Q(}xv3G_A?hJH@Efl5;`_awY zTv@oK5xIu8z^OE2zIot4a+F-7vvjw-Km_6jGK?NidI;jO*M;tmG!u54=qEQ=fmexqehaK~GZ6qHExs0iA|VLrcBQmAyB;n-flGOKmx=2P<8Vg~xFn z=n)OB&+~fMY4C8wJF@1dGDalqi5A1<$dT!W82~_wqi|z~V7q({G1`^=k5FcGwGO2jp6TA|5g&c_bP zG!sXT3G7k3S&!w#sOG7XhpB7>>g%j#R?Y)4UB4&>fjo?o1UO?V3 zTVb^0SA>s{0sulEj^~CK)h(l(L-sut^t&?=l%5B@5xMvd+rp8U-k{!*Ae5TtBlY-u zDFjNTG)e|}2LzCa@tg@~Qx2xgrDE0}g{L68RJ@E43rF1dg5bO*TroUDl!>lA zPm^ZxJug>p*##~O$J&Nm|9x-DvQ&Mz5q3-DVjOKW^QaY<)|dM?35^4u{2f#g3=XLj z5d~r!3+rMa4OD9sevYIt*079%&EcQa&ju^y$^(z0SwOq0P+ zUZNZpl`D+GSHIr6aEm0hArO}>s3VdZXol~kks9#GO)M^8?Q4upiLK!#_HYlmu z^|%GA^gIX=)stv$8$7h~Zl6y698N3&0E=(s68gh;CQM|1R&sw#Xh;V2k31MBB~&n* z+QX!I`hrY2e7qWkDuQ0kGkUwyBF&IBGZ{Krnc0{OG}WfwoRQ(r`q_fAU3C=spxU#2#ch|3{-?PI%&%nC zoxa{qrd<_(kZP|HRSav2Y9uXj;a(BVuQI`O(~%Bqk|_XRE@B;qO&1^(k_fAme?z^5 z+;VBfSd(fKOZPE{0=T>ZzZ?UP8;Q7*_~d$^6Ugf6MgAC;&8mskC#Y~4$C(Wz_n2x7 zNAj~ySLW79wT`-2&=^BkcCn~#32*$w7@MDxdn*O62%efe+h0=DEK&M(i7$}Qz2N&v z9YURZo?$tq(lt(L7%JP}Wa^oPtr|N*tMF|)+J0eaahy2WoQBR6D4(7fMkr!*6L<9W z{xAtbNevXxnBuPS03$eB1DQ<5p=O+y2ZClS#0g3`8snyKTQ35yaUZaWj;B$s?_$g6 z6fXPQi|*`0kpvV6Q7xW~P>cu?p@N}{A|NBGp&i%nEUT2C@bb47v1=I?wuRJnRoj(u zG)k47hgo~YXcq)}n^C{a3^No;IWgFVGVjByp5X4I@#BtD_71jPIFa%L0iW zqXvhIW7Z~;&=-j5`=zPkO;IsX4Ux$5{Q*}&j!tmRVFtrgi9(>QpaQGbLzDG z-H&g{$EgrdIS{;5jl20|)5);F?)}=b2A6zr8k)r?iDr3u2_59rv{Q4o&OUnRlT?0r zjqD9%I_hT&q<5MsF-cVj#XmxyJ|C5mN)TOj(Jhk7(f)(M7n&lH&l)=+catY~oZzbZ z`LH*p0SM?9&o8|{nFodROJ{R|^YOKV$k*6+48t=iWis-96W01%-kaF7lhT%2Iu(!U zt;TwhX|w8z(&|yz1=S=*#5j?eszw8;h{;9Id5Vc7m#r!q7BmU4&7^hXXLnmC#j#W+ zmrV$v>q4O<0q?>+QBC0!$ss;)jLd6v3Z=S{>q2?BII>4g55odnB*J)Vf${yI`t&az zxfWoY9Np#nd-PJzF%s-JKBd=cNAeSQoimFZUP5=?-4sE?*tHbUVOK=xsgr1FDVp#= zUV_?@EtAM9(;RZ1An-&3s*N- z8_D)D9bdTlXSRW@jhFEwY_E8c;`r*Ls){$-4yxW!ufd0)(+@E*P^Nm zauk4SiPuiR^^uSn%HDIL+`g(dbndB@Pu(eU>{F?G<@*}cnKi-qWU8Okth5dfOl1&n z7CN9j7S^qYZW$BfghoHn>c&(>m>TtTIuj-QxJth-t6u@tTd;Jt^t+k@`*r5LY93rM zi<)yCp0;6_=mtudo5-ZJ{SPtBfbBbHpkn)QfM&4z$%zZWijewZaTBb;=hU~hb7 z;F~H6jX9r6Ufn{tI6j_kkK9k|LhW!*jj~9^ak}?L+17YJaUCagjFi6St^`)8>%-e= zOh;AfJF`^r!QCkm+w#jEatR-m0{vIUR54?I=KCBq8K02`yR8-~5)d#qzqrz~KkuWk zaVDbr!F(*N?z3QrnCbe!3bzf^@?ULSWmW|A^9OI7=hlksIV0Z0w4AR`T@r|ISaO4) zj;^YkzkdBr#);#LP5iO^V{5Mk#52HYx)(7Qrc{t!*-{oA3|EsKeX7Q3ulZWcm>o0V_Y5u8ShAjoq&D@b_sB>IPzUG(Z#)(TR2I*9)rMus zEEfC&d$o|fWLczDoxTdqa_i#SqC#W*KmUmC#Dl61aIaTI*Vdea#ELSIl-gFr5tZVr za7-6jCncUt$zrtPs`{#S$|CmF?A2^G(cB=OKc53NP77pFN01_DOL&n$22?$lTaER< zVcm<`6L3i5FWPZY`)@Q}bgQc96eEZ`6szC_8Zid4vb4+M`L^yySH|wstuN8PV0=MZ z3TLsbv+`5WliX-FB4ZzdsL}H7h3<3n)y&0pVC?Hh3g zf5=g%f=xY>9XM?Ky zcQlPc6q7j&|b_nhhPKQ<#157OH?|?J)59Z>+r0#*GD% zB{(;9_rcBYdFr<~a<~{SvK5`5&CdlGWLak9F+m^l@M!Ua|P!=+=*o@Z@2UD(Wk6u2i$0z{SY$}>U^ z0sOx6`D&(bAH^M3x6_5^e8vJgKfsxtW~&UByWCfM?~W`5e+pg_$LOvm7NhJu-vj#J zOj@U3*B4=NGqpr>K1TiMW7IKn{JFmb6#Df`b<>a2k*GGRE4BaSf@>%qxEJ7_ zr@lnR2FMXJ-K2X22F@hL4fOF4=kQFF4)&jB?|y2Q#1bab4S6KRQYh6pVsVwEY(7zf z2a+LbF2$u^$JO3Ys{%*+{pLu=P9x#zF_CYDLhU74=#LzB+o)~PA=y(RqdgrAcDU{1 z&2t*Ky%CnE{8LN?)z8>HdyVodoZIDBx2?~QO$EgYvN1C7aTnwFgvOCBKgGnRh~yz{esUS?H>!?%IPp?t|z@bPt`{V7=1eONUEeFo@CXmq^m|YUMGToUvXr z{RnB)4^Mh|+E24y(cq9z*;!dXx5JuyvtA`wt$ht9F+A6F?wdzFqz0sAU4Ftw-jhoZx^Xzp~@Doo-ar36VTZpgbrI&glcvHjfBU~9Zo zI500k2cf+v;E2;B$1rsLv~;IW5?%CoqhSY?V;$J}h1^e@;d~$=!e7>_WHDt$_i>G; zN`}NR-J;YE0y5B)HxhCJg)I5mJAdAa$3{dV)pDUeq7rszAIVd-CK!8QwZ=qIy|iH33!rfv>1b2oaif(#Bjn? z6?sudFiNqUTrrTm!to~P$rLY$Zqx>xr~f0S7+FK;dNmeFzm1$;t_K)MSe`-5Lg${N zMgxxKhn^?u-ma^L~E6pR(kvB);ABPd9h|siDeyT7r?}MsrFF zDO*8rdBaR3D{=7gDBlL^@(#CDo6;}dK^mW68QrDz$O8mT=*KNR&rQCRpFF+c$J_sm z{50uV&v@PFD$XID_JJm;EfHkrHNA;@|LIvMp$%7w63%isg}LKqW$v#7l6cRSR3m~R z59NC@YC6(T7E+U={{Aj>vZFNSiE3%Qe#+zA+- zzO9fnnad-3qWt1@GgADkVt%~49&WK!swkEzCSU+?CBZO__|@SM`T|%lp^cT ztx6#Vk|{)e8@(dM9L3*wiqA()mmDfj9Bxtm?$|EyXNrryvA!ma7DBDCp=} z^hjXa*(W$S=&4^BZu%!lCZ^bE(h;TQE{aCnQX;b`0Yi8XRP7z-%_rIM^j5SDjXeS2 z%;}$=VFpW!Ta3Vn6HD1VK!Z%1oogYew}a^Dd?zla1|vTP=i0Rw&UHvR>PSw|;-gg5 z>9z37rU?78r6cMH{h%u4dMR^hPwEiMzM8 zd!+0zZlbTUK`1zc#5XmvctBDI6zP%F;1edPgzb&jskx@<#PIs#3s$K4ihkV<&TSsY zjeaveOX`3pT5A^6!)9grJp@bFFjO~_58=4ju9R+QWrR{)%_C4c&`1u=F!glW?7dk6 z9p9CHDtI@L%3#ZDT`_`pJTHgD_eo+%UOr(%!1BwdzVVfJYZon!)r_&zR2 zMT5<80F?{F%BUP8jq`p@mmgZ0<=>V^XxnABwH!__KShX%5yRS>3SLPOkm*2Su*K(< z8pO?~VyNUY(V5p)4bV59bD6jfJW2QwCA=NK{zixzcW6qT^xTGM4^(Nm39s%aNNock zoYBC;1;Vx`Fi$koQsUx3U`}(yMc80-1=H8--ga?7CpQ=VFt3l{#gKU}ACv?!$cJ}i zDI>IDb08ALx!IhHfbQ-k;7oJ&3E(^Nw0*-lPFl0!^KhO_B}Ok|yd+{`X1RD8m*YNF z`}JCu;+<49GWeJhZW3ODlno>>6DtX_#aD@2oH#s@X@?KEMcLfy7s23~tbA=ux+|of z9Akq53?9$mfL;p2RLrQq@f@}Ei}te`;jrl?tcdL@Y9x2bQgl?~XZvjOK2(AD>hkVR z*)fN%;JkOsa|?)vMfREy$>A&(nV4t}C^vllMeBInIyHEQ)q&~BQFXW(mfNgC>CNw| z6T_c_t-N%0?u<(Ym~dbmsG4AFoI%3;Y>(i+T0M6$x|=r#D}r{w)W8Ey+=O@R+*}Vu z(HrYwqQOtD^CWv$+4 zRgfZ1ZW-q?*CU$AvE$?L#4Vt^W$Vt)M`>gQL|iP!S0uQ~A{sijX6F{eVf<$Echifzl zN4L)>^kR~$26dPUxikw27;J;vxcLV%r_~oTvY*^m%;(LyLd1)-R zKc;Gl13ZBAGK29Byw&kW4QJ za>sRNp4ChjXsM2>zU_(r8Wqm46G$Hk=9m2wwK80z9%OGBHL4GJ$uXYx{!Lq@Wh_qz z>RaWrR})16x2#-ec)V^@z=IFD1g@X$$rWb=RdmZznk8OQq@0&yO0gp^bu(vT0KF2}^4&G}=VH z8{2INP>uofsV=#x=;&5s8x9iWy zhszcq(<+<#FgP%px=^Z0s2g5ZV7D>5a_JY!rM!Ej6pq{C=vrgPvN}6~p7=vYFi!q6 zZBw*PvF#F%LPfg-QcfCz=rq#?j()Bn&6MHidJU^=t}?QV0o?rI)cxNTU(sr&1-z*= zPh*(>V=c#&^0W=HHk*iBS0H^~gt=1;t)ufZN7_h2zq0F3H}0S62%d(sEK-`i1MdZu zkf7x<=*sUZaQ8LsQnIShmS`?>tFZgW#}d$G*GPD*=oW66>zVASLe+XTHakyH`-flk zla!%{I%*c$d=@i$`Tz%`&NWBnU!XOS4m3Gl-CsS+Z0DV za&@y_-hHWN%f7_7AFk57d*Lm2LmhSD-Z=hr@#x;(tcX>xpAL zXb?ZgLxWp#NPM~NXFO;lrxjJ0Mj^{h;);a*vK=x{x}y~IY6<#%b52DScSvy9}e$+xU6bWOrZ8w4b=YW>GA$l=4KUb z2p57VzCUu`J35g0MSIXyk@Bkv9Re<^bAyCFP@`83{D*ZHJw=tv_S*K0=c?lp^!i%B z9t5$m(uM7-z!>0b`rGyV6&{OC3y*AITEAFhgXRLcGpBECd92fSF zs^D9jLn`>D7TM9zP6Seiw&`yp5YcJv=t6omswdorTP7v~IGvA)pUdo)vGt{tHOfPypkL4^ z2M}Tifg=Z-ajLC#QCSVx`5-sK%9yA4_T81vn?l&@;%cZC8A<}cbHJTasT`i4H9eh{ z@;xcu$T8ItpIU+>vqpQ^!Y;Ahv-!mXcUV=9BHA*w>wS+t?|p6@u*D{Hc%j(e=!XK? zd7iYUNQb-c|7PIX-^;R&vOA*Fdn^%rjKGs8^jK2yi}kg@?Fif{O<;zAo|r|=!JC3i zSAw^?>-xgm`WY(l@wlRKr33uf3QIsmur5twpCjQm+20U7;^T9dn-g<21R&#<>&W{{r&whgI7cq z9Ya@1uR!;z-}HR$vF@jV%X4KXR)L%91_f^dQ2dwX=#2+ih(E8{oR{yUkSx1&QB zuC;4N_jo!Y5R}rD;~fHi=bvJsfS`(i+vtbu>E(x>jIK`hy}c23Esq+B2fkUcCOFW# zGx~Zib^`CWu>3|QRix|FSaR&v_6+MFs#imeGNM7n$bMT*O;xoFhuq^k-9iLJgbo-9 zRFYH9Y?DR;IOYzoHp8lFnI^u+do5zSoUHGyR<#7V&CaFTk#2dL=f(n2i;+-3=Jbrn z{*`Jk+MuvfbknbTtSsKW_;i7H=Oy;8cbnIPd`o`B71)o zM;(vZok>zV{{K+~$-0!AAgr{xgqdx8;%+DNvk2dyp_)#vPvmQuP6Sj2NMG><5VN7F z*>yr+bw2_AL7?be?Yx`EGBi_wTI9x}CzWs({>cEHCCa-gVFA87Tcgw$5|Heq{=hU| z5A(%elOYvDn4~m>wmB;RRp>~!7<;u^4=wo2W<)a1K5vi>zT1{cX_CrnM%>0g>0pAL z@Nzv@M`+%)!a9hoHTqxOg_5Wir&G^-HmC-Lshm+!)~ljCnArf6qWY^y1+JE3;}koX zo`C=HrBH}(U{ZxUH()8dVg$K)_+U3QMP|x5V>3YIP!Dw!AeaQSSIuF*f%z4wq#2Qs z1RJmvU|Dy)#m+e$A(**lmInEbZ+PQh*ug$y5e65%39rz`N#b zlnD*g3A|d*Q<}PA@xcpQfyvSwFjV-Dq_BbUKU$B{qCC{2%Z2vP|LOsAZ{Z$ig=?cs zfFIZR^~*n)fI%%hlGrJ&hLZX0f;e9R+CByW8zm6#)_PqIJ8$>`$FZ z5a3j>tSE;tnXDQtHmzy?V+_Kl!*xq6;iz81HV31X%s%mLI->%4(J&} zu~i6|+dm7PZ{0XhD78sy(ny(}5(rDcJh+16XAwYLt;52n6n;#Iga$DAE%e0;m|#kD}LO!lmd2aLQGWBTl+%(joj4=d*m!_|1_nby}|ArJ?SMysPAprR3_P!2_{~Pnv z;1ZU{|MPz9@qeiaOo>2ka5}(5rs;NR3G`*nz1wFAZ#h65q#R~O=3(O3?PqGTcHsp4 z6&P3XSAl8AN)Dulb?!P4Sd_^n9@sEOKn_i%euP&54JH%3 zj_W(yDDt12M|bQ9snZmfy$e8?p8s>fc0?Sl3#v+Iu;TF%XFvLhADYQEA&; zoK>qvK)=rt@X``P7P;`(!@bmk=I3vthp@-~(FK>-4xR7VO0Hyp+SZ)y+c5KxN9dc1 zy)5T&LHYx9=Xd13Ce1Y6kFPEK#L}qxdX$>&9XD|nf6*LwxD$HR;DS6QAMNYDF}CNN zuZ{Ou%+ygRU=ToIl|rmN_5z@#))??uAOyNv{mo_IT|+A9b_=6~-AaTe5CIR+FT+fw-0frUFl=P9mS{KXOP5LtDc#+xh zePf;#076qmOAVzIhmfm3S34p@^g#L|`~!BdfT`<_f&tMKw#$9-i|{V726I@wRYxJS zBM7t2jmMFj#}^Rc`-|nkQcn~N`=aPw!s^!*veiWYv$6_%B4$;07*{-AS z7G&L6;Q@yLC@HjY@Y3=zI*JsHx?3G)B)Wgg|HoxOO&)wF3)GN^L<%6BV5#if|1;?3 zGsXXgcd4NF-;7=W%;R;UkXvw|xUx#F)ql_mL-v(QXukX8-yIIQfp&VP3)<$iE)Tza8&lxau&7opH7dOd!xA|oH6ksmhQg;G83q_ zlYJ`xB&r@kg4c^M(|-%E>H_b&7%xy=La9wfWCslO3(4&GVKsI;VAk6Xpf*D#;=JR$ zF!JI23jA-GX9mg-fVN5UgFj(-YoFR=H~7jH77`UwhfAgdm+$IJVJ?ByS=f1a|1G(I z&AHEx&i_(xfWr>1?c@E;4LkmWh1w+;C^cla0LRXM2e1wEwbc=d>Zxz(XEqzuA-QQV zaixv@!nTjwgkalRfo*TT_hCR`*qGuW0068odF%BOz>K!+yUg4F(7l6@VD)_FBmD%b+j`)&EY^QF z5s!BJnGef;cr9q^K_aYGFf@}t0K%+?ar{Lg%|C@2=1%*jN45emSi%eW_t7LUuB5b1 zNeHtRnVmKaMX&%A0m>GR#{Z+30wf1UlRV#l>(hp%(4x>Bf&HN2*xz#bUp#4wB9S^D z1_x{X!m5BN#bFisG5i$WMos|)d1Pp99(rhe>y%?_O+?Y4COiDc+X|BZC1|kvXo%Bf zR~-H?c)9%pUNDLZbC2`HBNHB+NZ@}5^Z6rlBeF9AG?=LqhY?=lH@fz!u#dr9flIvA z0vu7H@BEQN6+;u7wyHm#C;9@$JQ#Ozfwc)d`#{t2c-Yz3`|ryCQmxblc5l~4;Cl$s z4QuBeFx;vj*=m917cV#IWwU{07R^=(OdgQfC;-{*PWiw0ZSIUNVI|bit*+CdwB_T) zv^#^}95iFVoko8gGmJfE{ePq&sfvkpw|_e@oCbQd=Put}9`*aLrmCdbnjY4Q!NLnA z%f~*U0p_&7WcqJ6MNaXbgTd^v^0#cjdOr1B-T$n?nxbxiVB%}TnlyPLyS)GE6tiGo zr`U|P{qKQ@lpJJ9l!r$vPdvaENvZ%`EyBa|ZSQ7@>-gIIyI;&)^u{tUF4-omjSOl| zJe-v1syl3{_+u5|diPD^a~?$q%^~Die9@ql1cWPw97;7qaTxMjq^5s;^{yf(J~4saL)+k+ws*^ z+R)uY%T}-AEJ>P|DiS&9{4n_LB`B93b2X+t;20Q^TrNCxz)2os?7@NI??N7qR4%_i z4R(0(rr@e0Bsw`{cZ)f>;AS}pHr(<21D?oPHp>UF-V0!TF zvVyN&5UVx$H*%#kdx^UpoQ~VMkUIhzq|p2Oxyhv<4CN2Z-$hNk?rp?w-O19G*+J&+ z?ueUz=ylkt!P`YehY}~F2O8Ev;;%LtfHz7vH~kh*vujC zcl$=fL&? zh*99?@&jh@S_<{3nNBZ4^65(K!_lj%1B3Xhcx8i=AIM^N0T0+TrdK?`y7Y(l7!wl5 zIW3(*Ok7=%9cNf)>$~so7VFxD_nYO1+bAV~PAzxml@~Q;56@I1uZL4Ea$wA)p%L;2 zWX?S4(yW7jfA+4gc-AS|6#$i_gAkM)8na}xt2;OIy(4B=KL?C%)-fKgCb&fIKP+FV zK(3ay>VWndwhg80XfMm_6~TBOML`bCx9RLg*s^B?ZPH`G?|%SAYVteZ)@>A)Uyn}t z(BL>fOxK|(DH@>mG41UQ^ETLxy}OK6@IOKU24F=D+x`r@U6+4qSFX>~?b|D=fgQXS9f4ULE=m!uES?Fh zx_s}mc*2QmPu#0g%Vs}H%z0mNS3V9~Z}HY`68 z0GJjZ|CsPHkRz&*^FA8i*b}8}^?+%f|J1k9*57|lfI-F}&26Q|Dll9tVsLywy$%D+ z{r&J3o!E3nM7UTYsJMw6)`0f}BPcR|KP`Ez${$`3s7CtFe4Gt>ELSjb{;+9f!QC1S zy|b)|8{*P&LzloS2#_+rOnD0O8^g`J;QC1%{ig<^MpQZ6PsTooyLFuX6abaa4DxW- zQHhIJ(k9Qqc;vukNOt!GLwKJjX%B<^dNk2a% z$^+Y?r}x*f4SU8|-m)0Z$uE;IKi}WI#Tk`!$;;P0lEvTW!{s>jExr{;{y3Y(HujlG z&+2yO=sYL?J%2Ay)7ljA`rvr|=7QRsa{#ya8cQl~CE{*yG?c*s15tdZgkksi_{8fI zEE!;}6EG_*nGth;_*@IQ`oeNNVAetY;L`WvcYq((?NGG5WPJ1T?|zbyGsnOy!I$-R z7j{RlcRMbZ5t>&rgDtqu6K?sMg16ExFBP^r(qWlU9*i(aqxFxEzMqt&p`O*490J4( zsMc95RQiQONi?g{IjGiUwO>4m#Aa7U7_Wj+$-D_1c?{fc4;{S!?bdiQIAOf z%`X`Ue_gTNu~11W$vYS+fOk2q`=Y#zM7ImBJTm)E9>(}qKF=jBf_~sR;*bBGw{7>` zC`^9xCL@`b#Z}dRr&E6MTX;nG6z}qcJqw+ml?wYVU45Eb*Nw^pqTRA&;{I|7HEhUwEcviL zVKL{lj4gh@8*-&jdLO^smx6GpN?qW8xE&5lXxVhl1orQY|7@6rMP;f40ny0DjRgPr!9#wH$>C%fy1f~xl-A(!&e8qFF2^GPo`G4tjP9n#`WA4DMehLFVsSFfMG7#3}01Ej61mixdi411z;cj zKl|W=HRff@bFzQs!f0QT!;*=5-^bNHO_7K*3Z;l_LK2{D3eRth)g<6Al^%02;149t z4hCW$oj(OYW6g;)b^7nm5b!FJB_&GfrN19p*!ydpQTz{FYiVSyqo4iR)K9VnfOW8- z2u+cTJ~t{c3P)zAXjkxm9BRHbz8(+fRaDACrN%8o@J3Zu2#W%XCTi+=^YY5ejCU}K zjexKKmxaL}J6Lj1Ble5WnJ9m_??_C*)% zu()8z>H!o zT`c%h*$A@iVwfoy_XW3=CO81hC8ZFZW7KcRMt?j7S0(aJPjrzt_YQ0wRv8VNoLB!k z>2BUh!F3pxn?LxR@$u-s0hYPuc@%^RBq8bnjvD>CR3ir3I6}MaLCY|L!I160+Xpv4^XZRRu@+;GF07!@Dr6_4&41V?>}I$ z7Qh|u4sYHmV69YCx^W+HM_?~^CsD9X|)zv3oT%VHS|etikx zh51~~MTYFmOrlp`ioNsqrA@$3XUb2Q@f;;E<<{HEUiMz2om0pH?&$uK&WTd|7vI6R`ihqidO8M^Y?jmd>ntC6AQPFkLMU2vF#Xsv- zxd7O?NwiEeD4`l+6;7I72fmtm={9m)&8LMUL+BVV9LxX!OIkh|Bfi6s2OG#2VH;T!ZW6(~S>v}z zzf*d5@`I>H1!Nc-=A5XIZ2Pc2bJB6G7$&_m!0yxQNFPXnd^EGy=z zZvTp@=BR@s1LQ<>k^SwRJU$TR;=GC?sRwYH2zB$oD`ucdxVI~46f3c7b$e}jVFeZ- zU_OVFFN7U8q0ecWMwyQC;9Cps(8%3$oclj@eR(`p-S__`N&}LqjLDQ_DD#vGNs%E# znKhy05;9yG`~Vzdt_DC+u_1K6|gd z)_bk>UhBBDP!3V)IC7}tmjxWhf3Lb=N;qnS?ySB^5PYJa>N>7i(RzVn+ntc`G65(nPikF`{eP)!}y zV0~PT&(rSkO^cXK*jaU0@vFjsY`uLi;h|@xu=*F~R!+4|U7L<;taL(&vE8n-t-ia5 zp)zgf?^3;<)fPiPKXGmOOr5NAj?FxEIs2IoA0a6ljnf&wtxm>gd-G<(rieL zXyj)UlN5)G$QuG`%#r(cRX{Ag%2(A{cE>=_b5YHmW|J)EG<*uTxh3h7NPnXx44SoM1XXrIwRytsocT@rm)Fd7>{jWSXP@UO^k~)I z3C{Z%^_}@@%&ogk>hEjxQ{`yWk!)A$Ok3?|GpI}?8a}OW&7JHNb^(Kgn@iLOi=Jrc zOh^f}={9!4OG%HVbQK4>*W|zU;pZ%J4sO?G-N(hpcfjg9!$?k!la|Js(v@v%FOGqb zRcS1VGcmydOVgcYlqVwPY8Km^!zO32^IN)vDP#o-jA~ncvZN;LHo`PLWeE+6a2wKo z;sm6*-Zwr+%%Y9sr2W#Ehjwm%S8|Ad!C--Ff3?d$48M==tkzF`Bx2(BX)Hds_0Orl z_!~X#>3e;Io*$pvCqiW;x@&a1rm|7M&bfQ@1;IkYP-C3a?HCt~eNkfKY2BGN8Zh(L zB)gz}zEF6w%yH4XcD^Ex_Wt%C9u&^^oVpKWx|lWT%lO`t_d1It>^qV$w#GY&r=Bre ztn{fbh$r`amI*gYY1d8#lk}*FcbLI#Ik=_wTiPH4X6JH*(c8!SGOO4fv*toC{Tgfu zj{lN2T0~<}#@Qw5&Ws-_SYG+f;Ap?8Z>7y^IfduMAfd1Pn49^s@5qFxh4(KT{o*5L ziO&m}w9Yjr6z!`b8@Y#Kl{cJMmmcLW2ZNi8HjQ#_$ew(c;BZl z$zZG-L4*9k7L;0APW9=(RY(6tpdq99tB2A`>uk>{fkma8T?*wpPAEx{y{Kym*6xeY zP(n$sVWyqjeZV3y2IBobSw39~g(FOaRl0a*4ulp$k)m^AGX$*goO|P3AwYRR+t>*7l6rp@+oc2h_&5Of z_@AjvUFaZzcR?~{5lRGFuo$v`{Dtd2N*WqHg{@mq3I< z)@;w6W#H*WZ)R+gwR`%}tmc7@=~?C9>!(3j#;-A3Lye3?Jg7#of(&HOl44k+9qf3h zZD(g@=Hd55XPDC+0jnN>jYy@T+%H}-&|~Qk^HoDyJt7)v#qu!|ohm>~M2m5G^FBRp zA5TrSfB?n;jcKb7qQuQ3gUmopZb-|LjBni|3%OW?C4FuS>@7xMGWL3VSOLFC%!fOP zZcU>yecNUUg%;&9qoaEpx4;(MP)<$LkpXWq9{Depa8(TPs&;IpN0!F*i(u#(3Lf%j zZrFGP8+g&CtDW|PW7ujF3$Rq%xILO)Wl1&Ta*%LSauGnCz4QJ66~H^#fO%1Xmbq`w z=w@90y5V+jH83;dB%(?E=DkNP5Sk2KxB8%;sv-H8-%m~SFaF@YY8UOCdPgdNbjd_p zeaRD+k>mzh+Kri>)i>*uMHVEjMI_uP<7($4B+MT4Du9#0kbn;j+Up&HqL=Gfxb7Wj zNdc1!fCaf~2-y%{m~l)dazh#dUSwv5AumSYV^kP#)2A;?T1{_tOhE-wEyMXCM?ChH z^PeS37tx>?(gvTH6;Bp!0_-O65?OV!1LcjaMQj^+Je**jbtH9HPrwa{ zvr2$ecq(zb_i%%E@>GMgvZ+h7W{KW}d|sJLULwT9cih0pNS*AY^)!6Bw66)ksX+bP z`=^M@r=p(9X-yIJo@qA-SXJ4xqT`(iZlX|;-=0`_yA-RWpk;9;0W`-Gg?}9;P!d5e z(%B&xA}KxM;sgT3N3SwVGQ;lQh-MgWX$VWpYLKZ-Jid&EQu~8+qHH3oW%m>B_e*2+ zf7UN5OyI@Uoq=s!0`w}!VDVCywsMl4U6%V@kBNTb2?MEv++vD2$SU?C6*IEpd7@5U0n+gwr8qWiOQl8W?7O>5uG zT+ODZh-%;6WgNvK=#@=RJo~l{vzvqCSaI3@88=Hy1umbJ&T?Of8v zQu$t;r8#lpOhS5kh(}F+o9yPz1xYE^&HOd=N|H03*7zT*K`*levpT%yU9c}r-J3u1 zs(01YMh#B((u+}HsKOiYO5ufF@?mhAO2((IWpR3q{(7X2twrp!a&MJpZL09}X5Mp! z<0c!%NaZ<{k%%O|NKaRID$t9(M@~Y7!_6|M7qgQuuZ$@=&`mTv4(Bd&VwnvqoKOz0 z%M$8UXg~8^$-74B>__lr&VUg({Pek?8N?A_X(Q8!hz)7GhScX)pw#frlB}Yp(L)Pk zrxX{qr;B?JC~|p^i8d*CFXWVpudT^>F0mRmD^)-NpDL1?4|wDB(xc3x_KY~5o~_Vs$1o4wB!9#LZKE1ZN@YFw4wc%@g@<{jq!W#D%u zPC^yap-+m=mx+aoVM{_9vl)nogd9N80Ygi|U=TpTPg&a>W`^j<;xwa^KJSo4mqYvw zRj;)bszOU4omqAg90RSjb3o)B^9X2$!@Ay7Z+2*D5LZjf$n@rBo$v$$-m8mB$wx|7 zZ=WK^B&iFoE3~nYg)&lq>o%gwtB;c}b?6RDGUgr1W|*zuRCgR``Vok7J3+sb`^oW_ z4)5h=ToQOz+sMxFJEZczno7IY|ETo2qA*_Mx(72$Bk$%}Q(iYY->||f`dHx(-kFMM zGf5#)YNUIO(3@U`G?c0VZ7@Vj;(~&<=azDnZdcJRztt$%a5XgH*VJ&7^`V(2_92h4 zi_y_VN41S_y@>y^VWfgmOG-+JAsvI}y=IA5niFhdGmFMPtpqwb4wiFQMOBtlKZ!48 z_Vu9ZTr4)3`iaH9c{1-!GtFB6OfVy8!~NU-A@Jm1h$T+W=TWXxp4iaGzv{9$LjTgouuu|`YK@W;uNrSqWk@`w3rqOr z*jQF~1`)vEsAFn3^nOB`ZM{3umfY+-3tI5Ci*`PEKrpwOg&6t)2!-kLj(G6exGQ=O z_%EV;0PptbVX!>tj6A#xXASo@!??Lb3`M8Y4fDtEMO}D+cRoBI5{)}7>z66>l#cya zUBCPG0sJR)!piv+Ko8=bQS9yY1cbrBa@dPHdwiP$JF;4XIj#<2#aWtB8Q%pm5;~?W zpQcl{LpAvb@|v`?vti6VZ*KX%Hd2Gm8O9gyd|V%{6GYd2ez~1hbicABZf?vT z6CEABduNRo=0G2z3rN8uV?T0*F-z6g{)Zuagl~Pkx>_~GO=PU}y62s_#Xzr`@YU8? z+>zO|f`WNyz}mgs#nXYyD$IS_H6>G#Ll5CGnVw)_(M_1L$tnuZv5WT{-!VVx9kuY8 z|LIyGu~cq*_iU@wnbUFmD%yt>#w#lc59-9?UG9XiVpa&<<8QAQ=Jza(mQ!bDReJ9b z5z$JlS)LE?4P>Po>ntF&h7H>wTCtHuY0ty>MVa`bwdx!RH>+urg*0exm*3sUZ8;=n zu0Gh%aa(x6;%4n3%t|w3R$rHt%thTg#g%sx6?yBEw;02SJmhOWHQ_5CBE!7~B)E-J zN5F?!4(ge5Y9Qd27Y92Z|3ylU$fKr3MI6jR;%W`{eZ-=O+bm|19vk;2Bi=Y@yUf8S zilSaWRibMA8$!=;O?dbmshF!eWLjh;RNwHphL4>g&%AK*%W{hH>w4~xr_hPdrdB0I*{+4+6nKPK^j{1E@b!$_VK5NTam3b4H z$b>->+o4_1-mQ4GfS>!p6lIyvbtu^8s9dVz@*KT}EwAX$_b)`&QLRvEOS@S6^d~=F zA{P*jxTD*;hR12(8VnCku0@^Jq%Lt8rj!)B+DEBG4J@EJQ%@9@Q#jORcYs6eNa?22 z%FfM1f-*d`CtO?aRz>go^xTS&j-;y%@_)uFJ7w1D%dYUU>T$o9GXoZM>wbwLHHK!q z3w11kNQU_G>e8D-uT5o-$kwce>B7x`d689-y+vy}T2df5drtF~fHxd)apo}4gBRAv zjZMC{M7W}9_LSBDEfpm=T6YflPjFR=u2~3#QUV^H1*}vxHJerDH`;z0N9Oaqm1+^= zASG`62y#J3>6KlP3xW{gBhrZF1@W$E$wcUede|U*5k}y#NudaEQxrx^*XsG&`)0z|@RX z{ZkdtTwOUlU+fx4oGZru&qVqk;#lb7^D>Wiul9r@V#Z=m$0A~1k=8wxlkM6c8Ys;) z7@(Qfps4|b2Nons`8OyVzPV+Wxh)0ABvn|*<(nKpuyN8t_A?@;N|TDb=vFEuj!{ZH zxF)Iz;T(!m#p1yM7gZ8~B!Zke`L4)QWE#?&2@T+c6zAU!fSvUR(r!2?T=r#f6fK8h62cBYjsgvm9{LeyItbWJW{*>Tj1C|R; z6y@kV9mZEUHVDpJDk1cjaCymaSLh!G&fG+1eheZFBi6)+h3IWFTL$Zie4daOZM)kg zDE5CdEUGBaYeB3QL=C1td+UhT?}Ky`%Ef%hR_QjsZ^Z`%fH#x<;+W$G8N~raZnrq& z5MVLKuT}eR2o<#DuFUUy2jjLj5(!-Jxae(pKGztw*5oPG?V%^S$vW<*xTdYb8necQHx%?9?)TM)4; zz}67V5U|Ao@0QF=(nlLZ56q@L#KR5iJ_@&iHAn0xaD7xaHtsNV+8E}y$I-R3)XbS8 zB#Vj#0{spA5fKN-|MV5Ol@j27HB{R+yzm+&`JCW{Z4k7Fusm@0VieuoEaKJ+CO~l9+Vo(i+j)Co7$FZf7OsI!j)Hb zHVE7F>}a84Q|CBjND{jJ!wy4qG>QD<_h&+6P@>FM%8DSDiO5$dnN244P0yTFko*oT zNWcr?%tKxZVNY+`giv4HwgW{+uo{CKgh*@mx+AYOVZqY377Ql;I`*rQRY(A24Csv| zlx!5P^Q#VGJf(yu{PdcHcQlf9D3xbuYMNC3S~NEU)4dZbCNWpEq&jOzx5NT z7rR|Sz@=u!C}I4nEwaR`knm!W5u>6%hE{IYfu?hmx$Tqkg6>gjw!o)xs8F@a`Orxa-tSh1D-N zTv-$M)q=uNV0Nd?Z=WsG!{h`YG1{O5?mKcFs!X(vwIb(CvCT>#s2U?0MmrBLF+Xjk zu)Go6>e87&V+37WXt8xR|mD-Lk?l?C?z5RBP3n zr&L#09M2FbfYMY;twTbY@?m@JabayN>cRJFPkK~4aulDsx%;&|Pl4HFfQ<9Dc9+g$ zQQa%gjd*Uk(7h;*fIKJn>tGj{~Aecy4vfs?BzdAJ^gR92mFCi&I^ z8KDU5kOPrdAa+=l=-@T$-;8(h^HYQSFW2oldYn*l#RPE zCBj`^e-c&fq!~R4ax4gk>Hnk@NpaMWl_#7RoMuPMemGJzh+tsVH2tU)vxCY!N;os= zdFe8Q+02!N*cz)aRQzBW<)b&dWLv*l63P2i_*?(9Mb61z`kqJ zFYKW2q{eRB(^pbfibk-&8;sGE6P;ZJEv4Dyl$DOw@P0Y8rn_ZcTm9FmQ(k<(bpIP-F zs+Xx4^@-~E^_B^7m9tdHQjHiqCu%%P+!S2;ZA7hKRijKldTR{?BQUA-tlvw@>XWiVTMH8!0T zMK&hn$!_eKePAXew!;uaP-iEAs)(DR)bIj+>whr~*l7d_+MTArIYT96i1P{d6^e0Y zz!g8_JzxqU{|M{2Un0vJyijSxsD%P0vN?cQhSm}N2nR5&sHn`(z(2BsOaUbsgFOC| zvSflNq+ugB_@R6IQSDiB5}C2NH~mZW#&slPNg|Fu%yvvFE)DBMVRlX=0@5BYfD;V1 z@XZ{tc21|kBp8?EwJMyIkXID*#2Br)ATx{Npu!*Ysw-QQ`Mzct?Ca9Qh#i0S$0>oIj#TCHS?EQ+X0T3LdUh;Pu9E1?1 z`v$EXP-Xlu;xvDrB|IyKTfsJvJce3+ge=ionLeD|%CxmPR2oGn7RE9X#X!}_{znD$C&D)5?yktU5^ZP3-fSKV) z0QDI$r432zY0Xvc^k2;5aatn)stwp&TXNWm1Bs-dyfW64vX1Retp&1VlA72 z+@5=i>^2}yI93B;!!)A9f$L_bGa$JJRtGr1@lgY+_Ps)3@lc6$i8Q&g^`D4& zESLf2k<6V({@yL`SreQ<`B6RYM%i5XD6Ydt1dt4AWh)6Xka_S#C-2J*DUq!;eDyEu znP8;|N`%xeA`A-kyy-$Gp2am)FqebT=pGi0zL6q0Svt@Hta>I#cnXKumZXeCPdY?!o=lJ}_HQd6Qc>FAJ z#Ddt=L5S$ufV&9wD3ORY?(}(fyE^F_w@D zlqvBh;Msr-Otd>P^;;ueE2sO34VRrm!=&mqse(vb8Nz*;M=KY`=s4uuuA);0(3~V+ z;x{jp{n)6|LGCew4JPV{Fj(XKLr|rFT)DQy9F{*`113r^2-6UuQBZ__STgdm1N- zvHVrTfq88;Pi4N)X0l9IcJ^B12JbzqY~bQO;{g##f=q#r&w4wy{6^2faM^pxf_c*! zFL{UFXzGVDg?-(gnVT$DOEVT&&0H_OSrM>go~&*$uWnd#qZV#s*_7^!Eq}#$8e8^F z;r+9kRb{vFDjGU#FS+9P`Uw=4&fP(Vi0W>&p4p}qzxK|cwRCznOTni*rK9^Z?P6SO zr*Jj(7z`ajnF3=I>ncZA+O=NkI<~`s>mCXeP4^C@i_LhBQ}lcH!>2ixUy224@*4Bn z?)scEvG}3R2RPinlQ3V_e>x4BJJFxIH?RrwUAbZpeYerEUk{}*uX`Mre~o5AV}LKL z>kFj~it8T*Usd4CZu|+$Nihs1rgXg?Z&9yYUpl=!{t@vIn(wjVS3^^{N;ok1j>+pQ zAGAelf(=?AOnY{nz%OaxPVGGs-${Yl?L9&Z-}6ENe3Qz_{Haa2YetAXAGo?1wHctk zS$yYV)0w!I?2n9b7ITZ&zit!#>`p3Kbz_?VX&Cv(V{GKYs|RnYHAjMA;=cTxt% z7PdH@iC4U5WYOyH^ZSd4+k6IZs_P0EM2EMbMsIDgZvb)SK=8!0+?mJtk!Ur_g^Z6a zo=5(;Rpl{MZftqkeVAgM_TkKJg&UqBz_h#H>H*vOCt$GX|k0<{n+U?L~CzGD5&%Prn>()yK)nR7I4KAuZI&&R!d^5)kgr_OUsV={M9 zCHFrppYIIQsso53DpZNVS3g({UzV_-DDjS!n#vZgW3JxghAVF^>}!-)e%*A2A3VOX zp;e0hR81-4<2H;h#II>z$LJ_Hvccp-Xu~sNG&Of>>hFA5)Kj?L(k(7@jH;!;EPt?j z?M+IXCKZYdE6xvQ`OG(C%j;zCd?}j@rnf)${qSIUo2{m%>FXE?>zto$!c%Xy41%b* z$GNW`G&W{2y`gxGWzw9mqyC99(?OMt+WjtR zzT2yIYJM%1HK}2ZE^L)Bae>h3+V}j~>6=#&qWx>RU9*a}JeP@6FitXJyY-0)>IA{) zeiGuOXm^rt4VF+EMo~^x+pnI{T@LB?H=qEX(;W=-WjVE7v}nj(@{re$dADS38tZW2 zi$g)^`rAh>`VY$g15AuP1%EnqP?bMRmWlgxRCR=r7z*yB?}|o+ivNying(CMp|nGV zm?Z*+2;jgR#(q;9Ml>N(QwRQ$PD@8@(Su*lwnVPZ#`ym+V~&ZTl%!iBXJ?K>@jGQ} zaqmf4uwg+*x#irpiJWF}NYnzq0d*Qd(rtrpkvt-6N^N=bEfkEtM=@T=@R&XP@1*j~ zn-vW!*4SVP07C$jBck2ggg+%<4@}xDJ>U*{M6YQUoJpupLuFv#wMpB{?yArEdRBxC z%QxNvyuZj=kXWJE=kAZf?`)k(7aKd)9-RJwgx>!Q1AO+4-M;-yp$H1Sf4&#{9y9w- z;VJN+phS<2yrIC6a>|p_YfqECfnk6s}&W2+*zm{#dM0N!@AHNJ%86oWw zeZ6C%787aMor(1Qfl$c^q2tr3@`s;GmrrCZH&eD!ia%7x9YxUg0?K%X(I8RZ?Vosp|(W7Rx_cC@^oWZXaI@y z56bK_YI@P`bNF>#IQhs|YL1-NkQIjAdACd=)O3Kfo0iQ3Qk3Y(A zpUWH!ptFmZi8M(RJLP&A&gaBG*Zf zJp(M%u>0|^7BrvL;~u=8htKyk4$7W+lo931k#ca&Fs&ZPJT+gU(nYVt?fAK~Zlxsh6PTbCh*okAaD$5eAd7E)zv&5pU! z2%3_vWav~WmI0l}!(cl4iO=^7uj-WYRKGV6+<#Jl{2V^7?fUlOwaJ<3<}a;Vpz7@_ zVexK(Z?=*;Bb3Mzh4$G=AQktN*r%*eWE+-ALEJeTn_T>TU9QCi_ioe+a|&a4uX}<@ z=g=Nqg|OVMN4NMgV47d!Ix!_)KG=M>bU}wB1fIcMhNGY7mq#%(AI|%i9AeAJ9}4y# zxy{@A@}QVf59>hkaS;V|J?`}*wr44LAjPJ7hcCTUhSC?>mQ_2#9a~(~{kWN`zSAaH zR7{j)+8Hp1Rx+Pxb}pLG(p z=?Ar2lH+4RLWzClLie_^uX&1<+Y}GcF5__gN83D1IQE!Z-kY&8i6S zF*oM!c;ed)HuoD$(6IQkOXqp$T>_0QEJC7jhotgNZ)1~%?{m%Z!zB7xtHZB4o>zO-YIeX_HJ7I&`tlCYm39M-}7W*jm{s; zgJ>9%Zy>e0G4~SQVS9}KF&-fJN1e0fTp7C?t;hAI>7>CKT~}V!F-QZuH{UW!6|cVF#@oNGDX`oE7*=T#o1FB*%<82cO6X)y;ud-&f z^lY_D!73O`ddja3^48-rCMI6aW$@7lYDj;x0+XWuy+Mk|s#S%cIpy&xB(48a<~m&m zPrGKUOxBdPbohUJx`^>z6O8v?cZ1EFn)G5F>=v=``K(qGj$OLbVv$5J2yKmn=$_yg7WI39!G+#x~z4}sKd)K2`tG)uYcu2jFs_=~$O6`fE1oU^P zi2^RAYzg5!|6XwIxDoel2dRaTFT5!aWcHY^{*a4Yg2V_EU;}9;)=^*=Sj^~4yqvKE z^>W-y1Hi2T?Zt;Mf@8bH_82}D199jrg;8tU29MyH4UmI+2T&+PMl>HG=Gwu?#9-lm zO$c{bMurZnUom_I#b~jlzksZS5F0R{$0xd*Hx_*|utOzWNr!A3W=_pUIBxz#X8fWl z7-wDOnr3At61I~LT}LA(YxbNJCEp~W%kMeRv%|he@3%OXy`;S>5fQ0F=QOkQwZ!Xz zy+m%Mp!w?1>PhXk+rWyn%d_@)cGb^ubKiC;XmkW{119M_iTfsSoEZ6pQ29NHofahf z{%4Iwt}ApBE7ibKs_4g2cx@PNv;xPhfZ<)3YjUbWEyFFhRuR<>%2oSyDV{g<( zHI1o{VKEVVnDfcpN{z*mXELb!blrD*bT+2pMz}t|86Noj@qSc666F95Gmd681?-=Q zCv8a_|Fet8_d#PA9|&d$H?8sQ`l;8KXYFfltEblV&SOSPhFhrE28>^bQ$7#} zZu<@EDIj4)Soam6p}U;ReBPsu9=gTPXhH8l1~%qAP;Mb?p4Rd0w%FRrrW2G^b1zS% z#9PK$qa6JLt7ObJsjb#ypWLEgyAZodHc(Tad@N8PnN#v%{-;ul6NF-k zWjuud&r!u1Cp)?iQf(+k=NM_Hx)ss3pArjks(uvDf4SfrsQH|o+$XJrk6jn^xh^CG zK&XsA1E|8B`*5iQo`KwlNWLGbf}qHnQUlJbBW8at2`p}CXMnbR=U+o*3M}Y{Pp3+b zbADg^8fDm6KEQ>!f_8@7$4*@1RWH;PpZ+>_gI*+}IM^kCjb>uyr0grD$J=Pyog~SP z!y5V6y%PmUHec*O!mg3*`fyST{J>s54$_o;L;2Vt#Pa4EwiJ=4-38CB3I-APuoWM* zji&q?QnHW|*d-llCwux-4RjY$G^ z(TI)QGwsh|x*DVE*ChA?Y42Th|0*M1Kna1dhy+stFUdi^$`7PHp12uxe7pGKl7$C< z;QKc^#%R4_K!w>!f0^j7W)~B8k*xyCVK_#i&DuIR$6Vh9J_Zzr)YV5c)o=n9DzJFq zd*6ZdVrX0dvHToqy~^GkDRJhJ9B>3s6y`jDnTg@S^ftmHjf;`+aMa50u0p4=HWc2a z35mqgpG`C`$O~Nl$_$#ER+Bn_fEvZ27LWHOYfnr@qnvbZBrWX1hdH|gd+i#?@=mQ8 zX{Vz&bAJ3F*R0XH#{b%VSWe9-G{7cSVbC*M(}JkqQ0$y}r$A`XMuWc+{=0pN>g66K zh;Nw}Yxya18x8s+>UEg?;FABJ({knrxmWECddBwTkI+6;2Kt2cYQ zQ2o0C7l?1=j8Pj;C8+y={)0lCCSV6KSCA5pl||<2k7)WFW4v#Q zDjXXDO}Doam5-(vY} z8)p9+s;a}`yCA~>(T_{{@c@;yjk*8z?W-c~FjG)}dS<%px#M}|Lby@I!Zf|=HmBA~ zYB=+)_>JzGJ0D!y^rf;Fyf!qk=*lnU^)w}2|1&7z)GLaiPKR-u(0A?3ACAxeu`);} zAwZF^@&)Lor4*JIE?hXuF;9_tw-6^NAn+pBx*-!Y#F%9OKX`amw)2~Fhb+-eE!|%! zg@(3I`tq|U3Z^*Q23hAg>DPt`r9On|`3KrhxN35{y*cB{eWWmZ-F9G5l-F%i9sAaA zQR@-Myz*iPuAy?}`+UmZG3j%w>z+e|jP=3OGL%PMmdzoA8Inmi(Q;Mr%tgy-76^6v z{$5B~|IA&!W@J(jdPd!3|8bnfbbmU_%I~-H-W_>L>)|OD3*!y7Fk)`+mT5C^V6&mZ_M2uhKUeC|EJq*L`<= zozj;*exdCiW|A4BK+zMFRo*(WJj*q;mV)U(4FvWKYq+}T>wxvb>iV>Ec5#S9trA7vHi%52KlNqwue zKASM{9(~uuKKl5tLs$ts(`$TMC%N2)Dikc1zKMB{^sg~EM>5O`uAN`+D(sSj_CcfY ziCE^DWDOMt$h*2V)jfHaRJ$h66Og=JEmm@e!aSGpj_BxfLwOX&?lt#TE>x z4+CZ0JOPrF@I6rRmo_GzG|&UX$7<7Ynpz3;iBM+UaO#1G4j4f#_4U@NF*J-F|pyy~Hl3$=#=y8IH42!msxGRmDLa{h(B+wHIf= zPSE=BMpk%d(?t977w!i~ufeN>{!6xUP`@T(RmEYH3SX;bKo7D)I5Z_bs|sNW#^Y@8 z5i{HU(uxI}|9-r<_uL8x1V;Xh*u4=`w-jxHw$C)^J+ohx&AbMy2|D-*j$Q;3%rrcJ z58!qlT<#d)&u>hh2ouP_mONxP(fTkRKgkOpdZSf=&gKW24m<6EPegP6dApkW5E%$J zlHZQwIea|8wl3#o)~%h#A>1ByyZ3X1e@pGU%wOvaDw^Q)ySTO=DU6N%_#)|Azc?;6 z|N8GYhMdLQzkc2gb5!}WMv%>VOd|!2t4}8VJH2t641UtOiq^rmu7gx&Gn9&Kf9zVB z?p{$|7+X!J=!xwAC~|$FZ$jCBi_eVo!RT)1+Nl?h~H@k572#AQXqCG-zT%1a)SKP)98E-5cAu3R@c^PdZ@ zV=viWasU5caOUAj^x<%B^aKYd`IEMKmz|DdUF;9J_?Qa9Kn`RD7?R&e;>he{8$Nx3vleVl>pwsAe9w^!q=2r_zQX1004R|O*I1m;5zV% zrK!lSYJGHt3;;lV$3<0D-^Jb@0KooCNpaGDW6s{SO@ggPD6yzS>8rj53x*~$S==2d zDOCFNUSHiky7szgY_ya*(_^JpW06V|wNDjy@7gKLbn|~G>={_uBy(o?j z`PJ0&cT?wyBfU!nR3q0yD^txLh{`*rGfLc-hh2zfT0Lzw2-CsYr{^^9{`sEaC3f<5 zVf~v`;7u9;P~?P?;WF$pXuHx#_Qn2WVYhoe5^WS{lE#~qUNc1vO6}Wcy94hKQoGRK zY1Y1=%!5BIkbrxldv7SaV-Qf%CBs8Wh!pvZ)9%E&TOyyv<>i2r`Ap+tAYifzz#G&Q zG?xOvYyp%`cKwkaMIV!K0O;1v1M>74*w|j=4RWPx7+x~AFz^uh*A}X@YdSWxz);#9 zKoc)jDH)E5j3o@<5QeItXA9-8pxmG|{%}hwRA`7s53lD9LS6gb>fy~!hZk0&5^53U zy}sQFBUd(yrW@wv)TB0u!te@2)0?X@yy9A8G&6N6mu9yKoi!K{$_7wb?Wnvu}g7^Cr&6^ehB!Zgpq>MT7}+&%Vhr= z{yx<9hk8GV?nYZn1>`A*6M{*Ip^DisPA=M|4b@3!*LlCnOw8m)3H~!N(59_ob@IBq z-~ssoiyzgau=7xX`DAYmSvU%Ek2*C@4N#W;FdO~OvhNA%MRY%CTpt8u28om8!b zGN2N{SHHu=Vbb@b_y0i%xxMn}XBx%CD0c_l0I&qyeid?=jF zeU+_WeT{Rld3QYfNHlgzM9S+khUPPE8$W4aaa=!$R6jTfumjp6SLKr6)}T_Z==auCbeO z^aqTdbOft~il>)vW0>B&dhv?V&D^c(`)`zR{XN_1h`-5yZPy3GvsKAd-l?(^VZ^nz zr)>dq59jz5QbR&7!9?k7F+1o8xKau&*$Xa~as~ec2R_-y#OWde`VnF>HQ)o?jAtm1 zDO3ae?9)+4cxZTy8p_UjzV42z0nWgy(hlvnXA`TZtsDy(=l2aM)|>@9Bk&^=`M!Ch zsR|tzbl>8*IqlOY(#O+zD$lA09CRFn9S^;Repk4)IV9K*496G7pT*19_167DcS-w5 zk9$N)M|o}i$o?uA(${WMP}cQ@xu(PGRD6*4qty2FLFnf6G0m|hsmot>(^-T`l!+~F zqk78xly$9OGV_LC*_3hJx{j^dTp2$;doMe}?)h$6CVfcxQGk#C*j5wwkt6d${{ss474I>G5*Xty;kx#}l>gn0nz&q~GVYnu*W%f5z&C8FCg_U$C6? z;YZ60f(9xEs|m)6mZN%Z5*wjVX-qWK9vWe#Y=vJz4gAGe`hB@9K3xAq>Rl<3(6qi= zr&*g8;0?Fr|$`3z~%b$0x;(YlfBJL=6J`%kUCU%`6$@$^=!F`hljJ+geNs@eZsVf$ipZoCc8 zcm{$2sX3Vy;;^!R@6HXs>q<{t$#mydPrRjYGsNNoabbSGZPEONS~r5dbEE4d^Ar{$ zIlvS07ygrnmi1VoTcTGuMd)3AT>b#^k={&j-SSdyNMf#!I{?;bO)ToZP>vIZoyj@ z^Cym&`#wIqiubZ#HK0YtQm(#DUit_l^~=l4YuB!Mc!*uq|98RR-2bhKj+UvfM_gVy z4G)_UiH37?RuK^jxw%?9J9d74a$H zngwH^X9!U9HRgYizBL{F;#>X9=qs0{Qi01qSJP{DI&d{WGzGWj)rQ(z)BM?0|2Jy_ zGcdUyM;_AKXr@x$<_Y?+sq!^qVW4B(fc!H1xx^aW^XdGLn zWl}-nZrxX4iu4w}?)`AB=d?R#`gT|W{I3L7# z^10_N7K*dwXake@wEb)FGbRG?95}4~OEp##b1_6(80f+cx#7MuY-6E|KF$1$4Mv8f zGXof(!WKGPVa|n5Y6pcHSUCIOmJp1~04BgOx1LmPW82ly?Q8yG2@^*U=MwK8NDP&o zb3-sPh+(s4hB?rgw1t^--+xDB4kQDS@7AtFOsMkQ2cT43>q^y zO7eOE=V0eC2aW^d@3mEkuo7O+r&kQ}WgT@iuS@Ee3{3x;-zz@#=FTf?y;a+JsZ|_swG67ny5IDF`OEIap2eOu`cv5GS@@3kMtw6-7iN72mB4_MznLaB8mv zViLKCVFj%_Qe!rR8r3!8M$&N0TWGIF(8E6B(f#=)NO%lH_jkeb{>)BQHx^LsI0H<6 zXO^*H1i6%zzsm{JOOTF1GOllpz()NvE%!%y!eGX;pPmKe_wf0$Cxun1XAIjWfBiv* z$*AUJrzywtkZ4=vks<91j~H`^q(XC;6KzWoOgZ5vsrofMe6!6~gy4cE2cgEN9F|OB zW-H6POz3d>x%)zGPw09$;KH={VJJ91Da|?qG0zFj5!G1 z0w-jN_4d24zmr4$qGS2M&x5v%X5M z)+#%l+nK8$AApe z`9%S9_{K+>2bQr+VSzi!6&4>v?PSOCsL^15X_kt=#L+-i(&ckuYxbas71P<*;7)kP zQrJ;91-9^Q4v3Y&6fKVwKcwQ|k$iG?3H&wRvw10JvWfHT>5h+DvH|^LPy3%F0%Eay zgLVsjN^G3Mnw(#6xy8pgZ^KBs!ar9A(7zMBpFuTU3>=@Sc=>udxOh6L zcwqq#r6AacQlb!1h(HnLyDKI4|0ofTeva+{Wd|=iNAU;uAiV1Oh8q6`G!A#DiC02^ Nmb#u=xr)t;e*sW&xfTEb literal 0 HcmV?d00001 diff --git a/docs/source/methods/criticality.rst b/docs/source/methods/criticality.rst index 5e141cd56..7959a5884 100644 --- a/docs/source/methods/criticality.rst +++ b/docs/source/methods/criticality.rst @@ -17,6 +17,8 @@ increasingly common with the advent of high-performance computing. This section will explore the theory behind and implementation of criticality calculations in a Monte Carlo code. +.. _method-successive-generations: + -------------------------------- Method of Successive Generations -------------------------------- @@ -43,7 +45,11 @@ distribution over some region of the geometry or simply a point source. Fortunately, regardless of the choice of initial source distribution, the method is guaranteed to converge to the true source distribution. Until the source distribution converges, tallies should not be scored to since they will -otherwsie include contributions from an unconverged source distribution. +otherwise include contributions from an unconverged source distribution. + +The method by which the fission source iterations are parallelized can have a +large impact on the achiable parallel scaling. This topic is discussed at length +in :ref:`fission-bank-algorithms`. ------------------------- Source Convergence Issues diff --git a/docs/source/methods/parallelization.rst b/docs/source/methods/parallelization.rst index cd6773630..7f0b9bd6e 100644 --- a/docs/source/methods/parallelization.rst +++ b/docs/source/methods/parallelization.rst @@ -4,3 +4,648 @@ Parallelization =============== +Due to the computationally-intensive nature of Monte Carlo methods, there has +been an ever-present interest in parallelizing such simulations. Even in the +`first paper`_ on the Monte Carlo method, John Metropolis and Stanislaw Ulam +recognized that solving the Boltzmann equation with the Monte Carlo method could +be done in parallel very easily whereas the deterministic counterparts for +solving the Boltzmann equation did not offer such a natural means of +parallelism. With the introduction of `vector computers`_ in the early 1970s, +general-purpose parallel computing became a reality. In 1972, Troubetzkoy et +al. designed a Monte Carlo code to be run on the first vector computer, the +ILLIAC-IV [Troubetzkoy]_. The general principles from that work were later +refined and extended greatly through the `work of Forrest Brown`_ in the +1980s. However, as Brown's work shows, the `single-instruction multiple-data`_ +(SIMD) parallel model inherent to vector processing does not lend itself to the +parallelism on particles in Monte Carlo simulations. Troubetzkoy et +al. recognized this, remarking that "the order and the nature of these physical +events have little, if any, correlation from history to history," and thus +following independent particle histories simultaneously using a SIMD model is +difficult. + +The difficulties with vector processing of Monte Carlo codes led to the adoption +of the `single program multiple data`_ (SPMD) technique for parallelization. In +this model, each different process tracks a particle independently of other +processes, and between fission source generations the processes communicate data +through a `message-passing interface`_. This means of parallelism was enabled by +the introduction of message-passing standards in the late 1980s and early 1990s +such as PVM_ and MPI_. The SPMD model proved much easier to use in practice and +took advantage of the inherent parallelism on particles rather than +instruction-level parallelism. As a result, it has since become ubiquitous for +Monte Carlo simulations of transport phenomena. + +Thanks to the particle-level parallelism using SPMD techniques, extremely high +parallel efficiencies could be achieved in Monte Carlo codes. Until the last +decade, even the most demanding problems did not require transmitting large +amounts of data between processors, and thus the total amount of time spent on +communication was not significant compared to the amount of time spent on +computation. However, today's computing power has created a demand for +increasingly large and complex problems, requiring a greater number of particles +to obtain decent statistics (and convergence in the case of criticality +calculations). This results in a correspondingly higher amount of communication, +potentially degrading the parallel efficiency. Thus, while Monte Carlo +simulations may seem `embarrassingly parallel`_, obtaining good parallel scaling +with large numbers of processors can be quite difficult to achieve in practice. + +.. _fission-bank-algorithms: + +----------------------- +Fission Bank Algorithms +----------------------- + +Master-Slave Algorithm +---------------------- + +Monte Carlo particle transport codes commonly implement a SPMD model by having +one master process that controls the scheduling of work and the remaining +processes wait to receive work from the master, process the work, and then send +their results to the master at the end of the simulation (or a source iteration +in the case of an eigenvalue calculation). This idea is illustrated in +:ref:`Figure 1 `. + +.. _figure-master-slave: + +.. figure:: ../../img/master-slave.png + :align: center + :figclass: align-center + + **Figure 1**: Communication pattern in master-slave algorithm. + +Eigenvalue calculations are slightly more difficult to parallelize than fixed +source calculations since it is necessary to converge on the fission source +distribution and eigenvalue before tallying. In the +:ref:`method-successive-generations`, to ensure that the results are +reproducible, one must guarantee that the process by which fission sites are +randomly sampled does not depend on the number of processors. What is typically +done is the following: + + 1. Each compute node sends_ its fission bank sites to a master process; + + 2. The master process sorts or orders the fission sites based on a unique + identifier; + + 3. The master process samples :math:`N` fission sites from the ordered array + of :math:`M` sites; and + + 4. The master process broadcasts_ all the fission sites to the compute + nodes. + +The first and last steps of this process are the major sources of communication +overhead between cycles. Since the master process must receive :math:`M` fission +sites from the compute nodes, the first step is necessarily serial. This step +can be completed in :math:`O(M)` time. The broadcast step can benefit from +parallelization through a tree-based algorithm. Despite this, the communication +overhead is still considerable. + +To see why this is the case, it is instructive to look at a hypothetical +example. Suppose that a calculation is run with :math:`N = 10,000,000` neutrons +across 64 compute nodes. On average, :math:`M = 10,000,000` fission sites will +be produced. If the data for each fission site consists of a spatial location +(three 8 byte real numbers) and a unique identifier (one 4 byte integer), the +memory required per site is 28 bytes. To broadcast 10,000,000 source sites to 64 +nodes will thus require transferring 17.92 GB of data. Since each compute node +does not need to keep every source site in memory, one could modify the +algorithm from a broadcast to a scatter_. However, for practical reasons +(e.g. work self-scheduling), this is normally not done in production Monte Carlo +codes. + +.. _nearest-neighbors-algorithm: + +Nearest Neighbors Algorithm +--------------------------- + +To reduce the amount of communication required in a fission bank synchronization +algorithm, it is desirable to move away from the typical master-slave algorithm +to an algorithm whereby the compute nodes communicate with one another only as +needed. This concept is illustrated in :ref:`Figure 2 +`. + +.. _figure-nearest-neighbor: + +.. figure:: ../../img/nearest-neighbor.png + :align: center + :figclass: align-center + + **Figure 2**: Communication pattern in nearest neighbor algorithm. + +Since the source sites for each cycle are sampled from the fission sites banked +from the previous cycle, it is a common occurrence for a fission site to be +banked on one compute node and sent back to the master only to get sent back to +the same compute node as a source site. As a result, much of the communication +inherent in the algorithm described previously is entirely unnecessary. By +keeping the fission sites local, having each compute node sample fission sites, +and sending sites between nodes only as needed, one can cut down on most of the +communication. One algorithm to achieve this is as follows: + + 1. An exclusive scan is performed on the number of sites banked, and the + total number of fission bank sites is broadcasted to all compute nodes. By + picturing the fission bank as one large array distributed across multiple + nodes, one can see that this step enables each compute node to determine the + starting index of fission bank sites in this array. Let us call the starting + and ending indices on the :math:`i`-th node :math:`a_i` and :math:`b_i`, + respectively; + + 2. Each compute node samples sites at random from the fission bank using the + same starting seed. A separate array on each compute node is created that + consists of sites that were sampled local to that node, i.e. if the index of + the sampled site is between :math:`a_i` and :math:`b_i`, it is set aside; + + 3. If any node sampled more than :math:`N/p` fission sites where :math:`p` + is the number of compute nodes, the extra sites are put in a separate array + and sent to all other compute nodes. This can be done efficiently using the + allgather_ collective operation; + + 4. The extra sites are divided among those compute nodes that sampled fewer + than :math:`N/p` fission sites. + +However, even this algorithm exhibits more communication than necessary since +the allgather will send fission bank sites to nodes that don't necessarily +need any extra sites. + +One alternative is to replace the allgather with a series of sends. If +:math:`a_i` is less than :math:`iN/p`, then send :math:`iN/p - a_i` sites to the +left adjacent node. Similarly, if :math:`a_i` is greater than :math:`iN/p`, then +receive :math:`a_i - iN/p` from the left adjacent node. This idea is applied to +the fission bank sites at the end of each node's array as well. If :math:`b_i` +is less than :math:`(i+1)N/p`, then receive :math:`(i+1)N/p - b_i` sites from +the right adjacent node. If :math:`b_i` is greater than :math:`(i+1)N/p`, then +send :math:`b_i - (i+1)N/p` sites to the right adjacent node. Thus, each compute +node sends/receives only two messages under normal circumstances. + +The following example illustrates how this algorithm works. Let us suppose we +are simulating :math:`N = 1000` neutrons across four compute nodes. For this +example, it is instructive to look at the state of the fission bank and source +bank at several points in the algorithm: + + 1. The beginning of a cycle where each node has :math:`N/p` source sites; + + 2. The end of a cycle where each node has accumulated fission sites; + + 3. After sampling, where each node has some amount of source sites usually + not equal to :math:`N/p`; + + 4. After redistribution, each node again has :math:`N/p` source sites for + the next cycle; + +At the end of each cycle, each compute node needs 250 fission bank sites to +continue on the next cycle. Let us suppose that :math:`p_0` produces 270 fission +banks sites, :math:`p_1` produces 230, :math:`p_2` produces 290, and :math:`p_3` +produces 250. After each node samples from its fission bank sites, let's assume +that :math:`p_0` has 260 source sites, :math:`p_1` has 215, :math:`p_2` has 280, +and :math:`p_3` has 245. Note that the total number of sampled sites is 1000 as +needed. For each node to have the same number of source sites, :math:`p_0` needs +to send its right-most 10 sites to :math:`p_1`, and :math:`p_2` needs to send +its left-most 25 sites to :math:`p_1` and its right-most 5 sites to +:math:`p_3`. A schematic of this example is shown in :ref:`Figure 3 +`. The data local to each node is given a different +hatching, and the cross-hatched regions represent source sites that are +communicated between adjacent nodes. + +.. _figure-neighbor-example: + +.. figure:: ../../img/nearest-neighbor-example.png + :align: center + :figclass: align-center + + **Figure 3**: Example of nearest neighbor algorithm. + +.. _master-slave-cost: + +Cost of Master-Slave Algorithm +------------------------------ + +While the prior considerations may make it readily apparent that the novel +algorithm should outperform the traditional algorithm, it is instructive to look +at the total communication cost of the novel algorithm relative to the +traditional algorithm. This is especially so because the novel algorithm does +not have a constant communication cost due to stochastic fluctuations. Let us +begin by looking at the cost of communication in the traditional algorithm + +As discussed earlier, the traditional algorithm is composed of a series of sends +and typically a broadcast. To estimate the communication cost of the algorithm, +we can apply a simple model that captures the essential features. In this model, +we assume that the time that it takes to send a message between two nodes is +given by :math:`\alpha + (sN)\beta`, where :math:`\alpha` is the time it takes +to initiate the communication (commonly called the latency_), :math:`\beta` is +the transfer time per unit of data (commonly called the bandwidth_), :math:`N` +is the number of fission sites, and :math:`s` is the size in bytes of each +fission site. + +The first step of the traditional algorithm is to send :math:`p` messages to the +master node, each of size :math:`sN/p`. Thus, the total time to send these +messages is + +.. math:: + :label: t-send + + t_{\text{send}} = p\alpha + sN\beta. + +Generally, the best parallel performance is achieved in a weak scaling scheme +where the total number of histories is proportional to the number of +processors. However, we see that when :math:`N` is proportional to :math:`p`, +the time to send these messages increases proportionally with :math:`p`. + +Estimating the time of the broadcast is complicated by the fact that different +MPI implementations may use different algorithms to perform collective +communications. Worse yet, a single implementation may use a different algorithm +depending on how many nodes are communicating and the size of the message. Using +multiple algorithms allows one to minimize latency for small messages and +minimize bandwidth for long messages. + +We will focus here on the implementation of broadcast in the MPICH2_ +implementation. For short messages, MPICH2 uses a `binomial tree`_ algorithm. In +this algorithm, the root process sends the data to one node in the first step, +and then in the subsequent, both the root and the other node can send the data +to other nodes. Thus, it takes a total of :math:`\lceil \log_2 p \rceil` steps +to complete the communication. The time to complete the communication is + +.. math:: + :label: t-short + + t_{\text{short}} = \lceil \log_2 p \rceil \left ( \alpha + sN\beta \right ). + +This algorithm works well for short messages since the latency term scales +logarithmically with the number of nodes. However, for long messages, an +algorithm that has lower bandwidth has been proposed by Barnett_ and implemented +in MPICH2. Rather than using a binomial tree, the broadcast is divided into a +scatter and an allgather. The time to complete the scatter is :math:` \log_2 p +\: \alpha + \frac{p-1}{p} N\beta` using a binomial tree algorithm. The allgather +is performed using a ring algorithm that completes in :math:`p-1) \alpha + +\frac{p-1}{p} N\beta`. Thus, together the time to complete the broadcast is + +.. math:: + :label: t-broadcast + + t_{\text{long}} = \left ( \log_2 p + p - 1 \right ) \alpha + 2 \frac{p-1}{p} + sN\beta. + +The fission bank data will generally exceed the threshold for switching from +short to long messages (typically 8 kilobytes), and thus we will use the +equation for long messages. Adding equations :eq:`t-send` and :eq:`t-broadcast`, +the total cost of the series of sends and the broadcast is + +.. math:: + :label: t-old + + t_{\text{old}} = \left ( \log_2 p + 2p - 1 \right ) \alpha + \frac{3p-2}{p} + sN\beta. + +Cost of Nearest Neighbor Algorithm +---------------------------------- + +With the communication cost of the traditional fission bank algorithm +quantified, we now proceed to discuss the communicatin cost of the proposed +algorithm. Comparing the cost of communication of this algorithm with the +traditional algorithm is not trivial due to fact that the cost will be a +function of how many fission sites are sampled on each node. If each node +samples exactly :math:`N/p` sites, there will not be communication between nodes +at all. However, if any one node samples more or less than :math:`N/p` sites, +the deviation will result in communication between logically adjacent nodes. To +determine the expected deviation, one can analyze the process based on the +fundamentals of the Monte Carlo process. + +The steady-state neutron transport equation for a multiplying medium can be +written in the form of an eigenvalue problem, + +.. math:: + :label: NTE + + S(\mathbf{r})= \frac{1}{k} \int F(\mathbf{r}' \rightarrow + \mathbf{r})S(\mathbf{r}')\: d\mathbf{r}, + +where :math:`\mathbf{r}` is the spatial coordinates of the neutron, +:math:`S(\mathbf{r})` is the source distribution defined as the expected number +of neutrons born from fission per unit phase-space volume at :math:`\mathbf{r}`, +:math:`F( \mathbf{r}' \rightarrow \mathbf{r})` is the expected number of +neutrons born from fission per unit phase space volume at :math:`\mathbf{r}` +caused by a neutron at :math:`\mathbf{r}`, and :math:`k` is the eigenvalue. The +fundamental eigenvalue of equation :eq:`NTE` is known as :math:`k_{eff}`, but +for simplicity we will simply refer to it as :math:`k`. + +In a Monte Carlo criticality simulation, the power iteration method is applied +iteratively to obtain stochastic realizations of the source distribution and +estimates of the :math:`k`-eigenvalue. Let us define :math:`\hat{S}^{(m)}` to be +the realization of the source distribution at cycle :math:`m` and +:math:`\hat{\epsilon}^{(m)}` be the noise arising from the stochastic nature of +the tracking process. We can write the stochastic realization in terms of the +fundamental source distribution and the noise component as (see `Brissenden and +Garlick`_): + +.. math:: + :label: source + + \hat{S}^{(m)}(\mathbf{r})= N S(\mathbf{r}) + \sqrt{N} + \hat{\epsilon}^{(m)}(\mathbf{r}), + +where :math:`N` is the number of particle histories per cycle. Without loss of +generality, we shall drop the superscript notation indicating the cycle as it is +understood that the stochastic realization is at a particular cycle. The +expected value of the stochastic source distribution is simply + +.. math:: + :label: expected-value-source + + E \left[ \hat{S}(\mathbf{r})\right] = N S (\mathbf{r}) + +since :math:`E \left[ \hat{\epsilon}(\mathbf{r})\right] = 0`. The noise in the +source distribution is due only to :math:`\hat{\epsilon}(\mathbf{r})` and thus +the variance of the source distribution will be + +.. math:: + :label: var-source + + \text{Var} \left[ \hat{S}(\mathbf{r})\right] = N \text{Var} \left[ + \hat{\epsilon}(\mathbf{r}) \right]. + +Lastly, the stochastic and true eigenvalues can be written as integrals over all +phase space of the stochastic and true source distributions, respectively, as + +.. math:: + :label: k-to-source + + \hat{k} = \frac{1}{N} \int \hat{S}(\mathbf{r}) \: d\mathbf{r} \quad + \text{and} \quad k = \int S(\mathbf{r}) \: d\mathbf{r}, + +noting that :math:`S(\mathbf{r})` is :math:`O(1)`. One should note that the +expected value :math:`k` calculated by Monte Carlo power iteration (i.e. the +method of successive generations) will be biased from the true fundamental +eigenvalue of equation :eq:`NTE` by :math:`O(1/N)` (see `Brissenden and +Garlick`_), but we will assume henceforth that the number of particle histories +per cycle is sufficiently large to neglect this bias. + +With this formalism, we now have a framework within which we can determine the +properties of the distribution of expected number of fission sites. The explicit +form of the source distribution can be written as + +.. math:: + :label: source-explicit + + \hat{S}(\mathbf{r}) = \sum_{i=1}^{M} w_i \delta( \mathbf{r} - \mathbf{r}_i ) + +where :math:`\mathbf{r}_i` is the spatial location of the :math:`i`-th fission +site, :math:`w_i` is the statistical weight of the fission site at +:math:`\mathbf{r}_i`, and :math:`M` is the total number of fission sites. It is +clear that the total weight of the fission sites is simply the integral of the +source distribution. Integrating equation :eq:`source` over all space, we obtain + +.. math:: + :label: source-integrated + + \int \hat{S}(\mathbf{r}) \: d\mathbf{r} = N \int S(\mathbf{r}) \: + d\mathbf{r} + \sqrt{N} \int \hat{\epsilon}(\mathbf{r}) \: d\mathbf{r} . + +Substituting the expressions for the stochastic and true eigenvalues from +equation :eq:`k-to-source`, we can relate the stochastic eigenvalue to the +integral of the noise component of the source distribution as + +.. math:: + :label: noise-integeral + + N\hat{k} = Nk + \sqrt{N} \int \hat{\epsilon}(\mathbf{r}) \: d\mathbf{r}. + +Since the expected value of :math:`\hat{\epsilon}` is zero, the expected value +of its integral will also be zero. We thus see that the variance of the integral +of the source distribution, i.e. the variance of the total weight of fission +sites produced, is directly proportional to the variance of the integral of the +noise component. Let us call this term :math:`\sigma^2` for simplicity: + +.. math:: + :label: variance-sigma2 + + \text{Var} \left[ \int \hat{S}(\mathbf{r}) \right ] = N \sigma^2. + +The actual value of :math:`\sigma^2` will depend on the physical nature of the +problem, whether variance reduction techniques are employed, etc. For instance, +one could surmise that for a highly scattering problem, :math:`\sigma^2` would +be smaller than for a highly absorbing problem since more collisions will lead +to a more precise estimate of the source distribution. Similarly, using implicit +capture should in theory reduce the value of :math:`\sigma^2`. + +Let us now consider the case where the :math:`N` total histories are divided up +evenly across :math:`p` compute nodes. Since each node simulates :math:`N/p` +histories, we can write the source distribution as + +.. math:: + :label: source-node + + \hat{S}_i(\mathbf{r})= \frac{N}{p} S(\mathbf{r}) + \sqrt{\frac{N}{p}} + \hat{\epsilon}_i(\mathbf{r}) \quad \text{for} \quad i = 1, \dots, p + +Integrating over all space and simplifying, we can obtain an expression for the +eigenvalue on the :math:`i`-th node: + +.. math:: + :label: k-i-hat + + \hat{k}_i = k + \sqrt{\frac{p}{N}} \int \hat{\epsilon}_i(\mathbf{r}) \: + d\mathbf{r}. + +It is easy to show from this expression that the stochastic realization of the +global eigenvalue is merely the average of these local eigenvalues: + +.. math:: + :label: average-k-as-sum + + \hat{k} = \frac{1}{p} \sum_{i=1}^p \hat{k}_i. + +As was mentioned earlier, at the end of each cycle one must sample :math:`N` +sites from the :math:`M` sites that were created. Thus, the source for the next +cycle can be seen as the fission source from the current cycle divided by the +stochastic realization of the eigenvalue since it is clear from equation +:eq:`k-to-source` that :math:`\hat{k} = M/N`. Similarly, the number of sites +sampled on each compute node that will be used for the next cycle is + +.. math:: + :label: sites-per-node + + M_i = \frac{1}{\hat{k}} \int \hat{S}_i(\mathbf{r}) \: d\mathbf{r} = + \frac{N}{p} \frac{\hat{k}_i}{\hat{k}}. + +While we know conceptually that each compute node will under normal +circumstances send two messages, many of these messages will overlap. Rather +than trying to determine the actual communication cost, we will instead attempt +to determine the maximum amount of data being communicated from one node to +another. At any given cycle, the number of fission sites that the :math:`j`-th +compute node will send or receive (:math:`\Lambda_j`) is + +.. math:: + :label: Lambda + + \Lambda_j = \left | \sum_{i=1}^j M_i - \frac{jN}{p} \right |. + +Noting that :math:`jN/p` is the expected value of the summation, we can write +the expected value of :math:`\Lambda_j` as the mean absolute deviation of the +summation: + +.. math:: + :label: mean-dev-lambda + + E \left [ \Lambda_j \right ] = E \left [ \left | \sum_{i=1}^j M_i - + \frac{jN}{p} \right | \right ] = \text{MD} \left [ \sum_{i=1}^j M_i \right ] + +where :math:`\text{MD}` indicates the mean absolute deviation of a random +variable. The mean absolute deviation is an alternative measure of variability. + +In order to ascertain any information about the mean deviation of :math:`M_i`, +we need to know the nature of its distribution. Thus far, we have said nothing +of the distributions of the random variables in question. The total number of +fission sites resulting from the tracking of :math:`N` neutrons can be shown to +be normally distributed via the :ref:`central-limit-theorem` (provided that +:math:`N` is sufficiently large) since the fission sites resulting from each +neutron are "sampled" from independent, identically-distributed random +variables. Thus, :math:`\hat{k}` and :math:`\int \hat{S} (\mathbf{r}) \: +d\mathbf{r}` will be normally distributed as will the individual estimates of +these on each compute node. + +Next, we need to know what the distribution of :math:`M_i` in equation +:eq:`sites-per-node` is or, equivalently, how :math:`\hat{k}_i / \hat{k}` is +distributed. The distribution of a ratio of random variables is not easy to +calculate analytically, and it is not guaranteed that the ratio distribution is +normal if the numerator and denominator are normally distributed. For example, +if :math:`X` is a standard normal distribution and :math:`Y` is also standard +normal distribution, then the ratio :math:`X/Y` has the standard `Cauchy +distribution`_. The reader should be reminded that the Cauchy distribution has +no defined mean or variance. That being said, Geary_ has shown that, for the +case of two normal distributions, if the denominator is unlikely to assume +values less than zero, then the ratio distribution is indeed approximately +normal. In our case, :math:`\hat{k}` absolutely cannot assume a value less than +zero, so we can be reasonably assured that the distribution of :math:`M_i` will +be normal. + +For a normal distribution with mean :math:`\mu` and distribution function +:math:`f(x)`, it can be shown that + +.. math:: + :label: mean-dev-to-stdev + + \int_{-\infty}^{\infty} f(x) \left | x - \mu \right | \: dx = + \sqrt{\frac{2}{\pi} \int_{-\infty}^{\infty} f(x) \left ( x - \mu \right )^2 + \: dx} + +and thus the mean absolute deviation is :math:`\sqrt{2/\pi}` times the standard +deviation. Therefore, to evaluate the mean absolute deviation of :math:`M_i`, we +need to first determine its variance. Substituting equation +:eq:`average-k-as-sum` into equation :eq:`sites-per-node`, we can rewrite +:math:`M_i` solely in terms of :math:`\hat{k}_1, \dots, \hat{k}_p`: + +.. math:: + :label: M-i + + M_i = \frac{N \hat{k}_i}{\sum\limits_{j=1}^p \hat{k}_j}. + +Since we know the variance of :math:`\hat{k}_i`, we can use the error +propagation law to determine the variance of :math:`M_i`: + +.. math:: + :label: M-variance + + \text{Var} \left [ M_i \right ] = \sum_{j=1}^p \left ( \frac{\partial + M_i}{\partial \hat{k}_j} \right )^2 \text{Var} \left [ \hat{k}_j \right ] + + \sum\limits_{j \neq m} \sum\limits_{m=1}^p \left ( \frac{\partial + M_i}{\partial \hat{k}_j} \right ) \left ( \frac{\partial M_i}{\partial + \hat{k}_m} \right ) \text{Cov} \left [ \hat{k}_j, \hat{k}_m \right ] + +where the partial derivatives are evaluated at :math:`\hat{k}_j = k`. Since +:math:`\hat{k}_j` and :math:`\hat{k}_m` are independent if :math:`j \neq m`, +their covariance is zero and thus the second term cancels out. Evaluating the +partial derivatives, we obtain + +.. math:: + :label: M-variance-2 + + \text{Var} \left [ M_i \right ] = \left ( \frac{N(p-1)}{kp^2} \right )^2 + \frac{p\sigma^2}{N} + \sum_{j \neq i} \left ( \frac{-N}{kp^2} \right )^2 + \frac{p\sigma^2}{N} = \frac{N(p-1)}{k^2p^2} \sigma^2. + +Through a similar analysis, one can show that the variance of +:math:`\sum_{i=1}^j M_i` is + +.. math:: + :label: sum-M-variance + + \text{Var} \left [ \sum_{i=1}^j M_i \right ] = \frac{Nj(p-j)}{k^2p^2} + \sigma^2 + +Thus, the expected amount of communication on node :math:`j`, i.e. the mean +absolute deviation of :math:`\sum_{i=1}^j M_i` is proportional to + +.. math:: + :label: communication-cost + + E \left [ \Lambda_j \right ] = \sqrt{\frac{2Nj(p-j)\sigma^2}{\pi k^2p^2}}. + +This formula has all the properties that one would expect based on intuition: + + 1. As the number of histories increases, the communication cost on each node + increases as well; + + 2. If :math:`p=1`, i.e. if the problem is run on only one compute node, the + variance will be zero. This reflects the fact that exactly :math:`N` sites + will be sampled if there is only one node. + + 3. For :math:`j=p`, the variance will be zero. Again, this says that when + you sum the number of sites from each node, you will get exactly :math:`N` + sites. + +We can determine the node that has the highest communication cost by +differentiating equation :eq:`communication-cost` with respect to :math:`j`, +setting it equal to zero, and solving for :math:`j`. Doing so yields +:math:`j_{\text{max}} = p/2`. Interestingly, substituting :math:`j = p/2` in +equation :eq:`communication-cost` shows us that the maximum communication cost +is actually independent of the number of nodes: + +.. math:: + :label: maximum-communication + + E \left [ \Lambda_{j_{\text{max}}} \right ] = \sqrt{ \frac{N\sigma^2}{2\pi + k^2}}. + +---------- +References +---------- + +.. [Troubetzkoy] E. Troubetzkoy, H. Steinberg, and M. Kalos, "Monte Carlo + Radiation Penetration Calculations on a Parallel Computer," + *Trans. Am. Nucl. Soc.*, **17**, 260 (1973). + +.. _first paper: http://www.jstor.org/stable/2280232 + +.. _work of Forrest Brown: http://hdl.handle.net/2027.42/24996 + +.. _Brissenden and Garlick: http://dx.doi.org/10.1016/0306-4549(86)90095-2 + +.. _MPICH2: http://www.mcs.anl.gov/mpi/mpich + +.. _binomial tree: http://www.cs.auckland.ac.nz/~jmor159/PLDS210/trees.html + +.. _Geary: http://www.jstor.org/stable/10.2307/2342070 + +.. _Barnett: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.51.7772 + +.. _single-instruction multiple-data: http://en.wikipedia.org/wiki/SIMD + +.. _vector computers: http://en.wikipedia.org/wiki/Vector_processor + +.. _single program multiple data: http://en.wikipedia.org/wiki/SPMD + +.. _message-passing interface: http://en.wikipedia.org/wiki/Message_Passing_Interface + +.. _PVM: http://www.csm.ornl.gov/pvm/pvm_home.html + +.. _MPI: http://www.mcs.anl.gov/research/projects/mpi/ + +.. _embarrassingly parallel: http://en.wikipedia.org/wiki/Embarrassingly_parallel + +.. _sends: http://www.mcs.anl.gov/research/projects/mpi/www/www3/MPI_Send.html + +.. _broadcasts: http://www.mcs.anl.gov/research/projects/mpi/www/www3/MPI_Bcast.html + +.. _scatter: http://www.mcs.anl.gov/research/projects/mpi/www/www3/MPI_Scatter.html + +.. _allgather: http://www.mcs.anl.gov/research/projects/mpi/www/www3/MPI_Allgather.html + +.. _Cauchy distribution: http://en.wikipedia.org/wiki/Cauchy_distribution + +.. _latency: http://en.wikipedia.org/wiki/Latency_(engineering)#Packet-switched_networks + +.. _bandwidth: http://en.wikipedia.org/wiki/Bandwidth_(computing) diff --git a/docs/source/methods/statistics.rst b/docs/source/methods/statistics.rst index 19cbf7320..387c9a675 100644 --- a/docs/source/methods/statistics.rst +++ b/docs/source/methods/statistics.rst @@ -33,6 +33,8 @@ X_n}{n}` `converges in probability`_ to the true mean, i.e. for all \lim\limits_{n\rightarrow\infty} P \left ( \left | \bar{X}_n - \mu \right | \ge \epsilon \right ) = 0. +.. _central-limit-theorem: + --------------------- Central Limit Theorem ---------------------