require(sparseFLMM)

# (skew-)symmetric smooths ---------------------------------------

# generate random surface 
dat1 <- data.frame(arg1 = 1:50)
dat2 <- expand.grid(arg1 = 1:50, arg2 = 1:50)

Bskew <- Predict.matrix(
  smooth.construct( 
    s(arg1, arg2, bs = "symm", xt = list(skew = TRUE)),
    data = dat2, knots = NULL ),
  data = dat2 )
Bsymm <- Predict.matrix(
  smooth.construct( 
    s(arg1, arg2, bs = "symm", xt = list(skew = FALSE)),
    data = dat2, knots = NULL ),
  data = dat2 )

set.seed(934811)
dat2$yskew <- c(Bskew %*% rnorm(ncol(Bskew)))
dat2$ysymm <- c(Bsymm %*% rnorm(ncol(Bsymm)))

# fit sum of skew-symmetric and symmetric parts with corresponding smooths
modpa <- gam( I(yskew + ysymm) ~ s(arg1, arg2, bs = "symm", xt = list(skew = TRUE)) + 
                s(arg1, arg2, bs = "symm", xt = list(skew = FALSE)), data = dat2)
# predict surfaces
preds <- predict(modpa, type = "terms")
dat1 <- as.list(dat1)
dat1$arg2 <- dat1$arg1
dat1$predskew <- matrix(preds[,1], nrow = length(dat1$arg1))
dat1$predsymm <- matrix(preds[,2], nrow = length(dat1$arg1))

cols <- hcl.colors(12, "RdBu")
opar <- par(mfcol = c(2,2))
# symm part (intercept missing)
with(dat1, image(arg1, arg2, predsymm, asp = 1,
                 main = "Symmetric part of y",
                 col = cols))
with(dat1, image(arg1, arg2, asp = 1, 
                 main = "Fit via symm.smooth",
                 matrix(dat2$ysymm, nrow = length(arg1)), 
                                    col = cols))
# skew-symm part
with(dat1, image(arg1, arg2, predskew, asp = 1,
                 main = "Skew-symmetric part of y",
                 col = cols))
with(dat1, image(arg1, arg2, asp = 1, 
                 main = "Fit via symm.smooth",
                 matrix(dat2$yskew, nrow = length(arg1)), 
                 col = cols))
par(opar)

stopifnot(all.equal(dat1$predskew, - t(dat1$predskew)))
stopifnot(all.equal(dat1$predsymm, t(dat1$predsymm)))




# cyclic (skew-)symmetric splines ---------------------------------------

# fit the above example with cyclic smooths
modpac <- gam( I(yskew + ysymm) ~ s(arg1, arg2, bs = "symm", 
                                   xt = list(skew = TRUE, cyclic = TRUE)) + 
                s(arg1, arg2, bs = "symm", xt = list(skew = FALSE, cyclic = TRUE)),
              knots = list(arg1 = c(1, 50), arg2 = c(1,50)), 
              # specify arg range to specify 'wavelength'! 
              data = dat2)
plot(modpac, asp = 1, se = FALSE, pages = 1)

predsc <- predict(modpac, type = "terms")
dat1$predskewc <- matrix(predsc[,1], nrow = length(dat1$arg1))
dat1$predsymmc <- matrix(predsc[,2], nrow = length(dat1$arg1))

# check cyclic margins
opar <- par(mfrow = c(1,2))
with(dat1, matplot(arg1, predsymmc[, c(1,10, 40)], t = "l",
                   main = "symmetric smooth"))
abline(h = dat1$predsymmc[1, c(1,10, 40)], col = "darkgrey")
abline(v = c(1,50), col = "darkgrey")

with(dat1, matplot(arg1, predskewc[, c(1,10, 40)], t = "l",
                   main = "skew-symmetric smooth"))
abline(h = dat1$predskewc[1, c(1,10, 40)], col = "darkgrey")
abline(v = c(1,50), col = "darkgrey")
par(opar)



# 1D point symmetric B-splines --------------------------------------------

# generate toy data
dat <- data.frame( x = 1:100 )
ps_obj <- with(dat, s(x, bs = "ps"))
B <- Predict.matrix(smooth.construct(ps_obj, dat, NULL), dat)               
set.seed(3904)
dat$y <- B %*% rnorm(ncol(B))
plot(dat, t = "l")

# fit skew-symmetric spline
mod0 <- gam( y ~ s(x, bs = "symm", xt = list(skew = TRUE)), 
             knots = list(x = c(0,100)), # specify x range to determine inversion point 
             dat = dat )
lines(dat$x, predict(mod0), col = "cornflowerblue", lty = "dashed")

# or a symmetric spline to first part only
mod1 <- gam( y ~ s(x, bs = "symm"), 
             knots = list(x=c(0,50)), 
                          dat = dat[1:50, ])
lines(dat[1:50, ]$x, predict(mod1), col = "darkred", lty = "dashed")