From b4eeb55e29b5a6d60588d639ee97c9f6f6ccd09a Mon Sep 17 00:00:00 2001 From: Nicholas Clark Date: Thu, 16 Nov 2023 09:02:33 +1000 Subject: [PATCH] bug in summary when multiple random effects are used --- R/add_MACor.R | 90 ++++++++++++++++++++++++++++---------- R/summary.mvgam.R | 4 +- src/mvgam.dll | Bin 1064960 -> 1064960 bytes tests/testthat/Rplots.pdf | Bin 22890 -> 22884 bytes 4 files changed, 68 insertions(+), 26 deletions(-) diff --git a/R/add_MACor.R b/R/add_MACor.R index 984db65d..10ff6a23 100644 --- a/R/add_MACor.R +++ b/R/add_MACor.R @@ -294,10 +294,17 @@ add_MaCor = function(model_file, if(add_cor){ if(trend_model %in% c('AR1', 'RW')){ - model_file[max(grep('= b_raw[', - model_file, fixed = TRUE))] <- - paste0(model_file[max(grep('= b_raw[', - model_file, fixed = TRUE))], + if(any(grepl('= mu_raw[', + model_file, fixed = TRUE))){ + insert_line <- max(grep('= mu_raw[', + model_file, fixed = TRUE)) + } else if(any(grepl('= b_raw[', + model_file, fixed = TRUE))){ + insert_line <- max(grep('= b_raw[', + model_file, fixed = TRUE)) + } + model_file[insert_line] <- + paste0(model_file[insert_line], '\n// derived latent states\n', 'trend_raw[1] = ', if(drift){ 'drift + '} else {NULL}, @@ -328,10 +335,17 @@ add_MaCor = function(model_file, } if(trend_model == 'AR2'){ - model_file[max(grep('= b_raw[', - model_file, fixed = TRUE))] <- - paste0(model_file[max(grep('= b_raw[', - model_file, fixed = TRUE))], + if(any(grepl('= mu_raw[', + model_file, fixed = TRUE))){ + insert_line <- max(grep('= mu_raw[', + model_file, fixed = TRUE)) + } else if(any(grepl('= b_raw[', + model_file, fixed = TRUE))){ + insert_line <- max(grep('= b_raw[', + model_file, fixed = TRUE)) + } + model_file[insert_line] <- + paste0(model_file[insert_line], '\n// derived latent states\n', 'trend_raw[1] = ', if(drift){ 'drift + '} else {NULL}, @@ -371,10 +385,17 @@ add_MaCor = function(model_file, } if(trend_model == 'AR3'){ - model_file[max(grep('= b_raw[', - model_file, fixed = TRUE))] <- - paste0(model_file[max(grep('= b_raw[', - model_file, fixed = TRUE))], + if(any(grepl('= mu_raw[', + model_file, fixed = TRUE))){ + insert_line <- max(grep('= mu_raw[', + model_file, fixed = TRUE)) + } else if(any(grepl('= b_raw[', + model_file, fixed = TRUE))){ + insert_line <- max(grep('= b_raw[', + model_file, fixed = TRUE)) + } + model_file[insert_line] <- + paste0(model_file[insert_line], '\n// derived latent states\n', 'trend_raw[1] = ', if(drift){ 'drift + '} else {NULL}, @@ -426,10 +447,17 @@ add_MaCor = function(model_file, } else { if(trend_model %in% c('AR1', 'RW')){ - model_file[max(grep('= b_raw[', - model_file, fixed = TRUE))] <- - paste0(model_file[max(grep('= b_raw[', - model_file, fixed = TRUE))], + if(any(grepl('= mu_raw[', + model_file, fixed = TRUE))){ + insert_line <- max(grep('= mu_raw[', + model_file, fixed = TRUE)) + } else if(any(grepl('= b_raw[', + model_file, fixed = TRUE))){ + insert_line <- max(grep('= b_raw[', + model_file, fixed = TRUE)) + } + model_file[insert_line] <- + paste0(model_file[insert_line], 'for(j in 1:n_series){\n', 'trend[1, j] = ', if(drift){ 'drift[j] + '} else {NULL}, @@ -448,10 +476,17 @@ add_MaCor = function(model_file, } if(trend_model == 'AR2'){ - model_file[max(grep('= b_raw[', - model_file, fixed = TRUE))] <- - paste0(model_file[max(grep('= b_raw[', - model_file, fixed = TRUE))], + if(any(grepl('= mu_raw[', + model_file, fixed = TRUE))){ + insert_line <- max(grep('= mu_raw[', + model_file, fixed = TRUE)) + } else if(any(grepl('= b_raw[', + model_file, fixed = TRUE))){ + insert_line <- max(grep('= b_raw[', + model_file, fixed = TRUE)) + } + model_file[insert_line] <- + paste0(model_file[insert_line], 'for(j in 1:n_series){\n', 'trend[1, j] = ', if(drift){ 'drift[j] + '} else {NULL}, @@ -475,10 +510,17 @@ add_MaCor = function(model_file, } if(trend_model == 'AR3'){ - model_file[max(grep('= b_raw[', - model_file, fixed = TRUE))] <- - paste0(model_file[max(grep('= b_raw[', - model_file, fixed = TRUE))], + if(any(grepl('= mu_raw[', + model_file, fixed = TRUE))){ + insert_line <- max(grep('= mu_raw[', + model_file, fixed = TRUE)) + } else if(any(grepl('= b_raw[', + model_file, fixed = TRUE))){ + insert_line <- max(grep('= b_raw[', + model_file, fixed = TRUE)) + } + model_file[insert_line] <- + paste0(model_file[insert_line], 'for(j in 1:n_series){\n', 'trend[1, j] = ', if(drift){ 'drift[j] + '} else {NULL}, diff --git a/R/summary.mvgam.R b/R/summary.mvgam.R index 34b8d30f..67c6acc4 100644 --- a/R/summary.mvgam.R +++ b/R/summary.mvgam.R @@ -197,7 +197,7 @@ if(all(is.na(object$sp_names))){ } else { if(any(unlist(purrr::map(object$mgcv_model$smooth, inherits, 'random.effect')))){ - re_labs <- unlist(lapply(purrr::map(object$mgcv_model$smooth, 'term'), + re_labs <- unlist(lapply(purrr::map(object$mgcv_model$smooth, 'label'), paste, collapse = ','))[ unlist(purrr::map(object$mgcv_model$smooth, inherits, 'random.effect'))] re_sds <- mcmc_summary(object$model_output, 'sigma_raw', @@ -601,7 +601,7 @@ if(!is.null(object$trend_call)){ } else { if(any(unlist(purrr::map(object$trend_mgcv_model$smooth, inherits, 'random.effect')))){ - re_labs <- unlist(lapply(purrr::map(object$trend_mgcv_model$smooth, 'term'), + re_labs <- unlist(lapply(purrr::map(object$trend_mgcv_model$smooth, 'label'), paste, collapse = ','))[ unlist(purrr::map(object$trend_mgcv_model$smooth, inherits, 'random.effect'))] re_labs <- gsub('series', 'trend', re_labs) diff --git a/src/mvgam.dll b/src/mvgam.dll index 2161ca821a5edcb6456981b031a7acf4b0321970..a8ea964988923fdde9d8056c6717ba9b91155529 100644 GIT binary patch delta 78 zcmV-U0I~mofI)zOL6C?8kwjIIi%kPb6cK^jh1&sz+X98#1BKfJh1&&%+XjW(2Zh@R kh1&^*+X{u-3x(SZh1(5<+YYtc4=LpXkVI9t2j&m?Su9Z^g#Z8m delta 78 zcmV~$y$wK66b8_6as9J~G2Q^GSJ~z&zC@y-lGsCpRlb@Xbj~qzW^S3}R=>`@YsK+m-C-mm3rUELd5ctu)`6R)aK~3Cm6YjG zu)#$%tgwJ}I)oIqF)>+&Bt|YFfba+&3j!!V80F#zBi9BDFj+r(m`oS1BR_m3{J2$- zbJ|4ohP;`80TyklU?>Ru7{|X*ILQy3D9V4eiR*uE6vp)J=CiS?%iSdOS5qAxZdP?X zo{^7zL|WX9x1tjI*_e*U#q8asxUaq~C!Gs{vTu!5tYl|Doey8sS!J4ZDO?K{6^36q zdQR<9a9aEtuAlRM!ce3C65x?l-8B<*`59&4!X#G1JeC$wY*mM_hv02=L_ApSsiy%spD%fIMx1&t z`7zhRr$6uL+CyqZ(5`+emSCM_p+~jL-^F}Ltu5G4AgA(deD!*0nplztt3x92j#`8y z;xpTbiC|p*vPL7*r`A(6&mnrEyEf?Dda^ny`39myaA5(MRV_doG}aL*Av zlEPA!7LmStb-w`-8*cb?3%19fKU`UDL{aH3i-sNw?Q@uBQOi3`g66$%0)e0QK#wCPvKFz#WVWS;M zoPySfS8V#XY62FRe1r`%)5hW0zL-9{=2S@YM}gJ5GSsa?=ar5i|RP^Ixc7}QdF(H`hVLMOBrHA*>|Cl52lIrnFG1P?M zWcNoRlZ?=w+FI~d3kP1R(7`qvKJygo=(3rN0iC+w!yEg}pWEEy<9uB8*0c343;IX` z*crRCgt!H1*i+n+1AZx0RdZ9gruDn^S7oSU3tjGN-3!+dp3J7qe0v976@7%6>pZ2y z)r*3e?{MgWGQ&$R42l!oI7&0X_^WPHO`v6vNE{p z?s9o{hKfrW0TLKK%V1|nWW|R~HUXZ!L_GV+p|i}8PdaO2OF3ZDddk80a*~@T7AZPe zYhQO~SX_ow$eX8Pm+af07>zDc9)7M2ga1WoXMBb-VR@iMTQ$8?5@f#T;(rmIvdG)8$7p0o z-!%LR3O-;nzKqw|mUrPgu(Ob;}eeeH{}O-!>+FI8^tnQ`=4!jPv?`|X4JBIEh!dwz~0>p_xa#E`V}0}_R5C#7OnM?K{S6w! zBY(fkDUJ`4P)kE1PyjrI6Q2CJ{XK>HN3#u% zF#mD^Rxq{P1G}0VlT;j<(VUqd1`pkHv@8tHRqltXWU3E9!;~9|cspAtv|A$d!VF`@ z`EW;dFzHxZT`)s8tXlIB(-)3F!9$gxVeoskWT=X<0SMk&(+ZT!;LB}RM zhnDa`P=E@Tj`@;QDlf$Fq3(v_k_84gLi2ZTWfqU$)63dgFBRx-9+;o+JZhi zKecORz@p2-kyLz`KD2GUA{smC*m*Lucy^7JTTbFQOb?d zz7@P+Slkst4e%=>Uz0SyCD(VQD5W^jqziwy!$NVsmz6GZ_su1ZrWQE71D$UtjP+Qe z3=e42zeik`Ga~+IVT1y~cxh6q+vj(*tqKf)(y>In%77I+ScLiwd{7LVY{pk}c1U&p zUS~u^%u~eE= zgSH-fP2=+*W^QDir{m}y4qZ@FAikYoSNEcmQU#Bm2BC4RRic1vmg48cIh&Eu7lujf zr@G%}#4{|DABDQ0OkwdhS&aVkTy_QmCCH=T&O+eu&~|q@&J6D)s6+-~_=aYDG*c!X zIL7Sfig7$Q?&*?9;q~1!GWIFoL^tclR7+Ou8dicRtaHBe}Y*0x# z&<hC7ryH(oqWz=JFgc(%Ln?NC%9HtY`g(1rIP>yD- zh=2R>%>s49vOChD{o7h(BHpt;!?@o{3801bLy{d{8+5fJmMBRzsyzAIKSGiT8x<*Q0l-)dpvALIRlJts zmoDuvA3QZk8cNb~QGtY2uZAcake?ueGl}0IB_!& z{W<0M9dwRRbYLTOlf`%K2smlHUaXSSDsLXh<8GK8+3TE9XX2|19ItLoH{;crh@s(f z339e?R>WIRX*yGDTNxKj{JY0hS#I-pR;tma`zn2|0!6gL30Q zquai-bEdo9j>&Ahi^&(7n0eRiKG5&3QM)g=0M%I1`Uvr6?VGHN*Ibe_eTvJW>g<*s3Dd6Wv z{Dbzg>o;%ga~f|D@2q!Lv&C(2-zq!L=!t!cUW4t_>otxxUn264NT6{SP5~2bhzh4D z6?%8EN-s~~Lr0i*UauRidHWp|3H}U%zQ^Q9>=BI2eH*Jj;W0iHiP*abnoQp^3%vHd z^^QpJp;;<3y}v8kzmZW5K)me~urf@a7b$*LO6zv<54-4jB{ib^)4 zExhagw_k)zNNg*#&l7oM8&<33&OVjw|?Jx+AP>veR4(wCwp1!_-W^ShGbE+b)Cder| zHlaiM^g+uhi+9dWeH8NCDL0c>enlM#rmjJllBZ|i`$ges)T)JoBklm6p`ZLH7EB<2m7 z)D=Izl<#Lr)mW3uLRoFaJoScEebOcUh)RlgR?Mvv7a|YCzrTvHdLxBPY09dFTuMGf zBhZxizbU$#2l>`5uu&30nJ9z!jjK)x^(ZrgoQU+AKo`T z-jGt2t#rD=uM*aME=C#6);I?69B~(Xum9|{Y2E2w@t5lwo0hI*mJFCtU7?3OPgI<+@Eh z9 z1!^Jw_bkruhZ^)V!jYBs#So5Vap*7fm)5L3-(Ifq6Ae{!PAce#DG*3uD^gDwUdVRS zyjSk`vuH*gpGs{Vi4Sm_QBlE>@7H_u*7uFsu&Q9s8ox7d$KU+GN*!4&tJ|3|?} zMjOS&>H%itLO5n4fACD>@cnER%QsT9E+7)C4F5)A0c!F9&Iqme|FU37I>%pDyM{P< zz%(~l(kuy$&{p08T;PK~t6LoclG0w;tkc_qZ0r`s76Q!rUm)ZAq_GcaFIeKuSN_gG z!-u_Iki{jeQn;#N`@350!^y{s-X68qdwJ-aK(o)K_qhj}hmzN3i`YV#pu*5U4#@UV zxQW|~)oIF+n{lLA4$8XWr54D?gOeWF6Rc;FvfBX**OB&;iz*S~?bQEbPZeXjPLC z<&2vN_v@N-ce%GC;Wi?SuD(R#`GtRJbAzDhc8l*G?r*!0CMnwzFdHlBWO(NL>961W zcN)nK@P_F%pWXAXTxd5AX?--+;jVHSacqYjjJ}GDa-~1+MHqeHyyfSAF>`1ezwj{0mvd#8^UW`| zPAKdbB3+ShFuT|wUF2c>J=b4EGu&`$5h#`N3+|M?{EKMX$=;CDXZS3M`3Wdkzyo~7 z4%402Y6ln9=KaLr?}tXe9mVZ4hHbI5O+W_SXRE!c*nlp81_9Zf|H;<~rHc^e?!byx ztJ|8}q!6X2gM4<^d+$ywMw~2+KMs7&&sAXl@pVk$zLN;->o-rR@7+aneXd&n6NY44 zW9;iUhA9?t2kIHcPoz%SVwE@S3p1t81-Kqtk3HJ#ZjvE2Bx`s+{{+3S)R?ZPRLuiT zfgE3InrjNupDGT@i<){}LDg&2MCndor1`BF%^Zq{4pJ7BO}+m_I(HAtI37Ow+_UZx zF!xLD*X#oPso4m#xJR-Vtrc>D?a$@Clm`lSCvu;huUj|*S~3r9X~XH%@jbA2`d>8^ zeEw0E2|Jqw$BA-4OZ#fT!M7hXwCy}Kb%omY+;Ql+S3acD*>njRx;OdQh(Z(j+NHyH zyeg#b_$-=7E=Vf5fyiuoEYX*DYgByUK(;mQqqggkGD zn)7WwDNJ+QA&{H@ajt*5`ifdpm~u+Fv?c{82;mkHNBSO&!=E8X(w?;9rf&O%8F z)2XqboNA|g)`iB&um0VhSgu$JH{z>TQVR#X;lZ7^Qw|l5WgYsO`D7X+8=j8X1wo5= zqon1oKK`h6VmB*?#SVQrl}5Xv(xT5LU>Dm;(cJBkF=vT+{Bb#h8=uwbv%-!TAf2qx zLn*~(AT+-z`y^I^EedyNbr^@zBmBcg@2maAn&j&<--;(K>1$Tk)nn+aqW!7%^ii?0 z+^9=HGi8QSwNU>{NS7<3+MGjbn&=v6l)j_2A9+uft``DX#tN@nK*B(gE9y9)gv2BT zxRrX+6^i^KpjoNcHJD&*b>sT$o3dBmP@EA+6&7p1c;G4%S!Nzt1x0+)QQ!Ki)}z~y z;3cPrRy2?XN<#FsoSmK7owX6}a?35|e<~FWs{)CEcEXtEqQcL}g_j2ZH*b=S{*et3 zS*gMMQppBiL$o{ZOn{ucwIVd>;Q!9uM7?1!J|P-gtE4Y{$nQu4{jV(vOPMi5#JQVY zN%3VzB)153Nu8UCbGldCV^D7LX%!DZQS+1IJC7#Gy8Go4FunsjmFZFJB|v|S`s~+mwxb=)?3aXA{3r9gEJ-bbUa4Ew*kjZg zOpsn9c4e2AhKLpw+USsag1t<;E)drzD2cvzd3J8ZrzHW!N3O^lcF*t+Iv?eF_F};6e)FR%AzrF(CVEk>3N?IbcVOZ7 zRB1c7BI>{W5-HJuXPQ2{v>qRLxnn!R+Q;I)x-6-7xG^odYFD~SZ-iUM)kf81rM^aY zZkqfqEyDV0`x^3ZbrI4c81Rz`{A9WFJ49&HVf%%&s8ZWBv5X9$p%+<# zrR&c@L(zaezvc)1M%5Ggd}PZ9oKdm0eNP!q9or{SAmq-0DeElK6oR#u*7n;~3>Ds!j9&1?QF-S0J#5RRahpe4hUe0#-nIOw<_5Hv; z+|#oM17^UJwxkf0vo${Qs{D{$Rh|)y57G2?G%c2&EzX`5Ti=ek7^#l<4R(E(*2jc5 z3>~*3)@=)GQ}DtAAHEFTXp2!dJ~ z?WND|xXnf16o@-p>cUA2|8CS@%;S&a&h1B!T(Spdx-cCF0(@sLefc%3vRq-}yLF>n zZLaUrR(OrX9>2SD=PpKhw~go!O|GIOOhuAtnQt|7^rYM!C|!&EYiI zy}{y=uYi&6$mZXkW!?DYP^lnWz@29P)q0Y?9ov?$1(ORKOU5lDZ1t_73-~pI|C7Jt za@EjkrEZr<$n@aAO#qzMEF$qYfmNd5TcLnGvdeO=EgG>`L5HQm1Zt3&?l@7Rl-A?m zR!F0&VgY8mX|mH(L+~595V3ccNCjp;9+PTaZ=Om;B zXi$I7JwKW4)kSt!3u4Kl^vh+y;G*K-_UFvv zQ#pR}(Y~oZZNOc{(Hn$sS!X?BQZwHEcp6u~5xWV}=s-3-fN-NL&KUu%R^%&RF5c7$ zPnR&%cnj*2+jsO{*s1AUR`uSe|3ffTwz3xebmhSS-?hr;wYP{oq&_9~XKD;vw~e*9 zKXHD2nf~4-PzZz97ZNP!>XKp~8?9`OrlC^Ai2@10Iaj8(ij7}rAaxdh_i=XF7?2z;sZyAZf zfcUz2^gp;|(WpwQ(^WcJhy)Jzgo66wA|wTla>4KIZzTUA0LDi6@4a-`auloudwwfG z5Hg^nA@O6Nfc(WTPA(6Qt^&fg50&ycdj9qC1VAAhxwZVH-4I5`5IMD?LzCVb^Z#DD zAwyp!A@^O;+Jk#%ZK~D-5y!M9j^}v6k(2aL#^30VaxIdCh5sePPiBL@{8HEmB`r;Z zEo4kOCU)0TK8OVh#8nDP_SElNX;VIVnzx ztd{pweDHK?CFeS*c6Wx3QfitSHq32W5_B3>RlHV(4O>QDCPcaJ9Rm{t{QTbnWawpS zo&1)5>@AbC4eDTou+H<7CEt2v$?PTFf9td*ls0>`G)iS+*8K$HN(3dDtl3}NNZoIg zpBRwxuL~;;d8N)z{>~pa=U46Bh*fpc09EPPxo%RtuTbw8A9*J{nq+kV_1jXxmcoK4 zNWA!s@0b&E|G>WZ9RW%w2mK9B?hCIQM*eIh2i7SuZ8qVdno+K8E+JfS@!i^+L3(U} zjZj2OBY&gXaRDuBqRcOY{=?)BAP~2|W$;L&gZQA=ytl{DQH$9DQu*Y%e=s94=nOyE z?g7I1kT9mIk>$WdT$L8O{B;Jn2TP$CV_TIHg7w zNXf0Ct6^nLwIsZY-ei+3_&b=-rn6O>9;MP|8qUDi-y5;)M%6Y0Er^Qw6QwN$+Z9=& z=60toBg(bj@F(9uy8XgbkEN3ED=8$|zu_<@%5{{F9Q8zXvtljmQvynW9L-PWOQ@3) zl`#cTL^rrN@@tYbAF^B7%od8k+y}?1Jk5=b8FBG~Dj? x^$Y6%_2+78Y8TG^&$ILA|NAI4wewp4{aVoUyarkVX1q&FU07Du$ii6o{{RfFP80wD delta 8463 zcmZu#c|4Tu*Edv>w5TkZNJWaM?8Z`gGjIIwRAtIEngb5%J@1oKtKi{+FuH?70viaryUo|;%u6Ushcq*b%E**RFdwW&B;kymRGMBnVv z2W*&K7AO%WlhIuS1aIm|hpdkw2_K7q&Dl=5TK3byjcF!bP7`K}?DoYk4$6gYR+v>a z(zUcUJGE*HyXQ3DbR)-VwSv!E4>6U~^0V&`AQ1#j+=v~KLM_rF*x3hvL_^L>0PFn5 zj`>UKg!S2=RHDEA=2WNLc2zTaYmp@t_$C(s?rdx8$ss+O_vUHN%4yYxu0*XYbteF; z%a-=(7%@pai2{|h1xg>JTy~{ol3vDp*~S}&FEcY92fexDzgM>{cpYCmfaTO2D_bKP z+ir>9Hwq1Zv!<04lc{f360>6EcLDXTirc(Eu6AvM`BEr;<7TNWD;Sf|yxrNhHA^W1 zSS^-=8^pBj{WFtnHI3!_f*1oxIdyP4UEQ{JMi!NWcsa;Lh zL~orn;H)i18%Sw$ZLdz@;-pSNfP41HW9u@QXpdJ)Tvmu$DFjYDs$6@kyf^j^wW@FZ zac}#=&=ia0qAlq*Jw4grqV0mr3oh!f`(dyf(J8)(2qEBmXJd7!?V-4}*a>F_YY6AT ze9DMnhf0emjAz`aaE4Y1#=K-!Pa>HGct1HC!3ct+IW%9TX9uG$#+G(#n4X;(s- z76xf^u*$PTt6LwvVxdLPY046D!*^jy__`gBWP}e8nYmr$5VDWpS`s-3(~^3YrJc}I zdE2#C0Solj^+bp(3k>hdijiuQx%*QNET<>4eKiWLIqX_7{ ztMs|};H3uuYTf11`3nW*mrn}hnHG$^vG91)XaRYjHx9 zYJGnih_(J@x{2*iMQT7ECMa@a0a5Tla=K_i@IdU#XGaqG0g;!PH z1QgoOSl5NV)ATitD$6MvhZl#ch}srC(@~{X341J|%zGs@f?NT+;8a9AWeKGmJ;q9W zDM=C{r!I9Z-;GVbz-*M)adaQBTxoTXlj*}T)7J;(3(lUfw?_Kjo-Jq*dA)YK2*d8~ z#=YI$bA;8JP-ICtIBiJ755P!{ z!#$1WZVN?s*!^soe7KUSCr=i|<0d^iRPalGKVVQRW&2=(t{PFV!N`P{Q)Y%-*rf@+ z7Xvdiy6%3fk-+2LGW*aPgn_M2n0=^4x<)rqEzVy+kv7*a_VG&|=q9ut2(I=Q2=Q2H zO%m=Sigiu7>IS+FSd6XW)HY}IfU_NH>2j>6${FD5CHsrn;|3Tju9-c7uJ;5kj*Hdc zzI~f#jShm|HT0K9sPF{o?+`pLUjW8)6Vp5}+}C#|wc}qU4Fb@DXK>HY+4nk*SNVJ0 z=9^8C5Xd%Ka-GVeHo>9};|9JDSO`T``3s!BeN)EvYZ8W)?H^CPpPk{8co?IF z4-cA{B}Y)##(RI_agDk1DAx}EK`>}{VsCvx;*8hVMlVKWH)sN^wYmzosLyyB-tv_2 zkbbF}_9ik$%Z5JfdLqzu!2$*kijFUQV39X$sJqfi%HR!jwUyLBxtelyWCvw@X40rP z;&-Ql$V7QQ+tUvtU9RqgAM3bbNL6jkT5f&VQ5@vj!QBxVG%>c`4i?>`7tFg<4qlQ} zoQ0Fms+0@?B)g-6T%V5@Q~P21%BO=Srbfn=VkdvpXy9?zixg3=-9$W2#wRLuGDQbi z8eeX$x?on?<^F*+8eHVC(t!85h0rLIMa_SVqbOu^&!SC2CTJRPft6p(lo~(Nn&wX>xp+A(^)&B5eq&B7d3wS&L@rAT@tk^TG z5mYXaB<#_Y^eprx4ku&?lkPljTT1pSIHPuPPS-=Sy7Mxo)bB@Hime+Tpo+ikhs$ zzlvcQQE&JNHo@{LigA+ofm@=9hyrAyYoI=S+^;;#f8};WkZbRgj5BjHuz*Xm#AN7+HPk>VPnLfnL31)i zy3K4hMbSqc?JxU6vmX|fn1~_}H|JZ&TprBPb@JtF9u}SLl?vd(*Ywb*`c(WHr}$vL(xq5Vs9_tb9%?-6(LRWO=V}O{db*D z8nL*%u55hl8vLAK_IQ8IC5E97TRX@cb%eHZ!$-R}Rq^^$uQ{*omn$Y4568|=rO>g= zu@&<)K#vMfkJdnIkSEVj=v|W`Zf(wRE6L;)yP2FZi>ZSOOBg)vJMoQiQ1v>zZ>4j5 z>pn0FRI0Enld5kdB9cXI7X63uVpRA~gG9ux?Isvql0k&o8_=u$L!TB4Va-PGMK{() z7E8GsZJqiYA$vOzJw5dKSb5xmp1}UvxaKd}kF?!uZZC?{olT5ZJgO?QbuHY?Heb1i zrHNIxj3)B;>PgrddW7{XL2Qb=I=f1d^04+DkAQOUXl@RS*XH5xb1;CViqeWh9o4BC7dA_D`p((CT*;58|XtS zts_18e}1jctZtfrkd8U3VUo#JzJlcTuJT67NPnp9`sxr>t=AbVmCegt)H-lRN%RoV z?_`w49>|Prl0Tw{hHYF%OIF&}?9XhXD(HFr7`w!9TJ*S`)hY3;4{;_zob-Ci`-wTU zs5;q$!hBJEruU3rqJF{$?)I4<3d4p8mD8L)>moeQV3h{Ker?>IPeniR4u)RI(iaE` zRQuZ_ly78vn0T__T_?chwsAP1Y+GF2#O3)~GXO*oH+H>T)7t%VjYo6gNbvBsQ1Eby z4+{Ai*moKHZj0?sw9l-|I$6~WD5j$DPK8cvx~XkV&xE|3OsiuTz!|z)RXO!uguTO0qkYpaXsK-?%vGUm6LT!8Xj>D-xw>d zYyE!R@xri-)7Ty{wxpAu?Uob%ua~D_PKxKO(-9f@A2-cIx3 zDcl|D>1CVXkl}>WZQgZs9bIf`S+;afod<4qm*S2a{b4oWEbTz@7J2u}ep;tDzXEva zJd~0j`LA`K^M{4+J`Iw3BQ9BBXY^i6^+Y*rb-nfLez+OAyvwD@Onw@BELH8csFX+) z=-k~Mdl0W0AIX7M)STeI?FF$%v|ab}^2-VK^UFRF{77D0e;>jI5pNK^e$DT$fe_*$ zwr>U#yyH`u_+$QZw4@!Yv?M?eusmXarst<3EOze;k5eX-#_>gks73bU%`5Qib|v>9 zy?V4C+-Yw*!i22EmP`vYRid|)OlBSPk(|IZpcdkCRw{)Au8p|f%Mz^;a_$qy2yf)D z#hvutBl24K>+sjgr(!`{3rq;w7%zoku9WmhNAk-xonB4Wr$t7oAB76}zC#UTkAm6a z=1RB{Dv!MSJ6pRF*y89}C1Ieg(#fqSb8oLvwTd0jYG;j>6`MBLb5QnoeEao%O~{__ zDndLnS%#X-bWbwIF>$5v~1sInr-*GCD2EPm)mEc;2~(e$1j;`0UX;kao;Yo9eQ-Vvv@*tgKv8X zo7oOlgoopTIPbhisjXux9=Vt!_zMhwd9RzzbJmRc$kvn2etLevB;xI2pz7l-CGY+9 zVpU?-;*a)~_sP^VwZ9E7Uh3Q9T-K|+_A?3#M~Xx_Oy}Z@`(+-z;D=I{Ep9T>8@*Z5 zg%N-k3le$Nta>_FC^1&zae0A~ewHZwnt)(YhmurY3!&-PuGv z=?hJbWNbkk%tyWC4cYqU$oU;WE9Jc*)qt$t$NJ=;<#FZfaLv#^&MzKsnizVEJZ`{#P$#7=)U4Miar7eOXGZK5ocL>3{Trx9Ij z&PY!!T#*LqO66C!`1-x)VV;GgxakrWKJ!l3C(j1zO1Z~O$lRaGm6uFD<~|+orlfi! zc8DgTTJ^;I=gs0S+~LR@RdZh|d@@KzFYmsMo{!`0z>w^oQ7`W=#)5U*m^S5g!cjE^ z*a5Bh<(ehc1kY&SGHhe6@__Y%$u!UGdj1IXtTO&ZuuVhn#Z?vMCx`MufPsXIt^4Lg zDt$LqgU*5Z0%@n8m0yJ;R(^8~g~+*F5@F|t&iLxJSB-l6<}_I#pi4jrKKPfp*ww9K zV_j2nPH~edXZ3Mxab<^J&*+YS=UMa8{ScixYzn%1Uqd3OEoi z?6BFo*Qb~!*H(Yc4+Z`omU3%vbv~SMopQM1`rxQCxV|$gJU^Kl9#<;7BnPqnfJnZn zIHa;3)uU-PiTpAI+in!@`Qvy``#pG#!C{aZ%Ku@+n6XCq&H5iTBT;cYbl~1da%2Ku zpS0ZttI)=~%B_|4L}3*tJ-%A-a$)bYSVYp7=jE7*fb>7x))c6|N)|cZcS-KOpC3ym zUhn;lExm9C>TM|eGjT0`{!;I|tcg46PIq4%!M4lx&89wmPRP+xF> zzK!RU1fTcOVBKum-zUcDp=(6ci_46e5i_fnr|4|(~yjStEVK$yU&rXm^R0{%cFRjqM+dfr9 z!7Hr|%+{Qv|J=DWm85qN9;_yVKC2}PbP3X&TGmM<(k8s3==(0Q%IC2Xit6Jp;|}N# zVFPX{_IZge>HeD);<eRWVrmmq4D)N5GTN?_JVyuGW+$7{5>EeZ9{;*5 zy3Y>ld-{_nR~3b#)ZAB2NHI1%hhdUSMpGcwNbss0x`+Mbg;R;f$e~OmNKv z2T{ly5v)^RX#*@_Wp#%+i1~jHS08-cX`;$S^y46w8s9e*X?C{dQYjzgB9{I=bgWGz z;JnVmrb|8Hjt$05{9@6^`ol>Y&B6;DF(p^N9~qPUVSE5e5w*k`-`Cc7L#h&`^=xoJ zF?*o&6SNC0hv*{Z2WeIsb(d6icZnh2Ka$De_BPNzPUr$guzZ71qq<9?ku`r<+pqW- zE%6{tNI%!Pec^Dl`?)rbm@GWKuyd$Jb{^2kd zBxb?o0)%M7+4+uE;LbQ#+!*d?$$>mimGmZ=3EI~$ka2mo<+L+96(w7 zxx<`dw4N;vC}j^6^Ib^cZuXoGlxF3BsQ+uAzosTh7@H|tcj>iIev9*~2!pzkkTA#p z_QMkw+q4jXb)#45= zmBzT-?y01?_WqIR1XE5GGN(7}UhVN<1Jt*;8+DiBPnz)dlBjJC6w9!*nm>GV0E**X z+z)kZYG?7kHWSf`{g-iFNnn`c^y%ckL=-)6%^Hv)w7vzCAXxQpb|(F> zqCudzxQK68&y_{(ex=Ppl&uZxoO#xFYmqc&x%OcXb63p!U2-DTGi5WUPeGt`Zk4PM zSXz>q5*G+C3Zfk(E!1lvQ^GC8}Tj(VL)3my!7&`X41nq6YMIU!; zVDvgi41zlb4$AASSH1-~2utpmsm&ipRG!4$xvEH=Hg`Wu>zV&1n0D3nDYDKnVF&le z*kghN#*G5z)cZ$L zfQ}`s;UEUj4F-s4tor%P{08@@EFO5xDM3U<5BndaaN8)pwBl50 ziw zyb9DxKaJ{27~R*VpD?MEHBc)3j_*k9>4X!3Z}<@hfg$mA@Qet6heSqQ$@HgHju;)D zh$G1i1%*WsH0j*3C11raX{YoNB-F2sH{hEg3HQPr|HpZli+IRO9V`<33P~xM@qo|v zAg#wY%+ce7rtB$D5rXO8uR?Bq)gUN$#QS^r#XpIeT{^0 zRK9wNp+2`mP9+X?CAVT^2wz`?If7*W{@;=&86yVL56wsRhMS%1@IEyNu@4gCX7_Bv zS@6|s-)V5?d!RRW)dHlxfoVw0TIfG`6Qfyeb3h#rSlE_^Tt6%~p&91L;2<`uZA__J zeKTsS-{b8Wct9Y=G4KE=v!P4Jrb5tz>fA6241Yr932uBRN~mJr6UPfHw$rC4(-t-+zL(cHpbk)jYu+FP!! z8K1B9%7@BGaUi=yLWqsNPOoYW>q@{H<4^oLzf*_~%HM=RJnk<2g79xUCRZx?Yo-$T zEcsNFvPgSCf;fD=kM)WXmIS2<4HL&yexnqI z?yd-=a=)zp&RC2XMBRhGzrh)E0T}9H`~7kV^2!P~tw)L0lgxq4sl}pC>d8Vh(8hZA zKHYiwUb9~Tu!_1sBW?|PLgEy#k&TqmHp1m`GR$a`$4bzgt9O)cHci2YQY~GAE z7&-+6BPGnSi;$~|*!}}y#}k_OkBIAoRr2sNcN_c0*`FHb_+l*sxrnJ1>ud&d?dhx< zwhOw%-+!Nu9G|Xt8UY?h?}F07(X#PyY2n9%-_wW2l`9@)GYfY2iEDr678t>Jz(FYH ziG%501}MAOWuJEKpm>k^i-qig!JF|vkEmGeSJ@NW&j=8LuS1}eEeQV{s0m(!Zyq>_ z%`5|&{W|v`WXB=kLsm;5)wMb*Z2Jw=u1wv