library(tensorPredictors)

# Load as 3D predictors `X` and flat response `y` and `F = y` with per person dim. 1 x 1
c(X, F, y) %<-% local({
    # Load from file
    ds <- readRDS("eeg_data.rds")

    # Dimension values
    n <- nrow(ds)       # sample size (nr. of people)
    p <- 64L            # nr. of predictors (count of sensorce)
    t <- 256L           # nr. of time points (measurements)

    # Extract dimension names
    nNames <- ds$PersonID
    tNames <- as.character(seq(t))
    pNames <- unlist(strsplit(colnames(ds)[2 + t * seq(p)], "_"))[c(TRUE, FALSE)]

    # Split into predictors (with proper dims and names) and response
    X <- array(as.matrix(ds[, -(1:2)]),
        dim = c(person = n, time = t, sensor = p),
        dimnames = list(person = nNames, time = tNames, sensor = pNames)
    )
    y <- ds$Case_Control

    list(X, array(y, c(n, 1L, 1L)), y)
})

# fit a tensor normal model to the data sample axis 1 indexes persons)
fit.gmlm <- gmlm_tensor_normal(X, F, sample.axis = 1L)

# plot the fitted mode wise reductions (for time and sensor axis)
with(fit.gmlm, {
    par.reset <- par(mfrow = c(2, 1))
    plot(seq(0, 1, len = 256), betas[[1]], main = "Time", xlab = "Time [s]", ylab = expression(beta[1]))
    plot(betas[[2]], main = "Sensors", xlab = "Sensor Index", ylab = expression(beta[2]))
    par(par.reset)
})


#' (2D)^2 PCA preprocessing
#'
#' @param tpc Number of "t"ime "p"rincipal "c"omponents.
#' @param ppc Number of "p"redictor "p"rincipal "c"omponents.
preprocess <- function(X, tpc, ppc) {
    # Mode covariances (for predictor and time point modes)
    c(Sigma_t, Sigma_p) %<-% mcov(X, sample.axis = 1L)

    # "predictor" (sensor) and time point principal components
    V_t <- svd(Sigma_t, tpc, 0L)$u
    V_p <- svd(Sigma_p, ppc, 0L)$u

    # reduce with mode wise PCs
    mlm(X, list(V_t, V_p), modes = 2:3, transposed = TRUE)
}


#' Leave-one-out prediction
#'
#' @param X 3D EEG data (preprocessed or not)
#' @param F binary responce `y` as a 3D tensor, every obs. is a 1 x 1 matrix
loo.predict <- function(X, F) {
    sapply(seq_len(dim(X)[1L]), function(i) {
        # Fit with i'th observation removes
        fit <- gmlm_tensor_normal(X[-i, , ], F[-i, , , drop = FALSE], sample.axis = 1L)

        # Reduce the entire data set
        r <- as.vector(mlm(X, fit$betas, modes = 2:3, transpose = TRUE))
        # Fit a logit model on reduced data with i'th observation removed
        logit <- glm(y ~ r, family = binomial(link = "logit"),
            data = data.frame(y = y[-i], r = r[-i])
        )
        # predict i'th response given i'th reduced observation
        y.hat <- predict(logit, newdata = data.frame(r = r[i]), type = "response")
        # report progress
        cat(sprintf("dim: (%d, %d) - %3d/%d\n", dim(X)[2L], dim(X)[3L], i, dim(X)[1L]))

        y.hat
    })
}


### Classification performance measures
# acc: Accuracy. P(Yhat = Y). Estimated as: (TP+TN)/(P+N).
acc <- function(y.true, y.pred) mean(round(y.pred) == y.true)
# err: Error rate. P(Yhat != Y). Estimated as: (FP+FN)/(P+N).
err <- function(y.true, y.pred) mean(round(y.pred) != y.true)
# fpr: False positive rate. P(Yhat = + | Y = -). aliases: Fallout.
fpr <- function(y.true, y.pred) mean((round(y.pred) == 1)[y.true == 0])
# tpr: True positive rate.  P(Yhat = + | Y = +). aliases: Sensitivity, Recall.
tpr <- function(y.true, y.pred) mean((round(y.pred) == 1)[y.true == 1])
# fnr: False negative rate. P(Yhat = - | Y = +). aliases: Miss.
fnr <- function(y.true, y.pred) mean((round(y.pred) == 0)[y.true == 1])
# tnr: True negative rate.  P(Yhat = - | Y = -).
tnr <- function(y.true, y.pred) mean((round(y.pred) == 0)[y.true == 0])
# auc: Area Under the Curve
auc <- function(y.true, y.pred) as.numeric(pROC::roc(y.true, y.pred, quiet = TRUE)$auc)
auc.sd <- function(y.true, y.pred) sqrt(pROC::var(pROC::roc(y.true, y.pred, quiet = TRUE)))


# perform preprocessed (reduced) and raw (not reduced) leave-one-out prediction
y.hat.3.4   <- loo.predict(preprocess(X,  3,  4), F)
y.hat.15.15 <- loo.predict(preprocess(X, 15, 15), F)
y.hat.20.30 <- loo.predict(preprocess(X, 20, 30), F)
y.hat       <- loo.predict(X, F)

# classification performance measures table by leave-one-out cross-validation
(loo.cv <- apply(cbind(y.hat.3.4, y.hat.15.15, y.hat.20.30, y.hat), 2, function(y.pred) {
    sapply(c("acc", "err", "fpr", "tpr", "fnr", "tnr", "auc", "auc.sd"),
        function(FUN) { match.fun(FUN)(y, y.pred) })
}))
#>         y.hat.3.4 y.hat.15.15 y.hat.20.30      y.hat
#> acc    0.79508197  0.78688525  0.78688525 0.78688525
#> err    0.20491803  0.21311475  0.21311475 0.21311475
#> fpr    0.35555556  0.40000000  0.40000000 0.40000000
#> tpr    0.88311688  0.89610390  0.89610390 0.89610390
#> fnr    0.11688312  0.10389610  0.10389610 0.10389610
#> tnr    0.64444444  0.60000000  0.60000000 0.60000000
#> auc    0.85108225  0.83838384  0.83924964 0.83896104
#> auc.sd 0.03584791  0.03760531  0.03751307 0.03754553