R/ossa.R

Defines functions owcor rowspan.ossa colspan.ossa decompose.ossa fix.iossa.groups fossa.ssa fossa filter.vectors prepare.v.filter dim fmask wmask iossa.ssa iossa get.orth.triples svd2LRsvd make.ossa.result save.oblique.decomposition summary.ossa print.ossa summary.iossa.result print.iossa.result owcor frob.dist.series.lowrank traj.matrix.ssa traj.matrix.cssa traj.matrix.mssa traj.matrix.nd.ssa traj.matrix.toeplitz.ssa traj.matrix orthogonalize pseudo.inverse high.rank.rate clone.with.new.series remove.ossa.pssa Cond

Documented in fossa fossa.ssa iossa iossa.ssa owcor print.iossa.result summary.iossa.result

#   R package for Singular Spectrum Analysis
#   Copyright (c) 2014 Alex Shlemov <[email protected]>
#
#   This program is free software; you can redistribute it
#   and/or modify it under the terms of the GNU General Public
#   License as published by the Free Software Foundation;
#   either version 2 of the License, or (at your option)
#   any later version.
#
#   This program is distributed in the hope that it will be
#   useful, but WITHOUT ANY WARRANTY; without even the implied
#   warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
#   PURPOSE.  See the GNU General Public License for more details.
#
#   You should have received a copy of the GNU General Public
#   License along with this program; if not, write to the
#   Free Software Foundation, Inc., 675 Mass Ave, Cambridge,
#   MA 02139, USA.


Cond <- function(A) {
  # Condition number for basis
  d <- svd(A)$d
  d[1] / d[length(d)]
}

# TODO move this functionality to `clone.ssa()`
.remove.ossa.pssa <- function(class) {
  class <- class[(class != "ossa") & !grepl("^pssa", class)]
}

.clone.with.new.series <- function(x, F, ...) {
  # Get conversion
  conversion <- .inner.fmt.conversion(x)

  s <- clone(x, copy.cache = FALSE, copy.storage = FALSE)
  .set(s, "F", conversion(F))
  class(s) <- .remove.ossa.pssa(class(s))

  s
}

high.rank.rate <- function(F, rank, ssaobj = ssa(F), ...) {
  s <- .clone.with.new.series(ssaobj, F)

  # Decompose for computation leading singular values
  .maybe.continue(s, groups = seq_len(rank), ...)

  1 - sum(.sigma(s)[seq_len(rank)]^2) / wnorm(s)^2
}

pseudo.inverse <- function(A) {
  # Moore-Penrose pseudo-inverse
  qrA <- qr(A)
  solve(qr.R(qrA), t(qr.Q(qrA)))
}

orthogonalize <- function(Y, Z, sigma, side = c("bi", "left", "right"), normalize = TRUE) {
  side <- match.arg(side)

  if (missing(sigma)) {
    sigma <- rep(1, ncol(Y))
  }

  rank <- length(sigma)
  # Check parameters' dims
  stopifnot(ncol(Y) == rank, ncol(Z) == rank)

  if (identical(side, "bi")) {
    # Low-rank SVD
    qrY <- qr(Y)
    qrZ <- qr(Z)
    dec <- svd(qr.R(qrY) %*% (sigma * t(qr.R(qrZ))))

    list(d = dec$d, u = qr.Q(qrY) %*% dec$u, v = qr.Q(qrZ) %*% dec$v)
  } else if (identical(side, "left")) {
    qrY <- qr(Y)
    u <- qr.Q(qrY)
    v <- Z %*% (sigma * t(qr.R(qrY)))
    if (normalize) {
      d <- sqrt(colSums(v^2))
      v <- v / matrix(d, nrow = nrow(v), ncol = rank, byrow = TRUE)
    } else {
      d <- rep(1, rank)
    }

    list(d = d, u = u, v = v)
  } else if (identical(side, "right")) {
    dec <- Recall(Y = Z, Z = Y, sigma = sigma, side = "left", normalize = normalize)

    list(d = dec$d, u = dec$v, v = dec$u)
  }
}

# TODO Make emat really polymorphic class. Such wrappings are very ugly
.traj.matrix <- function(x, ...)
  UseMethod(".traj.matrix")

.traj.matrix.1d.ssa <- .traj.matrix.toeplitz.ssa <- function(x, ...)
  .get.or.create.hmat(x, ...)

.traj.matrix.nd.ssa <- function(x, ...)
  .get.or.create.hbhmat(x, ...)

.traj.matrix.mssa <- function(x, ...)
  .get.or.create.mhmat(x, ...)

.traj.matrix.cssa <- function(x, ...)
  .get.or.create.chmat(x, ...)

.traj.matrix.ssa <- function(x, ...)
  stop("`.traj.matrix' is not implemented for this kind of SSA")

.frob.dist.series.lowrank <- function(F, sigma, Y, Z, ssaobj) {
  # This is efficient implementation of |T(F) - Y diag(sigma) Z^T|^2_F
  s <- .clone.with.new.series(ssaobj, F)
  TF <- .traj.matrix(s)

  sZ <- Z * rep(sigma, each = nrow(Z))

  res <- sum((sZ %*% crossprod(Y)) * sZ) + wnorm(s)^2 - 2*sum(crossprod(TF, Y) * sZ)

  # Fix possible numeric error, force dist be nonnegative
  if (res < 0)
    res <- 0

  res
}

.owcor <- function(Fs, LM, RM, ssaobj) {
  mx <- lapply(Fs,
               function(F) {
                 s <- .clone.with.new.series(ssaobj, F)
                 TF <- .traj.matrix(s)
                 as.vector(LM %*% tcrossprod(TF, RM))
              })
  mx <- do.call(cbind, mx)

  # Compute covariations
  cov <- crossprod(mx)

  # Convert to correlations
  cor <- cov2cor(cov)

  # Fix possible numeric error
  cor[cor > 1] <- 1; cor[cor < -1] <- -1

  # Add class
  class(cor) <- "wcor.matrix"

  # Return
  cor
}

print.iossa.result <- function(x, digits = max(3, getOption("digits") - 3), ...) {
  max.abs.nodiag <- function(mx) {
    diag(mx) <- 0
    mx[which.max(abs(mx))]
  }

  cat("\nIterative O-SSA result:\n")
  cat("\tConverged:             ", ifelse(x$converged, "yes", "no"), "\n", sep = "")
  cat("\tIterations:            ", x$iter, "\n", sep = "")
  cat("\tInitial mean(tau):     ", format(mean(x$initial.tau), digits = digits), "\n", sep = "")
  cat("\tInitial tau:           ", paste0(format(x$initial.tau, digits = digits), collapse = ", "), "\n", sep = "")
  cat("\tI. O-SSA mean(tau):    ", format(mean(x$tau), digits = digits), "\n", sep = "")
  cat("\tI. O-SSA tau:          ", paste0(format(x$tau, digits = digits), collapse = ", "), "\n", sep = "")
  cat("\tInitial max wcor:      ", format(max.abs.nodiag(x$initial.wcor), digits = digits), "\n", sep = "")
  cat("\tI. O-SSA max wcor:     ", format(max.abs.nodiag(x$wcor), digits = digits), "\n", sep = "")
  cat("\tI. O-SSA max owcor:    ", format(max.abs.nodiag(x$owcor), digits = digits), "\n", sep = "")
  cat("\n")
}

summary.iossa.result <- function(object, digits = max(3, getOption("digits") - 3), ...)
  print.iossa.result(x = object, digits = digits, ...)

print.ossa <- function(x, digits = max(3, getOption("digits") - 3), ...) {
  NextMethod()

  iossa.result <- .get(x, "iossa.result", allow.null = TRUE)
  if (!is.null(iossa.result))
    print(iossa.result, digits = digits)
}

summary.ossa <- function(object, digits = max(3, getOption("digits") - 3), ...)
  print(x = object, digits = digits, ...)

.save.oblique.decomposition <- function(x, nosigma, Y, Z, idx) {
  sigma <- .sigma(x)
  U <- .U(x)
  V <- if (nv(x) < max(idx)) calc.v(x, seq_len(max(idx))) else .V(x)

  ynorms <- sqrt(colSums(Y^2))
  znorms <- sqrt(colSums(Z^2))

  nosigma <- nosigma * ynorms * znorms
  Y <- Y / rep(ynorms, each = nrow(Y))
  Z <- Z / rep(znorms, each = nrow(Z))

  for (i in seq_along(idx)) {
    sigma[idx[i]] <- nosigma[i]
    U[, idx[i]] <- Y[, i]
    V[, idx[i]] <- Z[, i]
  }

  .set.decomposition(x, sigma = sigma, U = U, V = V)
}

.make.ossa.result <- function(x, Fs, ranks, idx, initial.Fs) {
  IBL <- pseudo.inverse(.U(x)[, idx, drop = FALSE])
  IBR <- pseudo.inverse(calc.v(x, idx))

  tau <- sapply(seq_along(Fs),
                function(i) {
                  high.rank.rate(Fs[[i]], ranks[i], ssaobj = x)
                })

  initial.tau <- sapply(seq_along(initial.Fs),
                        function(i) {
                          high.rank.rate(initial.Fs[[i]], ranks[i], ssaobj = x)
                        })

  names(Fs) <- paste("F", seq_along(Fs), sep = "")
  out <- list(cond = c(Cond(IBL), Cond(IBR)),
              owcor = .owcor(Fs, IBL, IBR, ssaobj = x), #MB just call owcor(x) here?
              wcor = wcor(x, Fs = Fs),
              initial.wcor = wcor(x, Fs = initial.Fs),
              tau = tau,
              initial.tau = initial.tau,
              initial.rec = initial.Fs)
  class(out) <- "iossa.result"

  invisible(out)
}

svd2LRsvd <- function(d, u, v, basis.L, basis.R, need.project = TRUE, fast = TRUE) {
  rank <- length(d)
  # Check parameters' dims
  stopifnot(ncol(u) == rank, ncol(v) == rank)

  ub.L <- if (is.null(basis.L)) {
        basis.L <- u
        diag(nrow = rank, ncol = rank)
      } else {
        crossprod(u, basis.L)
      }

  vb.R <- if (is.null(basis.R)) {
        basis.R <- v
        diag(nrow = rank, ncol = rank)
      } else {
        crossprod(v, basis.R)
      }

  if (fast) {
    dec <- svd(t(ub.L) %*% (vb.R / d))
    dec$d <- 1 / dec$d
    # Reverse order
    dec$d <- rev(dec$d)
    dec$u <- dec$u[, rank:1, drop = FALSE]
    dec$v <- dec$v[, rank:1, drop = FALSE]
  } else {
    dec <- svd(solve(ub.L) %*% (d * t(solve(vb.R))))
  }

  if (need.project) {
    basis.L <- u %*% ub.L
    basis.R <- v %*% vb.R
  }

  list(sigma = dec$d, Y = basis.L %*% dec$u, Z = basis.R %*% dec$v)
}

.get.orth.triples <- function(x, idx, eps = sqrt(.Machine$double.eps),
                              do.orthogonalize = TRUE) {
  # Determine the upper bound of desired eigentriples
  desired <- max(idx)

  # Continue decomposition, if necessary
  if (desired > min(nsigma(x), nu(x)))
    decompose(x, neig = desired)

  sigma <- .sigma(x)[idx]
  U <- .U(x)[, idx, drop = FALSE]
  V <- if (nv(x) < desired) calc.v(x, idx) else .V(x)[, idx, drop = FALSE]

  if (do.orthogonalize) {
    dec <- orthogonalize(U, V, sigma, side = "bi")
    sigma <- dec$d; U <- dec$u; V <- dec$v

    if (min(sigma) < eps)
      warning("Decomposition isn't minimal. Some singular values equal to zero")
  }

  list(sigma = sigma, U = U, V = V)
}

iossa <- function(x, ...)
  UseMethod("iossa")

iossa.ssa <- function(x, nested.groups, ..., tol = 1e-5, kappa = 2,
                      maxiter = 100,
                      norm = function(x) sqrt(mean(x^2)),
                      trace = FALSE,
                      kappa.balance = 0.5) {
  if (!capable(x, "iossa"))
    stop("I-OSSA is not implemented for Complex SSA yet")

  # Get mask
  mask <- .hweights(x) > 0

  if (missing(nested.groups))
    nested.groups <- as.list(1:min(nsigma(x), nu(x)))

  # Continue decomposition, if necessary
  .maybe.continue(x, groups = nested.groups, ...)

  nested.groups <- lapply(nested.groups, unique)
  if (length(unique(unlist(nested.groups))) != sum(sapply(nested.groups, length)))
    stop("Intersected nested groups")

  Fs <- rec <- reconstruct(x, groups = nested.groups) # Get initial approx

  idx <- sort(unique(unlist(nested.groups)))
  triples <- .get.orth.triples(x, idx)
  osigma <- triples$sigma; U <- triples$U; V <- triples$V

  # Replace (in nested.groups) ET's numbers for their ranks (order numbers) in set of signal ETs
  nested.groups <- lapply(nested.groups, function(group) match(group, idx))
  ranks <- sapply(nested.groups, length)

  # If we use reordering, reorder components
  if (!is.null(kappa)) {
    cumranks <- cumsum(ranks)
    for (i in seq_along(nested.groups)) {
      nested.groups[[i]] <- seq(to = cumranks[i], length.out = ranks[i])
    }
  }

  converged <- FALSE
  for (iter in seq_len(maxiter)) {
    lrss <- numeric(length(nested.groups))
    triples.list <- list()
    for (i in seq_along(nested.groups)) {
      cur.ssa <- .clone.with.new.series(x, Fs[[i]])

      .maybe.continue(cur.ssa, seq_len(ranks[i]))
      lrss[i] <- sum(cur.ssa$sigma[-seq_len(ranks[i])]^2)
      triples.list[[i]] <- .get.orth.triples(cur.ssa, seq_len(ranks[i]))
    }

    if (trace) cat(sprintf("LRSS(%d): %s\n", iter, paste0(sqrt(lrss), collapse = " ")))

    osigmas <- lapply(triples.list, function(el) el$sigma)
    Us <- lapply(triples.list, function(el) el$U)
    Vs <- lapply(triples.list, function(el) el$V)

    if (!is.null(kappa)) {
      # If we use reordering, force separability
      for (i in seq_along(nested.groups[-length(nested.groups)])) {
        div <- osigmas[[i]][ranks[i]] / osigmas[[i + 1]][1]
        if (div < kappa) {
          mul <- kappa / div
          osigmas[[i + 1]] <- osigmas[[i + 1]] / mul
          Us[[i + 1]] <- Us[[i + 1]] * mul^kappa.balance
          Vs[[i + 1]] <- Vs[[i + 1]] * mul^(1 - kappa.balance)
        }
      }
    }

    basU <- do.call(cbind, Us)
    basV <- do.call(cbind, Vs)

    # Oblique SVD
    dec <- svd2LRsvd(osigma, U, V, basU, basV, need.project = TRUE)
    sigma <- dec$sigma; Y <- dec$Y; Z <- dec$Z

    # Clone is used for hankelization
    # We copy storage for saving quasi-hankel circulant, it is used for stored fft_plan
    # TODO Separate circulants and plans just like in 1d-SSA
    s <- clone(x, copy.cache = FALSE, copy.storage = TRUE)

    # Add `ossa` class for prevent redecomposition
    class(s) <- c("ossa", class(s)[class(s) != "ossa"])

    .set.decomposition(s, sigma = sigma, U = Y, V = Z)
    Fs.new <- reconstruct(s, groups = nested.groups)

    if (trace) {
      svddist <- sapply(seq_along(nested.groups),
                        function(i) .frob.dist.series.lowrank(Fs[[i]],
                                                              sigma[nested.groups[[i]]],
                                                              Y[, nested.groups[[i]], drop = FALSE],
                                                              Z[, nested.groups[[i]], drop = FALSE],
                                                              x))

      nonhank <- sapply(seq_along(nested.groups),
                        function(i) .frob.dist.series.lowrank(Fs.new[[i]],
                                                              sigma[nested.groups[[i]]],
                                                              Y[, nested.groups[[i]], drop = FALSE],
                                                              Z[, nested.groups[[i]], drop = FALSE],
                                                              x))
      cat(sprintf("SVDD(%d): %s\n", iter, paste0(sqrt(svddist), collapse = " ")))
      cat(sprintf("NHNK(%d): %s\n", iter, paste0(sqrt(nonhank), collapse = " ")))
    }

    deltas <- sapply(seq_along(nested.groups), function(i) .series.dist(Fs.new[[i]], Fs[[i]], norm, mask))
    if (max(deltas) < tol) {
      converged <- TRUE
      Fs <- Fs.new
      break
    }

    Fs <- Fs.new
  }

  x <- clone(x, copy.cache = FALSE) # TODO Maybe preserve relevant part of cache?
  .save.oblique.decomposition(x, sigma, Y, Z, idx)

  # Update class for x
  if (!inherits(x, "ossa")) {
    class(x) <- c("ossa", class(x))
  }

  # Save call info
  x$call <- match.call()

  out <- .make.ossa.result(x, Fs, ranks, idx, rec)
  out[c("iter", "converged", "kappa", "maxiter", "tol", "call")] <-
      list(iter, converged, kappa, maxiter, tol, match.call())

  .set(x, "iossa.result", out)

  # Return to real group numbers
  nested.groups <- lapply(nested.groups, function(group) idx[group])

  # Grab old iossa.groups.all value
  iossa.groups.all <- .get(x, "iossa.groups.all", allow.null = TRUE)
  if (is.null(iossa.groups.all)) {
    iossa.groups.all <- list()
  }

  valid.groups <- as.logical(sapply(iossa.groups.all,
                                    function(group) length(intersect(group, idx)) == 0))
  .set(x, "iossa.groups",  nested.groups)
  .set(x, "iossa.groups.all", c(nested.groups, iossa.groups.all[valid.groups]))

  # Save nested components
  .set(x, "ossa.set", idx)

  if (!is.null(.decomposition(x, "nPR"))) {
    if (any(idx <= .decomposition(x, "nPR"))) {
      .set.decomposition(x, nPR = 0, nPL = 0)
    } else if (any(idx <= sum(unlist(.decomposition(x, c("nPR", "nPL")))))){
      .set.decomposition(x, nPL = 0)
    }
  }

  invisible(x)
}


.wmask <- function(x) {
  wmask <- .get(x, "wmask", NULL)
  if (!is.null(wmask)) {
    return(wmask)
  }

  L <- x$window

  if (inherits(x, "mssa")) {
    L <- c(L, 1)
  }

  array(TRUE, dim = L)
}


.fmask <- function(x) {
  fmask <- .get(x, "fmask", NULL)
  if (!is.null(fmask)) {
    return(fmask)
  }

  N <- x$length
  L <- x$window

  if (inherits(x, "mssa")) {
    n <- length(N)
    N <- N[1]
    K <- ifelse(x$circular, N, N - L + 1)
    K <- c(K, n)
  } else {
    K <- ifelse(x$circular, N, N - L + 1)
  }

  array(TRUE, dim = K)
}


.dim <- function(x) {
  if (inherits(x, "mssa")) {
    2
  } else {
    length(x$length)
  }
}


.prepare.v.filter <- function(x, filter, r) {
  wmask <- !is.na(filter)
  mask <- .fmask(x) #TODO: Cache it
  circular <- x$circular
  if (inherits(x, "mssa")) {
    circular <- c(circular, FALSE)
  }

  fmask <- .factor.mask.2d(mask, wmask, circular = circular)

  if (!all(wmask) || !all(fmask) || any(circular)) {
    weights <- .field.weights.2d(wmask, fmask, circular = circular)

    ommited <- sum(mask & (weights == 0))
    if (ommited > 0) {
      warning(sprintf("Some field elements were not covered by shaped window. %d elements will be ommited", ommited))
    }

    if (all(weights == 0)) {
      warning("Nothing to filter: the given field shape is empty")
    }
  } else {
    weights <- NULL
  }

  .multilayer <- function(a, n) {
    if (is.null(dim(a))) dim(a) <- length(a)

    array(a, dim = c(dim(a), n))
  }

  list(mask = .multilayer(mask, r),
       wmask = .multilayer(wmask, 1),
       fmask = .multilayer(fmask, r),
       weights = if (!is.null(weights)) .multilayer(weights, r) else weights,
       circular = c(circular, FALSE))
}


.filter.vectors <- function(x, vectors, filter) {
  stopifnot(is.array(filter))
  stopifnot(is.matrix(vectors))

  r <- ncol(vectors)

  args <- .prepare.v.filter(x, filter, r)

  F <- array(NA_real_, dim = dim(args$mask))
  F[args$mask] <- as.numeric(vectors)

  w <- array(NA_real_, dim = dim(args$wmask))
  w[args$wmask] <- as.numeric(filter[!is.na(filter)])

  h <- new.hbhmat(F,
                  wmask = args$wmask,
                  fmask = args$fmask,
                  weights = args$weights,
                  circular = args$circular)

  res <- as.vector(w) %*% h

  matrix(res, nrow = length(res) / r, ncol = r)
}


fossa <- function(x, ...)
  UseMethod("fossa")

fossa.ssa <- function(x, nested.groups,
                      filter = c(-1, 1),
                      gamma = Inf,
                      normalize = TRUE,
                      ...) {
  ndim <- .dim(x)
  effndim <- if (inherits(x, "mssa")) 1 else ndim

  if (!is.list(filter)) {
    if (!is.array(filter)) {
      filter <- rep(list(filter), effndim)
    } else {
      filter <- list(filter)
    }
  }

  if (effndim > 1) {
    for (i in seq_along(filter)) {
      if (!is.array(filter[[i]])) {
        d <- rep(1, ndim)
        d[(i - 1) %% ndim + 1] <- length(filter[[i]])
        filter[[i]] <- array(filter[[i]], dim = d)
      }
    }
  } else if (inherits(x, "mssa")) {
    for (i in seq_along(filter)) {
      d <- c(length(filter[[i]]), 1)
      filter[[i]] <- array(filter[[i]], dim = d)
    }
  } else {
    filter <- lapply(filter, as.array)
  }


  if (missing(nested.groups))
    nested.groups <- as.list(1:min(nsigma(x), nu(x)))

  # Continue decomposition, if necessary
  .maybe.continue(x, groups = nested.groups, ...)

  idx <- sort(unique(unlist(nested.groups)))
  triples <- .get.orth.triples(x, idx)
  osigma <- triples$sigma; U <- triples$U; V <- triples$V

  Z <- V * rep(osigma, each = nrow(V))

  Y <- if (normalize) V else Z

  fYs <- lapply(filter, .filter.vectors, x = x, vectors = Y)
  fY <- do.call(rbind, fYs)

  newV <- if (is.infinite(gamma)) fY else rbind(Y, gamma * fY)

  dec <- eigen(crossprod(newV), symmetric = TRUE)

  if (normalize) {
    U <- (U * rep(osigma, each = nrow(U))) %*% dec$vectors
    Z <- V %*% dec$vectors
    sigma <- rep(1, ncol(U))
  } else {
    U <- U %*% dec$vectors
    Z <- Z %*% dec$vectors
    sigma <- rep(1, ncol(U))
  }

  x <- clone(x, copy.cache = FALSE) # TODO Maybe preserve relevant part of cache?
  .save.oblique.decomposition(x, sigma, U, Z, idx)

  if (!is.null(.get(x, "iossa.groups", allow.null = TRUE))) {
    .set(x, "iossa.groups", .fix.iossa.groups(.get(x, "iossa.groups", allow.null = TRUE), idx))
    .set(x, "iossa.groups.all", .fix.iossa.groups(.get(x, "iossa.groups.all", allow.null = TRUE), idx))
  }

  # Save nested components
  .set(x, "ossa.set", idx)

  if (!is.null(.decomposition(x, "nPR"))) {
    if (any(idx <= .decomposition(x, "nPR"))) {
      .set.decomposition(x, nPR = 0, nPL = 0)
    } else if (any(idx <= sum(unlist(.decomposition(x, c("nPR", "nPL")))))){
      .set.decomposition(x, nPL = 0)
    }
  }

  if (!inherits(x, "ossa")) {
    class(x) <- c("ossa", class(x))
  }

  # Save call info
  x$call <- match.call()

  invisible(x)
}

.fix.iossa.groups <- function(iossa.groups, group) {
  if (length(iossa.groups) == 0) {
    return(list())
  }

  touched <- sapply(iossa.groups, function(g) length(intersect(g, group)) > 0)

  if (any(touched)) {
    iossa.groups <- c(iossa.groups[!touched], union(group, unlist(iossa.groups[touched])))
  }

  iossa.groups
}

decompose.ossa <- function(x, ...) {
  # We can simply implement continuation if we store reference to initial decomposition
  stop("Continuation of decomposition is impossible for ObliqueSSA")
}

.colspan.ossa <- function(x, idx) {
  qr.Q(qr(.U(x)[, idx, drop = FALSE]))
}

.rowspan.ossa <- function(x, idx) {
  qr.Q(qr(calc.v(x, idx)))
}

owcor <- function(x, groups, ..., cache = TRUE) {
  # Check class
  stopifnot(inherits(x, "ossa"))

  # Get basis
  basis <- .get(x, "ossa.set")

  if (missing(groups)) {
    groups <- x$iossa.groups
    if (is.null(groups)) {
      groups <- as.list(basis)
    }
  }

  if (!all(unlist(groups) %in% basis))
    stop(sprintf("groups in `groups' must consist of components from current OSSA set (%s)",
                 paste0(basis, collapse = ", ")))

  # Compute reconstruction.
  Fs <- reconstruct(x, groups, ..., cache = cache)

  LM <- pseudo.inverse(.U(x)[, basis, drop = FALSE]) # No colspan here!!!!
  RM <- pseudo.inverse(calc.v(x, basis)) # No rowspan here!!!!

  # Compute oblique w-correlations and return
  res <- .owcor(Fs, LM, RM, ssaobj = x)
  colnames(res) <- rownames(res) <- names(Fs)

  res
}
asl/rssa documentation built on July 6, 2018, 10:06 a.m.