R/multiple_align_functions.R

In fdasrvf: Elastic Functional Data Analysis

#' Group-wise function alignment to specified mean
#'
#' This function aligns a collection of functions using the elastic square-root
#' slope (srsf) framework.
#'
#' @param f matrix (\eqn{N} x \eqn{M}) of \eqn{M} functions with \eqn{N} samples
#' @param time vector of size \eqn{N} describing the sample points
#' @param mu vector of size \eqn{N} that f is aligned to
#' @param lambda controls the elasticity (default = 0)
#' @param pen alignment penalty (default="roughness") options are
#' second derivative ("roughness"), geodesic distance from id ("geodesic"), and
#' norm from id ("norm")
#' @param showplot shows plots of functions (default = T)
#' @param smooth_data smooth data using box filter (default = F)
#' @param sparam number of times to apply box filter (default = 25)
#' @param parallel enable parallel mode using `foreach` and
#'   `doParallel` package (default=F)
#' @param omethod optimization method (DP,DP2,RBFGS,dBayes,expBayes)
#' @param MaxItr maximum number of iterations
#' @param iter bayesian number of mcmc samples (default 2000)
#' @param verbose verbose printing (default TRUE)
#' @return Returns a fdawarp object containing \item{f0}{original functions}
#' \item{fn}{aligned functions - matrix (\eqn{N} x \eqn{M}) of \eqn{M} functions with \eqn{N} samples}
#' \item{qn}{aligned SRSFs - similar structure to fn}
#' \item{q0}{original SRSF - similar structure to fn}
#' \item{fmean}{function mean or median - vector of length \eqn{N}}
#' \item{mqn}{SRSF mean or median - vector of length \eqn{N}}
#' \item{gam}{warping functions - similar structure to fn}
#' \item{orig.var}{Original Variance of Functions}
#' \item{amp.var}{Amplitude Variance}
#' \item{phase.var}{Phase Variance}
#' \item{qun}{Cost Function Value}
#' @keywords srsf alignment
#' @references Srivastava, A., Wu, W., Kurtek, S., Klassen, E., Marron, J. S.,
#'  May 2011. Registration of functional data using fisher-rao metric,
#'  arXiv:1103.3817v2.
#' @references Tucker, J. D., Wu, W., Srivastava, A.,
#'  Generative Models for Function Data using Phase and Amplitude Separation,
#'  Computational Statistics and Data Analysis (2012), 10.1016/j.csda.2012.12.001.
#' @export
multiple_align_functions <- function(f,
                                     time,
                                     mu,
                                     lambda = 0,
                                     pen = "roughness",
                                     showplot = TRUE,
                                     smooth_data = FALSE,
                                     sparam = 25,
                                     parallel = FALSE,
                                     omethod = "DP",
                                     MaxItr = 20,
                                     iter = 2000,
                                     verbose = TRUE) {
  if (parallel) {
    cores = max(parallel::detectCores() - 1, 1)
    cl = parallel::makeCluster(cores)
    doParallel::registerDoParallel(cl)
  } else
  {
    foreach::registerDoSEQ()
  }

  if (verbose){
    cat(sprintf("lambda = %5.1f \n", lambda))
  }

  binsize = mean(diff(time))
  eps = .Machine$double.eps
  M = nrow(f)
  N = ncol(f)
  f0 = f
  w = 0.0

  if (smooth_data) {
    f = smooth.data(f, sparam)
  }

  if (showplot) {
    graphics::matplot(time, f, type = "l")
    graphics::title(main = "Original data")
  }

  # Compute q-function of the functional data
  tmp = gradient.spline(f, binsize, smooth_data)
  f = tmp$f
  q = tmp$g / sqrt(abs(tmp$g) + eps)

  tmp = gradient.spline(mu, binsize, smooth_data)
  mf = tmp$f
  mq = tmp$g / sqrt(abs(tmp$g) + eps)
  k <- 1

  if (verbose){
    cat(sprintf("Aligning %d functions in SRSF space...\n", N))
  }
  outfor <- foreach::foreach(k = 1:N,
                             .combine = cbind,
                             .packages = 'fdasrvf') %dopar% {
                               if (omethod == "expBayes") {
                                 gam <- pair_align_functions_expomap(mu, c(f[, k]), time, iter = iter)$gamma
                                 gam <- gam$y
                               } else if (omethod == "dBayes") {
                                 gam <- pair_align_functions_bayes(mu, f[, k], time)$gam_a
                               } else {
                                 gam <- optimum.reparam(mq, time, q[, k], time, lambda, pen, omethod, mf[1], f[1, k])
                               }

                               gam_dev = gradient(gam, 1 / (M - 1))
                               f_temp = stats::approx(time, f[, k], xout = (time[length(time)] - time[1]) *
                                                        gam +
                                                        time[1])$y
                               q_temp = f_to_srvf(f_temp, time)
                               v <- q_temp - mq
                               d <- sqrt(trapz(time, v * v))
                               vtil <- v / d
                               dtil <- 1 / d

                               list(gam, gam_dev, q_temp, f_temp, vtil, dtil)
                             }

  gam = unlist(outfor[1, ])
  dim(gam) = c(M, N)
  gam = t(gam)
  gam_dev = unlist(outfor[2, ])
  dim(gam_dev) = c(M, N)
  gam_dev = t(gam_dev)
  q_temp = unlist(outfor[3, ])
  dim(q_temp) = c(M, N)
  f_temp = unlist(outfor[4, ])
  dim(f_temp) = c(M, N)
  qn = q_temp
  fn = f_temp
  tmp = (1 - sqrt(gam_dev)) ^ 2
  vtil = unlist(outfor[5, ])

  dim(vtil) = c(M, N)
  dtil = unlist(outfor[6, ])

  dim(dtil) = c(1, N)



  # Aligned data & stats
  q0 = q
  mean_f0 = rowMeans(f)
  std_f0 = apply(f, 1, stats::sd)
  mean_fn = rowMeans(fn)
  std_fn = apply(fn, 1, stats::sd)
  mqn = mq
  fmean = mean(f0[1, ]) + cumtrapz(time, mqn * abs(mqn))
  gam = t(gam)
  gamI = SqrtMeanInverse(gam)


  fgam = matrix(0, M, N)
  for (ii in 1:N) {
    fgam[, ii] = stats::approx(time, fmean, xout = (time[length(time)] - time[1]) *
                                 gam[, ii] +
                                 time[1])$y
  }
  var_fgam = apply(fgam, 1, stats::var)

  orig.var = trapz(time, std_f0 ^ 2)
  amp.var = trapz(time, std_fn ^ 2)
  phase.var = trapz(time, var_fgam)

  out <- list(
    f0 = f,
    time = time,
    fn = fn,
    qn = qn,
    q0 = q0,
    fmean = fmean,
    mqn = mqn,
    warping_functions = gam,
    original_variance = orig.var,
    amplitude_variance = amp.var,
    phase_variance = phase.var,
    qun = 0,
    inverse_average_warping_function = gamI,
    rsamps = FALSE,
    call = list(
      lambda = lambda,
      penalty_method = pen,
      centroid_type = "mean",
      center_warpings = FALSE,
      smooth_data = smooth_data,
      sparam = sparam,
      parallel = parallel,
      optim_method = omethod,
      max_iter = MaxItr
    )
  )

  class(out) <- 'fdawarp'

  if (showplot) {
    plot(out)
  }

  if (parallel)
    parallel::stopCluster(cl)

  return(out)
}