## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ---- warning=F, message=FALSE, eval=TRUE-------------------------------------
library(KOFM)

## -----------------------------------------------------------------------------
simple_tensor <- array(1:24, dim=c(3,4,2))

## -----------------------------------------------------------------------------
cat('Reshape along modes 1 and 2:\n')
ten_reshape(ten=simple_tensor, AA=c(1,2), time.mode=FALSE)
cat('\nThe same results:\n')
ten_reshape(ten=simple_tensor, AA=c(2,1), time.mode=FALSE)

## -----------------------------------------------------------------------------
cat('Reshape along modes 2 and 3:\n')
simple_tensor_23 <- ten_reshape(ten=simple_tensor, AA=c(2,3), time.mode=FALSE)
simple_tensor_23
cat('\nReshape along modes 1, 2 and 3:\n')
ten_reshape(ten=simple_tensor, AA=c(1,2,3), time.mode=FALSE)

## -----------------------------------------------------------------------------
cat('With the time series "simple_tensor", reshape along modes 1 and 2:\n')
ten_reshape(ten=simple_tensor, AA=c(1,2), time.mode=TRUE)

## -----------------------------------------------------------------------------
cat('Recover "simple_tensor" by reshaping back "simple_tensor_23":\n')
simple_tensor_recover <- ten_reshape_back(ten=simple_tensor_23, AA=c(2,3), original.dim=dim(simple_tensor), time.mode=FALSE)
simple_tensor_recover
cat(paste0('\nTo see that "simple_tensor" is indeed recovered, the sum of their squared entry differences is: ', sum((simple_tensor_recover - simple_tensor)^2), '.'))

## -----------------------------------------------------------------------------
K <- 3 #order-3 tensor observed at each time
TT <- 100 #time length
d <- c(10,10,10) #spatial dimensions
r <- c(2,2,2) #rank of core tensor
re <- c(2,2,2) #rank of common component in error
eta <- list(c(0,0), c(0,0), c(0,0)) #strong factors
coef_f <- c(0.3) #AR(1) coefficients for core factor
coef_fe <- c(-0.3) #AR(1) coefficients for common component in error
coef_e <- c(0.5) #AR(1) coefficients for idiosyncratic part in error
data_H0 <- tensorMiss::tensor_gen(K, TT, d, r, re, eta, coef_f, coef_fe, coef_e, seed=2025)$X

## -----------------------------------------------------------------------------
cat('Test "data_H0" along modes 1 and 2:\n')
testKron_A(ten=data_H0, AA=c(1,2), r=r, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE)

## -----------------------------------------------------------------------------
cat('Test "data_H0" along modes 1 and 3:\n')
apply(testKron_A(ten=data_H0, AA=c(1,3), r=r, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE)$decision, 2, min)
cat('\nTest "data_H0" along modes 2 and 3:\n')
apply(testKron_A(ten=data_H0, AA=c(2,3), r=r, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE)$decision, 2, min)
cat('\nTest "data_H0" along modes 1, 2 and 3:\n')
apply(testKron_A(ten=data_H0, AA=c(1,2,3), r=r, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE)$decision, 2, min)

## -----------------------------------------------------------------------------
TT_large <- 300 #larger time length
data_H0_large <- tensorMiss::tensor_gen(K, TT_large, d, r, re, eta, coef_f, coef_fe, coef_e, seed=2025)$X
cat('Test "data_H0_large" along modes 1, 2 and 3:\n')
apply(testKron_A(ten=data_H0_large, AA=c(1,2,3), r=r, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE)$decision, 2, min)

## -----------------------------------------------------------------------------
H0_test_all <- testKron_auto(ten=data_H0, r=0, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE, all=TRUE)
H0_test_algo <- testKron_auto(ten=data_H0, r=0, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE, all=FALSE)
cat(paste0('Testing "data_H0" for all possible set of modes. The estimated rank is: ', H0_test_all$rank[1], ',', H0_test_all$rank[2], ',', H0_test_all$rank[3], '; the decision data frame is:\n'))
print.data.frame(H0_test_all$decision)
cat(paste0('\nTesting "data_H0" with the algorithm to identify each mode. The estimated rank is: ', H0_test_algo$rank[1], ',', H0_test_algo$rank[2], ',', H0_test_algo$rank[3], '; the decision data frame is:\n'))
print.data.frame(H0_test_algo$decision)

## -----------------------------------------------------------------------------
H0_test_all_large <- testKron_auto(ten=data_H0_large, r=0, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE, all=TRUE)
H0_test_algo_large <- testKron_auto(ten=data_H0_large, r=0, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE, all=FALSE)
cat(paste0('Testing "data_H0_large" for all possible set of modes. The estimated rank is: ', H0_test_all_large$rank[1], ',', H0_test_all_large$rank[2], ',', H0_test_all_large$rank[3], '; the decision data frame is:\n'))
print.data.frame(H0_test_all_large$decision)
cat(paste0('\nTesting "data_H0_large" with the algorithm to identify each mode. The estimated rank is: ', H0_test_algo_large$rank[1], ',', H0_test_algo_large$rank[2], ',', H0_test_algo_large$rank[3], '; the decision data frame is:\n'))
print.data.frame(H0_test_algo_large$decision)

## -----------------------------------------------------------------------------
noise_H1 <- tensorMiss::tensor_gen(K, TT, d, r, re, eta, coef_f, coef_fe, coef_e, seed=2025)$X - tensorMiss::tensor_gen(K, TT, d, r, re, eta, coef_f, coef_fe, coef_e, seed=2025)$C #noise series
reshape_H1 <- tensorMiss::tensor_gen(K-1, TT, c(d[-c(2,3)],prod(d[c(2,3)])), c(r[-c(2,3)],prod(r[c(2,3)])), c(re[-c(2,3)],prod(re[c(2,3)])), eta[1:(K-1)], coef_f, coef_fe, coef_e, seed=2025) #common component for the reshaped data
data_H1 <- ten_reshape_back(reshape_H1$C , c(2,3), c(TT, d)) + noise_H1

## -----------------------------------------------------------------------------
H1_test_all <- testKron_auto(ten=data_H1, r=0, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE, all=TRUE)
H1_test_algo <- testKron_auto(ten=data_H1, r=0, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE, all=FALSE)
cat(paste0('Testing "data_H1" for all possible set of modes. The estimated rank is: ', H1_test_all$rank[1], ',', H1_test_all$rank[2], ',', H1_test_all$rank[3], '; the decision data frame is:\n'))
print.data.frame(H1_test_all$decision)
cat(paste0('\nTesting "data_H1" with the algorithm to identify each mode. The estimated rank is: ', H1_test_algo$rank[1], ',', H1_test_algo$rank[2], ',', H1_test_algo$rank[3], '; the decision data frame is:\n'))
print.data.frame(H1_test_algo$decision)

## -----------------------------------------------------------------------------
noise_H1_large <- tensorMiss::tensor_gen(K, TT_large, d, r, re, eta, coef_f, coef_fe, coef_e, seed=2025)$X - tensorMiss::tensor_gen(K, TT_large, d, r, re, eta, coef_f, coef_fe, coef_e, seed=2025)$C #noise series
reshape_H1_large <- tensorMiss::tensor_gen(K-1, TT_large, c(d[-c(2,3)],prod(d[c(2,3)])), c(r[-c(2,3)],prod(r[c(2,3)])), c(re[-c(2,3)],prod(re[c(2,3)])), eta[1:(K-1)], coef_f, coef_fe, coef_e, seed=2025) #common component for the reshaped data
data_H1_large <- ten_reshape_back(reshape_H1_large$C , c(2,3), c(TT_large, d)) + noise_H1_large
H1_test_algo_large <- testKron_auto(ten=data_H1_large, r=0, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE, all=FALSE)
cat(paste0('Testing "data_H1_large" with the algorithm to identify each mode. The estimated rank is: ', H1_test_algo_large$rank[1], ',', H1_test_algo_large$rank[2], ',', H1_test_algo_large$rank[3], '; the decision data frame is:\n'))
print.data.frame(H1_test_algo_large$decision)