From a0133428adc5ec46fbe2cbe950f8e4ded8f00d4e Mon Sep 17 00:00:00 2001 From: Laura DeCicco Date: Tue, 21 Jun 2016 16:25:59 -0500 Subject: [PATCH 1/4] Updates --- inst/doc/EGRET.R | 54 +++++++++++++++--------------- inst/doc/EGRET.pdf | Bin 777347 -> 777428 bytes inst/doc/rResid.R | 37 ++++++++++++++++----- inst/doc/rResid.Rmd | 77 ++++++++++++++++++++++++++++--------------- inst/doc/rResid.html | 57 +++++++++++++++----------------- vignettes/rResid.Rmd | 77 ++++++++++++++++++++++++++++--------------- 6 files changed, 181 insertions(+), 121 deletions(-) diff --git a/inst/doc/EGRET.R b/inst/doc/EGRET.R index 9df5a92b..c0df8c81 100644 --- a/inst/doc/EGRET.R +++ b/inst/doc/EGRET.R @@ -316,7 +316,7 @@ eList <- setPA(eList) plotFlowSingle(eList, istat=5,qUnit="thousandCfs") plotSDLogQ(eList) -## ----Merced, echo=TRUE,eval=FALSE------------------------- +## ----Merced, echo=TRUE,eval=FALSE------------------------------------------------- # # Merced River at Happy Isles Bridge, CA: # siteNumber<-"11264500" # Daily <-readNWISDaily(siteNumber,"00060",startDate="",endDate="") @@ -324,7 +324,7 @@ plotSDLogQ(eList) # INFO$shortName <- "Merced River at Happy Isles Bridge, CA" # eListMerced <- as.egret(INFO, Daily, NA, NA) -## ----Merceddata, echo=FALSE,eval=TRUE--------------------- +## ----Merceddata, echo=FALSE,eval=TRUE--------------------------------------------- filePath <- system.file("extdata", package="EGRET") fileName <- "eListMerced.RData" @@ -345,7 +345,7 @@ plotFour(eListMerced, qUnit=3) ## ----plotFourStats,echo=TRUE, fig.cap="\\texttt{plotFourStats(eListMerced, qUnit=3)}",fig.show='asis',out.width='1\\linewidth',out.height='1\\linewidth',fig.pos="h",cache=TRUE---- plotFourStats(eListMerced, qUnit=3) -## ----MississippiData, echo=TRUE,eval=FALSE---------------- +## ----MississippiData, echo=TRUE,eval=FALSE---------------------------------------- # #Mississippi River at Keokuk Iowa: # siteNumber<-"05474500" # Daily <-readNWISDaily(siteNumber,"00060",startDate="",endDate="") @@ -353,7 +353,7 @@ plotFourStats(eListMerced, qUnit=3) # INFO$shortName <- "Mississippi River at Keokuk Iowa" # eListMiss <- as.egret(INFO, Daily, NA, NA) -## ----MissDataRetrieval, echo=FALSE, eval=TRUE------------- +## ----MissDataRetrieval, echo=FALSE, eval=TRUE------------------------------------- filePath <- system.file("extdata", package="EGRET") fileName <- "eListMiss.RData" @@ -365,13 +365,13 @@ load(paste(filePath,fileName,sep="/")) plotQTimeDaily(eListMiss, qUnit=3,qLower=300) -## ----printSeries, eval=FALSE,echo=TRUE-------------------- +## ----printSeries, eval=FALSE,echo=TRUE-------------------------------------------- # seriesResult <- printSeries(eListMiss, istat=3, qUnit=3) -## ----tfc, eval=TRUE,echo=TRUE----------------------------- +## ----tfc, eval=TRUE,echo=TRUE----------------------------------------------------- tableFlowChange(eListMiss, istat=3, qUnit=3,yearPoints=c(1890,1950,2010)) -## ----wrtds1,eval=FALSE,echo=TRUE-------------------------- +## ----wrtds1,eval=FALSE,echo=TRUE-------------------------------------------------- # #Choptank River at Greensboro, MD: # siteNumber <- "01491000" # startDate <- "1979-10-01" @@ -384,7 +384,7 @@ tableFlowChange(eListMiss, istat=3, qUnit=3,yearPoints=c(1890,1950,2010)) # Sample <- readNWISSample(siteNumber,param,startDate,endDate) # eList <- mergeReport(INFO, Daily, Sample) -## ----wrtds2,eval=TRUE,echo=FALSE-------------------------- +## ----wrtds2,eval=TRUE,echo=FALSE-------------------------------------------------- siteNumber <- "01491000" #Choptank River at Greensboro, MD startDate <- "1979-10-01" endDate <- "2011-09-30" @@ -408,20 +408,20 @@ plotFluxQ(eList, fluxUnit=4) ## ----multiPlotDataOverview, echo=TRUE, fig.cap="\\texttt{multiPlotDataOverview(eList, qUnit=1)}",fig.show='asis',out.width='1\\linewidth',out.height='1\\linewidth',fig.pos="h"---- multiPlotDataOverview(eList, qUnit=1) -## ----flowDuration, eval=TRUE, echo=TRUE------------------- +## ----flowDuration, eval=TRUE, echo=TRUE------------------------------------------- flowDuration(eList, qUnit=1) flowDuration(eList, qUnit=1, centerDate="09-30", span=30) -## ----wrtds3, eval=FALSE, echo=TRUE------------------------ +## ----wrtds3, eval=FALSE, echo=TRUE------------------------------------------------ # eList <- modelEstimation(eList) -## ----wrtds5, eval=FALSE, echo=TRUE------------------------ +## ----wrtds5, eval=FALSE, echo=TRUE------------------------------------------------ # #An example directory name # savePath <- "C:/Users/egretUser/WRTDS_Output/" # saveResults(savePath, eList) -## ----wrtds8, eval=FALSE, echo=TRUE------------------------ +## ----wrtds8, eval=FALSE, echo=TRUE------------------------------------------------ # loadPath <- "C:/Users/egretUser/WRTDS_Output/" # staAbbrev <- "Chop" # constitAbbrev <- "NO3" @@ -429,7 +429,7 @@ flowDuration(eList, qUnit=1, centerDate="09-30", span=30) # constitAbbrev,".RData") # load(pathToFile) -## ----getChopData1,echo=FALSE,eval=TRUE-------------------- +## ----getChopData1,echo=FALSE,eval=TRUE-------------------------------------------- # Sample <- getSample(eList) # Daily <- getDaily(eList) # INFO <- getInfo(eList) @@ -497,14 +497,14 @@ plotContours(eList, yearStart=2008,yearEnd=2010,qBottom=20,qTop=1000, plotDiffContours(eList, year0=2000,year1=2010, qBottom=20,qTop=1000,maxDiff=0.6,qUnit=1) -## ----tableResults1, echo=TRUE, eval=FALSE----------------- +## ----tableResults1, echo=TRUE, eval=FALSE----------------------------------------- # tableResults(eList) # returnDF <- tableResults(eList) -## ----tableResults2, echo=FALSE, eval=TRUE,results='hide'---- +## ----tableResults2, echo=FALSE, eval=TRUE,results='hide'-------------------------- returnDF <- tableResults(eList) -## ----tableResultshead, echo=FALSE, results='asis'--------- +## ----tableResultshead, echo=FALSE, results='asis'--------------------------------- print(xtable(head(returnDF), label="table:tableChangeHead", caption="Table created from \\texttt{head(returnDF)}", @@ -517,13 +517,13 @@ print(xtable(head(returnDF), sanitize.rownames.function = addSpace ) -## ----tableChange1, eval=TRUE, echo=TRUE------------------- +## ----tableChange1, eval=TRUE, echo=TRUE------------------------------------------- tableChange(eList, yearPoints=c(2000,2005,2010)) -## ----tableChangeSingleR, eval=TRUE, echo=TRUE,results='hide'---- +## ----tableChangeSingleR, eval=TRUE, echo=TRUE,results='hide'---------------------- returnDF <- tableChangeSingle(eList, yearPoints=c(2000,2005,2010)) -## ----tableResultsShow, echo=FALSE, results='asis'--------- +## ----tableResultsShow, echo=FALSE, results='asis'--------------------------------- print(xtable(returnDF, label="tableChangeSingle", caption="Table created from \\texttt{tableChangeSingle} function", @@ -608,27 +608,27 @@ plotConcHist(eList, tinyPlot=TRUE,printTitle=FALSE,concMax=3, showYLabels=FALSE, showYAxis=FALSE, customPar=TRUE) mtext("Custom multi-pane graph using customPar", outer=TRUE, font=2) -## ----helpFunc,eval = FALSE-------------------------------- +## ----helpFunc,eval = FALSE-------------------------------------------------------- # ?plotConcQ -## ----rawFunc,eval = FALSE--------------------------------- +## ----rawFunc,eval = FALSE--------------------------------------------------------- # plotConcQ -## ----installFromCran,eval = FALSE------------------------- +## ----installFromCran,eval = FALSE------------------------------------------------- # install.packages("EGRET") -## ----openLibraryTest, eval=FALSE-------------------------- +## ----openLibraryTest, eval=FALSE-------------------------------------------------- # library(EGRET) -## ----label=getSiteApp, echo=TRUE,eval=TRUE---------------- +## ----label=getSiteApp, echo=TRUE,eval=TRUE---------------------------------------- tableData <- tableResults(eList) -## ----label=saveData, echo=TRUE, eval=FALSE---------------- +## ----label=saveData, echo=TRUE, eval=FALSE---------------------------------------- # write.table(tableData, file="tableData.tsv",sep="\t", # row.names = FALSE,quote=FALSE) -## ----label=savePlots, echo=TRUE, eval=FALSE--------------- +## ----label=savePlots, echo=TRUE, eval=FALSE--------------------------------------- # jpeg("plotFlowSingle.jpg") # plotFlowSingle(eList, 1) # dev.off() @@ -654,7 +654,7 @@ tableData <- tableResults(eList) # dev.off() # -## ----label=savePlots2, echo=TRUE, eval=FALSE-------------- +## ----label=savePlots2, echo=TRUE, eval=FALSE-------------------------------------- # postscript("fluxBiasMulti.ps", height=10,width=8) # fluxBiasMulti(eList) # dev.off() diff --git a/inst/doc/EGRET.pdf b/inst/doc/EGRET.pdf index aa79d57d752ee0990ba2b060e8361556a9610145..6a68278d01f88bf059fc1856b22b2b4718f8a86c 100644 GIT binary patch delta 38381 zcma&OWl$Yk`aO)hOK>N+JHcIoI|=UYZjHOU6P)1g?k<7g?(P!&&Al^oXXaP+{^x!; zb!u0yUVA;V_vv2Ubpl!#=UN!63y>1F3~3T3kq7~Q>BW5&Fe1P0-Hnu)Nb=9boYu2k z_Bl~O{MG9Z`hpL=9Rb){i2a+pVbw7L(K%k8<-L9~?MzX-dS2Y2%JPT?0+;d4-uwhF z$C=3!y{f(6i^imsRSJ1j?6=vHyxS)ydJlK4*f=6Mmqe>BCHAGav|9S^ZO`90-oC$9 zn#BShVKQpEAE87Fc~PJCJ$+vE=K0Hfrm)VLP~;4C+|JX6Ob7C^wKDd~Gt+}|>h^>I zS(EHHJpoT6RX%!(1Zm*mt?9J@li~9-GY(+*Au%>CY1sg_gZK`;j6JlST*>6DSaoUF zSoxt@Nr!5sD!hYSW+KG+Gcl4s=CDs_{sSv#-RR^Wo@F_RLc2 z;B!atW^^`l@m);Z7UWd3W3RdSQu6*j!>t#c{<0n2m{QK6DU*B^YS98+uj~7bn5}O` z3Alb{4|=wlA0^>t*ax@mt=ay)BN_&xdm#Hwt*G;eH)CV``n3T(&q18`9|+KYfMKTq zm59WGPO^cmFXRpC5py|7BjhH{%&XJ?oUb38w`(;5MXD5^+VddYTx8wad3xYP9hy`J zwbps*>t-1ApR#_3KxD{vdNDWp5YlHqnk@0pL;!z;ckoW4gw*8cPz{%uvu~_T@4y54 z$Y?{CY1KG>aaAXY4SlM9`FrlIM;JijC!efdr#Tpp zxBW{k&kYV90lk~xpZtGiue_7x-A*^HYeC)Ri%#Xa$%fQ@!Bnu(kT39O=6Jh*XRg1G zvIT4wXmXy~IP>Mf$as#|5B)=Y_8Q;*Wla{uac`f}LVV`Y*}RQ@eYco6;hcl!Tmp5L5awlJ>+V{p_IdEBbvMHf z2Nj_FL3S+-YSJ{8dH8py{HgSJcm9U{?oPjjuar%Sx)Q=>v&%c0d!y4W>n9li( z;*Y(L)8AFYTgoT`kdT2tTjF;Nf5a&|1|&2OtWbY>V2}>>2v*aU_+_ESzS5YG!r6LM z;F$EoG+oUFM#dFhhIPspaYSD5-(HXB+)n`BqJv37q!*TtD83BI}g&nEh# z<9;ul7gyNu)?#Sd+5wA$u}9voDa94>K;rQaUY%y4j4n10Qij*tJC+R0sMm4Dc_Wjf z3YY=w87=j(!S?EM&ILJ!$y)^No&1}7k=ogVBo>&(?|(*pqdx{Vw^i#KS_|^KjM;xS z#qSA!#{(qEgue3nw}S9~`**E=uaO)9k*mF}R^P>WNG;&-^Mkv+^V;t0oh)YMDhHss z>)au&CQ~D5PB=8cyO6N%Wq1M}?%|MSeLfVMbPr@6Xm;p`sUxj^CHWSgXdUv$VzSru zxmT4Zld9!MsX~rB$+zpsbT*w(Qyg{a2Y=@HqvU_a|GPL{wOYdh8sQlV;SqnfmK*S| zrF&Bq4_lKX;WenFB@XlFDSuYv&w6I$LOjUKBrU^L&pcBw7H6CgJROHrbb2?W=an_G z+7j%gw5rgZlp!aW_X>9ZUK9T)r${u;6J0%{s3@CeO)4DWB8X#zyFu_A@6;=4`D^?UI$Fy{VG}2d*o_ zPBs01>>OK<08xKsN5iAiGHF5aPp{7GnI|_*Ij4te_xlJucnPeJBh%FLvV&{AFAB^a zi+EfAnE8*MB{aX1acFiEuZ%E_HcjeGh^)wVSDTYHmO}y-@4IYH-sHzJjAuI>5l;WL zGyb_-z)m4hc0_%5(4M!eNlz)90mKkLiFG{Yd`ExXUl3d~x-rzqZ1&;{aC9E3?vybFdP`RMk@^zG8~+e|$Yv zOKCEJ_v{-P6JzK6ykydQ{IieEvi43wi#pRP+Z7H{bl6w_@{bM3ty~seaoRI(s!2fX zAnonUy~)+Hy6dO$xcK0I^uZsSeln0 zzgcx?go6{yicI2B-hGV9RmP z2Ks;TvNHXbh~LBi77_IXI{ug9{{Nq@z9O2Q!+XzmH(Y zM=~AWveOB>Y|OrV^B`G$U~Iywsx;aUA{`bZ^k~9koB%&6q(O|=@@NWrXnVIoh1;O( zT?91fR5Kq#ipqJ$?38x^&bziOeNn~bJGrI^rkf`4Qd?(r7UTV@=X<2IG`$*`#B0Wh zuOvoiW)l3$BpF4dU$Ir;T-KHSSC@iBBE$CujEe8%=jIj#+czL6wEyVw(&EV;PX!Nwd~}t z1oe?y`Nt@buX$g0ZFZaSI(vPPP3sx4KI6das;^ZZHSrL}MWbR@!fmljQtItdz`J%p zA|KJntAW4QIx}s25VjSggR2voo%R0N!-QZrg+;1hr*W}6Yt!L)XvuU-Ebo9X&oBCl z{W>rwu%^@A1tX*jd^2T~1;deMRH7L+I<9!Y3XM_1+%?}Hk%2PDaId)H!0Ph3s_Z=N zpOQk11`Vm{-Qd+<=DuF$PS)Uep(6z~Irpn9c=XyS6L`7tUdPa8%14h698Jli#TtD}5bWoE84Qq` zjH4s=I!E6O0BR!-a_^HTZHR**!t7?Q;3D0GIC;`pqK0QkE_W)qomzc{@=vZO>)n_l zIam-yvnL)gD7~b|Hj)|-MNsj@ZdC_R9w(@kA9@+Mu?>mEYYAnhYGy0?UKguI*#p!E zoe&5jS6-q3O-}q|-ZTh(3VD&;&i-%FfmcKzC?J6bzPq6j5?Y_OuTfKCnmgVE^OEsr zJvUSQi5j`Muf@S6f ztt|}H2H$d}d|C_hxrNZ(V~|VOeD8pV&5cGCN|d++GMM+M!_&G@xVdt%-m%nK=0G-X zPZYL2qjNjuqCvMht-@!$0&?zn_nMhQCQsExJi1-k7&CPjY%rH|4JG#0wdRYg?-4rA zn>w6FLg2KWI<2vRVm!1=CTzznH;AiKBpK`iQ3+!*xJsKV< z-8|mmt&!O{<@MyRU~ja+i0_RGHMXN?mHCZARIiqG)=>?5HOZ)EvA;U#AMCnvFKfg| z3|Szoa24Qb?wi<5Q$OI*L9Ak1LqnxDzqoV&v3fQsoB=(8*?PPbeZ__e3kM;JT(arxhfOd#@4$ zVBoAhGQ&e`yD2=<%5gnuDRSiP!Y%Fbq z6kef^DV7}%NcRR>0NZM%JfyrnjKRtkhrQN7Rb|_h&GUO+sc?r536Ec)c%Q=SRsz0j z&N4ro1wb-Ny-xJMw-B0qrd+{o6`rmGMhzA-Rt!fx10N^xgW2b7cpPhIQjx6df>U~$ zX44$DKzfa=1aiox*AO}($-gyQ1J{#<_eb2r+W*rstoM^2iE z`z;nv=lzppIX~-t)nT1HR9grx5b?|GPV52_b_Tp!+kR};dLAz-9n`qIiuAGqQ>D-u z*f4`Gd7f@C{zHdllTe1;sfn_--KBcac8AbxCJ?I)@$Ae+7PvM(SiJLV%>~9K(^)2F zW{hRi)F5`+v#~z0(@PmUqH2)3JQny+uwVM^7H@dKvI3Hl@4BHy8*Jdf*KkB?UZEO- zWEY)bkzv=2U`n}tUG!jRjG}4*1nNhbU0LNWr(Ls_U&gxrqPNv?in?votso#pgo3Tl zf(InR?+p#F*;O*126o`!W?;Aova2)L{Jta7qccUEpJr^1#gYc}3t!C06b|=DsA2}8 zcI*0+@YBVo=82{L;`9JYhPQ&g*D-&4dL&Y})w$k&zfJ*zRn!#UQh50R-~a|#2cwtn z5U?Z=683E6T&?yPFw*rB?t$qRsKMl*Ng`*cklP7f^Tabh8_k7SpNTsC_*(+AilsO% z5<(O~zP#?Y?dE(9P)PBbSM6co)-a46Vi3uz9GwlXc(Kng59XPV1itV~4JvzT+6Qk+ z5KgO~81CWRu^)JsQ(}bxQ#pE5-eWA0H5hOty*(%4T82*!2vc#J7$0$k(MgdQa(a7O z1Q8FPv%ceiL96my+zexevkf~uid8ZNV(#SdQHO|~ksGCAn8q*G9!s*1v) zC)3^9G8Y@9P+b#Gm@wFH8>{vU$xv^oA%4UX!B+@4fyQmy;g%Q$&NUS8MId4FKM$dn z4zY)*VC!vb3%M0(VD%)gDGiATjjst-2e4OySRgLKV;Ll&0Q&#_c|E;A@;!e^Q} zGP2(Oq26y@W&+6!JoM~c@MznveVNRNeIiF0!U&zpo)j~cL2Gu4p!!6HfqX@1hkqdi z%Mxy?7tO2rfiboOs24M zrv74P`zY4EnJ#b(cDG&_WB#eEFE$VXztd|Cp^lkloUoF9rcaR}O^804?&%!8KeDx0 zD}($LH=@&?#<^-uV3e1^lqza)TWlp{sY|Fkf$M9IRuX#$i4Uo?Jws)0W@16Xw{WE3 zl_F|v(4@3JU@IBLQ2qou)jlZ#KY~oWYW#~bYz(*!O>;aP>fzE;^h;#rDn2Ks!UuPd z{bdb;^f}WUj}?s9^oWu`USw++45oSeF><+{1H7VU31`kT9?xk@ z38_}U??*lGLsj}{6McDBBt>jz;D&3|dj4ZV1x$^rjPpj4lFY`B@NX7W!lAKa4gK7o z3kD3zScVY8c5<_PN2Wq1jkj z8Ms(DQ#3z;CxbJwvZmxxgV6!I&1ZaFy^Dkj$kjUuVHd{%w678fapQq+&|sXFN^O7cXbyoM*ZVrU$kGbNZbr zk2sF=K?z1Wo8u%R@N)PeSQo<7u`)AQ+jxYn|p`X;vsT#U`zAT3%`PawH| z)!fM0(>>|R=|X%Z@$9TK${os&bh9iR-tca2?<|KGJ2W&e7>AEXQ8Pa^(6v2W1$|g+ zgS^q6nbF%ASq;bntrrKViLLTocqZ_X*e({2`bjT^nC|GKJR>TPzTX=zghqiTtu?uo zvTneWopa?EbvW;L(>7}}MzG=|<3@SI`7qsQ+$~L+B0zvCIma);Yi#zzKD=<*akQ}p zN=O|VOC5`j@9Sunz%iDyG7mGMPG0lC9dOia*pRS~I!hD`fSi2#6VRk__ z#-%M0z!yDw=*NL0+QZn6~MyCP9&x!tunF!GXns2S>L)B5ULd(OUL? zD)Igt8wZH%PLVzQmav+Dx{Z|l_PL3 zOiLAW5Bk{sskHo}yo;982{V}-FQn-~Q_*C>*$~rxxN9WZLGmo z*SNK)^r3u|C{TWja2B{r$uvWFLMHm^pU}r|QM#=GM1nA{*B@i8`PFcO>&Fd;L23jW z9cm`dn$xF=n#qXUOuQN>fmMyAejLJ*pB5KvJZNHYf?A8yi_(wV9&$Vg^V#k^7}Qc>Ztb?Uk!`C=o9|(sOjlz*vI0& z-#_vj`E`3kxR9aB$P}6=gHcoFxUtT4PDZ#7(TD0H?iXQILfUa^t2&5R_`s}?Jh$E}w^!={4%L`X`aGUi>*Qb^$KyXYTgK!P|VX8op-NH9Wq@(LUruHNHa=Dk`En4KzmL+Vm(0rQxr{^H6)sP zhTPXb&FE`*#wh8GygzZYyagky5p%0Q0NN2%!r@7Wc1YtrzG`W1_sQJyiAirfg2K}f zGpEQQVx0I@EC}%;BA9#xZVfz!5!@L5lXkNymeP9_=&}wYcZCHp1x7LF2unfo9 z=U5WoEGCY%?(@n|(>+9a@L|G`BG{l(U`C2($Wxma$g!*KY7I3-7E%k=EXBpdlWMt~ z@K6%fQmE`1E&bzsggBWXFPDh#@aM+TRNjM+Z)R5yN?bsnwC2OJMo4)H6Mp&wj{yK#YfU0ncp1h(ZSAF z?^U%4%(R76t$f!S^On#Qa1>(r$2Sql2hWf-GKW=P4n8b`BDx(mI$LKN3wh`u>8Ib z;bQ%L#qj&<_jLwS0`vC|3B(ZaDS~3)X6W2Z3|t(XZI_Px+s6AJq@iWtTO|LR^bZdfw`pMle3$a?f~oJpdH#mn ze=z^Rfx!OT%o;)vY?yyn#*-QZ!e90->;K4xlMTe%hrhv`n+Ac`v=szl^_Q*u|Igg? z8(fY5Z-3T~4dySGo0a?DHU165ZR%=-Ap6TvW=;8jVim?92C@GomhHa{=HJ7=4CX(H z19u_hQU5Mj=pMo<;9rW_{>$ZNOZk_}-300ZxvKd;QOt9Y^^AW_{M$1B4=_W34D>&rair+sf&JCTY;6B0JiJM0530Wh1$%%t;rn~gpKLH!0)K_}-@ZND ze?L12VReT2{ef*_b%wjb2j^hoP65q@#e`sD<>XFK?56~rosBiLzBbZBxVqA=uB~bN zLpFiIQ5C?sqGfJogY>XBw|T;#Wq#|)+9cexm+LzBX?pIoeP6VUVqQqC9em^9YbZi$ zt7o4>D*y>3j<;6mP;nI`6L*M+r}|G4}J-E4X$`}FOanIU1j!@h4Vq5u8e2w&eJwpA zRKc4mh~i3sY-%dOH2ei`7wXQg0jwnq#7c+u`Jx;|v!pJDc5!}^p!-(tYff5H(%A+~ zQ*|{{_m6`d1VN--sKy`=Y#pp@Vu_L)-87w#7(ysO?1|;fJTPFf46;|za}G4X>(x5C zx_;Q&b!r6l*I`ED#6yA^&PsS}P;d`zpqo2l1m&FTGNFKWfz<5a*z0@Dadh+*>>&1e zhg|U`nFpGT&bG!Mk-%K8K&Z!FH#XH_?lw~ChY7l;2KqjZ4}gQTfC#Ruxnu%`b9eSo zud;x{J0v*CJgCp4f~ha!^b8^fNF!)Rr^4&1Z#J=STNyM}_r?3iw*#c43_Lx=Mh9pO ze=FUpt&QUuKw2OFp;{!bd1BjbXCyxuGw|Bx7xhN7jK!_zsmLP~6j+i%& zH`)CA>8q|H3(K^+l#7iJKL}3@%6plz4ki%v=>D#(hz|4Cg=zo#T+PxC1&8t`u~mis zY{D?sD@b-<4;*JiV;Qxk-TQ{VxxYjx41g@zfPIQ!`hC=-X? z0R@R68bTF4C|`Y}pUd8TE9c+A-tsDVwX}7IZXO^zvR@C@p~U8xVl5Z z#!irei1u)gm~M8&FE9ct_*tAf{x&sVM@9tV6(P<1IPumk%$x$m0X^^f99pP9JI=@B$E#Ho3M|_*4tgq;RZHO;Ln^3f zQtU6{9TCply#tLu5X5nD#~EM}B&7b5KyHjAi>-81 zjtA$lWh!7lMt*pNHgk5UF6OE~n%31EYHuRe7G&;>n$po|RXDBh0YiStEu$f!g-ttV z$}qSHnf$&KgaZziVe&ZPy8(7IpCb*bU-CtN2v>PMB7WN+%BV9}7!km4nIr3yD9&u|U zx%NIotEY0P1(SWu?{qe@B3X5)jNk5AHMw0gESOv#8bogW=vdb1+r63^V|Vx%iOUr9 zbEMzPb4nG>n+O7b><&;q+z~+f76cI=mQyV%j{JGv1f@jBb9?^sONrK_PlF8l0`{!$ ztYNvV$J5U}Mc{-Sb!)saoc|j-L>!HNO@~nBiCA7WExKKzQ8S4e9=mOKy13Ts6*V1! zKuyvx>9iU?tHl<~1D+`+g$aIqOq?c)i*q2`cW)b`zJ(FW_iO;yY0hAg??`HQctSF) zj4#{XJcbZQ1Qq9N&#}8TZzD%~K75|L*ZKbK%m*nFC-Qu9#%K-$hBrF@r}<@eE}{)6 zM%}I7FExEBj;@&4_ME|ALXFCK5PryESF`xV zwxZ*j`=F<%q$>Ras^F*mHKAl<&s|Q+`+34_XqIBCOxBt{jwIVNmWQe;S!xERSP_y> zn2^`>F(D-ZZG^4AVn|Kovwza1T+Jj%6eh}(HL8T^m@Iu0zaUD)ag#L;yBmG@(_w zHE^0s%`tIu@wj<)Rhma{5t@k1jN1rto;Za5^rSQoFuZ>GD(tplA%VkJfIb?7i-4^x&69^ZI1ftt;@}?4fLSnEpvyF@QvBOJQ zZ=lvye;ov{qt;zWCzcyy82Q#w+!Wr_aN`{B0UZomn2e(gPzLe-#qlNnkU&N!D-)ud+G29XSLjpQ4HW z07$inyfC&Y+$!ZMU#d$F6m@xV^k|v7j}d}~Bo>@AaQu+*PGdB4DH{^EHt;_Z>3>rB zyxfA8C}Vm)0FT5d4Mu;gHy?Y50fyaBY6jshVwMy$p1!V?z)^7*eM8}DN{Kf0Hq1Ol z?Fd;vUWLjZ@BNDDXfval#zEH<0N)1!Gw#Rd0bx4DPL7SoCD zsXtwqCwo+f@ew;r993FQm|0b430)r^@Wm`+~fWH zyf)vXPd^UdlbavX-A6B0meZ`}@&d0yT;f+y7TlIlgh!~9jHh2<;SF{%Xd|edIA}0d zTbHge`C#rm!j407Joxf|ib&h?mxuw|t4JeX#XeaN&rKom?&|^pr>W4zub8Dg^|t0r zbqtf--10|#OFj`(vIR2`z!)k2j1_7NOQQLrE5C7*;UL;r&2h3*w}v_@%mo!h`M%QE z&BbFaI>pzRLhMPB9gvQ<#Z%n)D3aPH39ZxJr#^&f;2SV zkJuV%>V$4b1gddeXb*0LpRfObJq z#`rg&4gH`{+GDPw^rqO9m~i$$P7HJsIvZ20z&H5e)n{)2Ph`#Ci|torCa@P3^tdt; z7ug~RzQAR>k=7&lOAug{gRjc`x?<6oQX3y35l{TY#rFO(eZ*~J_hi8H)1l5@_ilP< zJZj47sU?fIK}jIxD9Q;UBA+X~UoPCXbGo~Oq#lXrZV`Ktkm1zsV~1wetqqZ%a7TH;e@kzP8rz7=Ii3zDQIkACwt}p ze7)KBC^P)@^ogi?=!{9O&sYB32==1mYr_Ol6Rt2`C8ku^nv3mORFhl-PuhJAD+&B? zTN!MaTitb+BKHS*vi*KKSy7f`#d@H=4bnRn`C~E(R?&48=aWyh*F#Xrvt2vwxIInv zhX=#t^0BvaNml^(BPaun1OxJlklts9$c&tfV;WjZqB11#DD_k?0F`o^&69Kg(5A(T zwh}t7)Fb~WJ4s2BX9bfVuMw9tA-TNY)*cq1aS!F?#=qfzvnc=+_1zenZ36qUTx5e=t zh{#<8xf~Ob(PH+-m@Y{eV}Ce409r;6J&c-?@2JX*D{GgnIK5}f z1J5F-UV+qxJ++ylGaeb%iw!;5m18o4p5(%W{Dt2}pEMOjZzlH`1;Zy^#d9AevSqG4*-H%hURW;_pN?Mg+&3=iX8KS$TJu_Ts_H zG?cf&fMmEW zb+9Hj3+b#_)O0OvRLi!oiL3kjp=Eqk0(KIA0)~j% zx#$l+sgT%nptOv9kRLf0O)G-(GePd>Q8FtYfSm9|6(ZlayFC*AiP@NH)S}RgpvT>6 z0P68;1e3i!GmJ%H9*aJB*4`vLqX-hUyL~l8@Q60stvsk$f?5{me&$9)wWJP+A@vlGYvKr5&m8IA`7^YtOBz6-o)a8w=9J}jYB_8ck^8RgE> zZHG8Kic!hXK`*xbm1@D%u1Tw9S0=j(-x|6uW><0jiF^>X$oT4R>7!T459bfFNfza6 z$;qD;NFXeDa*s=f{MsH`D))azShxII1eB`q`{NAnDio(wZe9;Wm?s4xhVRNHv8~24 zu9I2?7N70-#4q^L+%~*j53!i`HZwR&i;1~1U!pu0ZANr=*85O6$dZf#QMP83Ky|Qt zJxX(IW#i|W+b^e*JXr7X8POM(jslyW+iY60R*ZQdJJKLJ^urUf7O zR-b-ljC*91?_Pc%Us?3k5z+quTGDP|h4op!%Pb_Fw%xh}6fZ2~UV+F;mn1lS2vh6lIst#L#k`ju!hAQ@ zDcj*M8y9a9dV?Lt`pWQiZWBD^gJR7l)l3i7S(8dh>#R~Oq&}nX3>^IAM8G#^uZo|A zEJ#}6g`lGT4^U+Mvi(7wL>1jB@;|m^T)#QD9$u6Uufm*0QKOF5-dBA0tdukQDgU!v zaBCHb$7;H9*LJtV5hwHRL{YsQLi`Jew^;N;uSE!K+N=JYm6|(|XgVAd^!-ml-1ttp z0+emz8zQgxHpk2;;?cg1;Zoq_2g4-S1Fa5tTkDdk$8cQB3|hoXq$92^&RgvdV@q5I zdKUNSs7NX`+n*&0h43GeqcRb2-f^} zvN2dfFOJhnrUY9Rvx%AXKtgw5bZE;NE|-=&P_A&Ht^CrC?q^EZ2wDNK6@(*D<_pPO zBZBcI2?lQ8HAGFD#|h5`txf1#?L^wI|9g;ZfF)`?mZ_XxsFb4%Zzp zx}MC#`~uQgv2Ib&BrPb$8Uwpw^(%(;Y)0KusoA-(zW%AC0ZjGDzde~-@f!_6)qX}j z=*e2-$W79sj+%=1bT2Pp*@bgYIZANHj;O+CdYpO-OW{r0^>aDyn=9={@rmSxsPG}a z(D|;-)YEJwB>PrdQ!yi|xhlr`vFRkr3!G8-4V&Po(TbB!SKBvxnM;U(D{|t)(&UYs zYiJfhVfs^1Cj=*2orW8VS=C%SbFWnu5w}k#lKYa&n7g0B58yI^*wG<3S?%5UR|^GbEL6aGJ}VrAY*WU0iWb9PO`JJU_NsN2!AwMf6g&DW zoMKSSY&*LX-Cy(q0;!TWEyu8K&#|hngHgWO1xX*AijZYE(ihJ+yj4QTcr*3K$odi| z!6r?=6mJa7k1k2NTwzFZ=eG-~6KM~s$IG%pn;ctfro!*DUrCPJT9xsJ-F;w)W#4?W z*p)%A5efIZIVK+L{Yc{DMTSWSIF%t6gg?&MeA}fps`eWI;Gv@q*NB%b^fQ)eNNL!G zG^x}(@sf}cRIbX>(^c&1YvB+(>%tSTU-Lm3G3rC5uz8=(Ji9-a`nr?J=-tPP$j-T? zO7L(D)BpBleRiyL9lzAWYccsW%~;#3pG-mesS&2 zs^a)lkXM>12LwS{dhg%!&08fwPwus)i|xNbQpmi=T+95h`g}y7H2NUuf$LiQir)M z0d|3`%44TFwB)Ta#-20*kgIwRhSOJ$?7#p_e7V)5cT<9^-8NB4ExwEbSbg(4cA zz*c&v{d(zr&Fr#atf4kqe12g9@65f$*OJodaE?Nz5FH-zi{fWG)t=xwzBT$Yj}Hij z;-WuxXqyA?mo6rflvqb_9sLuYP%6puHp}azs*Zl3e;V6Wx)LNM9V7a({@rqfUFGQt zLaFh1n|3JV8h?G}m7U!tSk=lbzZH!!A0Ob2`)16qjoVu zR$f;Iqd}Zul!r`Xw5ktAPVC^)C@QgCusD7j>tH2e9FKx9%O}~w6E1aBfFZ!_4S;=9 z?;6gI*p&Cp#HV0=3Hq%Z>vowheesF9d1v5_M^Zhi%Ee+}=vRp)2`gTD@d)oMljetv z(@DnWh>JxDI+OK}@Mf)(9$BgD+8Ff=-&zSLXofF_e}!`7V3s#FnJItlG%OL{3me?6 zQ+ny-SN5zyLop(o`T*gkNR%k`6W}$?680-<3?id|g6N8h4ck#G88bB_jY398jg6#m zGv$4@PryI#&a$IS@M3K>^k6eLa0gLLG9|1#D+3>^hvG*sjS5bx2(LYzzi~nM^jyNF z(^3gtlm=h)T0TwDNcI}XQ&JsouA`j(Y+56($jGd()yw%3qO>ZH;4epK1_;&Y?hYK} zpi_1laVJWJNq!;o{E*>zD)fxJNg=Asrce<#?tSOM{}$rPN;-8T(7A+8&k8F(}}dX zd_)^UPi}g%+ED~JhW{v415zy6uAdh4`WzP5MV?ZnS8TZDzn;7}Fo;3x(kuDgRPS#- zvrTSJl*xZv-~3sq9ZeM}m8dr4P>L+ypbD>a(!R<_!7h(Fo}+H0E#9Y`>{&MyBzt|EIOeDpMIgryKmhsADt6HI#LxN`;ALz!iTHQ4TDH}S(=u|7Q z%TY#$H|iJM!&%HG3Vf^cwAom&gvK%DZk-R|1=3~ZLn$Y(i&#?nL&fGLeN@<~yx#T= zY#ZQ8{0po-O6K(vTxWFx-`U9g<+#$Q^Ehr%TrZ!|a_{31YpwDxaO3NJ0KJ}?{2}Zg zsEKpJ?WX%^RuQ_U;d@f$ywQ#^PG6(n79ZSnFch|bNNa+O04J!xW5I}%ft2Rifs~c3 z^1Af5A*P#lqmVM2`=v!Ni1buYL{FokrQP;o;mu2&iQ*XjxAY%pI^n!rdD0eZ%nw8U zN@ku#UX%Ab^&|cE0jy%W@wGH~$H-+HDG=iaK_kV|m*nR>t&fcCbj(irRznhV9%Db_hm&(a-Y zO6@}$HVxh{HNX5w(eq!Yn1q;!Dw&N5EYm|jZU5q%4W|S(#6qY;mP&t&{W`Y5h;zA} zqv~QaSV0c7D-LP-@$W5sOZGU35|PX+^cm&4pppsC0N}zCw8M`LRjLzDnu=Y*tWzvx z&a$FmEFgFHvWvE^sF!$ap%AVbBP6N`)fFvpAd?WAzFdB!WXpQxTg5Z4SL#7CSb8g zsI#79Za~`}7yB#5XrIbb%rd!D+oHi-k!CdpBK)&BB`Y%ZL+L^8L3rfQiZqHUp);ND z&WVYW9=f>*3Sx=^OJ{sWtT+kA%?{hyI;M_&UThY-8qycQPL>5IF9v4ol$#0(ryEBUhXugky zS0(7VeA}kf5%$SPL4qFLGZA?`mCBx^r-25C9&KkxenJVqK8)&6B#hK=qv~$@rT<}M z;Qnk%Ps1j~r*SAP0wr|-)>b&OKJ4D!G~Mg2fl>cpK`p(^N*N<=UnwkyTpC}R0c$Ct z9~B_HCH6Fp_q-#bM)h{FNpQ{~FbSvA6wSxBCYssVnUhXlGVel{q}N7^DY7hR+E#AW z4%C1Ms%}_4hZ0K+(vSl>P492^_jI(EQMrlDk}@W8th$iV*evMRL?LjT*;Oq?aN zdTZD<@4pB)o{<`)O7!9iq(okjuuto5@d2=~+9*t+RFWv%;T~{u*$?GMI4#+A_%Zm| zhPo)@!Lw|4iZ?k_jSYa#`MFU0^^b!6N$sC&BL4U+c!bu>1kutZ&m@q;GeHecXsnxR?hQZ6bGN3K$W?umUNQHZ#a+$CfIyhiM~}o}}>C ziJGNcHS3GUpzRzHV5Ef?u5aW~FW52y0(jT-sev{!Ka$cf;=|i|4fsyvHTX$KEb!|E znLR({y|}YKt8@j!8fvv58T4JW?gH&ioq81YBa7{i;lKqssqw}Mdj_|*Hw~~# zT77YLtHAcKplj>3kl(JU+c}RVXjKXDnU4!n9fBxa}ai!Yr%sR4p%r4)PX(rkoh zGY`?e)uJvcTu@SN`PL2^zW#oTV&{@G1$A7m{Mm(rP54-k=8@D6zbzTBbt31Tz9W}- ze9NI{YdpZbs2bQc3ezMVS{{H#+;a;MCHvC-;yt>C_U3uobO=Q-v~Hs9MVSwGN;hO4 zMBN^+hMsKutddL+aR@Z5OQ>-yU&N6N7dl)(fAYHdWsMh3jz_|BD)@%LbCa9Yui4GY zZ%-c7_y0>oUXrZ28Qqv0WUT+nxJ?pB?JmsgvG1G_{4D&U zysW_`BT|{@QPKWAxC4!+IEzyJ?Tk}QsHi@5uSz1|2J;?wK$ zu>zAdFwf%kWKjr2420E}-_NwN#4&wU4K7qwtb8C2Bn*c@vuyp$u;M+pI0{3W{K3B~ z!sw$^X!Kn zEqPa!`0|AopxSkLd5~>B{UEjT>x^%Nf?pVSrI;DFlxTjm7dYvUX_2wjU5PZ<@@amI z{1ps`h)L#M0Cn(nu!1MVlW~M?QbrT&M@}SlS#IVzDZBbO81Ea)PlULjp)^*0)XjmC z+AGKg2-vPlVn$1b5A3=&Xah$T?8}lZ|6gHe9&c4y?s3twI9nn(poGRGm?`O+2T&{o z@`_?NxR8L?05cVbWROH*DTpG-c!0OUB0+LM1qH`bPRP|1JStdu1rc!u2}Bb=q~7nl z_u2cL^*sCYx&Fb1-+IUAecoqY%RXn^aQP?0&Tz8dY=7vR?>_V9sfVs?^-=8+yIXeO zy5XsbJ2$NE^`nEIwHrNj@!oG_e>yPnJ6EoHvRk)_-Re56pLapWFd4JZy=u@~S08wH z-`8IHL&Lw1i6&((i%p#GTl`&upe_4((XdGI3djSK4@>pbWm-F`CnrUk(@ zH=mmfed^l34!L4R{qXyzJ-Ve~&$qITWA8~P96NgTWpgfQGi2r8Dad*Bsh-R3ed*L~ zYx=i1x943~cJKAzx-+)jy5oTZBi^csM@$;|`Pr|ud}-5&<=>h3HSsjwb`3C zjKA~Js|rK+f4cU#uhzc*$VEY6+Xu^!82{DAqfbBUuG6P|y7%$TqYnT2=qY!{*|bmU z@2vk|h%hmZTOIf);u_RcKYFhEkECR z+D8}9So*u4o;)_#{J^xsrbm0~&HJq${M*^9FTP{wr`w-swXgS#IU`oQx#oz&=YD_7 zzwY04U$$-gh7JpFo!QX)cekC{eR#JMj(4Bw{PnBM$D3PyH1Y@jx~u>B_twu(+%WjF z<#%3Q==|q>=9(XOea;(3ob$b(ez||-Gzw1UAwPnw`=nK3ud%B`rJ$I znR(%m+j@7pzV__SfBxUP{hJQm(sK6WgF8HS+kodzJhA?rcPB5}y>MCKm)XaJ# z>+H#OEBox&bH~?D=d;Vk7YB~)$q}7(0 zZw@-;)fcAx^;f?bH(>9oHOJrjT)%;Xw%+y92Y)>4$FG0;g-b@v|8>{Z{}?i`{++s$ zYc?<1eBC!jj_dQmyvJ4#Z)o+yRlBqI|N6y_*YBINJBvQN;hgMeceV{~=+$q~m<>JK zFQ5Fnw& zbLb&`zuR!<#ORngk3V@2(7?9q}K-I<+|S@SE1p>oN4kE~AhC z?SfIW-rd&gwO8-|``iw7oR$2@SbI-ocQth9=D&=+HJGuv&+_;@qVY>d*1lgc{Aq^_}ieXTiv&H z&4YD6tL?n&(eB+J?!0&P5wAS{%*cCQ-F>j1d-ui{KYMOj_w73`+*$w1rj;#M|KRh# zbl7>``|Db+`~GcTYc;;pHRoo-r+>WU=_z0BZ(sjtn`Q6TUA}4LjxT1_+<(~{7u{Ln97`>MsPj;$vZzffWBd#vW9)_s3GtnaYlCl0#!ve6UA4<2QH ze{4-B)3Uy1s(GSS;TTglvgRoB!v!_nAKqGDv-{@5OU)mCv45d^Tq(0m%)x659h?2L z)0D!k?UOE!{jbaAU25!a90XAsnT9Ea-t9yE$aPKUsfFIPo_a*#=)6KZbN$r9@wJ|# z9-Eh@@{~VLsP)xD-?W%ts5M9X=A0Rfs+So!o>^uesoivypRej{9 zaM8&(RqN#@VQMn_NHMG%v-gJTk0K7I^GtfFo9BjMV6K?SxrK_BZVe%`?3ltqPMlR}-@-Gq zXBEEJA~0+1Ph?)_$9SgD*0h^lIMRGPn;UoP>_U4p`Z}-+EY|lHusL8YLbGBv&0m?y zqpdt@5u5!MShy6R=PdxCIsF#Kjm=p5lbX6)n0Kc~2y z5pxR19qKt|I{!4I=kfe6b67>+C5-$wEzJXVJx9zf_BnQLvClw$%ma&PJ7X@Z^Ucz^ z<-xa<8t$?!4qZ&o?(>S%_nuc=))@N(DN6rZ%TH;XkLEE>U=Evafmw7rWj@_PGZ+A3 zU1M=b)Yhyk;V%GW;=@yKxC!QnS-OiA}-q6B?f%vwS{tUTT=#okPP#rf{S_ z%rKCL$pNRn5E_n$@RE0wn=QpA?K~;goz(N<#zW7_`rmb;gtygm>dW;|vSJd~_mWUz zfgA9Kd%kbv@>wl`#iufncV2HGL1gRm6lSnM!WkIo#EDYly76V33?8$TIiwB7m{9vggc84+HNb-ZZ zfovCIKm=W~9Md>q!p!;ISxlZCT;#-JcVO^lzHi_-F4kCX4on41-7?mxWc>!rA{mKx zu6fE@ghQ5cqkS`WX`u)2AvvaQDVKS)71tYiV&0m@-#+1QnOwRT27!>cjtAQ*%!^!c zFb?KUDwa;)G*W#7k5@s9TdZ7WI>S%pJyy1ux&jL$-;oCxMj=8Y-#3h6Zh^$ou)@YA zj?J7f_Cm?VA$Kno6FXkQ8uJ5!xry3w0=ZPX&=m`$ZxkisGCN2C_|$8$DiBIpXqVj(sSt&NJaZxsxvMFtcDJRtOe*l7uLD z$qQD_`C>q*l`8W4ShAAQecy{hd6|7as3I2Q@J$tU!RhQ&E@X*W9Y5s5dV*SfO+`AG zv}_vzq*)H^=B8EgO>r_7K;wxqDN3j&v z@Qc!)$kxb{lOWG{FFYsCXcy$Hx}rq&V_K?)aamAQ78v-bc+oEKgG8Qh;7i4g+;Eh< zb@--DCK-^A->^M7NrVBTFbo)lg;J2AqNCi7Ejz6 z5Xro#HmTs9APw`k%hDguUfM%Wh$@oTFu+|@JEp5zV4C`pTJW2?4J;+hB&kDmxf~D? zugDtWMu;7=lvtWR4i-a+)0;7Rh$+=5JDQXbnm}3>4jM;BUKlV^!YK4zc}wVv4VdQ- zZcZr7A?P2b4qrvcUlZamN*INaR?lGq)2O}?(k(w83R~IQksHR2R6-$=R-QY)nd0Cm zfHrax?x#2shf{atak7t!NmEurPJ}5XC&o6Do;~eCWqe=+E4dH_h&{Q$7%7m)G6(sl zu1@U0ula_wixl^B!_<|H{XnVR*sZoX-VN`PigSvi{QcxH)}yqmF()IlV`ZT>Lutu% zqRiFW5hH>n;6=s#dciq zT;HPlh!ZS1xs${Vq^x2oRWM(lR)Phi#tUM}o=NIt@`^i2>MCuUq}XR-K#WZ3@z6I= zLNdW_6-P&GBWx7dP3i;BfqRvo1bY(Z&SP{A)CI!Il<^ll>^K!v1xxdbf&LccyK;Fx zMx!X>IE|6*svX9T!s03z=c(EyVJ5D#lcK%FjuFjDD}zOntHHQ^`ANW%Dp*C?R@4Wf z2j!bsOgc{-Nvtpr@uf(O+7dT{MxLMo0SV7ZZ5KI3~edv{SKqFr^&9 zl<(@|eM)7U0V1Bos{H^TQeqLYBZ&ox6yGk|Lt1w(vQY)Q5L@WWqWr6vkZu66<-~!F zyvvD+Qpk2Bc&hY8e&_q5H;}<)EGs!cs2q5`3|25#icuTQuYxH#PrE=WV6afy!fu8q zqNG$1AiKn3Fr;NL32LF!Rl8J?9}8YOU&OGAil|8cvVQ`Ir78yuTxpHbE>6|NTu({C z5LOo|+Qs;6=AN`HW<%W|=KELoVh&i7A}5Jyj~G%IaaN5;DKf%1Q7jf>?>FC^Jr06~Xe2pF;#R>tski9syA|ec z~F8H9~3 z%~dzTCYGcIlMx8o;irg+qd3Z^BB4JGW4#khFUgp37|I9}p=;fjc$`udtVLsy<7J98 zM4cor#L%@l(^#HKpsG9)tg)P}S!{<^kW@e$l^ug7a{x0$iF-^rb)E4aq}zk5Mg~QY zI*F7dwl1BJ&Ne@RoYS3^^G&-rZ^}W16vJK|6AqZVO<&g6LV z7&?Wi^O$GFf?!OrR8q?$R8cX)%jN(eVp0xBlvzze7(qy4`l`E)fkqf!j7I=J&;5>* z(NG0l!BlY8T72?)LKzjjqa9*JPE2Y^Ehg~MTSawL5y2RX$^Zf*`XmRW5I`Q2iy0zr zljqt9yy%N#+hqSVngawrs*n)ZQd}_>lgPM zlnA_Pmuk#YkwjD{XD9o5p4~{1{$&m5;hW3QlSiOP%5RHF8wSj+g7N80F|m(v*u>OX z(Hx*NQBm12B$h2Oj3RaJTnL51T2iKPBq~(kq5lXlib;s1s2zhz_2~t;WKBCMZ3HnR zMKJADeiJMz2L3GlE94IlU6QN8q!B@9a;sv|*uqk6vo#zF)y{9uTgCIIp$cJIy2z*u zmxh_8!151F7h(mZUN#lO6$dzCMixQ7z^#7Uo#5ER!b->RTIYu<|4nb z{|foT6c0so4pWsm^g?zlkBJkmf~D%(k(O|6)dsU71cD6S^KCBr-<6E#^Al0}jiVp5S(ipCU3lqlF_|JfP>Pm(pk_m!QQ zM4meqAqho6nV4L+?2E?F`!>*@NSE@UaM;n(azL(_Goq!q6ECBL}Ro#$aSK`N@DPbPk4`P@Pi|6tuZ<>=Dt9e7TJI_zvmf=8LxS(^?+# zVlrSXTn(z0p{Jq(;Q_fjyt-y4N>Q;hUL-^Ic8*L20BufQ8t7n3eqsueD-t_)QmT(M zCS@vpy~f2f$xlouZkOmx#7OEs+EHPq+Noj~+67Xkf`$1w4D_d7L8%*xf^pr&BwnJ@ zD8;C%Tnx1)((Y!CSc*Hyp9Y!k>{ppsMF(t^opj77lOSO%7wstfoOW<#sa-{W0wUv8 zXBD-iok~9V_zKee1ygAmFqxJieui-*ZCx;xmjP3G8Jnw@ULr@LVz*EW_7|JOL-%Aj z9UHGnwl#`2FwGw}g)FZfA(CpCG+}L7Va4Dt`_F(ZQm)WOhy_N$pbTkPpF;5r;jOZc zEE|Y~bIGOW$HRUi*C{|OZ0MD5NE}(;cP@mI0eId@27oCUkWlCr3;)@sew4U{-J~M6 z)yvX`Vt2S^V$o#9q)|XU48pqHu6)MmNF@nv-cC{h>sK0O3`k(DJYXq=2^PtqEg_f0 z?E4*M%RAN=6`wp>5+zby(x(#YDjy>SimGk`4F52{1~4MM62sgpuC$E=7$jDGDGO>QehPounuU`k(DTnvQSnpi&iyJ^i}2r z#k652?JgD*R>2%a8Jjp4N5_Cfh??!rZVpUpa0+b8&SsP11H>2>oiPcd!9%JET_xM0 zl2l*fFY=gNZeWf_p+6#4`BY$7JYo_nDCH!qr3zND%PY)kqlLMNi}OaxaUl4mf_wNJ z>P-V|SqkLD0z(H?Dlqhq%yMufQMp}3@iDP15puv(sRoz|l_NYoS(HPNOl2>?Lg9{O z)WV^fh2#ZOMV4SHgkmvi)Pm_-g843G7{xLaJPOAFkZL$VLJGg7)}$R#f7OnH1c4!- z6=NbgE7!fO`0zQf%0y=yIg^VCxrvf7v9r|mM9h@WxWa^FvSgShqp}t~EGO1w$M~*P z&8l{)B$ala484P?lw3ez5MtiSzoys?U?i0sn5qZ{!{!l-A)%o-oxajpz{)oVBX%Yhh%;EuP*vuzA@HURmjcIlBc3$- zzo8$U?TuAzWoE+vbwrk?gAU zmcmw5^1}2g?!kCz$)kbED#y=zG>W%5E@^U5I-=ZLFqPc_Lyn8hkaTp>V_1)#in&bVT{ZiQ@yL)W}E#|uhUkg(POVgjoa$;aF)d+Sg@D+B<2+ zAec4kKnZC|YmCk)J#*O0g)#C?;V?p@@mA=|grs#hRsmr93^MXUEa!PI*S78AV1t6uC_x+{>z&$}TU6NN(345lys zU7iY-ta5=IlFAlGSjLen%v(Se1g0zQg`6a*mYh*9pdzVM*%DY&@VOipw-%~&%!#rq z4-zajZC*uX;i;)LB~mcZT)l+DRIzl=Y1gcFu;e*|;c<#Nc(cD$qSH)eVVP50Bpc%1 zhy^0dWLS^)YLMKzuga3sj%(U%4y&>-Xi4NZ1521HW~Qu2vji1oQe;qStzKckkfzvA z%z?#$Mi3w_nRGl*OT2L-+F=Y!OPY5!H2F~+((pxjrj60MCczmIiBt{oS)?*Y=zCJy zgORt)FOck$vNdhBdOnyA{Sm~TECvMQ28jXjrX?ex7F5;}ERp6N7}+5)9`>g%UwG#o z2lYA;7B9m+5pcl@SN4RBBfQC;qSsZh|R842%)wV$kqvfuyTdJ9K_t z0rPqcKD#Ct)}!vMURtLwS!lT!m`3sf{*KJ|qxVU53&9DfGLK-9)IDIzGqN>}GNu5g zlFMLCFYK1%0)Q$K0;XbvG5If9AP7dNM^J%~POfqlnL;~Nh?3!LDkBa|#rMFJIRK_V zp%L;ft-Lt|I1~j-(Ge^EnhYfl6W>Jc1olw7Bef#Wn8>UwyJS2S9RcIJ5al_{d#@Kx z`FR}m3(~%6-vPWl5%dp|e*I4G+dmqR{3!MYB>nmh=+ott|F_9`UHbH?89!mvxCy@= V_lpZ_tg>~ynjSs+pEIE5{{UD4u0;R< delta 38241 zcma&NWmH|uwl#{oy9Rf64-UZ{f&_v)1a~G9+}&M+ySoK1z1Es-X7X=ls>(+KWsPI;qKO+v`T!Wv>Zhr|dP0Q80eTYtLG;c`Edq6A@KY_ucsPk*;*F%hKE39v;fd*@(X^h= z)D+j|R=D%iFvX7r8Zg1CJnJt#MEQRWs3zJvd>vyctD|hNwL!>Sbg`+pmNM3BW*<(& zq%VJ2EB->zT9eZ;7;669-AIa2x;9P!iOjinCyo#W7nn&O^doliyXDb+(#5GO(w)n6 zzwqSvQZ|vn>iFhL_g$}LY})z&U=6(Nst$*C?LlX_{7%I`1vXJhd4JW~JP}i|qBj(| zM(lWcA$?u8{M>y+Uux|gZ%V1@eo!(w#rC{Df`1LKR?kExwKPXCVCGBpVG~xt^`j6D zMr98%Cg4t*;87h-wWhx;#wW>}<+Z4{gbhlD!f@FKZU}~aDF)Af1)|{Ea~*a$8Mozx z#vixcA5%McbJO2{WN`)~l9e=7v*cWjwvz6Tl9e0H&|d*oZ3UFvBd&vFS2 z?%-D9gfKe@ng3?1>q1A)qGVLiuzRcm?hkB@PtwM3UrZ9uDhW{%*tUPK_XpxXOTXo5 z<&)nM(jYC6{%F_Ae=C!vWa*9A3jFtaztef#)BXNQdGFH_--0aqFPl<-M%_+cl&jF_ z_K_(t1^%u-!lL5edc?`o2a41?1PLNl?)bGfL;g?`$aM9g;(teuv2Q|*_gBk)VW%sx zkEutGe%~JGnW{Bb_c}m&CH^*^bp4p2y_531chu8S)mZ+C1&O=jPljt$@2C=w)Gp5z z2kp7?QPVj(Fll>)|XH$V{j}$k6s+IM&pi z+c$`3N=MGnDoAd%U_HD9{w37E)2w^p3ej%M$aT8NH8L|X<8b<3X=2=N&FAK3m@3b! zLE}|uz6bJqzy7tYziYJjy|jd@!Tyg*f7a9Y zgH9R6Spk91O4?=49qALdE&SFAKz@AE-Bk}_>PrD9x?M_9Ws-a1!kKBLWBJ*%Hfrcf zTL)uOy_deGbw!3wdJ1mKfaJkgkWtYDoiWNc#Vyd^DIfGELKgQLgQ6OP{!vPQO7*WM zeWy)9d3J+v0(BMocNhIjqK`n}_DHW2sGT~y@GmS#%C z8hMfGW|1zKZdvf0jHj08#=izbM_T{U0KdEW&ldGIZjJ`3Ft~Cserxl)w61{P&Ha13 zey5>Ffo>oe^$=dr6pr%els_u(cUS)|W9vY)RW@;rWp$s3^#x6qHJuA)4Z-E(6^&Aw zmJ;Xv^1RdO_|T({C*~AN1c?8h>%CMjBzz%@r_63wEsJ(Q-l_uWsGFu~mG;ebrZTb7 z%fuAO(1jgJ-o)|$t8)GjQIM1h7Fhw4F~D5$vbW1ECM7mkI}MmMP7CtfprJMR^skor zOKVU5TRHtzwtBTVYMM0ajHTH_xr8Nc6z_!>M6iRpVhe#ZYsT|%`EiMzabtozNUc8( zz<*N!@BWMRjhuSJ4KPRMJXw{dbsu%H&}Rg!vR;2504(BCtd75&;_Mok1y#d8L`wWu zvr205#}P1Q#++=Y>}Z;wlh3G6@Q;+X)CnR1i+7ziCvOU2X~t{q4hW~etM$)w`;R?% zgavTr-SJiSTFSNzIGAJr>rgz|@F({lrf#oSTctOm12j7rf1XKy&Z?zCl0)vx@x~82 zU*~gHoqN==rNejv_#zhM7D11DS(www}l9y5c2)BV@e{~vJUZ?FXM ze{igKLx+VgMxAW2h8?&=HJ7=Guw=U)-gc;JGZ2-pP(@RGn1X| ze?g<4Ku7+Hn4SG^5&s_kop%TyOoj0;D(*kQF#leLJqc9@l;J;;;WmaygVBEYhgf2B zu&?C*5PKE@hW9@K|G#5!tN@#({)fEv3t$eP{vkGj7@QOdjFmY~sgD-WT#H(5Km|HG zMVJ*;;a!bGF>X=okZl+>$f*r>&qK3=pR-r?&f%t{EtPnm-=4>z3lfaxt_bag(H@LB z2#|tyzj0shlwICsUcR}JtlkqfVpmoe9S4vOieb4o@-Yp89Z}H2Zp*qiYF|*jAE;4# z*22}lS835@JQcMRYmESoU$(feR{Zq`8miJn*UNflH}1g%)xSDUD&rx#ZENOXAjFPsD&*G`IRM%IJl{@!rbt z6;k(Tckvow7B$2pv3fcY1nRw#i-h`7cY78L^_SN*0<&vxT?{oWXs=2n?Y>JPO%55_ zE^XR2(@OFZ{%^T<>{(dKwgFh$wW13ntu+BAIGd?VV%cA7*LP=5+bQlixZR`saT_Fw zV&9Tq&5mqsnAU+W-@h_d)h6}%A`}#DKJ^W+5HLj(QqT^@&Z^KN)llEy?VK@syWZq} zXxuu84{Y+0E5GCl({vthbG8#7iP{Mb7EWo^YpwAHmVyF{2pG_IHiVsD;A5{#qw{xU z!ZS_Y)E4jdX0UR;ZP>|jf}|Ro%0}OKYyA{?xOn@5c~%RYV>R5|Q^}zeeT@@3MDTXH zh;^V;K>2u&blD;37I_eJH$QGm8~{=MZRiRv-Bp~5uZgvudYt6yHs@o#+bvsTMVQCT z<^c23AUTRA??P?HAq_onVIM|~vq*7ZFxpE!?|Z`IQh+moe<=Cl@GG^r z|7(p_5KT#(FNy$RmU+RrV@XOPdEPi;3ar&8JAAECuU6zEw608WH=ha1n6lC+U6)*8 z$RXyN{PrpTU(&dcc7sqv!BpNiHI_N~EWAWcr2NtS0oqmvSYG&!1Hm}DzEI^3OF?9c3g z-|6AI!bo9_Qh=X9OcS4VE8-wDeL^NbUyNaig9ugSom2)B!c*A(HjV*4##v+t@$sCJ ze`m!>J;U9Jq8_?6j5VRGpjCch68z&(Gt)*>1QlUHnG*@)jK$p=@-}b*v;T1F0Zm= zY_+dbv~a42var8chWr~p?b+giyP=NgC@VX`yZ&^yl-^h5f&Z*C4z@i`V4r^3ksZ7E z38D?%c306{)GM(yef^7u^>~KqTm;ML^6XEQFZ%5VrE|Q~)f=$^Ri)Ol^p%}H3O@>0WeFowbC8VG z^4PRPR3OL64=F3bb%wKYl#*wmfKRP;OstvXFZWP-nyQ|T>E+dFg-s5S12q+QmN4HX z&&frbYbXLBc5nL!myU!N%5+m6fxB*zm;7n7Q$YK+Js@obk5{CsR+^)1th{=YLjQH#J>MUql@Lq)i%0Ye4+Mvl%EUQUvROfN4Y-ta`8DaH!Nm4 zfLM+hdHJzp5$_;IduCsJsrBH?T0Od8lqcs5oDId`o=buD)l^qQgAD;A@J-+2KA_v~ z;>WU#+Q{(ShmZ;B5v257Y#iyB^L-Lnk~}(Xe1v8QJyq0+OQ==2M37c*AzO)_47z4? zhvAm}p2x%$H$DXO?JywE4@yB#aUOnMjZBm{bmxEm6Un*xSBdgr(8!k#fjaC;-l)*9!YCU{xj{T6ZH6G1epYFLF?^>K%#;pw)!+6+0P;It-;T#*KD+>L z#?jF&#u#@O2Fp8K1UswnUYQ~ygB_^cS@?F5JKQ8VQ<(>3^NXRA`M8tZ@f(=3M5eV` zpLEYppy(Z5j*%4&k)5y($RuS~!hFWjROPE%#rwSC2Wj{B>ob5mEUlXaqVHg921c zK{o`jfz(P_DbaTJu5kD=A$W@y$b|-Q+bVm06Ou8C#lHxrm5(w|fXn}oFFyY;J zE>4nam6Q>LJE(B9G2 z(B}@h1mT#zMD74{>o>eD&E~D|aEN_t6y&DEj;Ljeh)%p1MyrJR7D9xHe5?6BLH0p` z3B8t|2u=oWYDv(M%-Z9BEsjErgIY8IUr6Fj;&r-2v&|i}HbMFgT{&#!47H()eSkp^ zSdn7(*_twNKn*WbI}&Hfk)6}|5Yz3sro%t{d2{obr1H*&K94d!GhbiM+||y7k&8@h zvr${D8JA|w)ML}G0h(+OA=YJNw5^UqrS~~}@#P1Oo;`T@A_fswJ0F^_Y@Fna!v+U} z2ImWD3SFOH)aXaO4bNwu3$L1&6q3Zk0OBN?+prlVhJd$)8I!iZgy5;c4jU{W-fMORyJlv&NxL17&+OJzVLx3g0b@a z`PPJM8MV`N_WS~8sGW^eb@lMdjX_=bOkHT^a=^-S8q}JrXh31NQ0>{k=J0q=Ax>+^ zx86RU&DPF$6ULx~A+Ag=u2fN9uXDnL&G`P_*4E^CQIvSVo08@>Z*NGEvZE~(&26`@ zJ+k_=-f?c?AQ*C%TkxL4TL-PAY%1K@q6Dl+oZgGqiwhZ4?#hpPE_BeTic)9fzPH5W z#?7uLS{PHWjlB$F$&siRSm}Y&Bc5W$aytt8%Vv8lG0)w~*Vo`xv|wkMW?GIc2_^!N zQwO@6V)_zbSWQ9@{<#9mmU$_G14@JhZ2k9yTyW;Y&$(#EZt49oh-Fa1xgij`g+M(9 z$bw&sc7Z(@MzaylukX=`?CjU1n=(;d_J|Q2H3IQpWHu9Y!bC^7EktqFB0(+dE0Hq} z;nZ^6Z$4KablXtz`knQK;L%Ffs*Z{D3at#KQxq4*L-3^&&AS8s+yiqQqS+uOb17px z)YVC;3|&o{ysY1>YP2?67hdh*>a;k+v0`wN&sF zm7Blb)ESJD>SKs;&uh2StNk<@ZMEq+VrlcK5x9x zSsfdh?!3<;=9;j2Mvf^tIo|-*i@@_w@M1qxP-m4-x?A2)U$s_l^E5QA3jMb*| zW01l^eJ+v%EVP1h>7{>SdJF@#vznD)0B4=cw#c287-qVXHOT@y;hU&{biw| zV(_Ao!p#ROaAEOhM0=+#USS-~`OO8Za+*spu#~ud+&d4!2rndW_Pf@-{wVAZ)DU$R z6xBN(9he>+=zSga>w3ievH%#Bj56I#`XHp%^ss!IohpTuD&V&*U#+!3%Y z2s)R|sF0=y5(H54RtV=V8eavf;q7jzsw=y<9E%Q?7Uw?P$9w&hc=zZz13<6gbx~2H z%A4Z2aA-4oz)p#7gb**s!fOJVyJ|ZY{NC)|@;Ugq1HMp=Dn8QaBLTw5bZzd}V|y;; zR?9NY)P3=sSP> zW|L06TI08JGd>>=CKUG%;%V2lcHx=+><&aCN`uWy_C>WdQ=o&ygNjwC$MM|@Ka2AN zwO--Vqo+bctEeA^1_5Gz!{R8B=T!4Ua2)(r;X76n5~@^i>9+MkB@YXRLrJ#OM0tgD zN^upe8Wx?fp>XnVI_w_7Q*P)+cV0OhAfmYyzLa@K{^;7|@Wds{K0@Li z7MKW{C~##GUhHK3`O<@ecy&aov1yC&JVfbk!UMd>6BBJH3V^c6^DkO-gRJ>oHs>NF z3Zgdcp%g0>X0eNxr6Cz!Jjw=8rFZ#S_87T`A93U{R|PR#Szk!y<%?Vn^8-$G;Az(A zzM5Hm4F2*JE9lb#(K_%Vu&b<5uoV77&+N_vF1nI&2d=wpS#mCO16kd$=-=eis%{lv9zr?fvN!0_I~@kTR6Y z5j_Xau&hucbXc(qj%zdyszpwz451{tX}}|+k3^Cvy$kpUDDy#ES(sygFVKd;$sF^v zLR!nwpH*)t8~+IUN#a1o&meF))$?|jo$k@{6Tj$X1q{-%(sqgnFPv$2<;SZR1qGK3 z*GME+S^i!OQqw^e1nc8>sWm`+!&H+1x}bF<(1t{$z3lA?^szBV6y}ntv5Ah%yEYwsH#2`cE zf?QsRrDR!rGqJxm@u;dKDKg!3sVmpQdnAnfB79iLv^YiQH$F03`EADQP6G_I7DAqu zyarHS4TD?ircFyUZTxGs+@WtS0jHuavWlvhA0v;_C{_$7*j<3 zW!QnCCHSNdyLzIRSg=Vu0q>r}+n##lu_#|57$uFtr#U~HOJ{1y8>6jCwi?hp zu2xP^{mK}L$ z=AcqivEwo|`}TQ;Y^v`0A!Wy;tedM0Z-VRIZKD`U3%SHMavkxZ7c%LCNlrTts9JBW z!%|2tyF&@+S#}Y*0haj<=hOsR138Xs75SlSw^@I#A7yL4re%*9R6-4Rm8D`lXA(0| zT}kJ6#g<3(tRpnV1pDO1O-UtD_MDsU<}5|nL_vth_GM~S`Nu?&(_l#WF36^z=&1t@ zJG1D=1qohSaj3hP&|v0lL>S`2??2-x4b*YOUW4k@ow+nrRs^|%l+r!vR#&GP9+L2VZoo4aPY~|EGU$x5Gbsb;{Fzg z>M|K|5*?W<;*XizlnP0Pq0)RrRWiN|p@_{e;F09fSn+--L_1-5xF}Tzh2*YqAK`esdfcPPV0Y(7447>1rIj*FWaqBYsBWh#V>5#6CJb{)9&I~iJ(AH8 z4BG;ROi*@m&M4&%3t#3@5USQ;7?|D=XLK!ilMuzgpHNwTJB`dd#BAJbj4Xe=ma)g7 zGEv1MLdC5z!6nIygR5b&a{qq9L(KmBO^xL@gq4klk%xmVPKo(L+y)|4(yBPPITQ;s zBL@#>(o_`;dXlaLI0-Z>&mSnx#y|;h1N{GZTO2HZGx>fG|J&rlZB#P=&n5e(V530r zs=u7O|BKnT6nvZHpE6I5f$vfMQ!wcrIR8J6c?1^%_a8GW2}5xF<(~du30|v1ApYg( z{$F5g8;EzTe~=)_9RyzEasY(&U-tC>gZWQzeF{V(!C%aO_jEb_%hUZ6+}r>`hVf6q z(jyQ9IR9~cIg%6?!2V|{$N$vPpTyp~5DKXOC${l7Hqih7geCn+iRlj6qxn~?oPV!l z&i}Q(_$*{9(_i4f8~hwx{|hV!0fiv;7x?e4GROaNwg)~zk^bdCbNs)tYUWTFe>u?q z2X;ONNwr~V-?<~X$bKjnS7gEkTPicARF z_uAT8&Qu{jo=ebXCd9%f+A27N3(Ps*sz3Us{rGB1^Tw#l!gF59o%eC|a=m(cUA#qV z*|s_gd^Q5a>G1H-011ds63S;z}MkhVsM=I+zE&8#)qL$kf(l z;V$LNlX~XC0{)n zn&5&oGL+Mjqy1fo8gMOhkT1kGiZa<@AYD%!weOK|(cihFAs`&}C0WBsX!2#&_4o^q;qRur&yFC)J2!+HX7HU3}>dSSl|-X}VEd`Lcg z&3gAMKkHHll8pBD`zjC{9Gt<2WnT4mG+_2Nk`X5{yGFYEyC=FKL3qJ^G7}q4fq>ax z_F!*vfc+a}nt%C=XF9Pwh^0d`{N%(q-UDy;LsN1qIB10;S9aFKyUfc3CL9%2*YNl> z=2y@`a64}`SjD_nV&%%L0xxKFF^sl1Ad|hBu zoMegq`ZlrrZcIS}>gF^w{OAkaHQfL6xg=f#hASiZ`}Nfl#^;?V zq|^B{Yp4}Ghq00tg87BvvqabXO+J`ieJunS&Km<78Io?KhF2R>`g^JC8}xgp^r4{P zV3V&9urP5M)dO0Y&<~3^E*@XtL#z5XExtOLf|jMXrN4KZFD+~G`uH+3Rvtg6 zf%%ock>#^5q8S)NH-MV8bR16*SllvgGA&Od!g*OwY(U?%+-F%@wk|;U3Mph$3~v+k zDKpP@aMrz^r*Z#)KB=@j&%9seVzL1)Tonz3t5&SO391GYNK3i>s(4hc&7YO95-&m6 zO4!oxMgLS_DCo+4DdSj0vb`C&)iop@OHexp`vWFBgN#L1_uEZ0`2dR7I8XG=aS~r1 zedkmWW(+~^6hhD@bo2)9=^0;X2SoujQ?Z!(T9)B1l`G6gwFL0TQuwwE+$Tn$U!REJ z3Rf+s*6uh{w;QDX>H({<2K=X$bA<+22G>_g9jyNPyH!8c&Kv)(`)c_L=oKXK>Wiu` z^m=*rNB7>fQMn&&5m?5u+bN{w$@Rp0iAhsow$u&XNsqM0T{lVnU}8 zYN~oDikG`3nH(siH`9;|+gLFG8yr#g5u+~jLS?`T&d={ioAq4x>qRG{JZKef>>rk# zl3uR9*g1G$6b=M~No3H2C}$9yG9sAs-y|%spD}*&V8Y)RHk@G6rFK-e4Yb6*Y^Nom zgFRZ)=X_lMzMnEw0K0~dRsSnZH_I*dOfcc%kwiM+%o4$k$-d^AcC+~jfTL&~EMSD4 zcNw0-&}ASFo3i2>7Z`i4Mf&)Zy_{A>M7@bG%8J*RiPkcMPG*)8y*1e`Aw`{m$u@V2 zvxFD>ow;9zPcj+Ggg4kvcY&6*0mYj~L;mioGu_iUi?YL+p{S4g)oBSLF*Dab8Nb%Y zb~ssmCsx#nb%%t(D4+XmzDPKGo!(nSta$}5uI_mQin!|2x#(0GE6CyjcV`9Q9` ztDZkN_hBKqNjvF#htK!1pkp3HokaDLJeeLI{_(=P7!#$B>|*}p$PMbGNbB@3Es3cd zu1pv1$hL1YGX(6s{ z!x74rd_Hq4q$$=cK+c4bv-~zzC$JHpEqBJOzRIdL)yKt`A!lv9)8}$@3?l|N?FjS! z3B(8xv?a=5r+=#;LMzepw@?-gm9)ldWX~zu^HOatrrj6XPj=7IO~pW~h2~QkV;uBB zH+@M=tI*SMyM72oeC}3H`8enpPE=^&#a~M9rnTl~q$fuVs36L*TFAen-GS^O1r+Zd z*&$IP%~>8g$eg6~$I@>1*@lWH(~F$QaP_J~6crU^J@(H!Dg|r~{`_MR_a`F~m2A~`qNW$XVK0HnpZrb^1;jr(}&T6HW6OC!@LBd)9dSc6c zpomNch;r37%a;ejUA^U@oX--}K(KD->zlkc{nk$)D!zW>Bf}WUY6+ea6tX$QG?}}X zdy<&8Y9t?{G+pmsuag!KTRoFxkI&#^>p_G3pRr(k7Dl_e?W{p;{{4iDXG_ zHiJcSFjYL*;eLL_BFVJ1|9JrTn0FMU4Lh>9sz%gZC)}Zjx2}HEho~dcF;LkcrIh=I zm$+U(tf!H&Ef6kt!5wRKw)ivYgKY;8O++674D#G#_nrMTz?tn8j>5_Kxbd^@jz<|b z!bfJYk_PF0k}WUA;7Jy}jQs`C_{Gxe5<4EjQ7!`)=|tp}kTEdftG7__!sCdvBz!QedV9BG`s&sHDhOSD$me{|5O z1E?Z)pphUGVC(GqGT|sj5Kg%ZL@Io(zqYues zTMg{+CH>N`ZrZ-%Y{t|X3GDaepU@F*qT9UZ1L=gOBhp#wT6s4k;kLa4KCdZs6%pS|^O^Z4 zF1jeniHC0jx!OjI4(gZ5hYB*qg{g(~jptjUBCLOfJy|Y)OFlPzT^lxH_r6x%tEV55 zk4^D{knx_d&xn)|NQ8iNihUJ(DQ+;g=6kVG_YARyg0#6|e)YQ@Mc&sgW2_K8BH$U81kksv?k>=j{1FaPDLGD%&aLQ`=P1evBQw*K~~I2;ff*ARkS2 zbo+V`E5m5c)?nV3NKu|Zq$b3p4R7Z!;GA06u#9aph>V#wxqHdxQ!N@yjT?Nf0$b5_ zFJaSo5qT>9>0a!I+2~#5=L&i{9e_kfKzI@A%3Yer(%r@!m5QnC81`DkXqPUGhJ7Ax zH(hIK8uN?ou_IanvY=C|BwY}iq_Ls3!PZGW8uNt2w|&0vZucYKNyj&?l z5nFfBWH0;3p0M^atf)A}0@STUDc80kdiw%FmG__h@n=lvg`XSLlJdcymP3UgP_o=O zf~GxRTyn)~b(F^3uM9Qi6FUmI547}tc{@J)P((&lmz;-Sf#+y%lFL3d4SD($5>=3eU zNq&1V5~h(3(jW)Dd3E9Bnm{Ie%~8#*$Ge+Fja%bH+R*31y1LFn39edm6h2kiIx7n(=Oep(j@rv^_AY z>SRRi#i@vlcW@-I24DxHnDyA~$2eH3jQvtp_k8?BAch}woJREIkT}C=U46D@AGZ~C z{1%mxmZjUhRd+OVcux?H5E=1xIAJa=L~~Z50%|S|;z_0dvFt zRA9Fff>&;GTk3vGPmexit9(qXcN-J>C^aN@(#)ot8z&Pbm;x+A=dWeTd@^CqrgEpUWyUmeMqR;nZ`AkKoHp1J7N%n^-4xqta~@soR(TN_kQdy$N2rf03VTk_MKK-?mhiyA44l56AY z__Uvk#D!gMXT^iN=8f`w;UMLqKifD6o(Hu9IFeshCW7aNoL7ENQ>QP{8O`)+*yfcq zm>3Y%4}vWbES8C9M+VHLcqUUCoBVutcLu?>dnpOb2L`WGe@32>RWFiV(N@B7OWS7> zq(%E2xqZ`*YE{^9y-AE$;yazl$Ca%u>Av}z8oz! ze&H>lD5CoOS)e|gCbLowzp23FEy(oxzWwG<7;lve@g&~=-%S_To znEx8Y9K^jP(4Sey>+g=|H+HoSMD}2>q;e9kn*=_SVd`!dwSO;{D029bdKsxKJuX&d zqa!a=c0j+L|I~ce7Nt`#ns(t0yJ0_&6n#*385P`WmP4jbMDH^=CbDC`z@w}tpZ=!Cxn*Mbd33M?@9cG>&LC_ z1a#u4w_?I+?*V=2a9Tx!9D&5DZy;E{Jkt08M0PZU?HVU_1H{iBy4Yq_n)FP{H~6cG zdnl#Cz4&|hdvY6X`Qp0E@{!*@!)mMjB5mbhd}?PvmLoN;C)!V~G_>fkmD$a~F$ zP#40(5fzG4^Gz!I+ZqkU)b)hoCQ6slum|$!tucjf?Aaq=(gr)pXJ zW>-IItx9`Zoxjj}CdVY8q-x53$nUMp4crMOA)1>yJ3WmD`OpQD?_MSKW5|<;HzXG* z%S_bN^j8^aD4iGG$cD$Z+O`2!IXK0xzb+87cbdBtyzg ziYPa@3qNEJ;+5-l3UjR=uY7ijr@CimuSzvLMCg~tgfBUMi8 z`RLG63Kn@?$+jrmlMr`f?-s+L(pcp3iqkAcuHEFc3Znd~N5Fe4;W&%@wbgTUus$Kc ze1wKabZ&;AtPrYSW?gM4xy#qgrM!nYIHj%pCb#Fk%7nNVhI!5JXg6LXnp`{5T2 zglUP*2L;962Ul_WU!Km$I-i*ULkt{~-xp5QiOW+|)18~lbB{eI1yV^c%vH%A!yn)xTYf_2 zDJZ`k{m77vniimQVUeYL{Kb=igc{;O1X|RuZY_{TTO@El?7IiKxfesZh*L z7#f0bW)eT^UR=YLoLB$8?SpVNW}<&nH;>E(OF=97&?bymxWeI*d+nV&u^&in9G-&* zDC?JO*>7WC6kmDj8`7IJZ7sA!x}_nUoK$aXamTBF?*HQM;x4ag_fgc#{yFgytmLyn zR00HUWH|{}hiLq~B!IvLQR~Twk{qKg8IBa@)p`AGSfe@!B|dYzU)W7$)F~?jqu}-1 zg&W%Thxv{QH%Ue#R_->NOWdmuD+`xNJxZaS88KCcS~z}CWI2&!@OCC8!cF#DQzDp^ z+!5U$e=ez1UO=S~=h#r98<(DxuN6xv27M9C4~0U+_3R5)&;eF+*BvC1Wa6D;&KLX> zZ*91E6!t!{1~Y@XDTsuWNpCQsIwc8&rYHLJ1&w2HCaDlQ++-4Vrf3=EMMv#(!6<_) zdh6DFEW#*tr_53avdWBPs$HeBTLv__CEW_8t4Xu2YW*if2z6fFw?*g&=Q>z27aqR| z=Qptp4j5QuegzaQ#2B_p7TNoumx;YbnhWqDmnDuyg@|Kd<}~UXXA%c=^Nr0@yk%1M z?6{|^*e-=0M#@8^GID1IxLNLKdOT+ezq*CJG~Kp1)>OjgwPG;ECJpBIXJvQCQPsK;7_sv^)$M6M`$X|@6ZM$@z*SyXoOXva| z?S5Z`gV~1ZJ0DfaAzb+oJK7?$XP_Jrs&%kUHM+(&O6PX(th$wl*2PS#5mo*eRu%=) zdFn!Wt0|z}e>_tuVCBtfZ^Y%tB z;`~DB9ob-{zk2^l4js*>cfeE<#TWHW3GIk2GW# zCRu!{2DDe`eMXADgB@m9qG5V+lSHJ-myOO~Urb2Eu~@A=x@pQgTjgf21afs0OR?M0 zv>5@n2ma<*V2@FAqc3kUMl|1TEs?b&#fT_2LZh+v)ddC9?)`3m#x+0Zu1vW(qRF&j z69hFDa7sNdQ|$>>0}_z5RQsD=G!n6J7G#op#Vf`4{Zc;elU%Di@bCr&x5yTJ_Xmx{ zsRds(qh2Q)3|)VC4X371CunuOp~10npBDmC^)Du9PIlo;dP5vk*nU>gm}5;ykWL9Ye}?(%w9!w9YS$9UlLtc!k4Ub-+)n$QA078s?H2Yz)9mPM?>Y6 z=4YY3+GaFjq;pVupzye%R%stdI%Mg1cSbii=h_3?Y{|p^0a%N4<&fp<3jQY^yhniX zq&ySTa_FbrQb;N&%HVFj@49F=ih75fSfc%bReUI5^G*k-TYVcdo$hhE;SDO-|!g(IHyEKd&(> zPAMmQHp{@Di3Ym7WVX_(^nv)|jsali*<+g-=DWgZRM1<1U!IL;SE{DyO}!FSR~vTf zBc5OK$Y90qm4+!ix)MtTQek#BVtjxBS}0@mYna&csGQv$6#NxzWi#R+{Qj=`@Aj(# zvW%Fi33y%A&|(-vYP_lS0_Uky^^rRq=R)2M@vGlQDEK6juzi_KJ4okkZ-A^vcL9ZE z>EjfoMTrA0u}35MV_b$$^(R%zl z@jj1q-@>P!qIuO2o^pE$^xS4us+qFz4j`jboT~^8U36`>+pv)llsi(pAaIef-65;# z$#Wrhb`sl9U1VIFBsjcZmw_ysxyeYop(C!ISe&J2YxJ_fC;I1@)O9ArO`r`OrT`A} zoEqmSad_q0#0V#j!NWP23NF5v%2K@d1KgUx4z%Do`EY)<`45iI9D+mc6dZkX#%rn8 z=j(1>Oo3)n+lTzhUPuB7kecYqrn=bbqB{B=8brsVD>krMzmg~@lR`2|r&i*^^mOHO6${hvO z$+x;~zkho9VLQytgbZd|7o(Bpt0Zqo42^!v0o};qb)jk>nrorfh(fGz9Ajf0HkKsm zH7nqWF#Gtng{&YbG66hAnzCLUY zu0OHUgXwV;aOI4(n3ZxlQ2c&D@_2?G7ZF}HlT^lAB7grJ?7jEQY};n%SZ+f3?3Dp( zj3HVl>Io^45tcO;1YfE{Z}I>|X2|RxkP#V<9qXirS5!61m<040l7D!DP2Qo#TtVw-RroHFH{p}!uC1A5rr6($ApI~NTJvKlg5Tl1 z9v8MzOYhg?3dH#b{vQEGV<$|c)k?Jbt*VQnzI<_LxoNwQ5?dXuzS_jT3gqFZ1Oj74 zSyO9Z4$UQ;L3Kb1`p=D)x#tyW7Lcjyl7NqU!9nC2mU`BfZ|r5zG7)&(SC-|cgdy#+ z8S^uN<+^9WJ1&y+?XYd^v$#)+E{Qnqw0&T2^s6Ffpi<#_Ud#c-@%c7Mh;?Xo0r-YC z-${8#LsQDsQaXp(ptb7Ypr}t_;a&2F*8m(&oKXjIg+L&xdG+qgD%I9Wdud8f(aum| zpk5Ei4=2A89^5HQS4cdlbc(4kXVqW`vV3ycbZUSsc1k{Tb~#&#RBUMJ_Ek~7e<m}2z&Lqib7b+ej}=Tl&hz2Xj0Qd`C{yPwAN5XnNyXR zv`WGB7d!;4DoW6F5N|K@Q!ijqQ#m%`cs32=-{XP3VJ?=y;`tv@lRM$X;WIh<1`C8{ zI-5+is%2*wE3ElhS0i=!4hlHHBgM-OWnAT5+_?yo0fH`AV{ z4pF~Yu4&39$!v&>iBXL^d{Z$w4KvG~hTRIOBO1Wb)f9*p zBe4L6SJmKU-?aHdd=hPgMCRnDqWg)Cv?5Hu>J&4%xTk9rSEu%4(mrAaWgQ!3V#b(r z%~#d?W66aviT;DMWMuyZL%-?|gYnX-I~GvI{*&9Q^t*XMJv#jRVQO=2trsuE9t#;J zGU8}t&DfCm!E@!7jNo>?nUfH+ufL+GzH|W%8v(vM5?!EZTZ|VS^m944GoWRe8Re6$&i$ILEy7_&Tl?mxcH)&Y(jZcwN;}`7bmPhFW=f1j~_(XaK zIZg!Bwq5Vf!5XK5X_6;8pE_jS2SvWT+N!+;FQ3|5?kLR8Q{7PIb6J=~adv z<{K{&LuTjC6}IIz@-UnPvRo_M>(pTm@p$9RjXnx5cnl6%S{L{0oWya80$}^fOmnfp ziY5hdhW)zJwi@;Y1u_tDH&@bfX?3wqW-h;G8rYhONp0w4xspe453huu&KJZ1K?WHa zG0ENN3sw^c^w+~{s0h)s+Y}nigkKO`d$rjVL=BJ&4!06CjA^u!cGqYo6tdANywux;lC(bSm`WSgz(o(_kZ{GO(aOy6;c8H7YJinl^Ybh<&jfdbuW#p?v3YyK4He z7m@9u6n|I}=G8*2XsIkw)d0*jK}u<`HgvxNUmcrvv3jcDyx`XXU}!b&`U+KIH7a4r zz0vFq)o&)7^JvrW8u6w}jHYukcv!_{`$@yBEs$qm-RH&|N9$HojDKGFlU1!G{XxoH za#CSYIjm(RMje{6N8o!=r*`zOox3&&+@$!oLF7Pt-%>C_)2$7@dujpES+8ELCAoNK ziKPbZ7D4(le2Yp%;DrjO`?hn)KJVs^>`EgWa~dl}=cx@wE^g(gFlz9Ab6@7iQ1Ats zTu+?T3v$s36IopG!qeEnpU;k9Gw;9yy5INZ2&pGOn!y^->?mNOhhyfV$Wk#|8+0TCPSeXil4{)j% znyxDVWM@$+*+7rW_-h@b*4Cx0#*?^feaMgZ+Q&A@hQ(a>8VRoMy>E@yIwU+V2&(U} z4wm?7z0W_!?B{P4dpfAj=>0GlL&UQDul-NT(kpIk(2VjknC{D=Ol?P`?UC$ZZO{t! zrJ+Z9zZ9~D1HC#73e!wLjlgNr&uzg3_JKZ;Lti7K&w_9-A5;@%rBcj}i3Yy6guzb+ z>1&J^++%ueKqOi?<5R6|YiWY(<`>*FL;u*ja}1h>|4w@6o}ts3==#K)pW>mA=K zF%)*nXraAxT~P4#l>gi#M1rqOjHws)(vP*msVL>P^_DM5HV3t*4`NO(~6V6d55>_A+a@gthGU_j*V;9Ih}POiH(by{wer?>)T?(=7RBlXKLN znV#V2Xkf;%)FZFq%%h20bds+wYqf6l*?q+rDI7xn&Q*P>Vc%m?x%)nvXh&^rJZYS! z7k4te8G6EEnFV>B{0z+skweQOB)POIUDb+75A3HsBN$UuiW?U2FUX%j ztMS{054%`!`Tq(#7icZZDv#%d((_>%prV0-fM`1O-Y=r4h}KX!qUeaABY9!GE0+`n zUr0cu!$nlynUO{ta4c_;qDEPXk(g1U=B>!fl2BJbpma~1-}8Ov`_6g)`>ZuC)_T_2 z?{nSHzVH2fzWnz;ZrJ-PtB1ZAEdAM*N6&ch%!OYZ{_0In{ZG#x?>zm`$Busfgw2Pp zyQc5N6W)Jl%4OH58*kb=X3_m;%w7KZgSISt@s;;i-#ljB1y8NbzxAv4PWa!^ix0nN z>r0!b|M=NCy)Rra@WQ@-yKd0hV#euhi~hFvD^FfC;gB0P_xS9pQEg+kZ+K<;J_EC* z%P-#X@LsE@oVMWeyLtr=Ubx|->Dz`pyzS}bYbTF;e$#=w{^_p=p1b)~~*??@y0?cEYgXJHOE9qJzIPddgn=o&VeoE3#9rI`QDQwjcY}py@le z-Z1Z)f1f<5_lwJaQJ!Cnn-)!7bbN8^UGpyAef`;w{A90pj$QtZt8QC;^=s$;`3oCf z-usuMo*e#Rdph9q<1V~$(VEfkE_h_?Z>DdJiajrR^n1Pkl1w{o=xJ+Tf9awze=I+^ zYy3MKr(ZeZ)n5G{Xur1qqK!*mXdg0Z{XGM(?=#_|UwrDYwcmd1%w@x_Su^mQ4|;5! z`LR#!b?(Z;E-eNvnEdjf32o!wKP|j_-8WAkaoCt)r`$K^(63J4d!I+vY<~2o&(7QX z$Uk3q?u;=*XFmSo_+MV}yG`SETye&wi|<_b*wUd79&_E8)kl7_*Im7*Z9D#?Z+-8} z!ET$MJ^Y}9?|U=eD}no*RI6r&VlV7$+^ltX<%X0_(_2BD| zTyya1ONuQQjCg!X&(R0`_|w0c9epr&%iu5dAOE8jS3dp2!*{)YTiePJi|>5r-b3zs z;_w}>3?JS12b(vIop{RVK^LCV@9jnH`=4~k)#L8_(TAIUHs+suzH#s9L0fme{&wvf zCw{n~7&~~{*N26(?`@x&|Msrizq9hXy+)mT;<5K``~F^6TyWNtfBe(u`tCY#;xAwQ z{i#LtiI=}|#2)$VN3WXP>rczSJ9qaV&w2Lwb^Tu3u=J+Ex4*D_>8z;_KD=~8zi;3A z-haLM#x*xy`}V7&{_xD9pZw0a%XS;I>yo3=i@y89aX)|dyJx&H?5z*q{d|wl9KP+X zHUB-YSHDT84_v?G+FH3{+RTTSJk{g1pG>;^_%Hu{`w>U=y>w1~z`BLe*uPzI$9+!~ zZA*s#?Cn2Kh*XZ#KY?0v<7Yu1l{>o3PX^vKI+J#@shos*WW-SOMmz20ws>5@0kx$cDD95Q^tabJ6H z$IKT-ulnMzpLt~b{Tt6;{n}@j|9#h`qu)L7m1(zZxqQ!2ZAaWV^^T8!VbYlY*zv(v z_nrUVeSbZ1|LcPHx4bj;sa?-~@}(EXtr+yf!4ICBpK|##UmLXh*Z02pqaL65)TMo= zmIL>x?N@&E{kEe%e%9He&pdl%>HnbZp0d{`YX?34{tw!Q&hG;5T@D#l+xPL2C)ajs z{)R*Dh2@uK)(+|@iFAKnIkUEGzhr>t{u^+`fHv2rQIz;`IpBuc@clybDGt5zD)-4V zpCT_QpS+=V^!|Qkz6E()?prX-x1R#PJkNcK%qKre%a!g^YCd^ExO7%`^M!eom#23B zDNBO#?pf?G%5}mpP2+OA8)mid4?(jFUiX{=XJN9 zpNCxCP2GOVsPdlG{iig|$_?(5!85-G1H0Ya?WdUQxypSqHwSNo<&!satBFCgA7);; z?=1{7DDgvx`Mg^g78oS^fd{|2Z@#&2-}l3%^BoNtT=`xYmZ#3I4c{#-+vnE~+nu8G zRr6~@ck{}Pt_jLLzfKdCN4dr?$GRpgXSgOVS2Td@8^E2eiOPQM^>3qS%F~~(eY%{{ zUORv>w>RKVxF#((wX;~5mpvEor(Z@38f#B*e|qH{*97GQ4hpUQ%KaIa|8PxG4qixu zICPCyPIXOCE~6=bQq+2tFD_FKxIQ8~&retFr_#=^HZpij`0 zQ*PlVb}p^$w?~viLnR&4GIC8^PH;_F&T)-jK2S2+gM|3gUQ`_EMR9p~C>P)b0eaP3 zsGkY*(l1a^%q<0JhrTVEYTEabEEWp!1K-m{{4mI3S=f*KNGZureEwB8=%+byumUSQ zKg*L)w141bK_-g?=(t>zLl6gXqGMAeY6lrM4;N%TwH%5J=n^QN*$=%<;TOSgA*_1_Z6xzFU&czf#sf{+KW;z(6LeK^BW8mYq()wxDIw9XoH!SKzv52Bt6x5G zXKlzHkW9JhPS`jpdoH7i%7~_D&DWUD1AUQ@2Nq_m^b2<;EmG*Xb;crat9xc_o|yR} zh_KPG#l;K`R8(-b+*k5|<$OU_oWx!r3&%+s26F3hlE<>?!i%98rQaB)%Vl9!;4peR ze!Wf>re7GF6T?liY5Gw#C|v9Xxvw@KuplwllxNCdUXXSJi{nT%CG#b%D{gXe8s=ca zoGw#Y#^5NH2~oYCaxpN2g}ML`+Akb16|@3NRVsU-pZhWxSQ^S$Hkg}(gnrh^I30u@ zh5`1XaB&!VnVilG1BhKP!(bJxUKk_dv|o&!Z(v;vnxeUz*ec2}@eDg5s@1{r!X#Ab z4J@`YFUe97DP9rA`~Z#aFbs_*4-=ffYCmWuRE3JLO#+?I2{Yv;ywePNtK+|^>^vMy zth}mV4fC`qEL!~{(64BX!YI>gibB=eo)>vtg>6yy^t}ijCil)y>EI>HW}*bE#K53` z1B+Fdp`XFO_A8J+<|GMnPdkQ*L1C0+ksX`sL7ZMQ##NNZf!mD3ewjXi_%tZ+nw|6-ydLR@r=QpEyd-k_Cpr1t7kUpI}pU3H+&>hJ1}iGtLQ|3>V&3Q7&0`l#PB;DW>_{2-Oqegg<0y$;V$MoLeFTsBu(o4*_Ayx zoUg(Mt}{{(ooKL?`sLO4)_NCCe^mDbD&nT|+akp@iVeMIzORgiKNgvOgiY0&g`paj za7EnU%ozrP9o3R!5q$gVroA+^)&NVWY2u}sA9WGV4RM2=)D>92qBINbZDpZFZWgJV zg}Z{0AREkboG$|_f^LAJdNE8G8M7fnM9S6i+lVF(;=={%u^_>{uujHZ#Z4Y{I&)16 z9!parUbrD?E_4hG@lZXTgQbcE6iUTDuH!#X`Bj1HU@ym0)jQ79v>RAvc)%G5swwi6 z<5kxnIu{qJJ^&&}OhJ&HJNvcj2Q43@J(=Vx@5NZV$ zC>C(mbLA2kHm$Nd_y?G&M+OXgKu+w#P8Q~?aqSat&lS}^Q4qNj4wr?65)C1Keo?2+ zflc8nJ~z!%F&@3cfw#R%)1bz&u>mms^uP{`T-5|2R zg@q1x(}g&Vl)~%QkOKYSG1(wzQ?cyhC8z-pEC@vd`bA8NM1ipfa?B!;k$Pmq)G$XH zK{PrxB{-rm7z{$M>z|gOi8J5ucAN#O#r!xUb7&@v49#F$sSgBY5!Y&sZMd6A;4D_w zNwKAL5xh+`rMMD;tE~y^?yHMV8&ezj6vkRzf+Fuq)b;y!QE{WiQt+xo zLpsz58fTg*v62TcSV;3+y#%hQFeVpCgxK(Qh83n_8P%jw9(*|?mVreos)2znhgHts zLOxfD5XLzSP-?8C&xKvow=qfOox4hks6-DZjh4iTXQ;GL%G>7Kpl` zVGIm!T44wpxz50g{h}LSpnfhL1C4(HGrUcX34=ivcIjc=j1MO^IIPqD!_*@%po0n2 z2yuZK6Uk$&NudR3v0Sm`l3SIl1OIe14lyv0qnkl8D608wX?3sz&Y>cm6dyM)^M&Z> zLV*}1YHYBXOvSBBp{bk)M&#ZuqH+N*lYm5{U(5ajn5rZMMiNIR3`pG>807(a5a-K^ zDkoHkQC^v^Kqq>cAubm+(z^>$(Te^ME8oJHv9gE@k{M?f81|6Ze%u z$q91eDhD}OZa58?kzWZ;nu>=6Ctl|xrWq?!C0SwY${@+OS2;}Fge+&KiXd0izNX?9dB?GP2x~m@u^u(lrjKXA%idpIK56|ZYg@{en)Fd7= zwPAqq>$h?`?9^^xrrLs2_=Qns7|3c@l12{H#RCTK%EG|Z(V`!uY5FB5Fic0Htfr&Y z)R_$pvg5=E#EjIa(ovR@OY6D@_;+<4MiNuDUG2~D1svFKnK3te_`i0iH; zV@;L9kpNBl$}OUyW2~QY#}84M6y@lFJh;Wee6_Lh?!m6^lj#)Un8_+QF{>y{V5}3E zO^FTAIBUQN7t12pWE!rd9|cC*&#E4(TPlHJw>~~khQTU>*@|Wu=$EM7O_2mntBiGl zVc~klR)hMDb&JTHj(WSQ`{6d|vcj?m$vqVX^h>*enS{411ZdDK=Wn6k#loxX%)-b+ z3-h(MfPM&k-5`ZcX45gcr%OMFe^-_sD7q3nB(pH0SiB0*nyRdFg?Z|+x{$5PIbeCA zwh4xzmNs5J(l~g38q(ANVJzZ229AsN9pgCV~w{1^_cMAcMzL1`vnQa0rp`DAjt{5Du$Y zRP}`^3sFLHMHeb0@jW#BNjy3nXe!l=7MMyATVYWX+|BRrveEWE?0FdFjF;0^-iut$OiGjL}4z|eg2AHb>L{Rm!qy6qXvz5vdJ`Q#+EmQM8f((EzQBs;Q~T zeGrb-`CNsZz(|WYiCM>gQv)&~sU_19`>qN@w00)vn#J2#(i$qDqJwg$>R4A?*To}f zq&5IjBW^mH_zL}SgG@imB}6Ei60UVP=a1=Tq9UchP9}^}6I-Wrin2=6k0hDEU>rl3 zDEAh>H?c(?u5MgA{I_;USWT4xFS84;IMn0rCc zn2+amSWPoDM5tyMFmqyHRb}V!50DU;E`oO>QjdP7Aej}6H%I6c z)aZPatJ{NQg&LkN?%I_=94$b=#)}7L(ijxznA&4tzUJ70p(W}Y?6m&?qG4p$E}kkx zh~kaZx!|P-qh%~gQCXJ5;nsHT2makN5#{5U44D!fFj%7NV8#Md#{{b%+eiBmqR=aW z+>k$}W6-Ic8YxQWGN?yHb_eTB#_~WyEn~()Hrfwj?buAq{-dLzewWTQC*hTYifT4v z0w&CttM}vH-Uy-dl7i{f00Cw^VPHHmmXo+=CrZ@>1yn9D7V%@o65SU2FUI>Z_etE@ zPKb;WS>Z~Bo6pReGKKj}&J7rOk&exjlb)+hD3mfNhF3qEbQv)>IRd2v8dqX=N-j-5 z6AC7c$V=}=zjEY;+W3MBBOQh&R2n1XrIxv?EH?5cAuynvo7nuGdvMdx=+2bs(2f8- zbFB%6jb1-R2p-QwpRb+YOR0eO$Mwg{8=tRzw$Mnud(xxU6vfPhT@?c<=C;CeYG11` z7^UjgJE579QEANAMaW^7s13+wj1qLui_FbZsi2-0%$TZHWvp0iPRLPmp=KN=164!S zlK`{qjqPKU8_5NA4Y3%^%QZJLD3q{x6Q%J>>=GVinayP4FLDrDHUW%#(y{PC2fvJO z5f(K*_9XVb0TTk02_aC!4ZIVd{bH?8VXAhqs7!@`dm*4nTyWOa9Igj885=dX-rTvv zqwXldpJ+(LE>;X?Dq^U*bzt1S!pL9Qe${jiWDIGOY?rZEI~DoltQd;JOi&tFuKqd3 zt8pZOnHLp|HKj=IX_E>D#zsXhl99pSOLuK6KN-w$)6W#`Vy>I#Tn^SD_Bt7y8#b?Y zT}O2qfst_1g|Rx!W-x(7j#B1>8Pjzv_{XfY1VZ9jB?hlC%;N;q2u4fS37HS_Rep3T z$j}yH@Vrn|^AJGIC#*J8c#8-zSruTEaLZ<}aaGLWskTCDEh~nf&cMWz!A{)@T5;T< z)@nasrW$}(C8kW-y)+T!QNY zCS|sCDlqfJ9GH;-z>GtIk4#8T4}t;FJX|x1MMsmicMi4+5MbufDKKL&5+LPA7Ib0v zQl$#Ri%&btMLMhqVAQ2mI|3u!sxV}=u$X%zu67_WQ&j@Y)nhsN*J-{0GC+1yR~aMeU2E0_=l zD7}rfQRlf@^OLiYxLzkmClTfMvyr1R|HKn`y*koAab!GVz%l=CmvaV;7|}Lm>cq)Y Xubq71MQx~DD(TvW3_0mvPHy`@m7h2F diff --git a/inst/doc/rResid.R b/inst/doc/rResid.R index 8f83c1d9..1ecc4b14 100644 --- a/inst/doc/rResid.R +++ b/inst/doc/rResid.R @@ -1,14 +1,38 @@ -## ------------------------------------------------------------------------ +## ----echo=FALSE---------------------------------------------------------- library(EGRET) eList <- Arkansas_eList + +## ---- echo=FALSE--------------------------------------------------------- plotConcTime(eList) -## ------------------------------------------------------------------------ +## ---- echo=FALSE--------------------------------------------------------- plotConcQ(eList,qUnit=4) -## ------------------------------------------------------------------------ +## ---- echo=FALSE--------------------------------------------------------- plotResidQ(eList,qUnit=4) +## ----echo=FALSE---------------------------------------------------------- +par(las = 1, tck = 0.02, xaxs = "i", yaxs = "i") +library(truncnorm) +x <- seq(-8,1,0.01) +z <- dtruncnorm(x,b=1) +x[902] <- 1 +z[902] <- 0 +area <- 0.01*sum(z) +z <- z/area +plot(x,z,xlim=c(-4,2),type="l",ylim=c(0,1.6),xlab = "log of concentration", + ylab = "probability density",main="Truncated normal density function\nfor the log of concentration, mean 0, standard deviation 1\ncensoring threshold is at +1 log units",cex.main=1.1) +polygon(x,z,density=40) +x <- seq(-8,-1,0.01) +z <- dtruncnorm(x,b=-1) +x[702] <- -1 +z[702] <- 0 +area <- 0.01*sum(z) +z <- z/area +plot(x,z,xlim=c(-4,2),type="l",ylim=c(0,1.6),xlab = "log of concentration", + ylab = "probability density",main="Truncated normal density function\nfor the log of concentration, mean 0, standard deviation 1\ncensoring threshold is at -1 log units",cex.main=1.1) +polygon(x,z,density=40) + ## ------------------------------------------------------------------------ eList <- makeAugmentedSample(eList) plotConcQ(eList, qUnit = 4, rResid = TRUE) @@ -20,12 +44,7 @@ plotConcQ(eList, qUnit = 4, rResid = TRUE) plotResidTime(eList, rResid = TRUE) plotResidQ(eList, qUnit = 4, rResid = TRUE) -## ----eval=FALSE---------------------------------------------------------- -# eList <- makeAugmentedSample(eList) - -## ----fig.height=8, fig.width=8------------------------------------------- +## ---- fig.height = 9, fig.width = 8-------------------------------------- multiPlotDataOverview(eList, qUnit = 4, rResid = TRUE) - -## ----fig.height=10, fig.width=8------------------------------------------ fluxBiasMulti(eList, qUnit = 4, fluxUnit = 9, rResid = TRUE) diff --git a/inst/doc/rResid.Rmd b/inst/doc/rResid.Rmd index 728c71a5..7e76d81a 100644 --- a/inst/doc/rResid.Rmd +++ b/inst/doc/rResid.Rmd @@ -15,25 +15,29 @@ vignette: > \usepackage[utf8]{inputenc} --- + +```{r,echo=FALSE} +library(EGRET) +eList <- Arkansas_eList +``` + # Introduction When censored data are present in a water quality data set, depicting them in any type of scatter plot is a challenge. This applies to the problem of plotting the actual data values or plotting the residuals. What follows is an example that shows the problem. -```{r} -library(EGRET) -eList <- Arkansas_eList +```{r, echo=FALSE} plotConcTime(eList) ``` The solid vertical lines, which represent the range of values that a given censored value could be, is very distracting in terms of getting a picture of the overall behavior of the data. We can see this even more if we try to look at the relationship of concentration to discharge. -```{r} +```{r, echo=FALSE} plotConcQ(eList,qUnit=4) ``` It is difficult to see the relationship between ammonia concentration and discharge. If we look at residuals from a fitted WRTDS model and discharge we also find it difficult to see if the pattern looks reasonable (a horizontal cloud of points centered on the zero residual line) or if there is some substantial curvature to the relationship. -```{r} +```{r, echo=FALSE} plotResidQ(eList,qUnit=4) ``` @@ -41,11 +45,36 @@ Here again, the plot is not very informative. What solutions might exist to res # The concept of randomized estimates of censored values -If we think about a fitted WRTDS model, what it is telling us is this: For some particular date and some particular discharge the model estimates a unique conditional probability distribution of concentration. That distribution is always a log normal distribution (that means that the log of concentration is normal) and that distribution of the log concentration has a particular mean (which in the EGRET package is always called `yHat`) and a particular standard deviation (which is always called `SE`). The `Daily` data frame has columns with these names, showing their value for each day of the record. However, for days on which there was a sample taken but that sample was reported as a "less than" value, we also know that there is an upper bound to what the concentration could have been. It is the reporting limit. *As an aside, all of the following discussion would hold for an interval censored value. EGRET allows for these but they are not common. For simplicity of the discussion we will only consider left censored data, but the code being used here considers interval censored data as well.* +If we think about a fitted WRTDS model, what it is telling us is this: For some particular date and some particular discharge the model estimates a unique conditional probability distribution of concentration. That distribution is always a log normal distribution (that means that the log of concentration is normal) and that distribution of the log concentration has a particular mean (which in the EGRET package is always called **yHat**) and a particular standard deviation (which is always called **SE**). The **Daily** data frame has columns with these names, showing their value for each day of the record. However, for days on which there was a sample taken but that sample was reported as a "less than" value, we also know that there is an upper bound to what the concentration could have been. It is the reporting limit. *As an aside, all of the following discussion would hold for an interval censored value. EGRET allows for these but they are not common. For simplicity of the discussion we will only consider left censored data, but the code being used here considers interval censored data as well.* + +So, what we can say about some particular censored value is that the log of the true value has a mean of **yHat** and a standard deviation of **SE**, it is normally distributed, but is constrained to be in the part of the normal distribution which is less than the log of the reporting limit. Such a random variable is known as a truncated normal random variable. Fortunately there is an R package entirely focused on the truncated normal distribution. It is called **truncnorm**. We can use **truncnorm** to generate a random number for each of the censored values and this random number will be drawn from the lower portion of normal distribution with the correct mean, standard deviation, and upper bound. The idea is to create these random values to substitute for the censored observations. We call them the **rObserved** values (the "r" denotes that they are randomly generated observations) and they reside in an augmented version of the **Sample** data frame in a column called **Sample$rObserved**. The figures below illustrate what a truncated normal distribution density function looks like. They show the density below the censoring threshold for two different examples. The **rObserved** values are samples from these density functions. The area under each density function is equal to 1. + +```{r,echo=FALSE} +par(las = 1, tck = 0.02, xaxs = "i", yaxs = "i") +library(truncnorm) +x <- seq(-8,1,0.01) +z <- dtruncnorm(x,b=1) +x[902] <- 1 +z[902] <- 0 +area <- 0.01*sum(z) +z <- z/area +plot(x,z,xlim=c(-4,2),type="l",ylim=c(0,1.6),xlab = "log of concentration", + ylab = "probability density",main="Truncated normal density function\nfor the log of concentration, mean 0, standard deviation 1\ncensoring threshold is at +1 log units",cex.main=1.1) +polygon(x,z,density=40) +x <- seq(-8,-1,0.01) +z <- dtruncnorm(x,b=-1) +x[702] <- -1 +z[702] <- 0 +area <- 0.01*sum(z) +z <- z/area +plot(x,z,xlim=c(-4,2),type="l",ylim=c(0,1.6),xlab = "log of concentration", + ylab = "probability density",main="Truncated normal density function\nfor the log of concentration, mean 0, standard deviation 1\ncensoring threshold is at -1 log units",cex.main=1.1) +polygon(x,z,density=40) +``` -So, what we can say about some particular censored value is that the log of the true value has a mean of `yHat` and a standard deviation of `SE`, it is normally distributed, but is constrained to be in the part of the normal distribution which is less than the log of the reporting limit. Such a random variable is known as a truncated normal random variable. Fortunately there is an R package entirely focused on the truncated normal distribution. It is called `truncnorm`. We can use `truncnorm` to generate a random number for each of the censored values and this random number will be drawn from the lower portion of normal distribution with the correct mean, standard deviation, and upper bound. The idea is to create these random values to substitute for the censored observations. We call them the `rObserved` values and they reside in an augmented version of the `Sample` data frame in a column called `Sample$rObserved`. What is most important to understand is that these randomly generated values are never used in any of the WRTDS computations that result in estimates of daily concentrations or fluxes, or annual average concentrations or fluxes, or any trends. These randomly generated values are created strictly to provide a more-easily interpreted set of diagonstic graphics. Here is an example of a plot of concentration versus time using this approach. +What is most important to understand is that these randomly generated values are never used in any of the WRTDS computations that result in estimates of daily concentrations or fluxes, or annual average concentrations or fluxes, or any trends. These randomly generated values are created strictly to provide a more-easily interpreted set of diagonstic graphics. Here is an example of a plot of concentration versus time using this approach. -The solid circles are the uncensored observations and the open circles are these `rObserved` values. One thing to note is that if these values were to be generated a second time the graph would look slightly different. There is no unique set of values that will be plotted. The random numbers used will be different if we go back to the same data set and generate them again. Below we show such a figure and then redo the graphic with a second set of random numbers. The code is shown with the graphics. +The solid circles are the uncensored observations and the open circles are these **rObserved** values. One thing to note is that if these values were to be generated a second time the graph would look slightly different. There is no unique set of values that will be plotted. The random numbers used will be different if we go back to the same data set and generate them again. Below we show such a figure and then redo the graphic with a second set of random numbers. The code is shown with the graphics. ```{r} eList <- makeAugmentedSample(eList) @@ -59,11 +88,9 @@ Careful examination of these two figures reveals that the black dots are exactly # Randomized residuals -Many of the diagnostics we use in EGRET are plots of residuals (on the y-axis) versus predicted value, or time, or discharge (on the x-axis). Once we have this `rObserved` value, we can compute a randomized residual which is called `rResid` and is stored as `Sample$rResid` in EGRET. The computation is this: - -rResid = rObserved - predicted value +Many of the diagnostics we use in EGRET are plots of residuals (on the y-axis) versus predicted value, or time, or discharge (on the x-axis). Once we have this **rObserved** value, we can compute a randomized residual which is called **rResid** (again the "r" denotes that they are randomly generated residuals) and is stored as **Sample$rResid** in EGRET. The computation is this: -All three variables in this equation are in log concentration units. +rResid = rObserved - predicted value. All three variables in this equation are in log concentration units. We can look at our residuals plots in the following manner (using the second set of rObserved values computed above). @@ -74,29 +101,25 @@ plotResidQ(eList, qUnit = 4, rResid = TRUE) # Details for how to include random residuals in your computations -In order to use the `rObserved` and `rResid` in making graphs in EGRET the process is the following. +In order to use the **rObserved** and **rResid** in making graphs in EGRET the process is the following. -1. Bring in your data and create the `eList` in the usual way: `eList <- mergeReport(INFO,Daily,Sample)` -2. Run `modelEstimation` in the usual way: `eList <- modelEstimation(eList)` -3. Then augment the `Sample` data frame that is inside of the `eList` with the command: +1. Bring in your data and create the **eList** in the usual way: **eList <- mergeReport(INFO,Daily,Sample)** +2. Run **modelEstimation** in the usual way: **eList <- modelEstimation(eList)** +3. Then augment the **Sample data** frame that is inside of the eList with the command: -```{r eval=FALSE} -eList <- makeAugmentedSample(eList) -``` + **eList <- makeAugmentedSample(eList)** -4. To produce any one of the following graphics, using the `rObserved` or `rResid` values simply add the argument `rResid = TRUE` to the call to the graphical function. Note that it doesn't matter if what is being plotted is an observed value or a residual, the argument is always `rResid = TRUE`. The original option of showing the vertical lines for censored values remains available as the default for any of these functions. The call can either say `rResid = FALSE` or just not include `rResid` in the argument list and the graphs will appear without these random values and without the open circle/closed circle symbology. The functions where this approach applies are these: `plotConcPred, plotConcQ, plotConcTime, plotConcTimeDaily, plotFluxPred, plotFluxQ, plotFluxTImeDaily, plotResidPred, plotResidQ, plotResidTime`. It is also available in the two multiple plot functions: `multiPlotOverview`, and `fluxBiasMulti`. +4. To produce any one of the following graphics, using the **rObserved** or **rResid** values simply add the argument **rResid = TRUE** to the call to the graphical function. Note that it doesn't matter if what is being plotted is an observed value or a residual, the argument is always **rResid = TRUE**. The original option of showing the vertical lines for censored values remains available as the default for any of these functions. The call can either say **rResid = FALSE** or just not include **rResid** in the argument list and the graphs will appear without these random values and without the open circle/closed circle symbology. The censored values or censored residuals will be shown as the vertical lines. The functions where this approach applies are these: **plotConcPred, plotConcQ, plotConcTime, plotConcTimeDaily, plotFluxPred, plotFluxQ, plotFluxTImeDaily, plotResidPred, plotResidQ, plotResidTime**. It is also available in the two multiple plot functions: **multiPlotOverview, and fluxBiasMulti** +For example. -```{r fig.height=8, fig.width=8} +```{r, fig.height = 9, fig.width = 8} multiPlotDataOverview(eList, qUnit = 4, rResid = TRUE) -``` - -```{r fig.height=10, fig.width=8} fluxBiasMulti(eList, qUnit = 4, fluxUnit = 9, rResid = TRUE) ``` -# A final thought +# Two final thoughts -In the EGRET User Guide (http://pubs.usgs.gov/tm/04/a10/) the distinction is made between the graphical methods that are used to simply descirbe the data (and these graphics shown in `multiPlotDataOverview`) as distinct from graphical methods that are used to describe the WRTDS model of the system (such as the graphs in `fluxBiasMulti`). If the `rResid = TRUE` option is used with `multiPlotDataOverview`, or other graphical functions that normally don't depend on the WRTDS model, they now become a hybrid, because they are using the WRTDS model to generate the random values used in the graphs. Take a graph such as `plotConcQ`. With `rResid = FALSE` it is a pure representation of the data. No assumptions are being made. But, when `rResid = TRUE`, it is now a representation of the data which is partly based on an assumption that the fitted WRTDS model is indeed a correct model. The fitted WRTDS model is partly determining the placement of the random values that are less than the reporting limit. If the analyst wants a "pure" representation of the data without any assumed model for graphing, then the `rResid` option should be set to `FALSE`. +In the EGRET User Guide (http://pubs.usgs.gov/tm/04/a10/) the distinction is made between the graphical methods that are used to simply describe the data (and these graphics shown in **multiPlotDataOverview**) as distinct from graphical methods that are used to describe the WRTDS model of the system (such as the graphs in **fluxBiasMulti**). If the **rResid = TRUE** option is used with **multiPlotDataOverview**, or other graphical functions that normally don't depend on the WRTDS model, they now become a hybrid, because they are using the WRTDS model to generate the random values used in the graphs. Take a graph such as **plotConcQ**. With **rResid = FALSE** it is a pure representation of the data. No assumptions are being made. But, when **rResid = TRUE**, it is now a representation of the data which is partly based on an assumption that the fitted WRTDS model is indeed a correct model. The fitted WRTDS model is partly determining the placement of the random values that are less than the reporting limit. If the analyst wants a "pure" representation of the data without any assumed model for graphing, then the **rResid** option should be set to **FALSE**. -Also, with figures such as shown in `fluxBiasMulti`, there is a kind of circularity in the logic. The circularity is this: we are using the graphs to assesses the adequacy of the model fit, but we are using the model to estimate some of the observations. This circularity is not a fatal flaw to the approach, but is a reality that the user should consider. On balance, the authors think that using `rResid = TRUE` in all of the plots to which it applies is a beneficial approach because it enhances the ability of the analyst to interpret the figures. But, we reiterate here, the choice of using the random approach or not has no bearing whatsoever on the quantitative outputs that the WRTDS method in the EGRET package produces. +Also, with figures such as shown in **fluxBiasMulti**, there is a kind of circularity in the logic. The circularity is this: we are using the graphs to assess the adequacy of the model fit, but we are using the model to estimate some of the observations. This circularity is not a fatal flaw to the approach, but is a reality that the user should consider. On balance, the authors think that using **rResid = TRUE** in all of the plots to which it applies is a beneficial approach because it enhances the ability of the analyst to interpret the figures. But, we reiterate here, the choice of using the random approach or not has no bearing whatsoever on the quantitative outputs that the WRTDS method in the EGRET package produces. diff --git a/inst/doc/rResid.html b/inst/doc/rResid.html index eda770b6..1b4d81b1 100644 --- a/inst/doc/rResid.html +++ b/inst/doc/rResid.html @@ -69,7 +69,7 @@

Using Random Residuals for Censored Data in EGRET

Robert M. Hirsch and Laura A. De Cicco

-

09 May, 2016

+

17 June, 2016

1 Introduction

When censored data are present in a water quality data set, depicting them in any type of scatter plot is a challenge. This applies to the problem of plotting the actual data values or plotting the residuals. What follows is an example that shows the problem.

-
library(EGRET)
-eList <- Arkansas_eList
-plotConcTime(eList)

The solid vertical lines, which represent the range of values that a given censored value could be, is very distracting in terms of getting a picture of the overall behavior of the data. We can see this even more if we try to look at the relationship of concentration to discharge.

-
plotConcQ(eList,qUnit=4)

It is difficult to see the relationship between ammonia concentration and discharge. If we look at residuals from a fitted WRTDS model and discharge we also find it difficult to see if the pattern looks reasonable (a horizontal cloud of points centered on the zero residual line) or if there is some substantial curvature to the relationship.

-
plotResidQ(eList,qUnit=4)

Here again, the plot is not very informative. What solutions might exist to resolve this problem of the graphical representation of the censored data?

2 The concept of randomized estimates of censored values

-

If we think about a fitted WRTDS model, what it is telling us is this: For some particular date and some particular discharge the model estimates a unique conditional probability distribution of concentration. That distribution is always a log normal distribution (that means that the log of concentration is normal) and that distribution of the log concentration has a particular mean (which in the EGRET package is always called yHat) and a particular standard deviation (which is always called SE). The Daily data frame has columns with these names, showing their value for each day of the record. However, for days on which there was a sample taken but that sample was reported as a “less than” value, we also know that there is an upper bound to what the concentration could have been. It is the reporting limit. As an aside, all of the following discussion would hold for an interval censored value. EGRET allows for these but they are not common. For simplicity of the discussion we will only consider left censored data, but the code being used here considers interval censored data as well.

-

So, what we can say about some particular censored value is that the log of the true value has a mean of yHat and a standard deviation of SE, it is normally distributed, but is constrained to be in the part of the normal distribution which is less than the log of the reporting limit. Such a random variable is known as a truncated normal random variable. Fortunately there is an R package entirely focused on the truncated normal distribution. It is called truncnorm. We can use truncnorm to generate a random number for each of the censored values and this random number will be drawn from the lower portion of normal distribution with the correct mean, standard deviation, and upper bound. The idea is to create these random values to substitute for the censored observations. We call them the rObserved values and they reside in an augmented version of the Sample data frame in a column called Sample$rObserved. What is most important to understand is that these randomly generated values are never used in any of the WRTDS computations that result in estimates of daily concentrations or fluxes, or annual average concentrations or fluxes, or any trends. These randomly generated values are created strictly to provide a more-easily interpreted set of diagonstic graphics. Here is an example of a plot of concentration versus time using this approach.

-

The solid circles are the uncensored observations and the open circles are these rObserved values. One thing to note is that if these values were to be generated a second time the graph would look slightly different. There is no unique set of values that will be plotted. The random numbers used will be different if we go back to the same data set and generate them again. Below we show such a figure and then redo the graphic with a second set of random numbers. The code is shown with the graphics.

+

If we think about a fitted WRTDS model, what it is telling us is this: For some particular date and some particular discharge the model estimates a unique conditional probability distribution of concentration. That distribution is always a log normal distribution (that means that the log of concentration is normal) and that distribution of the log concentration has a particular mean (which in the EGRET package is always called yHat) and a particular standard deviation (which is always called SE). The Daily data frame has columns with these names, showing their value for each day of the record. However, for days on which there was a sample taken but that sample was reported as a “less than” value, we also know that there is an upper bound to what the concentration could have been. It is the reporting limit. As an aside, all of the following discussion would hold for an interval censored value. EGRET allows for these but they are not common. For simplicity of the discussion we will only consider left censored data, but the code being used here considers interval censored data as well.

+

So, what we can say about some particular censored value is that the log of the true value has a mean of yHat and a standard deviation of SE, it is normally distributed, but is constrained to be in the part of the normal distribution which is less than the log of the reporting limit. Such a random variable is known as a truncated normal random variable. Fortunately there is an R package entirely focused on the truncated normal distribution. It is called truncnorm. We can use truncnorm to generate a random number for each of the censored values and this random number will be drawn from the lower portion of normal distribution with the correct mean, standard deviation, and upper bound. The idea is to create these random values to substitute for the censored observations. We call them the rObserved values (the “r” denotes that they are randomly generated observations) and they reside in an augmented version of the Sample data frame in a column called Sample$rObserved. The figures below illustrate what a truncated normal distribution density function looks like. They show the density below the censoring threshold for two different examples. The rObserved values are samples from these density functions. The area under each density function is equal to 1.

+

+

What is most important to understand is that these randomly generated values are never used in any of the WRTDS computations that result in estimates of daily concentrations or fluxes, or annual average concentrations or fluxes, or any trends. These randomly generated values are created strictly to provide a more-easily interpreted set of diagonstic graphics. Here is an example of a plot of concentration versus time using this approach.

+

The solid circles are the uncensored observations and the open circles are these rObserved values. One thing to note is that if these values were to be generated a second time the graph would look slightly different. There is no unique set of values that will be plotted. The random numbers used will be different if we go back to the same data set and generate them again. Below we show such a figure and then redo the graphic with a second set of random numbers. The code is shown with the graphics.

eList <- makeAugmentedSample(eList)
 plotConcQ(eList, qUnit = 4, rResid = TRUE)
-

+

# now do it all over again
 eList <- makeAugmentedSample(eList)
 plotConcQ(eList, qUnit = 4, rResid = TRUE)
-

+

Careful examination of these two figures reveals that the black dots are exactly the same in both, but the open circles are different between the two. In looking at these kinds of plots we are not necessarily looking to see what actually happened on a particular day, but rather to understand the pattern of the relationship between the two variables being plotted.

3 Randomized residuals

-

Many of the diagnostics we use in EGRET are plots of residuals (on the y-axis) versus predicted value, or time, or discharge (on the x-axis). Once we have this rObserved value, we can compute a randomized residual which is called rResid and is stored as Sample$rResid in EGRET. The computation is this:

-

rResid = rObserved - predicted value

-

All three variables in this equation are in log concentration units.

+

Many of the diagnostics we use in EGRET are plots of residuals (on the y-axis) versus predicted value, or time, or discharge (on the x-axis). Once we have this rObserved value, we can compute a randomized residual which is called rResid (again the “r” denotes that they are randomly generated residuals) and is stored as Sample$rResid in EGRET. The computation is this:

+

rResid = rObserved - predicted value. All three variables in this equation are in log concentration units.

We can look at our residuals plots in the following manner (using the second set of rObserved values computed above).

plotResidTime(eList, rResid = TRUE)
-

+

plotResidQ(eList, qUnit = 4, rResid = TRUE)
-

+

4 Details for how to include random residuals in your computations

-

In order to use the rObserved and rResid in making graphs in EGRET the process is the following.

+

In order to use the rObserved and rResid in making graphs in EGRET the process is the following.

    -
  1. Bring in your data and create the eList in the usual way: eList <- mergeReport(INFO,Daily,Sample)
    2. -
  2. Run modelEstimation in the usual way: eList <- modelEstimation(eList)
    3. -
  3. Then augment the Sample data frame that is inside of the eList with the command:
    4. -
-
eList <- makeAugmentedSample(eList)
-
    -
  1. To produce any one of the following graphics, using the rObserved or rResid values simply add the argument rResid = TRUE to the call to the graphical function. Note that it doesn’t matter if what is being plotted is an observed value or a residual, the argument is always rResid = TRUE. The original option of showing the vertical lines for censored values remains available as the default for any of these functions. The call can either say rResid = FALSE or just not include rResid in the argument list and the graphs will appear without these random values and without the open circle/closed circle symbology. The functions where this approach applies are these: plotConcPred, plotConcQ, plotConcTime, plotConcTimeDaily, plotFluxPred, plotFluxQ, plotFluxTImeDaily, plotResidPred, plotResidQ, plotResidTime. It is also available in the two multiple plot functions: multiPlotOverview, and fluxBiasMulti.
    2. +
  2. Bring in your data and create the eList in the usual way: eList <- mergeReport(INFO,Daily,Sample)
    3. +
  3. Run modelEstimation in the usual way: eList <- modelEstimation(eList)
    4. +

  4. Then augment the Sample data frame that is inside of the eList with the command:

    +

    eList <- makeAugmentedSample(eList)

    5. +

  5. To produce any one of the following graphics, using the rObserved or rResid values simply add the argument rResid = TRUE to the call to the graphical function. Note that it doesn’t matter if what is being plotted is an observed value or a residual, the argument is always rResid = TRUE. The original option of showing the vertical lines for censored values remains available as the default for any of these functions. The call can either say rResid = FALSE or just not include rResid in the argument list and the graphs will appear without these random values and without the open circle/closed circle symbology. The censored values or censored residuals will be shown as the vertical lines. The functions where this approach applies are these: plotConcPred, plotConcQ, plotConcTime, plotConcTimeDaily, plotFluxPred, plotFluxQ, plotFluxTImeDaily, plotResidPred, plotResidQ, plotResidTime. It is also available in the two multiple plot functions: multiPlotOverview, and fluxBiasMulti

+

For example.

multiPlotDataOverview(eList, qUnit = 4, rResid = TRUE)
-

+

fluxBiasMulti(eList, qUnit = 4, fluxUnit = 9, rResid = TRUE)
-

+

-
-

5 A final thought

-

In the EGRET User Guide (http://pubs.usgs.gov/tm/04/a10/) the distinction is made between the graphical methods that are used to simply descirbe the data (and these graphics shown in multiPlotDataOverview) as distinct from graphical methods that are used to describe the WRTDS model of the system (such as the graphs in fluxBiasMulti). If the rResid = TRUE option is used with multiPlotDataOverview, or other graphical functions that normally don’t depend on the WRTDS model, they now become a hybrid, because they are using the WRTDS model to generate the random values used in the graphs. Take a graph such as plotConcQ. With rResid = FALSE it is a pure representation of the data. No assumptions are being made. But, when rResid = TRUE, it is now a representation of the data which is partly based on an assumption that the fitted WRTDS model is indeed a correct model. The fitted WRTDS model is partly determining the placement of the random values that are less than the reporting limit. If the analyst wants a “pure” representation of the data without any assumed model for graphing, then the rResid option should be set to FALSE.

-

Also, with figures such as shown in fluxBiasMulti, there is a kind of circularity in the logic. The circularity is this: we are using the graphs to assesses the adequacy of the model fit, but we are using the model to estimate some of the observations. This circularity is not a fatal flaw to the approach, but is a reality that the user should consider. On balance, the authors think that using rResid = TRUE in all of the plots to which it applies is a beneficial approach because it enhances the ability of the analyst to interpret the figures. But, we reiterate here, the choice of using the random approach or not has no bearing whatsoever on the quantitative outputs that the WRTDS method in the EGRET package produces.

+
+

5 Two final thoughts

+

In the EGRET User Guide (http://pubs.usgs.gov/tm/04/a10/) the distinction is made between the graphical methods that are used to simply describe the data (and these graphics shown in multiPlotDataOverview) as distinct from graphical methods that are used to describe the WRTDS model of the system (such as the graphs in fluxBiasMulti). If the rResid = TRUE option is used with multiPlotDataOverview, or other graphical functions that normally don’t depend on the WRTDS model, they now become a hybrid, because they are using the WRTDS model to generate the random values used in the graphs. Take a graph such as plotConcQ. With rResid = FALSE it is a pure representation of the data. No assumptions are being made. But, when rResid = TRUE, it is now a representation of the data which is partly based on an assumption that the fitted WRTDS model is indeed a correct model. The fitted WRTDS model is partly determining the placement of the random values that are less than the reporting limit. If the analyst wants a “pure” representation of the data without any assumed model for graphing, then the rResid option should be set to FALSE.

+

Also, with figures such as shown in fluxBiasMulti, there is a kind of circularity in the logic. The circularity is this: we are using the graphs to assess the adequacy of the model fit, but we are using the model to estimate some of the observations. This circularity is not a fatal flaw to the approach, but is a reality that the user should consider. On balance, the authors think that using rResid = TRUE in all of the plots to which it applies is a beneficial approach because it enhances the ability of the analyst to interpret the figures. But, we reiterate here, the choice of using the random approach or not has no bearing whatsoever on the quantitative outputs that the WRTDS method in the EGRET package produces.

diff --git a/vignettes/rResid.Rmd b/vignettes/rResid.Rmd index 728c71a5..7e76d81a 100644 --- a/vignettes/rResid.Rmd +++ b/vignettes/rResid.Rmd @@ -15,25 +15,29 @@ vignette: > \usepackage[utf8]{inputenc} --- + +```{r,echo=FALSE} +library(EGRET) +eList <- Arkansas_eList +``` + # Introduction When censored data are present in a water quality data set, depicting them in any type of scatter plot is a challenge. This applies to the problem of plotting the actual data values or plotting the residuals. What follows is an example that shows the problem. -```{r} -library(EGRET) -eList <- Arkansas_eList +```{r, echo=FALSE} plotConcTime(eList) ``` The solid vertical lines, which represent the range of values that a given censored value could be, is very distracting in terms of getting a picture of the overall behavior of the data. We can see this even more if we try to look at the relationship of concentration to discharge. -```{r} +```{r, echo=FALSE} plotConcQ(eList,qUnit=4) ``` It is difficult to see the relationship between ammonia concentration and discharge. If we look at residuals from a fitted WRTDS model and discharge we also find it difficult to see if the pattern looks reasonable (a horizontal cloud of points centered on the zero residual line) or if there is some substantial curvature to the relationship. -```{r} +```{r, echo=FALSE} plotResidQ(eList,qUnit=4) ``` @@ -41,11 +45,36 @@ Here again, the plot is not very informative. What solutions might exist to res # The concept of randomized estimates of censored values -If we think about a fitted WRTDS model, what it is telling us is this: For some particular date and some particular discharge the model estimates a unique conditional probability distribution of concentration. That distribution is always a log normal distribution (that means that the log of concentration is normal) and that distribution of the log concentration has a particular mean (which in the EGRET package is always called `yHat`) and a particular standard deviation (which is always called `SE`). The `Daily` data frame has columns with these names, showing their value for each day of the record. However, for days on which there was a sample taken but that sample was reported as a "less than" value, we also know that there is an upper bound to what the concentration could have been. It is the reporting limit. *As an aside, all of the following discussion would hold for an interval censored value. EGRET allows for these but they are not common. For simplicity of the discussion we will only consider left censored data, but the code being used here considers interval censored data as well.* +If we think about a fitted WRTDS model, what it is telling us is this: For some particular date and some particular discharge the model estimates a unique conditional probability distribution of concentration. That distribution is always a log normal distribution (that means that the log of concentration is normal) and that distribution of the log concentration has a particular mean (which in the EGRET package is always called **yHat**) and a particular standard deviation (which is always called **SE**). The **Daily** data frame has columns with these names, showing their value for each day of the record. However, for days on which there was a sample taken but that sample was reported as a "less than" value, we also know that there is an upper bound to what the concentration could have been. It is the reporting limit. *As an aside, all of the following discussion would hold for an interval censored value. EGRET allows for these but they are not common. For simplicity of the discussion we will only consider left censored data, but the code being used here considers interval censored data as well.* + +So, what we can say about some particular censored value is that the log of the true value has a mean of **yHat** and a standard deviation of **SE**, it is normally distributed, but is constrained to be in the part of the normal distribution which is less than the log of the reporting limit. Such a random variable is known as a truncated normal random variable. Fortunately there is an R package entirely focused on the truncated normal distribution. It is called **truncnorm**. We can use **truncnorm** to generate a random number for each of the censored values and this random number will be drawn from the lower portion of normal distribution with the correct mean, standard deviation, and upper bound. The idea is to create these random values to substitute for the censored observations. We call them the **rObserved** values (the "r" denotes that they are randomly generated observations) and they reside in an augmented version of the **Sample** data frame in a column called **Sample$rObserved**. The figures below illustrate what a truncated normal distribution density function looks like. They show the density below the censoring threshold for two different examples. The **rObserved** values are samples from these density functions. The area under each density function is equal to 1. + +```{r,echo=FALSE} +par(las = 1, tck = 0.02, xaxs = "i", yaxs = "i") +library(truncnorm) +x <- seq(-8,1,0.01) +z <- dtruncnorm(x,b=1) +x[902] <- 1 +z[902] <- 0 +area <- 0.01*sum(z) +z <- z/area +plot(x,z,xlim=c(-4,2),type="l",ylim=c(0,1.6),xlab = "log of concentration", + ylab = "probability density",main="Truncated normal density function\nfor the log of concentration, mean 0, standard deviation 1\ncensoring threshold is at +1 log units",cex.main=1.1) +polygon(x,z,density=40) +x <- seq(-8,-1,0.01) +z <- dtruncnorm(x,b=-1) +x[702] <- -1 +z[702] <- 0 +area <- 0.01*sum(z) +z <- z/area +plot(x,z,xlim=c(-4,2),type="l",ylim=c(0,1.6),xlab = "log of concentration", + ylab = "probability density",main="Truncated normal density function\nfor the log of concentration, mean 0, standard deviation 1\ncensoring threshold is at -1 log units",cex.main=1.1) +polygon(x,z,density=40) +``` -So, what we can say about some particular censored value is that the log of the true value has a mean of `yHat` and a standard deviation of `SE`, it is normally distributed, but is constrained to be in the part of the normal distribution which is less than the log of the reporting limit. Such a random variable is known as a truncated normal random variable. Fortunately there is an R package entirely focused on the truncated normal distribution. It is called `truncnorm`. We can use `truncnorm` to generate a random number for each of the censored values and this random number will be drawn from the lower portion of normal distribution with the correct mean, standard deviation, and upper bound. The idea is to create these random values to substitute for the censored observations. We call them the `rObserved` values and they reside in an augmented version of the `Sample` data frame in a column called `Sample$rObserved`. What is most important to understand is that these randomly generated values are never used in any of the WRTDS computations that result in estimates of daily concentrations or fluxes, or annual average concentrations or fluxes, or any trends. These randomly generated values are created strictly to provide a more-easily interpreted set of diagonstic graphics. Here is an example of a plot of concentration versus time using this approach. +What is most important to understand is that these randomly generated values are never used in any of the WRTDS computations that result in estimates of daily concentrations or fluxes, or annual average concentrations or fluxes, or any trends. These randomly generated values are created strictly to provide a more-easily interpreted set of diagonstic graphics. Here is an example of a plot of concentration versus time using this approach. -The solid circles are the uncensored observations and the open circles are these `rObserved` values. One thing to note is that if these values were to be generated a second time the graph would look slightly different. There is no unique set of values that will be plotted. The random numbers used will be different if we go back to the same data set and generate them again. Below we show such a figure and then redo the graphic with a second set of random numbers. The code is shown with the graphics. +The solid circles are the uncensored observations and the open circles are these **rObserved** values. One thing to note is that if these values were to be generated a second time the graph would look slightly different. There is no unique set of values that will be plotted. The random numbers used will be different if we go back to the same data set and generate them again. Below we show such a figure and then redo the graphic with a second set of random numbers. The code is shown with the graphics. ```{r} eList <- makeAugmentedSample(eList) @@ -59,11 +88,9 @@ Careful examination of these two figures reveals that the black dots are exactly # Randomized residuals -Many of the diagnostics we use in EGRET are plots of residuals (on the y-axis) versus predicted value, or time, or discharge (on the x-axis). Once we have this `rObserved` value, we can compute a randomized residual which is called `rResid` and is stored as `Sample$rResid` in EGRET. The computation is this: - -rResid = rObserved - predicted value +Many of the diagnostics we use in EGRET are plots of residuals (on the y-axis) versus predicted value, or time, or discharge (on the x-axis). Once we have this **rObserved** value, we can compute a randomized residual which is called **rResid** (again the "r" denotes that they are randomly generated residuals) and is stored as **Sample$rResid** in EGRET. The computation is this: -All three variables in this equation are in log concentration units. +rResid = rObserved - predicted value. All three variables in this equation are in log concentration units. We can look at our residuals plots in the following manner (using the second set of rObserved values computed above). @@ -74,29 +101,25 @@ plotResidQ(eList, qUnit = 4, rResid = TRUE) # Details for how to include random residuals in your computations -In order to use the `rObserved` and `rResid` in making graphs in EGRET the process is the following. +In order to use the **rObserved** and **rResid** in making graphs in EGRET the process is the following. -1. Bring in your data and create the `eList` in the usual way: `eList <- mergeReport(INFO,Daily,Sample)` -2. Run `modelEstimation` in the usual way: `eList <- modelEstimation(eList)` -3. Then augment the `Sample` data frame that is inside of the `eList` with the command: +1. Bring in your data and create the **eList** in the usual way: **eList <- mergeReport(INFO,Daily,Sample)** +2. Run **modelEstimation** in the usual way: **eList <- modelEstimation(eList)** +3. Then augment the **Sample data** frame that is inside of the eList with the command: -```{r eval=FALSE} -eList <- makeAugmentedSample(eList) -``` + **eList <- makeAugmentedSample(eList)** -4. To produce any one of the following graphics, using the `rObserved` or `rResid` values simply add the argument `rResid = TRUE` to the call to the graphical function. Note that it doesn't matter if what is being plotted is an observed value or a residual, the argument is always `rResid = TRUE`. The original option of showing the vertical lines for censored values remains available as the default for any of these functions. The call can either say `rResid = FALSE` or just not include `rResid` in the argument list and the graphs will appear without these random values and without the open circle/closed circle symbology. The functions where this approach applies are these: `plotConcPred, plotConcQ, plotConcTime, plotConcTimeDaily, plotFluxPred, plotFluxQ, plotFluxTImeDaily, plotResidPred, plotResidQ, plotResidTime`. It is also available in the two multiple plot functions: `multiPlotOverview`, and `fluxBiasMulti`. +4. To produce any one of the following graphics, using the **rObserved** or **rResid** values simply add the argument **rResid = TRUE** to the call to the graphical function. Note that it doesn't matter if what is being plotted is an observed value or a residual, the argument is always **rResid = TRUE**. The original option of showing the vertical lines for censored values remains available as the default for any of these functions. The call can either say **rResid = FALSE** or just not include **rResid** in the argument list and the graphs will appear without these random values and without the open circle/closed circle symbology. The censored values or censored residuals will be shown as the vertical lines. The functions where this approach applies are these: **plotConcPred, plotConcQ, plotConcTime, plotConcTimeDaily, plotFluxPred, plotFluxQ, plotFluxTImeDaily, plotResidPred, plotResidQ, plotResidTime**. It is also available in the two multiple plot functions: **multiPlotOverview, and fluxBiasMulti** +For example. -```{r fig.height=8, fig.width=8} +```{r, fig.height = 9, fig.width = 8} multiPlotDataOverview(eList, qUnit = 4, rResid = TRUE) -``` - -```{r fig.height=10, fig.width=8} fluxBiasMulti(eList, qUnit = 4, fluxUnit = 9, rResid = TRUE) ``` -# A final thought +# Two final thoughts -In the EGRET User Guide (http://pubs.usgs.gov/tm/04/a10/) the distinction is made between the graphical methods that are used to simply descirbe the data (and these graphics shown in `multiPlotDataOverview`) as distinct from graphical methods that are used to describe the WRTDS model of the system (such as the graphs in `fluxBiasMulti`). If the `rResid = TRUE` option is used with `multiPlotDataOverview`, or other graphical functions that normally don't depend on the WRTDS model, they now become a hybrid, because they are using the WRTDS model to generate the random values used in the graphs. Take a graph such as `plotConcQ`. With `rResid = FALSE` it is a pure representation of the data. No assumptions are being made. But, when `rResid = TRUE`, it is now a representation of the data which is partly based on an assumption that the fitted WRTDS model is indeed a correct model. The fitted WRTDS model is partly determining the placement of the random values that are less than the reporting limit. If the analyst wants a "pure" representation of the data without any assumed model for graphing, then the `rResid` option should be set to `FALSE`. +In the EGRET User Guide (http://pubs.usgs.gov/tm/04/a10/) the distinction is made between the graphical methods that are used to simply describe the data (and these graphics shown in **multiPlotDataOverview**) as distinct from graphical methods that are used to describe the WRTDS model of the system (such as the graphs in **fluxBiasMulti**). If the **rResid = TRUE** option is used with **multiPlotDataOverview**, or other graphical functions that normally don't depend on the WRTDS model, they now become a hybrid, because they are using the WRTDS model to generate the random values used in the graphs. Take a graph such as **plotConcQ**. With **rResid = FALSE** it is a pure representation of the data. No assumptions are being made. But, when **rResid = TRUE**, it is now a representation of the data which is partly based on an assumption that the fitted WRTDS model is indeed a correct model. The fitted WRTDS model is partly determining the placement of the random values that are less than the reporting limit. If the analyst wants a "pure" representation of the data without any assumed model for graphing, then the **rResid** option should be set to **FALSE**. -Also, with figures such as shown in `fluxBiasMulti`, there is a kind of circularity in the logic. The circularity is this: we are using the graphs to assesses the adequacy of the model fit, but we are using the model to estimate some of the observations. This circularity is not a fatal flaw to the approach, but is a reality that the user should consider. On balance, the authors think that using `rResid = TRUE` in all of the plots to which it applies is a beneficial approach because it enhances the ability of the analyst to interpret the figures. But, we reiterate here, the choice of using the random approach or not has no bearing whatsoever on the quantitative outputs that the WRTDS method in the EGRET package produces. +Also, with figures such as shown in **fluxBiasMulti**, there is a kind of circularity in the logic. The circularity is this: we are using the graphs to assess the adequacy of the model fit, but we are using the model to estimate some of the observations. This circularity is not a fatal flaw to the approach, but is a reality that the user should consider. On balance, the authors think that using **rResid = TRUE** in all of the plots to which it applies is a beneficial approach because it enhances the ability of the analyst to interpret the figures. But, we reiterate here, the choice of using the random approach or not has no bearing whatsoever on the quantitative outputs that the WRTDS method in the EGRET package produces. From 50f2aa43ca32a4109ee43f10127532c72edbb45c Mon Sep 17 00:00:00 2001 From: Laura DeCicco Date: Tue, 21 Jun 2016 16:39:09 -0500 Subject: [PATCH 2/4] Changed rResid to randomCensored --- R/boxConcThree.R | 5 ---- R/boxResidMonth.R | 6 ++-- R/fluxBiasEight.R | 16 +++++------ R/multiPlotDataOverview.R | 8 +++--- R/plotConcPred.R | 6 ++-- R/plotConcQ.R | 6 ++-- R/plotConcTime.R | 8 +++--- R/plotConcTimeDaily.R | 6 ++-- R/plotFluxPred.R | 6 ++-- R/plotResidPred.R | 6 ++-- R/plotResidQ.R | 6 ++-- R/plotResidTime.R | 6 ++-- inst/doc/EGRET.R | 54 +++++++++++++++++------------------ inst/doc/EGRET.pdf | Bin 777428 -> 777435 bytes inst/doc/rResid.R | 12 ++++---- inst/doc/rResid.Rmd | 18 ++++++------ inst/doc/rResid.html | 32 ++++++++++----------- man/boxResidMonth.Rd | 4 +-- man/fluxBiasMulti.Rd | 4 +-- man/multiPlotDataOverview.Rd | 6 ++-- man/plotConcPred.Rd | 5 ++-- man/plotConcQ.Rd | 4 +-- man/plotConcTime.Rd | 12 ++++---- man/plotConcTimeDaily.Rd | 6 ++-- man/plotFluxPred.Rd | 4 +-- man/plotResidPred.Rd | 6 ++-- man/plotResidQ.Rd | 4 +-- man/plotResidTime.Rd | 4 +-- vignettes/rResid.Rmd | 18 ++++++------ 29 files changed, 137 insertions(+), 141 deletions(-) diff --git a/R/boxConcThree.R b/R/boxConcThree.R index f926e190..9d06ebd9 100644 --- a/R/boxConcThree.R +++ b/R/boxConcThree.R @@ -45,11 +45,6 @@ boxConcThree<-function (eList, tinyPlot=FALSE, paStart <- 10 } - # if(rResid & !all((c("SE","yHat") %in% names(eList$Sample)))){ - # message("Pseudo only supported after running modelEstimation, defaulting to rResid=FALSE") - # rResid <- FALSE - # } - localSample <- if(paLong == 12) localSample else selectDays(localSample,paLong,paStart) localDaily <- if(paLong == 12) localDaily else selectDays(localDaily, paLong,paStart) diff --git a/R/boxResidMonth.R b/R/boxResidMonth.R index b91f923a..1e3ed5e7 100644 --- a/R/boxResidMonth.R +++ b/R/boxResidMonth.R @@ -22,7 +22,7 @@ #' @param font.main font to be used for plot main titles #' @param customPar logical defaults to FALSE. If TRUE, par() should be set by user before calling this function #' @param las numeric in {0,1,2,3}; the style of axis labels -#' @param rResid logical. Show censored residuals as randomized. +#' @param randomCensored logical. Show censored residuals as randomized. #' @param \dots arbitrary graphical parameters that will be passed to genericEGRETDotPlot function (see ?par for options) #' @keywords graphics water-quality statistics #' @seealso \code{\link[graphics]{boxplot}} @@ -36,7 +36,7 @@ #' boxResidMonth(eList) boxResidMonth<-function(eList, stdResid = FALSE, las=1, printTitle = TRUE, cex=0.8, cex.axis=1.1, cex.main=1.1, - font.main=2, tinyPlot=FALSE, customPar=FALSE,rResid=FALSE,...) { + font.main=2, tinyPlot=FALSE, customPar=FALSE,randomCensored=FALSE,...) { localINFO <- getInfo(eList) localSample <- getSample(eList) @@ -65,7 +65,7 @@ boxResidMonth<-function(eList, stdResid = FALSE, las=1, plotTitle<-if(printTitle) paste(localINFO$shortName,"\n",localINFO$paramShortName,"\nBoxplots of residuals by month") else "" - if(!rResid){ + if(!randomCensored){ resid<-log(localSample$ConcAve) - localSample$yHat } else { if(!("rResid" %in% names(localSample))){ diff --git a/R/fluxBiasEight.R b/R/fluxBiasEight.R index a3f791e6..a70031db 100644 --- a/R/fluxBiasEight.R +++ b/R/fluxBiasEight.R @@ -21,7 +21,7 @@ #' @param cex.axis magnification to be used for axis annotation relative to the current setting of cex #' @param col color of points on plot, see ?par 'Color Specification' #' @param lwd number line width -#' @param rResid logical. Show censored residuals as randomized. +#' @param randomCensored logical. Show censored residuals as randomized. #' @param \dots arbitrary graphical parameters that will be passed to genericEGRETDotPlot function (see ?par for options) #' @keywords graphics water-quality statistics #' @export @@ -40,7 +40,7 @@ #' dev.off() #' } fluxBiasMulti<-function (eList, qUnit = 2, fluxUnit = 3, moreTitle = "WRTDS", - cex = 0.7, cex.axis = 1.1,cex.main=1.1,rResid=FALSE, + cex = 0.7, cex.axis = 1.1,cex.main=1.1,randomCensored=FALSE, col="black", lwd=1,...){ localINFO <- getInfo(eList) @@ -72,27 +72,27 @@ fluxBiasMulti<-function (eList, qUnit = 2, fluxUnit = 3, moreTitle = "WRTDS", par(oma = c(0, 10, 4, 10),mfrow=c(4,2)) plotResidPred(eList, stdResid = FALSE, tinyPlot=TRUE, printTitle = FALSE,cex=cex, - cex.axis = cex.axis, col=col,lwd=lwd,rResid=rResid,...) + cex.axis = cex.axis, col=col,lwd=lwd,randomCensored=randomCensored,...) plotResidQ(eList, qUnit, tinyPlot = TRUE, printTitle = FALSE,cex=cex, - cex.axis = cex.axis, col=col,lwd=lwd,rResid=rResid,...) + cex.axis = cex.axis, col=col,lwd=lwd,randomCensored=randomCensored,...) plotResidTime(eList, printTitle = FALSE, tinyPlot=TRUE,cex=cex, - cex.axis = cex.axis, col=col,rResid=rResid,lwd=lwd,...) + cex.axis = cex.axis, col=col,randomCensored=randomCensored,lwd=lwd,...) boxResidMonth(eList, printTitle = FALSE, tinyPlot=TRUE,cex=cex, - cex.axis = cex.axis,lwd=lwd,rResid=rResid,...) + cex.axis = cex.axis,lwd=lwd,randomCensored=randomCensored,...) boxConcThree(eList, localINFO = localINFO, printTitle=FALSE, tinyPlot=TRUE,cex=cex, cex.axis = cex.axis, lwd=lwd,...) plotConcPred(eList, printTitle=FALSE, tinyPlot=TRUE,cex=cex, - cex.axis = cex.axis, col=col,lwd=lwd,rResid = rResid,...) + cex.axis = cex.axis, col=col,lwd=lwd,randomCensored = randomCensored,...) boxQTwice(eList, printTitle = FALSE, qUnit = qUnit,tinyPlot=TRUE,cex=cex, cex.axis = cex.axis, lwd=lwd,...) plotFluxPred(eList, fluxUnit, tinyPlot = TRUE, printTitle = FALSE,cex=cex, - cex.axis = cex.axis, col=col,lwd=lwd,rResid=rResid,...) + cex.axis = cex.axis, col=col,lwd=lwd,randomCensored=randomCensored,...) fluxBias <- fluxBiasStat(localSample) fB <- as.numeric(fluxBias[3]) fB <- format(fB, digits = 3) diff --git a/R/multiPlotDataOverview.R b/R/multiPlotDataOverview.R index 2c898f13..67c6b7b2 100644 --- a/R/multiPlotDataOverview.R +++ b/R/multiPlotDataOverview.R @@ -13,7 +13,7 @@ #' @param cex.main magnification to be used for main titles relative to the current setting of cex #' @param logScaleConc logical if TRUE y in concentration graphs plotted in log axis. Default is TRUE. #' @param logScaleQ logical if TRUE y in streamflow graphs plotted in log axis. Default is TRUE. -#' @param rResid logical. Show censored values as randomized. +#' @param randomCensored logical. Show censored values as randomized. #' @keywords graphics water-quality statistics #' @seealso \code{\link{plotConcQ}}, \code{\link{boxConcMonth}}, \code{\link{plotConcTime}}, \code{\link{boxQTwice}} #' @export @@ -24,7 +24,7 @@ #' # Graphs consisting of Jun-Aug #' eList <- setPA(eList, paStart=6,paLong=3) #' multiPlotDataOverview(eList, qUnit=1) -multiPlotDataOverview<-function (eList, qUnit = 2,cex.main=1.2,rResid=FALSE, +multiPlotDataOverview<-function (eList, qUnit = 2,cex.main=1.2,randomCensored=FALSE, logScaleConc=TRUE, logScaleQ=TRUE){ localINFO <- getInfo(eList) @@ -41,9 +41,9 @@ multiPlotDataOverview<-function (eList, qUnit = 2,cex.main=1.2,rResid=FALSE, par(mfcol=c(2,2),oma=c(0,2.4,4.5,2.4),tcl=0.5) plotConcQ(eList, qUnit = qUnit, tinyPlot = TRUE, printTitle = FALSE, - rmSciX=TRUE,logScale=logScaleConc,rResid=rResid) + rmSciX=TRUE,logScale=logScaleConc,randomCensored=randomCensored) boxConcMonth(eList, printTitle = FALSE, tinyPlot=TRUE,logScale=logScaleConc) - plotConcTime(eList, printTitle = FALSE, tinyPlot = TRUE,logScale=logScaleConc,rResid=rResid) + plotConcTime(eList, printTitle = FALSE, tinyPlot = TRUE,logScale=logScaleConc,randomCensored=randomCensored) boxQTwice(eList, printTitle = FALSE, qUnit = qUnit, tinyPlot=TRUE,logScale=logScaleQ) title<-paste(localINFO$shortName,"\n",localINFO$paramShortName) diff --git a/R/plotConcPred.R b/R/plotConcPred.R index e3299694..c613b75b 100644 --- a/R/plotConcPred.R +++ b/R/plotConcPred.R @@ -18,7 +18,7 @@ #' (for example, adjusting margins with par(mar=c(5,5,5,5))). If customPar FALSE, EGRET chooses the best margins depending on tinyPlot. #' @param col color of points on plot, see ?par 'Color Specification' #' @param lwd number line width -#' @param rResid logical. Show censored values as randomized. +#' @param randomCensored logical. Show censored values as randomized. #' @param \dots arbitrary graphical parameters that will be passed to genericEGRETDotPlot function (see ?par for options) #' @keywords graphics water-quality statistics #' @seealso \code{\link{selectDays}}, \code{\link{genericEGRETDotPlot}} @@ -33,7 +33,7 @@ #' plotConcPred(eList) plotConcPred<-function(eList, concMax = NA, logScale=FALSE, printTitle = TRUE,tinyPlot=FALSE,cex=0.8, cex.axis=1.1, - cex.main=1.1, customPar=FALSE,col="black",lwd=1, rResid = FALSE,...){ + cex.main=1.1, customPar=FALSE,col="black",lwd=1, randomCensored = FALSE,...){ localINFO <- getInfo(eList) localSample <- getSample(eList) @@ -80,7 +80,7 @@ plotConcPred<-function(eList, concMax = NA, logScale=FALSE, xInfo <- generalAxis(x=x, minVal=minXLow, maxVal=concMax, tinyPlot=tinyPlot,logScale=logScale) - if(rResid){ + if(randomCensored){ if(!("rObserved" %in% names(localSample))){ eList <- makeAugmentedSample(eList) localSample <- eList$Sample diff --git a/R/plotConcQ.R b/R/plotConcQ.R index d8095984..e4896140 100644 --- a/R/plotConcQ.R +++ b/R/plotConcQ.R @@ -24,7 +24,7 @@ #' (for example, adjusting margins with par(mar=c(5,5,5,5))). If customPar FALSE, EGRET chooses the best margins depending on tinyPlot. #' @param col color of points on plot, see ?par 'Color Specification' #' @param lwd number line width -#' @param rResid logical. Show censored values as randomized. +#' @param randomCensored logical. Show censored values as randomized. #' @param \dots arbitrary graphical parameters that will be passed to genericEGRETDotPlot function (see ?par for options) #' @keywords graphics water-quality statistics #' @export @@ -37,7 +37,7 @@ #' # Graphs consisting of Jun-Aug #' eList <- setPA(eList, paStart=6,paLong=3) #' plotConcQ(eList) -plotConcQ<-function(eList, qUnit = 2, tinyPlot = FALSE, logScale=FALSE,rResid=FALSE, +plotConcQ<-function(eList, qUnit = 2, tinyPlot = FALSE, logScale=FALSE,randomCensored=FALSE, concMax = NA, concMin =NA, printTitle = TRUE, cex=0.8, cex.axis=1.1,cex.main=1.1, rmSciX=FALSE,rmSciY=FALSE, customPar=FALSE,col="black",lwd=1,...){ localINFO <- getInfo(eList) @@ -85,7 +85,7 @@ plotConcQ<-function(eList, qUnit = 2, tinyPlot = FALSE, logScale=FALSE,rResid=FA xInfo <- generalAxis(x=x, maxVal=NA, minVal=NA, logScale=TRUE, tinyPlot=tinyPlot) - if(!rResid){ + if(!randomCensored){ yLow<-localSample$ConcLow yHigh<-localSample$ConcHigh diff --git a/R/plotConcTime.R b/R/plotConcTime.R index fddff605..18435363 100644 --- a/R/plotConcTime.R +++ b/R/plotConcTime.R @@ -29,7 +29,7 @@ #' (for example, adjusting margins with par(mar=c(5,5,5,5))). If customPar FALSE, EGRET chooses the best margins depending on tinyPlot. #' @param col color of points on plot, see ?par 'Color Specification' #' @param lwd number line width. -#' @param rResid logical. Show censored values as randomized. +#' @param randomCensored logical. Show censored values as randomized. #' @param \dots arbitrary functions sent to the generic plotting function. See ?par for details on possible parameters. #' @keywords graphics water-quality statistics #' @export @@ -42,9 +42,9 @@ #' eList <- setPA(eList, paStart=6,paLong=3) #' plotConcTime(eList, qUnit = 1, qLower = 100, qUpper = 10000) #' plotConcTime(eList, logScale=TRUE) -#' plotConcTime(eList, qUnit = 1, qLower = 100, qUpper = 10000, rResid = TRUE) +#' plotConcTime(eList, qUnit = 1, qLower = 100, qUpper = 10000, randomCensored = TRUE) plotConcTime<-function(eList, qUnit = 2, - qLower = NA, qUpper = NA, rResid=FALSE, + qLower = NA, qUpper = NA, randomCensored=FALSE, tinyPlot = FALSE, concMax = NA, concMin = NA, printTitle = TRUE,logScale=FALSE, cex=0.8, cex.axis=1.1,cex.main=1.1, customPar=FALSE,col="black",lwd=1,...){ @@ -108,7 +108,7 @@ plotConcTime<-function(eList, qUnit = 2, plotTitle<-if(printTitle) paste(localINFO$shortName,"\n",localINFO$paramShortName,"\n",title3,sep="") else "" - if(!rResid){ + if(!randomCensored){ subSample<-subSample[subSample$Q>qLowerBound & subSample$Qy+lZ1)dCG|DzVCC+Ju0J8R%jSs${vwD8^$q~@ zgGWUZugo;GeSP|l{gKi^s;F1%uWk(S?lj!(XUx$NA1pmdUXm+jl0wy-L^ zsAK=U{;hiYD_;z4v6puoymzjR$)~+`u31y^I(MDizGQfbbERM7Jc)I{1#I* zeJPVs`moWOYb`%&FzWq6w*TpF<+0`42IJsqo9~Ub^5BPc4Eyu><=IZZ)-!5cST@(j z(R}~hy#-OzVj2b{-*;$!WK7#N{?6$SCd~qobHr`o2#eR$CQau5|Hgry&Gi zzv6g5gZrwV=v_N+>9%y6uplx0oo?>A6US!9XUdP9n-{#)P?eYQwe04hkHtSe#osx7 zv~<+w@okQdIghnh<^56Cc4#TP{A*18g5pabGu~hKnLht`j#kC|Rn&|BcCGBfl1b&G z$7D_oNIAKENvF)P9<8^u<#j43_py$+#l|~x4P1XdebAjx)_LzX{r!L-Gja2{a}84A z>~Y^`9?af-qI;b#?Q+(d_`d+%!3@v8?lN2QD9bZHD?06Ne8aq_<4kwNbUxymm(%#} z)ZJNn36TLykFMxw&~k9WP7qhTQ~$K0%@;#|oZ%CezlmBBZF6|8QS0be%B}0$zhjNR z-g16gc4$Q-mzM4R`E~sw(|rASW0T(>3sYLhCu^Nrm$j>G;0ShsPb<|)Y>W!Wi^~7aK+xoi)PTT0SwX0L!1MkqlVfWUbZShKy=M+1v%j?728rbsoyDC_s za`exEFHB~|PI-E&bIZ3aM=ogl{NRZD)3XMR&d!esyE(M~&Z#3awKVtHSmd~mg+DFl zxCU-}zvi=@PZAD?avgB*#@b2r)uYT)wx<~KB z8`D0|eGtWYId{ow(zqZ?yz=Mb`EfaA*XEl29NT54V;AcR*5XALsu0~<=lXZNxU}8J zmILOrG^-d=eV;Renm%Eihx>L-b@{OU^p2T_{qhgCPx@AJ-#6s&)SeDO=9>F75|xRl z3a@xRGrrv2(<^c97Wvvck^RRk%dRWodK$m}CPEy}~df|1Z&CG6XMunR@#r{v((oa(4bSZz_ zEb(GNNl}gY@4duaeg2|ligmO)dg!-7@}$7mBixl+)#|cwSwBx6`rgxOR>Z*4o3f*Q zYOI&n$)4F!tcba`V8n$>zM~^}m*!l4i_f3s+ne3E5`5%Ns$Q(Yi^Ns^EgBC4+ zcO_i+ONiMu!`0PFu%N!ftL~%5_&)udn4T~rmCu}h>eRx?(Mw0WnhYH^N-=9_i%6k3 z!E03ZHmib$oot%lf6{Hz_WB2hH@+Kp^+EeN`)0=nWOp$(G}8$>9}1>^eoQg2x5EeL z9?!~+qkO(4=RM7ScfNd>AZy*SCLiFQMwSQMG5Oqa#`1xWH-rSrGeQE+=1=tP9%a(ng?>o6`;$pR%<;NH4?Heq4v%lnaZjjxpB#L z?-M>g(>Gcf*gQF~xt2~#eYs4gU^qrCxIff3uKJa{TdFm%Yt_%%F9p2pm+a-OWo24V zGn^IN9%?uHh0bysgwKn$b~OJlgz|7L9p$e<99we@{Hj=E=5?(^D-DK1S`$P3Re`bT zo7Q#H|HQc3T6>6fC5EbE6}^gGD+EJtt^Ubt`)hY?_aE;~r)rNg|Bv_Y)@!$}+QDz; zpMO~Ue&_!fpZHLFe&7FiztvQyLnCduG|c^?B|qZn#wSTl_;cqRrEWtc{*ywCKQ}qv zc&7Dm^UmXQ--qm4^L-C`%iPX{a*h4X1#q|C*J?qcCJIyFO#@7je>J$tBPL6 z6_!L7mIs(UDQ%zJ{K(;nZ&!Cn?`;y8Y-e*-`;CQV{THr*$zC5i{mdKMd2xb4;?JWa zyGDJBJ^A;@41WHo=gN;+E5}_H&WlVAeY?kelTW7ejIsar=07(LNF3#16yp%S z)I)N(Uun~l@Zq7p@unGD#(1x9Gwed2M{UP#EP3T?l~`C5-~Hax@>HvMS<@fxW@naX zu01g?%FFMx*UUTagIrdY@0r8yH1U|V+Qcit>0!pH_|=nJzO!EHU%aPX^7Fj#pcBgm9x7CZj#zl}^ujsL85>`QtZ|4uHRe?OQa&>{D7KAN z)0fX27M4l26dl~O_vL_f(UZ%EJ}P;2*JStXXLDQ?b+%3NaotqBbmSGsFCWTBHG7x9 zzc)!LF6}hXFm`NMulxo}7kt<5b)A}jGy17psQUcjEot?U$!d0Jl*$Z|LAoqrCaYd2)IAildim(T)(+ZxpAL6J)4&QlkBwCrgQE? zN3+aWcl|mm_GAxgQh&_9>dq;$*avMs1lBq4!}24}9e&ce*eP~H_R6A^g#(vNe7hj| zL)^+vmh;6#uiww)x`oa= z-8-i5&bzf}aH4Nx^Y@Q7%ffU+hS>h&q1~(1;1{!R8tdFp#7@xd@xG(|%d7>%q5>CK zB-T3>D&KnYn`?O8mbyA4ho(v9rWt0k1s|?#y>qrpNpOT^qus?vZk&}pSYhKedDp=g zrAex^tHEhk-Ui*v?eM=@yWa9!OG@wg{Lf?>)BDMvPf837@09tU*^;wHP%er4v-!}m zw=E5Zzt}T!H z>Yg^5mFqUuX-UMX9oeIX&$t=obAP?;`r~_7JV%W?lXvjG?VV%KwyZvSaE5E*l9XJ2 zM_6F{`)kMhShf8c5!YSs!s62AO z+9oACCRZRF86%S!dej2N_THINzTr`qlTx-zWLAEbl&-4!uRX;%6cn% zK1xgTNZeES=zGMT(9L&Nn+}vVe*E?9qgHvFz3#qt&ODVj>S>};&y9-}cMWD*q`&gp z=k@V+YEsh7HrD*?af?b5j7B^5G;h&=^ZVCk_09h3>eaeI zitd;pO>PzHJu!%=b82a?SuSZ&T~mzaT|4=*ebzq}f|!f}3ryZi7IYuh@em(6;eEK1Ec zbNrc@_;m9|SKA9O^{-|3Rka>;Z@oeKkhBQd3FlM(UA{>7_SKq~+o#yGi(KDhQojJ* zFG9}R#s~lBc~XDE&yGE#`-k|Cd$G6yzxV8q^0zkzv@7GwCvDnva*}b$hw+bJD%-q# zB3Tmc=$L-CiE3x-IU8Q8%v#)z&imuRqmY00w){sG{;Zj6d!K}neMi6lSFhW=vbT-h zCNy4Pl5f7FdFF-AUzqSmJxc=|ElNh5uREaNe_8|DIPM9|xEnn#^UCMgmLJNtm1TYF zH_tvXo!{qa-g6GWG{InagvY@7qn`(UO8EG7_hb zb#~twxoOT^e}mQo*0y@|p>fQ3gBPB!R@~~d+0EQx+tohKX@)Khe*W1s=wd(pb?008 zE$H*rX6^`iE1j@4qX!yj4QTZ(dCh)RhmdZ4J-P=pT0cH6M8&V}@+vF74C88p23)WY27+WaPQj~XtpHJ|GW-H#lhPg=u};{Cl}L!OW=9lr!O<9$nYHspXr*k1`9N9PXX!v}@|TkvVlPP4MSt z4V>Ur@9$-N`l9-df3x@E<8=6;t)jg)w8)&#Di-E93-}&?cICS7IXU0Wj+D&o_E*R7 z=)L0xPW7MTXRZ9BcTTcb@06BZp871xnUhrSkMssx{4%beKF)bw_`dE@+5Q$|=IX3_ zWZd55j9t^^u3>%3<~`#Y3|EwOi0pcC#m4JB#yng6UHxH6d|9MfM9gYaq08scyN#?? zmp?3aJTs6Pde%p}ByXk96^9N+vl{9>=rGRAqLat6*5&nVN4C?8hMo3c`+wN%NP z${DG(BFwEr3lS}BRz1D2vmFftFI%0VMht5$Q`?9d^vdLU?1VGeeJ2 zs;w1_Fv?CRLR+q`@S0jq+xYUtyYFQN6TEgcI{NX;50!D}nD3oqrugN_)dj()_PC9= znDg-1;v{C?*HQ~-%L#VbXKwCd&e#=)TDcoLdzd@BA5Kj4{@13t_0ONDPoG%)AvPk< zHgbEwgSitdB42v-+oYQF+AHJ9x9qd;C(rA-_T8R+KX=<6ujlUJ<8bbzJD=rJQ0koM z+-v-i+(FkDKRnxYh*RgNJwH!w-V=5&JHmHT=chNy(oF`M4cj+;j83SmZ_`)4;x;QqucO%s0u-*MO^{;uP)|s+KHqmQ(*BSlR-0l>xa?I^_`tC!%uAP+j z0rIlz5dK7U7p34=5C&zv+!&1;D^Z% z@4Gbl_RmMpTOY&fOrBbBWA{j-XZd=8Hn*1^)rt7iZi=P#3`JV1Ug^ja%S#Jm+}zeY znmD+)-EYE_iP1&|<`coQKm7bt%-=jtslhvSd-u^-CHpWlK)Iy0o!*yUx1D zeJzf?Ty(bCqO(R9&n9M%F|~K)wL1;)=vurZz3#|6kG&UvpAyrqQ{lYvS4>0y?LAfh znQDD~zXslgDNXA77H`V{quJtSo-1{ruUE-e&`w*4rf65vtbnb?#F>Jgu>B z;Q8Rk1K#_+*uJS`W$Tp5=f)&VDnB&EZrLIm!@zdV+lvzBtczaB&TH@Jo}Q99H#Q;u zwpO}ny4%qvZ^O?&yz=+nppic$J2u@82|3#NYp-8y&xIo0_lFqV*|<68U7a zPTO?)pU#ePU*v0*--p>2F=+6Pj?SJY@iOnwtLw~Xd%GGrwGY$l$;J-xyd7XPv1G|^ zu8z~5!^=DV^=74!m$}6Bz~g1ToEP~vYcRslAkJk8`^Tld{nLB4I^KTBmhF!u>h?PV z94^~7zIE%g;b7jw@%txBqwzZ*z?Z|Ja&Gl zwDTLO{*GZ$W{Eo19yjYO=&kqazE}R?T%DE^ zZ)%M+@93fH?vyiLJz!iTvxxq={L0lHcJ-$FPHnzwpO=1%=o^ctH+s5z+`VB*aW{Kt z^>xYfTrocA@wUaEbJ?z6r}edWDm<*;=&}2YNJm=(pYun0d@LJ0*e}X|LnD7z`wPY$ ze9~nNX1x8`?cqXp=o7UIl-A?+2HV4yM*)Q6k-F#o!rMDmU z^TT7V&l`VHz9K1V{+Qx;<6hm`j{eiN+ma=p${yaD;!)ByFLHf-=FZQ^3q2Fw&X4|F z@NLnj1z-1%*7|lY%3;OfqRYbkh9?#*dXvhv)4pwUV#PlP&3i{QUw%V(|1ir8Civ7e z<9VE!g_*y>#Z>+AIk9y|XPCZB4GzCoe(g|W*Pi-OeCmX|=h9XvdYsR1RX53Z=#U1z z>xK*r_A4{(=PFws-F)l`^$<_N&89eYm`8$S(bkMlZNEH>_AM~7pY&!*mi@?_(RX+& zOW$6Lr`e9zTauLie3;LP%r2+xpC+h6%KcUsIQhA{4QbY6!<+3EPx8&0HqdeWuqA7a z)&1U0ts|XQ4q*AKnG&Jh1*LY+4Yz!4>i2v*cAC|-&sw3nC;PTvRB~wIj$zxj^q6md ztA3wp6PoM{ns6fU{^0<%C#zF#4n1UPVf1bT*ZAjJlV^X&@3jhleaEj&lU8T__?#aR*KC&r zo3{THb-k!VlHbuQ7rE@0{+{LYJP#ax>392H`Ty?zar^l9FOe@N$DLQkwYsG&4S4pq z-pIxWmT$;z_~k-fq21)&XVh9lFC@&%T)T&zcR=giu2}Q&CR2+Qw{JB~GH>C3)I8t+ z{FF<3R&+eIcf_=`-g;Z(RY#|}g)0Opf9y!t6(5C+8iD8{ll_J>a8hkLVX22MZHI z`rEZXI$*g~NW&C;oek!jOu|=0wtpG=S4+dAc`=6$eCk`I7;;JPQBsr78Uo%$Vq3HTpBt_*E+Ml+v}hMLw)K^Xk!=a*@!D~8QRwOg-=L=(cTW7-ad6$b}wsu z!dCSV*0yWe#jlpj;}z_7v5Y<(-6@`iUu``QW(d+U5_B9mKlQne|gVytDNTgHYpvFWs> zpw6i&%=}AQznpt$*dSWYUU!n-I|?I$r0mKmH40Vvx?U%8c7@ zXIOqv=cUfMmRm2zf5j$h1?*E05ym_^5veExF8ZY}NCq?Q;w5i0Xl^>O{0xvGTtVng!OqFrspY zidyR{G;28DtGnH%5<^*=Up9${9H^_;uy<{Q>R%O?dAYCNAC)OsnPv(CW9U@9W|c|U z-(x)3tk=V#dQ?q9DJ!ljdD&gP6;*R+61`O?Z{m%N_{=gryW5ph{0^bn&3cRebB3~# zdOqDNXQ-KySIE?wP26JjKlS`Cj9{w;`uMAlKZ9he1(N3azp&+U!F7qGQZn&t))L8% zW=d13_|J5vse$-ZvaC|d2&PLVcFkDol|rTxW{a<6)GJoa3aLvawk9$)^$wawFyt{z zPQQ|>g|Xr*j(R0$IN_-HN=3br$z&1BBz8^d`2>e`5@TW9GKnR<7^CK-!Y+(!CM9QN zREQW0#VM(MFoM=5OzFP-w^J$QwcG9gENMeU#AP`!|@z$*o{5faJ#wmrIpx>@0mE+cefA>nJM$K#4 z?_a4D3gMq%P@kblEK@QNwKaIBrr#-WGqcxVHR{r294u+7I7~qiS;i_WqevoS zRFJi`7$&EPER(AtY~nkZdxA$~a;YpL1g|)XS~8hjCb)%2?6g$Ex)4c6ZB8v550P}z zk_k`52P5c);=u{s#D`oMB0gAQvG`C3do{r88ep;b;Dk2oD&F-LA4*}m_+W)_JOpN? zq=|5SouoMy%N5~Lp}~4QYM6gOQ4zGCFYaSEDhf-LFhY+@1(p0z+xSiADLoVcs53r7a zhe}|@heGhiL-1S<1`8L!U~(Z&kaQHrheO}K5+o+{RdFnNRk%LTn;^Gb4trV~hO3e* z3HvFP%b`Dm=c@JjjB4gw#j#ya6%2GPQz=QMB*Ih#38hLIwn$8|b*2ZIDY`N;qXq;0PtB{NYJkVMQ z-(Dfyr?*Qibign|uJ~IiG}wU$C$M<%)w5N2x;Q_}63iim&q;DGWK&Va2*2(eW1(D~g_co&GV$&F_P$&r*#xbDG$O+YO1u2+{6fY^f zOd7^WV96SjnhSFd`RGsYp$R4YC|DA#}G?O>in;Vs)z)1x4|f5uDG9K$Nve zL?x+{rHo8X4hBrg5M!}HHAQ|n`!&Qy!{`j$0G5(@7aQn7H-^r5H@ARVf2X z6Q=_VAu2hRtKrb_8W{#6Z6+ColoBUKA!kX(WE2dkyQK^qWs)+CLZu|GNy;-E=uhsD zVJJHe4Nt`djGCnvQ7h!+d{AYi7(@NCDo9C%4z0A~Qbq}-M{Y)~tW*|Ncr3_Il9gpx z6?si8L)vVql$F+0wY-KISXeOP-cdtEWK3)ZH%##Yt6(S{&nh@lnZY_M#1o7Est>|L z`%ouQNEOt=xOfVtP_gt_HF*%xs;Nj`VAZf;*o;X3>fNYmvEyKcs9mdnb7Iho+Mt^9 zW?(vNNrc9RatO3Y(c$E*ikuK?o+2xU_<+7T6(yc<7AQvNI9^TZbPg^ANiGh}pwwm{ z+tj8R87U=F4zUWsn9#O>ioAgHSCJn!M@B8IVkw2hDQFW77l1x73=BCRa#gi#Nd^6p z*d)UwSYAZN1h*<=GzTkS4oE5}AfLpQz#9}f4SzKbjC1g$k&r5sxOf5)s)uxifJfSL zDPS^6FDMZEQ=3sjB~`4!SCbzZK1y_yqBBZuMnTCtB#aVSm?p~fDrLmsMEXPNQPxhO zR8+`kO>*IIGJ=t$5LimWS_#_J9xDS4U`s0WUzH935Q0Y`D#`E)cvQ#*7-f2qf`Dw3^C3T5b&y~F zAMQk14b=FA2)?OU0C=E3IiXTcQUTdXH7O+!Sg2NFEH6{4$oLG1sB034NSFTxCL>*K znM@^9lcbiZs5}l12N5IuKjdfB=wK6VMD>-(WlD+(p(&{QL?azdL7am;v0Pm0dNwHfRNO@73M49TqU9BXRP3i*R= zlGg-cPoM(SExBt1m83Zq?HcJBBVj`;DVCLz<{fyv;x`FFSbdU&2v^9Hz@H~WF(k7n zLj)LnXmUPQrg6=1{wne#eMgB7%nKPG$e_2#1QIfq44Onm!qu2h2{(s=)%G<+`itpK zN=w1FC6fujLk@(Pki!JzHQ~IZJ1-Mmcao_&DZ~E;$WkXkWsuT(95OegY;&;NRFlajXmsc^Vm`oBJb7Zc-=rriLxM(HYu6&rK+lj9!X^I=OrfU7 zDk~3C<^N$q6&WUCM4Cg-62_Qi7?G!xuB;Bo$!sDQl;KNg%I8 zbj?VyLrs99f>Ny}V+LFaJTAijLqP%gNaDnBH7K2q0vMU|5hFWNEd$2>2F6wSn3et? z5NpNdq9Oqdxo=`KYSak`6OD>Tr8re7AHZmyK?Q}pP97%~+uOzUY7Ul@(GLs^b)Rx6 zEG4l)z>rxY=Yysp#2EQI7*N9h0}#;#HCRq&!+=N0HtJ1k@^tXvV4sNjP#&(XZv689 zfQS?Wi6~W(*$$AE%Ije)RZf#5&QTNM2^i9m#2UO<8Yh7=C^Pa(P<|oREhZE@=4xM^ zcn(5jQSMWdHUbqf1c8MAhc*L_bT+Vvid45^IgQlwfXPWyE-pgq09ZAo9jZ31@&5pX z60QPPQRbo2s&W)ZNhbo88p^T2SgK|R@dCQl;rDoIwjdZal<+MIIh0OM)qVxhW74o1jEDJeL{q-zF7XGl*Zl6J;QUqOM^3h~3!w%K2ALost5#BX%%U%eBU>QLjP#11e@nJG{P&|Gn~ z|0jnYq!vNNgd@p{94MuHVDV{~T3SQVXhc>4MOf>=pg(%ls8kxR1lb|t#E7gZvZA6x zrfxA-MX6gQ3@S-MoDKSiTtr!29Q>6BLI;A9ymv8yN`}^`KVoTOKIm17ft8|{qrt24 z{}8>yp(h8!ua&ZF@Rms@6HN~Iq>h*mHA&)%d6Z!hS(7SQ#d5!aAyG?Bh?Euu zqs_L8&J&`CL&2#C|4$CyqLd06o5;ML*Z@MQWi=$3JTZ8X(vBc>m3FN3{~+14%7nBN z}B{AUwKOsGb~LCh{EwS;d-DEzto3qNYY#O2DWhmI$M+ zQ>;^x`9I`_5SWt-)A2j-khh`m(8@`}$n-15BC}D6_sjo7zK9`nIe@Vwrvt{3wLHL3 z2&|k?lRE?qLPo);2n+HrQxhSPSXD=p6)_l{RYK+Bm{`fAe~<0a2o#)J`-ezL_`IQJc$$vi6SckB9i>@ zpvfQv9v4VWtRdrBDvkppl4&F&$|&?AFvuhxZHO2-lG-GmTm`okV@3O4>1>GE*Q&tq zm+nWi7)plJA}Uf7g8Y=3#aJreFE*dmJe;^bIxNt>M~IHtU_mV+2iK5PfcO0y7*!d^ zd@7288)wBU~|?OKLs9YC9?-{b6xvu|xKxcBZNV533Lh!^I>Gk=V;pRbZ;* zBLx!4Bp)H|#gkqMnptS3LFS8$1(5bdO}X}F;FM5_M={?_=>Qd4dNdMN?f)SMjaK8@ z69Psx9EdPV2dE%pBx52)LuN`K*KBPQ4f>0vB+|XY>AB?fp^S97khM`z_l{yqrIoGn|A43_Ef(4x(M(%= zVZf+zEnrl89asTA5;+zwLv70NsL-KBk+?g=lvJ)yj5EpnFkmbxC4g}x!-=H_k}~K! z#)T689|oh2Lhc&=DMf0u1X8C%LcI3r;QYZ)q{gB*tP)n~|6zwz4M6PcB?Ae-;8znH z1dP;#XxF8JRKTbHgixeKpk`#N8IfKqP@)SCe@IbljoQio#Qn@qyTxC&DW1f29QkkV?IXFw(mP zjFu7>PAkbDVkkgv1_ja@Fpd9*{3(@o2Mk|@5GO%0olLWc{p%ECa>z20Qi2j4>1neZ zUuB7^`~WejN+}wY-U!VlWzk1U-!VE;NWFjx82ObIO40DTa5$0vRpdt%g>)i-NHrIT zFtRuW7_xl?wUEfE)bmwqD3I2xL78cy|A($hWCW;B9OG`fKf^kdWDMS zkVl09AvcIlc0x#j1&5^?uR`F&BrM4iU|A}P2aKv|z;UJ8k^w{5SG|aMser-wP%xgU zEd5oHMog{}k_%JIl36=^0s_-dPKeI|NGd>=Rx4PIjuby1(FiWiU$roj$s~)Ru7e&i zVi7o`WGI3!9gw*r@sN=uHv@L76$y&;$6$N`MGS^nC1W5CCV-R+I2=@y24L`4$oasm zN$LbAfpnl08+6pU=?Md zT8oHZq>x?XMhPqGuuMPBsDDP5>paF~23mm*_3J+NZ zRv3FrGJz*=0lOqaY81_oYpc0WY!eLxWI_pBS5Y^OL>9%yP-&GaT5$YF;)On=T4WHX zhtVVD6R9rJ2^YJKDMf>?-bf!8cvO{7{0K*5yv11lS12gzDC#~H@L>oq5O~6|)%ezM z;%$is^(`PhJvAnv~j7AV4Nm8AJqXl(-%ANl-B) zoCvDMht1Tocp__y-yBg_i-rX<%7PdwNehaxgm(s&0z|4P0TQ6MkKFZd1RH)1CDG7@FjiT~%#{m^T|JfQ;Hnr> z1kxv=KSBVWnhg2i5K-N9fWed#8w8BtEqp_U6O#c8NJD)^BetnfdKHDc%9B){P9!V2 zLGhyyG6{k%S*3!kApl6_h%ppKr}jz2!DQY6oDScKz!SesB*&uXpz4Y%V@R}-Ass;u zh7zL&L+=yfgt1hE5;mwJ6$SiY)Lm<@fp5s*ToD!mh*U!&U^23{2pIA&wI_t9j{Brx zq#DEpiz&bgm#Rv=iG#^_QS8p6DFgR|(7*6+t3E`p;sto3q|1cvq#FG3D~1G+`UnAj zDD?pY%nJTvbK-R9RU^eyG)ZI}0=+Mqe^G~39~j?opzo0YB03>u8##Nbx&jy)wu!N* z8Iod$3^tWbg9Sitm53n$M4GJjgy?stTv73B4>BE(_GskAiLv-#gG}uqsziHU6`m%B zgy>M!Pr&f~P8DKB4Hz&~vq>1dQ$ixd|8s_H4FQbmLscS%myr<=iu?>&mIsWgtOEuO zLT*Ns9STPIWr(9;RDmb*es^DsF;Xf$tc2 z*C97UXG?&GS~!J=gc=3odQkF?@R8K%n2^er0H$IXLHb0}dn_9N7}l}3i_#H}oeSsC zn|0tE)h>*agA20DPFxRD`~NS>DAOK24Cc)rH+R1O+)0xSU~&-;8dz934|Xy5e>Rc? AB>(^b delta 27162 zcmai630O_r+n3UGBorA!nG&jV&fe$ja}8#cNe)oCqvz_05*YvJ;z0-Pkk51X$K4o`wXDux@ynniN_|00m zd_#VG{fM>eA|H-3NV525azL^2Y|8;&1$*2dy!0GclsI5in@2~2lUnw=5#-s^qSmr{ zmlv;$E#i({y?;OXRGj~~C#O!?wmvcM)j!`J6~rAq8$7qg8Py!myvgf?fj`bCG&en9 zlx5q^X3nv}-=5{DdvzQ@mfIdEmsPv7^m zA{2!V58@gx4^tmuGxCy5`gqS9|2ek#y(|CNd#vcAWIiOF&Gg>;Zc)STS9g1Ht`dh& z)%yCUe$h(Z=3h^j-*?Mo-{X)K4sUB^DIcn@20sWYiCk1iT0OIRUk7RHkJT`ciSftf zRt_op#vL1uvO1q-J<41i;t-XbsOtFP#@C%OM->gS{g+rdpoM>`iOt8eCZ1M(W+rTD z`B?Yz_+bXz!Tm>QF|ZRYq)(!4C?F*qPkoLe8+3Cp6IvZeCyA6>Dmq{k3(vA z-*0ZcAb0bmu1|dYZe7Up&}%-bPQ2fP;m=m-l}-ywXf~eB|8VC$Hyyc)Y!t&duVPol<5^aafdc_8@DQ zzsO@}1HF8Yj~F|rSf6*h)3NuMDm2uAD95|PrW_7B17N6jk@w8xxy=l~s^+##OdT`dJ0>^ez~l5l@8aH`FWl^nTz<{! z&+!Mg1@sTpD_?*pG>aE9FSA#Ud>@$abI$gKvJF45e(&r}x;MSo$}U4vN@oh`1SzC* z>Am4r3H?XAb(xs;@Ws$qS8ttC+}jad&o9=_wK%o-);{;(QwLfH@U@IfMFMLO*XewX z4^O=C$a+|jZI8hDfrCf7{X2QE;?c=)fjj+n=d8$oyf+V?Qw^kGmUbyiex~#=507qL z*e|2ifvlBr?INegKg}-u&_0FR)J^-@^FKl)N*s>P&$$M{r}yZ!i@z{A_gb;TaA6OB zXnn-He;h_7Iz^?h^PAh2eElYms<*o2djGoSs7mq+ZbIc6c`Ze;UjykEOA(Nx|+kbzk7wr z$8_%c=~d3@CHL?0HAcMYSu}5kU-rwaZ@udseyrq&zVu7#F}PQp;`-t(nyx8eU zTl~Cl4DNq*PP+CFrMTs#>&ND##E-b#c}Q?nzu*v2wD?jHDm}Vij;uF`i*Z|cali;y z{n93-J5suzQQLKjdSymFG(4L5qqW_hsJtfM-vu2?UUXvfjsqu|X-!@Urw6sqJ763> zDb?{D&P3UEPF$KRjFgyuy}8rUrLozymNWW%ws&{TT^)Vlq)`x8GI>P-!`tV5?a}R` z&W{GAyIFb)#4}Y#-QC}R|8(ow0n3`{-V7Nbm}YPO8M|+4)6a|c>|JU3XFr$!k>)R- zZgVSUWZWL?vtwtk9kaa##BzSVzE1y+{K4sO$%^-le=v2{4DZ=9m5t`wzBp=ANptWy z#m&+WC!HuB@UZ#LJ$Du+M`wQbAKGVj(b0K@A&n(&X05K=vcB89O5139hn>+dMJ9{;M=`D`_$s^7pC~;{px!y$}=lzR?Yl!%~+}} z3^Hn3C9U2!bL8FWahnBxSEh?m2iHYO*S%c6e9AfSwd+Oe7p8wV<(O5i-p1+MlUH}0 z`=N>b{-|-$KJ84uXKcN=%qV@yC2*=t&?Xi;X!*4K{_w(1vyJ|Rk&O-a?2pPBHRgw` zl~(e&m2&5tze;dHmwxm5-V=t+U1allYGZ!7^3nUlTQ}{AF$;Hb_YiIjujl%@RZ0F} zA;C+sAnnRE7X{bzfmiCgJ8A5c_-}3_thY>!o@-dc`SOynA4LV^f3UYn|62>zWh8WV z$=&p;U|z!91xHivcwTqCSV;~`ToSp1=vF+!c6 zx~7$WiLJjb)Y@WJ38N`WBDA!*KNNCVS5^d9c;8S8naU5_iTpU!O>(Z1dLPSq!5+7yN;%xigH&Y3=H)yI?3AAilad|YIg zV0Ps2gx4#YUv_CaKcTVR4ehtP+8V7NF+X8N(V5@Abqt@jPrc6%Y45(JIB54l`)kt& zzTmlUXSR-Y?*H&xpM((Umjvzo)BHT-p0AJmD6F+<)0ntk5ufdD47Z$6Fww2kK+CAz zo{q8j-SNhjCslx-?u8iH1QL=1O zpF82_7FLgNj#K)1KkS*k)zfU(Qip${U+JFibH{Oe9{=1tN-EirA3I`t-%pNXjZb&F z*4p^9L&;0umx+y&Qx_Y1AGb5>*-$dtBm3~%1pB7z!s;|QY?!0yW;NzjYsuWC8pdzi z_8qMKJKj(Fw0@UqF6|AEIb~%nx^rGyb9#nZ^M2EZe0Z7Vn>6dq6zNxszBA4Hc5*Ih zP;%jR=j>5uQ}mri^V_`bdRK4Q!)`#^%M;!mDl%KSWFS>s1G?W_*WO~`dy{7Ye0*r`_QDM}eJuiYR!z_RWEkr2#My6_ufCjV`SHbJ z_3YD&FU974zIohV8mg2VIY-|pHnm#Vj-NSlYTkw*kGR7$9a3e(DrFHT`oJW)anx^{r5?oo8tmoW$u09`1-kO zWnHfoR+k4jF&Wcl$K*NZHl1xCbI^+4wS1&@&RnbS*{c?R)fxTDwomx8-$@dK?O}DC zBHTXmQL_S{M-X&EZv|Ea+yw zN7wo1=CXEMM&I+T_i+1ow^vU_+}-U{STpH;LXzK=Aj6TflYJU|^2%0tJ9qNf!8kv- zv?^@+s<6;N=O7F9^XGMhu*8MwuUhzao9z`G?v$?|&s#LEyG7by)>eO=Px-ZaHQ4Fw z=GmhCF2fOCCRQ6CsCx1jug-PQQ#be__pQnH5)>U=q#BSSB(ahQK-Axn=>h3wD`fM|)onzB=zmK@^`i4I9w*&|I z_I+7&SL(H+k@*~>jK-2j=S-I#J$-n=$WJwk0>)W+o_lXJE3B{kxS(&oQ>LxxSEBc> z^}|`Y&0k$O#ZNmrU`MSTX{HwsKbraK^m6xs(}ZX53mv)#&rEWg|L@Ai9fFTN3fMO^ zOIB16P{Xc!$&K{irJyP5%Bq{bq}pJI9xJ7|fpUeyYg_Bl82(^QH%czqqjQ^qGI2#V>ewe?V^B zoRY4o44<)IFRFq5#W_0``0Fpb^ya`#-SOu(M)aO`ZvPm)>2LH{FP+CH)+=PYCpYq4 zZ02w`J8`Gy-dkI}CQh_FaLc;d^pRO>7q7kbAnJwl;2NXtC!6aU3`?JW|BHR~tpSl6 zzMa!^KH+Pqo!zi&T#4?96}|8K85kR6Ptf15UNUEle%lGFmHKOA!{yZkc@FH29`Q*GuVTw(=$x< zBD)sUsqs8~;a!W&J0~YOF+Z9mFYR<{`PMz%l^?=``d1rP^7_onnBf8agoStA&)cep zj9+`nw(E|=SF?f-{hSdsDbW6IC_k)I_lpA;44o19%&qAykG%1(8h%scE;?vs7Q|=7 zjLBGZO&%I&YG|2yJ}-F0`0q1oW~_T=+Nx0(ljgN-V^Ytbax{ALW9z<#)%D!%KYY1d zFPObNv)JXZba}lOF&(XgoGl|l&o||I_tlcl94p`4x^*3QXG4<>?Y~Zw?HMPu`D}5+ z%XY?dPxo}EUe-?s8?`d)YThp(Pp{FL7izbw{QDvO3%l0PkL`J>&i+e^pg|4BJRj-g zzk7(;`3>0vE_HRfp6z%v@NC-W*l;YI@6$W^&h(x0jzVXkxfR8j&~W)+_h8Y z@FT|9-k3DgAky*Cq=(NohFJA8H`{$;YtMbvQnIvb?Am+$w){}f1LO84m`gL?PVjmD zk4L)C0k0Y|zUcg`5B4LcWe>XTn=!?ed!D$qTb<9+SzhH$7~VnbzNSlaozY z`|s~f(nC-084=Yozfnl?jJAG9vO8aOu=Q@U`r>n!M?nu#oX&--hR2MER3uq0Zr)>t z?cvld5j%Ho3$7LTXt#9(t;kJje9ng2He zgYRKKLk~YT=j!#VKd)MP!Pm>Hp9i-5-NdNz(c9Ozo80z}ng7VNWOK;g^vH zeRNrdu~o^#EkBtEQ#)uc*OoG6{$s2S zn}2kYezRG|XJIU|$xps*o?LruH<3_6&87`I< zMtysb;ayXw`T)Q6_a#G~r5hCeR83c&VV1jff8*Kxo%8h)r1_U|hZ6oBP~&6mhTdn!JdV_Dx!=%!-2ubuYqbNzP8S;w+3wVCN@|C_i+ZkI z$GeWsXi@aft!nMxP1S20-PF=*bBi8cDLdTl`FDr*Zm;L#mSJTZzT;4j2M-5^H2UDZ zzyJ2!G{v$uBbKgSwEN-ardy=5KJ>8B+E}-@Hfywbu{zG}Y_%Pe+v__Yj~zA1qFtAz z*%L!2@2t76hok?J?*Uyto6Kz`#JU;ZU3Wb;%gruzPsFWGDc5Q~y43MTi!<&&uAV)c z5WlC-@J1_>i}~<9!i1MY-#I-i9&=^$cjNEfm>lP|^Sh@Wvgv)WxaIKizFW*OI6Sn)$m9h-el^zLZ5#Hf{@8xE$JXgJ!f)B_)qU12n_aSP{lgXk z(-)0CwtLt0pmk%;$v-!WuVb^aqg8F)rS}adH+rz4?WM*}USmBvm|XwL*FNWC)MVAd zl7Ft9jy-DV_~hM*H!B7=d73oX%Wvqd{D`3&jGjv!hSXkul^m6smVV^H;nr$`7qUN-h9UD4q(L27f_fGwTuG53RHD94v zKE!(68fA}c^P~Ir^sasK{kgHyds|}f7x6#tZkrR&JU6eqWa4b~r8a?kw;Byu<2tJa zFIi_5xY=n#uNE2X673o(C6~YOi?aP*s^>kEeSTS|Xnh-b)4H?Yt=abeck>V7ZEmV| z^-VE;uq=Bcj!FXxS_6RMt#rN+5zwx{5)@=Cw*}(nkqMVL#^0Cg}%^zO%&M#TjYEK&M`_tHCuP!yJv2vhs>!GzT zuFja$XszF|*%KNSg}w0@;nZMI-NZCio8B8I?`yPp!`d^S?+-E=FfO&dRyUd7MX);$^9{Kxf8@v~?8genumZ~Hp;Zt`)OPv0P?Sesdw z-{ns_X&Ao8J}qoV;Xv!4TZb&$yqj;=_G){}SiRLv+P$~9;&Z{~j&t3COi-A0X3HhZ zHcc2F9CfvMX5VENVSOjt8BE}1r<+lc&Xs!h{K+|#q8G}=Ut)&Reg8AA?YLT)sPxLx_aYM_f_NC9SEEs zF>}tGb0zW7&kj%R;4ydq`oM|&#Oq&2PyEuwz23;}$Gblr+hyX|L+ve3YS~j%} z|51K1F>~&%OX&%BU&N_@oZ3I=>A%EOTv!}SaxJJgSW;d=KF<59b!*0s5b&H)iX${`Aq}}!f-FwWpGRQl8u>SQ!6RQpV zC;7y+dWQ~0eD{AG6nc4Pt!e|OuhbrBXjCw0%;XjwBb#?Ea&&Rtc*{_y_V83D#B1ST z)6Y3EC;Ak6n9aJEa(v~nIXkB9Co&yj5nb4Fe5_5A6;voDKh_8L0n zMZk;o<6p$Jx|EsIt@D{*NqlX$VJK9niw)#uU)k;7w-)H!nSCI9JMjT2@? z>vO-h9c2@LE=lU*sdIR>k;OHedI$PExu7ZxTVrmRAGv>lX}ZVcW6aLPOM2G&Yv+%w zf5GhXTfGe(LtCzWbnx2NvrnvDJU3grU#x5A`PS@;ZLI&rG3`q}bnd&*+$40-*?>&1 z1}|c_I&Ld?#2e-4=iT_8pYpVJbfTfBu)K~|ZK0*J_8SvUVJn9>qGI4V-~nGl8XjlI z3WHp<-A&X|Tb0s+RVr;|;#;ZQRvCWOu0eS1F#YhSM!Mk=nW0eNqCHv{zKUGU3H@U9 zO$9GkZHpR=QlwoMWs(r>sy(1_d1{N(NWyVK&NuDu2_wg8?`ZL#6jKAWk5(k%s_bLc z0qy6O|Jld=8``g|{s%<@=CR_xr;yc9r{VuTU((G=xgvRr0q!0a1ik`P?Y_Th+ zke7YdMz_Xr!T8{1PX6qnmDv^2Dr&{Y^ucQ;e9TTdzDZTQx=Y zI^9lA|2axfmhQ6u?!)h=u6<^CAFRea@JpPj`AxkA6@64mpto^3i*pA^-Seb zsv#q2g-J|=HoT;*2~)+}PvTpDysf9C-pb`FA%U0J*HcrkWKvepTY^_~pHkcscvSg{ zQDgV_h_7U1hm1<06pF-G)HE1QB^TV5Ry7u*RHy_&d_{MtP_aU;_==i5!$OC0%kYY$ zroqajLXh}MMZIDerEpJtMbA*KlnYh?5__R;0Jc%aF~ZycNr(DMYC9My%L$i2NlzGx zoe|8ISJl5-d5lCH~3Z#0Nq^`h* zO6p@PdWNEjN|+@+6hf5v;Dj{sp%(JQhfFY8jRzyxuP&Rzc=4xP2oWEwa7;wyQsI&K zlM}SV@K6cu@DNfrV3-MElKMhw6#f;8f1yFFNTC$2i4R8j5?0n}<2Cpx7rNph{IeDl zU9_fjMd2c@7S4!xyntT5u0bzqL9({A7e;)r!f^2+7Xt7QmswXhyB1VRAxp#+!jG~? zi**t!o;U?k8O{!Qlo^#$Mq`|c+*(GZR+2Kqs2NfOq%x&c{pWxQSJp`^G2=YU1v55T zj~UMt`P6GY_*A;W+4z^pU-8vH;$qUvfJYUyrs1G{z6MkIyud@~5`jLNM5236U?Ksm zkHA!>Mq+N95Fi5BUI-W{1o>LLdMUnQg&&b+8`@$6NX7{u8M&c!CueOan@1F0U~AIE z2P5Q(51H_2v81tJD@bYz&Bb1O3#H?Z7tlxZP;|aUK|9AHujf=;d`w4YQhmvevioG-V`GN}@VOH+BsrJATma-Qji!c zW06f14Sa>Amnmf^79cQAE#)*PQ$s3)01BDMBqls`j3&LDVP)jxg@R)_qCcgSlk?BG7vCLI)zdJ z4MAcY7Cc`R-4mn08CHEZ|OZEp1 zpp+6EwG;C?ujS&nV0 z5pw0qqJoZPVABaYmSI(dUXZd3X{n(*q?NTTUmDPZ5bM+I5I2iho+bt>3y+1u;IUAc zf>dT`t0etaDa#`qqE@1W`6N5U##0!49lEWWT!?7cNb@db)r!&+Ci4F;TM_aST#zbY zh{&VJBJ`tHjHc8C)oO<17R0Uci6}dNn2UmBwL;FS$PS_ADYhyQHV~7>R5;S{k-}l1 z1YM!v)s#+$H6R@eDXflynv6n8DI^6WBTot1!UNRYPMjP?nKOj-SFm59mLlD(c}$Qs zp%jurNtv5Yc*rZN+haHPnA;h2#DyHv?3DE|X4Fr^lhEVLduK@Liw5>_gNP;b=9;iM4P zN`xXWIfNLb@Z*3{CzfMa%JOousum~ZsZ*!ql#~I4>qt7)@FkU$(&I!UMot>WomvP~ zN||{4^;sd9O4>hV6_oI+WHR#9As}NZI;2WT-xiZ5{9=*+h~}u%q{e9{1rr!0ggEU~ zTTTTpj^u($#!_3SlCwMwGW6c5;L=bXvK%R_Dz%h6?oyRnMp<5!8va^Ye}A1|c%5YM zgT`!Sjaho0Rd8gJm z7BxErrm6y?EHCm7RIDOmWH7~JPdE}p+JXuN04PZo0+W+mkjbTL@|1`Nmb7)iREW+> z^e6B%Z za%w+ea%c|1tEd>Hyg%#`AzLyYxqc0+2)RYJO4Cx7VN@_HLpVxWUSQDG#H1M(F(rv9 zkvJlpU>W>aQhx(OPK9jDz_lPTmLsjF49Ng84+^Yu`l-wbmN7_JP%~uJRWybI*Blsj z13Fj}gUNxcQ<#j@qB4e8AuB}@QhFX8vZOf%rXU?-F^xwlDVBxmBF13F`zm0f{FiGf zum%>MB2-iChdWP3Vn}CEh6otERAia6qFu`t9ZJ9C*>#tzHLJE(8 zBE*rPFeGV6j8jnpr$8h?E>59>d8*V_&j}23|iI5J9q# z-mEmRI0mOu@d{EvID{qSAP7B4R|;E#Y$@3x!U(!Q=p>2`N$XPVum2~5`%QYXz>sDn z8UrJ96QDx`L(wVKbX!hyA>cn6tH_`M>1zd|m7(A?yQX?#fV2T(=`E_bI>IAD1n`s^=v;{`X46aS6?Ns4(GKBw!#<=d$ zm?2^YWb6q{j%yuKFoCH_Zpo3vAYF1`P?>+mDEI#WB3-RHu!xb#FksZ#hPQ_FJkb`3 zeNy^CM`dJrfK=*MbpJ?)0i2=IdZ2?^AqU2l6Ezu@A&nRuN6{N9wLcXq7$8zcWQVx! z*AyD0^x%^aibgII^YcVo#Jr?8gp5cf4*v1~&={9aL}MBDm>d`vYMKfJ0>`kpI3n6& z5K3m`xX&XO2mW)UWkWcMGeL9+M~cSeRE`m*6MILtmCNacAa4M=QvR1)Ksi!Kr2ZB$ z(xU+;s{%u{A~7%{V<@3!;3NH&k23!cw>-#n6OG{_QTvVqL17RrLaB&bLgfu`$%D)z zF@}PN%^?Ry;)OJysCHpU6AcUpndlH0nf?QvvI?wBwnY9{u%8uG7u6x|TuBT;FF9#e zEhk+#te2r=3p%PYR=WQpW5^s0Fe-%y4EZo>GL*-Q%Q)N(5_C#z9VwW&Gex*g5TfRY zl5}Kg!6(r){N)xWS_ukMQv-8Q%zp*0H0dr^s1dD*ZL$6e_AAuXy&?#yy9Z#jUcfCS zb?ktV7a_>FYX&BbF$dg&P2Rp1wrUx&^b{c`NnnT$$$SLjN+=DAjwSsN&{3`?Fy$X7 z>#u_hkfRJ7a`dEIk17z<4iE#Y5pR;13V9H6aR|dT{U7jOtoop|t{P%X`YR}GpcEo7 zO8J5NRbcc%6Xh4{$8J|Tr8sAlLPQFY%twfKv6?Hvy(_#nYTt3wMx7k!&GN9c!a)}6 zTu4(P)^Jua2<2-+HH$aC#6pxv$I%$!Jn>+^66%T~L`f%wAtz5^h#568WQdUGp)eQ_ z3d40V;UJ?im8S$oWYtAhi|IJ3Lp7N$#5Fpc6M_y`>lN9catE1_;#!%UHQW&L%%Ve% zj62X)^VS6$feaR=O%04P5>1Rphe|S?FAhSc2F2@p(pv{csaE)fh*2~L0frJK8iTeK z_1|CV2cevv44h_ikC96v4IvyaRGtxoz!@V$M$oZ;we24V8K8p9%AmA?GDIl-MI28I zB3DrtET#EYaxy4OA)Rn$pbRRBP@`Z!Wm!>O##Avdl`0Uc?}Sh1U|8WDk@J*N4myf0 zAdD0p^gWGL$n2lO7PO^qA{c}UW$te!!wL{5i0lx)E{#!lG-ylR(Ey|FXv9l=4L3*> zEbIc6$OVHS@Wf;gp_A!PF=b2IGGMeshwU&9L4$Vr2V+%B%|H9HiqZb^p@4mQ}ONrwsr`h?>MXr6lTAB1x1 z*aIV-EYP8Hha3bOO6V%FN|H3nVi~iBLnUsz3dI6EDgC%3BV7bk*+4%N1H)iax=*~r zBe{ShN;+uh4|%n+{)8)O5{ajhOxGgWq9VN0ERTMX1CaOH7RhSM@kAS+TwFQvMurorJSyG-Vr2HsZL2P*Bb=W=)k%p)uu3pxjnLMyH5utHPudXK;v-=_6JG=93bPL8z=D zDsZW&6Bt}%qCa>8q_2YuC{*f|()~T};;j_vD2W)UEPzo3YS6V*pal%h4lx-~VQNee z{0AYU36^4%w+=cgw+@U75|IR>l6(k6sBaHSOUy_|8jWE^$yuYah44+VB}$4;RE~sp zM1NHMhB!ufs9NDjLygMG!X0GGsA~*h3esB#hH@xk(x|kn0wYr?4aJf?xM);kwNU4Z zSUgWEFx1d!)|ZP12x$ZIRty-bfQfGL4FZ`d0Uat62s$X|GEuCo9-ud=BaUk}(szb_ zrR)qijw}W>2z4`#N?OfbEmaQ+4wWWEW87kqfja_Lsva2_%4x~Au!L1$qd0CWz3l{Bo*7Js`{|zZ|Sc%>ly)s8nNyEQp0MLZXADLmd)vJYi~~ z?ltI8jSsCvo>{C8^`RPsqMfEP;<%Kda`k~n|9ilDn%1P~v$x!JK z9s|SkqZ-4vAu&LtAc;#MCyh}mR`(|&L*)sSI@KT4nT2rtrldAL2_X4~8#7Ae5z$ba zh~zOP-iU_inL<#g6~TQjT~Y17k-h?vkCGF-GbH258bMU`QgO0)_eJ%KxQ0v;pl1n;@3iuUe zc3V}6Klmts+Dha&NkdZkV+fX66Bo(YTjiso;?kgDPucbto316slP6Mrak0t>gr;Nz{$85>7Gdurvn&i0vaW>^s3~u`!ulfTSx) zIfA~TN|%8B%{iBwAAr<{8%Vn&olHyyMUsRC$CoQ8gdp}E6I79E0{v0n>p=6vsUl1R z2&vLSV2C}b$xt78!CyzJoTNjTLrf5=P0aI&?)P8iCN?Gm#nL)Ea*tt`$s-0g7gzh_ zz=##7wosN8Tl^=@i{c~XhRDXES|uB5(R#K4GTNc#>8AQx4CDgVp09za}cYBt906y<^f!*@^wogB$? z8pClS=nz#`;;QV=Uw08JstH1|DuBk+f)@T71y3a-@YMtLu@Nv7@{`@-vrvLT9Lo5~ zD|KMhtvE2Tb{hf;{KYVgVC>SfW>XY$Jo?| ziNKH|A_u`)p$6e_FHO-==|0$a1nMN6O3sr$9|)=1H-rURENlB(ioZiLR=Ps_DEFAt{l93jjIvLKIJ$B}7|Cyd%hN7vJ N$_-n!a_Q@8_Using Random Residuals for Censored Data in EGRET

Robert M. Hirsch and Laura A. De Cicco

-

17 June, 2016

+

21 June, 2016

@@ -100,12 +100,12 @@

2 The concept of randomized estim

What is most important to understand is that these randomly generated values are never used in any of the WRTDS computations that result in estimates of daily concentrations or fluxes, or annual average concentrations or fluxes, or any trends. These randomly generated values are created strictly to provide a more-easily interpreted set of diagonstic graphics. Here is an example of a plot of concentration versus time using this approach.

The solid circles are the uncensored observations and the open circles are these rObserved values. One thing to note is that if these values were to be generated a second time the graph would look slightly different. There is no unique set of values that will be plotted. The random numbers used will be different if we go back to the same data set and generate them again. Below we show such a figure and then redo the graphic with a second set of random numbers. The code is shown with the graphics.

eList <- makeAugmentedSample(eList)
-plotConcQ(eList, qUnit = 4, rResid = TRUE)
-

+plotConcQ(eList, qUnit = 4, randomCensored = TRUE)

+

# now do it all over again
 eList <- makeAugmentedSample(eList)
-plotConcQ(eList, qUnit = 4, rResid = TRUE)
-

+plotConcQ(eList, qUnit = 4, randomCensored = TRUE)
+

Careful examination of these two figures reveals that the black dots are exactly the same in both, but the open circles are different between the two. In looking at these kinds of plots we are not necessarily looking to see what actually happened on a particular day, but rather to understand the pattern of the relationship between the two variables being plotted.

@@ -113,10 +113,10 @@

3 Randomized residuals

Many of the diagnostics we use in EGRET are plots of residuals (on the y-axis) versus predicted value, or time, or discharge (on the x-axis). Once we have this rObserved value, we can compute a randomized residual which is called rResid (again the “r” denotes that they are randomly generated residuals) and is stored as Sample$rResid in EGRET. The computation is this:

rResid = rObserved - predicted value. All three variables in this equation are in log concentration units.

We can look at our residuals plots in the following manner (using the second set of rObserved values computed above).

-
plotResidTime(eList, rResid = TRUE)
-

-
plotResidQ(eList, qUnit = 4, rResid = TRUE)
-

+
plotResidTime(eList, randomCensored = TRUE)
+

+
plotResidQ(eList, qUnit = 4, randomCensored = TRUE)
+

4 Details for how to include random residuals in your computations

@@ -126,18 +126,18 @@

4 Details for how to include rand
  • Run modelEstimation in the usual way: eList <- modelEstimation(eList)

  • Then augment the Sample data frame that is inside of the eList with the command:

    eList <- makeAugmentedSample(eList)

    • -

  • To produce any one of the following graphics, using the rObserved or rResid values simply add the argument rResid = TRUE to the call to the graphical function. Note that it doesn’t matter if what is being plotted is an observed value or a residual, the argument is always rResid = TRUE. The original option of showing the vertical lines for censored values remains available as the default for any of these functions. The call can either say rResid = FALSE or just not include rResid in the argument list and the graphs will appear without these random values and without the open circle/closed circle symbology. The censored values or censored residuals will be shown as the vertical lines. The functions where this approach applies are these: plotConcPred, plotConcQ, plotConcTime, plotConcTimeDaily, plotFluxPred, plotFluxQ, plotFluxTImeDaily, plotResidPred, plotResidQ, plotResidTime. It is also available in the two multiple plot functions: multiPlotOverview, and fluxBiasMulti

    • +

  • To produce any one of the following graphics, using the rObserved or rResid values simply add the argument randomCensored = TRUE to the call to the graphical function. Note that it doesn’t matter if what is being plotted is an observed value or a residual, the argument is always randomCensored = TRUE. The original option of showing the vertical lines for censored values remains available as the default for any of these functions. The call can either say randomCensored = FALSE or just not include randomCensored in the argument list and the graphs will appear without these random values and without the open circle/closed circle symbology. The censored values or censored residuals will be shown as the vertical lines. The functions where this approach applies are these: plotConcPred, plotConcQ, plotConcTime, plotConcTimeDaily, plotFluxPred, plotFluxQ, plotFluxTImeDaily, plotResidPred, plotResidQ, plotResidTime. It is also available in the two multiple plot functions: multiPlotOverview, and fluxBiasMulti

    • For example.

    -
    multiPlotDataOverview(eList, qUnit = 4, rResid = TRUE)
    -

    -
    fluxBiasMulti(eList, qUnit = 4, fluxUnit = 9, rResid = TRUE)
    -

    +
    multiPlotDataOverview(eList, qUnit = 4, randomCensored = TRUE)
    +

    +
    fluxBiasMulti(eList, qUnit = 4, fluxUnit = 9, randomCensored = TRUE)
    +

    5 Two final thoughts

    -

    In the EGRET User Guide (http://pubs.usgs.gov/tm/04/a10/) the distinction is made between the graphical methods that are used to simply describe the data (and these graphics shown in multiPlotDataOverview) as distinct from graphical methods that are used to describe the WRTDS model of the system (such as the graphs in fluxBiasMulti). If the rResid = TRUE option is used with multiPlotDataOverview, or other graphical functions that normally don’t depend on the WRTDS model, they now become a hybrid, because they are using the WRTDS model to generate the random values used in the graphs. Take a graph such as plotConcQ. With rResid = FALSE it is a pure representation of the data. No assumptions are being made. But, when rResid = TRUE, it is now a representation of the data which is partly based on an assumption that the fitted WRTDS model is indeed a correct model. The fitted WRTDS model is partly determining the placement of the random values that are less than the reporting limit. If the analyst wants a “pure” representation of the data without any assumed model for graphing, then the rResid option should be set to FALSE.

    -

    Also, with figures such as shown in fluxBiasMulti, there is a kind of circularity in the logic. The circularity is this: we are using the graphs to assess the adequacy of the model fit, but we are using the model to estimate some of the observations. This circularity is not a fatal flaw to the approach, but is a reality that the user should consider. On balance, the authors think that using rResid = TRUE in all of the plots to which it applies is a beneficial approach because it enhances the ability of the analyst to interpret the figures. But, we reiterate here, the choice of using the random approach or not has no bearing whatsoever on the quantitative outputs that the WRTDS method in the EGRET package produces.

    +

    In the EGRET User Guide (http://pubs.usgs.gov/tm/04/a10/) the distinction is made between the graphical methods that are used to simply describe the data (and these graphics shown in multiPlotDataOverview) as distinct from graphical methods that are used to describe the WRTDS model of the system (such as the graphs in fluxBiasMulti). If the randomCensored = TRUE option is used with multiPlotDataOverview, or other graphical functions that normally don’t depend on the WRTDS model, they now become a hybrid, because they are using the WRTDS model to generate the random values used in the graphs. Take a graph such as plotConcQ. With randomCensored = FALSE it is a pure representation of the data. No assumptions are being made. But, when randomCensored = TRUE, it is now a representation of the data which is partly based on an assumption that the fitted WRTDS model is indeed a correct model. The fitted WRTDS model is partly determining the placement of the random values that are less than the reporting limit. If the analyst wants a “pure” representation of the data without any assumed model for graphing, then the randomCensored option should be set to FALSE.

    +

    Also, with figures such as shown in fluxBiasMulti, there is a kind of circularity in the logic. The circularity is this: we are using the graphs to assess the adequacy of the model fit, but we are using the model to estimate some of the observations. This circularity is not a fatal flaw to the approach, but is a reality that the user should consider. On balance, the authors think that using randomCensored = TRUE in all of the plots to which it applies is a beneficial approach because it enhances the ability of the analyst to interpret the figures. But, we reiterate here, the choice of using the random approach or not has no bearing whatsoever on the quantitative outputs that the WRTDS method in the EGRET package produces.

    diff --git a/man/boxResidMonth.Rd b/man/boxResidMonth.Rd index 201f5d6b..b681684e 100644 --- a/man/boxResidMonth.Rd +++ b/man/boxResidMonth.Rd @@ -6,7 +6,7 @@ \usage{ boxResidMonth(eList, stdResid = FALSE, las = 1, printTitle = TRUE, cex = 0.8, cex.axis = 1.1, cex.main = 1.1, font.main = 2, - tinyPlot = FALSE, customPar = FALSE, rResid = FALSE, ...) + tinyPlot = FALSE, customPar = FALSE, randomCensored = FALSE, ...) } \arguments{ \item{eList}{named list with at least the Sample and INFO dataframes} @@ -29,7 +29,7 @@ boxResidMonth(eList, stdResid = FALSE, las = 1, printTitle = TRUE, \item{customPar}{logical defaults to FALSE. If TRUE, par() should be set by user before calling this function} -\item{rResid}{logical. Show censored residuals as randomized.} +\item{randomCensored}{logical. Show censored residuals as randomized.} \item{\dots}{arbitrary graphical parameters that will be passed to genericEGRETDotPlot function (see ?par for options)} } diff --git a/man/fluxBiasMulti.Rd b/man/fluxBiasMulti.Rd index 4028e07a..518bb0ae 100644 --- a/man/fluxBiasMulti.Rd +++ b/man/fluxBiasMulti.Rd @@ -5,7 +5,7 @@ \title{Produces 8-panel plot that is useful for determining if there is a flux bias problem} \usage{ fluxBiasMulti(eList, qUnit = 2, fluxUnit = 3, moreTitle = "WRTDS", - cex = 0.7, cex.axis = 1.1, cex.main = 1.1, rResid = FALSE, + cex = 0.7, cex.axis = 1.1, cex.main = 1.1, randomCensored = FALSE, col = "black", lwd = 1, ...) } \arguments{ @@ -23,7 +23,7 @@ fluxBiasMulti(eList, qUnit = 2, fluxUnit = 3, moreTitle = "WRTDS", \item{cex.main}{magnification to be used for main titles relative to the current setting of cex} -\item{rResid}{logical. Show censored residuals as randomized.} +\item{randomCensored}{logical. Show censored residuals as randomized.} \item{col}{color of points on plot, see ?par 'Color Specification'} diff --git a/man/multiPlotDataOverview.Rd b/man/multiPlotDataOverview.Rd index 55da8b4c..9836c682 100644 --- a/man/multiPlotDataOverview.Rd +++ b/man/multiPlotDataOverview.Rd @@ -4,8 +4,8 @@ \alias{multiPlotDataOverview} \title{Produces a 4 panel plot that gives an overview of the data set prior to any processing} \usage{ -multiPlotDataOverview(eList, qUnit = 2, cex.main = 1.2, rResid = FALSE, - logScaleConc = TRUE, logScaleQ = TRUE) +multiPlotDataOverview(eList, qUnit = 2, cex.main = 1.2, + randomCensored = FALSE, logScaleConc = TRUE, logScaleQ = TRUE) } \arguments{ \item{eList}{named list with at least Daily, Sample, and INFO dataframes} @@ -14,7 +14,7 @@ multiPlotDataOverview(eList, qUnit = 2, cex.main = 1.2, rResid = FALSE, \item{cex.main}{magnification to be used for main titles relative to the current setting of cex} -\item{rResid}{logical. Show censored values as randomized.} +\item{randomCensored}{logical. Show censored values as randomized.} \item{logScaleConc}{logical if TRUE y in concentration graphs plotted in log axis. Default is TRUE.} diff --git a/man/plotConcPred.Rd b/man/plotConcPred.Rd index d1e8ebbe..1426f1a0 100644 --- a/man/plotConcPred.Rd +++ b/man/plotConcPred.Rd @@ -6,7 +6,8 @@ \usage{ plotConcPred(eList, concMax = NA, logScale = FALSE, printTitle = TRUE, tinyPlot = FALSE, cex = 0.8, cex.axis = 1.1, cex.main = 1.1, - customPar = FALSE, col = "black", lwd = 1, rResid = FALSE, ...) + customPar = FALSE, col = "black", lwd = 1, randomCensored = FALSE, + ...) } \arguments{ \item{eList}{named list with at least the Sample and INFO dataframes} @@ -32,7 +33,7 @@ plotConcPred(eList, concMax = NA, logScale = FALSE, printTitle = TRUE, \item{lwd}{number line width} -\item{rResid}{logical. Show censored values as randomized.} +\item{randomCensored}{logical. Show censored values as randomized.} \item{\dots}{arbitrary graphical parameters that will be passed to genericEGRETDotPlot function (see ?par for options)} } diff --git a/man/plotConcQ.Rd b/man/plotConcQ.Rd index 9c063cf6..19498afa 100644 --- a/man/plotConcQ.Rd +++ b/man/plotConcQ.Rd @@ -5,7 +5,7 @@ \title{Plot of Observed Concentration versus Discharge} \usage{ plotConcQ(eList, qUnit = 2, tinyPlot = FALSE, logScale = FALSE, - rResid = FALSE, concMax = NA, concMin = NA, printTitle = TRUE, + randomCensored = FALSE, concMax = NA, concMin = NA, printTitle = TRUE, cex = 0.8, cex.axis = 1.1, cex.main = 1.1, rmSciX = FALSE, rmSciY = FALSE, customPar = FALSE, col = "black", lwd = 1, ...) } @@ -18,7 +18,7 @@ plotConcQ(eList, qUnit = 2, tinyPlot = FALSE, logScale = FALSE, \item{logScale}{logical if TRUE x and y plotted in log axis} -\item{rResid}{logical. Show censored values as randomized.} +\item{randomCensored}{logical. Show censored values as randomized.} \item{concMax}{number specifying the maximum value to be used on the vertical axis, default is NA (which allows it to be set automatically by the data)} diff --git a/man/plotConcTime.Rd b/man/plotConcTime.Rd index 2f987538..2f181944 100644 --- a/man/plotConcTime.Rd +++ b/man/plotConcTime.Rd @@ -4,10 +4,10 @@ \alias{plotConcTime} \title{Plot of Observed Concentration versus Time} \usage{ -plotConcTime(eList, qUnit = 2, qLower = NA, qUpper = NA, rResid = FALSE, - tinyPlot = FALSE, concMax = NA, concMin = NA, printTitle = TRUE, - logScale = FALSE, cex = 0.8, cex.axis = 1.1, cex.main = 1.1, - customPar = FALSE, col = "black", lwd = 1, ...) +plotConcTime(eList, qUnit = 2, qLower = NA, qUpper = NA, + randomCensored = FALSE, tinyPlot = FALSE, concMax = NA, concMin = NA, + printTitle = TRUE, logScale = FALSE, cex = 0.8, cex.axis = 1.1, + cex.main = 1.1, customPar = FALSE, col = "black", lwd = 1, ...) } \arguments{ \item{eList}{named list with at least the Sample and INFO dataframes} @@ -18,7 +18,7 @@ plotConcTime(eList, qUnit = 2, qLower = NA, qUpper = NA, rResid = FALSE, \item{qUpper}{numeric the upper bound on values of discharge for selection of data points to be plotted, units are those specified by qUnit, default = NA which is equivalent to an upper bound of infinity} -\item{rResid}{logical. Show censored values as randomized.} +\item{randomCensored}{logical. Show censored values as randomized.} \item{tinyPlot}{logical variable, if TRUE plot is designed to be plotted small as part of a multipart figure, default is FALSE.} @@ -66,7 +66,7 @@ plotConcTime(eList) eList <- setPA(eList, paStart=6,paLong=3) plotConcTime(eList, qUnit = 1, qLower = 100, qUpper = 10000) plotConcTime(eList, logScale=TRUE) -plotConcTime(eList, qUnit = 1, qLower = 100, qUpper = 10000, rResid = TRUE) +plotConcTime(eList, qUnit = 1, qLower = 100, qUpper = 10000, randomCensored = TRUE) } \seealso{ \code{\link{selectDays}}, \code{\link{genericEGRETDotPlot}} diff --git a/man/plotConcTimeDaily.Rd b/man/plotConcTimeDaily.Rd index 6378c4da..3c3ce0f4 100644 --- a/man/plotConcTimeDaily.Rd +++ b/man/plotConcTimeDaily.Rd @@ -6,8 +6,8 @@ \usage{ plotConcTimeDaily(eList, startYear = NA, endYear = NA, tinyPlot = FALSE, concMax = NA, printTitle = TRUE, cex = 0.8, cex.axis = 1.1, - rResid = FALSE, cex.main = 1.1, customPar = FALSE, col = "black", - lwd = 1, prettyDate = TRUE, ...) + randomCensored = FALSE, cex.main = 1.1, customPar = FALSE, + col = "black", lwd = 1, prettyDate = TRUE, ...) } \arguments{ \item{eList}{named list with at least the Daily, Sample, and INFO dataframes} @@ -26,7 +26,7 @@ plotConcTimeDaily(eList, startYear = NA, endYear = NA, tinyPlot = FALSE, \item{cex.axis}{magnification to be used for axis annotation relative to the current setting of cex} -\item{rResid}{logical. Show censored values as randomized.} +\item{randomCensored}{logical. Show censored values as randomized.} \item{cex.main}{magnification to be used for main titles relative to the current setting of cex} diff --git a/man/plotFluxPred.Rd b/man/plotFluxPred.Rd index e3f914b5..27a444c9 100644 --- a/man/plotFluxPred.Rd +++ b/man/plotFluxPred.Rd @@ -7,7 +7,7 @@ plotFluxPred(eList, fluxUnit = 3, fluxMax = NA, printTitle = TRUE, oneToOneLine = TRUE, customPar = FALSE, col = "black", lwd = 1, cex = 0.8, cex.axis = 1.1, cex.main = 1.1, tinyPlot = FALSE, - logScale = FALSE, rResid = FALSE, ...) + logScale = FALSE, randomCensored = FALSE, ...) } \arguments{ \item{eList}{named list with at least the Sample and INFO dataframes} @@ -37,7 +37,7 @@ plotFluxPred(eList, fluxUnit = 3, fluxMax = NA, printTitle = TRUE, \item{logScale}{logical if TRUE x and y plotted in log axis} -\item{rResid}{logical. Show censored values as randomized.} +\item{randomCensored}{logical. Show censored values as randomized.} \item{\dots}{arbitrary graphical parameters that will be passed to genericEGRETDotPlot function (see ?par for options)} } diff --git a/man/plotResidPred.Rd b/man/plotResidPred.Rd index 79f52f60..f65cd00d 100644 --- a/man/plotResidPred.Rd +++ b/man/plotResidPred.Rd @@ -6,8 +6,8 @@ \usage{ plotResidPred(eList, stdResid = FALSE, tinyPlot = FALSE, printTitle = TRUE, col = "black", lwd = 1, cex = 0.8, - cex.axis = 1.1, cex.main = 1.1, customPar = FALSE, rResid = FALSE, - ...) + cex.axis = 1.1, cex.main = 1.1, customPar = FALSE, + randomCensored = FALSE, ...) } \arguments{ \item{eList}{named list with at least the Sample and INFO dataframes} @@ -31,7 +31,7 @@ plotResidPred(eList, stdResid = FALSE, tinyPlot = FALSE, \item{customPar}{logical defaults to FALSE. If TRUE, par() should be set by user before calling this function (for example, adjusting margins with par(mar=c(5,5,5,5))). If customPar FALSE, EGRET chooses the best margins depending on tinyPlot.} -\item{rResid}{logical. Show censored residuals as randomized.} +\item{randomCensored}{logical. Show censored residuals as randomized.} \item{\dots}{arbitrary graphical parameters that will be passed to genericEGRETDotPlot function (see ?par for options)} } diff --git a/man/plotResidQ.Rd b/man/plotResidQ.Rd index 5c7297ee..cb624b72 100644 --- a/man/plotResidQ.Rd +++ b/man/plotResidQ.Rd @@ -7,7 +7,7 @@ plotResidQ(eList, qUnit = 2, tinyPlot = FALSE, stdResid = FALSE, printTitle = TRUE, col = "black", lwd = 1, cex = 0.8, cex.axis = 1.1, cex.main = 1.1, rmSciX = FALSE, customPar = FALSE, - rResid = FALSE, ...) + randomCensored = FALSE, ...) } \arguments{ \item{eList}{named list with at least the Sample and INFO dataframes} @@ -35,7 +35,7 @@ plotResidQ(eList, qUnit = 2, tinyPlot = FALSE, stdResid = FALSE, \item{customPar}{logical defaults to FALSE. If TRUE, par() should be set by user before calling this function (for example, adjusting margins with par(mar=c(5,5,5,5))). If customPar FALSE, EGRET chooses the best margins depending on tinyPlot.} -\item{rResid}{logical. Show censored residuals as randomized.} +\item{randomCensored}{logical. Show censored residuals as randomized.} \item{\dots}{arbitrary graphical parameters that will be passed to genericEGRETDotPlot function (see ?par for options)} } diff --git a/man/plotResidTime.Rd b/man/plotResidTime.Rd index e87198a0..f8273a16 100644 --- a/man/plotResidTime.Rd +++ b/man/plotResidTime.Rd @@ -6,7 +6,7 @@ \usage{ plotResidTime(eList, stdResid = FALSE, printTitle = TRUE, hLine = TRUE, tinyPlot = FALSE, col = "black", lwd = 1, cex = 0.8, cex.axis = 1.1, - cex.main = 1.1, customPar = FALSE, rResid = FALSE, ...) + cex.main = 1.1, customPar = FALSE, randomCensored = FALSE, ...) } \arguments{ \item{eList}{named list with at least the Sample and INFO dataframes} @@ -32,7 +32,7 @@ plotResidTime(eList, stdResid = FALSE, printTitle = TRUE, hLine = TRUE, \item{customPar}{logical defaults to FALSE. If TRUE, par() should be set by user before calling this function (for example, adjusting margins with par(mar=c(5,5,5,5))). If customPar FALSE, EGRET chooses the best margins depending on tinyPlot.} -\item{rResid}{logical. Show censored residuals as randomized.} +\item{randomCensored}{logical. Show censored residuals as randomized.} \item{\dots}{arbitrary graphical parameters that will be passed to genericEGRETDotPlot function (see ?par for options)} } diff --git a/vignettes/rResid.Rmd b/vignettes/rResid.Rmd index 7e76d81a..b158d9e2 100644 --- a/vignettes/rResid.Rmd +++ b/vignettes/rResid.Rmd @@ -78,10 +78,10 @@ The solid circles are the uncensored observations and the open circles are these ```{r} eList <- makeAugmentedSample(eList) -plotConcQ(eList, qUnit = 4, rResid = TRUE) +plotConcQ(eList, qUnit = 4, randomCensored = TRUE) # now do it all over again eList <- makeAugmentedSample(eList) -plotConcQ(eList, qUnit = 4, rResid = TRUE) +plotConcQ(eList, qUnit = 4, randomCensored = TRUE) ``` Careful examination of these two figures reveals that the black dots are exactly the same in both, but the open circles are different between the two. In looking at these kinds of plots we are not necessarily looking to see what actually happened on a particular day, but rather to understand the pattern of the relationship between the two variables being plotted. @@ -95,8 +95,8 @@ rResid = rObserved - predicted value. All three variables in this equation are We can look at our residuals plots in the following manner (using the second set of rObserved values computed above). ```{r} -plotResidTime(eList, rResid = TRUE) -plotResidQ(eList, qUnit = 4, rResid = TRUE) +plotResidTime(eList, randomCensored = TRUE) +plotResidQ(eList, qUnit = 4, randomCensored = TRUE) ``` # Details for how to include random residuals in your computations @@ -109,17 +109,17 @@ In order to use the **rObserved** and **rResid** in making graphs in EGRET the p **eList <- makeAugmentedSample(eList)** -4. To produce any one of the following graphics, using the **rObserved** or **rResid** values simply add the argument **rResid = TRUE** to the call to the graphical function. Note that it doesn't matter if what is being plotted is an observed value or a residual, the argument is always **rResid = TRUE**. The original option of showing the vertical lines for censored values remains available as the default for any of these functions. The call can either say **rResid = FALSE** or just not include **rResid** in the argument list and the graphs will appear without these random values and without the open circle/closed circle symbology. The censored values or censored residuals will be shown as the vertical lines. The functions where this approach applies are these: **plotConcPred, plotConcQ, plotConcTime, plotConcTimeDaily, plotFluxPred, plotFluxQ, plotFluxTImeDaily, plotResidPred, plotResidQ, plotResidTime**. It is also available in the two multiple plot functions: **multiPlotOverview, and fluxBiasMulti** +4. To produce any one of the following graphics, using the **rObserved** or **rResid** values simply add the argument **randomCensored = TRUE** to the call to the graphical function. Note that it doesn't matter if what is being plotted is an observed value or a residual, the argument is always **randomCensored = TRUE**. The original option of showing the vertical lines for censored values remains available as the default for any of these functions. The call can either say **randomCensored = FALSE** or just not include **randomCensored** in the argument list and the graphs will appear without these random values and without the open circle/closed circle symbology. The censored values or censored residuals will be shown as the vertical lines. The functions where this approach applies are these: **plotConcPred, plotConcQ, plotConcTime, plotConcTimeDaily, plotFluxPred, plotFluxQ, plotFluxTImeDaily, plotResidPred, plotResidQ, plotResidTime**. It is also available in the two multiple plot functions: **multiPlotOverview, and fluxBiasMulti** For example. ```{r, fig.height = 9, fig.width = 8} -multiPlotDataOverview(eList, qUnit = 4, rResid = TRUE) -fluxBiasMulti(eList, qUnit = 4, fluxUnit = 9, rResid = TRUE) +multiPlotDataOverview(eList, qUnit = 4, randomCensored = TRUE) +fluxBiasMulti(eList, qUnit = 4, fluxUnit = 9, randomCensored = TRUE) ``` # Two final thoughts -In the EGRET User Guide (http://pubs.usgs.gov/tm/04/a10/) the distinction is made between the graphical methods that are used to simply describe the data (and these graphics shown in **multiPlotDataOverview**) as distinct from graphical methods that are used to describe the WRTDS model of the system (such as the graphs in **fluxBiasMulti**). If the **rResid = TRUE** option is used with **multiPlotDataOverview**, or other graphical functions that normally don't depend on the WRTDS model, they now become a hybrid, because they are using the WRTDS model to generate the random values used in the graphs. Take a graph such as **plotConcQ**. With **rResid = FALSE** it is a pure representation of the data. No assumptions are being made. But, when **rResid = TRUE**, it is now a representation of the data which is partly based on an assumption that the fitted WRTDS model is indeed a correct model. The fitted WRTDS model is partly determining the placement of the random values that are less than the reporting limit. If the analyst wants a "pure" representation of the data without any assumed model for graphing, then the **rResid** option should be set to **FALSE**. +In the EGRET User Guide (http://pubs.usgs.gov/tm/04/a10/) the distinction is made between the graphical methods that are used to simply describe the data (and these graphics shown in **multiPlotDataOverview**) as distinct from graphical methods that are used to describe the WRTDS model of the system (such as the graphs in **fluxBiasMulti**). If the **randomCensored = TRUE** option is used with **multiPlotDataOverview**, or other graphical functions that normally don't depend on the WRTDS model, they now become a hybrid, because they are using the WRTDS model to generate the random values used in the graphs. Take a graph such as **plotConcQ**. With **randomCensored = FALSE** it is a pure representation of the data. No assumptions are being made. But, when **randomCensored = TRUE**, it is now a representation of the data which is partly based on an assumption that the fitted WRTDS model is indeed a correct model. The fitted WRTDS model is partly determining the placement of the random values that are less than the reporting limit. If the analyst wants a "pure" representation of the data without any assumed model for graphing, then the **randomCensored** option should be set to **FALSE**. -Also, with figures such as shown in **fluxBiasMulti**, there is a kind of circularity in the logic. The circularity is this: we are using the graphs to assess the adequacy of the model fit, but we are using the model to estimate some of the observations. This circularity is not a fatal flaw to the approach, but is a reality that the user should consider. On balance, the authors think that using **rResid = TRUE** in all of the plots to which it applies is a beneficial approach because it enhances the ability of the analyst to interpret the figures. But, we reiterate here, the choice of using the random approach or not has no bearing whatsoever on the quantitative outputs that the WRTDS method in the EGRET package produces. +Also, with figures such as shown in **fluxBiasMulti**, there is a kind of circularity in the logic. The circularity is this: we are using the graphs to assess the adequacy of the model fit, but we are using the model to estimate some of the observations. This circularity is not a fatal flaw to the approach, but is a reality that the user should consider. On balance, the authors think that using **randomCensored = TRUE** in all of the plots to which it applies is a beneficial approach because it enhances the ability of the analyst to interpret the figures. But, we reiterate here, the choice of using the random approach or not has no bearing whatsoever on the quantitative outputs that the WRTDS method in the EGRET package produces. From df387a5d438ff040fbbd1930956f8c90018f2fa5 Mon Sep 17 00:00:00 2001 From: Laura DeCicco Date: Tue, 21 Jun 2016 16:40:18 -0500 Subject: [PATCH 3/4] Update version --- DESCRIPTION | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 9ec9c49a..bd451123 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: EGRET Type: Package Title: Exploration and Graphics for RivEr Trends (EGRET) -Version: 2.5.2 +Version: 2.5.3 Authors@R: c( person("Robert", "Hirsch", role = c("aut"), email = "rhirsch@usgs.gov"), person("Laura", "DeCicco", role = c("aut","cre"), @@ -10,7 +10,7 @@ Description: Statistics and graphics for streamflow history, water quality trends, and the statistical modeling algorithm: Weighted Regressions on Time, Discharge, and Season (WRTDS). License: CC0 -Date: 2016-03-07 +Date: 2016-06-21 Depends: R (>= 3.0) Imports: From d0cd64cab724e236d4172304d62ea055a67f3e5a Mon Sep 17 00:00:00 2001 From: Laura DeCicco Date: Tue, 21 Jun 2016 16:41:54 -0500 Subject: [PATCH 4/4] Adding news --- NEWS | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/NEWS b/NEWS index 8449a60c..0a90ee4d 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,7 @@ +EGRET 2.5.3 +=========== +* Argument for randomized censored values names randomCensored + EGRET 2.5.0 =========== * Added an argument to the concentration and residual plotting functions to randomize the censored values for the plots.