diff --git a/NAMESPACE b/NAMESPACE
index c913ff05..091e1d20 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -108,14 +108,12 @@ S3method(summary,mixo_pls)
 S3method(summary,mixo_spls)
 S3method(summary,pca)
 S3method(summary,rcc)
-export("%fp%")
 export(auroc)
 export(background.predict)
 export(block.pls)
 export(block.plsda)
 export(block.spls)
 export(block.splsda)
-export(block_tpls)
 export(cim)
 export(cimDiablo)
 export(circosPlot)
@@ -123,22 +121,15 @@ export(color.GreenRed)
 export(color.jet)
 export(color.mixo)
 export(color.spectral)
-export(dctii_m_transforms)
 export(explained_variance)
-export(facewise_crossproduct)
-export(facewise_product)
-export(facewise_transpose)
-export(ft)
 export(get.BER)
 export(get.confusion_matrix)
 export(imgCor)
 export(impute.nipals)
 export(ipca)
 export(logratio.transfo)
-export(m_product)
 export(map)
 export(mat.rank)
-export(matrix_to_m_transforms)
 export(mint.block.pls)
 export(mint.block.plsda)
 export(mint.block.spls)
@@ -169,10 +160,6 @@ export(spca)
 export(spls)
 export(splsda)
 export(study_split)
-export(tpca)
-export(tpls)
-export(tplsda)
-export(tsvdm)
 export(tune)
 export(tune.block.splsda)
 export(tune.mint.splsda)
diff --git a/R/tens.block.tpls.R b/R/tens.block.tpls.R
deleted file mode 100644
index 45bfd50e..00000000
--- a/R/tens.block.tpls.R
+++ /dev/null
@@ -1,48 +0,0 @@
-# ==============================================================================
-# Tensor block pls generalization; developed @ Melbourne Integrative Genomics
-# based on Kilmer's tensor m-product algebra
-# ==============================================================================
-
-#' Tensor block PLS
-#'
-#' Developed @ Melbourne Integrative Genomics.
-#'
-#' @param x A list of tensor inputs (termed 'blocks') measured on the same
-#' samples.
-#' @param y Tensor input.
-#' @param ncomp The estimated number of components. ncomp must be explicitly set
-#' as an integer in tpls.
-#' @param design Numeric matrix of size length(x) x length(x) with values
-#' between 0 and 1; the i,j entry in the matrix indicates the strength of the
-#' relationship to be modelled between the i-th and j-th blocks. A value of 0
-#' indicates no relationship, with 1 being the maximum. Alternatively, one can
-#' input "null" for a fully disconnected design, or "full" for a fully connected
-#' design, or a single scalar value between 0 and which will designate all
-#' off-diagonal elements of a fully connected design.
-#' @param m A function which applies an orthogonal tensor tubal transform.
-#' @param minv The inverse of m.
-#' @param mode Currently supports tensor analogues of canonical, regression,
-#' and svd PLS modes. Defaults to "regression" mode.
-#' @param center If set to false, the data tensor will not be centralized into
-#' Mean Deviation Form (see Mor et al. 2022). By default, the mean horizontal
-#' slice of the input tensor(s) are subtracted, so that all of the horizontal
-#' slices sum to 0, analgous to centering matrix data.
-#' @param bpparam A \linkS4class{BiocParallelParam} object indicating the type
-#' of parallelisation. Does not have any effect if transform functions
-#' explicitly set using \code{m}, \code{minv}.
-#' @author Brendan Lu
-#' @export
-block_tpls <- function(
-  x,
-  y,
-  ncomp = NULL,
-  design,
-  m = NULL,
-  minv = NULL,
-  mode = "regression",
-  center = TRUE,
-  bpparam = NULL
-) {
-  # bltodo: add docs link to existing block.pls for user reference
-  NULL
-}
diff --git a/R/tens.mproduct.R b/R/tens.mproduct.R
deleted file mode 100644
index a68055e6..00000000
--- a/R/tens.mproduct.R
+++ /dev/null
@@ -1,321 +0,0 @@
-# ==============================================================================
-# utilities and functions for Kilmer's m product
-# ==============================================================================
-# bltodo: use classes for m and minv?
-
-#' Apply a function across the last dimension of an input vector,
-#' matrix, or tensor. This function defines both a parallel algorithm using
-#' BiocParrallel, and also a simple \code{apply} algorithm. Even on Windows
-#' Machines, setting \code{bpparam = BiocParallel::SerialParam()} offers a
-#' notable speedup for larger 3D array inputs.
-#'
-#' @param x Numerical array input.
-#' @param mat Function which defines the tubal transform.
-#' @param bpparam A \linkS4class{BiocParallelParam} object indicating the type
-#' of parallelisation.
-#' @return A tensor of the same size under the specified tubal transform,
-#' denoted \eqn{\hat{x}}.
-#' @author Brendan Lu
-#' @keywords internal
-#' @noRd
-.apply_mat_transform <- function(x, mat, bpparam) {
-  if (length(dim(x)) == 1) {
-    return(mat %*% x)
-  } else if (length(dim(x)) == 2) {
-    return(mat %*% x)
-  } else if (length(dim(x)) == 3) {
-    n <- dim(x)[1]
-    p <- dim(x)[2]
-    t <- dim(x)[3]
-    if (is.null(bpparam)) {
-      # apply algorithm --------------------------------------------------------
-      transformed_slices <- apply(
-        aperm(x, c(3, 1, 2)), c(2, 3), # permuted tensor with mode 3 in front
-        function(slice) mat %*% slice # left multiply by transform mat
-      )
-      # permute back to original orientation
-      return(aperm(array(transformed_slices, dim = c(t, n, p)), c(2, 3, 1)))
-      # ------------------------------------------------------------------------
-    } else {
-      # BiocParallel algorithm -------------------------------------------------
-      transformed_slices <- BiocParallel::bplapply(
-        lapply(seq_len(p), function(i) t(x[, i, ])), # a list of t x n matrices
-        function(slice) mat %*% slice, # left multiply by transform mat
-        BPPARAM = bpparam
-      )
-      # unlist and cast into array with p facewise t x n matrices
-      # then appropriately rotate to get original n x p x t matrix
-      return(
-        aperm(array(unlist(transformed_slices), dim = c(t, n, p)), c(2, 3, 1))
-      )
-      # ------------------------------------------------------------------------
-    }
-  } else {
-    # error: some array >3D has been inputted
-    stop(
-      "Only order 1 (vector), 2 (matrix), 3 (3D tensor) arrays are supported"
-    )
-  }
-}
-
-#' Validate appropriate 'null-ness' of m, minv inputs, and apply default
-#' transform as needed.
-#'
-#' @author Brendan Lu
-#' @keywords internal
-#' @noRd
-.stop_invalid_transform_input <- function(m, minv, t, bpparam) {
-  if (
-    xor(is.function(m), is.function(minv)) ||
-      xor(is.null(m), is.null(minv))
-  ) {
-    stop(
-      "If explicitly defined, both m and its inverse must be defined as 
-      functions"
-    )
-  }
-
-  # due to above checks, the code below executes when both m and minv are NULL
-  # and returns the calling state the default transform used in this library
-  # which is currently the dct-ii transform
-  if (is.null(m)) {
-    transforms <- dctii_m_transforms(t, bpparam = bpparam)
-    m <- transforms$m
-    minv <- transforms$minv
-  }
-
-  return(list(
-    m = m,
-    minv = minv
-  ))
-}
-
-#' Create m-transform functions as defined by a matrix
-#'
-#' Returns functions \code{m} and \code{m_inv} which apply tubal transforms
-#' defined by the matrix \code{m_mat}.
-#'
-#' @param m_mat Function which defines the tubal transform.
-#' @param m_inv_mat Function which defines inverse tubal transform
-#' @param bpparam A \linkS4class{BiocParallelParam} object indicating the type
-#' of parallelisation.
-#' @return
-#' \item{m}{A function which applies the matrix m_mat along the last dimension
-#' of a given numerical input array. For 3D tensors it performs the mode-3
-#' product.}
-#' \item{minv}{The inverse of m.}
-#' @author Brendan Lu
-#' @export
-matrix_to_m_transforms <- function(
-  m_mat,
-  m_inv_mat = NULL,
-  bpparam = NULL
-) {
-  # error: non-matrix input
-  stopifnot(length(dim(m_mat)) == 2)
-  # error: non-square matrix input
-  stopifnot(dim(m_mat)[1] == dim(m_mat)[2])
-
-  # invert m_mat if m_inv_mat not specified
-  if (is.null(m_inv_mat)) {
-    # error: non-invertible input
-    m_inv_mat <- solve(m_mat)
-  } else {
-    # error: specified matrices are not the same size
-    stopifnot(identical(dim(m_mat), dim(m_inv_mat)))
-  }
-
-  return(list(
-    m = function(x) .apply_mat_transform(x, m_mat, bpparam = bpparam),
-    minv = function(x) .apply_mat_transform(x, m_inv_mat, bpparam = bpparam)
-  ))
-}
-
-#' Create m-transform functions to apply the Discrete Cosine Transform 2
-#'
-#' Returns functions \code{m} and \code{m_inv} which apply tubal transforms
-#' defined by the Discrete Cosine Transform (DCT-II variant). This is equivalent
-#' to Scipys DCTI-ii algorithm with \code{norm='ortho'}.
-#'
-#' @param t The length of the transform.
-#' @param bpparam A \linkS4class{BiocParallelParam} object indicating the type
-#' of parallelisation.
-#' @return
-#' \item{m}{A function which applies the dct-ii along the last dimension of a
-#' given numerical input array. For 3D tensors it performs the mode-3 product
-#' with the DCT matrix.}
-#' \item{minv}{The inverse of m,}
-#' @author Brendan Lu
-#' @export
-dctii_m_transforms <- function(t, bpparam = NULL) {
-  return(matrix_to_m_transforms(m_mat = gsignal::dctmtx(t), bpparam = bpparam))
-}
-
-#' Compute Kilmer's facewise product. Note that the for-loop implementation is
-#' relatively fast, and very readable. There's also a BiocParralel
-#' implementation here, but it lacks significant benchmarking results.
-#' bltodo: the parallel algorithm is probably stupid remove sometime?
-#'
-#' @author Brendan Lu
-#' @keywords internal
-#' @noRd
-.binary_facewise <- function(a, b, bpparam) {
-  na <- dim(a)[1]
-  pa <- dim(a)[2]
-  ta <- dim(a)[3]
-
-  nb <- dim(b)[1]
-  pb <- dim(b)[2]
-  tb <- dim(b)[3]
-
-  # error: different t for each input
-  stopifnot(ta == tb)
-  t <- ta
-  # error: non-conforming facewise dimensions
-  stopifnot(pa == nb)
-
-  if (is.null(bpparam)) {
-    # for-loop algorithm -------------------------------------------------------
-    fp_ab <- array(0, dim = c(na, pb, t))
-    for (i in seq_len(t)) {
-      fp_ab[, , i] <- a[, , i] %*% b[, , i]
-    }
-    return(fp_ab)
-    # --------------------------------------------------------------------------
-  } else {
-    # BiocParallel algorithm ---------------------------------------------------
-    # bltodo: benchmark / investigate preallocation here?
-    return(
-      simplify2array(
-        BiocParallel::bplapply(
-          array(seq_len(t)),
-          FUN = function(i) a[, , i] %*% b[, , i],
-          BPPARAM = bpparam
-        )
-      )
-    )
-    # --------------------------------------------------------------------------
-  }
-}
-
-#' Tensor facewise product
-#'
-#' Compute Kilmer's facewise product cumulatively across any arbitrary number of
-#' tensor inputs.
-#'
-#' @param ... Arbitrary number of numerical tensor inputs.
-#' @param bpparam A \linkS4class{BiocParallelParam} object indicating the type
-#' of parallelisation.
-#' @return Cumulative facewise product.
-#' @author Brendan Lu
-#' @export
-facewise_product <- function(..., bpparam = NULL) {
-  return(
-    Reduce(
-      function(a, b) .binary_facewise(a, b, bpparam = bpparam),
-      list(...)
-    )
-  )
-}
-
-#' @describeIn facewise_product Custom facewise product operator
-#' @param a Tensor input.
-#' @param b Tensor input.
-#' @export
-`%fp%` <- function(a, b) .binary_facewise(a, b, bpparam = NULL)
-
-#' Facewise transpose an order-3 tensor
-#'
-#' @param tensor Numerical 3D array input.
-#' @return Facewise transpose of \code{tensor}
-#' @author Brendan Lu
-#' @aliases ft facewise_transpose
-#' @export
-facewise_transpose <- function(tensor) {
-  return(aperm(tensor, c(2, 1, 3)))
-}
-
-#' @describeIn facewise_transpose Alias for \code{\link{facewise_product}}
-#' @export
-ft <- function(tensor) facewise_transpose(tensor)
-
-#' Kilmer's tensor-tensor m-product
-#'
-#' This function works cumulatively across any arbitrary number of tensor
-#' inputs.
-#'
-#' @param ... Arbitrary number of numerical tensor inputs.
-#' @param m A function which applies an orthogonal tensor tubal transform.
-#' @param minv The inverse of m.
-#' @param bpparam A \linkS4class{BiocParallelParam} object indicating the type
-#' of parallelisation. Does not have any effect if transform functions
-#' explicitly set using \code{m}, \code{minv}.
-#' @return Cumulative m-product.
-#' @author Brendan Lu
-#' @export
-m_product <- function(
-  ...,
-  m = NULL,
-  minv = NULL,
-  bpparam = NULL
-) {
-  tensors <- list(...)
-  if (length(tensors) == 0) {
-    stop("No input tensors provided.")
-  } else {
-    t <- dim(tensors[[1]])[3]
-  }
-  if (
-    xor(is.function(m), is.function(minv)) ||
-      xor(is.null(m), is.null(minv))
-  ) {
-    stop(
-      "If explicitly defined, both m and its inverse must be defined as 
-      functions."
-    )
-  }
-  # use dctii as default transform if user does not specify an explicit one
-  if (is.null(m)) {
-    transforms <- dctii_m_transforms(t, bpparam = bpparam)
-    m <- transforms$m
-    minv <- transforms$minv
-  }
-  return(minv(
-    Reduce(
-      # bpparam MUST be NULL here to prevent double parallelisation!
-      function(a, b) .binary_facewise(a, b, bpparam = NULL),
-      lapply(list(...), m)
-    )
-  ))
-}
-
-#' Tensor cross product
-#'
-#' Compute the equivalent of ft(a) %fp% b
-#'
-#' @param a Tensor input.
-#' @param b Tensor input.
-#' @return Tensor facewise cross product.
-#' @author Brendan Lu
-#' @export
-facewise_crossproduct <- function(a, b) {
-  na <- dim(a)[1]
-  pa <- dim(a)[2]
-  ta <- dim(a)[3]
-
-  nb <- dim(b)[1]
-  pb <- dim(b)[2]
-  tb <- dim(b)[3]
-
-  # error: different t for each input
-  stopifnot(ta == tb)
-  t <- ta
-  # error: non-conforming facewise dimensions
-  stopifnot(na == nb)
-
-  fcp_ab <- array(0, dim = c(pa, pb, t))
-  for (i in seq_len(t)) {
-    fcp_ab[, , i] <- crossprod(a[, , i], b[, , i])
-  }
-  return(fcp_ab)
-}
diff --git a/R/tens.tpca.R b/R/tens.tpca.R
deleted file mode 100644
index 7ef9e249..00000000
--- a/R/tens.tpca.R
+++ /dev/null
@@ -1,192 +0,0 @@
-# ==============================================================================
-# Tensor pca based on Mor's TCAM algorithm
-# ==============================================================================
-
-#' R implementation of np.unravel_index. NOTE: currently only works for 1D to 2D
-#' column-major conversion, and returns a list of 2D indices. Returns a matrix
-#' output of length(indices) columns, with two rows. The first row corresponds
-#' to the sorted k indices, and the second row contains the sorted t indices.
-#'
-#' @author Brendan Lu
-#' @keywords internal
-#' @noRd
-.unravel_index <- function(indices, dim) {
-  nrows <- dim[1]
-  return(sapply(
-    indices,
-    FUN = function(x) {
-      c(
-        (x - 1) %% nrows + 1, # transformed k tensor position
-        (x - 1) %/% nrows + 1 # transformed t tensor position
-      )
-    }
-  ))
-}
-
-#' Extract tensor columns specified by .unravel_index output.
-#' Effectively achieves:
-#'
-#'     self.loadings_matrix_ = hatV[
-#'         :, self._k_t_flatten_sort[0], self._k_t_flatten_sort[1]
-#'     ].T
-#'
-#' type of indexing in Numpy.
-#'
-#' This is the 'tensor compression' based on the ordered indices specified
-#' by each column of k_t_indices as described in Mor's TCAM.
-#'
-#' @author Brendan Lu
-#' @keywords internal
-#' @noRd
-.extract_tensor_columns <- function(tensor, k_t_indices) {
-  return(apply(
-    k_t_indices, 2,
-    FUN = function(index_column) tensor[, index_column[1], index_column[2]]
-  ))
-}
-
-#' Helper function to convert compressed matrix form of the singular values into
-#' a sparse tensor.
-#'
-#' @param mat Matrix s.t. each column contains the f-diagonal singular values
-#' @param dim Dimension of output tensor
-#' @author Brendan Lu
-#' @keywords internal
-#' @noRd
-.singular_vals_mat_to_tens <- function(mat, dim) {
-  n <- dim[1]
-  p <- dim[2]
-  t <- dim[3]
-  k <- min(n, p)
-
-  tens <- array(0, dim = dim)
-  for (i in seq_len(t)) {
-    tens[1:k, 1:k, i] <- diag(mat[, i])
-  }
-  return(tens)
-}
-
-#' Tensor PCA-like dimensionality reduction
-#'
-#' Tensor analogue of PCA introduced by Mor et al. (2022) based on Kilmer's
-#' m-product algebra and tsvdm.
-#'
-#' @param x Tensor input.
-#' @param ncomp The estimated number of components. ncomp can be explicitly set
-#' using an integer value, 0 < float value < 1, or left as NULL which will
-#' default to the maximum ncomp value possible.
-#' @param m A function which applies an orthogonal tensor tubal transform.
-#' @param minv The inverse of m.
-#' @param center If set to false, the data tensor will not be centralized into
-#' Mean Deviation Form (see Mor et al. 2022). By default, the mean horizontal
-#' slice of the input tensor(s) are subtracted, so that all of the horizontal
-#' slices sum to 0, analgous to centering matrix data.
-#' @param matrix_output Collect the top singular vectors across the tensor,
-#' organized by the magnitude of the corresponding singular vector, and place
-#' them into the columns of a matrix. Corresponds to the 'matrix compression'
-#' type of scheme described in Mor et al. 2022.
-#' @param bpparam A \linkS4class{BiocParallelParam} object indicating the type
-#' of parallelisation. Does not have any effect if transform functions
-#' explicitly set using \code{m}, \code{minv}.
-#' @author Brendan Lu
-#' @export
-tpca <- function(
-  x,
-  ncomp = NULL,
-  m = NULL,
-  minv = NULL,
-  center = TRUE,
-  matrix_output = TRUE,
-  bpparam = NULL
-) {
-  if (length(dim(x)) != 3) {
-    stop("Please ensure input tensor is an order-3 array.")
-  } else {
-    n <- dim(x)[1]
-    p <- dim(x)[2]
-    t <- dim(x)[3]
-    k <- min(n, p)
-  }
-
-  # use dctii as default transform if user does not specify an explicit one
-  validated_transforms <- .stop_invalid_transform_input(m, minv, t, bpparam)
-  m <- validated_transforms$m
-  minv <- validated_transforms$minv
-
-  # bltodo: add scaling as well?
-  if (center) {
-    mean_slice <- apply(x, c(2, 3), mean)
-    x <- sweep(x, c(2, 3), STATS = mean_slice, FUN = "-")
-  }
-
-  # NOTE: passing in bpparam configuration is not needed, as any
-  # bpparam specification will already take effect above in the transforms
-  tsvdm_decomposition <- tsvdm(
-    x, m, minv,
-    keep_hats = TRUE,
-    svals_matrix_form = TRUE
-  )
-
-  # flatten out Fortran column major style, then get sort order
-  singular_values <- as.vector(tsvdm_decomposition$shat)
-  k_t_flatten_sort <- order(singular_values, decreasing = TRUE)
-  singular_values <- singular_values[k_t_flatten_sort]
-
-  # get explained variance
-  # bltodo: check relevant tensor theory for this!
-  squared_singular_values <- singular_values ^ 2
-  total_var <- sum(squared_singular_values)
-  explained_variance_ratio <- squared_singular_values / total_var
-
-  # process ncomp input
-  if (is.null(ncomp)) {
-    ncomp <- k * t
-  } else if (ncomp > 0 && ncomp < 1) {
-    # pick ncomp to be minimum number of integer components to explain the
-    # inputted variance
-    ratio_cumsum <- cumsum(explained_variance_ratio)
-    ncomp <- findInterval(ncomp, ratio_cumsum) + 1
-  } else if (ncomp %% 1 == 0) {
-    stopifnot(ncomp > 0 && ncomp <= (k * t))
-  } else {
-    stop("Please input an integer, 0 < float < 1, or NULL for ncomp parameter")
-  }
-
-  # convert the argsort indexes back into the two dimensional indexes
-  # corresponding to the singular values matrix
-  # NOTE: these are the collection if i_h's and j_h's in Mor et al. (2022)
-  k_t_flatten_sort <- .unravel_index(
-    k_t_flatten_sort[1:ncomp],
-    dim(tsvdm_decomposition$shat)
-  )
-
-  # compute the values of the projected data
-  # bltodo: benchmark this against \eqn{new_data *_M vhat}
-  # bltodo: zero out uhat first? svd truncates it to rank already, but can
-  # further truncate to ncomp columns?
-  x_projected <- tsvdm_decomposition$uhat %fp%
-    .singular_vals_mat_to_tens(tsvdm_decomposition$shat, dim = c(n, p, t))
-
-  loadings <- tsvdm_decomposition$vhat
-
-  # extract the columns in compressed form
-  if (matrix_output) {
-    x_projected <- .extract_tensor_columns(x_projected, k_t_flatten_sort)
-    loadings <- .extract_tensor_columns(
-      tsvdm_decomposition$vhat,
-      k_t_flatten_sort
-    )
-  }
-
-  # BLTODO: (to add in) zero out and compute rho as well
-  # see _rank_q_truncation_zero_out() in `tred`
-  # need to use this for tensor output to be appropriate
-
-  return(invisible(list(
-    ncomp = ncomp,
-    x = x,
-    loadings = loadings,
-    variates = x_projected, # bltodo: maybe just adopt a different name here?
-    explained_variance = explained_variance_ratio[1:ncomp]
-  )))
-}
diff --git a/R/tens.tpls.R b/R/tens.tpls.R
deleted file mode 100644
index 16cb771e..00000000
--- a/R/tens.tpls.R
+++ /dev/null
@@ -1,284 +0,0 @@
-# ==============================================================================
-# Tensor pls generalization; developed @ Melbourne Integrative Genomics based
-# on Kilmer's tensor m-product algebra
-# ==============================================================================
-
-#' Convert a singular values tensor (in compressed matrix form) to a set of
-#' indices corresponding to the (column,face) pairs of the top `ncomp` singular
-#' values. NEEDS singular values to be in matrix form.
-#'
-#' @author Brendan Lu
-#' @keywords internal
-#' @noRd
-.obtain_k_t_flatten_sort <- function(s_mat, ncomp) {
-  # bltodo: if ultimately only use once just place it in function body directly
-  return(
-    .unravel_index(
-      order(as.vector(s_mat))[1:ncomp],
-      dim(s_mat)
-    )
-  )
-}
-
-#' Get the k, t index corresponding to the largest singular value. Mirrors
-#' .unravel_index() but more efficient as we only care about the largest
-#' singular value in tpls.
-#'
-#' @author Brendan Lu
-#' @keywords internal
-#' @noRd
-.obtain_k_t_top <- function(s_mat) {
-  flat_index <- which.max(as.vector(s_mat))
-  nrows <- dim(s_mat)[1]
-  return(c(
-    (flat_index - 1) %% nrows + 1, # transformed k tensor position
-    (flat_index - 1) %/% nrows + 1 # transformed t tensor position
-  ))
-}
-
-#' Run some sense checks on x and y tensor inputs for tpls
-#'
-#' @author Brendan Lu
-#' @keywords internal
-#' @noRd
-.validate_tpls_x_y_dim <- function(x, y) {
-  if (length(dim(x)) != 3) {
-    stop("Please ensure x input tensor is an order-3 array")
-  } else {
-    n <- dim(x)[1]
-    p <- dim(x)[2]
-    t <- dim(x)[3]
-  }
-
-  if (length(dim(y)) != 3) {
-    stop("Please ensure y input tensor is an order-3 array")
-  } else {
-    n2 <- dim(y)[1]
-    q <- dim(y)[2]
-    t2 <- dim(y)[3]
-  }
-
-  if (n != n2) {
-    stop("Please ensure x and y tensor inputs have matching number of samples")
-  }
-
-  if (t != t2) {
-    stop("Please ensure x and y tensor inputs have matching number of time 
-    points")
-  }
-
-  return(list(
-    n = n,
-    p = p,
-    q = q,
-    t = t
-  ))
-}
-
-#' Tensor PLS-like regression
-#'
-#' Developed @ Melbourne Integrative Genomics. Intuition: Fourier-based m
-#' transforms compact the data across time points, but preserve the information
-#' association for each sample and each feature. In doing so, pls works by
-#' calculating the singular vectors associated with a 'global' (across time
-#' points) top singular value, and deflating just that relevant face.
-#'
-#' @param x Tensor input.
-#' @param y Tensor input.
-#' @param ncomp The estimated number of components. ncomp must be explicitly set
-#' as an integer in tpls.
-#' @param m A function which applies an orthogonal tensor tubal transform.
-#' @param minv The inverse of m.
-#' @param mode Currently supports tensor analogues of canonical ("canonical"),
-#' regression ("regression"), and svd ("tsvdm") PLS variants. Defaults to
-#' "regression".
-#' @param center If set to false, the data tensor will not be centralized into
-#' Mean Deviation Form (see Mor et al. 2022). By default, the mean horizontal
-#' slice of the input tensor(s) are subtracted, so that all of the horizontal
-#' slices sum to 0, analgous to centering matrix data.
-#' @param matrix_output Note: ONLY AFFECTS THE OUTPUT IN "tsvdm" MODE.
-#' Collects the top singular vectors across the tensor, organized by the
-#' magnitude of the corresponding singular vector, and place them into the
-#' columns of a matrix. Corresponds to the 'matrix compression' type of scheme
-#' described in Mor et al. 2022.
-#' @param bpparam A \linkS4class{BiocParallelParam} object indicating the type
-#' of parallelisation. Does not have any effect if transform functions
-#' explicitly set using \code{m}, \code{minv}.
-#' @author Brendan Lu
-#' @export
-tpls <- function(
-  x,
-  y,
-  ncomp = NULL,
-  m = NULL,
-  minv = NULL,
-  mode = "regression",
-  center = TRUE,
-  matrix_output = TRUE, # FALSE only takes effect for tsvdm method
-  bpparam = NULL
-) {
-  # allowed modes mirror sklearn's 2D PLS models
-  # see https://scikit-learn.org/stable/modules/cross_decomposition.html
-  allowed_modes <- c("canonical", "regression", "tsvdm")
-  if (!(mode %in% allowed_modes)) {
-    stop(paste(
-      "Please ensure mode is one of: ",
-      paste(allowed_modes, collapse = ", ")
-    ))
-  }
-
-  dims <- .validate_tpls_x_y_dim(x, y)
-  n <- dims$n
-  p <- dims$p
-  q <- dims$q
-  t <- dims$t
-  k <- min(n, p, q)
-  # maximum number of non-zero entries in the f-diagonal singular values tensor
-  # of tsvdm(XtY) is k * t
-  max_rank <- k * t
-
-  # use dctii as default transform if user does not specify an explicit one
-  validated_transforms <- .stop_invalid_transform_input(m, minv, t, bpparam)
-  m <- validated_transforms$m
-  minv <- validated_transforms$minv
-
-  if (center) {
-    mean_slice_x <- apply(x, c(2, 3), mean)
-    mean_slice_y <- apply(y, c(2, 3), mean)
-    x <- sweep(x, c(2, 3), STATS = mean_slice_x, FUN = "-")
-    y <- sweep(y, c(2, 3), STATS = mean_slice_y, FUN = "-")
-  }
-
-  # project x and y into 'hat-space', expensive computation so do it once and
-  # save into variable
-  xhat <- m(x)
-  yhat <- m(y)
-
-  # process ncomp input, much simpler than tpca as we only accept integer input
-  # bltodo: investigate non integer input options? explained variance semantics?
-  if (is.null(ncomp)) {
-    ncomp <- max_rank
-  } else if (ncomp %% 1 == 0) {
-    stopifnot(ncomp > 0 && ncomp <= max_rank)
-  } else {
-    stop("Please input an integer or NULL for ncomp parameter")
-  }
-
-  if (mode == "tsvdm") {
-    # simplest algorithm - just uses everything from the tsvdm call of XtY based
-    # without any deflation steps
-    tsvdm_decomposition_xty <- tsvdm(
-      # bltodo: change to facewise_crossproduct when this is improved
-      ft(xhat) %fp% yhat,
-      transform = FALSE,
-      svals_matrix_form = TRUE
-    )
-
-    # loadings are just the u,v tensors from the tsvdm call
-    x_loadings <- tsvdm_decomposition_xty$u
-    y_loadings <- tsvdm_decomposition_xty$v
-
-    # project the scaled / transformed x, y data onto the loadings
-    x_projected <- xhat %fp% x_loadings
-    y_projected <- yhat %fp% y_loadings
-
-    if (matrix_output) {
-      # bltodo: tpls explained variance?
-      k_t_flatten_sort <- .obtain_k_t_flatten_sort(
-        tsvdm_decomposition_xty$s,
-        ncomp
-      )
-      x_loadings <- .extract_tensor_columns(x_loadings, k_t_flatten_sort)
-      y_loadings <- .extract_tensor_columns(y_loadings, k_t_flatten_sort)
-      x_projected <- .extract_tensor_columns(x_projected, k_t_flatten_sort)
-      y_projected <- .extract_tensor_columns(y_projected, k_t_flatten_sort)
-    }
-
-    return(invisible(list(
-      ncomp = ncomp,
-      x = x,
-      y = y,
-      x_loadings = x_loadings,
-      y_loadings = y_loadings,
-      x_projected = x_projected,
-      y_projected = y_projected
-    )))
-
-  } else if (mode == "canonical" || mode == "regression") {
-    # preallocate output arrays which we will fill during the iterative process
-    x_loadings <- array(0, dim = c(p, ncomp))
-    y_loadings <- array(0, dim = c(q, ncomp))
-    x_projected <- array(0, dim = c(n, ncomp))
-    y_projected <- array(0, dim = c(n, ncomp))
-
-    for (i in seq_len(ncomp)) {
-      # compute tsvdm
-      # BLTODO: tensor crossprod to speed up?
-      # bltodo: minor as ncomp is usually small, but we technically do not need
-      # to recalculate the whole tsvdmm, we only need svd re-done on the face
-      # that was deflated as everything else is still the same
-      tsvdm_decomposition_xty <- tsvdm(
-        # bltodo: change to facewise_crossproduct when this is improved
-        ft(xhat) %fp% yhat,
-        transform = FALSE,
-        svals_matrix_form = TRUE
-      )
-
-      # get indices corresponding to largest singular value
-      k_t_top <- .obtain_k_t_top(tsvdm_decomposition_xty$s)
-      curr_x_loadings <- tsvdm_decomposition_xty$u[, k_t_top[1], k_t_top[2]]
-      curr_y_loadings <- tsvdm_decomposition_xty$v[, k_t_top[1], k_t_top[2]]
-
-      # note that only one face of xhat and yhat is relevant per iteration
-      # using k_t_top[2] we can basically just work in matrix world for a bit
-      curr_x_projected <- xhat[, , k_t_top[2]] %*% curr_x_loadings
-      curr_y_projected <- yhat[, , k_t_top[2]] %*% curr_y_loadings
-
-      # perform deflation
-      # the calculation of the x regression coefficient and deflation step of
-      # xhat is the same regardless of pls mode
-      curr_x_reg_coef <- crossprod(xhat[, , k_t_top[2]], curr_x_projected) /
-        as.numeric(crossprod(curr_x_projected, curr_x_projected))
-
-      xhat[, , k_t_top[2]] <- xhat[, , k_t_top[2]] -
-        tcrossprod(curr_x_projected, curr_x_reg_coef)
-
-      # the calculation of the y regression coefficient and deflation step of
-      # yhat differs depending on the pls mode
-      if (mode == "canonical") {
-        curr_y_reg_coef <- crossprod(yhat[, , k_t_top[2]], curr_y_projected) /
-          as.numeric(crossprod(curr_y_projected, curr_y_projected))
-
-        yhat[, , k_t_top[2]] <- yhat[, , k_t_top[2]] -
-          tcrossprod(curr_y_projected, curr_y_reg_coef)
-      }
-
-      if (mode == "regression") {
-        curr_y_reg_coef <- crossprod(yhat[, , k_t_top[2]], curr_x_projected) /
-          as.numeric(crossprod(curr_x_projected, curr_x_projected))
-
-        yhat[, , k_t_top[2]] <- yhat[, , k_t_top[2]] -
-          tcrossprod(curr_x_projected, curr_y_reg_coef)
-      }
-
-      x_loadings[, i] <- curr_x_loadings
-      y_loadings[, i] <- curr_y_loadings
-      x_projected[, i] <- curr_x_projected
-      y_projected[, i] <- curr_y_projected
-    }
-
-    return(invisible(list(
-      ncomp = ncomp,
-      x = x,
-      y = y,
-      mode = mode,
-      x_loadings = x_loadings,
-      y_loadings = y_loadings,
-      x_projected = x_projected,
-      y_projected = y_projected
-    )))
-
-  } else {
-    stop("Unexpected error in tpls, check 'mode' parameter input")
-  }
-}
diff --git a/R/tens.tplsda.R b/R/tens.tplsda.R
deleted file mode 100644
index 30074f69..00000000
--- a/R/tens.tplsda.R
+++ /dev/null
@@ -1,151 +0,0 @@
-# ==============================================================================
-# Tensor plsda generalization; developed @ Melbourne Integrative Genomics based
-# on Kilmer's tensor m-product algebra
-# ==============================================================================
-
-#' Check if y input is a vector, matrix or tensor that can be inferred as fixed
-#' class labels for each sample across all time points; if so return these
-#' classes in consistent form as a vector.
-#'
-#' @author Brendan Lu
-#' @keywords internal
-#' @noRd
-.y_to_vec <- function(y) {
-  # input is already a n x 1 vector
-  if (is.null(dim(y))) {
-    return(y)
-  }
-  # input is a n x 1 matrix
-  if (length(dim(y)) == 2 && dim(y)[2] == 1) {
-    return(c(y))
-  }
-  # input is a n x 1 x 1 tensor
-  if (length(dim(y)) == 3 && dim(y)[2] == 1 && dim(y)[3] == 1) {
-    return(c(y))
-  }
-  stop("
-    Please ensure y input is either a vector, or a matrix / tensor with one
-    vertical column representing outcome categories for each sample
-
-    If specifying categories across time points, please specify
-    `multilevel = TRUE` in function call and pass in a matrix or tensor with
-    the number of columns corresponding to the number of time points
-  ")
-}
-
-#' Check if y input is a matrix or tensor that can be inferred as repeated class
-#' measurements for each sample; if so, return these classes in a consistent
-#' form as a n x 1 x t tensor.
-#'
-#' Note: We do not have to check that the dimensions are compatible with x
-#' tensor input, this check will be done in the call to tpls.
-#'
-#' @author Brendan Lu
-#' @keywords internal
-#' @noRd
-.y_to_tens <- function(y) {
-  # input is a n x t matrix
-  if (length(dim(y)) == 2) {
-    n <- dim(y)[1]
-    t <- dim(y)[2]
-    return(array(y, dim = c(n, 1, t)))
-  }
-  # input is a n x 1 x t tensor
-  if (length(dim(y)) == 3 && dim(y)[2] == 1) {
-    return(y)
-  }
-  stop("
-    When `multilevel = TRUE` in tplsda, please ensure that y input is either a
-    n x t matrix, or n x 1 x t column tensor, representing the classes for
-    each sample across the t repeated measurements.
-  ")
-}
-
-#' Run tensor PLSDA-like analysis. Note this always returns a compressed-matrix
-#' form output.
-#'
-#' Developed @ Melbourne Integrative Genomics
-#'
-#' @param x Tensor input.
-#' @param y A vector / column matrix / column tensor with n class labels, or a
-#' n x t matrix / n x 1 x t tensor if multilevel = TRUE.
-#' @param multilevel Set to TRUE if y contains repeated class measurements
-#' across the t timepoints specified in the x tensor.
-#' @param ncomp The estimated number of components. ncomp must be explicitly set
-#' as an integer in tpls.
-#' @param m A function which applies an orthogonal tensor tubal transform.
-#' @param minv The inverse of m.
-#' @param center If set to false, the data tensor will not be centralized into
-#' Mean Deviation Form (see Mor et al. 2022). By default, the mean horizontal
-#' slice of the input tensor(s) are subtracted, so that all of the horizontal
-#' slices sum to 0, analgous to centering matrix data.
-#' @param bpparam A \linkS4class{BiocParallelParam} object indicating the type
-#' of parallelisation. Does not have any effect if transform functions
-#' explicitly set using \code{m}, \code{minv}.
-#' @author Brendan Lu
-#' @export
-tplsda <- function(
-  x,
-  y,
-  multilevel = FALSE,
-  ncomp = NULL,
-  m = NULL,
-  minv = NULL,
-  center = TRUE,
-  bpparam = NULL
-) {
-  # most parameter checking is done in tpls call, only the checks that are
-  # unique / necessary for tplsda are done here
-  if (length(dim(x)) != 3) {
-    stop("Please ensure x input tensor is an order-3 array")
-  }
-
-  # BLTODO: may need to vendor unmap.R into big copy-paste script for Saritha
-  # or library(mixOmics) at the top would work neatly too
-  if (multilevel) {
-    y_tens <- .y_to_tens(y)
-    n <- dim(y_tens)[1]
-    t <- dim(y_tens)[3]
-
-    # we need to determine all the unique levels that show up across all
-    # repeated observations for all samples
-    factors_list <- lapply(array(seq_len(t)), function(i) factor(y_tens[, , i]))
-    unique_levels <- unique(unlist(lapply(factors_list, levels)))
-    if (length(unique_levels) == 1) {
-      stop("y should contain 2 or more distinct classes")
-    }
-
-    y <- array(0, dim = c(n, length(unique_levels), t))
-    for (i in seq_len(t)) {
-      # unmap conveniently allows us to manually control the total sets of
-      # groups (class measurements) our input sample is drawn from
-      y[, , i] <- unmap(factors_list[[i]], unique_levels)
-      # bltodo: preserve names, see below
-    }
-  } else {
-    y_vec <- .y_to_vec(y)
-    y_vec <- factor(y_vec)
-    if (nlevels(y_vec) == 1) {
-      stop("y should contain 2 or more distinct classes")
-    }
-    y_mat <- unmap(y_vec)
-    # bltodo: WANT colnames(y_mat) = levels(y_vec), and be able to propogate
-    # names into tpls - need to work on preserving names for tpls and tpca
-    # when mixOmics2 is rolled out to end users (see existing plsda)
-
-    # fill y matrix into a tensor with the same t as x tensor input
-    y <- array(y_mat, dim = c(dim(y_mat), dim(x)[3]))
-  }
-
-  # note the rest of parameter validation will be done in tpls call
-  return(invisible(tpls(
-    x,
-    y,
-    ncomp = ncomp,
-    m = m,
-    minv = minv,
-    mode = "regression",
-    center = center,
-    bpparam = bpparam
-  )))
-}
diff --git a/R/tens.tsvdm.R b/R/tens.tsvdm.R
deleted file mode 100644
index e13c9041..00000000
--- a/R/tens.tsvdm.R
+++ /dev/null
@@ -1,87 +0,0 @@
-# ==============================================================================
-# Kilmer's t-SVD decomposition for order-3 tensors
-# ==============================================================================
-
-#' Tensor SVD-like decomposition algorithm
-#'
-#' Return the t-SVDM decomposition from Kilmer et al. (2021).
-#'
-#' @param x Tensor input.
-#' @param m A function which applies an orthogonal tensor tubal transform.
-#' @param minv The inverse of m.
-#' @param transform Set to FALSE if the input is already in a m-transformed
-#' space, and you just want a facewise svd returned (this is used in tpls).
-#' @param keep_hats Set to TRUE if you do not want minv to be applied to the
-#' output (this is used in tpca), i.e. return outputs in the hat-space.
-#' @param svals_matrix_form Setting to TRUE will return a compressed version of
-#' the singular values tensor S, where the singular values of each f-diagonal
-#' frontal slice becomes the column of a matrix, with t columns in total.
-#' @param bpparam A \linkS4class{BiocParallelParam} object indicating the type
-#' of parallelisation. Does not have any effect if transform functions
-#' explicitly set using \code{m}, \code{minv}.
-#' @author Brendan Lu
-#' @export
-tsvdm <- function(
-  x,
-  m = NULL,
-  minv = NULL,
-  transform = TRUE, # differs to tred, control initial transform to m-space
-  keep_hats = FALSE,
-  svals_matrix_form = FALSE,
-  bpparam = NULL
-) {
-  if (length(dim(x)) != 3) {
-    stop("Please ensure input tensor is an order-3 array.")
-  } else {
-    n <- dim(x)[1]
-    p <- dim(x)[2]
-    t <- dim(x)[3]
-    k <- min(n, p)
-  }
-
-  if (transform) {
-    # use dctii as default transform if user does not specify an explicit one
-    validated_transforms <- .stop_invalid_transform_input(m, minv, t, bpparam)
-    m <- validated_transforms$m
-    minv <- validated_transforms$minv
-
-    x <- m(x)
-  }
-
-  u <- array(0, dim = c(n, k, t))
-  v <- array(0, dim = c(p, k, t))
-
-  # bltodo: less ugly way to write this?
-  if (svals_matrix_form) {
-    s <- array(0, dim = c(k, t))
-    for (i in seq_len(t)) {
-      facewise_svd <- svd(x[, , i])
-      u[, , i] <- facewise_svd$u
-      s[, i] <- facewise_svd$d
-      v[, , i] <- facewise_svd$v
-    }
-  } else {
-    s <- array(0, dim = c(k, k, t))
-    for (i in seq_len(t)) {
-      facewise_svd <- svd(x[, , i])
-      u[, , i] <- facewise_svd$u
-      s[, , i] <- diag(facewise_svd$d)
-      v[, , i] <- facewise_svd$v
-    }
-  }
-
-  if (transform && keep_hats) {
-    # make clear returning in hat space
-    return(list(uhat = u, shat = s, vhat = v))
-  } else {
-    if (transform) {
-      # minv will work on s regardless of what form it is in
-      output <- lapply(list(u, s, v), minv)
-      names(output) <- c("u", "s", "v")
-      return(output)
-    } else {
-      output <- list(u = u, s = s, v = v)
-      return(output)
-    }
-  }
-}
diff --git a/man/block_tpls.Rd b/man/block_tpls.Rd
deleted file mode 100644
index dc32f29b..00000000
--- a/man/block_tpls.Rd
+++ /dev/null
@@ -1,57 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/tens.block.tpls.R
-\name{block_tpls}
-\alias{block_tpls}
-\title{Tensor block PLS}
-\usage{
-block_tpls(
-  x,
-  y,
-  ncomp = NULL,
-  design,
-  m = NULL,
-  minv = NULL,
-  mode = "regression",
-  center = TRUE,
-  bpparam = NULL
-)
-}
-\arguments{
-\item{x}{A list of tensor inputs (termed 'blocks') measured on the same
-samples.}
-
-\item{y}{Tensor input.}
-
-\item{ncomp}{The estimated number of components. ncomp must be explicitly set
-as an integer in tpls.}
-
-\item{design}{Numeric matrix of size length(x) x length(x) with values
-between 0 and 1; the i,j entry in the matrix indicates the strength of the
-relationship to be modelled between the i-th and j-th blocks. A value of 0
-indicates no relationship, with 1 being the maximum. Alternatively, one can
-input "null" for a fully disconnected design, or "full" for a fully connected
-design, or a single scalar value between 0 and which will designate all
-off-diagonal elements of a fully connected design.}
-
-\item{m}{A function which applies an orthogonal tensor tubal transform.}
-
-\item{minv}{The inverse of m.}
-
-\item{mode}{Currently supports tensor analogues of canonical, regression,
-and svd PLS modes. Defaults to "regression" mode.}
-
-\item{center}{If set to false, the data tensor will not be centralized into
-Mean Deviation Form (see Mor et al. 2022). By default, the mean horizontal
-slice of the input tensor(s) are subtracted, so that all of the horizontal
-slices sum to 0, analgous to centering matrix data.}
-
-\item{bpparam}{A \linkS4class{BiocParallelParam} object indicating the type
-of parallelisation. Does not have any effect if transform functions
-explicitly set using \code{m}, \code{minv}.}
-}
-\description{
-Developed @ Melbourne Integrative Genomics.
-}
-\author{
-Brendan Lu
-}
diff --git a/man/dctii_m_transforms.Rd b/man/dctii_m_transforms.Rd
deleted file mode 100644
index a1c2f56d..00000000
--- a/man/dctii_m_transforms.Rd
+++ /dev/null
@@ -1,28 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/tens.mproduct.R
-\name{dctii_m_transforms}
-\alias{dctii_m_transforms}
-\title{Create m-transform functions to apply the Discrete Cosine Transform 2}
-\usage{
-dctii_m_transforms(t, bpparam = NULL)
-}
-\arguments{
-\item{t}{The length of the transform.}
-
-\item{bpparam}{A \linkS4class{BiocParallelParam} object indicating the type
-of parallelisation.}
-}
-\value{
-\item{m}{A function which applies the dct-ii along the last dimension of a
-given numerical input array. For 3D tensors it performs the mode-3 product
-with the DCT matrix.}
-\item{minv}{The inverse of m,}
-}
-\description{
-Returns functions \code{m} and \code{m_inv} which apply tubal transforms
-defined by the Discrete Cosine Transform (DCT-II variant). This is equivalent
-to Scipys DCTI-ii algorithm with \code{norm='ortho'}.
-}
-\author{
-Brendan Lu
-}
diff --git a/man/facewise_crossproduct.Rd b/man/facewise_crossproduct.Rd
deleted file mode 100644
index eae42324..00000000
--- a/man/facewise_crossproduct.Rd
+++ /dev/null
@@ -1,22 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/tens.mproduct.R
-\name{facewise_crossproduct}
-\alias{facewise_crossproduct}
-\title{Tensor cross product}
-\usage{
-facewise_crossproduct(a, b)
-}
-\arguments{
-\item{a}{Tensor input.}
-
-\item{b}{Tensor input.}
-}
-\value{
-Tensor facewise cross product.
-}
-\description{
-Compute the equivalent of ft(a) %fp% b
-}
-\author{
-Brendan Lu
-}
diff --git a/man/facewise_product.Rd b/man/facewise_product.Rd
deleted file mode 100644
index 4ad53883..00000000
--- a/man/facewise_product.Rd
+++ /dev/null
@@ -1,36 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/tens.mproduct.R
-\name{facewise_product}
-\alias{facewise_product}
-\alias{\%fp\%}
-\title{Tensor facewise product}
-\usage{
-facewise_product(..., bpparam = NULL)
-
-a \%fp\% b
-}
-\arguments{
-\item{...}{Arbitrary number of numerical tensor inputs.}
-
-\item{bpparam}{A \linkS4class{BiocParallelParam} object indicating the type
-of parallelisation.}
-
-\item{a}{Tensor input.}
-
-\item{b}{Tensor input.}
-}
-\value{
-Cumulative facewise product.
-}
-\description{
-Compute Kilmer's facewise product cumulatively across any arbitrary number of
-tensor inputs.
-}
-\section{Functions}{
-\itemize{
-\item \code{a \%fp\% b}: Custom facewise product operator
-
-}}
-\author{
-Brendan Lu
-}
diff --git a/man/facewise_transpose.Rd b/man/facewise_transpose.Rd
deleted file mode 100644
index 8210cb6e..00000000
--- a/man/facewise_transpose.Rd
+++ /dev/null
@@ -1,28 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/tens.mproduct.R
-\name{facewise_transpose}
-\alias{facewise_transpose}
-\alias{ft}
-\title{Facewise transpose an order-3 tensor}
-\usage{
-facewise_transpose(tensor)
-
-ft(tensor)
-}
-\arguments{
-\item{tensor}{Numerical 3D array input.}
-}
-\value{
-Facewise transpose of \code{tensor}
-}
-\description{
-Facewise transpose an order-3 tensor
-}
-\section{Functions}{
-\itemize{
-\item \code{ft()}: Alias for \code{\link{facewise_product}}
-
-}}
-\author{
-Brendan Lu
-}
diff --git a/man/m_product.Rd b/man/m_product.Rd
deleted file mode 100644
index a591be68..00000000
--- a/man/m_product.Rd
+++ /dev/null
@@ -1,29 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/tens.mproduct.R
-\name{m_product}
-\alias{m_product}
-\title{Kilmer's tensor-tensor m-product}
-\usage{
-m_product(..., m = NULL, minv = NULL, bpparam = NULL)
-}
-\arguments{
-\item{...}{Arbitrary number of numerical tensor inputs.}
-
-\item{m}{A function which applies an orthogonal tensor tubal transform.}
-
-\item{minv}{The inverse of m.}
-
-\item{bpparam}{A \linkS4class{BiocParallelParam} object indicating the type
-of parallelisation. Does not have any effect if transform functions
-explicitly set using \code{m}, \code{minv}.}
-}
-\value{
-Cumulative m-product.
-}
-\description{
-This function works cumulatively across any arbitrary number of tensor
-inputs.
-}
-\author{
-Brendan Lu
-}
diff --git a/man/matrix_to_m_transforms.Rd b/man/matrix_to_m_transforms.Rd
deleted file mode 100644
index fce4e337..00000000
--- a/man/matrix_to_m_transforms.Rd
+++ /dev/null
@@ -1,29 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/tens.mproduct.R
-\name{matrix_to_m_transforms}
-\alias{matrix_to_m_transforms}
-\title{Create m-transform functions as defined by a matrix}
-\usage{
-matrix_to_m_transforms(m_mat, m_inv_mat = NULL, bpparam = NULL)
-}
-\arguments{
-\item{m_mat}{Function which defines the tubal transform.}
-
-\item{m_inv_mat}{Function which defines inverse tubal transform}
-
-\item{bpparam}{A \linkS4class{BiocParallelParam} object indicating the type
-of parallelisation.}
-}
-\value{
-\item{m}{A function which applies the matrix m_mat along the last dimension
-of a given numerical input array. For 3D tensors it performs the mode-3
-product.}
-\item{minv}{The inverse of m.}
-}
-\description{
-Returns functions \code{m} and \code{m_inv} which apply tubal transforms
-defined by the matrix \code{m_mat}.
-}
-\author{
-Brendan Lu
-}
diff --git a/man/tpca.Rd b/man/tpca.Rd
deleted file mode 100644
index 599af9ed..00000000
--- a/man/tpca.Rd
+++ /dev/null
@@ -1,48 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/tens.tpca.R
-\name{tpca}
-\alias{tpca}
-\title{Tensor PCA-like dimensionality reduction}
-\usage{
-tpca(
-  x,
-  ncomp = NULL,
-  m = NULL,
-  minv = NULL,
-  center = TRUE,
-  matrix_output = TRUE,
-  bpparam = NULL
-)
-}
-\arguments{
-\item{x}{Tensor input.}
-
-\item{ncomp}{The estimated number of components. ncomp can be explicitly set
-using an integer value, 0 < float value < 1, or left as NULL which will
-default to the maximum ncomp value possible.}
-
-\item{m}{A function which applies an orthogonal tensor tubal transform.}
-
-\item{minv}{The inverse of m.}
-
-\item{center}{If set to false, the data tensor will not be centralized into
-Mean Deviation Form (see Mor et al. 2022). By default, the mean horizontal
-slice of the input tensor(s) are subtracted, so that all of the horizontal
-slices sum to 0, analgous to centering matrix data.}
-
-\item{matrix_output}{Collect the top singular vectors across the tensor,
-organized by the magnitude of the corresponding singular vector, and place
-them into the columns of a matrix. Corresponds to the 'matrix compression'
-type of scheme described in Mor et al. 2022.}
-
-\item{bpparam}{A \linkS4class{BiocParallelParam} object indicating the type
-of parallelisation. Does not have any effect if transform functions
-explicitly set using \code{m}, \code{minv}.}
-}
-\description{
-Tensor analogue of PCA introduced by Mor et al. (2022) based on Kilmer's
-m-product algebra and tsvdm.
-}
-\author{
-Brendan Lu
-}
diff --git a/man/tpls.Rd b/man/tpls.Rd
deleted file mode 100644
index 2877f787..00000000
--- a/man/tpls.Rd
+++ /dev/null
@@ -1,59 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/tens.tpls.R
-\name{tpls}
-\alias{tpls}
-\title{Tensor PLS-like regression}
-\usage{
-tpls(
-  x,
-  y,
-  ncomp = NULL,
-  m = NULL,
-  minv = NULL,
-  mode = "regression",
-  center = TRUE,
-  matrix_output = TRUE,
-  bpparam = NULL
-)
-}
-\arguments{
-\item{x}{Tensor input.}
-
-\item{y}{Tensor input.}
-
-\item{ncomp}{The estimated number of components. ncomp must be explicitly set
-as an integer in tpls.}
-
-\item{m}{A function which applies an orthogonal tensor tubal transform.}
-
-\item{minv}{The inverse of m.}
-
-\item{mode}{Currently supports tensor analogues of canonical ("canonical"),
-regression ("regression"), and svd ("tsvdm") PLS variants. Defaults to
-"regression".}
-
-\item{center}{If set to false, the data tensor will not be centralized into
-Mean Deviation Form (see Mor et al. 2022). By default, the mean horizontal
-slice of the input tensor(s) are subtracted, so that all of the horizontal
-slices sum to 0, analgous to centering matrix data.}
-
-\item{matrix_output}{Note: ONLY AFFECTS THE OUTPUT IN "tsvdm" MODE.
-Collects the top singular vectors across the tensor, organized by the
-magnitude of the corresponding singular vector, and place them into the
-columns of a matrix. Corresponds to the 'matrix compression' type of scheme
-described in Mor et al. 2022.}
-
-\item{bpparam}{A \linkS4class{BiocParallelParam} object indicating the type
-of parallelisation. Does not have any effect if transform functions
-explicitly set using \code{m}, \code{minv}.}
-}
-\description{
-Developed @ Melbourne Integrative Genomics. Intuition: Fourier-based m
-transforms compact the data across time points, but preserve the information
-association for each sample and each feature. In doing so, pls works by
-calculating the singular vectors associated with a 'global' (across time
-points) top singular value, and deflating just that relevant face.
-}
-\author{
-Brendan Lu
-}
diff --git a/man/tplsda.Rd b/man/tplsda.Rd
deleted file mode 100644
index 06983acd..00000000
--- a/man/tplsda.Rd
+++ /dev/null
@@ -1,49 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/tens.tplsda.R
-\name{tplsda}
-\alias{tplsda}
-\title{Run tensor PLSDA-like analysis. Note this always returns a compressed-matrix
-form output.}
-\usage{
-tplsda(
-  x,
-  y,
-  multilevel = FALSE,
-  ncomp = NULL,
-  m = NULL,
-  minv = NULL,
-  center = TRUE,
-  bpparam = NULL
-)
-}
-\arguments{
-\item{x}{Tensor input.}
-
-\item{y}{A vector / column matrix / column tensor with n class labels, or a
-n x t matrix / n x 1 x t tensor if multilevel = TRUE.}
-
-\item{multilevel}{Set to TRUE if y contains repeated class measurements
-across the t timepoints specified in the x tensor.}
-
-\item{ncomp}{The estimated number of components. ncomp must be explicitly set
-as an integer in tpls.}
-
-\item{m}{A function which applies an orthogonal tensor tubal transform.}
-
-\item{minv}{The inverse of m.}
-
-\item{center}{If set to false, the data tensor will not be centralized into
-Mean Deviation Form (see Mor et al. 2022). By default, the mean horizontal
-slice of the input tensor(s) are subtracted, so that all of the horizontal
-slices sum to 0, analgous to centering matrix data.}
-
-\item{bpparam}{A \linkS4class{BiocParallelParam} object indicating the type
-of parallelisation. Does not have any effect if transform functions
-explicitly set using \code{m}, \code{minv}.}
-}
-\description{
-Developed @ Melbourne Integrative Genomics
-}
-\author{
-Brendan Lu
-}
diff --git a/man/tsvdm.Rd b/man/tsvdm.Rd
deleted file mode 100644
index 30be0b9f..00000000
--- a/man/tsvdm.Rd
+++ /dev/null
@@ -1,43 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/tens.tsvdm.R
-\name{tsvdm}
-\alias{tsvdm}
-\title{Tensor SVD-like decomposition algorithm}
-\usage{
-tsvdm(
-  x,
-  m = NULL,
-  minv = NULL,
-  transform = TRUE,
-  keep_hats = FALSE,
-  svals_matrix_form = FALSE,
-  bpparam = NULL
-)
-}
-\arguments{
-\item{x}{Tensor input.}
-
-\item{m}{A function which applies an orthogonal tensor tubal transform.}
-
-\item{minv}{The inverse of m.}
-
-\item{transform}{Set to FALSE if the input is already in a m-transformed
-space, and you just want a facewise svd returned (this is used in tpls).}
-
-\item{keep_hats}{Set to TRUE if you do not want minv to be applied to the
-output (this is used in tpca), i.e. return outputs in the hat-space.}
-
-\item{svals_matrix_form}{Setting to TRUE will return a compressed version of
-the singular values tensor S, where the singular values of each f-diagonal
-frontal slice becomes the column of a matrix, with t columns in total.}
-
-\item{bpparam}{A \linkS4class{BiocParallelParam} object indicating the type
-of parallelisation. Does not have any effect if transform functions
-explicitly set using \code{m}, \code{minv}.}
-}
-\description{
-Return the t-SVDM decomposition from Kilmer et al. (2021).
-}
-\author{
-Brendan Lu
-}
diff --git a/tests/testthat/test-tens.block.tpls.R b/tests/testthat/test-tens.block.tpls.R
deleted file mode 100644
index 90ff0082..00000000
--- a/tests/testthat/test-tens.block.tpls.R
+++ /dev/null
@@ -1,8 +0,0 @@
-context("block tpls")
-
-test_that(
-  "block tpls",
-  code = {
-    NULL
-  }
-)
diff --git a/tests/testthat/test-tens.mproduct.R b/tests/testthat/test-tens.mproduct.R
deleted file mode 100644
index c74a74b9..00000000
--- a/tests/testthat/test-tens.mproduct.R
+++ /dev/null
@@ -1,208 +0,0 @@
-context("m-product utilities")
-# bltodo: add tests for matrix inputs
-# bltodo: run on linux system to test parallel algorithms
-
-#' @description Use for internal testing.
-#' Performs a DCT-II transform using the stats::fft algorithm. Produces
-#' identical output to the scipy.fft.dct implementation in Python.
-#' @param vec A numeric real vector to transform.
-#' @param ortho Logical, if TRUE the output is orthonormal.
-#' @return A numeric vector representing the DCT-II of the input.
-#' @seealso docs.scipy.org/doc/scipy/tutorial/fft.html
-dctii_scipy_equivalent <- function(vec, ortho = TRUE) {
-  # https://scipy.github.io/devdocs/reference/generated/scipy.fftpack.dct.html
-  # https://stackoverflow.com/questions/11215162
-  # NOTE: this is dtt output scaled by 2
-  n <- length(vec)
-  p <- exp(complex(imaginary = pi / 2 / n) * (seq(2 * n) - 1))
-  unscaled_res <- Re(stats::fft(c(vec, rev(vec))) / p)[1:n]
-
-  if (ortho) {
-    return(
-      c(
-        sqrt(1 / (4 * n)) * unscaled_res[1],
-        sqrt(1 / (2 * n)) * tail(unscaled_res, n - 1)
-      )
-    )
-  } else {
-    return(unscaled_res)
-  }
-}
-
-test_that(
-  "gsignal's dct2 is producing the same result as the Scipy algorithm",
-  code = {
-    test_vector <- array(c(19, 3, 6, 11))
-    transforms <- dctii_m_transforms(length(test_vector))
-    m <- transforms$m
-    expect_equal(m(test_vector), matrix(dctii_scipy_equivalent(test_vector)))
-  }
-)
-
-test_that(
-  "facewise product works the same as the naive algorithm",
-  code = {
-    n <- 2
-    p <- 4
-    t <- 3
-    set.seed(1)
-    test_tensor1 <- array(rnorm(n * p * t), dim = c(n, p, t))
-    test_tensor2 <- ft(test_tensor1)
-    expected_result <- array(0, dim = c(n, n, t))
-    for (i in 1:t) {
-      expected_result[, , i] <- test_tensor1[, , i] %*% test_tensor2[, , i]
-    }
-    expect_equal(
-      test_tensor1 %fp% test_tensor2,
-      expected_result
-    )
-  }
-)
-
-test_that(
-  "mode-3 product result matches naive nested for-loop algorithm",
-  code = {
-    n <- 2
-    p <- 4
-    t <- 3
-    test_tensor1 <- array(1:(n * p * t), dim = c(2, 4, 3))
-    m_mat <- gsignal::dctmtx(t)
-    expected_result <- array(0, dim = c(2, 4, 3))
-    for (nn in 1:n) {
-      for (pp in 1:p) {
-        expected_result[nn, pp, ] <- m_mat %*% test_tensor1[nn, pp, ]
-      }
-    }
-    transforms <- dctii_m_transforms(t)
-    m <- transforms$m
-    expect_equal(m(test_tensor1), expected_result)
-  }
-)
-
-test_that(
-  "different mode-3 product algorithms produce equivalent results",
-  code = {
-    test_tensor1 <- array(1:24, dim = c(2, 4, 3))
-    t <- dim(test_tensor1)[3]
-    transforms_default <- dctii_m_transforms(t, bpparam = NULL)
-    transforms_parallel <- dctii_m_transforms(
-      t, bpparam = BiocParallel::MulticoreParam()
-    )
-    expect_equal(
-      transforms_default$m(test_tensor1),
-      # suppress warning "MulticoreParam() not supported on Windows, use
-      # SnowParam()" when running tests on Windows
-      suppressWarnings(
-        transforms_parallel$m(test_tensor1)
-      )
-    )
-  }
-)
-
-test_that(
-  "forward and inverse dctii transforms invert each other",
-  code = {
-    test_tensor1 <- array(1:24, dim = c(2, 4, 3))
-    transforms <- dctii_m_transforms(dim(test_tensor1)[3])
-    m <- transforms$m
-    minv <- transforms$minv
-    expect_equal(test_tensor1, minv(m(test_tensor1)))
-    expect_equal(test_tensor1, m(minv(test_tensor1)))
-  }
-)
-
-test_that(
-  # note this test does tend to throw warning messages on Windows
-  "different binary facewise algorithms produce equivalent results",
-  code = {
-    test_tensor1 <- array(1:24, dim = c(2, 4, 3))
-    test_tensor2 <- array(1:36, dim = c(4, 3, 3))
-    expect_equal(
-      facewise_product(test_tensor1, test_tensor2, bpparam = NULL),
-      # suppress warning "MulticoreParam() not supported on Windows, use
-      # SnowParam()" when running tests on Windows
-      suppressWarnings(
-        facewise_product(
-          test_tensor1, test_tensor2,
-          bpparam = BiocParallel::MulticoreParam()
-        )
-      )
-    )
-  }
-)
-
-test_that(
-  "cumulative multi input facewise product works as expected compared to 
-  custom operator",
-  code = {
-    test_tensor1 <- array(1:24, dim = c(2, 4, 3))
-    test_tensor2 <- array(1:60, dim = c(4, 5, 3))
-    test_tensor3 <- array(1:90, dim = c(5, 6, 3))
-    fp12 <- facewise_product(test_tensor1, test_tensor2)
-    expected_cumulative_fp <- facewise_product(fp12, test_tensor3)
-    expect_equal(
-      facewise_product(test_tensor1, test_tensor2, test_tensor3),
-      expected_cumulative_fp
-    )
-    expect_equal(
-      test_tensor1 %fp% test_tensor2 %fp% test_tensor3,
-      expected_cumulative_fp
-    )
-  }
-)
-
-test_that(
-  "cumulative multi input m product works as expected",
-  code = {
-    test_tensor1 <- array(1:24, dim = c(2, 4, 3))
-    test_tensor2 <- array(1:60, dim = c(4, 5, 3))
-    test_tensor3 <- array(1:90, dim = c(5, 6, 3))
-    mp12 <- m_product(test_tensor1, test_tensor2)
-    expected_cumulative_mp <- m_product(mp12, test_tensor3)
-    expect_equal(
-      m_product(test_tensor1, test_tensor2, test_tensor3),
-      expected_cumulative_mp
-    )
-  }
-)
-
-test_that(
-  "m_product() throws appropriate errors",
-  code = {
-    test_tensor1 <- array(1:24, dim = c(2, 4, 3))
-    dummy_m <- function() 0
-    # only defining m without minv
-    expect_error(
-      m_product(test_tensor1, m = dummy_m),
-      "If explicitly defined, both m and its inverse must be defined as 
-      functions"
-    )
-    # accidentally adding () to an input meaning it is not a callable
-    expect_error(
-      m_product(test_tensor1, m = dummy_m()),
-      "If explicitly defined, both m and its inverse must be defined as 
-      functions"
-    )
-    expect_error(
-      m_product(test_tensor1, m = dummy_m(), minv = dummy_m),
-      "If explicitly defined, both m and its inverse must be defined as 
-      functions"
-    )
-  }
-)
-
-test_that(
-  "facewise_crossproduct works as intended",
-  code = {
-    n <- 2
-    p <- 4
-    t <- 3
-    set.seed(1)
-    test_tensor1 <- array(rnorm(n * p * t), dim = c(n, p, t))
-
-    expect_equal(
-      facewise_crossproduct(test_tensor1, test_tensor1),
-      ft(test_tensor1) %fp% test_tensor1
-    )
-  }
-)
diff --git a/tests/testthat/test-tens.tpca.R b/tests/testthat/test-tens.tpca.R
deleted file mode 100644
index 321aef1a..00000000
--- a/tests/testthat/test-tens.tpca.R
+++ /dev/null
@@ -1,63 +0,0 @@
-context("tpca")
-
-#' Unnames and ensures the top row in a matrix-type output is all positive. This
-#' allows for comparison between pls results that are the same but just differ
-#' by a negative sign due to svd solver or other implementation detail.
-.make_signs_consistent <- function(mat) {
-  mat <- unname(mat)
-  ncols <- dim(mat)[2]
-  for (i in seq_len(ncols)) {
-    if (mat[1, i] < 0) {
-      mat[, i] <- -mat[, i]
-    }
-  }
-  return(mat)
-}
-
-test_that(
-  "basic tpca sense checks",
-  code = {
-    n <- 4
-    p <- 5
-    t <- 3
-    ncomp_input <- 2
-    test_tensor <- array(1:(n * p * t), dim = c(n, p, t))
-    tpca_obj <- tpca(test_tensor, ncomp = ncomp_input)
-    expect_equal(length(tpca_obj$explained_variance), ncomp_input)
-    expect_equal(dim(tpca_obj$variates), c(n, ncomp_input))
-    expect_equal(dim(tpca_obj$loadings), c(p, ncomp_input))
-  }
-)
-
-test_that(
-  "tpca agrees with mixOmics pca",
-  code = {
-    n <- 3
-    p <- 10
-    t <- 1
-    k <- min(n, p)
-    ncomp_input <- 2
-
-    set.seed(1)
-    test_x <- array(rnorm(n * p * t, mean = 0, sd = 5), dim = c(n, p, t))
-
-    mixomics_pca <- pca(
-      test_x[, , 1],
-      ncomp = ncomp_input
-    )
-
-    # bltodo: we do not even need to specify identity transforms here right?
-    transforms <- matrix_to_m_transforms(diag(1))
-    tensor_pca <- tpca(
-      test_x,
-      ncomp = ncomp_input,
-      m = transforms$m,
-      minv = transforms$minv
-    )
-
-    expect_equal(
-      .make_signs_consistent(mixomics_pca$loadings$X),
-      .make_signs_consistent(tensor_pca$loadings)
-    )
-  }
-)
diff --git a/tests/testthat/test-tens.tpls.R b/tests/testthat/test-tens.tpls.R
deleted file mode 100644
index 5a8b3a1e..00000000
--- a/tests/testthat/test-tens.tpls.R
+++ /dev/null
@@ -1,244 +0,0 @@
-context("tpls")
-
-#' Unnames and ensures the top row in a matrix-type output is all positive. This
-#' allows for comparison between pls results that are the same but just differ
-#' by a negative sign due to svd solver or other implementation detail.
-.make_signs_consistent <- function(mat) {
-  mat <- unname(mat)
-  ncols <- dim(mat)[2]
-  for (i in seq_len(ncols)) {
-    if (mat[1, i] < 0) {
-      mat[, i] <- -mat[, i]
-    }
-  }
-  return(mat)
-}
-
-test_that(
-  "tsvdm-tpls mode agrees with tpca",
-  code = {
-    # it actually takes a lot of coercing to make these two outputs comparable
-    # because their default configurations are aimed at completely different
-    # problems. that being said, this test here is probably still useful as a
-    # sense check.
-    #
-    # problem 1: squaring the singular values means that the k_t_flatten_sort
-    # compression is different and picks out a different compressed matrix. Sort
-    # order of the squared values is not necessarily the same.
-    #
-    # problem 2: loadings (and therefore projections) may differ in sign, think
-    # this is dependent on the svd solver implementation?
-    #
-    # problem 3: if you use centering (on for tpca and tpls by default), the
-    # faces of each tensor loses one rank, so only rank - 1 columns of each face
-    # in the loadings tensor ("v") will match. for tsvdm mode this means we get
-    # one less column to compare.
-    #
-    # problem 4: tpls computes the svd on XtX, which is rank-deficient if X is
-    # non-square. this means that the last few columns of "v" are absolutely
-    # crap and should not be compared. these columns may not even exist in
-    # tpca depending on the input dimensions.
-    n <- 4
-    p <- 5
-    t <- 3
-    k <- min(n, p)
-
-    # this test fails if `test_tensor` is rank deficient, as is the case if we
-    # use array(1:24, dim = c(3, 4, 2))
-    set.seed(1)
-    test_tensor <- array(rnorm(n * p * t, mean = 0, sd = 5), dim = c(n, p, t))
-
-    transforms <- dctii_m_transforms(t)
-    m <- transforms$m
-    minv <- transforms$minv
-
-    tpca_obj <- tpca(
-      test_tensor,
-      m = m,
-      minv = minv,
-      # turn off centering to get over problem 3
-      center = FALSE,
-      # return full tensors to get over problem 1
-      matrix_output = FALSE
-    )
-
-    tpls_tsvdm <- tpls(
-      test_tensor, test_tensor,
-      m = m,
-      minv = minv,
-      mode = "tsvdm",
-      # turn off centering to get over problem 3
-      center = FALSE,
-      # return full tensors to get over problem 1
-      matrix_output = FALSE
-    )
-
-    # compare elementwise absolute values to get over problem 2
-    expect_equal(
-      abs(tpca_obj$loadings),
-      # only compare k columns to get over problem 4
-      abs(tpls_tsvdm$y_loadings[, 1:k, ])
-    )
-
-    expect_equal(
-      abs(tpca_obj$variates),
-      abs(tpls_tsvdm$y_projected)
-    )
-  }
-)
-
-test_that(
-  "the first component of tpls regression and canonical match tpca",
-  code = {
-    n <- 4
-    p <- 5
-    t <- 3
-    k <- min(n, p)
-
-    set.seed(1)
-    test_tensor <- array(rnorm(n * p * t, mean = 0, sd = 5), dim = c(n, p, t))
-
-    transforms <- dctii_m_transforms(t)
-    m <- transforms$m
-    minv <- transforms$minv
-
-    tpca_obj <- tpca(
-      test_tensor,
-      m = m,
-      minv = minv
-    )
-
-    tpls_regression <- tpls(
-      test_tensor, test_tensor,
-      m = m,
-      minv = minv,
-      mode = "regression",
-    )
-
-    expect_equal(tpls_regression$y_loadings[, 1], tpca_obj$loadings[, 1])
-
-    expect_equal(
-      .make_signs_consistent(tpls_regression$y_loadings[, 1:3]),
-      .make_signs_consistent(tpca_obj$loadings[, 1:3])
-    )
-
-    tpls_canonical <- tpls(
-      test_tensor, test_tensor,
-      m = m,
-      minv = minv,
-      mode = "canonical",
-    )
-
-    expect_equal(tpls_canonical$y_loadings[, 1], tpca_obj$loadings[, 1])
-  }
-)
-
-test_that(
-  "canonical: tpls agrees with MixOmics pls",
-  code = {
-    n <- 6
-    p <- 7
-    q <- 9
-    t <- 1
-    k <- min(n, p, q)
-    ncomp_input <- 4
-
-    set.seed(1)
-    test_x <- array(rnorm(n * p * t, mean = 0, sd = 5), dim = c(n, p, t))
-    test_y <- array(rnorm(n * q * t, mean = 0, sd = 3), dim = c(n, q, t))
-
-    mixomics_pls <- pls(
-      test_x[, , 1],
-      test_y[, , 1],
-      ncomp = ncomp_input,
-      scale = FALSE,
-      mode = "canonical"
-    )
-
-    # bltodo: we do not even need to specify identity transforms here right?
-    transforms <- matrix_to_m_transforms(diag(1))
-    tensor_pls <- tpls(
-      test_x,
-      test_y,
-      ncomp = ncomp_input,
-      m = transforms$m,
-      minv = transforms$minv,
-      mode = "canonical"
-    )
-
-    expect_equal(
-      .make_signs_consistent(mixomics_pls$loadings$X),
-      .make_signs_consistent(tensor_pls$x_loadings)
-    )
-
-    expect_equal(
-      .make_signs_consistent(mixomics_pls$loadings$Y),
-      .make_signs_consistent(tensor_pls$y_loadings)
-    )
-
-    expect_equal(
-      .make_signs_consistent(mixomics_pls$variates$X),
-      .make_signs_consistent(tensor_pls$x_projected)
-    )
-
-    expect_equal(
-      .make_signs_consistent(mixomics_pls$variates$Y),
-      .make_signs_consistent(tensor_pls$y_projected)
-    )
-  }
-)
-
-test_that(
-  "regression: tpls agrees with MixOmics pls",
-  code = {
-    n <- 5
-    p <- 8
-    q <- 10
-    t <- 1
-    k <- min(n, p, q)
-    ncomp_input <- 4
-
-    set.seed(1)
-    test_x <- array(rnorm(n * p * t, mean = 0, sd = 5), dim = c(n, p, t))
-    test_y <- array(rnorm(n * q * t, mean = 0, sd = 3), dim = c(n, q, t))
-
-    mixomics_pls <- pls(
-      test_x[, , 1],
-      test_y[, , 1],
-      ncomp = ncomp_input,
-      scale = FALSE,
-      mode = "regression"
-    )
-
-    # bltodo: we do not even need to specify identity transforms here right?
-    transforms <- matrix_to_m_transforms(diag(1))
-    tensor_pls <- tpls(
-      test_x,
-      test_y,
-      ncomp = ncomp_input,
-      m = transforms$m,
-      minv = transforms$minv,
-      mode = "regression"
-    )
-
-    expect_equal(
-      .make_signs_consistent(mixomics_pls$loadings$X),
-      .make_signs_consistent(tensor_pls$x_loadings)
-    )
-
-    expect_equal(
-      .make_signs_consistent(mixomics_pls$loadings$Y),
-      .make_signs_consistent(tensor_pls$y_loadings)
-    )
-
-    expect_equal(
-      .make_signs_consistent(mixomics_pls$variates$X),
-      .make_signs_consistent(tensor_pls$x_projected)
-    )
-
-    expect_equal(
-      .make_signs_consistent(mixomics_pls$variates$Y),
-      .make_signs_consistent(tensor_pls$y_projected)
-    )
-  }
-)
diff --git a/tests/testthat/test-tens.tplsda.R b/tests/testthat/test-tens.tplsda.R
deleted file mode 100644
index 9f7d9c59..00000000
--- a/tests/testthat/test-tens.tplsda.R
+++ /dev/null
@@ -1,115 +0,0 @@
-context("tplsda")
-
-#' Unnames and ensures the top row in a matrix-type output is all positive. This
-#' allows for comparison between pls results that are the same but just differ
-#' by a negative sign due to svd solver or other implementation detail.
-.make_signs_consistent <- function(mat) {
-  mat <- unname(mat)
-  ncols <- dim(mat)[2]
-  for (i in seq_len(ncols)) {
-    if (mat[1, i] < 0) {
-      mat[, i] <- -mat[, i]
-    }
-  }
-  return(mat)
-}
-
-test_that(
-  "tplsda agrees with mixOmics plsda",
-  code = {
-    n <- 6
-    p <- 7
-    q <- 9
-    t <- 1
-    k <- min(n, p, q)
-    ncomp_input <- 4
-
-    set.seed(1)
-    test_x <- array(rnorm(n * p * t, mean = 0, sd = 5), dim = c(n, p, t))
-    test_y <- letters[1:n]
-
-    mixomics_plsda <- plsda(
-      test_x[, , 1],
-      test_y,
-      ncomp = ncomp_input,
-      scale = FALSE
-    )
-
-    # bltodo: we do not even need to specify identity transforms here right?
-    transforms <- matrix_to_m_transforms(diag(1))
-    tensor_plsda <- tplsda(
-      test_x,
-      test_y,
-      ncomp = ncomp_input,
-      m = transforms$m,
-      minv = transforms$minv
-    )
-
-    expect_equal(
-      .make_signs_consistent(mixomics_plsda$loadings$X),
-      .make_signs_consistent(tensor_plsda$x_loadings)
-    )
-
-    expect_equal(
-      .make_signs_consistent(mixomics_plsda$loadings$Y),
-      .make_signs_consistent(tensor_plsda$y_loadings)
-    )
-
-    expect_equal(
-      .make_signs_consistent(mixomics_plsda$variates$X),
-      .make_signs_consistent(tensor_plsda$x_projected)
-    )
-
-    expect_equal(
-      .make_signs_consistent(mixomics_plsda$variates$Y),
-      .make_signs_consistent(tensor_plsda$y_projected)
-    )
-  }
-)
-
-test_that(
-  "tplsda single time point y imputation is consistent",
-  code = {
-    n <- 6
-    p <- 7
-    t <- 3
-
-    set.seed(1)
-    test_x <- array(rnorm(n * p * t, mean = 0, sd = 5), dim = c(n, p, t))
-    test_y_vec <- letters[1:n]
-    test_y_mat <- array(test_y_vec, dim = c(n, 1))
-    test_y_tens <- array(test_y_vec, dim = c(n, 1, 1))
-
-    tplsda_vec <- tplsda(test_x, test_y_vec)
-    tplsda_mat <- tplsda(test_x, test_y_mat)
-    tplsda_tens <- tplsda(test_x, test_y_tens)
-
-    expect_equal(tplsda_vec$y, tplsda_mat$y)
-    expect_equal(tplsda_mat$y, tplsda_tens$y)
-  }
-)
-
-test_that(
-  "tplsda multilevel y transformation is appropriate",
-  code = {
-    n <- 6
-    p <- 7
-    t <- 3
-
-    set.seed(1)
-    test_x <- array(rnorm(n * p * t, mean = 0, sd = 5), dim = c(n, p, t))
-    test_y_mat <- array(letters[1:(n * t)], dim = c(n, t))
-    test_y_tens <- array(letters[1:(n * t)], dim = c(n, 1, t))
-
-    tplsda_mat <- tplsda(test_x, test_y_mat, multilevel = TRUE)
-    tplsda_tens <- tplsda(test_x, test_y_tens, multilevel = TRUE)
-
-    # both multilevel y inputs should lead to the same underlying tpls call
-    expect_equal(tplsda_mat$y, tplsda_tens$y)
-
-    # the y tensor that goes into the tpls call should have the appropriate
-    # number of columns corresponding to all unique categories found across
-    # all timepoints specified in y
-    expect_equal(dim(tplsda_mat$y), c(n, n * t, 3))
-  }
-)
diff --git a/tests/testthat/test-tens.tsvdm.R b/tests/testthat/test-tens.tsvdm.R
deleted file mode 100644
index e358598d..00000000
--- a/tests/testthat/test-tens.tsvdm.R
+++ /dev/null
@@ -1,37 +0,0 @@
-context("tsvdm")
-# bltodo: add tests for permutations of different configurations and input sizes
-# like the `tred` library
-
-test_that(
-  "reconstructing original matrix from tsvdm",
-  code = {
-    test_tensor <- array(1:24, dim = c(2, 4, 3))
-    tsvdm_decomposition <- tsvdm(test_tensor)
-    u <- tsvdm_decomposition$u
-    s <- tsvdm_decomposition$s
-    v <- tsvdm_decomposition$v
-    vt <- ft(v)
-
-    expect_equal(
-      test_tensor,
-      m_product(u, s, vt)
-    )
-  }
-)
-
-test_that(
-  "reconstructing original matrix from tsvdm without transform",
-  code = {
-    test_tensor <- array(1:24, dim = c(2, 4, 3))
-    tsvdm_decomposition <- tsvdm(test_tensor, transform = FALSE)
-    u <- tsvdm_decomposition$u
-    s <- tsvdm_decomposition$s
-    v <- tsvdm_decomposition$v
-    vt <- ft(v)
-
-    expect_equal(
-      test_tensor,
-      facewise_product(u, s, vt)
-    )
-  }
-)