L4 data, p = 22 (# composition covariates), m = 6 (# control covariates) and include intercept. Time ranges from PNA = 5 to PNA = 30 and are mapped to standard interval [0,1]. Z (integral of compositional data) is scaled to sd = 1 (column-wise). Select from k_list (degree freedom of cubic splines) in $[4, \dots, 9]$, and lambda.factor is set to 0.1.

library(compReg)
#try$Funcomp.CGL.fit[[1]]
input <- get(load("C:/Users/zjiji/Dropbox/Shared-ZheSun/FuncCompReg/Real Data/L4/graphcode/leaveoneout/data/input_K=4.Rda"))
#apply(input$X / 25, 2, sd)
try <- get(load("C:/Users/zjiji/Dropbox/Shared-ZheSun/FuncCompReg/Real Data/L4/graphcode/leaveoneout/data/[0,1]10folds_Weighted_cv.FuncompCGL-0.1.Rda"))
#names(try)
#as.numeric(names(try[[1]]))

try2 <- get(load("C:/Users/zjiji/Dropbox/Shared-ZheSun/FuncCompReg/Real Data/L4/graphcode/leaveoneout/data/[0,1]cv.FuncompCGL.Rda"))
# try2$lam[1]
# length(try2$lam)
# try2$lam[1]/tail(try2$lam, 1)
# class(try2$Funcomp.CGL.fit$`4`$W)
# dim(try2$Funcomp.CGL.fit$`4`$W)
# head(try2$Funcomp.CGL.fit$`4`$Z)
# head(try$Funcomp.CGL.fit$`4`$Z)

k_list <- 4:9
GIC <-as.list(1:length(k_list))
a <- k_list[1] - 1
n <- length(input$y)
p <- 22


# ## with n 34, p = 22, p * dk >= n, p < n
# ## GIC1 = GIC3, GIC2 = GIC4, GIC1.1 = GIC3.1, GIC1.2 = GIC3.2, GIC2.1 = GIC4.1, GIC2.2 = GIC4.2, GIC1.3 = GIC3.3, GIC2.3 = GIC4.3
# ## GIC5 = BIC1, GIC6 = BIC2, GIC5.1 = BIC1.1, GIC6.1 = BIC2.1, GIC5.2 = BIC1.2, GIC6.2 = BIC2.2, GIC5.3 = BIC1.2, GIC6.3 = BIC2.3
# ##
# ## alpha = switch(GIC_type,
# ##               "GIC1" = log(max(p*dk,n)) * dk,
# ##               "GIC2" = log(max(p*dk,n)) * dk,
# ##               "GIC3" = log(p * dk) * dk,
# ##               "GIC4" = log(p * dk) * dk,
# ##               "GIC1.1" = log(max(p*dk, n) * dk),
# ##               "GIC2.1" = log(max(p*dk, n) * dk),
# ##               "GIC3.1" = log(p*dk * dk),
# ##               "GIC4.1" = log(p*dk * dk),
# ##               "GIC1.2" = log(max(p*dk, n)) * log(dk),
# ##               "GIC2.2" = log(max(p*dk, n)) * log(dk),
# ##               "GIC3.2" = log(p*dk) * log(dk),
# ##               "GIC4.2" = log(p*dk) * log(dk),
# ##               "GIC1.3" = log(max(p * dk, n)),
# ##               "GIC2.3" = log(max(p * dk, n)),
# ##               "GIC3.3" = log(p * dk),
# ##               "GIC4.3" = log(p * dk),
# ##               "GIC5" = log(max(p, n)) * dk,
# ##               "GIC6" = log(max(p, n)) * dk,
# ##               "GIC7" = log(p) * dk,
# ##               "GIC8" = log(p) * dk,
# ##               "GIC5.1" = log(max(p, n) * dk),
# ##               "GIC6.1" = log(max(p, n) * dk),
# ##               "GIC5.2" = log(max(p, n)) * log(dk),
# ##               "GIC6.2" = log(max(p, n)) * log(dk),
# ##               "GIC7.2" = log(p) * log(dk),
# ##               "GIC8.2" = log(p) * log(dk),
# ##               "GIC5.3" = log(max(p, n)),
# ##               "GIC6.3" = log(max(p, n)),
# ##               "GIC7.3" = log(p),
# ##               "GIC8.3" = log(p),
# ##               "AIC1" = dk * 2,
# ##               "AIC2" = log(dk) * 2,
# ##               "BIC1" = dk * log(n),
# ##               "BIC2" = dk * log(n),
# ##               "BIC1.1" = log(dk * n),
# ##               "BIC2.1" = log(dk * n),
# ##               "BIC1.2" = log(dk) * log(n),
# ##               "BIC2.2" = log(dk) * log(n),
# ##               "BIC1.3" = log(n),
# ##               "BIC2.3" = log(n)
# ## ) / n
# ## even version of GICs and BICs multiplying log(log(n))
# ## odd version of GICs and BICs multiplying 1

GIC_typelist = c(paste0("GIC", 1:2), paste0("GIC1.", 1:3), paste0("GIC2.", 1:3), 
                 paste0("GIC", 7:8), paste0("GIC7.", 2:3), paste0("GIC8.", 2:3),
                 paste0("BIC", 1:2), paste0("BIC1.", 1:3), paste0("BIC2.", 1:3),
                 "AIC1", "AIC2")

MSE <- matrix(NA, nrow = length(k_list), ncol = length(try$lam))
N_zero <- matrix(NA, nrow = length(k_list), ncol = length(try$lam))
GIC_min <- vector()
scaler <- vector()
knitr::opts_chunk$set(cache=TRUE)

\newpage

10 folds corss validation

#try$Ftrim$cvm
#plot(try, k_list = k_list)

par(mfrow=c(1,2))

plot.args_CV = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                    y = try$Ftrim$cvm[1, ],
                    ylim = range(try$Ftrim$cvm),
                    xlab= "lambda Index",#"lambda", #"log(lambda)",
                    ylab="MEAN-Squared Error",
                    main = "10-folds Cross Vadilation",
                    type="n")
do.call("plot",plot.args_CV)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = try$Ftrim$cvm[i, ], 
         col = rainbow(length(k_list))[i], 
         pch = seq(length(k_list))[i])

  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       try$Ftrim$cvm[i, ][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which.min(try$Ftrim$cvm[i, ]), col = rainbow(length(k_list))[i])
}

CV_min <- apply(try$Ftrim$cvm, 1, min)
k_opt <- which.min(CV_min) + a
lam_opt <- which(try$Ftrim$cvm[k_opt - a, ] == CV_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

abline(a = 0, b = 0, col = "grey", lwd = 2)
title(paste("Groups:", paste0(paste(Non.zero1, collapse = ","), ","), "df=", k_opt), 
      line = 0.5)
text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
text(tail(seq(0, 1, length.out = 20), 1), 
     beta_curve1[dim(beta_curve1)[1], Non.zero1],
    labels = paste(Non.zero1))
par(mfrow=c(1,1))

\newpage

Leave-one-out corss validation

par(mfrow=c(1,2))

plot.args_CV = list(x = seq(length(try2$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                    y = try2$Ftrim$cvm[1, ],
                    ylim = range(try2$Ftrim$cvm),
                    xlab= "lambda Index",#"lambda", #"log(lambda)",
                    ylab="MEAN-Squared Error",
                    main = "Leave one out Cross Vadilation",
                    type="n")
do.call("plot",plot.args_CV)
for(i in 1:length(k_list)) {
  points(x = seq(length(try2$lam)), #log(GIC_m1$lam),
         y = try2$Ftrim$cvm[i, ], 
         col = rainbow(length(k_list))[i], 
         pch = seq(length(k_list))[i])

  text(length(try2$lam), #tail(log(GIC_m1$lam), 1),
       try2$Ftrim$cvm[i, ][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which.min(try2$Ftrim$cvm[i, ]), col = rainbow(length(k_list))[i])
}

CV_min <- apply(try2$Ftrim$cvm, 1, min)
k_opt <- which.min(CV_min) + a
lam_opt <- which(try2$Ftrim$cvm[k_opt - a, ] == CV_min[k_opt-a])
beta1 <-  try2$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

abline(a = 0, b = 0, col = "grey", lwd = 2)
title(paste("Groups:", paste0(paste(Non.zero1, collapse = ","), ","), "df=", k_opt), 
      line = 0.5)
text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
text(tail(seq(0, 1, length.out = 20), 1), 
     beta_curve1[dim(beta_curve1)[1], Non.zero1],
    labels = paste(Non.zero1))
par(mfrow=c(1,1))

\newpage

GIC1

$\textrm{GIC}1(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \times df \times \log\big(\max(p\times df + m + 1,n)\big)+ (m+1) \log\big(\max(p\times df + m + 1,n)\big)$

GIC_type = GIC_typelist[1]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste0(paste(Non.zero1, collapse = ","), ","), "df=", k_opt), 
      line = 0.5)
text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
text(tail(seq(0, 1, length.out = 20), 1), 
     beta_curve1[dim(beta_curve1)[1], Non.zero1],
    labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

GIC2

$\textrm{GIC}2(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \times df \times \log\big(\max(p\times df + m + 1,n)\big)\log(\log(n))+ (m+1) \log\big(\max(p\times df + m + 1,n)\big)\log(\log(n))$

GIC_type = GIC_typelist[2]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

GIC1.1

$\textrm{GIC}{1.1}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log\big(\max(p\times df + m + 1,n) \times df \big)+ (m+1) \log\big(\max(p\times df + m + 1,n)\big)$

GIC_type = GIC_typelist[3]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

GIC1.2

$\textrm{GIC}{1.2}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log\big(\max(p\times df + m + 1,n) \big) \log(df) + (m+1) \log\big(\max(p\times df + m + 1,n)\big)$

GIC_type = GIC_typelist[4]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

GIC1.3

$\textrm{GIC}{1.3}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log\big(\max(p\times df + m + 1,n) \big)+ (m+1) \log\big(\max(p\times df + m + 1,n)\big)$

GIC_type = GIC_typelist[5]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

GIC2.1

$\textrm{GIC}{2.1}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log\big(\max(p\times df + m + 1,n) \times df \big)\log(\log(n))+ (m+1) \log\big(\max(p\times df + m + 1,n)\big)\log(\log(n))$

GIC_type = GIC_typelist[6]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

GIC2.2

$\textrm{GIC}{2.2}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log\big(\max(p\times df + m + 1,n) \big) \log(df) \log(\log(n)) + (m+1) \log\big(\max(p\times df + m + 1,n)\big) \log(\log(n))$

GIC_type = GIC_typelist[7]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

GIC2.3

$\textrm{GIC}{2.3}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log\big(\max(p\times df + m + 1,n) \big)\log(\log(n))+ (m+1) \log\big(\max(p\times df + m + 1,n)\big) \log(\log(n))$

GIC_type = GIC_typelist[8]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

BIC1

$\textrm{BIC}1(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n) \times df + (m+1) \log(n)$

GIC_type = GIC_typelist[15]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

BIC2

$\textrm{BIC}2(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n) \times df \times \log(\log(n))+ (m+1) \log(n) \log(\log(n))$

GIC_type = GIC_typelist[16]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

BIC1.1

$\textrm{BIC}{1.1}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n \times df) + (m+1) \log(n)$

GIC_type = GIC_typelist[17]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

BIC1.2

$\textrm{BIC}{1.2}(\lambda) = n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n ) \times \log(df) + (m+1) \log(n)$

GIC_type = GIC_typelist[18]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

BIC1.3

$\textrm{BIC}{1.3}(\lambda) = n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n ) + (m+1) \log(n)$

GIC_type = GIC_typelist[19]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

BIC2.1

$\textrm{BIC}{2.1}(\lambda) =n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n \times df)\log(\log(n)) + (m+1) \log(n)\log(\log(n))$

GIC_type = GIC_typelist[20]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

BIC2.2

$\textrm{BIC}{2.2}(\lambda) = n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n ) \times \log(df)\log(\log(n)) + (m+1) \log(n)\log(\log(n))$

GIC_type = GIC_typelist[21]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))

\newpage

BIC2.3

$\textrm{BIC}{2.3}(\lambda) = n\log(\widehat{\textrm{MSE}}) + 1{s_{\lambda} \ge 2} \big(s_{\lambda - 1} \big) \log(n ) \log(\log(n)) + (m+1) \log(n)\log(\log(n))$

GIC_type = GIC_typelist[23]


for(k in k_list){
  pre <- cbind(try$Funcomp.CGL.fit[[k-a]]$Z, input$Z) %*% try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ]
  MSE[k-a, ] <- apply(input$y - pre, 2, function(x) mean(x^2))
  N_zero[k-a, ] <- apply(try$Funcomp.CGL.fit[[k-a]]$beta[1:(k*p + 7), ], 2,                                  function(x, p, k) {
                          A = matrix(x[1:(p*k)], nrow = p, byrow = TRUE)
                          a = apply(A, 1, function(x) sum(abs(x)))
                          return(length(which(a > 0)))
                          },p = p, k = k)
  # scaler <- log(max(p * k, n)) / n * k # *log(log(n)) 
  dk = k

   alpha = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) * dk,
                   "GIC2" = log(max(p*dk + 7, n)) * dk,
                   "GIC3" = log(p * dk + 7) * dk,
                   "GIC4" = log(p * dk + 7) * dk,
                   "GIC1.1" = log(max(p*dk + 7, n) * dk),
                   "GIC2.1" = log(max(p*dk + 7, n) * dk),
                   "GIC3.1" = log((p*dk+7) * dk),
                   "GIC4.1" = log((p*dk+7) * dk),
                   "GIC1.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC2.2" = log(max(p*dk + 7, n)) * log(dk),
                   "GIC3.2" = log(p*dk + 7) * log(dk),
                   "GIC4.2" = log(p*dk + 7) * log(dk),
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)) * dk,
                   "GIC6" = log(max(p + 7, n)) * dk,
                   "GIC7" = log(p + 7) * dk,
                   "GIC8" = log(p + 7) * dk,
                   "GIC5.1" = log(max(p + 7, n) * dk),
                   "GIC6.1" = log(max(p + 7, n) * dk),
                   "GIC5.2" = log(max(p + 7, n)) * log(dk),
                   "GIC6.2" = log(max(p + 7, n)) * log(dk),
                   "GIC7.2" = log(p + 7) * log(dk),
                   "GIC8.2" = log(p + 7) * log(dk),
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = dk * 2,
                   "AIC2" = log(dk) * 2,
                   "BIC1" = dk * log(n),
                   "BIC2" = dk * log(n),
                   "BIC1.1" = log(dk * n),
                   "BIC2.1" = log(dk * n),
                   "BIC1.2" = log(dk) * log(n),
                   "BIC2.2" = log(dk) * log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n



   alpha2 = switch(GIC_type,
                   "GIC1" = log(max(p*dk + 7, n)) ,
                   "GIC2" = log(max(p*dk + 7, n)) ,
                   "GIC3" = log(p * dk + 7) ,
                   "GIC4" = log(p * dk + 7) ,
                   "GIC1.1" = log(max(p*dk + 7, n)),
                   "GIC2.1" = log(max(p*dk + 7, n)),
                   "GIC3.1" = log((p*dk+7) ),
                   "GIC4.1" = log((p*dk+7) ),
                   "GIC1.2" = log(max(p*dk + 7, n)) ,
                   "GIC2.2" = log(max(p*dk + 7, n)) ,
                   "GIC3.2" = log(p*dk + 7) ,
                   "GIC4.2" = log(p*dk + 7) ,
                   "GIC1.3" = log(max(p * dk + 7, n)),
                   "GIC2.3" = log(max(p * dk + 7, n)),
                   "GIC3.3" = log(p * dk + 7),
                   "GIC4.3" = log(p * dk + 7),
                   "GIC5" = log(max(p + 7, n)),
                   "GIC6" = log(max(p + 7, n)),
                   "GIC7" = log(p + 7),
                   "GIC8" = log(p + 7),
                   "GIC5.1" = log(max(p + 7, n)),
                   "GIC6.1" = log(max(p + 7, n) ),
                   "GIC5.2" = log(max(p + 7, n)),
                   "GIC6.2" = log(max(p + 7, n)) ,
                   "GIC7.2" = log(p + 7) ,
                   "GIC8.2" = log(p + 7) ,
                   "GIC5.3" = log(max(p + 7, n)),
                   "GIC6.3" = log(max(p + 7, n)),
                   "GIC7.3" = log(p + 7),
                   "GIC8.3" = log(p + 7),
                   "AIC1" = 2,
                   "AIC2" = 2,
                   "BIC1" = log(n),
                   "BIC2" = log(n),
                   "BIC1.1" = log(n),
                   "BIC2.1" = log(n),
                   "BIC1.2" = log(n),
                   "BIC2.2" = log(n),
                   "BIC1.3" = log(n),
                   "BIC2.3" = log(n)
    ) / n

    if(GIC_type %in% c("GIC2", "GIC4", "GIC2.1", "GIC4.1", "GIC2.2", "GIC4.2", "GIC2.3", "GIC4.3",
                       "GIC6", "GIC8", "GIC6.1", "GIC6.2", "GIC8.2",  "GIC6.3", "GIC8.3",
                       "BIC2", "BIC2.1", "BIC2.2", "BIC2.3")) {
      c_const = log(log(n))
    } else {
      c_const = 1
    }


    scaler[k-a] <- alpha * c_const

    part2 <- scaler[k-a] * ifelse(N_zero[k-a, ] >= 2, N_zero[k-a, ] - 1, 0) + (alpha2 * c_const) * 7
  #part2 <- scaler * ifelse(N_zero[k-a, ] >= 2, 
  #                        (N_zero[k-a, ] - 1) * k + 7, 7)
  #plot(log(MSE[k-a, ]) + part2)
  #which.min((log(MSE[k - k_list[1] + 1, ])+part2)[N_zero >= 2])

  GIC[[k-a]] <- log(MSE[k-a, ]) + part2
  GIC_min[k-a] <- min(GIC[[k-a]][which(N_zero[k-a, ] >= 2)])

}
cat("\r\n")


par(mfrow=c(1,2))

plot.args_GIC = list(x = seq(length(try$lam)), #GIC_m1$lam, #log(GIC_m1$lam),
                     y = GIC[[1]],
                     ylim = range(GIC),
                     xlab= "lambda Index",#"lambda", #"log(lambda)",
                     ylab="GIC",
                     main = "GIC",
                     type="n")
do.call("plot",plot.args_GIC)
for(i in 1:length(k_list)) {
  points(x = seq(length(try$lam)), #log(GIC_m1$lam),
         y = GIC[[i]], col = rainbow(length(k_list))[i], pch = seq(length(k_list))[i])
  text(length(try$lam), #tail(log(GIC_m1$lam), 1),
       GIC[[i]][length(try$lam)], labels=paste(k_list[i]),
       cex= 1, pos= 4, col = rainbow(length(k_list))[i])
  abline(v = which(GIC[[i]] == GIC_min[i]), col = rainbow(length(k_list))[i])
}


k_opt <- which.min(GIC_min) + a
lam_opt <- which(GIC[[k_opt - a]] == GIC_min[k_opt-a])
beta1 <-  try$Funcomp.CGL.fit[[as.character(k_opt)]]$beta[1:(k_list[k_opt - a]*p + 7), lam_opt]
beta1_C <- matrix(beta1[1:(p*k_opt)], nrow = p, byrow = TRUE)
Non.zero1 = (1:p)[which(apply(beta1_C, 1, function(x) sum(abs(x)) > 0))]

B1 <- splines::bs(seq(0, 1, length.out = 20), df = k_opt, intercept = TRUE)
beta_curve1 <- B1 %*% t(beta1_C)
matplot(seq(0, 1, length.out = 20), beta_curve1,
        ylab = "Estimated Curve", xlab = "TIME", #main = "ESTI",
        ylim = c(-max(abs(beta_curve1)), max(abs(beta_curve1)))#,
        ,type = "l"
)

 abline(a = 0, b = 0, col = "grey", lwd = 2)
 title(paste("Groups:", paste(Non.zero1, collapse = ",")), line = 0.5)
 text(0, beta_curve1[1, Non.zero1], labels = paste(Non.zero1), pos = 1)
 # text(tail(sseq, 1), beta_curve_temp[dim(beta_curve_temp)[1], Non.zero1],
 #      labels = paste(Non.zero1))
 par(mfrow=c(1,1))


jiji6454/Rpac_compReg documentation built on May 31, 2019, 5:01 a.m.