From d774a9e9c7affa68f314076bc5a5fb225aef6a7a Mon Sep 17 00:00:00 2001
From: Nicholas Clark <nicholas.j.clark1214@gmail.com>
Date: Tue, 2 Jan 2024 15:11:57 +1000
Subject: [PATCH] better monotonic doc; allow for by variables

---
 R/monotonic.R                   |  31 ++++++++++-
 R/mvgam_formulae.R              |  11 +++-
 R/stan_utils.R                  |   3 +-
 man/monotonic.Rd                |  29 ++++++++++
 man/mvgam_formulae.Rd           |  11 +++-
 src/mvgam.dll                   | Bin 1061888 -> 1061888 bytes
 tests/testthat/Rplots.pdf       | Bin 22911 -> 22896 bytes
 tests/testthat/test-monotonic.R |  96 ++++++++++++++++++++++++++++++++
 8 files changed, 173 insertions(+), 8 deletions(-)
 create mode 100644 tests/testthat/test-monotonic.R

diff --git a/R/monotonic.R b/R/monotonic.R
index b1573a96..5d95634a 100644
--- a/R/monotonic.R
+++ b/R/monotonic.R
@@ -82,6 +82,35 @@
 #'                  points = 0.5)
 #'
 #' plot(mod2, type = 'smooth', realisations = TRUE)
+#'
+#' # 'by' terms that produce a different smooth for each level of the 'by'
+#' # factor are also allowed
+#' set.seed(123123)
+#' x <- runif(80) * 4 - 1
+#' x <- sort(x)
+#'
+#' # Two different monotonic smooths, one for each factor level
+#' f <- exp(4 * x) / (1 + exp(4 * x))
+#' f2 <- exp(3.5 * x) / (1 + exp(3 * x))
+#' fac <- c(rep('a', 80), rep('b', 80))
+#' y <- c(f + rnorm(80) * 0.1,
+#'        f2 + rnorm(80) * 0.2)
+#' plot(x, y[1:80])
+#' plot(x, y[81:160])
+#'
+#' # Gather all data into a data.frame, including the factor 'by' variable
+#' mod_data <- data.frame(y, x, fac = as.factor(fac))
+#' mod_data$time <- 1:NROW(mod_data)
+#'
+#' # Fit a model with different smooths per factor level
+#' mod <- mvgam(y ~ s(x, bs = 'moi', by = fac, k = 8),
+#'              data = mod_data,
+#'              family = gaussian())
+#'
+#' # Visualise the different monotonic functions
+#' plot_predictions(mod, condition = c('x', 'fac', 'fac'),
+#'                  points = 0.5)
+#' plot(mod, type = 'smooth', realisations = TRUE)
 #' }
 smooth.construct.moi.smooth.spec <- function(object, data, knots){
 
@@ -351,7 +380,7 @@ add_mono_model_file = function(model_file,
                                    '; // monotonic basis coefficient indices\n'))
     mono_idx_data <- list(coef_indices)
     names(mono_idx_data) <- paste0('b_idx_',
-                                   mono_names_clean)
+                                   mono_names_clean[covariate])
     model_data <- append(model_data, mono_idx_data)
   }
 
diff --git a/R/mvgam_formulae.R b/R/mvgam_formulae.R
index d9889046..a16ff2f0 100644
--- a/R/mvgam_formulae.R
+++ b/R/mvgam_formulae.R
@@ -13,11 +13,15 @@
 #' \cr
 #' \cr
 #' The formulae supplied to \code{\link{mvgam}} are exactly like those supplied to
-#' \code{\link{glm}} except that smooth terms, \code{\link[mgcv]{s}},
+#' \code{\link{glm}} except that smooth terms,
+#' \code{\link[mgcv]{s}},
 #' \code{\link[mgcv]{te}},
 #' \code{\link[mgcv]{ti}} and
 #' \code{\link[mgcv]{t2}},
-#' time-varying effects using \code{\link{dynamic}}, as well as
+#' time-varying effects using \code{\link{dynamic}},
+#' monotonically increasing (using `s(x, bs = 'moi')`) or
+#' or decreasing splines (using `s(x, bs = 'mod')`; see \code{\link{monotonic}} for
+#' details), as well as
 #' Gaussian Process functions using \code{\link[brms]{gp}},
 #' can be added to the right hand side (and \code{.} is not supported in \code{mvgam} formulae).
 #' \cr
@@ -31,7 +35,8 @@
 #' \code{\link[mgcv]{jagam}},
 #' \code{\link[mgcv]{gam}},
 #' \code{\link[mgcv]{s}},
-#' \code{\link[stats]{formula}}
+#' \code{\link[stats]{formula}},
+#' \code{\link{monotonic}}
 #' @author Nicholas J Clark
 #' @name mvgam_formulae
 NULL
diff --git a/R/stan_utils.R b/R/stan_utils.R
index a005ff87..60698f76 100644
--- a/R/stan_utils.R
+++ b/R/stan_utils.R
@@ -2613,7 +2613,8 @@ add_trend_predictors = function(trend_formula,
       model_data <- append(model_data, idx_data)
 
       idx_lines <- c(grep('int trend_idx', trend_model_file),
-                     grep('// gp basis coefficient indices', trend_model_file))
+                     grep('// gp basis coefficient indices', trend_model_file),
+                     grep('// monotonic basis coefficient indices', trend_model_file))
       model_file[min(grep('data {', model_file, fixed = TRUE))] <-
         paste0('data {\n',
                paste(trend_model_file[idx_lines],
diff --git a/man/monotonic.Rd b/man/monotonic.Rd
index 90668ffc..f65b39e7 100644
--- a/man/monotonic.Rd
+++ b/man/monotonic.Rd
@@ -97,6 +97,35 @@ plot_predictions(mod2,
                  points = 0.5)
 
 plot(mod2, type = 'smooth', realisations = TRUE)
+
+# 'by' terms that produce a different smooth for each level of the 'by'
+# factor are also allowed
+set.seed(123123)
+x <- runif(80) * 4 - 1
+x <- sort(x)
+
+# Two different monotonic smooths, one for each factor level
+f <- exp(4 * x) / (1 + exp(4 * x))
+f2 <- exp(3.5 * x) / (1 + exp(3 * x))
+fac <- c(rep('a', 80), rep('b', 80))
+y <- c(f + rnorm(80) * 0.1,
+       f2 + rnorm(80) * 0.2)
+plot(x, y[1:80])
+plot(x, y[81:160])
+
+# Gather all data into a data.frame, including the factor 'by' variable
+mod_data <- data.frame(y, x, fac = as.factor(fac))
+mod_data$time <- 1:NROW(mod_data)
+
+# Fit a model with different smooths per factor level
+mod <- mvgam(y ~ s(x, bs = 'moi', by = fac, k = 8),
+             data = mod_data,
+             family = gaussian())
+
+# Visualise the different monotonic functions
+plot_predictions(mod, condition = c('x', 'fac', 'fac'),
+                 points = 0.5)
+plot(mod, type = 'smooth', realisations = TRUE)
 }
 }
 \references{
diff --git a/man/mvgam_formulae.Rd b/man/mvgam_formulae.Rd
index 929b7794..de3d1019 100644
--- a/man/mvgam_formulae.Rd
+++ b/man/mvgam_formulae.Rd
@@ -21,11 +21,15 @@ are non-identifiable.
 \cr
 \cr
 The formulae supplied to \code{\link{mvgam}} are exactly like those supplied to
-\code{\link{glm}} except that smooth terms, \code{\link[mgcv]{s}},
+\code{\link{glm}} except that smooth terms,
+\code{\link[mgcv]{s}},
 \code{\link[mgcv]{te}},
 \code{\link[mgcv]{ti}} and
 \code{\link[mgcv]{t2}},
-time-varying effects using \code{\link{dynamic}}, as well as
+time-varying effects using \code{\link{dynamic}},
+monotonically increasing (using \code{s(x, bs = 'moi')}) or
+or decreasing splines (using \code{s(x, bs = 'mod')}; see \code{\link{monotonic}} for
+details), as well as
 Gaussian Process functions using \code{\link[brms]{gp}},
 can be added to the right hand side (and \code{.} is not supported in \code{mvgam} formulae).
 \cr
@@ -41,7 +45,8 @@ extensive documentation for the \code{mgcv} package.
 \code{\link[mgcv]{jagam}},
 \code{\link[mgcv]{gam}},
 \code{\link[mgcv]{s}},
-\code{\link[stats]{formula}}
+\code{\link[stats]{formula}},
+\code{\link{monotonic}}
 }
 \author{
 Nicholas J Clark
diff --git a/src/mvgam.dll b/src/mvgam.dll
index fa4f2bc453e333247240d6d702357d582ff8c023..85f8f1edec9c32d64c263a30562ab9d06267b30a 100644
GIT binary patch
delta 76
zcmWN|yA4240EW@H-tQR07A&B*%~kw~sA%mWX5s7Q&Y<G-$GXnCZk<iPZ_D66K{Fvd
fVT3ON5sFB}A`z*`L@o+Zit3<lufLwrUBBZ8O=}}5

delta 76
zcmV-S0JHyq=s<wzK#+(6HGPqaP6N)$5P{o;+X02!0)^WHh1&#$+XaQ&28G)Ph1&>)
i+X;o+3WeJXh1(2;+YN==4z=44DcJ%weYXbM5BgbZ7bAoK

diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf
index c63ab130434970a2731cdf47b4f5c2cab91ae0db..1b6f133121f9d6df6fb594dbbb63f9b0aaa76949 100644
GIT binary patch
delta 9634
zcmZvAc|4SD`##w!N{gjJ>M2TDBKue>Axk0a*t2F|!eqV46RGTF?1~TukHJi1Fg#+&
zmTeeg#$;zKiLo!gTlKu}^Lf8NAAj6)f9AZd>pafmJkH~qeI(;|Bx8aii}VdC8R<k-
z-d_(D7K1{JvdM~!!a%%un-94jXhvV=bv3E4M4fVHWOg^pW~@EIETN$m!<_L<o^@wk
zz-R2PgVFkl^pjOd&OTT#2O|fK`mw~o;+op9?Q7+#wR`l1CF}jIbL-m!HYIiIdvp1?
z@-5er+W4)lmG$y{`c|^a>RJYFx1fG+HX1#;^Idg+t=tA!0Xj+aFss4q`#V!!d)PN+
z%Ud`vV4t&`MvwQJq+4vRN27B&YxhVQ)kfNZqcnPBRB-9FRs0|=LxQt)>B0)rc9`IL
zU_<Nbtt?vI+o$0ULGS3xVbuo@;pUr6&^Fv2Po5;<Z2VH#M&k0uH>_&i?i9}C?5r0M
zooiaO1+0|M^J)hL+N>ChOYbk>X{NwH*zVREr`P_Zh4f0_2RXk;zg=Rlzx7gRn%^d|
z_kLI>U+-Wy&*Eg*63gE`QkuB+_27CH%8%ilaSQnVi_GLvTgN?Y!KK+GBSc=EjSj2P
zdU!F2j%xk>S#r|;OTuLK(#pe4-2TA+#>zUdzp=d(wh!!0>Z#U-t<y}g*0?RWDlTke
zuN%kHxV$y6gxmiHz-{)MXrxis)BBJ3Jm~9TKA9@kq4R2XM^ra8IybG`#E~<!Z^^ox
z{WaAWfdde%c0+U6?oF`;^xS^Q_F9)y=}nx<2ev(?X)6R-sNJWiGZ@jLG_5KQm^<Xh
zr8ok9?0e7YChfuPe^>8wd<(|GMMtLecs-@2VLXf5h3qgW=XLh(p>5T+oSm7f1>%R;
z0UltkYC*^JzRb<#WfWYz=^CB4V$BQICuAs55jG~^u#oXmZRU8;+)EKd*-3r#H>29{
zb<fq(f!P(_3EJ1Kmwq%>vimb2|B_UPwmnNqaNSa_A3nWofnWoy-Zz|n>h%~9?T6Vc
zp^Gh?@KtU4(n>h-L9g*mnfIe9HjlmJW=tJZ<!1QkS~vd+b9FLOon_rBWO`2*gYxjs
z#BO;)N3Pklhlkk(M-$3l>EZRF^2NRmeXeeQdb_?+)MyqpFzO+jK6DSrlN3-#O+@Wb
zLps*qkLSgy&sAy}_7<pBwR|Ty_Fb%WikEnR|BF~kDzLw;J3}%F^)gy&@iX$tG9Gvz
zQxr<_Lp{`0E?Q%ao&EZ)+<5lo>{sYn#mrKoUu(LB2Fb6=k?ch-Q0o^5pYKUXW_k?y
zS+++N-haTxQ{I~y#|I>~SXjsk<N|DU3Ia9Uw@7}Lls%8W_f=~+o7Ce1IWo6Neys@U
z(sEv1%|z(uja64tUfTxGS6M-1$hjG(m5WW`dvvE2Ni4k_mca>qZiO%B%;4<Vr@Zq}
zad?|n+hb+&NUIwL^-T353@$p+A{VNOLxO#ph0R-TTe3v}fa(*R><-#0@xyjMt!kX?
zT6V1|`Es^fy#4c3dkfa)qwtg(U2wd^6h0Zd=c{7uUbyK=CwN~YsBy16qIDB=lC8}9
zf?OSk1EzvXkh^XWT1AWkX0g&b<BXAo1S%Vw8jb7@@{)Q@Q`x?a4TV?Qq%*(r$lU#~
zl^41w_f`T}na2y0$SB;x;I$RbmCw6~aW~Q^Dp;U?x07EwEr8VRByP0uwj-vDEmbo;
zgJ<eXhr#LXiZw~o`yFUyuPNO@y&U{VwRn?tyn_3jQklJYhInLHN6`Z{5`M(wY)+uy
z69nU;OgzSOjy@$hUB_*KlU~1Hd#FQ%%t+Al<pn%o^sH|WhG!<^6sPXM92Y-Qnx^23
z#jbK&KF-5}L#J!7!rMMw<H!A>E}!pu6*WaJI5=T!9u~Q}Wt$A%U>(n~!Rc6=i;fW(
zY_0chDmeGp%=u4`I1Iku_Up1!7VLr|Y)Ujq*RaP*nzB}uGi7O)21dZIK?Q_I&)l*?
z0Tv4Ba|rW*F03KNDX-0ab-TmlZB<zT(OrA#CdpW~N4pDS5+niZ`on2!$hEllz8M78
zb*d!wN%er4BfW{_T|Ku^m6lmK-$G1&!jtzXMA4b}!8D*tYEB5^?7!!NZX~TvWI9$a
z1qGB85Dtrf%BWT*-Cyn6m2#V+K<i+@k1k_`9!l6aiUr$V_XI|}w<Ow_0KFOH+`~zE
zK|pAsw^ayb)>V|MDC3mdE^QFEH_&*Q0sWIDSp<`5&hc*)4EKt76+v-+U5QY;_6SR?
zTe0<?r9eM-1b>B;B)hdi*kj4sUU}p8lvoP1#cOFI%kj+M+uBO^YY5KuE+CoGR97jg
zRx(<fufc8>zglqa^&cC4OQPf`8Q79i!Z2$op+A=?xc*oWWW6h$vObZ;Udm1Hp-0j0
zhn5_-SJ*9jMabASrXBTaU)Fa+IsfQFb=WIqoPY4)@x8*6<M_?5E~0G^Iq1~kGBRY{
zt|y@z{d&C#GU!Aj4z^pa0@;>&$m~v1ALLbNk%+nbxy}!dAEJB6eL+0iEWr(QYnd8R
zXVt34$BlJQoc*X2D8dcuf+Mdwjl|Ezdf>pYurqq=U3T3qy4zOi5p}*xqE3e?k!yGI
z#<S<Ol3h)Za>ZUKO%F>a2{NQz+QLt73DIMzt#QNJ6dCP`dk9Dw=Bun^R~3$NmZjHL
zYDCHq?W;9q0-b%Xr)2g(N_Ov2FcdNYx6K62Nzv_ra5}LqX31bbCBt~O?J00d;M3sM
zuE@2lJ>6(>Zp`(yX!T1|oF$rlt6k?p2;aMk1)P#W)X%i6yn>|&T_VFZ7eIs*frY2e
zJ?Hu~m)wECogNKgkM0j1kG3g%F+ArhUbv4X+P`&Ec{qJ7gS8?v*qP{XDSoYah~%8t
zJG|-cBI6?4bPfYcnf9pmmN|@bSVBbe3zZg(ZS_Ieaa~pC922|2>Yit_c9F4pEj|Ki
zpIcP$Pa!)}Qn2a;m8gj<p`h77G`D&1dKWT%7<d~{vf91+9!4O}w?X2Ec2`+Dnn%oA
zZwzMFR+@2|f!OA%A|fpEN}hTeSmKg=rwOU;nCk1Vjme<Ycy?t8Z6QaZ2u8j~%-GSW
zr0rC6K8%M3e<DI(8&Y1(x{{~1-X+DMU~0JZY^J5IQbjfyma=IL`_UzJ?2OI2#&ck+
z2A#D0sS;MNY+3*_@k)9!9ZXydSh7sJj)C!bStOwcM!X9xp>?xweO;cGgP8^8q}4u|
zd>R;hP8lR)e%<xSu5+=a^I48*)gv~?_6@6YDcJnS1#Np?M2K5emalWC0A7~@%_~Z;
zeX?SiLV@PyL3=~=&dY!pdO(36L1!arZ)sXnPY_^dC?QnJG>KEKo*{KFVQJS(oh6e@
z^=MWtCHW&~J4f0e?Wl60Kxb%q%JnC>qPJuh)v_1a9xBDbIzi5oB~L!=)+?QRr=0Mp
z{-Ja^#zbe088h=>t>C7j&ywi)v89<cZNkzQ3N+;{ktt~UmF6o?PDdbQ;qe;^G-_a8
zW?}Y+dmIIGWB5rh!uc53=>hh8QEQgSi0n9R|2i@eLc1|J0*yB}{c+`w_~j={A<Lqz
zyNQ%P908#fq9LlN$(8uWZWcnKBz{l7;8@Q|a=t7&9+qNuk`P1Lw7a@bgmf2Bpno{;
z=~4o}A#T)G`f-wSfX+5oei4|7Td%L^`((?9UJm}sb|<YSh>%y+WS1&GnSmq?!`ikl
zd14xm9F{~e#j)bfI=qh*hdiZ%z6Yp~g1j^9D$iRM4&+iAABcZScQ~05w%S!w)T<RY
z<rMJa?vj!+e9D!_E3mU_*tBAEne4r}Nc9w1^a$3a-eV&f1Lcgc;n<*<%qP=?-FKuZ
zd$gvgwL^fLqE=I5U{@BKDca)ndck7Wz}DJ)mRHbC^O7$EO9aHtbziqDO}1Css-mP+
zQ?`2{uly@4UlOWT90p1~9Np(AL|C+iTZQ<gN%P_TY(11Cx`_`1tH;!mA){5|i&;X7
z=CckKd8f015%Kt4ALJyZ_-hST?(~ue#xb@sUm;1ZKe^>?ZN6@0mTn)!O*&N%<)xtt
zc0TMLrDn6OB-r`vl4VZWm1O5&H!Qj=j9`0ZupcdHy0Yh|QAinmurxL<DKSOhFDG^1
znaNTY%(ik@C8ZNV{+KDi#FaMvmn_dsN3V#H0E+M6$?cVIS=e&UvMVtJUBMvd&RyTC
zPl3Mkv}zkmz_R;#Zuc0ik8pE6GANI$)3u7cSe0Y$Sk|#SB_~ubNO%7rv?MyP#^sS$
zp9!Hs)f9;YSBto*f!Fg6dQ}->M84RRGKcujqUG~h%^c_M2MhPG<Qk)4NLA6Gcgk;o
z-H2sSGgpkCEQ*e+>ip>H6>N`s96#;R^Pc3k>4=eCz#s-ylI!4OZy%sCxGVyuDVVQV
z)C(s`%A!|5aFBC_po!}4TMk1NVqB%cI;>Y$g{V=M6_@U4@-1#g7w>^#v?}$aq}8rm
zdC{DZ*!!xyCFf|xRmCjeq+Mf>Tqa=s$y*tBRqM)mo0>0D`KF5xR~SB9ZFt~SRW85P
z>5oSIkm`>%Vwwus5s)sti&Ba+4E>_W<eAgkIDS3nhk_ui@7y1&{9k7uq_#aD)%6Z-
z)4GZJJo9(+M0c&j0L$B|4|4I-A#FP7#gLsk3EWPF*d<|LLd-a=VKe1)yU%nufLWFr
zL?g@;+0TB;pRC9hul!pG!#>#lF<n=~`Nf#9`sp9;zMqq-LljCGJr(yn?T8i+{>sN)
z&rrJGTIO(y4_%v~n>Lzwe{K?<CXKP{HMuVgO=)gDvgUAzCS0di^SH*_V7u!`#`&D-
zal)7|pm}fL;j>b-Y5+zTUbX0a2RM5~&dU5cK@b**+*?Yqv`f>qwV3A|@5|Sgn}3=v
zZh883=8ac-0Z?7I)u-(GO$Ghy{1^B5Qd@Eu&9-4<_9>#DcUE{qIwUqy_boS4ZN8tz
zQLpSJf4ia*2=A8NfnQw#|M&B*u4Hf#x&AcWEM%1~A3VgM(OqAY{W|;54v@r{^wxM`
zFz$@ui=~O4V|`PWw#y>Z%&yIOqNzcj$9zB0=;pU)@+iVbE%{oN{dN}@VCE8w(K4AS
z+`EZKJ+1DHWt2UOJT?#mtG)29kFdf%J(f$Jy!XyT>B}>g*!vG)#&q-8hCG5WroUwZ
zSyr$YJEWz@>J+&pY9?=<G6@W_cf^X~?SoTd`(oQ=rpFLHldu0WU%tvtRTR0R7&D3+
z|7d(z$<-rChpsD5bDe{;^JdFy&3#qXei>vH+O1)NJa_L;6D3myx<Iz|)*lAn(KW>h
zLfk8pLSxI0dse7N@=dEPZvshOa&?xnw9b026E%S}$I#_dnNLLgbOwR^p{vlNK2C>@
zE^yzncox5N*+J{QXK=cZx~Wd9^_n6Njnz@F(&I-FuVW>J0)LK1_AeW)HEW+S;<0pm
zd~8f*T<I{^a0vrl(PlI~^`g#~IaDF^ko&9lmnU&f{T^@V&BFmIGBU(y%Dj&HSk#q{
zJL+wvp2aT=7G!0m$UX#s<x5qh4Cmh8Cf?A8A4%CxN4y`3Uw^v@nUnpJjY4JV)m_b?
zJK2bSGGMEFMQ_#--GH~=9z^ITC4W{vd9`Zr`&ib>_Ve9a&_aE=?6ss^h^&LptCRkb
z$Qq*TX8$WB)N3t&T=ZS_dAAN7ZylS@fo<u=IR$R5-J;Hj3KhRjfDM=nI}ClOOh|=G
z&ojW)soeg9X7-fd!Aq~|KFIs*?<VSa3tKtp+xBA%dVMU89$jc_Xb?&m?Xb{K;X26p
z*<B?EpL6`+zJ)~8vXw+~=uBt*J?R2k*+sqWR-sRNq5;PW0!|XQe<axTN{KiutL*S?
zH#%y)2M(}q4N^n7g`fTz{m-c#he^3Fho$ysgo0!&=MlLV*pl-F&?E1@W;2>S<MeiF
zJ?H38>Jx3O%Wp4^v@t*%tW{Mngt%%G5pvr|#rf5l9$tyYIwR;Om0K!7xUu6f4I7?Z
zD%ER)6D~B6Ct=-t10`<pz62Pw%9Hpf+cot+<jLCXk0;_SkuqF^at;GBT34@N?A|rq
zdatZ?Z9q&>@G2%j=y}b8ONy>~#&gkvD|7#%2Z8r_f@!ehCNeN4=~(1b#y@j?;`@1F
zx6$aro|<9fY1oh@{8unX?Ue<jRrlpwPe6Z7nk1YPWFSjh2u=ITW9zv7BO9}4k&jqT
z6kEMXc<&lo3#$zgxJuPuxVCvTAp0VXYp$iw`jGC&_PK4XOCU?;{>zf19o5<Rd-FSq
zI_VmIpTflD;}c$ppO7WrJg(esf}IPo^n5N_G^eHchHK!0gT(#*JWqY#U)y8E5Ns+_
zm*BQm*Qyv2RlZ^(3K^cSxmKIy>^L>5fW<nfRp)K<1QrvCUwMyXt;HVHY!~r`7MX>W
zGQ{j9>V#{kZ5SqS={yvKmA^wMyK_y`XkIrRvhI+d4(oWk&LCgjy;|-U(qPt23Q(cv
z0`7#2(4Kdg@hS8Z+|wgVItekfvuO@J5nG2GhmVnQjXA?HzY(W?mUED(MjFUn4Nmt+
z{@oagy*jpkVp?jM;nG3Z-ua;2`vz>n-q>*s)9$%N#=n;P-OY@Ce~O+s{GZT*rq^+l
z`Sz<75TFBh_Xq`Lr1@NiJ;a-zDdJ9WBt~*mH%<g%x9<I+?gRGpb40#_^{-8B?4mpL
zf`V;QLXDW*eq4f!0_qODAJmcJ<I3+VnYca&394L35U2GM8LbcT1FeUOPW!{U19Ck;
zoP@W(W`w4-W<<P3@)LDl3Wc6lKxaqJ0%P9&j?4UKPawCyCw5083{23<ou_ophy+7l
z&ZH(zDj?DSVYQel>6TdNMAGs78ZmS-)t)3_r=Myy+Z<hXdfI|AgZz;0Y-;e3zN+dF
zL48qko7m@`5i+MvM4;CvqdfZhSj!|OW`m4AUP*D~<)gG5sqOl%CI~#2Dzfqqi|+J?
zJ`@Xp!;PA4o{RoPFo4C7Cq4K^=~AF_5qRso?LJF+#0JJ#3y-dKMU_hKo<ACI8@qK}
z&%$#iq@B+yAn0zB`xR}JxiGx$dx|!d0qzW;pjJvW#S%I}yyjOF<jTBNecL=d7o@N+
z<o}M3FMn~9>n{#0pj|KYMqI~rbN1F7TOR|H(eERZ+)Mk-TaKhISe@?Ld%5k02|2UT
zqE7=(3iXk_CEnc?zTvp?IC|v9Mk}9MijnGY1#@ai0!tWVuB;xL`&+5!b*j6$T>2jc
zkLp}2<6&5pIGkb<C=@cC8^GPB{ykAH>Nn74fH!;~d?l_X7Fk&8F>Pr0#i@=6kJ{0l
zE;ipgC4>^^^M%yvJi?1@`yO48nTTo}UDFcyl5~f>wz{Ub_&9-vRFg2v_55!Wch;|t
z*(4YDZH)=ZKO=-sBQ15SpNxzz_^^uKTxt?rbP7^E5X@+-$Uv84^xxdTsRHEJhfhfp
zyfC$&P}yD%k13yL*C`nI>(MwwKQv;N`S}{FZpun@sM1U1gP%QkzdZX6YFIAM40eof
z_MjDh;w1Gw-KPjW-D;g%Sj?@J(sJ_HH;tfmxB~JHugr0uJ&o66NDU}?O<E_#x*g-r
zf@q?FBL+B7kNYAAGp1m&vqr=PU+dtK{mRr&43oewPkB9@eTi0)JvHBj4xu<=3vYEK
z<}v4;MK)Z$OYVF!@%SHtu>{_vT$%%q9X{D|{WQ<sYZ|wge*|(Zp%MaI%z;4*i*tN6
zjwl7TTBlp7li}b8(uAP|2sS{35A2Bcc_87^Ljwc593~#(j88Pu0{Xgb8*Rc^*@=_R
ze|aX?F?)qPpLdjd*3@lF)(H9?ON7T-8pfaS0(lJ`!tn*VrZuph9ML?wKy+*A!X+|c
z5r!Z6^nzn<lPpmpo<g8x`TsM?=i9;1$a<{|U|xtRcwQs>rCfVHD}c7$+Ps5JoMshN
z@NPY*vb4_k0_WB-v8j*PmMuD}t<d_-Pd>8}Js+$19`SS;_JDP}O3Yoa7xCVP*49%H
zZj{Y{A`j&4pM_HUcIqODUi#`I-y8aKkG#vP&*lxiXgaugMa)=<tssWk3lNx#e&W|j
zAem3Ix4qMaU|Ry}$=+B|k$w-oN{2BfS8;Z>bN-y&EB|G7MJ<M*YX#JJaW)pKBHXox
zmgeoN@4i}BWRf>iy0M_b%jp>$^5Q~%*DzS0l+)|^XQbx;D8S5~2muA;u;FCg>eoK(
zQ(k}tbT74o!Z|R?n_shui)?i<cy{sSzXU{XyLa9b!=$#JsLgso`lxVI9b*bAGAi~?
zu-I{}hIqA8{|sY07&P(*9-{Pa-o=fy+5ZN#A4{(UIQ&xI#9K<yUzf;1uFW)T=fhGS
z%LSwDkN&_40FU$uMc>de26lpVzR=nVZhtvkY0Lr8KLvAEp6RRHch8k=11afp6F)y^
zyWG;1g~{0p_ob961wrT27j=#v*nlF!lP^3hKI&v}XViUguDBCO+jh<jeD6NSTRvNf
zQ<>jsh*Rc&?8!CE7qCYSfn<PotJ&*99>6@>>_ha=n?{hX1`iJl|5Cj9_S873Zwfg0
z5d74ZSzNRTS+zpz%P>}eGjAtQ;k<^U7W818oe07pM0v@bU~CELz~|+=6D+oqY(DN*
zaE2Tz$aju9BHwc80}neSo&n3mKi)=Z$?)w`r`}b#BIV7&Vz}gC2^`N6k#DFt+aOsL
zb7+v&m(z+#o}kKvuVKVboMn=6Yp}TKjbQagAgG?N2JXjwW{j&~=!)}Qi%nqSw{T(}
z`faNj`AInu_M~uo<l)P!X`KnX5JvpWz9+Wy0&K?kWBb`<hPd~(ZXD_gJuc4Fq2tt4
zl$Q9wO$>qtFvfLCG*r0osq1mfUyrRhxdrae>8-2Rl!fy_fivZNgQS?tVSh^#_+D0o
z{qea+zn%1g?e!4jNk)Eqt;))$RAMn>AJYi8R?ah=pD_ey-w4u>&6`7rg+UnJ{A%~e
z!J6*qrl0L2NKK`rs6^PC_a?N=0%CWox46P58Y-?cXw(1By1aRZteIlO8du@K^!OGu
z=F#u-fM(&Ri)NEY;TZ5JM*QkT<(Q*Ti1dUawG7e_GbZ3)2N6t-ND!ifsWRhbE<)(1
z<iFL@7~|gmo(XL7vUq3e8SV5MD2vx@y1-EG5qVkQxVP@0(x$+Xymmz}lAFDMO&<BQ
z;=-IbBkITReuoSdZ*JN!;LYq3j_o>GSH`<srVe?d3#sMR$A7tJ&VyyocRqI?;F4s-
zKRtQsw?9#7sMuqXCoFu;^$SPa_Gi><@>0e1AM1_ELZfXuE%Z@afQSLO@H4~Qgu2zg
zZZhJx4$k|_E^IU4*LbNzveZ;_=#f(eM0J5#CjKYX*syj_fxXyW>d)G|p?(o@fFS0<
zGBWY2{T4<6`6zS>Sh(w$lQFbe`~U-oVY^?ry-M%t!|hbg%(tHq2504EcM$pA2>|(?
z?kR6~MemN^{3S)$2i(q%uopeMromg^7coc^Y^d-Jx1YO|882}>>Y%Ldsa8D|{=}C)
zstuAkyk5U%z&c|befg|>)TxlGdW{xy={zciY2o&h40ta5WXR>$7j>>Nf_e7UB}~HJ
zGRaRhX#R==5MW@DNAMxMS#DZMC|M`+LIsj>aX9;`$OgOI*2E>rBgA*7yb<4#cISzY
zZiy3=k0-qPjqX}KRjiTC`?V1I{OxP`0R|tHIX}_X^2Fb1n0#-aG~FCj>}LWG6nPUf
z|L5)RHHF=FRllSGM0RU~gP#__rw?3t5YsPyoJpCVI^-X2?`oRm#rs+O*I{AVcX$NT
zbs8kk8>#<B7GyElG4DUy`q7d11O<!Py$3n;<HO>81;Q$AbADDN<toEx$BgGz%@H>S
zA}rlW8=8G3VZ9)%r%_@|cr)L;eQx%L5s!`@1KWP=79bNcdUZq2y{tV+@vpNI&@dp<
zV%o(S!N{L+oWJt-{ZFMqvd$&!5MJMoukzbv%_Xk(j)scAcmF6WDQ%LYR<tJtG2mY^
z;$Og{5FX=9{26TWQNO7xNHr4#g@KI9t$7yr*$C|t%zz(ejB^=T3}HGC_9!#JZ{KD_
z^=(b(QvmZk^3PcQm8S<&9JiXhedk!-9^>L6(1tv*|B=RTm;%DuvipHS4PG7WP>$Ue
zdRTQoY&n$Uz-c}&x_>B69m$kN8coAYGr~YWovbzd!vgE?kfq%qnZlvQZAfUSWRstw
zntPEp$^dwq$|~3<n=JAg`g)Q|Rb@oGBk!5;Wel8p6Z4IlDuPH2<p}%c+B+T>-QHy}
z^V8Vm3p^k=$VWCLbcu<dWD{<G_w{VeND7IFzjoRi0jkSCb{On{H&mIWgyv|e)7~9}
z<>+74&j`gth<7$fN;F<Sw*?&g(C%A7W8}y3%SU;Iu#7HpS)~{qimPBAtiNb6-y#ui
zA4wbh=Z0omGIdgi0gt%b6wBf&$0jIDHG?!r#;Rqi3S`Qg=<89><ed{ekzB_S;2myX
z`*azVHLBzq;Z8Up952;OPQTlOq+xHo+u*=}zxTfG*~ksxmx~24i=hl-=9!h6XxZ5Y
znxq!hbnD4v$}hRomy-VRywBrPq#!4C2&7HLalt_Q)`P0+Khvv=5>+ui<1eR*eesW8
zd2?(MndF^>jyjjG-M`D8G7M@unBFhTO!&wF)|~anTjW;a%yobX|I~WY1Co>e5WIml
zRJ^Axf|}kIApNmEx}?NO(0g$0aRqTS{jiGoA8H4Mm*tQ*nNCz^%?Y>PWx#*^<t@NM
z(7-~3z@Vls<xX8bL2xjA4fPAC<2Y<a1V}S}+G4?6sKat#3o}0-o4g6>+b_Q>2vCRi
z4)K!$SU$-$Sh$W0NvG4ou?0C27cgdxgA;atA`B0-mBV)e*?H&zfvn^i``a6U4MmO%
z0PCtfJUMz}KtA`Z(qM!=95-Gj^z?K}h~-_-;-8goIThKdj~crf$H;FH#y#^p&(LuD
zbNZ?g_Q}rm0K29exi=M54!6L#*Xy?jzgf&I9rSU99W}8x-aXo9r3O3Q`j3ZNYDdZ)
zIt5BjXSjPSMrH@I$ia`5mpH9}U4dKOeMmS^N3+Al5vlFg=g?OUiY56k&Ech<>9!px
z6K}`RCr8-45du^*02-~j%fW*L{Vs^@E=RzKCPBpr)LYQ>di==9X^h)WrX1&2E8Ap<
ztN2Wf<qIsn(>SF>E#v-P6ssh+!C(+ATExV^k`(u;!4dMzxf~meI{lyNBWVQ#IjQ~q
zn?=ZxBE~ZWP|Cx??Q5p8ye_!H4<;oyz!b-*qE0N{x=b~@H8u2~pMA&|z)`6)?>&dY
zI#-+r2h*)@KX#EBZquf0v5Ns8!8C^`xbjP{R4?=T$N$Aek=4E<(W*lKJ@Agn6THIF
zew}kxAoWZ@J26<hvJgdb)Di+0Js}<~`u8L=zXdb>DGB&>4fw4$5OX&Hsu<lgsc<f5
zQK><)-u69Jf)4w6`Os1hx4$Zqj(IMs@+>>te&yff%NV~2*5uS9k*6gjfmZE?iu3lS
zf#vu2P<cXMSl@=*H>+snQq|#N=c4bDn=fTQ&e&Ph)EV6PfQ|%8=RFme@Cb(?)Bbd_
z4yW1QKM9H=eB*mD>=+w#2p`3WUpq>SQ(ZjKP_fVu<n_zf{$@P$k9fxKCe=;Xtz_81
z9#e5#AD)0g{(grr1fAchOtwEod-2Ttq@Rpu^Rsaj3=;fFcr~r)(ph=c$7?|JtwST+
zC-_)^SpN<A4Dgq@1HLEBM2jx+_vQ3o{jZcmhQ4$$obv9_3EvsOrN0W6{h$u+@z_Ae
zTTfHPI{i5G;lY+6N4p&9Ost&LR10O)!q4kJIAA-$$t-)U`H%M7jTQ|}RL**4Z^VWI
zfW`$IPw%*}*}oB@GsjijjYxe{ozDfPnU;ZG*!j9zoU<zIjo8^(Z&@7=Zm4d?ub;uU
zjDC<!$_Z#5&0O=^+Rf)|9xmLmesSaj?{*<!8|JmodD~0leC4XP)-hE{g&P-c9Q<AI
zym0jAgQS9#oWlQkCoT2=c_%BC9AakoP)6>5zAhytB_s8J-rbb>->syiZYuomvQkpg
W((=hSpz24Z737Wz390L89Q}WTB)YBu

delta 9607
zcmZvA2{@GP+kTcPq0%CI<+W6n>^m)%LMGX_k%W{nmKZaZ$5MGE6j>r<$u5myCd^nW
zTeiqPj8XQ9!I+t>{U6o)e&6r+{T&X+Gtc3Ep6fdA^E%J-p5;oGpOq|6HTNnhC@7zQ
zrpN#L0_`v)!K$X3q{Si*I9*Hdt}!P~IIqYh?allK&xRd3X~4>&#CPC<t*PWwX%T&n
zm<Gsz1tghrAM%m=)?UoJaol%{z!ABj=<SzZrs<_N$n7OM)p-Z}Q{mg-v|T;wK&^mo
z5KAhspNK^jdSROz(L-WuIr0rZN-Gv9&h?wCvDBTNim>$!S_M$=2p<X{toBswY{B#I
z#|;5t>#YFFslpjJQ{mj+FcErYqmfe4U<-Nu97ljszrPzQF>L?{^&4y9TfB;TZ>c!n
zY|ba92m5PC_A#7S3OMb>;5M<l7t+akPakeg52eMSyF)@~<IK)h#KhKQhx|_ECsSZs
z=qSN>AR7SO8?N_t5YsFCayBa_ob@&_4BF68!{#c#b6$p63*nT|J4(XL8ybqz@vb9b
z#(VW1iCiD1>@UxBkiQw9rLjq97^l&9Bq%e~hWqOY%Npy)lxaYf@YR&RZiy@|+AHSN
za$1&E0Gm6D2`5svZ1UxpAsMNwKfm?tFcKzqfDPQl7DI1m@!$j@j7GDmgX(SK^y=$3
z83}s)Eh`%%@P_TB?2ewF^K^PeRqW5|IP8W)gF1v-q02ax^S%d_Hxp7?CYQTL?@t0j
zsI3+U!f&HqIc>wTIj_uc2KHWkyo-+ZgLs8}&K>?Cr?PR0g?IfkI#mi4>rBjXptrwQ
zk;?RL(#wHY<EuLzHxq8wGacGHVLLuElJ~+qX6gZTfBsP07swXU&x<=#!Kttw$<Q+F
z4~tZobBW48=f}-df+bIBmGj3@peCcll<;vQ261Ry?QryGku}S>GZw1T*X<un7^mXA
z*UG!*Y5YGZpBaz2DIDWo5pE{|#V+H^Y#jtrpn2`(!h3b}A!_KCehre~hS-6{@gZXQ
z*OU>om*O#=o{ElXVo=z%HrzD-29amg>ctGthGB3YPg!=lKr-Y}#LB$>onVvxudIY#
zlI*AQt!nNzDJ<g7W=V)9pePHAXpwy5G;?5!$35rM+wyJi(%z4Tf^pM;bftva56&Ui
zdY(+N0GYll+4MxL@M)C-oSXx%ol0<sug0Tt?{Tc~0k2CF1$2TrUZ}&oIl%f@v}6GT
z)-htNu#!I_q?~fOu)+HExs=N`Zgcms?j5HzaZj+*cRC(T)XGNx?F3#wHNT}vd^`Ov
z<r*ZzkypB+kNDi~>;5(X<E3%CcE)m4=H{^-Se2}iXKloTdYHx=#^8O8&4<(1k>`Q!
z%IQP&#w~i}v~;~A!H+Z6|8i(G{i$aS*&(Cc=-VnYaC1wgY67Mch+7<tE;pJC>AP?)
zwC6G+r$e^~@eFV(BuH-fc71UcBDR<y)3@F;^6=Lj9pN6@*Q$W~OG;2}=bjMd-e}+K
z^}F|q7ay;G&()*p(#pMt<s1{j##Qwk<zXDeRLY1x4>fGCZouI~@#pBXGauG^)?T;E
zm0!)ae}T;`hFf6)n2}gUS`tNfWAx*zTU-6Jax*M6lTf%~oUQW`XEIaHqd|H9aiaR>
z)#$jLnu+HPoftqYV=e&(5G$<84s{-3`U^cIPGEyY>bJXumF)tFJ??Ur-;+D94yDcO
z&J;YXQS{H>Kr<9<u65#c;4@`0{ByLlx<R*5ibC|Y%+`QTHw%@y)|U&4RiT}*<n;*M
z&%*b@4{YbxcrVar6lMv>FmPq&H9TwQnemJKXhT2#YXB4CCk4AdsHQn{3E{R>MEWq8
zx>T}oe#6&8GpketF=UxN>l<BSEvV4zu&{PBIK$3*@y*$Cr2T-4X-h+WSQ5U<_G39q
zBYw-nulsaE6{WkTozb*B|FRjS`;sO*m@n!%kbY5XOk?7MCq}oAyDUb_qnGy^ocon!
zmQuz-8DJmSU6(Xro!5SMZL`cex#~^9fR`~9Gis&UYusI9jgm)nzjt3vQ?(M1EKAVx
zh_ULdNmsS*4GR*wTiOlDXxHfQeyOA3p&)1tnU%>{`_ygcPPjc=)~&awnhC!X3%-rv
z)Z;zSFW9sffnj*m#gJ;8-MhxLN||P)n#XjYeK7x1)lEbV;y9%agGjmARwh8ewAp+@
zyI*@W&O`^&Z@s;A&(?s)*o0)4*NG}Ce)RZrpuK~WUT1Xh?4yrLWz}YzI!TB@_Vd9W
zE=pl24?75kvDSTKVreQ<0@bwi8fU&e+AidxW53$1SdLhAHyFZ6QMEG+*p@QgLaD%%
ztcqqDCFslOoX(h5teWgOqz5{yTG+nOrJEL#v{K#=wWzjl*NSDTP+_Y@^%nu@3XX<h
zh;V|v#hg@YX-@{O0FmOEE{pWAV_0>Es9JqF%2Q0o(Z?rwWb7;HH{k5SGhpRSiGvO^
z%*o7qN}++FFJeha8-HztP^U71=KB7lvZ@LW;s!~zai)4jJ+~Te9xS3-<#qTmqPlU5
z@~GMFxkduU)}}0*ipwB9*)%7h7ZE6rX9sDMB)c_uXVt`Oq2Y7P?pu#uB~<YXiZE|O
z(v>lZ1al-t=k7Pf;CYf=%1CZ7YPS7O4*}CXl>9wswBmWErDqRbY7r>@5mhQlmH84j
zbUR@{bYLp8x!l&zkXptNN6w-Uv2x>LSCsFxA0^OkDpN>ym{=|3?5FK{5;6!_>yBb%
z?wA&i=kDLAs}YEE923Kh$fG=bqcM*V*7I-E*SZy3&zoiAF-1?fokgp&+qEbg)!#A=
z?#hS}Fw-kWf_7CpKvLDi?ae!X*#~_*T^VB(Z}n#VO6GvzsT6;zOorFgNXo!P@2{oi
z0_QgGdMbIU_zn*zAyQ_qITZ2SaIQoQUXD+!Dt@;G?zgi?pJA7T^=kJ5&zt-%{=5Us
z?%Z*m!tTntPNE;K&2KOTkI|Rcp{k+_#DhKPk~X8|$%Eqr0Q1HjLPCCADpWh+UOr~u
zraY8|$Jizik`Tw&U3kU>otph)zpf@f*HW`<Ni27ADO(`fB?R~{$@Hd9t5y{xnl#$K
zt4<!{L}P+I&|zBB-3yDdYIz7~Cl5#`CR_tyO{oQ+Ro!U*B-q1=nuy1A9?TU%<=KDK
za)Spr0i6c(UANQUq~kF&Tp;AZy&1h)YCQw<nFi(eK0O=|vezMCWWE{KU{GwfM$GO@
z;(~tEvcb^&ouyl(HUe#JEQ4eh`7@C~n;LnXU_+o>q2+k+pGof}hvj4TjVGi*nBAwH
zp$9)pAU$3rJj<1}5}M1aU{0AQ1$&I7NKk-SeH&S$;L0pl391ZD!KBY<Wac=}b?4@~
z&g3<_wPwHIg}NsVug3(&$J(`Q;)e*BH$A4CF#pg{y*Gbn!E{Q3O;EF!wpX_6wFO`4
zJoT)<p<Gd8ZNgx$nKNS)C~5jq#~vY#a5idrS~!nfYxje5new8twK~t&tMklJ9xniE
z%urQf2Sx)m``&SU5uK&-7qx68Ze3{9znuf6dA9s3b?JL%K)T;jnNj(<^&^w5TQ!Km
zJC%|^KJLnhb32Oh9JWwpkG%p&ZMZG=ihZg)rI}mTzr<MS<~#RZvscTcDE#_t9wO@+
z!9OX^OD_olRX@i~mH7~btaDY1e-3C3p)vw%tCG>>35Hda!?8t)9#omZz`~Imuo8)`
zf|ea5wpseg_e=?!>XGTY#e4-(i;@^l_-tu8XX}U-r758AIGs1p?;?y(LUnyEGr;0?
z(#x}kJ!V;RnT4}`{q8F%L#CtkkININk2(;X2oqm8(UrwWbJ!YwQwlN`2cTkJ`p*t-
zzakF0;*4aVw*&5JUj5L7?9Sw@CYHwpTQIxd!$}<k8f>wqx}-l0=bnBq_4Z=PtYOuo
zO5K#rKCRdQ1&YbN60_G5_Z$Q=Y7va_$wB+Fjw@0*sKLlsPG)z8l=D(1d?7D*wr(Xq
zi)2^OWR?ih>!&P#ff^rxb`cqOsHygaOH!<oa2D_}i>0!=ivCOAJw%o=RqRbks$~_y
z9-i>2F*~dkqz2okQeD_9#z<n~-T{uHNbXG{0aG-;<f*+q)2`A=MD|fY0z&3*J9K?9
zlOmoynCrsqPL|4B$~;G(Aox$VgAY?U*ZCSoa3eY_Y2E?=UA#YH`;~6j#qAt~+!Jd@
zF*QiMle1P2t7Eo}es(kc%i4?0iv;t>?g910Ool~iXL8NNNA7aG7~QK7<ZKt?C@hcg
zQmaGx`TffBF061D!cfT@A_zY_Mw%O$O8)+=tmhNSuJcuIl?R^_S?2v>rat1h&Sz2M
zwfJ>EDj;+B1~|+#JrmF#L^~ecNtlI=#1ZeUyVa;J*60_8N(Sp8Y2S8kE*JcT_s>DR
z);k>5iQCK=2Z8(^T@3*wzC4Fuzj^g<XWmbOV4vM9%yYN1dMB>=#h2KA-l<%2M+rH5
z@KEe2C5Z)9LE!Ai+&RCUtAFn$j{Be{-_y5%+`nD~h0LGsNQ5G<bg>;mOyc&EbtRI|
zE?8I#yxn|$xAA2{L)1*-l?zY4Svo8BT3YOvn<uX%%dfvosPp7GjN+Z!Qtx*9(ZxPw
zWmc<jEbl;SxW9X{6~ASoG4;Kye)3R9Ysl8({!F-l?}#|OxHj#ie@D#pJO9S9pcfPe
zYd|LRA@^9{sx2hd?)>_?pS^tET}6S+Q;#!_UTed&A2u0Wovz||H~*t9#tp8~l(O>1
zYqh{xuJ12ZZ=~9b%6jn(OP!n2?Ng7hVmfZ6-!_kZcEqcpGa%{Beiw+GwxhDs`{tE*
z9f#01W<S>a?n|axq@8?cvU;e+X>j=)w6>}RU=u<h<A#Km_*B@tTDy{4vzKnrS9i?h
zLX`9aQ)O0N2xD%PY4oSffrJSA$$cvgP0Bm_8b^N~62s^pruAemE9gzS?!YHq!`?%-
z9%TQ}dzAfy?4&!ooSoYV{{A#FpRqYpd1hb5Dum|_v|{g`Cmiuk5JLz{SUY=+yVWKz
z9bI&?{+W@QX~qz1EJRo8^I1{1rhaju*$++?IIc@#bX{|DUpVy?Ma?-ZZ}DrG?FY}^
zIi{sgDOU`<mWx-_Qj)Y18E)zjWmC~?P)B=qp)~0E-UEC>*<gw*7mP3}hYa2j{MJ#K
z{zL7Sl#W0e&sWUN$vbxa75r&J=B=(ka7tH_`!J&aYgzC5leg9@(g(<z?NW*{6Aj;r
ztoCW&4XsT1D&U+l@KK&$PTM!J{}k~2D7|(gEACFE;FSl~+BSRWJXz3mi2k;<%`trb
zEdoA&GNltKRG?%EmwCZlq)VpK+-o#aSTCu;30^YLNPoG9-o?<zaeP*PD;>EG%vJ0e
z<)gPjB2$l!VhinVSUN;9D~Hg-s?x#fp_nYcXKEdbt}JaQ;tTy`tlz(?B%i(PaEK`|
zd-2uhZ2Ky$1$udm06xDA7Mw28h*gL0?N)abl@^MnFVdg)2Z{Ol%N&al$~u4Gm;|o2
z__Sd)ajea+)s9~}iFZ(Yfgc+Q&=M~Z8OkOFS)rO0%C*!xMCsts4F>9y0ZB%!^;T1y
z=A08)^{u7`eclqxrK@IqO#N4%oj)9|9%3HYjKG`az_b@73=m$+UyqBAa=5`T@Iec?
zRj}jZBUQPG$D9ZeT4GMpjgQEf?h=TF$f{C@faYNILP$Sk{9uzw2@oBRy!wIH2g(lZ
z4dfA_)!RX8+A^jY=wXw`N0<}T8g`F91IfLHgnMU2aIe9hO}QOeDH+HEtH0}h`ZuS%
zSyUeEIzxFAi%RYI)iPgPF-c|06I{$Z)%>QZU+#3Uv*-;|`vNqq;jQSQ7_!7e?4oDN
zNeT}Mc*9~=#q%l4eWI+;ajFZ%>kH|)Re#F5k$i9I<Smb5do|LE0>}KKE&)f{gZtUy
z591B(b^-H&&*K-!HdV0aI<U}S9SI3|N`|YQNFZwML~JUVR4wBuc2bPs>>I7+Ds#Zx
zr=KlN5TB2F(rh9)Gl=y4a1=IlB9_^O0)7Ks>1$Qgr8;_Benjx`x$8zwNe}!I&v9sS
zAKf6nGmb^Qn)*WzMf$R05L&zUf3XARYW~^0+9pNbAys<h{NMZ+qy?E*f_?e#2LpoT
zrB3t>-V6#p|1(ynD=~tZpNXTq(Encz1Ad#fl3p3?#;<*3WNK9Qn(xUOHP|muf^7|M
z?Vw*^iN0%NC*RONlQ?vyDz1>go_2kf5Z=c7;DGc!;=eSxaq|sb+wj&b)3HxEx0LT<
zv6Cf1xZ=X_{T}tBVL5>V6yF+fG)kmI%x_pmk_pP@3D}DxU%qG1pT`1tJ8!WLn2XM`
z0=QkM;?J!58y99dvRIF9ln&eq9aKBOXgA@45x>ZI^7gUHKVA9Vx&xtO;A?(eHdY~p
zi++!<rP|GPB~hrviKTZS-!^gMUGQGFZI`?t!UMOv{_T{bqb*;ckrtO{t%O~XufMUi
zph&o(-z)`=^Q3EehiNfmQ<v-6U%bnDcJ;VOk^FVZ@z2{rbPLYIzoFxM`GbS{+FBK4
zIOJc3r^K)Eop4!5OZ6f={cOP!SLc@TW@k0;Xx?XMBYPowTTs;RnO5AcKC;XjaR?Sw
zSvaCbIWu%Wt@l}mtAIpc^g56#ZHMa5jpx!l5aJ`u7I6uPs!^c45KcmTT$(fuzKW3I
zf||7BfAQNrG)6>GuWixQj$zI-Hx(>rcRZ5K(9%@?n@$dfO%i%O4t3#V=1?()4bQH=
zLtade(kG5%GUtd3PIBKu+q%%EQ(YGf?c%>nG$+=_ymaVkveW>C1MIYoI39f0rMJyS
z4S~Cw<sKvPWQrhTaKE!{Rd{5RfigDU)q+Y{J&Wh;0i7nATNLb=S!4pN;qQwsfv@?+
zN~EB0WF9T?ZkCh5)qyWR$$KMrrQlznLS@vvGO|>!x8n0D=BHj6H&>R@x$xb~m&G1F
z0*e1(sQq_2`KFA9CI`#eyhXmcj(^6PXuI}4=5U9y2#qh4e*ue7HEmrBme!6`SrF2O
z)#MC0kmm4T#0A~8%qs^JHK-49maLrbw23SAgnCpzD@g{I(5##7|BZ{OlR7ptD2Ivt
zTYT%skH_klPr1A;Edy{5`jA&;^+f*WljcX#d3v0nD?83}@?hJ4#ANUKi?Fsx4}NpS
zaULeCaS*G8YQp_Co;*>c>vH0Pk!H#s|6}+@xDjyxVYAltM>uaRDj5g{43n6Nk4yIw
zTzaoLu{2!y(QPn)HZiuyh;-tQh}>)fGHuwpwMJDxWnt?HS&4~i9bIVUDUVv~fz?!L
z1%g-G;<Dsh5#0BVGiyUng$%c`kcFqeScL=pjG`Q8MNmzA|L}?*i5U|%k>61%)`WpZ
zBYuY`NUo#VQnOWuB<ev{r0LQ>jGI!ULA_t|QG4B|LB0O&lJR@X<Mv$C40`kUp&-nd
z?(q}*Z+<=3`aaQ9)@%3&uch-%&5Jn=EJ@S!54xfT(86CjQeXpVa#+kT_o{mox3;U;
zPnf6FT9XjkBEZT<+=*snT;;C5$`DfgNnEzkxc|_|7b$9#f_%^wUL(CRHi0TNPleS<
zPe!xEDd?=3?@;(Uo)!p+NE0{O@2(USVg5HZW45gqoSNP4bdOJU39eOlqu=Ix20pZr
zl-Uv(>bU4(Xa}^75W-GnBLHY&Nw%Nier@D3G~LO7m<oRQSJ_L~?Fxvu<vrh)mTn4^
zjXJW<JW-r4$ZpCFTNjWCO3iMny~apnty{pSj_>2uiwY2$|6{;iXs@Ys)Po?4R{U+3
zpWai&ZmdH3K4NFqJSqf&x8sC|k(rl>L8B(zqreRQn(sMKqwZz=yLPNiMP)9dsvWBM
zwiR6am^R!O+2C+~J0xek%|B06!)N)ol`R)|EPs`qStCVKDnMu0H3pEo-cJ`F&FG|(
zQwqidnPrUJ`mY(Rxa?*MrK^Mng{O|T2>r?pibv%R<DHxVhsU6_+77E3Uy@taoyMZ)
zqVHQ$&yo&@^hh`$E;M{UT}v9&a=yah);FSp9U}7J{`kM7?yDKR&#PX>vbOw&yUz6K
zT{}CIqQ>8dB>C?r-L(9iyiv;P<%myU<=Q^hF+w1%B4;J6wr9>5lj9=!Nz$YBTQg@i
zfPa=2<2}7K_fl)}+1^_s<>PP$!a?^yNyiVl>rnB36#YzXT_m$_{Mdi;7gIOQ%RL@^
zRv`SdbQ(94<2sC!N9f^>y>Q|yUEZsW_nue7AQ<*fOL&o+uY#gE|Hpobt+zM<HWA5b
zom5kuZwp<ZrY?uAbocqyjns))ysal5=0Ii&hB+LOVDz4wg`uMPU}_EJG5_+n&zIvB
zoIF0m%)5ODT0S?|{ZT(yY-i`N>j+7<i-vl6ja~js75MJ_(YAJnqaXjHY~ABpY{PpD
z2*51gOT?PXEN9^N$qW8@pC{hZI4f|N{UB6Z9E!%vPQtTJYKGm|xXQwE_+Xu5%jQ$i
z2)>w0YR@wHHg8-K>sUGedQ(dEsvUs)=28^~jgK~m-WXKPV(>kHcAs0wt$P6za)UWc
zKp|J#Go}a8Q~Y-usBGEcp7SqUSmV*ap5o^B-QtZ)ds>4krThyDKb0Af^>w8beYp2p
zP+>0*W}i`SL<=-3ti`aVa)cbTV^bg0V{NTvX`K>aQ@PjZpyX4=;Zs&tG=98PNGyg;
z*p=m_$ohAoQ=ID9&V6tHWN8i`FMvHB#ugvLcV$i~<J2aCHNN-fc%Va*?P!ioPPfNL
zHuXo5<qAjF^$Rj1T(!boN2Tta$y#x<JLuyT#u|T(D2=-1GwB2w-?z;oOyN{7D4c{g
z_+sbq`J#;qvYz__UfNjw>~2)xq`!q%sMWkw>qr{AAOKjF66!BcW~TOok56RXlSKR4
zyu-e`9!o6#mEOe2X`~e02v_^Qr!Dic`qhl1;Zu$1b9;>G?=%0RJghw^U(7A++RpJg
z9Dl6huoT%=uTf!<NI6ryqC2JAh~D0>j^$A=kjmy;E^4vH%O4@X(Q1tODq=N35KW69
z0}}K-@j2|RDot8+KJ`gFg(W`cPud0hice@%c&;3P;UeA06dBIludd4LV;Yb}7ulaj
z(*zeDJJ8Qi#|NQLj>@rMSO4VPz_(U3TEm`r^v9Dz8(nyJu`gE|(UabgCg*#Z8^xQa
zg0AX>yGnbal=<;i7Qjb#;co_4{+JAtC7upxM0-acf6a#9n@3v1__HjrKHOP@B?T;$
zf&tlDPlAhl!o<h;C_Ys*Q1Vx^YcRT8>VGtAZ5{c*Gtk^)V^6Dakoc)q+rpieFdyiR
zgX9y&KGR3>Ki+w<#&@v^U!EufeqpTT!wr`OtFrbJ=NtNxP|9+mAZ30}Ye<3Ap%imH
zAM3ApLpSp9GA2bg3L&m=wGlnvs6YWRJ_G7EuQ2Ao>h2D_RjXaj<PFnQU}^ot5?`FK
zyutHxa(DO9E`bKKgUE9;+u$Re@c+VFmXOd0C_E=$MuPYbRV}~x#4daeB#ydVh0*&H
z_<TN}Clc!82cTg#8sV<IA-`M+JV)p@Du8kUqM+<wT@mbg=r3#4MknGCi4Z$fCh9NA
zCLv4K_+wD(&JSALZOzzL7C-c#H-^dtnJlw_m6jVl{b3x4VFu5#0EqFh{=NNek2n0*
zt#+z;ciHffY}&2matOaTUwYz?;)jO2hK8`m2h{Nt9Vyhm%O4#f%Iou86Re~Z3g`dW
zZE}{S)ix5570EA_ZqleAYH9Gtp~T3;Ae}GAG|erzJ=4{-<R`xocqO`hYebXY6&bpk
z2%z?~+UD*#Gz&W(c#QMCB_~*Isa_(K9t8rM@fIea-jXQvt>XHXy^qc)llQbfW*b-z
z%`0kf@drmHqARF<byH|e6QbqHLms|QPTM^q#^UXk*o41i%t-lgOJEg*%>35O(4E}E
z=f>o}){bh@IkW?JsbH%>KC#jVkgH)K;jXkjt)Kt+F{k=qs}b22S6=kymAJwHY!uwU
zz;IVVTrhzpAu|;O7k&Gzk48{eoG|%B0(;)S+vx?c?&`;}%j%asQtQnH$a2@`;OEuc
zb&_0H6Tvp$WgopchYw8n4ns8PKd!fP4iSq|1HhB8%iW}K*Jt0uXBi0hC1j)1tIrd+
z6R}%^-XhS#T@HM>-aMtU{GBCU<K5st2@=G0i(Oc}sPKfSeE`<>(i|uWjUSF^1ow+)
zAB5&*HKK1pbFgHU&t|OgbPx?NIcgPg{z8VB9qy|Cq!{YmvB>ZNwoXL2g1oeZ;7O~7
zv%l4#0?7%F2ooiz{cGB4t$u7y=@G8bNE_ufp)0A=w}#l%12ufg-oHT+#?NTV7qbhy
z?tc>+g&U4@-qXrKAKLuwn#3mSF7cE-@pdtfgbwt3^Xi)ZMl(8S`@YXt+lzosgzJ;L
zuTHiE3w#H+H)!Q+XmBX4Xz#B|DPX1J7j$qip*Sl#TM~<~_`@A8zYo8{CO923#%qJ!
zNGY-aiLv%Xa|62}K5wPl22yqN2)Of?#hum1A8Gwq*%KefCM-nh7l(8ImH<KbQRGVT
z`#W{(BL=`9<-xJFKF!`$fLtHO-3k4wo5kRVt_asZp;$0E!d0dz+tL@hePs>2)~;Cn
z6Kxk8(K#q4sCKut1T7~%Se1`(j#cZ}dogg5xuX2fnoJg0goQYSOKHoeqd;j33wNE{
zJtO|pAK)5-u&b`JQJPgJfP(|gJAKHZu)Z#XWfxJ(Qo%xK&<e7JU08@y-I`b|zx1ea
zR{(yg`XZMPJRH~GVsb}$hy`&?&-sEMIv1K0dl+v9)-J8k_Kc=_Kz$^T_D*Ft{V0;6
z!I7c_R%L))_R0NxRNQj7iwxB83fFCkAY_{|-j=l*++P0GslqTf@|-T8a3gx!CSvzo
zc?10R<Dicz+k%vj*J>wX@=m#fAZC^YA=~q8MR!h<etbhqzL?Xe>_kxh+6SV`7=Cd{
z!|xwx6Wm6GSoVX1?0!LYYD{p)5dnT*=9n7ua-5P1IO!)2fh)ai_mxDX(V&7jDe|y7
z=;gn|;%`Gx$Z9Haa!T>1bku-ztJR3EG}~TNvRyXQivq`fO=;Ly47yWZx+lKo8l&EB
zVv9R)Ol&HUe8R7b6C587_*k38F)mfG-12Y;Y^uim;Gux!L{Plso1y{RD4ci@16${z
z^-mxN6&S^22Tf6|zB%=&V7SX4Ki8zvIdcD%4zk5m`B?~S{NCH~e>8H@6)iLc8bOSR
z85nD223;k|C!`(godf+B=>r$@>O^*#fMrrC%e><ay^MV(lc9Nk0aL703sA0pQ9ays
zn<L01z?#|zY7Q#gbvF&Gjr|c<P;w5-f#94jYW*`?n|oRr;0<hE^&kuB1|Wef*_rQV
z6aIGKIM>b%V3wZJ^J`PMB@C2s<YuGKM|`|Q$tV6KVQ}%xkA1>`YmB}vVHGNv{gB7U
zCZT3G&2D4!Bvv_+!(E$oujY~sW@QDV^@+{mSs-!86H?x+T=eF27#a8Dh_`EhTT!G^
zmn!=9!LMlGRc--4d3RQ2VgD}WiRXgoAA-U6VA&`Hctv48_*Ik`+4d(WTTcP*a9oK9
zw0IB@CfoY|cs_m)^E=BLH!7r%wNC~;)0HF-gG3)S3%$ThW}X6xum}Oe$o=8;cv_S@
zu<x%%Cl<8%@{nzzsK=3CQ+9bj-!{29({&?BaQc7boEKZE6wWS8<2wn^D`viC@pziE
z#w*zO^NoE$i^^XEIDbo*?%rAa21^vDPiMI(cuw9usDQ8RS~;w)%F}#$;T2hy?6S=1
z%&Km^^|FCDL7;7Ta=K1Q(AneB`VnbUDpp+NVMb@yCS&92T8Y?{q)%DMPiD*7x8TO=
z=2qs3)|QEU)x_+;=82aK_{Mhe!RFDzE$H9-Z}S}8jK*#voOcd<Ib$e|nKr)4qo;UY
z<>dL@pOfAvxqe+JXeeqZ|DSux%Kx8xYHCR^+gmr(6#whr6%`d#)&9>tCDs3am7=1u
c!hb&&yr-&^bjr?vOGQDEOG3f`V#xLX09{R)lK=n!

diff --git a/tests/testthat/test-monotonic.R b/tests/testthat/test-monotonic.R
new file mode 100644
index 00000000..a38a3a46
--- /dev/null
+++ b/tests/testthat/test-monotonic.R
@@ -0,0 +1,96 @@
+context("monotonic")
+
+# Simulate data from a monotonically increasing function
+set.seed(123123)
+x <- runif(80) * 4 - 1
+x <- sort(x)
+f <- exp(4 * x) / (1 + exp(4 * x))
+y <- f + rnorm(80) * 0.1
+mod_data <- data.frame(y = y, x = x, z = rnorm(80),
+                       time = 1:80)
+
+test_that("k must be an even integer for s(bs = 'moi')", {
+  expect_error(mvgam(y ~ s(x, bs = 'moi', k = 11),
+                data = mod_data,
+                family = gaussian()),
+               "Argument 'k(bs = 'moi')'  must be an even integer",
+               fixed = TRUE)
+
+  expect_error(mvgam(y ~ s(x, bs = 'moi', k = 1),
+                     data = mod_data,
+                     family = gaussian()),
+               "Basis dimension is too small",
+               fixed = TRUE)
+})
+
+test_that("monotonic only works for one dimensional smooths", {
+  expect_error(mvgam(y ~ s(x, z, bs = 'moi', k = 10),
+                     data = mod_data,
+                     family = gaussian()),
+               "Monotonic basis only handles 1D smooths",
+               fixed = TRUE)
+})
+
+test_that("monotonic for observation models working properly", {
+  mod <- mvgam(y ~ z + s(x, bs = 'moi', k = 18),
+               data = mod_data,
+               family = gaussian(),
+               run_model = FALSE)
+
+  # Monotonic indices should be in the model_data
+  expect_true("b_idx_s_x_" %in% names(mod$model_data))
+
+  # The smooth should be an MOI class
+  expect_true(inherits(mod$mgcv_model$smooth[[1]], 'moi.smooth'))
+
+  # The coefficients should be fixed to be non-negative
+  expect_true(any(grepl('b[b_idx_s_x_] = abs(b_raw[b_idx_s_x_]) * 1;',
+                        mod$model_file, fixed = TRUE)))
+
+  # Repeat a check for decreasing functions
+  mod <- mvgam(y ~ z + s(x, bs = 'mod', k = 18),
+               data = mod_data,
+               family = gaussian(),
+               run_model = FALSE)
+
+  # The smooth should be an MOD class
+  expect_true(inherits(mod$mgcv_model$smooth[[1]], 'mod.smooth'))
+
+  # The coefficients should be fixed to be non-positive
+  expect_true(any(grepl('b[b_idx_s_x_] = abs(b_raw[b_idx_s_x_]) * -1;',
+                        mod$model_file, fixed = TRUE)))
+})
+
+test_that("monotonic for process models working properly", {
+  mod <- mvgam(y ~ 0,
+               trend_formula = ~ z + s(x, bs = 'moi', k = 18),
+               trend_model = RW(),
+               data = mod_data,
+               family = gaussian(),
+               run_model = FALSE)
+
+  # Monotonic indices should be in the model_data
+  expect_true("b_trend_idx_s_x_" %in% names(mod$model_data))
+
+  # The smooth should be an MOI class
+  expect_true(inherits(mod$trend_mgcv_model$smooth[[1]], 'moi.smooth'))
+
+  # The coefficients should be fixed to be non-negative
+  expect_true(any(grepl('b_trend[b_trend_idx_s_x_] = abs(b_raw_trend[b_trend_idx_s_x_]) * 1;',
+                        mod$model_file, fixed = TRUE)))
+
+  # And for decreasing
+  mod <- mvgam(y ~ 0,
+               trend_formula = ~ z + s(x, bs = 'mod', k = 18),
+               trend_model = RW(),
+               data = mod_data,
+               family = gaussian(),
+               run_model = FALSE)
+
+  # The smooth should be an MOD class
+  expect_true(inherits(mod$trend_mgcv_model$smooth[[1]], 'mod.smooth'))
+
+  # The coefficients should be fixed to be non-positive
+  expect_true(any(grepl('b_trend[b_trend_idx_s_x_] = abs(b_raw_trend[b_trend_idx_s_x_]) * -1;',
+                        mod$model_file, fixed = TRUE)))
+})