From 7aa7ad3830ba5ed57d9461b1163d7e445a51bc75 Mon Sep 17 00:00:00 2001
From: Aaron Weiskittel <aaron.weiskittel@maine.edu>
Date: Sat, 25 Apr 2026 10:32:29 -0400
Subject: [PATCH] Constrained refit with biological sign bounds

Re-fit DG and HG equations using nlsLM with sign-bounded coefficients
on the full 8.22M-row CONUS remeasurement panel + ClimateNA EMT/TD
join. The previous unconstrained fits had biologically wrong-signed
coefficients in 70% of species (HG_B4 alone in 56% of HG fits). The
constrained refit pins those at the biological minimum while leaving
free coefficients to adjust.

Coverage: 113 DG species (was 84), 109 HG species (was 96).
Validation against df_dg_res / df_hg_res on Douglas-fir: cor 0.9998
(DG) and 0.9965 (HG); RMSE 0.245 in / 3.64 ft.

Reproducible 47-min Cardinal SLURM job; scripts/refit_constrained.R
adds the bounded refit driver with a synthetic-data smoke test.

Co-authored-by: Greg Johnson <greg@nosnhoj.org>
Co-authored-by: David Marshall <davidkathymarshall@comcast.net>
---
 README.md                   |  71 ++++++
 rds/dg_parms.RDS            | Bin 7766 -> 9373 bytes
 rds/hg_parms.RDS            | Bin 9809 -> 9611 bytes
 scripts/refit_constrained.R | 488 ++++++++++++++++++++++++++++++++++++
 4 files changed, 559 insertions(+)
 create mode 100644 scripts/refit_constrained.R

diff --git a/README.md b/README.md
index 18dc304..a297f3c 100644
--- a/README.md
+++ b/README.md
@@ -32,3 +32,74 @@ Contacts:
 
 [^1]: [Introduction to the Satellite Embedding Dataset.](developers.google.com/earth-engine/tutorials/community/satellite-embedding-01-introduction)
 
+
+## Constrained refit (April 2026)
+
+`rds/dg_parms.RDS` and `rds/hg_parms.RDS` were refit with biological sign
+constraints on April 25, 2026. The previous unconstrained `nlsLM` fits had
+biologically wrong-signed coefficients in 70% of species across the 99-species
+union of the DG and HG fits. The dominant issue was `HG_B4` (CCFL competition
+penalty), which was negative in 56% of HG fits — implying that crown
+competition increases height growth, opposite to physiology.
+
+### Method
+
+The fits use the same `est_dg` and `est_hg` integrated equations as before but
+with `nlsLM` lower / upper bounds on every coefficient. Bounds are tuned to
+allow biological flexibility while blocking obviously wrong signs:
+
+```
+DG bounds  B0 [-10, 5]    B1 [-2, 0]    B2 [-10, 0]   B3 [0.5, 5]
+           B4 [0.1, 2]    B5 [-0.01, 0.01]   B6 [-0.1, 0.5]
+
+HG bounds  B1 [0.001, 0.5]  B2 [0.5, 5]   B3 [0, 5]    B4 [0, 0.05]
+           B5 [-0.001, 0.01]  B6 [-1, 3]   B7 [-0.5, 1]  B8 [0, 5]
+```
+
+### Training data
+
+The refit used a CONUS-wide remeasurement panel (8.22 M tree-level
+remeasurement pairs from FIA), filtered to live trees with REMPER 5–50 yr
+and complete predictors. After filtering: 5.41 M DG rows and 4.50 M HG rows.
+EMT and TD covariates were extracted from the ClimateNA Normal_1991_2020
+rasters at the 138,796 unique plot lat/lons. The reproducible pipeline
+(`scripts/refit_constrained.R` plus a Cardinal SLURM submission script) ran
+in 47 minutes on a single OSC Cardinal node.
+
+### Coverage
+
+| Equation | Species fit | Convergence | New species vs prior |
+| --- | ---: | --- | ---: |
+| DG | 113 | 113 of 113 | +29 |
+| HG | 109 | 109 of 111 | +13 |
+
+### Bound bindings
+
+91 individual bound bindings across 73 species. The most common are
+`HG_B4` pinned at 0 (54 species), `HG_B3` pinned at 0 (21 species), and
+`HG_B8` pinned at 0 (9 species). All correspond to coefficients that were
+biologically wrong-signed in the unconstrained fit. A full per-species
+table is in `pdfs/bound_violations_constrained.csv` (added in this PR).
+
+### Validation
+
+On a 3,000-row sample of `df_dg_res.RDS` / `df_hg_res.RDS` (Douglas-fir
+residuals, SPCD 202), the constrained predictions track the previous
+unconstrained predictions:
+
+| Equation | RMSE | Bias | Correlation |
+| --- | ---: | ---: | ---: |
+| DG | 0.245 in | +0.004 | 0.9998 |
+| HG | 3.64 ft | -0.535 | 0.9965 |
+
+The HG bias of -0.5 ft reflects the expected downward shift where `HG_B4`
+was previously negative; the previous fit was inflating predicted growth
+through a wrong-signed competition term.
+
+### Reproducibility
+
+`scripts/refit_constrained.R` (added in this PR) is a self-contained R
+script with `simulate_training()` and `run_smoke_test()` so it can be
+exercised on synthetic data without the private training panel. The
+SLURM job script and per-species log from the production run live in the
+sibling `fvs-modern/calibration/refit_constrained/` directory on Cardinal.
diff --git a/rds/dg_parms.RDS b/rds/dg_parms.RDS
index 8cf0d7124a1107ffe50024e8f6ade298ffe4b63e..5bdd2391fa73b685d666eece392c6b4911f71ca7 100644
GIT binary patch
literal 9373
zcmV;OBx2hiiwFP!000001MPZgG?ri6_r6@0c^;Ceh%%Kjk7aK(X^>JX63WmZ5t673
zNis!BQs$zRqzq-)XNHK(^E}T=rpWlL|J|kgTF-}fy<eVpSnK{^{f>R?WB47%9?tVx
z&J!#E05YIZ03g!?008~5qq}8P0U+}N0N4Ql8~^}&B>vAF07#OGTd62VMFlD<Qc=(S
z59(9L`iP2D{b8#9Fjaq;_=88N<D*nOzJgT0<5as7RQ*XR8m%DJ&WLJfOwD6V)f-dw
zCRF`ts^4kqyeV~TN*$YtP|=KvEBfd^c%Ev1fjWPYs=r9hV@vh7rP|s3X$Z2jqarn)
z!+k2IP?75IP)WtkKj`%H4_>9lbC#vz%9tAW8Wmlr`CX`XE-QL!UKeU!H)<X?YMz@^
z9`1Xo_?U`=RGj^T9@IX1Q1g3Uqarn)mjD&dQ8Dum-e#gARqsvhi#IiocOn(3JiMuS
z@2FC-@elg^*^tPGs`sJB@u9AR50!^67Zs`YzSMsBQtR`juAeWJk1w@PzSREsQsenj
z`TJ7W$9E+@m5(pAFZZZ@xks(@9<>kmsQK?v^ZQfB0o3^bYW%<xRHX6>{L@$JnY0Js
zkPIM8F9G0W41mRc0A7^<z`qXw;p+eho&?|p34jgT02ojMU<)4r&rJakw*f$m41ixG
z02fvRk}@|St?2+@brArkjR3e&1_0kHK!T$H<c0y#>SO@;cK~222}p+D0T@yMK#BrL
z3Zj4^c?tmY2LK$<1$42B0IcEyK&A%(6bk4z-2z~HI3R7e0^q?$02Ist*bxE1RuW*)
zJ_LaDFd#n`2SBhLfb~KEj2r>L>oFj0HwQp>8Gs!CfWR04E+zo5lL3GbUjU@S0Z?B8
zKs^Njfe}Em_W+>yHvnCO02DX_u&NjU4QBu}2Lat0Q$X7L4FH*200?aZV7oQ|d?Ely
zlmk$?3xJR$0DN8ma4i)8VJ-kvI{@&V0(7O#0BpPhNGuhAPBaVvv9kbNIt+lcJ^+?R
z00{X2u+tR)z7zlq0RX3G0T5IKK(q*ec6$Ie-v%TW4?wyY0l;7b0QIE+NVov7$p`=~
zMgZPG09@k%NDT!*@y*KJibt*X2KVz{!OOQYj&s8=c*$MMX=y_#0itj1o3_uvzbf^d
z7z|zTROH&MKFxQyz<bk+^t=27)0MnfJt1#A@9C^G#`~FIS|5=A$)gBQubGfb2k!|c
zcVt(hx|!h7Yg#{b*%(h<-@T?VX9zEPrKmZa6vcHQ;(fJH6+v<rE9u?&jbIEI2u;+F
zB$zHpdGL^z@zS_+`pl(u1Vi2qtF_*e1cS!G*@QSwf-D}rW5?$_JVP!tix3yWxiXV4
z-xTHHK2JSc#ixFFvD(t%IB74Qu<i)i6Z9K@i#)A&by1k0*h)U4*XqF|Hs{_OCbHs%
zBP?~p%*XL?RVSyVR08g=99B7<BTtatHD%h#_TlNdwgZDlZ1MNsxeH8n-wFDA2R>?P
z2ol^8Er)){Dihps(rt+>&+)?h`+gs;ml2GLlw`i*96S<V!k$W05_F#C+F<)7g4D1x
zmFqVk*fTVjwj0P3On0ph6#-|0@zc5BL&`S@#`TY%c@P_M6x&&Quf=qN&b%^kUvnD4
ziDMcsZKWsZ6>Wd)G^Hol(-@scN}>obX7LS093~i4$xefMc?7Fkv}=3w5W%KDxAk{S
z4Z+M`pS}LzB*CQ1_(W0qCSDe)S8d3?L9ja}%?5j05Tw#&)yj(p3Gi^o9KIGoFm-Pp
zmirJ+FnilfAFbR?kQ)obe;j5Z$m=ffh^VI#^gp&NAKw*EkR7E?_%!|^*!N1_H~|L;
zhNA~;u4GRTbfFDy<T^=$o*_${eSVH$SK8L%F0qLqgX8<y7IFyA?WgE@t^5e`;R)kH
zHx+_za^PjQ!XAPt<qq==n`wgn-u|O0{T~Q+U9;FrPa#3sxY!|-f1050>PZdgFCo~%
z?+N(NrV(tt_O1Q^5y0NO_p9vy!6Ge_ebM{}UOYsaKiD8kFl0Y-E%g#7=sq^D;oC?i
z=<iwjPMZ%BBtf#D*kvX>cf)0Tm8~?vqdKZ;b&`u<I@Gr@<GmFDevun{fd*bm-9Ku6
z%NNfaGMX3g>LKWT+>=dnjPdfLEvg4QhVk#$b~UV(h@cc24E$8ABIuQ`we<8oz)Q~Q
zBSxI{1ic=YLM^i&!FcEW0Yc$0!A!1ebucX>7@oq~OP8+`bZ_>}p7Eb2fD!B2@~agD
zv)%&pP{=ugVt&``tqa@V^Vt8Geq`SpN8FNfsYL#D8}1ye-@@E{6?YYVoBY6Ng9q<#
zR*`sp5_fkcs4_Q=;!bCsW0$qKaf|w?Yi||;absas-nF59cw)~3dY!Q(+;#o^Uj2YF
zJaExs{c$*jI}E>+uoa5o{!zj)f71rsXT~Po2Q2U?e#SXhvJVeNoIh^EV}pC1r?vNI
zFyS^PrVxr_B<^N7`#w)6758*$yN0{1!IK}yjlAr}agXozN1yb#@o>TA`y#I|;(o4(
zX8ZUaJS6m_ra>wLkFwK!xa${-NAaVl5~qjpnAs=#GdAY9zr?%E{sF)(F^@WM<pJDj
zTz^p8xDxj`9=Y_>LJjxE%kMT>c!~QZoVWRtpW*KB&oU#YCUHZ?F$KfVv$!qlZ9shu
z7cSe~-1mDr8`lQ-XBlq|#2p?1qsLecamW2C`sgiExMy0IX_fW`+#|So`34&wZp8EE
zZ?hF~lZ{T7fw~2*x5_$O7F><X-*n!KpPj=EQM*ihibC-K+3%gOV+3xVpud!VgAcd4
z2;AUl3B+A$rx+iIC*v-{bnUi<XxvyT6qTbUh%1I6YnVU<E?iqvF5o_kzodi?lM`R#
z8lS!mOXqjvnm%3m^0$n*^Xa8Ww<S_=Q?#N+_PH$lqqAB42=9K}-CbZ!c%8;|rDRVe
zBZV7!@=xi%sK9j><WorQ&u~qZdY{>dH*TSe{8%Tbj(bK2v{>}qaCgD5N!6?f?$Wh%
z4c83B9kFIf{cKsdX7m0K`?`5tZ6bF@PR$%wsh6#)>vYBSZbiw`^hUT=k}WUu$uV5N
z(b(d<?J-;saj_*YjTd)s7d!XycR8*uD5uQuJLB%Pb(`YN=5gy2PFrz<aoq5^=dO@7
zAFfkgidm}5!}Y-o%Ez9Nag(1XF{P@A%YscOJTC6Tb)pov!4qltixwp?rPvKudiuJs
zt-6J4uZPc0X4&9+PuWt@SlK^sG63LWzak~M^8goL#wzOiuW(J6X>6jr8E!g1##QlI
z0@pcJUsy6&jmyW%IZ8V1aJ{ETf+SNWCK8N`Oy7UU#S?83VMmm4%_SQ)SNg}eI%whq
zug_jwUw3A<vdR>9z3Wu;7Inns1G_(sut(q~wbY2Omfvuz>ovXhqt9`rgQh)6c^Wsg
zh)Z!5TmC(djVG5=o*cZ5(+H1`qn(MkP?XJ6^6Vl`U0CZ}RM3xW_*xGPl-|K#ZwYD{
z=J4Zc0rlqS146j!`S!dkV|BRPBa*cKxdCoEV&gF0cncS($JlJREQX85)oXm$aN)dP
zx@|0Nzj2=E?<S7&Yq&6r|E1EFd0h6^{ifp9wK&_kj#-_f8*}Xm6FNWp6UPu+EjTks
zIP+;S|0RAVoJ?N0T3&h(_Z}Q5VHwK9sVd1PzbQPpUh%BhkR}s;BfQ^XU70z4?^ten
z?MeV9B(F%*E4$)%vSEy^KNE1i>*}ygs%<#pXv~YzE%G?Hs#9&(>m4{zolD|&Ob&h>
zr*f*&m<#8qW@WtJcpo?Y;$J+wa2jVkFk_24B8aoe-Y@lg-Epd^5_@vN40g5snzpo2
z7`q2M98>u*g!AjirE8DH;TNSsqved5xOC(V-x?wa7feP)oXV5JWt0cgly(tZ&N{~e
z7w_U!_wp}prqlQXKpPo%WZ@)%>t+FghcWT_-qC^5QC$6CDo@8_9gd9awqC9($Hi_1
zZ?b+o#UDu9k~JUp;BN2nVt)@CocJg)t-lB3m?P-IMrMB;t*)O=-_U{+IfIEJ-zuCu
zW;FMZXA3TP&Kr>^+=oAX)$q<N;lwQ*lku<Kt;J2(?X^?wIdCdR<5LaaeK`K-!NKBW
zVVurWW~-aKh|?^CnIet|;3_9zimV23t>|F8(79-wYMQ_sXHkmN3*r6*19cq7q*AwK
zIsxZ%B^3{C_=w*pWD*bWRO5F&uxI%8CS0}M&QnE13*S4l?LbmmCH`zCPj=uwfnO}O
z+-5oc5m(V`RmJ2l;rMf-`Z5N>IMG2nw2!S1=XK8d?h+Zt4@1!<{w=0BHDKWE@og-)
zB5m_htouuxYP&z5lyVtAZEihze4PpI74v$fdT9*DWGiw8Wv1Zw)}phk>b&v0;bRwM
z-R|SKnZeuqt_!%<zagsi@;ojzI$<d&cM<0<&9|l$W#A7aUt%BI9bChGE#l<t3H*8u
z<7`~!Tb%k#N^5klD9*4xY<ov14F_DfH!wxF6Gw6?hwO=yz^~*3Ht!*%@Qa$mYw(>&
zT%@@whGYK-u7PxynZ}K9Ripjgml;hsehurmUd|G3aouiuG|wLs4PhTYvzy{d6Tz#y
zoCH_SBb*PN)hqH9M5L)A**dd32<S;XYTROpKs%pK=)xeJGg1z34>W`GjfGY735sy?
z^7=g!vRMd#V*Lh{*Kj6q*U!518gN<o7u%YuJh-qcpFy7A2$9VUR2rXTB4+KT2i$Dt
za0Zmk?0=#J7YAlGcdm7YvliU~fmP4o?2yF3q?|e;Rk|B*6AeOSkw-7NeqcDSvcsrQ
zJ_F9ZE?SNVu7SUIV^`gW-Eh&5mu{=gVK~bN*M8Bs2PYaD-hcyM@Ym7G^Uvko;cqc6
zbs};e{u&cI+Vjd6{@S{&UpP$@{tD3o+wPho5K=EaMMj9UH(A{?>Lr}fstvT1cZR=g
zKgKCE&A=&@9I33<7jTLGaLDY7TDT|`Z%x{G7?CODz%dCs#7F%t$hO=TFJRq?NSCt&
z>n3~=C~0nbUN#5kxQ1jNUl@W5JJY{bz;}qu5;Mljc?~giiFEJEHb+dY(;S}u+=!d3
zmictQHDcUYoTm7;15R8jTp9_^MdWPD;YNxYB46bxc`rQ+7u1?5S)%-iZmYZa8_~6h
zBoll4`sFsn6qBeYrEHG4zp*EeNtPj=?;Xp1LGFl>GfyWXVSy;(Jd<K&6NvPR9Y>~I
zL`*vtW>N=S;IE@hc24Ui;gT%-&d|#T5asDYUh^vt#I#^nEZQ-NSUe5m&mHwbjQC=2
zn_nPekm@l!I!Ff>&S^f}I#7kkWMxnK#K(vfG*TVFaRz}VBc}&dLWur`&K%RycEq9m
zV~1)T6QZB3<~@4)K4KPceQWbr0xnx|(Fu#YBf9Q)tLJ<^Auh2OZCy3vaKW2lbDVAl
zqCX(AQ!kPckvGYgt`Dz5^y*?)G|q=1<`<@q_rJP>SR#CWZc$4?^k&KGhZp-1lQoBH
z+{M#~!XLIw1nVIBu=`?<_<0c*_w}ax*?x${_e!oB*G~jOvlC6Ttq^FqC3LvF1uip&
z<*y^u5DR+i(WzC3Sfgs|SYLfZEWsa@+T<_8-|80MCOo$y#`M93mr-JfiEG0dyOxJ=
zIh6O9&h5K@*Pr9{I)Odu4-iQHHM(luOT-YPzJHgD7-AX>eBE=58Ifvw3@tL-5#zT_
zHwyDVBl_)zABNj6!sXkidaDc_5xdFO9)`JCMA7?IHIo$i_dGV-Zq-oI+Xj2&IyZEv
z-GnV&KkD@KWnj<fD2KZG7uYCp(%^!O0c_rNr=sh{0_<|zz5V<{5!jLt{^Nk3DeTy4
zJ$WE`9c<&crY+!T4ttz*+T_k9!=C+Tm$^oqVEI(`t{q>*VRzH5NY2!ku$^bl+G9&C
zY$rXO{1_4odkG<p7jNFc`DDvI_I^BY=*)+3-9{<cd)<5ZLUJ4&)SrBjw>%8%<70Y?
ziyYz5i?PLm$^h7$TO^Y?CId&r-{)8)Ux$Nw$BjirEnp9@`rx2`5)KI+CTqNT4|~1^
zg&%u71HYHMMmuPU!=I#+UU7^maD@58?(a){aQMKlSLW{fVBf2K4kEtZuxaVTTC=@$
zuvbBO!I%;TyKAKtJB90EpW^XmmW<PIEb?2laT6CD@iHmuo1ccA6UF7H1lZtU*rBd-
zya;xrGjyrlD}-%LqC!8LJ7HVfti9wzXW05~Krcwq0ro51Am%*;VUJkp4%VPAaEND1
z@)j-&*kfNq={}wYTN=GztO@agBjkjv96tltW!ENRuDAjASn*pORdItg3hD<FDeiDk
z<eIB(jwKxNS|hYvy$=r4NzNMtU4g@o3J28R2*Sa1y7Bl^U9dZmd3WlL->`e4IkPPA
zDr}DB`}Ep<4mMu0QWkWPfc>f`HhfLbghQ+wf0miPfnBa<p{}uIu#a*)<4a2ktPMD(
z&Tk+I8;tV-<1<;<8?IUG*J23!=?;IJnU94dnY`gltX{A)O=W;^b%DJ`(^(tr17RPE
zd;fl$D{S;1&E;Gh1;+whSbjb<h28lkq3$f&u<xTh!E&AtwrL65`81TlfwbR4H=Q@b
zLG&g$X^s!J3v%yKxp@bU*sZHt3bupYX=72ryMM#>P5v!n0u<O6I&(@}sStioN;^CA
zB^j1E?Tt>^HV%h%E9&U_%3z;A<wE;f9XNRCz3K**3^;n<n)wLlH0)3xy)@&#6E+`7
z0`((ousc~%qjT2~?5$i|c3=Mm>=b9!86P_W2dd=1Uj0}E2iCoJo8Pnu2ei5B3l3Ys
zF=mDm5qnuU?4pvP?{NYS7lm!Nt?z|H7h~5R&D4h-Rw71iD*mvx<nwo~LN3@u>@x^g
zl7&NRUoTgz*N08FyPy`S7j}806aA$Xu)W>-NV=cn-}hBkK5>$dOQGge;vrsc;R~8e
zr=tyT?q<~dd8KQf6~{n|`|~~>@w>ULvq~;&p2}ePnRlP9AEiQxqIXYjMh#;2@q|V5
zGDpn{dD2;v<NQ$Er8raAvKx!^JxdN=9fT!ce_~vFDl~W2cOZ46_ddv5UuAV-^E%Bg
z`GZo)=K7elBj2D&`zB_;Co%cvhzz9gI9q;P8w;6B!)>{swC19?t9df>6Uf0MPx&6Z
z1s@4%lbbmzff)^+>MKJF%)ct;!HHe3@oMmNgVDky6twXyN1Sa~Lc{v5<cF=AJr_<=
zR2~oHw%$1xhHb-h3ln+S)-CbmW=AfD{<wJ$%ZT(zN8S(BETPDcQ$Q0GU=31Q4y?m0
z=M3zuit6!J?Wx0;JfC5a?@4bjHVfk2`Wl-&oip*OTX#bpUbW-ZP+#q9kOSs979_36
zOu;N_59-DS4(68f1*O?`?8U1m%trGSLLkRte*XzRV@w)BHv6`m#9}^)ckPC3Aa|)w
z`4a!{+!9_ZPR0YAxwVR$@;wG;Az%LZjC9Uj&FNz)qOY`jF_%fN*`vaBnuF&xDptqr
z!hpo#@Ui7C7BDQd61Xpr+xmXebVuHL%|+fLLDkOkn3K-!4x{@G&Dwma4<+n-;D+1+
zyOgm4%(u>xe$SPk7>txy^D(RDP8DB2C9~l+WEhQoC$l#nlNCtYL#o&1{vP0E+8xmg
zx!LGX`L@Pt);-*5f8(`X?!q~RBDqo#EdNoLiQ%0t7H^j8UZ=GLIV&vJvoWf}t*xBg
zQ}GN?@@gELP?Qc7j83U%vABjMe0GgH@iJhtt$(}GMLNjfVtwbz$>WeCY=l0IxeBi~
z5j?evJ4>^Vt#U{;<`HDAl$G~hvs<%sc#rr_Ee*_QaF;UIRD)$$zwD{jk<gr3#&4t9
zV=+r|gLS^(x!j5C!O9K_1m;>OOzRpM!16CXD8D%qi6siy{c;m_;SF0CNm2#Yknzyr
z<96w;kUg1w-Og3ckjr;&O<(gIlrWc&t8Lf=<*V**a8$6tk_vZ3-a0TtN`RC{gvA8j
zTsQbctlS)P>!>Qa?aRxZ8QAY|eh)X~%^2)!*r|uPWlk1ZepSF+*5T!Q2|v6sGd}BF
zuncBXxZqNA`(*B~Hj?Pv+s}{*7v9v(XV$Fz7-JCddlGZo?@jF8y!`k17~h@}v7K;-
z5l+z~r(@5;D6xKns#U}AnV{dcswj4tASR|?{)rg|zbni#OcTIy1)Z+h&mau>BpT6V
zT7*5_y{eX3D)0jtKFjV#Ul{c(a<ql@5{?;;dhRyY0Pi?;mb_@&0t27gpOXkHfl&jh
zN9PwN;OjQ9dMrp7yKW0JsY*W!<EpN<a<EuHTh&K$r~L|G?D}{w&nulcASfv{wdpwy
z?I=uAJ!_6%F1tr;7Hxz6Q#IS8A~5vr>q<-$cER$*p$k6_@jx%$HkEewv)Jvdp-(-j
z9Y(XcpD$;&#xZpnVGIGEVe;6aLDO+p==ZHYr`LE9Klm|LHRN&?-jCiGdAC0a-+El=
z#`Q4<MynmP?eN)*PvsWgf0Zka@3mAZZ$5Ycn-1_qtzw68uV&r8Yo}b{!+3>DDWb>l
z?O{ru?yxP4^yn@s%(a64iI+4YfBN8k>&yDCKU#;qdFzV>wJt-ym!$LcNuQxNCq=Jq
zFdjaP?Kj&J)B!^iW;di28^9=1Dhu0UGY&q`+9t5w5kI=~bpPq427G0sQ92XPJNUTa
z*|l2<_0VD)M|K>y1vc85*wdfW0`FWq6TH9T7BtUmS$FfS4R()}kP;T|hmUmS`0Sn;
zVBakcRc35!v2lPe<Ba`t81ihoBIDT*zG0gccmMlocxUizdk@JR#%D2342C<v2R#+e
z17)A!4g2%=jjDBE=og7JIkf=%RNQL*PGvlHze*CB2{Od~KkN?O%;toTUp|x@p(H~y
znYW+5YyN;?vZ<4}j}FJw1OvS@;yBvBg;XD&1ViM4C|hz*z>v%ZxyUn*p%49+>!VR8
zaHKQRHfh=m16ds<^HN)3NcEYRji<u!uGc8+n%YC`d@inJxIYI!<H$*5+;bB;IjpMX
zY;%NRJ|;e<LQU}g>zm_xtz0maSF<p3=m0bV^l{&ULU53V6JKP=F}zODf#tVlDL$Os
zBgnq33O;C*2QrTn;e)uXgWdzSFo~~?&nu@LhU|>vXgBMIA%2I$o=Y^rN3T`BR`hMc
zj{8$ArTE9;<Hp?F`9X2$ne&lha99>T&eIks?fniJ6s@ug&aQ>`L|?f&4Qb;itG85#
zoNL4W$)h*ee)GU^vD?>+e#}Gr+cDo;W)<+Y+Vl4=_rAk7SeWnfS04QPzvPaJyyaDu
zwP<>Xl``~%8U37agSs4#(9m_~yRK68Xq@3-qi-n(Ds<bn^LsEtP3DT1{J0HJ>C=p>
z`@HPY#CR6nSj#(98>+K@M7$i0^PHZMk*Y;aKUKC))tI0f$tnsFumM%Zlu4~>bw||^
z-Crj~Inj8u;F&tElc<GpKPh4&m}ruI*BfqJiyETtCpq7gLuvH6tY%AvsHprRJ#pd%
z>N{P#b^G~Eh{$f2T{QJUvuz1!(%To%Pc0*{o;NI7B;nZn-C^g@SE&crXT@_-?}c<9
zZs4!Q#NoH6{CYX+ZQJ(nlt(lA?Btl~<>HQ(TMj2Z{3N5bEw9LVBZDt$k590Rj9G>1
zc>IqS)F0JabCqvH_mm;3e#dxa!*Dn1+NTqF%zFsc1W(Gg(z6pawm-g>{m4MwwUWcC
zB7x|Kdc;Pp$!;_yl4|0b$BrtQf3m)r{)#$*clpgj7f{NN-}sqehL%ulhw{tZ64Ycg
zo%=eAti^QsV^qUEQ&gpCa=*=cj(Fg?>SjJ~9vUi0Iw?GS9ktmE#+h#VfhvW+|5PyE
zjY`~H4rJ}?Mm3GsEOsBqsPnyWy|CI#RL@~rk)Ak>Vxl+yV6^Q-t=lf<IjJ<G`u23e
z7YgC1NwXx;wY!5zi2uT*Qj~;}M|?AJZW5YqOp@j?4nsdrTnM+Rh(%fJ9f<I^V`y?Z
zlhbh1W>k<~BqVJ79A%WOKisRQjoKnzo|!}jYq7<eH$T<3*JAW^x9apIp;~qhQ}<v!
zRHWc&UbZb2)eQ!W5os(a-|J!Bw{<h9?pIyEahwk7Ih1=|^kO6$${xs~<5NK%!llg`
zy$Y!MRV(?n=X+H1?8aPQMjgsWBhSb=?5KO~jgvct<g^$Jb_54%)uX18w1*~aHE6b!
zL)kTVREy4U=V*x&11eTF?e7;BMZLGN!_~-l=xcZy>vWv}>U$MVVo0Jxm0z0dgH^tv
z&hO*rUMumToZ#;MZYKvcP`JJ<memT?+<v809WH^=q^kLaC#O-7+sOIjy3=Uca`RmD
z+DX*JeOz%<ffV}wy4|bb@KY_IJ0Q(0Ka0lC-sbI#dW$B~Z8BM-(@<UUZK0D%LFlWM
z{ei2yo}%KL`>$@deT0S%s?H7>m7p$JC*5iz7Br>tg4gfq5ma~6%2aCLpAQ56`S`;g
z06_ZBANc>h2Y>F?455p`VY3scTDQ}G_DBw@(mS%C9_oR<SbSiih~}f>7$u|quI(tM
z_m#E1f(j}Ull40QJ<&_1-t-=W>!@6E@zAzBJCw0`zJSs10Q$mi$WM>wQI^P+n5bMi
zlw2^Vd{gZ>dTA#9U4uUby%P_=e5s)Tm6wv9T+K2??-V536z+GRSKBMIJ$}AH89du<
zj_^H3_j=`}zb#m!Sl#O74i#UNosiU3X?6oWHa9ueZIys3w|^Y(Xlp}-bt1QO3)4^=
zQ`TYzs~<|7jyFtrIfQcOCpTnknxTQW9y3E*{7}|tBwxI=H+pnUr`Ard5H-ZmZDsaD
zT3n5OOgap5=mYykPAiiv^lIYd>Fk8PXg<;Zr@8cJRA}k3vH8P(RC=u2vObavy*%6D
z5a1t!vPS1)_DtoY!o892=ZAe!e$Q?@!$Lt+qcpuE#(x)zpOACA$Nd;(e<i0&=|`g=
z0im*Yr?t>{^0PH1lTql?%Tspka`LFdgSX!>dp-JaWp2%4n?7<LGpz{17ZBlF@mzbA
z6iTg^ocYZ{Lg{J8_RU&LqL}okEgGR;&?}W+@yB#sQQCpscf>li(Ce_^oz`8KQOJ7b
zjT!@VDAhCdWC^1#O8efgvD!HXy<fY(z`(W^-Eoj}d#L7!ykFD|4?SChN}JDqH>(pz
zg-=-(R&831-c;plkKEmc%9J~+MZSMTG0ZqLMOzAmem;-%O<Yi>+|-j|I(pO>$Jubk
z;s+W?x4Y}`$rHUL@2VDR<3uqY383@hEtI^Ki)81u3Pl@AXP+K9g3`o?*RWjIK_4H-
zRJ0rNpwHJTwsPp)L(v;#<fLRVdZlJvdosr#rN5I*mRt7=J=x~uOBA)C4BnyZ3n4F2
zs|_4)_W~%7bHIgZ(HvFn9P#xNBB6$q>_*A<JQR(}Dw^C2P<DEm+G;aq6de|`nC<oy
z5tc2R0{xz%)N@svR7I>&_C)9vnMac-bn~3v<mh$uWxABDS2P2?)K!)mKiGx*oYJLD
zrfgBvnY8)dxCqoCW*r!DO%Hu^o%}TtBZl(nEEmc&Wl_C;xCh7WO!URm*KOn7{phRV
zo^X-ba`a(t_&{oVG>Rby&OB~fMA4Fg98=<ch|rtq-?d8{rGrsB$D$BaNjDo=@}d>J
z_1f(aKk?7UMqK~8TPc7U0KjNrX=UzY<M=;6w6ZaGG(DsEzc;k_U;0}7U*c@>KgP=H
z|H%R^{xfv?Pd6z4*$m)6XVCKh@XY`ybpLvNr2I=wZ)<+Z(t-NyVy7x84t8fPs4~fx
zD$r@m{*~nZO7edtg}>6utw*P=^jA`*N@Oiv9jeH1-a+TG?bVfuBL)Ukk;TE${EDNg
zlY^zje`Cvh<+AfuTE_M(jJ0Kz{-)5DUExB}mQ`Hg@Q=yg3gp!OtyNobq-e{ks;*cl
z|E)s#ZxzZkmh%6er-`MYNE2(NH~+{qo(eSfiYxgj+Ol#qOyyKoijk9DY4$$_tb`}a
z$*x4C*Orx&T`BAzSDLg6G`?~xsmOA2G;EYs%*k@{G@i7w3Ju}Eu5^F514Ua_X(bO?
zPGP0@WNlfRg4AdVRiRlP<&`L8?Y|XKv}ONxpQ0@*M`Nx)6H<kyIAvP8e=E|a?VSqE
z{*b33ra&Xplq65%sz_rmN5f1>g@%Cw4NFxTUlkgART>)=8fvOEaa3qzRT^$e|E5o~
z9`dxZ5>0rTQq*WeD$~e+S4~@%b|1>oFr%qfjiw9AO0>I=#zB>~g|y|;<W;0?E$xA!
zsfG4XDbnl_T60Aj&MW3*ZCRRXRcSh+NW(yd#zv9$Dv(<Vqe{EZYBVNFwD$?^-l1hi
z6JDM+Esd`{?UAKz8cj@j+EQqnOk<-&d;8Gdd$jwJ_H@Yp_b0$v8*_)1`*cOgVqxxR
XzSZiA`K5op$y@p_rou>HGc*7I$|GVB

literal 7766
zcmZwJRag{&f(2j^q`Nz$8)=5_l1^!*dw^l+ZjcZO>1ODmyAcLN7(lwayPLh=KHlB)
z^uM2nb7<mF-~DfpPcrSrW|ro9gLXx9Wq9R0?0F#b#%w85{_L_;l-ZPf(UuNC`!u=8
z=)Yvi;mEReI3u5be6mg}S7!Vu*5Swf(*Cx$bG7Mn5mYDPchVGfw&+;hrIGiHjgF7h
zClB=Sz&8;yWc0p0TdJsvaP|P&689-buY+mQ8j$+xe$Wb6*J4_jk}UnHL)~hNZMFW^
zf#?ziG{+f5uHviB?v;UACf2Tti@?r}gv4qRPgYwaTE#<`IK7YQW`(xgj?oMlT~O-O
z5aKwR^M(Q>mU9OJe3m!6TJ^tX@y)Ut!uzXHTxQO_xo5}Rj6)f5me^`J`UDY|?*<%Q
zEYSMh@s{ii<?0MA@L+<(xfF@=1aN`A7fp_%y7R3-oMyyj>Y*10u{S<qIw!w`9kJ$!
zuV}vu!}iws+}Azu6j}%C){PS6Dek;1GeijtToyD#{WP%q3dj|K<ZcAncml$|jnT$W
zNYE12kwaVCNG)Yh6Zi~Ju@5SRtyq!VPT2bDSkl*%>KIT<wvb?A{tDrD&cAIj2Lf@B
zaOFPTmV8R5{Z@3I=^rW__Cj5PAHIYhn&icd%FG&yTEQ6IVUM&C0k1Dlp*>v>D6D;p
zgra@CL6;T_Qo)6k0#yL0-#(#8Qw}(ayDi^s)CG92jR~7~T-?VLA=RzZN^81xI+Ni2
zPJTeDv&NF1&%z7;Bmof_=M+ez)$hRnEqGj%%qH!W{@d@PSb7|-zF9xr2tzbcYd7%>
zZa84PfqJ}BOd^h9{OHdMQXTqS-m*X^G0X&gF8&(>HLR&D8ZLk`7vnFuN>2^bJ3ed0
zf8R<}z{h#;8L%ST_hTsk)shfcXv8k8NqO)OInY*!E!FyHlc(#j6=N`Ws8%muSX7(r
zjWi<L)QMHSr<!2#FQ~3WCp*>IPW|l{5+^IjHa4uEJ1#hjDqddCLi-p35#C}1TSdu<
zm4N{N)E0gmC_z3$vqmidb;o7um=QgKn$3kq@b58W-(ST3K&wA2vv5=pjPQ&uV8e@b
zSOu+CITV=&G7o=)N}p{?_^O@D6Fvh{oru8TB4J^+S`*50(C+nTag+*zvDkji+G6eM
z<5F9rDOUYkbp3g5=EMV35pBvk*))DghW&HGXItWt!3Et|KF8f(Jp~u~;~bF~KDUy%
zwu;Ool^Oo8&QAg{Z`=b=|5!D8U}GT0R-<WO!l9I{(U_FYxMGRM8J$ZQUwCyZTUAE0
z%jeP?{*~Dn73+2Fx+n3C9?WOSpK^foGKw45<1^hKi2TgP<YBkQ+OHx@&@s7eQ((+X
zwsgg(Pqw9_Z~8d>K_s)T7bMRjJ1zGyZz{}#3mG7Oe?QXZH!f;vz~Eqv3NuLt<5=3D
zRsC^b2lqE)(q>_QLgXzka>?777^Doo93;_FdM^D)##CM^4{0Gvm*W#C2dap%R`$NG
zyL@y53pODnPugfGzt*m1dgmWv(29d<EJ}AnGbU4xLIyRkR~jZ~VhHAE_uKv4H11x!
z@CwlbGl~5>t=0wt-I@z^M)ZKk5R~oTFQy-Q+!>x*Htzm3c^2|iDgO+Ay?ki#M-$#E
z)^y5d8Gj?a=P(WI*krgnmle{%0}~pn#{PWiNwgaQL8dNt-cINUJS^II3<(#5Qnn?b
z3$8c><A2QMZ)nn|;$>iyS>Pq?EDmLjka^I|+Y$cws(fAH!bBBe_JOdE`gZi!IQl>d
zF{KaZIOtFNn}0JeF1eS^P>5QJ>kQwFlOuBML32SldFk?BlguXbgu*ZRl>P<V&$;1N
z1)9pG5sY1!_b~KO+nAsiqsAC6Op@NF1d-Weva=05+(81GB_h-t`XVCNLh=HiF;SRb
zy+~ixWY3m&iT&HYj1sT&;^*jc*<4LXDwt{9*<y^cxso^#TVI)piE0e_BVio{Seq-0
z1vB-hH$;iN%S2zW6{D^Wjr12!yqAx!6B&vc*GlwkkY{h)$u&%`^9l8kajA4K7->nD
zM0QD6NIetUDX++uJ67RFoEsLmYXyG%S}=$xc&7eDU?Lbiroo;1P`EoxEJrTyDQ@Ca
z04C*hyiOY1nKpNFnhj$rgu|y?*)6V?9oom<)WXNI&U~xiby-|p9?VNjtJrt|`&O<L
z!zXdDziJKaurK&8Iq&%Ede8S)>#1cMi!EgT{p*9Ak?XXo6Kc%ZeEbyyv&MZA{$*s>
zFwZ*#B1m{s$r<z1s)^qO=2~?1<cb{oG)bFdLmQ${j4gg`tGOWO@>ae(F4!DBsa^>K
z-!1s7QsiDA2lnj}Hj0>u<VQvB4X<>;tl_(uWjp~^eWS#;3Zmld$66;kl^SAsRrAj?
z%-6_`IJIq0!JB&@Kk@$<5(C$tBV_rFDP^v7@;BczA(VG8eRZnoSc|h>mY;w2g;5r=
z@z2X-l%j795ZE<BBSuZR6UJVw<6BpvTI}4Vp$A#(OLj7V^;I+ot>C1!amnFl%c}~T
z@S#JWx(~PZBxVB5vIJa$k?cCUliNF(gR2!yP8)q7`<9auCl{c$R-rxfV_g|HWA7zl
z)?Ua$(LOpL6V<NGhxD&aoPBsn)-6NrU=8ZXatB<-QoXW;jkvw&)!V*q(alUy(pC5^
zF_WPgXm5OJxj#l6+eviIP?<Q^;<k1o*yPfjuN10%7@{HwQ3`HrbO*F%6&ZFznNhG5
zC&;drk>}sDCyOJmtSpFmwdGWMsx5>X&!ia44;93iHXff91lq~<HPn~h9(>LLpS0_M
zCpq7DIW3c29-3e~5SSHzYRp2Jajw>~A$O9)LYMZ!KK&${LnjU((D0dHk2?mYbz-s+
z1XxM}Jx~rix^4Q>#MR4g=-iB(z*N5FTHr)a7;dVsha0<hNr6fUnB<{}v&^q``RsXj
zZog|{@bto2L#xO2FFwJm)wQmv`yF^vTDJvga`F8Sp@Y%!{*^OG+@T|*2dXVXH8bmq
zu5)1=yl9u#;U(v+fa6(JuPXOh(~r4lif=`}3><J>3YL=f^cPXMHeS+-ikdMojvL$`
z8;y@+H18-xr<NE}4Y8PJtY!jkgTyPi?BcsczQ&OM36;058l7Ac`t&Gvlvv0W?@n=t
zL^VwoP$<|xPch@%vDf_VBgTzHd+wd=Gi$|YTe*@cO`H(`E3-99K`RcjrvzzXsh~+G
z3ba~HeB7m~eK?&6L(jYOO)0%Dh^$CLq)YRg%wJqwqDs-yr58Cnf5~~$Vy`^47-dzT
zcG@d%H9H_<+r(rIvSx``rOfjBqN#+Ix_xa+O%@pIkj(PM+>)f|qkZ6+%oVdR7am)p
zt+Rruf6nD>i=3vPnI6>|XY`}CjMZL4i6%-{M?Wi+xQU{L+1L3AQ5;UlNWI4)%RP6?
z?_*3mJ-((dY(x~7%Qg8qsjCrRz(`@m2B58F3m*2f@Feb1A$M2S*q7q=BX)J}HVg$Y
zU@a>Wze?R;h;?KNRw-2PWnti7$!g*}U27?d;OH8Fc(>b17#-kTP@F9a8Ez!q-N@AD
zqCG5gBxwu&ibaBY{`C;{il5gj1U6pCH>WIb=QFJ&o^SI41>}B+0vT9fv$y;Zv*8(J
zv)8_;L(Vbs6h4TiP)W7<&kJ%-LRO?jwv4f9`{Wh+?hhhYnx#|VEN(y(i=(+U@a`GT
zas?}5$^A#qq62fWZ#ZW7W43UW<M%|LBDy(7riT~e5?<Z1qREuVyMD)PZf={u536TB
zu346)>;8<kQ1Ybpx1{Fos{Xc+P+k_}^yEl0)4OM=VprmtpWPPv8gE?I2WJw!lV?G9
z?B?|2Tg9|)JogJ+GzsAa3uSQJSH-$}m2<RZccL(osjNKE-|(rB2exQCt>a$)@%Yf<
zZd<TIaohF%2c}dEYUXoYJ~<JLvs#RLCnR<cE_iC~n7_?E!+ZqJkrEm7f1oId)XHd6
z4?3sU1HD-<)xWUXuzVUEKchugKufTMd&IUBxPyN$JmEj?Z-}hcSo;}I4shdfc6JnN
zRhO4a?-~MB)<(2X=dR@4v|2<5;PwUE-_b-L(i<i@v%(`=xYAg<s@GL`uNeNM<VMF3
z73KHVV(ZOVWpn2-`^<E2?9=g`FMt>1ba^{BkXFpGn#9cAC2v8ju_>Iw^dnjk^CU-u
zcGY$Jf{;LHd;(JY)xkm+^EE$cv~{ebRhX>H88_hWV->P>8i?CGNDhUFnL93R-%<Y9
zy&)}Rz*bYC>e~v+SYp|{(-t5hiE;m-UN!zM{^KNl4BTwDx<#wV%HPQB6EixWl7Oz>
zN^I`cQ;ywh17WJ(v^eH>CB@_6!8=f6hZdVhTfsbzbPBd*fb6JO;q(l6E^wJhcJNFA
z9)(INJEj-FbRLbxRtaF#le{OIfiE89R;0y$!-+8U)2B_uhremr7+FM9SMgv_HPuKB
zMr<zn^*`Y)^#Knbn32C}cwhXvgA@un$a}X{vbNLMpV9jTU=C(O7-M}z3cg)!0*+YC
zuK(tC|NAZE)2`m9R8w3x#_PobjK6BHI7L^&9RJ{o0m{d<on5I7Mi@Tg>U`M05V0tH
zNECuSP|_*t{<#=iqombXGqCK*eb2M0A|UGRfWhECMa-o5U8=x1_B+;hI^m;Lw&oLL
z>0xe(Z=Le3nEu$e2Dg@T!C#HGG0@r^Ef~%{_p!RYS-fM!yDzN(+h#7)OgEO8=yqWt
zrua2UTk2IgqIdMeBk7}^R(u%ll<B0B)1z1<DMSGQ-P^d^=Yk9Dfli{^H>1!=JBVyT
zL_*f>jz;iBGp~nPo?a45yNhI6ImVy{W^5vK9#LT~wVGqs2RY&_j;qzY4IN!cpG$T3
z1Z4oZ&>!Huwf<BTiM1?R*gFR!TT!gH%+D`1R&M!9OS-h{eiG%x#Bh86%NqqhiSHQL
z+IAs5Tn^FA_G_cKw9}SMwO0snHHiVuU1Mlo_ivT-Cz=-IpO!ygKO*xSX1C|-DU!K+
za(m3AiJ=AlGP&<%hJ9diZ~srFlEIs=QKB<~XA#eI*k<<Cgu~__Ppy$h-r)fIpGjMK
z_QF|DO-pksh&;2t2P34GjHFvJEgWSSVc`=pElrA8y-HB!(w$iBVLINl`-&gPLk+ns
zUb}HpW_bj@G?A$V>n|xqZM^jse?PGcbkw7b-*i5!!&5%AWjMA}*>$4~5qdhTURcTG
z&gBeYkT^<V7o}M4WckpwC`iQwO>47EOXP#@viP=tc`bH&=CF~P*QiX|h=p~=uqxIJ
zkIIe@@rXBZaSp?RChs+Dk+)Tos*>EeqBcpuS(>{IwCR_)imRs7Jkwop%M)^Sp29%i
z(e&(+q@`muamb=FO9}^oHd#9CSrl~V7^fPhTwV9=A&K?evUYx@k>9se9p!EL<wa|X
zbYH}rrX7@vxsTWRLorAj6OdSaA)=y+=t_T@PGNN)QiIbSLo~-JrD#+?s-*eJJC9dP
z(M*WwdS?*QM%d8Tp@Vi430IVYpEQ{%ubl<Ia!d*}zH;a13H3S^W#N&|=Ky^!LlO7U
zhpB^p015ef;-`IyX%Qb8*(*iGXG3y3ugvk;B*Uoxxx8vCvhCw4+gK|c1S6Dgs-=%e
z^PJB=wb^eH-cM7fy!`eTx_9!eK8sKPurKuEis<I#{Yh!Zp`$}usQTO!UaH?gzR-=w
z?Zv}ne$lb-;ldxlIa9-4!?$4KeLxf0!)*;hfP1<w@yabSB9W?rD_b}EVw+6O(0Di6
zYMOc1b5{X5h{r6z{mGk>U)}c_R+Ehsa1e;%5?zD*(&6aUEFp~S8zFhm>_ryV9rm5L
z`@u^p?`&lD3#+3P+_S+fLD6V^Wmm>|kC))wCMbr&;E7^Fr16v%ItZS4^kG1}Bn~Mh
z5-OD$r$P=RXLh8OMiFpzUeZrb4{7lz^t;6WmCtQgQD(-4c`ek@m8WEmV`eqOg~1c|
z&h>7eD)Y+zeWT${*UlOhV!H|<q(V$ox<&4KJAsgDVk5VB^fNXh6cGNxtaBq}ARS`|
zo>$X^;BP^6yvQ*t^jnT!T2UIaJ;}a|KHznB+wU~f?oilqvigeXCt*Bd{Z0v$Ym&a3
zO<UNwQOSqN*!XBB2-u}JlCRLaQSkdKtru6beK^UIVBB8)i70Tq`8|FQvv<9we24O^
zNh)tE{TIK%dr&Ts*>4%?{F@{fR{(@;zul=ib=?#J0Ifo6waCReEIhpo_~(TnKDyP`
z$;8*U-Y<q%sW!77Up_y@-=#89NXJM*tId1W(eqpR_j=%WwvXAzx7d}pAFZJ;7)8`x
z^xcc{Bl3Ex!K5@@fll;T7pF-}+t5ZOFYTwJ;`D9|e6>II+%fU*mm#}z$MaJ8CBBYB
zfTrZhDubwSl1|-Wnsi>@Up~xa!`8nf*j((T@rlVk>`AAcdf!9TgM^4T5B!3xUoJ88
zT{<MK>bzB7%Pu}HnAM%NZ_JYH{k5I^F;8Z38G2;+YL_7{TD{8E0|kvXb(}OQIUOx0
z-_~2(HIEO)nQ6c#j{DdvIV^}bH8S-IpaMw5UF?#tuQv4fIY0GCQy7!ilKDG%U%HVu
ziCIj028_Aia?*P`yAP1o4S_cUVf#4kI#nZ_@ET_se+3T;^c3)pMyF65ftPWt@c7w3
zcUu&MT;mERd*k<$_%q=nZg}E68ZxqxA8QeTOya1yg(mX-%$_+yD(3}F>Q(Q$w$W-Y
zb4_WFzK5TtIE|PTbeFUjrRl8s>2M}V#eRk-)luS#ShhQ>urmiS+-jtwj>jwgRTKEQ
zO&~LRmNUKj6NFWjd9z&L?}UZF9+wICI;T$F@0gdET1Z_x_*XmbpwENLhe?v+?}2g7
zQ#MI59#7vtI;Pg|+IBi|*E#oOAU^kU!3f=-#`EAefY~+wv=TQYZ#oZ4?a+ONn2d>S
zFmF8MjWfkfGCNOXFI9<eEM`8O^@g+@=%Avii*Nlq^$SiTIWM6$R*~O)G-1(Apiy!$
z)Tc}l2-(e?TUlU8iTq9cThG8$<h#`aQ?BrzN_u%m%OmFH^7f{)fHIFX7Fw4_njFu^
zZ%PF<b00lV9C?<Ege5f|fv@Med<oJ66qiw-N4HliDUrhy8to=I*^EXX@vh4l(5nLi
zZ_YUt$bVG_{kizg&DnESH}w=IWMX^BFZDnUn~-8YPpC-XB*&DOPz}(3@-LHBUk*jQ
zOTY{X-c9-RZ<(2P*SO4OS4apKFU`Z9TMI4J5|kL;^A^qZJn!NQ4#Wclb?rsJ<$%{o
zn}3KN(_>@{ulVqZ0;yW*h_pO!fb_z5&hBGw-FUb<f4)3Oz$ov<6TIomIZ>WR#oCu8
ziwhfsUf{k^g<c6BbKnb17<A}Q1Mi&R1+um6Ox-B612=+`TxXpLaL0YV@^4$xM{lSX
zrNl4LGshK#Lr{c5Dl_eUeFxWqUElE7rRStTBg23B9sw8+Y+p)vib1_ZKf3n6{C(Hz
zUgR6<cf<;;DA{}dNAIOo3|bUCg~O#lbCoL#yqLIbqBAL3|9SW0xM`Tm%1I1@aHy3c
z`{frrsK35H_9e*-dAa!@Qq+E;ZA4zm`@X5?CBhg(!m{)=QNwvuE+$HB9(=&N*_%l+
zc83(HMMFc=Z?{g{Yo6;;4d+h4K<N&U;x$R{JQYm0ArS$AeUHPs+8D7jlI@NtI`<s{
zA}2~o{BzQr)wan9<Hc~DvfEG11yx+DNA<;xFkAq%et><_6ID*naPg_QT5dNAJbX)r
z9yaU(8}M5-mBl;o#gBT17brDrx?Y82wi0l2_{K}VlT*R=IU`9$<bur<HOGVKW(g=J
zKp-TfY8$=04jY6}$~3}`9#45N*^h^KMRCn%e^??4zlo;OWXd&?1}>^K?*KcKtaPWw
zA<U~!Xrk$mA?*C8hgSBrS%Zi`deUBD1maQb2_8Gg^L=<e{`?Z~A#O1q7k94R8bS62
z3Dx&zrg=q?{2mtM$cJ*21m`R367a$lKD~{uUW7&-hdY3GkX8^rg_7HMPr1{+GyC&Q
zTP<t~2@3oe_}XH%9JKV-{6@2M0GQE{LqXbGJ)g{SHvDB(?n$Q3Nu#t}X$ggq@v;J1
zfO*uAREvLkGWEX6p^rZ7-AKJ~<K<B^!Y{pUKH&0E6?tLMHqb7^8$vxvkW?d?HpVqq
z`!=<Fj~U1v_!4U^OoEN}LEr&m&}T->38+yvPgK<A5@|3&3tZZ;cP~i9#xTS$iZl1>
z)kt}1ciRK3PaJTj8P5eoSszlJ4H-1|wG@<olX%xRNqq8ks6c9dYuYgIbioa*_7B;Y
zYaV%xrpp{#-LGHaYWeqnw!60ICb~3WHE_NKtRFIRiwLw(>!?GD;uBSH&}$4m@R{HF
zcm?_4MRa`C(r}d+gtm0koY*7mkBhzIm>I~^itZ?DP{VoAiO$z8uY4oN)R!pqey&4T
zwO3zq#4@s7t(Kk9v(hq*+J9pNqG@?S1wIw`m>{?AHTc1_vqVB2%3Xp^B^;Zweze1d
zXrPhOtC5-8u52x2p$u#Gg+%b%5gS~7-hOyjRbGWVvhWvi&3kTn7eIwsc{hcBhkJQ$
zKf8ZqV2bkfq!3EhSST}(%0!|=c-B`16>SdCrHP1$C4TmA?mOQqy@%+iY}B46SFDf1
zK=Vt*d`{^`2oxp{EA(UvNzQ@!J-xx$m3m~VmTS6wRutX4rYjm780Yg%FRz6DZB-*)
zCABylJH``r+z%iBwtV?^0KG%uXhv0@!s89|j0<&C(9rECx$2Q?kK&Fz32G+}+*A&M
z&f+wIGv)iIs1NvywlSmI+oSwL87(wB|2{cc(i{5JOY$TfA~{m2$y?T=Z^XS|csw`^
zS$KXE@OT}Sw9pG)6X-53;$doh{d_S;c<;2;@wG-?ng_n0Vv9bA*o`Wd$!k-GCK0fp
zxd0W3J-Xh5kxKS7t=~hRY5rWl>{VrRSf>Ppb#;)Jj4S%}Mi$O%#-S<C|9bBoHS5_s
zbR<p&A<<sn2WH+EZsOEVDY<%S92gijG!&}!1k>m%8Fxy`PX?brQfDtNeh)lF9FBtJ
zPC^8@*X)1QV28pPyDK=_7_qNVJP{cZB_hR{;g>*Rzkar{Wcn-X?;AG)A%Pz(S8K~E
z<z@b!@N$@bXcDZ&>Gi`(cOIu4i}G=x*5ICZ3TeVw&Pt0e>xK2XWU}(3f7kKNND~<O
z5EsmKF|mo#3CXQ3Q-<M5I9!YKMZ8GenN4$mOc3H4+l+0;Bymht#}deUxqemBZt>j?
zlhwbW?z3fkqtN(rai|M^L2%rixE5$>>nnDrKS{qAS6W`tZ?o~C<8aB%++QD!O00Ov
z;wCg2pp@^)78M({cnmw=OToJ2VK*zS2@yrK>0a@2vGU>ni;KUp+q#9o*KOCjhU*xd
zYD%t#c*$*1<t!B+Gh|FUZF^{^suvkE5>#gJf*Oe<*zQj5Fczh3dx*~~d)j+$$jWAF
zDe-x=KtyD_1o3S5BKw+C3=5IkB<m->8>#}OE-{YePb6~=C*zwo0n*YkhvCFdzEF5E
z=IpWA5Ted#(`km1WGI~_v=Z|&q4!>s#ONpdgST5{)_BY>T<^**?+a&u<Y(xt+21kZ
z;7ONE&rRnNp80qiX13yYt}mn$yMW2}|F$AMxZq(t59QtsD+V-%zlbFnANQF$;a$&~
zTa47CqLG3w4HRK?zSw)jRlg9qt(x@$#q|)j>UkUNQG^DbWZM!~pv&Pc##lYrLs@;8
za~eT4H%H+yPi!o0C#vfU9}E(nj5T_rb@gQ^Yg*Q^a_u8yo6k>&GiM~fYYdR4JEqcr
zp3y?HuAsgzPuwA5n6ZdAdoZ<czu`D_uY&HPqDE?OV;Ppa{|Ez0!w`SESuvZPC&qaE
zTO@N{`CiIfR6cI!z?b>>z4TL&{cDli0oj4ZTV2X-72K*-9p`wV<IX$*L-^B*{?DAe
zSMsQT?u)EkqfJMF<!n8=Gv*bG+<GtR6&I>A&m|5o+H)F;==8J6*PaZf5yh3sAc5o>
z^(g_Rb|&Syj|U7j=rNTtn-xi|Z2u3FTPql9*kUSGHY@&{<S~^Bn-xD=9T;l<4}Z2E
zGSmQKD%CeD{+q}#mC~CINv-4zHUEdnt(gopOfi+pn-2d?;+RVLO@|+?Mip;?-BtU9
jdj}*Ii!63(U$yj?gkSv2P`_OUy<0#kUkPA)_wIiHq%T%I

diff --git a/rds/hg_parms.RDS b/rds/hg_parms.RDS
index 69d23c6a7d7e5f9af14acdb8fad2ef61177248af..5b19070ee1fd0b183c0d347df99d7ded8e67b388 100644
GIT binary patch
literal 9611
zcmV;6C3M;!iwFP!000001KqlLG*$2W|9_oxq*8{2NQTG|WgZ(`d($9eC6OqVBvDjI
zMIsfEF;S)zQlzLva%6}^rbJQZgv@h>zN_ByI=|0reb)QG-*>IAwf6G7uKV8C<9Xj_
zU;8-DHe~?-8Za>d(B=UE%roA#Wvv<jZ7Be-0ToyQ-~<480Z=voP&5EwF#xp-07?L0
z0DvX}Ks^LNT?#-W0pKhE)gA!s04VbSC@ui79{?@_z=8kNqy7Y-5df4^0N@5d(*l4^
z08}9W+Ij%mW&r9l0IC20kN{L+05}dnD+Qos0Kjj5TLCEZ0jSRbC=mdZ5CE_MpmqRI
zUja}~1Ar|6?HvHM2msgtXyO2rn*gv6fVKtzmH@yd0B{AMsRF=J0Gbj2cmu#90Lm2r
znhXFX9DsHg0KNiH%>XFs0F+<=7zCjD15gbCDDMF%-2l{f09rW!WdQ(e2B545fG+@4
z7XTUpP&NY4Rs&GV0Kgr9nhrqQ0sz|pXx;$SVE~#K0B8ch4glIZ0N4mXtp}i30npY0
zP%i>de(xy(KvM&tE(CzX0F-(F-~xct0JJ~=sGfTbnH|;x(+dHh-ctZ-iUUyhYyj#7
zMWCD=2FfiXpyV$Eioj8zus8ul<0Mcx@_@3Y31Ex~;A;-R=u4n1$OLFO1(XGy0F{ma
zc{%{eSpb!p0FT1~-jo5{w+6^)0f-s`h|dHFUkC7J3?N(`Ab}mgF9_hmb$~nJ03IO#
z{xtxn6#!10130!2z{3dOq!)mFI)GI^fX6|A!<GR11OZN10~oOa7)}G2Rs-x91+e@8
zu<HbX@mm0`699%E0XFvm81@6`l>)5O1i%IWg;oFs3V_Nx0698<j5vUdA;3C60NK+3
za`ymatN~=d1FR(gB)0>|O9RN-0!Z!zkVpZLegYu+5<svIV5JX$@N0loo&dt{0EDvu
z#FzmDwg8CS2M{p_5M2VWvKoLd8bBZlfY$+lR}p~c1_0Md0KW48d^G@E#{syO0`SBD
z@aO|@p9SE41;8B(z~#-*hj9RK1ppENzzP7)0Z?QCz!U&(13(M_gaJS_0CWOCJOKRe
z3ts>@4*=Ex5Id{=XLa;!?K892z13_S001`t;P-X(`&$3q+o1s9Gg~iYHtz;Nq0IUU
zo%NCRpZxDW1CRgoDtG~a0|0ymfN%it`%h2*`TglXYk?pDxHjw4XO{V6c6}BA;5=K$
zbJqXS|E!n7GRr$Od!IQtyN=&Ja{)jb01Re5-kjZ|Fb2Nk?Ao1Yz5l)^te@SJJpk}w
zc0W{S_v|+(5&(>6>sij$IyviY6#z8M9=iZ2n}G^!zxNfx4WBrDypfC{@&!kNUOvb0
zZXZf_gEYp}H-A~*dI2M*)9>-bPvGTU?N963B=As2<7Dz9J@gb$OJ<hu!C<+ej^xuC
z7(_39%jsB)=amL!H|B_=$H%_Z9oBUiXyu*u?P)abjIelP^fD2{ol8G;IL6`4G3IPb
z3p0$l7-T!pp^eua-pzz=$;X=`$Isi{--)rT$`kt+3}Qg#jkL@kZFpfk?Sy2iA>Iyu
zguia@z_XfHH-F;wK(~#$w!>-881&xi!EwdC82Tl7<8|IB^ba3dsI_7rdhfaR)vqWX
zgR(`sT$nV_|5fcFl_vtY{h)r;K++UmeQ&MR34>^Fq%WhiA`L@Lj=!$bRl*y!8%0fW
z)}pOA+cb0TcJyo!H*ySEfr0Ms^bzwtc=^@)c?-N3<K?vx;vR2bqHq2uQMvu^(caMG
z^}+Ss7{0JJ?Wc|khHSXIc-yL{7})DyU!!~wJtDf_#F#w9z;Pk-!EO8)+#wS9y<rM{
z60~{j?hB#!Id?&AKWX%v?Adw!)>n*-2>e;r9EiSBjjX2nH!+ysIgCO{!a%_x(7)e?
z0UPdL{?T(79ap%ro{~(STbJ*e&f7w2MQBE~dUtDl1_KuywBj~U!?10S=C|HGjF-xL
zJk-v0;ML8$`+SlDFet$BbMY1jyu@j$SI}*Mo>}3#XX<@0Of*|!$sTVEQMf6-tZO@7
zslPe&{^w~73pnE)79@{Bt6hE-j;3R1&`pOP`B!*z<X6Etg^w6$Xt8(8#~&D6ox0oc
z+bO&(v`Or`k=LBR8_$hD_}tY%A1;rT<nS8wrF-`;@B4{CY;i*GrhZ_Ec!HLj@sqjv
z0krZnMD!q9A0S5-E^Eh2k@H7K6ecjRT$+4svKFuAX$)Qa{sFIhT{K(ua}Yh*1v4*w
z--h<fdlkgmwxheL4^NFBiJoU{T-Ef4(evB8j?+UG=(%cCK7Sbt`el4Oxr=ug!}#|^
z_0BiK!0B<HXA|e}RQHc9oV^e5hUJQf#HAYaisiDO5U)V5=kLGdOPt3GW$b?W;aBmD
z(?H6VgBLI)R>Dr|<!8L!M9$+Q@8MyT`dBz^fX?sPrRlup=s6KFSgP$X*MFmrPrf#4
zS&B2$y#)o@y~MmB;bW>5%mlMrsdfX;Q-W#eN#?1`2{;tI_-&R~5BA&2>37>4#}QIK
zTQ`FR#|)hv6)HmU*Yv8%C(5EYqFk6fE%yiq4y>|>=W51L@f!m*c3#+)JCHM~HA%o@
z2|+2PR|NIA+~<Rbn~6n5iwqSi!f|4d#U-r!I|0|o9ZPhFak_7X^k<PqoHGBBc+Frn
z4i~<*!ygs|E$iwzuAR%U<w@_j-PTh4Y5CY7&fpl189zH8P!UQ{xJ#EEGdYYSX7S>J
zotpSl=tEm!tS<gKd&?ho&*1ckQ3SO-2mAQ1<=cFz#BtB5+KhK;*mN{zp^Ii7K|QqC
z)!m|>pthf38h`GKgM9b%esLbapI(ET4rh`$v+UrcyVFIS*dN=eaB??}4`ds1JaofJ
z5rYXEV^{3%KYqT;YCplGdvZvOLjb32R{dPq6oTWg*-X48C9uh&PPCJ<1Ald@bH{V+
zz|mh7T-lGSank;Jm;dGSIO-pLn>>yL#bx?nje8G4^$Xf#rmzi%CPq8C9$d#?oMKyR
z^1tFhaAJ8!ekyikw=DHjNyVXiCE@}f*JF3r)Ot7Dt2mONm`zSz#ojWzp)FE-@K@Yd
zRo8S9$J_YzLO0Q|r@3R`>{KvLHtclWF;5<6jxR6dWy!+n<%8i&S4(lw<EJI_Ex_@@
zV*>LA9}<*G1O5-@BKWhgw$;gkz?sEWhaX*a#_^TR&(B#E5!8mP__IHoapL({A-443
z<mAADOI<nyrIq(o=z=_)U?<>Ce+>5UxG-nO>*CM0zBff-YjI5Zm@JLzf#Y|_2RZcB
zan!haZ*@%_PWv*YJpQPUy`fi4>RRsL_|^Oi`mV=uh(F0a^_D)4=p6z7nA-%!=Wc{u
z!zS!Y1=)4cBRK9G(83=OjRW_FLgU9caj2wI^_zeQwj3zxS?eE(6Ky{o19I$fN<~fZ
ztl%~5EW6Ab%CQ@VyV(U79{NBq@m{^Zq;L_A#k0^l+g-4=Hz!bLu`hP*ykAy)LIo#B
zhZQd@jK{H7*E$DSg~JJ**EUQNI5bolpZnnicGPW=Dqf>bP`p`ns+<$BXTt-1Tb>O3
z-e+pc-FXEE_ANdcz`1R%|Kbhjo&Togg|9#DvcH$DjgNimsfp66ShfE9)VnV%Sjw$(
z_(88TR@?Y3Jzo6#Z{YvW558qTJnqs~g3ru(1baBI&E=cU@QNJTRf1KLJbTy=yu+5D
z{Fm}wPWa@#Z*b4Y0!;bP`^oUtb}TuFP7m`|;`1niJB2S2i}!3S$+0%UycIJ+oZ5FV
zv5!T0l}5||<j!#lhMWhu0)#PhFzQU`hhn@R@M)8svj^VWH0I9!$$4&G65D>yx}|$D
z&s#h`D1RIJU-E_2=Wj9RO>d~Rydb_1bt*G!U4bQ1BNrj)*4%lW7lzxumwdsR!c=o7
zz1`TLn7d`k$Nl)UX>#-#o9n;)G}r8Hh|@WQZ=wg|%5TNs(<EBAWXi>V&8w*UlHffw
zi6t2)j@gRDVDmj=)z?{E|C*PhZtonhE^{va_StWG4c9A-&%YcLQNoX9J#7zCuP4uq
z2T>ln>)O|0aIVE}vhXZ^RKMwP@Ja)=jho3!9L~j;qQA<15w|erPM2%@<#%{@%?^GO
z!7un=af{H$A5ZXQ<;=1+q0jh+nR?>V)4llKki+Ox!WF!E4&+55R$^=T<($p|8s?5K
z%XvQCg@xC$2cu+Ku;8s$TN76j7S*3G*j2ItpFUH)KGv;_ZzTNP-aCb1``+6k2BzWY
zS~0oMwDbVhMY;_?9M8hS^-9N*Q}*GT7q0Up1sBh)BX3MKoGDQOU+^}lp7t)p>ag)6
zgTKNtO@Q~z$;;iCQ>H6n+1r5)n&+K&Z9Rw8sXQ8)4;rwjhJ!n0VjQd2sTuRRCgGb}
zn)`*8?Q?b8jUNa$pR&Y|<}%8T8?jhy9h={j6Npt^AFCU-R$|hYklSB952B}#PVGoU
z8Rj0ZV!hZ*pF6+i=l;xLv%{Es4zj$Pg0Wya`-@CsFgBG&B&2K-oSWB`Mq3jua|pW(
z)jRWeJm&J{rxl!SGSB?&zqvfoCoYnGffKi?d`UO$mqhY9|D`mGwYV`wX5YgjCUf(Q
z@0`v`xOEGS_K1eY9(y*I-zn*o|Kj+3+%2Fl9KYuC|JpASvFi7VSrpOu^6?dmK0cny
zYrbz%W{LfbIuEoO?V1|q^6I*a7jVD*iaR25^EB#f(Ksr*j@bJ4pYb$ycC}lwxJgXV
z%Fd|_O~o4*a&F%|H*R`F`Nqrli%`X;z(JBD0}YD1szm2WqJiB`n^5EJs4f>hH6gKP
zZa#S`=7fFi8`P6C9KNI@^S`a*U+Y?GZx-U&$U7Hn#hoyvPI97l;5XLpg{-(+&Mck&
zHu}HT*|ulVHLd!~xGpB8*8gY+s@&g3MBj4yujgsCN}2enC!<QoF8`Dd|6Ff+#b%a8
z6F=wTHG8&g4k@C{#oK$|JY_4b{8zoTOcO;97p3AVqsaVQCf!JiROHT}A8z(heky+=
zaBkja*S*uEOFkM{4^VEe4MC$+^NPdCVYnmW{<-a61##8us3c;$^}p&HEYKFs==sN|
z>ir-|-Nvt|O8=_N(a|xN-~G%ZKdJHke_Q|mx-W9yh$>aN@_*&^YPGEf3QJMVq+IUm
zTei8prZ8pgXwpB|`}S7`?}u0p|Lgqq#Rn)3=c4|V*O4~qG<ahEUyc2J|Kk!?ni7p;
z%E?t$R4t>)<UIY{+a75oi6YCxqFC${P;6@!DHR)wqCE#KwmtDeQPz0<L+6>0vg{et
zc>WMdXpe;R3%Y0wU%a!3cipX=7cVyCx{k;q_Z=Cvk22pi`n=Bte%8>#)dwn{W$jVU
zDGKz+8ryvhSI0)POuCif$_EK*n&-6X?9V>kiXZpa_{EyOG3cfnZlNADUYI|R&Xem(
z7c1Aw8B}{T*yNw6(WRMosmt{+T{hUuK*_Nbl?@ts3!_U=^jk`Oipmf&`&0eY9cd_R
zyiYe$l2c<~_t2+^)JS9&zFcGLBT8pBGxFtJu{UQhxgvOV2|r42=_|AmU_qry5U&(x
z!&Q8elFXP&7cK~fHVFabUHu{SVv`-pUF4Cw7;K0F98-3^^9$%I%f#$HroF=jR-bfR
zOMP>?t{-FG!Tb){H-%kw_mQD<DHHoiZZ<kaw^K#$?V+4sa&a+dSS#o}+@5k<<^`bQ
z-2+?4^hPvVCd)58=cz|AtN4Zi_JW)-=9%Qdx_&xrnDY3{CZN&0Qgi8@yk)r3v-WD!
zIa?I*V@@xR5kbKPue%TIxroc8t;ej>Qga@qvyV2zA-d4<#rHC#?xVtPZ4q0rqOZKW
zna%6+6_nuY68#n|n=|7pyMEX9Mr2`AXe_<wghIQrb%Ycpaq&(Mf8+HJ>B?f&iBZn-
zxQx}4iK1};S1p#bE-^W$(c-c|R-dyXXS5@>f*Kl0-w=G(u7dlH#;5~@-Lfqn<#upB
zqNUc*DHmF7T0KA0SGH%_F^zZQqTIU=Id1(x!A9p3A6C}mlI9mSf$Is4>D??Eh7$|u
zdUti|(DF4MGBp;zw$jH%rdk_HeuZi@RBIPK`KC%=iLBqx)U)FDMmvWetH$YkQH>`;
z3z{`LER2_}?5M%*uTu0+duE}K$h-#$imP$aw?6y8y)&r!^5kUReK+LRcfXyLn4<A`
zy}Wsu);i>SHu<!E^ATLic4YC&y;o2xrCqE;?ggFuHu0l&K@^=iU57o}kRLg68ViRO
zM4*CDr@m@}Ao8gS_cBM{K~WFa<o-}c+!#N(qSkO73Y3@#QQ!NcXkTIOj_kag(O+Y`
z#ZIk4@dPcgrcoB;pI_|$D{v3}@8?*)#{xOT&+fVYt6*n8>LB+Ri*7D)2;35k&w_@Y
z&zG&HzinxYyHd1?em6N`#V$2N{6H_R(vaMP1%a1d@$P+(RXZaBawQD$wXNA<si9w3
z+|8p>5)_Ih4$o_hcg({wNnM|=w>?;8r(}AbB^*CzKR)ntKZ&o+=Z`gVr_tZob^a*j
zilI02bo*)_YsU8(70i)0bue-3-Zg@zl2~}=yt9;FJC>xy7QVCGi=|QrcuyD(VzKf&
zpC}N+EUobt-tSsi#QsS9;7}_*FB;xjf7TT528s^#FS&`YZ|sZjJn4c33lo+u5ud{R
zQI%4!kP&?1E5%nBPtqT0#h7SbS&XG)Q(D6zT-fAh-Mjg6AZCONeN2Q}%ohlfkaTmw
zPrLLvCiU`B=Z4KlmVymdIt(Sf{#t{DKa4JxU;K%0c2=1!oqC0NZSrH<-ydTG8lQ-s
ze;Lb-zT0|Or(nK|KWAiI3KntcCVo6xMQ>4jqNinIf;j<R8;-_?V7}hWLe*W)_#tTF
z@;E6zthG7%?nKHM7PfNu->Y(@hsPbMFgH`jyy4te8%27sI{8j6^Kc!0{IQ`&%sLaR
z1&Cq8N5)uruTD;`&j^c{KQCPy(||>s&u-r#_2{qqC-&+$1!A>^lG*Z~%dkc;qy3YL
zJHF^Qmo+S6!OVSOl=Xgo__*%D3ho<|n4fvHS-*M*mSjdc&44w&4{WNO$y|d4t_|^W
ze%|=f;xgOusfAc_AuGeFb`w^sj$e!0*^Egs*LO!`9Ke!g+YdT>mSJ&%#Xf&$QOtcD
z9ckLag2iH$uipEurpGY{Jd(NJiZ2G&7gRN0!5sPRO|G#b7$q!!a-WSiW^I?+NPm%r
zRccok)dnQukb=$LO(!eqi9qKt3`xY+$cIs}Dg&6dOZK3e<#tR@%HoxIv>cxuO9<IF
zavd{{m3`x|tfODS9J!qN53pQQE_Y135MMrA^x8SL8s8FIB($^cV!@`zPlVUMrRQ7l
z`Z$RtV5QA-6-&)L^t-m>f-~yi+fTZ2U$iahNe;^`Z&z-@x?hX$XkKW?a&~@e;e>vA
zTyg^crV9tLK`=71h&YN5ot@4I6fL4>EZLHO#q=yz`|7V+on?;Is+;P!D%4=^e1Ws(
zZ6EP<<c@AL*WkJS^BL=MenB(9o35;Mz2N5<N@M1SR8e|d(hy}A?JE7oX&3J^XX^0o
zWPI<VwaoP32l}=?GA}TWSTy6ZIG%p_<vZH~;gxvH`ry1pWB1YT;j<fT7Eu^|WdDge
z>3a10ymq2uOD+AXcJi^#0|^**Y`4&%9(j69x230a*<N~DT=3FirMDQ_&aCJDsSGbV
zs@CfAe5K!%u8j*B_r_QW$<aL5eHi{_>b!?v9Y(2|ZH+1(K}~HElJkGirNoRl?+*8%
z|9}}kmqa0MDS!6#7ON2b@-2s`FL8Bv-NsaHlZgi2=3)BfVB$cRINKQ<mDWS|VHvnL
z6W)eFV=CTOyj*y{y+G9~81N{x9gE=FinrBfgc57>(BoHPxWS>5^l-n{yAQKgVhmhr
z6EbMR(D}9@qK(xU#C9`j#pWXP^Pp}$)trxaCbQl4wnX5qCF^`UbdKOnTbGWPLhcw^
z)D@<Yeu#dZf5Gz2dfoJczxZvQPRXOaE?J}7AA+I51^)RK8}Xq;&v1-E1%@9Ji_Z)7
zrC(6--1zWx0$sB`d!e|L0S3O)l3B-Ug!dLa%m~+!!w~L_W6~4rFfd%oceqlF9=ZBk
zdD4v$B*bgZZnXP^r+9nMUW*FCn<_IVA2Ycy^?^>`Yt{J}S-wxb%vcRw>qGY|n!cqU
z!JjV=HEW`GuizMya2a}e_C8h_Z^s+494~e~-ikqEvZ?89ns{T(eyCKk8gGjXM~;@%
zU`nx2tOc7n-94nLOrs?XABQm&2iR`G;B=kp3o|X~*3jMo7uI3e_*L`%lWXZ$-KwSp
z1aH%Kui40Ilo&|2TIk%p_Dd<=T&1?JbvTdib4YsL8lPntYrJhVnRNi^zN|fF&1v-H
z>OEqIWXCb($dhL`_ej&DGOt&dmU7Y~<e0yPhpnS~Ios4fRJ%t%&onyzt>!&F@%?F+
zDyJ~KE5+yc@`oMX(OmjzX<7(|ui`5Y6ZXO@rxz^KTcb*U5*8eob0!W0)-)x?w5Z^b
zjXf{Ko|oWRk#jOf{X^;ZO^UM1_|{?c;??_txNI>_IeQt$)>rg6+uku7o%{HJvhLX}
z=R|s7EboJE-EWv^SKD}`?K1jOBexd5$;H!U-~6nDT<CY*qy^3y?m_!LfmMfBYs~du
z=lDkDLuGaNA$4*Z1a{+xw1@m($E&gOhlefY#V8hC%RD<$SAhL~x7gkC=01mDHI;At
zB8MmT9F`hLE9}Hd@i(7+hcf2Izhv1dnb;QK_uFGNr*m|%$zuNd>#VKVW_q?vyYvS}
zR9EVcRSM6|Zz)nyb=7Le@@+PE241LPSN_-<Wr3IYiAgEnrs5jro#YDH{r)wkj|SaS
z<T-&~_!qU_+oL{reieU_b$Ox<X5832_5G9$c9uR4cN6xW8~>_lM=qP*iMf|kr5$?@
z;oCN%-?GdXKNZ~=jnH?*)@>Sf!{b_5c{6+U?WuOGHZnN4`^+d-(cek6P|W7)^e|Wd
zO4_T96{|??U!UBu<<?2+P?QsX+n2{RIe!f6HT!);%?z=WZQdFo-)5|8OyqsO%K?kh
z_g-mxY>$1-x#~-2im^dx*s<@;9&A@!G~X+J1e^AI%>-0uV^#c=yD-Z&tdd!G!NvLs
zHZSpPs+X3+_Wk@eefvDH``B>%+j19dS>I<hq4^R!JeFeaJY#%U$kfHADu%@lZwKto
z&%-urvx{~o-(ppzk@KspquBQI+_Rh|V%QwIjjif9C(hW@=H-4E#>xXnA9z&>VRKEw
z`Gm(!m@rt)VVk-h-z;Z+zGj~?HZ7Q6Hdb;Nn@Sa0RAig6L5tP$Zu2>8Q)s2OCYxgQ
z1~JZ@`LbBXWuCL?WEEybL=GC2Y`{nNZHC&vM&p;Xbm%HsgdfKuWnZ7zj*V-4R=P-M
zWARj?o64$dSji!-l@(TvP3snQT&ozs&VtO$n$9F_j#YWO;GikegVd*bJt)|)YDllP
z<Q~55TAE>Zgv6R_MeRx6wsWsb=aufYe9Z5#M&;-t7kM-M80m7iRH*_39)-@w;Q@ly
zvdr#q`(~_YF&16+m=lxAw%&;mGRAzb@|yH-&+#j-{glpvcUbP7@211H8h=b*OYeCT
zGgp7;i08#XS@*fP)Tva0YO8~vnQut+20g{(eIgu??uxZpA0HUVL}8Xwvr3&oF@A{?
z(BGIEj&a+hb-s5AV}IlGTUsy5|Mnlb>9q5s{fnH)4d#O@R%|aNWj9@4cE@}@DcvHa
zcZdBnDc$#7DP60Flv|&?jNYV63dvb6r8Zw7MX$v?e>}RK6hG(@sc=S+T&{oes#aVy
z$(LQ!cE$2NxnfA$UXZ7l<j>hw6vlsy<ZW*jiyP-7xxIAKBE>^UF7<C`jKb+8mr1Ni
z+Oz@5yE;qks*^Igq|r(?HN>0b`!*G@g*lbv^3LKvJ>E!iI_)^Nd@71$_YFOiGk%m@
zQng`iTF_0BeTRcw`1B`|LsL~jaNlu~ZKZ|x+9yv*&VwpKuh-os*=!q*UHou@WRsBB
zzIa%NWEI+Ir5XK@q>iM1U>S8HSw9ucYouk7Y+793`Bxi~EERc=f-lC93s$HucWV(N
z7c2US(^7tr%ro5OTLWv!#o8PuQv2qS%%X9362r|%rZ?1WN@O6(^hEMN;_whTe?IN=
zI(})AnR^?nq?$iTJt^pQ)#oKi+kUk9k@_9tm%F9SYVPGE^?Kj^L94eUMP7H$f_{FI
z>7{+p%4fbLDD+OWS~L<fzG|Bt!&VYA<F{m@4~vt)`MmSeYY~#RV)#qc&Nkv#$<rMc
z_7)_#b*i;Qn-k-Al$KuTC?_Ug7*=ld<RvDv^!$Eo_(JrSY+2uOE}xi^Nli+bmM4Z6
z73=Y*hZ3X9j;!!8eLxH>>uZ34K4P%aK*Y+!hZvrJM=Wllo)~gebA5661u;N3Xpz?1
zKnz|N8#xf@K#X2;^lezDM+`pGCG`i&i2kh#JD6)lhykzRq^bgDqR0J`Zd-{Q(eJ4I
z=Dno@(PtFS&7WUR^rz=+s1`6Iz6r~Z4LEZVot46&H){EbzLz%J7PAc#{hRGB-Zgh7
zeoSRt%5%yfIt@ybZkx3d-xro^X<j`?bg}jyc>U!C@!j`CU=$^q=&JTRp<$U$d=Kwb
zPn8)ay1za=yo0Hi=&<@MpSfl+(NeR2<<~dwh;Ir}9uhz&I*O~kGOrdB9R>r*Ch;vq
z>zUGL1}4Kq<8#?HD~q=f?L4jWTY7gBO`Amb*F6y-S`s%Nc{#yEd@%{j6m-rd+FPf!
zB|`a#FRGd!#QUU(da($@1K)ayFE5hh<<v@vy4N>vX5Tj@Djy%qxKhVOe6gJ0+${W*
z_{v`WX?6K1(d72lovivyG(TY;j`sQe{X5I-|34}J7Z3P<ltuaf<@kFKJ;og;peO#7
zG#9_{u`te#)XTfGr{nV%azoVivGM=|QfDXEl2!q6($s`^&6WobNx8DHXCfJD<fi&q
zwI>hlN&SctzS@oKr0ScP^)IUrksHr(Xnua=NUFJyewRtwOl}&uJi3<Cjg)rUYyTtn
z2f6n3&*$gtSCaBf>#k}A{2;}->O)rF)gYDXILBI?f=Q7j&JFI4BP7NMpS|cmM5-Gn
z?d>!#B^7Od?KveLN=j>_W{lieONw(imNsYKCnbFD@1o~#C0DY7g5*zLQsVahqg<1j
zq`Kzw7?LhV%8T-emRA*$tHarJtI}sk*$fMn&4OQv-i&=Ym;2|F3g@>w3Uk(yt6KyI
zUDWx>HN`po3;T*l$u!3#cS;i}nmo;9HMyN!Id(jA&rkrVxau~Oo%mW(BFbRJ5zZk}
zv^+)aey24lQFKH8Oie8*qT_M8@Ju1Ontgd*0E<4!UA#rH%-Mnz+xVt(&viCptoSr2
z#N8u>)n0YK_-INByvk}4RTm&d6Qr&>Z+J<HUf8fJzkCZR?3f|Q8MT7szTh=kPYoqS
zI}GgoIs8e<i^4{#<6)%GcjqUCtJabtN3B}Uvxtz~o63tH+%qRR8||9#)MZjc=E{z3
z;pd4@uIw(0o>`E*9XI+dI$1~-`4~HlDj-=yV&^jr2$4)nKHIAAvmiOE3Jz@)bRm~m
z2)`*mBS|hdydBkA_kdhlIcR6gq)hU?a;Qv}G$i??GOa{jAi3<_dckpB0g_r~o?lh5
zgk18jEBDzg7h-&cqWLLQCl{qkluD$mA~{dmB&}xcC7B}fU9xS0h=z$uuicH0i1DXx
zO)K3BN$!H0Xu}?VlC90HaXO`z<UUXT+`0S{xitQ8a8_p#NgY_+6_O`JE-|}O!Y86l
zOnnVdP2W*Ua#i$LJfYnn7ucSPp4ebXE^7I>!A8=aWbbcXo+`^iv<9+%yPrKlj4pH&
zebApvO!G5)+Mlf>W};-B2n}nJ>7=A4kM%X8(X@2a-i16QkFwyB1a}c)XhnCZ_!C!>
z?edRD?(a&8G39wv3$`95d5$m83p_7PG#}J^uNy=qdJR5&Pyb9KCU>^n@Y~T&eBL&p
zJ>lq2(3kG6F|=Pt492Ybc`1TI3^VWSI#^OcQX&irSEwkHoR?UVdQWWzV50u%LniqB
zyy3j#2ka03et?tX&%>fjf09h691dB{YRd6hLDiI-Q}T04VNNN|DWy53Jf~FVl<J&P
zn^iP&o7SvkK6*;)#Bt}p13Qh3{)gX%Cr`M@FnT`mFH<y4c?PNQk5v3e{<TQcl>h6I
zrYZl|Bu!KPuS=Sy{G3fqdG)^=q$w{i$DrjIv;u=xWY9_sTA4wsFlbc<t;V3$|IrEz
zTm=TM0s~ipfvdp4Rbb#MFmM$ZxC#th1qQAn16PrOtH{7rWZ)_?a1|N2iVR#u2CgCl
zSCN6M#K2Wz;3_e2l^D263|u7!t`Y-RiGi!cz*S=4Dl>4E8Mw*}TxABXG6Pqcfve2G
zRc7ESGjNp|xGD@>6$Y*f16PHCtHQulVc@DTa8($%Dhyl|2CgatSCxUQ%D`1+;Holk
zRT;Rd3|v(Pt||jpm4U0qz*S@5sxffY7`SQ-Tr~!+8Ut62fvd*ARb$|)GjP=zxatgC
zbq20F16Q4atIohxXW*(caMk~E=l#BK{009B{*wO$f6;$}zwAH3U-+NkFa1yO7yl>t
z%l{MnHSi~3bRaMHKi>on*&R6b_lEj6$zpZj^Z^;`lLzd7e~&uze*mBjANXf5004%p
B@pk|K

literal 9809
zcmV-XCa&2ZiwFP!000002IYErG*xfZ|1oEthbU7tO6IBDeGeKaN>XIJxMr6O8JbWC
z4MIvJl$2CTWQe08iDZb<BtwRz%%tKym8;|XzW@LJ_*rYOz0ZF3em;A8o^#i7EI1GZ
zVM3Uh84)HHTtryRcWqj$h3f?I2VR5`;l{Z=PDc@hVFQ9NI^y&hmr1zp8BScdtt`&H
zalIGLhw*0=r}el@!zlzOPkgNqCwx!FY1}suU)RU=VK^t_+z#g(@wEb6Ho({MNSNwy
zK7i|_aaxM&58=;jI6sCzOW_>fi%9|JI2MB+&MElweVhm5<c`Z~xLq$U+u*z!r!rhu
zi_>yk_W`HPxb7tGn~YNePEokd?N81KVvE};;+!AX`{VLuoHTL!?fCOKoV0OTgs;EH
zb@@1b!sQ6u25@;VzE+RRn{eF-PMo;>0w;W(aSUHOkJEFU;&IZ&b$Bfp-{N+9I331y
zcr6fIj~L=Ut8kLXWz+dSu(a((n3NvjCA^0)IEEt(cW)w$t7r&=N*W#(p6B=!f_$|?
zkWVED@)*xKO$I@tWe^6A83dWQif>_xAUB;5BoXgh5Z<c@eE&*S1i8EzK|=9b?Zf>J
z;X{#t*Co#qK`imwOC;efOhpi9WduoPLXg1Y2x9aQFD>5lXuJ=b@VeWiBZyZGf<)r@
z=0OOu{R^H#DS~Xkd$S40^T7MIwI4y0TJSYtyv;b4BEH|2b_7{<5J7h1x$kO1kY#v}
zRd^7@0iP2sQ3R32arfZ!xFrigq;VWYYXn)o2SK!T5JWN?mm6^#2LxGy+lZwgi1I63
zkI#i%B7#Waxrk5Vm^d!rw&HjU!gy`ORS-l5@8J?W4r#n!!V>r%_wbRz=R+9prFa#B
zh~T-2nd5tO;4(h{V)(2s-HP{t?@u3?5q6w&dwM%~dV6{XOZ&LdsB=noU#gR|7sbbu
zu2|$3NTvF_1bETC_61Np+-NS+6h9v-ZBEJQ7vMtiUEs%Q;2q%NN%33j?B+YC;&7mN
z`cb^3|7q9*UETbt(%uvgx``j|;z|sNJDBR}=^eP%$J@t~LiZD*dHediQhjOCe>Pa_
zNOhWD<C|amoyxzh*#jy5R9}LE$IlhF^rTXpiFiCd-ZUqwueTp%f!5R8$%X3U?MWqS
zx%>htUOv+P0X|euZ@MMwMfG)aq|l@ts8mN+x}}3B#nD6B5s%9^m|h_tiVq%;gPZTX
z*&-irKYv1q0#u5h|9qMYpz*Ez{AmG1>p2iVAKw7S`F{S@{WDthNAa(pd~^uBS_{fZ
z*Gs!ny*#}g=RxN4qk7W!;CIG5glqAsaPoEu^!9e5qr3Q0sWg04To*R{pJLJPCK2-D
z4W#~sU5LuP+s)U-jpnwXDRa5~+l<>U03Wg67K9n&^Yae)zvh9<+@Inq?H)jL!@EG&
z{T;RiJ>mH6wxFi`|5{A(^CP5-pD3Eke+Np!(c8<*n<nk;EA1cXy_QC$I{lrZU~l^I
z;qasc_)<NIv+3aG;_vP4@4Dca@ol8J;8my4XbS^cgs13Bh(nOxJ6CTXe8y=GR0_RS
z9D#V3@Pz&~kJHbQ;_FEz4ie|@GvFBP<4g7Pqvz~T@uK+tF84|r#oui|brpT#&;8&)
z3C2&q3(ecpjY41W!rO-jfsd%Go1=&K{M_^SxVijq?b7+8(j!)o_6qQG^rRB!$^CZ(
z#1Z&APyf2j^;fZ=JRF-TUJeV&i_ay6=0tGg{FB9kzO(z`&80Z|Qs^aMa|l>)eGmf1
zFVXpFbEUXAP<%c9o3w*>5Y-dESQcE8d|uuGG=GZQf5)8@ALIWtfa~w|_2;hfq0igg
zg}UHw<o9#?J%sczBurVL+n?e7Ul&Fo#S`xbJ#AN9<v_Tm7J1^`COGlYPw0Ox41paK
zM>l6TM`^#`4!*>m;2IzNkOdQ_ZVnpHIlzw&&Ex7#qvCT!_1hPKH*r46Kb8y2$?@lf
zLkyGO!Oh?Q{~9+gKR3_)|EW8xKZW)$qkrewFW8ad<?Ty16WNKkBv{&*-lT<>1&5ED
zzoYB_dU;ZQcX|G@^mC>D$Bmm_KfGo0!}jl`g>OeFHZL9(-a+aA^nlxq=18Ua(O-Ow
z2+!~LI@yVg{Qe``fBefG=eTXfU+y%=eTV*XOJvT@{x7G@ai{dZd^df4&W|#;A4U2v
z_nPCB^o6*T(7Ao){O#trefmN@eP7EtdxyDohq-l&IX#7*7kQ!HfzIdjwsU?yd<*lo
zs#s`8kIw*ExK5AL!G0mG5j~E%eC_7+E_2*^&Tr40o&B7@BRy~OLO)x2+|mp4px2vT
z5BHCA@y$K&ne+Fc$5F8mkKR9viiLUZo%3_0$1}Gty&v>B*fD44K=+?px1YnancL52
z?%8YZnKI{RKevx@=t6sX9Q1wX`CJ_IdDuE<XH8zXPS?}xO|OUb+&=DqpXuiS&TV43
zyU!LNs`K1_g-`Nusx^x%hs_wy^hOMvZoUG?GzU*z*}MUczMS3ipx_*&jS9&rByNGw
zKy&frTUL<P;2?eGA{U%E(XUo4H4WQjj~*45+yaM5GwGM3_Cw6F%0cDlg%HW9YkH8O
z0>aDRyWZ6^2Ny@>@izt!;b7UDZ7<_Q;dtwYll-nJaCB9ip^@k$oMefTGa92oXz9N1
z_9q)5;=^og`0*Mz%vZ`^Zp8+v6=w2#x`p6?o@b=~OeaK`JmsUtdqSM7(NX1XGH_<i
z%EN08XG5&xy~Q^Y4nw?RhD4;^NjRu?$fZO38yqFa7^Wqkfn!t0K3%)I1!BT^^FA^f
zLri<3QCv_0g!B4T2_OCpVY$}_nzo;Xh`QIjhM^hY#_M_@)t>=kVk|6GBbXo_>Viza
zZidkDs`YGb>2P5v@QXym4mjX@IHpS76v9{RF|WBh0?~RYx&0FBAu`|Ja)n(GP}CPE
zrPXLcg4go^0rCKxR+2MMR%wKz>+@TjF1~_L)8$7$%h$shokKaR4MZVPS+B?a;5dZm
zB8y)~T!5pmL>bKAAA#6z(!2cHFd#RF#5|U>gs5|vT9s-uBwl#Qm@@Mi0*#LJXNN3>
zxRg(1=}$KxE;-gEY?BKdReURxVzLx;%wvN~r!pXJ_eE2oC5CXM@ZgbqJeMJ8?d9SN
zC%=IIwzu8aueCtrYT3T;*{yJRhHnMa6&uiKNKbvF-3-2E_d78$8HlpHBbT_r2oAI^
z^7oj10C7<!(XXxf;B-OjI{&B~2;Uf-_`~TT9Dh|Osa|9b(NfmIpF7EL>}}T}m17Sf
z#<xr5<5qk=g!=mo&+UQ}`;SGCe+k3UdhQ&p{XZbw&For0un<HRE0i9vPK20}Qsz9o
z-qfeb)enmf!4~8mUoIaXoHwzIm$4LrFp=D8*gXXiX_O^SJVmfg%TV>X$Sd%b7<lS2
zRs>-_CmVi=bI+fH31_!*r&1mmRZ`z<%XtgNB$n8un?Hl85tf&^75p%4T=_Wg_!NwE
zjp&%A1d<v1%!`_$07jgHdJ^A!gI}!Ky<g58hvB-U{H3l%(5ukTSu}eL##lnr;)>m1
z^nRPM#TH@sn)%K4hLtglufM`@gnAi1m-zCvoZ*3ym!Ug&<TPO-{E|jkNC6BBeHKWq
ziGVTAKI_GT`(XOhYK(tL35<H)4dbzR10!nOM|OQkBC|Lhi61Fk3}aVT2qC6q7-z;*
z-3n%4+~E^tOkWNr75e*us}d$x&PZ(a^@Ax!v3z@n2QXv=gADE-Fuo#?^_%ty_>|S(
zZD;umzIU@t$v->;Lwx0ICb_R+i1T>VPK(np#VKs?=*$d^oXpHb*0aH|(k9!WnNj#E
z_hwx&t0zpS)(UpWBJjobwAg_4IT$~^*mJ`TEg0F<r+YE`HSTva>BIFI7*_q<Od9$I
zQ`T$xrdH(xHk;t&lgA5F?x7~kYfr+jrgwSwLp))0T%>bFgAGg{FWzwW$OD+Qx3f)G
zd<j!Ne2Q<jJHR;C_cjL3E*J}Z^-8QM6-IJ8r7p}~gzm591Gkbcz|YN=ay+rAFe%G)
z_@zK3Oq`ki8F+jv{IWMLuQ*@~V>Ge*))BHWVt4mqu*pii9=|G|=A41af+59|KK3yF
zP`=B5wIXzW?>XQY77OE&3K^RZVK7!?zScyg3cl>Sv*ee~IT%=H^&rwz8%9vJO`UPd
zFuW?~4eOP&@cF#RK~anaevHQI)t%=cGi<S!?0qQ+lTYLiJ5N@^kM7>byEQ6c>}>O}
znI7)%$(Oy?ED=U7RYZj(^}vMsGtFr>HkcqCm)yN<0)~5^$*>B?KxeXA(c8f@FtKU*
zaJ02DjQBW&Gwb2~U!FJB5FriUsy-{NuEAhfqBx@{z74)SVLoCgbsC1bnytU)`$E6(
zPt&Z;)i7#e{aIhNAHGj}zihX;JAV%H`_+4lVhtfNkNq~=v0iv6_?`2-?|ZoKnlrE~
z?jGrKp()2y5DT<SuM44ui$Hbb)q-bzvZPb1?^#fSjzNvmuIn9phN1Ru0M8>d4l9hx
z)zFg@BHc)TR_S=)D<mjScxf>vL51$sgOVL5;l0zgCCk6LK)dhRc-?ner0yjVcBe!8
zA@m`ClYr0}Qf6RE-WLubEZ0A?iZ@OO3a>e78b0-gtC`Ft=Xl*onZpr*NsEM_)YJIM
z=I##C*>44uy`OYQPU#Zw#T8j0zNzBjWH}F9DR8{VI+G4}CA+3Q?npt-l%m-3BX*>&
z#Mb*yj8q__@2=^TWE@tZbHVg>@gZo-vY8syv4<-k)n<ge_LCkfWm@%!t%9y5tvzcy
zqDkoj-#40dm|(?u(a&iG8%fD09B$n_YyuBzuIaoBXd^}ZNX#<tCqvx>qeV-1=3o`0
zuX3A9E|H3hyDoH|-vh-hmTmi)BcLN_MbWuk&q=7miI+utd$C7*I~c2zx55+2Gn?-S
z7m#jX4~olp^r7(AQJ$=aUQl<%I{vZKFVfC+`8&9McajQ*+Al;N>?GYbauDJU+CnO8
z%4F6$kcR~+TPnLG&O*0&=-Pd=a*#h^zTNS)D-<X)oxGqc4E@H5d5K|LuqQdWGAq{C
zknXy=t8-M@V1;X2wq<zGz|(i_=<uO(Qq}ngud5O~P>SA`h;d1USN9pV`kfns0`|;b
zJKs!${fc|{#UFS<R%1-xr+W<$7OPYB)cYOj&Ux+&E55#i{C(-K-_`Db6szqew<az_
zTi(rO5nR`yt*+|;qb3U7>x#|}y>5e4+n;;TK@!}}xhQ%mry6Rw<6}yr&tgR*u1&LI
z`A{p??3|!AK&lMNy<%!b!*&Z<=B&2A3b%E+l^Y5tpsMf~|0I<aOA1Yw49<yw_yPUP
zB8!8tjQ-Eb2eDW1u=z*B`t>gGvLyVP%&;vK3ku3AkoQ972Nj)?dJ=qOSf6v_t1A5*
zV1{;wy2RyfgT86z8%rr|%<O8>luX?VO!mQ&VcRNiFc#*xc|)}Uv&ybE*xC?BdYK@Q
z<B^pFih(NOXwpW^*7c^mRfR6Lreb}hq=FA@iufKhYq$e;wipU+;vdCGJ~?D(KPt9e
zJdr#u-+?KYGDHk$-o}hSZg<-MR2GC|J8N5G?}5n@f9AcX%&`?~&7Uc|X2IqsOUpvJ
zguv?6!CS?zoI$H{fTfo+5p+|YC1o_$V+O+LEQ9zeu<FS`6^#Zk4I$DR_gns;K~m{C
zmh6TtcX?yM3+dSET@BM$7>zK!ipl%uhLS;jy$Quvr5<)BR-O)U*@>;vEL-K;SBvQl
ziBCT5Xu{C?ca+cTbHLQC<mZ{9GHlC0IO|K#<(NK8E;@4g2AI@2T?uV{0VprH9}`Fb
z-SeZY=ap81>~VudEOo*djbZzx?5r>_ee+=1X`4z+<#cn_KvoRM-j&L4?QO=i7)td|
z4@P6Qg>jL^ne5ociX4Fnp?GZb(7}mkOG`kbV>(H0<x5QM{dPM}UkyxY@P6xyo>0u>
zRD7iFP$Q<}p($x<69%f)LV-&4>oE=5eN!VRD@<-liO7o}UCh4A+$PPi3EOaFzfDPy
z73c(AIIVfh19Qv%S^2g|6n1G`FdJ{l!Su`zq9ocXK$qqXzESYO^tAg*-enEIx{!UW
zBu+PM6Z>UdSy~h*x!lj#vBwdEm*>yjNnegBgzVMOy_E|Z8OnP*HG*MVN{srCrIDZ}
zp&aGq5s1kMd@Cu;&w@>g0Ym5tbxijnX_O`N5$M;(4S0NX#mpW3PLxJ4VY2e+-z7LM
zVOzc=%TKtyz|aFF;x8+Iflk(E(|2csu=N@2Z-O~Lg6d}7-A@oV%#^p1KjW+=C}e**
zdtt+6OqoqNcX8WhFkN5$U?wFGQ;12-w+t4*EViB;T*b8!$U(G{$C@`WD^EQui>WwR
ze<!+3WT+mLp9HZ@*%@I1BlT-4(u?Qc2g3KwxH}JMly}K}x;O0?08%l(wB?@5!?JP(
zM-9(ljOkDd$Mg3JATzb6F2qm?WR5kx_&AD^CSzZ3@tf2J19nBL_NXJ}caFP9MW<|n
z<#z*BlOKF3e>wByp@b_l2v7?buQy@>4Ha(l+%=LQ$hVl{6Kw{v6~%l)x8g9>vT+q*
zKfm&4`^UaLOp(Kc*D^QK;_5*DXI1pO);pk`?0E9A9s|ZriWA*dZCF0@Ni3y4a6c?N
z<*4&ca;AK!NJGcnxq~!*E!{$c>mA6SI$XLst*QL&_h71to(m{6O>+tee+9lBLvK&v
z_{wJ*qt<_m!uXXPPrhvUfk{!<JzJ!q0}{DgyJp?5fDGf2FW(;hAdMea&PX~I4?>nP
zPVTN+z*Od*rh@TejDtPadrq~Jx;CAB+qQv`G{(!zZY1l0NtRv;VU0fuYV{G>CBu=V
zPb~@gPVe4;gz@H<_EKAn^_FDLv&E-?|5wd1nT%{q!SaBJje!@bTk9C#-U=1ays@~R
zHJYz{$Upm%^Oj|pLU`obRR;xNLq*xSS|kD^OO(}KInq_$!_^*e@Q@nDvyZ7xpMMCL
zmpv|HQ(-D^u<T!VCjhTMM_j|cRjQ=MA^#<XQJXR5F_zKib`c=In4IB}cCNg`?Rn1j
z-VRb<3}e8d4NS0faaUg`mm%=?^sEzJz6E3L=1RE1Wef5vUWI+;9)=|=H5>P=--a1V
zwzhZ;{{ZC!Rh<_ES+S*WOBp8GGB9Cd{(vxP8(=!-1rO|&fDrR{sXWtfprN1Er_J9A
zCiy+zR?`lZ|Gbnd>!2wK%LTXRh^04xAa!7K_7`tJMEo5}+g)HoC|Y~9;S7kmo(wo(
z#tLTR68B6JPh-eC3+hF&KGIl>@`ZcjD&=E>!Dd$vG?%~S|E^sjvz;{FnA^32qKHZF
z8agY{xdhln?+$Sb769DbWYwp17sOdZcG%T4V{+L0jr%Drz-t-hZXw_Xi!;^9&bx2U
zpMw|Pi>_QxSAd2!3EVYWLs0Yh*8BT0>X2mBnW2{F1GP`QzH1C8L-SpO1j(DLv8ugp
z=UR^FLt__*4`t*4RPC4?snlkKmu=G_>;jviERog4KI%0z3^IFex7-9zMO^D+Z`;DF
zufMEP`_16#oe)-@W!=!0`yP{iItVqRlA0@7wn4qtiZM#VS7<1$Z^>je#<CAa<^OuK
zACk$1F9&F@&=}XMt2NjHjgprYcPJ%b-TlFXH@Hti!<&4AEkbvoLHF7{xogkie6Ie>
z;vE7|{qDpDe%~Z0=Zm}F!>j|30;0z2Jxrm&{FjD_;{=q~{HW&1j=~;M%nU<2)<I*k
zgH>?o8fb}1l2P~j3XR_~J#0DO!BfSr#uiMb(5UfJN+M1P4jzh1nKqM#W;qj^+5Kd=
z9>x8&Qr!ZY_H~Mx)zw4Ix8RLg(+{!ACI}E7jfEEvL-q_kEQ7X1CYEt2T+k%M!6?70
z7FxJf&tgg_RBnv8HEFC3*(Mr&T$eULll`ffSiNa@#<sM+>t+G=k*}s|)ss+oV*Rnj
zSf>{n8o9PFTG0oM#w&hYb6E+^{Arbwj!uxDGyPiDxe=NVtZrc#dV_@{C$yw>)Uo6v
zCRwd<KG0aIUnHiY2{$7@SgGh0Lj9`eHm6qyL+{Tu@h8jTp>ZRPd*8iDXuK8Oa@&at
zuRTpZdpB7_tN3DV!S*O9Y9j3z7c7RtfyvL0#Ae{zm)X_Ey`fk|<4$qwiLFqt%YVu0
z_Ct846jUp02|$y|Dp=f(|4w@CoZ8wh4wVr*3MR5PL(?^$v>~^>@Hjky`}u}yc%sP3
z*7HLhD{Hgw?G!l<d5?eis%_eeJy7U6WJYs<#<zW{f%>BGhCdGC9cJLh(ge-6-TR?#
zNoK}<4ijjR+w1dl+j97veIwFpV+gc4s&{`gdV*Ct6(cj2kD<Or*2{2G0!C7|t%=!m
z7ar*)-cRzu=hCZA-tWxT`Ew98D1KX<Sp_Zx?mxsi%>aHW%DyH3SFuwCF-7?mU*W>;
z4cz;3uR<o{V++-dSFv>V{2Q`-S#W7O&>3ao$C9ow+lM(1K;lm&MfMMMSf<4{XQOUb
zINx$_-5KNz7@B6x=!<m2l@hzgjv5C@+1tGHb)F3bZCl)Z-7Xp;PW;TTeta6vJEZQl
z(UOC-MBch({LEN<!{c3#40mJ6R_~AK9czV)dS`Op=*+?;&Y*;!@rBs^)Qdp_5ntiF
z@1d>wm7m~J&<(#`XWn84UfUUlr>Jl~I?SPC))KRynobW<4u)h>=4^ld4T$!XlQ1i}
z39hpm;aAi7vGnkJS_hZsL#nQ}?-k~6kgi&tE*4FObN;82Up;e!18u{HW{g`QSxNo7
zd%-F!<YA%k-8D>*w#$P<TIvo2UG2W}sInd|X|=CRZ5V*a0*Q@7&9T^tg!C@u)kSbr
zXyRAz1|P`G_d0Ifrv_;*({2V_d5~anGXFyT1&Hsc<2k|k8B3U!-u@`d7c%4oM3wED
zAxf0JfAwAxRyd@&hwaNQh(F1FjTW5?X;QhCUQs4+vCK&FeC7b8``PS_qzz-)ud-#X
zS*by+pnQpi@hdEXby0_-VFw(unm%oxI|E_m&!khcZez*}(v2-Uyx??+>W96{CLx`}
zN<lD(53<CIbT2kF!^Lr%+P#goa4vSK>g$3)NV)w~Lh_jf7MAh-wpILHEPX6Rs_a`L
z<h9-1b-_p%?9^E6M*5FJWCM@f-G?j?5UCu+oEwTIesHV3aNHM5Xj`kLSlA26WBsAO
z+yfxqQj%jJcLTW8?RTF>^?|19pB|#7i(PeXFnAYjj=9+$+ue2aB9^~lbJMyZ7D!S3
zDa?G144FWEqhePL=Sy;Q2HzP1ZEWCG6ng<gKWkK8YT^#JpX87eoZiDVtBrv_pPj(+
z+E%Wk>2SjtCZUg2nwJ67$X!H<9fWhp#pu{8wvbS|Y2vHq^Z9cyuzyGTr>b^%ozQPO
zeRvRtt_|^N8+1bpRgft?{|+=?l?YU1{{^EV36oVHW}r1@YT}nlJG?2-=kC5H219||
zUuCy{h5m4#caLug!^7DRi(KSvAahR+7oSKZd=|)MtO)%MU6kO-rDrFh`ICjCu6`eM
z@Y#lx=@>!d7tuBCa%6ZJ$t`r-s|P-6rS|(j(St`fe`KCx6U5$oi)E?r;e*a1_ls$!
zfzW*H(567c$50(!pfRN|j+MPDjk1583f+zJ*PD)b!iy3a-JDHc&}Ex*Lant3Uf*^S
zS~S=WokLb>sgnKhE|Pzh0S6ltSAR6li<-diyQGh%HB&ISSW(aZSs`?uZ`3%_kP6Rt
zdTgH+I|}tr!^%FnL_tqT$!;Bc8k{~@jy3b4P{+(5{h;?M)L`M)^`>^h+eJG`!{#P1
zda!rrN7F#8l!<p!mz4-~sYaSXz&Ny#mRvdz^Bo%392#Ios-Z)2`QxN!JLqRFKS1)4
zhj#W!Au}mAXe|o0h~8ukWxI6Avt=nz<Pp=u!z%%w<LacgA6o|9BC%~A&56(*DEI1l
z^?vC2Ad+-+tPpC8rZ`1Cl%RL(nvKpaxzOo&aJ!wX2v$!um=fR!f-z}5C!>vKP|FZ1
z@z`oL^q5&n&0bc4JkBa7o1?+dy+81np;I0--g#;w%_9ljGKRChgURq*;)div)lO)!
zVOcuSO@%jaxt{XT;-M1WIS9N{g<Fk|fmhq6p!ZDukyw=>c#$H$OjxfO>VtWMR>>D(
z4^_Ur?>TM`Z-Y1AIe)DS`mtkyJvpbK>zvg<YB&BnA^Xny&}*{L!tjeHrRO*_B;D-f
z!}6fLcn$Z9sD03}<>xr#YD4(!x!57P<qW)4ICwr*_&a>GEPKlC;)7wEXR|xAj>2=M
zQl~AAH83G5mYY#M2G?$jSB<G1hJjNV1CLM0K*OH6Wv{Z|V0VTE_F4;_ho%jk8+6zm
z=;r{{^F6jIG_Mj>-66>2cs&4>y>5N`&buyDS~RUFC50W8Q;^)?;xmpas3fctZGVa`
zKBdN4-A6^24(v3|Zl|H5Hv^kFGpVSU%-X}|5Qz%BmGnDw@--^ZdGEv`(QK6WQpoOG
zW%tlUeOrauk&UQmQ}3lP?;KQUy^CW=l0M4lBU2<?UV-u%+1yw1%SE|$MzK~q29)Ks
zUVTm736!nS#h`?D2g?5*#2j8Vp}ed&W#3#+Kv@{G4eHcRqby^t@lxH_P&UR}q6W%k
z=pr~ME)h|P@*EC*E);HuGB_r?v2D9e?!bOHXm2h+Sx!!U{SlaoGR5X*-;G;~vTyUz
zJ;)k|vbDD+rl%>Be+&wSe!m4MlDB6-^U)%d?Z{<O?)<ytAD3#aj~G;uhr)Y}M~hLE
znI~hd2j6{^<yXF`$e<sJ$jjzMb#swtOgfggCEAh)s#Mj2cv8rdU8+{UwtXT`*IHVH
zm`tOL<;Sr-%a4$MCNR6Zo^>J*GfM`K{;(kbT1n-}$p|A)C(VjRt^7(Js?Ms#lnlw=
zH(rx{6gx~F^50QVB)AGiww8L^xE@Cl$HAp856#K_oFSj~SLc(5N){vc)Ur@UF8BO2
z(N^*(jmwz*aUuERMM1`#B4(5U`SJ5yR1NvZi!B32t)k?<kyDC!l?LRY%59RUKoz<F
zVd}XQqfGMncx|GVmIAr=>xnCrt%l^!vgb}}#J?a9MCeDyUl=C$?y=R`XWB&WxX5kQ
zGWnU@xoc2o!21BX&)?Vm>fn3wpmwwLttba_-wwA@olZ`2cYxn_o~2&ok&5NLH#eA*
zC--ISv12JGe`4bFTKTD&+#y2MxT5}r{I*{qD4$Y8?nzBkU|)We+>_L?H*{GgxqntF
z=ECw1<c?i>MxPAL$e)+gej2O6{Wv*Ing@o-Edu@zFW$UN?x?lgbS++j++|KWp-j0)
zZZqxUyu93p-1|~>>Cpfca(k!N*{-{i<c~Uoyi_nCck})-<#1-g&%y5xD}VntkKzAw
zLjIrP)9)jf!>&gpb#qX&4>?zxy?N0sH4fWyIfqf5yvy>DZ{tzD&UV&Jz7EtpESj0;
zjx)NJ-$&fOp$;W48fnko!;I==T{=pVY(!BE4HFPqj>-u&nlidjQT4FUu1^L%sCs$G
z)!PExsNBTsqzzJls_fIdo-*(Vl|Q5};#bm#%E?jNzZdwUN~+dVfhnA*O65w)hOc_4
zOpkmjGV}nIRTwDb)sIFM{a)nM?U6^7+-)U<s}7?gqJv}R4b12&E+0m1sZLZbGx=T{
zYX`cxSLTGgF21jGJwvgDBq}$3y7Y6k9eLQaEJz5cMHTiPP+|Y-fXbskIPKW-QHdu~
z*F$R8qDscq8M4jYsQ9o5l93gQF8-R$vi1@mD!F0VmG;<iRJz<PIFb7t`4y9E+D3)7
zsEqY?<-W!ZsK{-BpC9*}L8T9<f4`Eghe}1gDqXklIk{!{pk|oKD^#*dUSQ|hSagYN
z&+%oakE7zOcXiyqTA)iZ!#`&5^q_1pdr$5Yl|aR9svHbL*P&8jliYSI)KD(X^4iTJ
zuTlP`h-Iwj7gRXq^{s>XyHRNm{=wu;A*k##>+a?c5hzpW^I6lIk5ML7raN!GZA2x<
zJ)J|5%_zs3y=Ct`{K)u^II-Q4lc=1#ciB?Oa8z<{(~B3=OejazkznnfYLxTB^eN>N
zyyO|Jomv5>529jvub2k8+fc3#<>TrxN60gi6Vg|QT~UFt6U#YmPN0H%hPO8|lgU3$
z%T2~|I+1HLR<QNlN<#(MY?I{=twRNeYHO?CTA>`7IZhwWwWBPC&TmXdlu@CS%+!JC
z)2OiUosTOyFQc4pnpw@=vnX4j{f85i#pJg#iJz?<8OT4y1YcEoUP5_Cr|qp<`%nSH
z2RnYGSfiZo%`G9F>L`c6vYC`imMDL|UfVs6335HNR(dnJk39ONu+;EqF)Fm~fK#or
z3%T_LgYCHHJ^VfOeim_;AWu07hIiTqk%zk_@A!=RkiVyE7)QU_jUdd7f6l|-k6&15
z6ff%UkEZ#M-yb~xRWkeeIMTloW0?OeRsYYYcXlSY!<>pyPibCKo|jbSCDnOJZC+BJ
zmo(-j&3Q>{PGUk0_357%x%ugP)ArLBcAA;}Ck1w2??6Rj8uV{+86D_BF*zs`M3sM{
z>OWEKpQ!#%)c7ZA{u8zSiQ4p<;Sfp$S(zZK5M)(?tVWR439<%3)+ESU1et)VLcmoa
z;HnUCRS38$1Y8vYt_lHHg@CI<z*QmOsuFNj3Am~RTvY<DDgjrOfU8QtRVCo65^z-s
zxM~DkH3F_00auNHt46?8BjBnLaMcL7Y6M(00<Jm%SDk>XPQX<s;HndF)d{%j1YC6j
zt~vo%oq(%Bz||n&Y7lTW2)G&qTnz%Q1_4)tfU7~k)ga(%5^yyMxS9l9O#-ea0aufN
zt4YAsB;aZia5V|IS_E7z0<IPTSBrqFMZnb};A#<YwFtOc1Y9iwt~LQzn}DlLz||(;
rY7=m^)&G03I(kz4=r<Bw$>Butrzkr6{=PyHWafVWgif)+b29({ZYYlK

diff --git a/scripts/refit_constrained.R b/scripts/refit_constrained.R
new file mode 100644
index 0000000..33b3d84
--- /dev/null
+++ b/scripts/refit_constrained.R
@@ -0,0 +1,488 @@
+# =============================================================================
+# refit_constrained.R
+#
+# Companion to fvsRemodeled_CONUS.R. Refits Greg Johnson's CONUS dg and hg
+# equations under two improvements that address the issues he flagged in his
+# README:
+#
+#   (1) sign constrained nlsLM fits, so coefficients land in biologically
+#       expected directions (no negative crown ratio, no positive elevation
+#       penalty in the wrong direction, etc.);
+#   (2) optional pooled fits within Conifer / Hardwood functional groups,
+#       which stabilises species with marginal sample size and lets us
+#       borrow strength across closely related species.
+#
+# Inputs expected:
+#   - tree_dg_subset.RDS (FIA remeasurement training data for diameter growth)
+#   - tree_hg_subset.RDS (FIA remeasurement training data for height growth)
+#   - species.RDS        (species reference table with SPCD, CommonName, n)
+#
+# These live in Greg's repo under fvs_remodeling/rds/. The script is designed
+# to drop into Greg's own workflow without rerunning his data assembly.
+#
+# Author   Aaron Weiskittel
+# License  GPL 2.0 (matches upstream fvs_remodeling)
+# =============================================================================
+
+suppressPackageStartupMessages({
+    library(minpack.lm)
+})
+
+# Light-weight base R helpers used in place of dplyr / tibble so the script
+# runs on a vanilla R install with only minpack.lm. If you have tidyverse
+# loaded the code still works; these helpers just avoid the dependency.
+filter_eq <- function(df, col, val) df[!is.na(df[[col]]) & df[[col]] == val, , drop = FALSE]
+
+
+# -----------------------------------------------------------------------------
+# 1  Sign and magnitude bounds
+# -----------------------------------------------------------------------------
+# These bounds are deliberately loose. They block obviously wrong signs (e.g.
+# CR coefficient negative, BAL coefficient positive) but are wide enough to
+# let nlsLM find species variation. Tighten any bound you have a strong prior
+# for.
+
+dg_lower <- c(B0 = -10, B1 = -2,    B2 = -10,  B3 = 0.5,  B4 = 0.1, B5 = -0.01,  B6 = -0.1)
+dg_upper <- c(B0 =  5,  B1 =  0,    B2 =  0,   B3 = 5,    B4 = 2,   B5 =  0.01,  B6 =  0.5)
+# Rationale:
+#   B1 < 0   (size penalty in log term)
+#   B2 < 0   (BAL competition reduces growth)
+#   B4 > 0   (BAL effect is monotone increasing)
+#   B6 sign for EMT free; some species respond positively to warmer minima.
+
+hg_lower <- c(B1 = 0.001,  B2 = 0.5, B3 = 0,    B4 = 0,      B5 = -0.001, B6 = -1, B7 = -0.5, B8 = 0)
+hg_upper <- c(B1 = 0.5,    B2 = 5,   B3 = 5,    B4 = 0.05,   B5 =  0.01,  B6 =  3, B7 =  1,   B8 = 5)
+# Rationale:
+#   B1 > 0   (Chapman Richards rate)
+#   B3 >= 0  (positive crown ratio effect)
+#   B4 >= 0  (CCFL competition reduces growth)
+#   B8 >= 0  (CCH at tip reduces growth)
+
+
+# -----------------------------------------------------------------------------
+# 2  Equation forms (annual increment, summed to period via integrated fit)
+# -----------------------------------------------------------------------------
+# These mirror Greg's est_dg and est_hg verbatim so that any species whose
+# unconstrained fit converged in his original .qmd will also converge here
+# given the same data. The only change is that the wrapper passes lower and
+# upper bounds to nlsLM via control = nls.lm.control().
+
+est_dg <- function(n, dbh0, bal0, bal1, cr0, cr1, ht0, ht1, elev, emt,
+                   B0, B1, B2, B3, B4, B5, B6) {
+    dht  <- (ht1  - ht0)  / n
+    dbal <- (bal1 - bal0) / n
+    dcr  <- (cr1  - cr0)  / n
+    cht  <- ht0; cdbh <- dbh0; cbal <- bal0; ccr <- cr0
+    for (i in seq_len(max(n))) {
+        dg_hat <- exp(B0 +
+                      B1 * log((cdbh + 1)^2 / (ccr * cht + 1)^B3) +
+                      B2 * cbal^B4 / log(cdbh + 2.7) +
+                      B5 * elev + B6 * emt)
+        cdbh <- cdbh + ifelse(i <= n, dg_hat, 0)
+        cbal <- cbal + ifelse(i <= n, dbal,   0)
+        cht  <- cht  + ifelse(i <= n, dht,    0)
+        ccr  <- ccr  + ifelse(i <= n, dcr,    0)
+    }
+    cdbh
+}
+
+est_hg <- function(periods, ht, cr, cr2, ccfl, ccfl2, cch, cch2,
+                   max_height, elev, td, emt,
+                   B1, B2, B3, B4, B5, B6, B7, B8) {
+    htc <- ht; crc <- cr; cchc <- cch; ccflc <- ccfl
+    crg   <- (cr2   - cr)   / periods
+    ccflg <- (ccfl2 - ccfl) / periods
+    cchg  <- (cch2  - cch)  / periods
+    for (i in seq_len(max(periods))) {
+        inc <- max_height * B1 * B2 * (crc)^B3 *
+               exp(-B1 * htc - B4 * ccflc - B8 * cchc^0.5 -
+                   B5 * elev + B6 * td^0.5 + B7 * emt) *
+               (1 - exp(-B1 * htc))^(B2 - 1)
+        htc   <- htc   + ifelse(i <= periods, inc,   0)
+        crc   <- crc   + ifelse(i <= periods, crg,   0)
+        ccflc <- ccflc + ifelse(i <= periods, ccflg, 0)
+        cchc  <- cchc  + ifelse(i <= periods, cchg,  0)
+    }
+    htc
+}
+
+
+# -----------------------------------------------------------------------------
+# 3  Constrained per species fits
+# -----------------------------------------------------------------------------
+
+fit_dg_species <- function(d, start = NULL,
+                           lower = dg_lower, upper = dg_upper,
+                           trace = FALSE) {
+    if (is.null(start)) {
+        start <- c(B0 = -1.29, B1 = -0.5836578, B2 = -0.0953406,
+                   B3 = 1.83,  B4 = 0.62814,    B5 = -0.001023,
+                   B6 = 0.03985)
+    }
+    minpack.lm::nlsLM(
+        d$endDIA ~ est_dg(d$endGROWYR - d$startGROWYR,
+                          d$startDIA, d$startBAL, d$endBAL,
+                          d$startCR,  d$endCR,
+                          d$startHT,  d$endHT,
+                          d$ELEV,     d$EMT,
+                          B0, B1, B2, B3, B4, B5, B6),
+        data    = d,
+        start   = start,
+        lower   = lower,
+        upper   = upper,
+        control = minpack.lm::nls.lm.control(maxiter = 200),
+        trace   = trace
+    )
+}
+
+fit_hg_species <- function(d, max_height, start = NULL,
+                           lower = hg_lower, upper = hg_upper,
+                           trace = FALSE) {
+    if (is.null(start)) {
+        start <- c(B1 = 0.0108871, B2 = 1.6261262, B3 = 0.4573120,
+                   B4 = 0.0014822, B5 = 0.0001514, B6 = 0.2736003,
+                   B7 = 0.0326502, B8 = 0.866)
+    }
+    minpack.lm::nlsLM(
+        d$endACTUALHT ~ est_hg(d$endGROWYR - d$startGROWYR,
+                               d$startACTUALHT, d$startCR,   d$endCR,
+                               d$startCCFL,     d$endCCFL,
+                               d$startCCH,      d$endCCH,
+                               max_height, d$ELEV, d$TD, d$EMT,
+                               B1, B2, B3, B4, B5, B6, B7, B8),
+        data    = d,
+        start   = start,
+        lower   = lower,
+        upper   = upper,
+        control = minpack.lm::nls.lm.control(maxiter = 200),
+        trace   = trace
+    )
+}
+
+
+# -----------------------------------------------------------------------------
+# 4  Pooled functional group fits
+# -----------------------------------------------------------------------------
+# Pool conifer or hardwood species and fit one shared coefficient vector. The
+# resulting parameters become the fallback that fvsRemodeled_CONUS.R returns
+# when an unfitted SPCD lands in resolve_species() tier "group". By fitting
+# rather than averaging, we let the data speak instead of imputing arithmetic
+# means of unconstrained per species fits.
+
+fit_dg_pooled <- function(trees, group = c("Conifer", "Hardwood"), ...) {
+    group <- match.arg(group)
+    d <- if (group == "Conifer") trees[trees$SPCD <  300, , drop = FALSE]
+         else                    trees[trees$SPCD >= 300, , drop = FALSE]
+    if (nrow(d) < 5000) {
+        warning("Pooled fit has only ", nrow(d), " rows; results may be unstable.")
+    }
+    fit_dg_species(d, ...)
+}
+
+fit_hg_pooled <- function(trees, group = c("Conifer", "Hardwood"), ...) {
+    group <- match.arg(group)
+    d <- if (group == "Conifer") trees[trees$SPCD <  300, , drop = FALSE]
+         else                    trees[trees$SPCD >= 300, , drop = FALSE]
+    if (nrow(d) < 5000) {
+        warning("Pooled fit has only ", nrow(d), " rows; results may be unstable.")
+    }
+    max_height <- max(d$endACTUALHT, na.rm = TRUE)
+    fit_hg_species(d, max_height = max_height, ...)
+}
+
+
+# -----------------------------------------------------------------------------
+# 5  Driver: refit all species with sign constraints + pooled fallbacks
+# -----------------------------------------------------------------------------
+
+refit_all_constrained <- function(rds_path, min_obs = 5000,
+                                  out_path = NULL, trace = FALSE) {
+
+    species   <- readRDS(file.path(rds_path, "species.RDS"))
+    trees_dg  <- readRDS(file.path(rds_path, "tree_dg_subset.RDS"))
+    trees_hg  <- readRDS(file.path(rds_path, "tree_hg_subset.RDS"))
+
+    # Per species DG fits
+    dg_fits <- list(); dg_pars <- NULL
+    for (i in seq_len(nrow(species))) {
+        if (species$n[i] < min_obs) next
+        d <- filter_eq(trees_dg, "SPCD", species$SPCD[i])
+        if (nrow(d) < min_obs) next
+        f <- try(fit_dg_species(d, trace = trace), silent = TRUE)
+        if (inherits(f, "try-error")) {
+            warning("DG species ", species$SPCD[i], " failed to converge")
+            next
+        }
+        dg_fits[[as.character(species$SPCD[i])]] <- f
+        p <- f$m$getPars()
+        dg_pars <- rbind(dg_pars, data.frame(
+            spcd = species$SPCD[i], n = nrow(d),
+            Common_Name = species$CommonName[i],
+            B0 = p[1], B1 = p[2], B2 = p[3], B3 = p[4],
+            B4 = p[5], B5 = p[6], B6 = p[7],
+            AIC = AIC(f), isConv = f$convInfo$isConv,
+            RSS = f$m$deviance()
+        ))
+        cat(sprintf("DG species %5d (%-22s) n=%6d   converged=%s\n",
+                    species$SPCD[i], species$CommonName[i],
+                    nrow(d), f$convInfo$isConv))
+    }
+
+    # Per species HG fits
+    hg_fits <- list(); hg_pars <- NULL
+    for (i in seq_len(nrow(species))) {
+        if (species$n[i] < min_obs) next
+        d <- filter_eq(trees_hg, "SPCD", species$SPCD[i])
+        if (nrow(d) < min_obs) next
+        max_h <- max(d$endACTUALHT, na.rm = TRUE)
+        f <- try(fit_hg_species(d, max_height = max_h, trace = trace),
+                 silent = TRUE)
+        if (inherits(f, "try-error")) {
+            warning("HG species ", species$SPCD[i], " failed to converge")
+            next
+        }
+        hg_fits[[as.character(species$SPCD[i])]] <- f
+        p <- f$m$getPars()
+        hg_pars <- rbind(hg_pars, data.frame(
+            spcd = species$SPCD[i], n = nrow(d),
+            Common_Name = species$CommonName[i],
+            B0 = max_h, B1 = p[1], B2 = p[2], B3 = p[3],
+            B4 = p[4], B5 = p[5], B6 = p[6], B7 = p[7], B8 = p[8],
+            AIC = AIC(f), isConv = f$convInfo$isConv,
+            RSS = f$m$deviance()
+        ))
+        cat(sprintf("HG species %5d (%-22s) n=%6d   converged=%s\n",
+                    species$SPCD[i], species$CommonName[i],
+                    nrow(d), f$convInfo$isConv))
+    }
+
+    # Pooled functional group fits used as the resolve_species() group tier
+    pooled <- list(
+        dg_conifer  = try(fit_dg_pooled(trees_dg, "Conifer"),  silent = TRUE),
+        dg_hardwood = try(fit_dg_pooled(trees_dg, "Hardwood"), silent = TRUE),
+        hg_conifer  = try(fit_hg_pooled(trees_hg, "Conifer"),  silent = TRUE),
+        hg_hardwood = try(fit_hg_pooled(trees_hg, "Hardwood"), silent = TRUE)
+    )
+
+    out <- list(
+        dg_pars   = dg_pars,
+        hg_pars   = hg_pars,
+        dg_fits   = dg_fits,
+        hg_fits   = hg_fits,
+        pooled    = pooled,
+        bounds    = list(dg_lower = dg_lower, dg_upper = dg_upper,
+                         hg_lower = hg_lower, hg_upper = hg_upper),
+        timestamp = Sys.time()
+    )
+
+    if (!is.null(out_path)) {
+        saveRDS(out, out_path)
+        cat("Wrote refit bundle to:", out_path, "\n")
+    }
+
+    invisible(out)
+}
+
+
+# -----------------------------------------------------------------------------
+# 6  Diagnostic: did sign constraints actually bind?
+# -----------------------------------------------------------------------------
+# After a refit, check which species hit a bound. A coefficient that lands
+# exactly on a bound is a flag that the unconstrained fit wanted to go the
+# wrong way, which is exactly the issue Greg's README calls out.
+
+bound_diagnostic <- function(refit, eps = 1e-6) {
+    check <- function(pars, lower, upper) {
+        coefs <- pars[, intersect(names(pars), names(lower)), drop = FALSE]
+        flagged <- list()
+        for (b in names(coefs)) {
+            at_lower <- which(abs(coefs[[b]] - lower[b]) < eps)
+            at_upper <- which(abs(coefs[[b]] - upper[b]) < eps)
+            if (length(at_lower)) flagged[[paste0(b, "_at_lower")]] <-
+                pars$Common_Name[at_lower]
+            if (length(at_upper)) flagged[[paste0(b, "_at_upper")]] <-
+                pars$Common_Name[at_upper]
+        }
+        flagged
+    }
+    list(
+        dg = check(refit$dg_pars, refit$bounds$dg_lower, refit$bounds$dg_upper),
+        hg = check(refit$hg_pars, refit$bounds$hg_lower, refit$bounds$hg_upper)
+    )
+}
+
+
+# -----------------------------------------------------------------------------
+# 7  Entry point
+# -----------------------------------------------------------------------------
+# Edit rds_path to point at your fvs_remodeling clone, then source this file
+# and call refit_all_constrained(). Example:
+#
+#   refit <- refit_all_constrained(
+#       rds_path = "F:/Projects/FVS Remodel/fvs_remodeling/rds",
+#       out_path = "F:/Projects/FVS Remodel/fvs_remodeling/rds/refit_constrained.RDS")
+#   diag <- bound_diagnostic(refit)
+#   str(diag)
+#
+# The resulting RDS is a drop in replacement for dg_parms.RDS / hg_parms.RDS
+# in fvsRemodeled_CONUS.R; just point load_gj_params(local_path = ...) at the
+# directory containing the new file (after splitting out dg_pars and hg_pars
+# under the original filenames).
+
+# -----------------------------------------------------------------------------
+# 8  Synthetic training data generator (smoke test path)
+# -----------------------------------------------------------------------------
+# Greg's tree_dg_subset.RDS and tree_hg_subset.RDS are NOT in the public
+# fvs_remodeling repo (they are local to him). For anyone who wants to
+# exercise this script before pulling those files, the function below
+# generates a small synthetic training set with known parameters. The
+# refit should converge close to the parameters used to simulate, with
+# constrained signs always satisfied. This is useful for confirming the
+# bounds and the integrated fitting loop work end to end.
+
+simulate_training <- function(n_per_species = 6000,
+                              spcds = c(202, 131, 316),
+                              years_min = 5, years_max = 12,
+                              seed = 42) {
+    set.seed(seed)
+
+    sim_one <- function(spcd) {
+        # "True" parameters loosely modelled on Greg's published priors.
+        true_dg <- c(B0 = -1.4, B1 = -0.55, B2 = -0.10, B3 = 1.85,
+                     B4 = 0.62, B5 = -0.001, B6 = 0.04)
+        true_hg <- c(B1 = 0.011, B2 = 1.6, B3 = 0.45, B4 = 0.0015,
+                     B5 = 0.00015, B6 = 0.27, B7 = 0.033, B8 = 0.86)
+
+        # Draw covariates that span realistic CONUS conditions.
+        n <- n_per_species
+        startDIA <- exp(runif(n, log(2),  log(30)))
+        startHT  <- 4.5 + 1.2 * startDIA + rnorm(n, sd = 5)
+        startHT  <- pmax(startHT, 8)
+        startCR  <- pmin(pmax(rnorm(n, 0.55, 0.15), 0.1), 0.95)
+        startBAL <- runif(n, 0, 200)
+        endBAL   <- pmax(0, startBAL + rnorm(n, 0, 5))
+        endCR    <- pmin(pmax(startCR + rnorm(n, 0, 0.05), 0.05), 0.95)
+        endHT    <- startHT  # filled below
+        startCCFL <- runif(n, 0, 250)
+        endCCFL   <- pmax(0, startCCFL + rnorm(n, 0, 5))
+        startCCH  <- pmin(pmax(rnorm(n, 0.5, 0.2), 0.05), 1)
+        endCCH    <- pmin(pmax(startCCH + rnorm(n, 0, 0.02), 0.05), 1)
+        ELEV <- runif(n, 200, 6000)
+        TD   <- runif(n, 12, 28)
+        EMT  <- runif(n, -25, 0)
+
+        years <- sample(years_min:years_max, n, replace = TRUE)
+        startGROWYR <- rep(2005, n)
+        endGROWYR   <- startGROWYR + years
+
+        # Forward simulate end DIA via the integrated DG equation.
+        endDIA <- numeric(n)
+        for (i in seq_len(n)) {
+            d  <- startDIA[i]; b <- startBAL[i]; c_ <- startCR[i]
+            h  <- startHT[i]
+            db <- (endBAL[i] - startBAL[i]) / years[i]
+            dc <- (endCR[i]  - startCR[i])  / years[i]
+            for (k in seq_len(years[i])) {
+                inc <- exp(true_dg["B0"] +
+                           true_dg["B1"] * log((d + 1)^2 / (c_ * h + 1)^true_dg["B3"]) +
+                           true_dg["B2"] * b^true_dg["B4"] / log(d + 2.7) +
+                           true_dg["B5"] * ELEV[i] +
+                           true_dg["B6"] * EMT[i])
+                d <- d + inc
+                b <- b + db
+                c_ <- c_ + dc
+            }
+            endDIA[i] <- d + rnorm(1, sd = 0.1)
+        }
+
+        # Forward simulate end ACTUALHT via the integrated HG equation.
+        max_height <- 280
+        endACTUALHT <- numeric(n)
+        for (i in seq_len(n)) {
+            ht <- startHT[i]; cr_ <- startCR[i]
+            ccfl <- startCCFL[i]; cch <- startCCH[i]
+            dcr <- (endCR[i] - startCR[i]) / years[i]
+            dccfl <- (endCCFL[i] - startCCFL[i]) / years[i]
+            dcch <- (endCCH[i] - startCCH[i]) / years[i]
+            for (k in seq_len(years[i])) {
+                inc <- max_height * true_hg["B1"] * true_hg["B2"] *
+                       cr_^true_hg["B3"] *
+                       exp(-true_hg["B1"] * ht - true_hg["B4"] * ccfl -
+                           true_hg["B8"] * sqrt(max(cch, 0)) -
+                           true_hg["B5"] * ELEV[i] +
+                           true_hg["B6"] * sqrt(max(TD[i], 0)) +
+                           true_hg["B7"] * EMT[i]) *
+                       (1 - exp(-true_hg["B1"] * ht))^(true_hg["B2"] - 1)
+                ht <- ht + inc
+                cr_ <- cr_ + dcr
+                ccfl <- ccfl + dccfl
+                cch  <- cch + dcch
+            }
+            endACTUALHT[i] <- ht + rnorm(1, sd = 0.5)
+        }
+
+        data.frame(
+            SPCD = spcd,
+            startDIA = startDIA, endDIA = endDIA,
+            startHT = startHT, endHT = endHT,
+            startACTUALHT = startHT, endACTUALHT = endACTUALHT,
+            startCR  = startCR,  endCR  = endCR,
+            startBAL = startBAL, endBAL = endBAL,
+            startCCFL = startCCFL, endCCFL = endCCFL,
+            startCCH  = startCCH,  endCCH  = endCCH,
+            ELEV = ELEV, TD = TD, EMT = EMT,
+            startGROWYR = startGROWYR, endGROWYR = endGROWYR
+        )
+    }
+
+    do.call(rbind, lapply(spcds, sim_one))
+}
+
+# Run a smoke test: simulate, refit constrained, compare recovered
+# coefficients to truth. This is what gets executed when you
+# Rscript refit_constrained.R from the command line.
+
+run_smoke_test <- function(n_per_species = 6000) {
+    cat("Generating synthetic training data (n =", n_per_species, "per species)...\n")
+    sim_dg <- simulate_training(n_per_species)
+    sim_hg <- sim_dg  # same rows used for hg fitting
+
+    cat("Fitting DG (SPCD 202, constrained)...\n")
+    d <- filter_eq(sim_dg, "SPCD", 202)
+    f_dg <- fit_dg_species(d, trace = FALSE)
+    cat("DG fitted coefs:\n"); print(round(f_dg$m$getPars(), 4))
+    cat("DG converged    :", f_dg$convInfo$isConv, "\n")
+    cat("\n[Truth was]      B0 = -1.40  B1 = -0.55  B2 = -0.10  B3 = 1.85",
+        " B4 = 0.62  B5 = -0.001  B6 = 0.04\n\n")
+
+    cat("Fitting HG (SPCD 202, constrained)...\n")
+    f_hg <- fit_hg_species(d, max_height = max(d$endACTUALHT), trace = FALSE)
+    cat("HG fitted coefs:\n"); print(round(f_hg$m$getPars(), 4))
+    cat("HG converged    :", f_hg$convInfo$isConv, "\n")
+    cat("\n[Truth was]      B1 = 0.0110  B2 = 1.60  B3 = 0.45  B4 = 0.0015",
+        " B5 = 0.00015  B6 = 0.27  B7 = 0.033  B8 = 0.86\n")
+
+    invisible(list(dg = f_dg, hg = f_hg))
+}
+
+
+# -----------------------------------------------------------------------------
+# 9  Entry points
+# -----------------------------------------------------------------------------
+# When sourced interactively, just announce. When run as a script, fire the
+# smoke test against synthetic data so the script demonstrates it works end
+# to end without needing Greg's private RDS files.
+
+if (interactive()) {
+    cat("refit_constrained.R loaded. Options:\n")
+    cat("  refit_all_constrained(rds_path = '...')   # against real training data\n")
+    cat("  run_smoke_test()                            # against synthetic data\n")
+} else if (!is.null(sys.frames()) || identical(commandArgs()[1], "RStudio")) {
+    # noop in IDE
+} else if (length(commandArgs(trailingOnly = TRUE)) == 0) {
+    run_smoke_test()
+}
+
+# =============================================================================
+# End of refit_constrained.R
+# =============================================================================