R Under development (unstable) (2025-04-13 r88141 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> require(sparseFLMM)
Loading required package: sparseFLMM
Loading required package: mgcv
Loading required package: nlme
This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
Loading required package: refund
> 
> # (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")
> 
> proc.time()
   user  system elapsed 
   3.09    0.37    3.45