From 82f48b6773b9f5028d422abc5bfb065fd8689606 Mon Sep 17 00:00:00 2001
From: Matteo Fasiolo <mf15002@it031161.wks.bris.ac.uk>
Date: Fri, 29 Mar 2024 15:22:21 +0000
Subject: [PATCH] mgks now based on distances and treated as "si" when inner ==
 TRUE

---
 R/I_prepareInnerNested.R | 115 +++++++++++++++++++++------------------
 R/I_prepareOuterNested.R |   2 +-
 2 files changed, 64 insertions(+), 53 deletions(-)

diff --git a/R/I_prepareInnerNested.R b/R/I_prepareInnerNested.R
index 2d56eaa..a2e7dd2 100644
--- a/R/I_prepareInnerNested.R
+++ b/R/I_prepareInnerNested.R
@@ -13,25 +13,26 @@
   
   da <- length( alpha )
   prange <- (sm$first.para:sm$last.para)[1:da]
-  Va <- gObj$Vp[prange, prange, drop = FALSE]
-  
   type <- class(o)[1]
-  if(type == "si"){
-    alpha <- drop(B %*% alpha)
-    Va <- B %*% Va %*% t(B)
-    se <- sqrt(pmax(0, diag(Va)))
-    edf   <- sum(gObj$edf[prange])
-    ylabel <- .subEDF(paste0("proj_coef(", sm$term, ")"), edf)
-    xlabel <- "Index"
-    out <- list("fit" = alpha, "x" = 1:da, "se" = se,
-                xlab = xlabel, ylab = ylabel, main = NULL, type = "si")
+  
+  # Remove scale parameter of inner transformation
+  if( type != "si" ){
+    prange <- prange[-1]
+    da <- da-1
+    alpha <- alpha[-1]
+    if(is.null(B)){
+      B <- diag(1, nrow = length(alpha))
+    }
   }
+  
+  Va <- gObj$Vp[prange, prange, drop = FALSE]
+  
   if( type == "nexpsm" ){
-    inner <- expsmooth(y = si$x, Xi = si$X, beta = alpha[-1], deriv = 1)
+    inner <- expsmooth(y = si$x, Xi = si$X, beta = alpha, deriv = 1)
     fit <- inner$d0
     Jac <- inner$d1
     nobs <- length(fit)
-    se <- sqrt(pmax(0, rowSums((Jac %*% Va[-1, -1, drop = FALSE]) * Jac)))
+    se <- sqrt(pmax(0, rowSums((Jac %*% Va) * Jac)))
     edf   <- sum(gObj$edf[prange[-1]])
     ylabel <- .subEDF(paste0("expsm(", sm$term, ")"), edf)
     xlabel <- "Index"
@@ -49,46 +50,56 @@
                 "p.resid" = si$x[ii], "raw" = ii, 
                 "xlim" = xlim, 
                 xlab = xlabel, ylab = ylabel, main = NULL, type = "nexpsm")
+  } else {
+    alpha <- drop(B %*% alpha)
+    Va <- B %*% Va %*% t(B)
+    se <- sqrt(pmax(0, diag(Va)))
+    edf   <- sum(gObj$edf[prange])
+    ylabel <- .subEDF(paste0("Inner_coef(", sm$term, ")"), edf)
+    xlabel <- "Index"
+    out <- list("fit" = alpha, "x" = 1:da, "se" = se,
+                xlab = xlabel, ylab = ylabel, main = NULL, type = "si")
   }
-  if( type == "mgks" ){
-    d <- ncol(si$X0)
-    if( d != 2 ){ # ONLY 2D case handled at the moment!!
-      return( NULL )
-    }
-    if( !is.null(xlim) ) {
-      xlim <- sort(xlim)
-    } else {
-      xlim <- range(si$X[ , 1])
-    }
-    if( !is.null(ylim) ) {
-      ylim <- sort(ylim)
-    } else {
-      ylim <- range(si$X[ , 2])
-    }
-    
-    xx <- rep(seq(xlim[1], xlim[2], length.out = n), n)
-    yy <- rep(seq(ylim[1], ylim[2], length.out = n), rep(n, n))
-    
-    X <- cbind(xx, yy)
-    si$x <- as.matrix(si$x)
-    if( ncol(si$x) > 1 ){
-      si$x <- colMeans(si$x)
-    }
-    
-    inner <- mgks(y = si$x, X = X, X0 = si$X0, beta = alpha[-1], deriv = 1)
-    fit <- inner$d0
-    Jac <- inner$d1
-    se <- sqrt(pmax(0, rowSums((Jac %*% Va[-1, -1, drop = FALSE]) * Jac)))
-    edf   <- sum(gObj$edf[prange[-1]])
-    
-    mainlab <- .subEDF(paste0("mgks(", sm$term, ")"), edf)
-    ylabel <- "X[ , 2]"
-    xlabel <- "X[ , 1]"
-    out <- list("fit" = fit, "X" = si$X, "se" = se, "x" = xx, "y" = yy,
-                "p.resid" = si$x, "X0" = si$X0, 
-                "xlim" = xlim, "ylim" = ylim,
-                "xlab" = xlabel, "ylab" = ylabel, "main" = mainlab, type = "mgks")
-  }
+  # NOT CLEAR HOW TO DO THIS WITH DISTANCEs
+  # if( type == "mgks" ){
+  #   d <- ncol(si$X0)
+  #   if( d != 2 ){ # ONLY 2D case handled at the moment!!
+  #     return( NULL )
+  #   }
+  #   if( !is.null(xlim) ) {
+  #     xlim <- sort(xlim)
+  #   } else {
+  #     xlim <- range(si$X[ , 1])
+  #   }
+  #   if( !is.null(ylim) ) {
+  #     ylim <- sort(ylim)
+  #   } else {
+  #     ylim <- range(si$X[ , 2])
+  #   }
+  #   
+  #   xx <- rep(seq(xlim[1], xlim[2], length.out = n), n)
+  #   yy <- rep(seq(ylim[1], ylim[2], length.out = n), rep(n, n))
+  #   
+  #   X <- cbind(xx, yy)
+  #   si$x <- as.matrix(si$x)
+  #   if( ncol(si$x) > 1 ){
+  #     si$x <- colMeans(si$x)
+  #   }
+  #   
+  #   inner <- mgks(y = si$x, X = X, X0 = si$X0, beta = alpha[-1], deriv = 1)
+  #   fit <- inner$d0
+  #   Jac <- inner$d1
+  #   se <- sqrt(pmax(0, rowSums((Jac %*% Va[-1, -1, drop = FALSE]) * Jac)))
+  #   edf   <- sum(gObj$edf[prange[-1]])
+  #   
+  #   mainlab <- .subEDF(paste0("mgks(", sm$term, ")"), edf)
+  #   ylabel <- "X[ , 2]"
+  #   xlabel <- "X[ , 1]"
+  #   out <- list("fit" = fit, "X" = si$X, "se" = se, "x" = xx, "y" = yy,
+  #               "p.resid" = si$x, "X0" = si$X0, 
+  #               "xlim" = xlim, "ylim" = ylim,
+  #               "xlab" = xlabel, "ylab" = ylabel, "main" = mainlab, type = "mgks")
+  # }
   
   return(out)
   
diff --git a/R/I_prepareOuterNested.R b/R/I_prepareOuterNested.R
index dbca5c9..8170e02 100644
--- a/R/I_prepareOuterNested.R
+++ b/R/I_prepareOuterNested.R
@@ -23,7 +23,7 @@
    trnam <- "expsm"
   }
   if( type == "mgks" ){
-    raw <- mgks(y = si$x, X = si$X, X0 = si$X0, beta = alpha[-1])$d0
+    raw <- mgks(y = si$x, dist = si$dist, beta = alpha[-1])$d0
     trnam <- "mgks"
   }