R/General.R

In manifold: Operations for Riemannian Manifolds

## S3 methods

#' Returns the Riemannian metric
#'
#' @param mfd A manifold object created by \code{\link{createM}}
#' @param p A vector containing the (single) base point on the manifold.
#' @param U,V Matrices with n columns each containing a tangent vector. The metric of each pair of n columns is calculated. 
#' @returns A vector with n entries containing the metric values
#'
#' @export
metric <- function(mfd, p, U, V) UseMethod("metric", mfd)

#' Norm on the tangent space induced by the Riemannian metric
#' @inheritParams metric
#' @param p,U Matrices with n columns. p stands for the base point(s) and U the tangent vector(s). The norm of each column in U is calculated w.r.t. the corresponding column in p. 
#' @returns A vector with n entries containing the norms of the tangent vectors
#' @export
norm <- function(mfd, p, U) UseMethod("norm", mfd)


#' Geodesic distance
#' @inheritParams metric
#' @param X,Y Matrices with n columns. The distance between each pair of columns is calculated. If either X or Y is a vector then it is recycled.
#' @param ... Passed into specific methods
#' @returns A vector with n entries containing the distances between pairs of points
#' @export
distance <- function(mfd, X, Y, ...) UseMethod("distance", mfd)


#' Obtain the geodesic curve
#' @param mfd A manifold object created by \code{\link{createM}}
#' @export
#' @returns A matrix of points lying on a geodesic. Different columns correspond to different time points
geodesicCurve <- function(mfd, p, h, t) UseMethod("geodesicCurve", mfd)


#' Riemannian exponential map
#' @param mfd A manifold object created by \code{\link{createM}}
#' @param p,V Matrices with n columns. The exponential map of each column in V is calculated w.r.t. the corresponding column in p. 
#' @param ... Passed into specific methods
#' @export
#' @returns A matrix in which columns correspond to points on the manifold
rieExp <- function(mfd, p, V, ...) UseMethod("rieExp", mfd)


#' Riemannian logarithm map
#' @param mfd A manifold object created by \code{\link{createM}}
#' @param p,X Matrices with n columns. The logarithm map of each column in X is calculated w.r.t. the corresponding column in p. 
#' @param ... Passed into specific methods
#' @export
#' @returns A matrix in which columns correspond to points on the tangent spaces
rieLog <- function(mfd, p, X, ...) UseMethod("rieLog", mfd)


#' Project data points in the ambient space onto the manifold
#' @param mfd A manifold object created by \code{\link{createM}}
#' @param p A matrix of points to be projected to the manifold.
#' @export
#' @returns A matrix in which columns correspond to points on the manifold
project <- function(mfd, p) UseMethod("project", mfd)


# # TODO: 
# geopolate <- function(mfd, ...) UseMethod("geopolate", mfd)


#' Project data points in the ambient space onto the tangent space
#'
#' @param mfd A manifold object created by \code{\link{createM}}
#' @param p A vector containing the base point on the manifold. Data X will be projected onto the tangent space of p.
#' @param X A vector or a matrix containing the data points in terms of the coordinates in the ambient space.
#' @param projMatOnly Whether to only return the projection matrix (FALSE) or the projected data points (TRUE, the default)
#' @param ... Passed into specific methods
#' @export
#' @returns A matrix in which columns correspond to points on the tangent spaces
projectTangent <- function(mfd, p, X, projMatOnly, ...) UseMethod('projectTangent', mfd)


#' Calculate the Fréchet median
#' @param mfd A manifold object created by \code{\link{createM}}
#' @param X A matrix with each column being a point on the manifold.
#' @param ... Passed into specific methods
#' @export
#' @returns A matrix with 1 column containing the Fréchet median
frechetMedian <- function(mfd, X, ...) UseMethod("frechetMedian", mfd)


#' Calculate the Fréchet mean
#' @param mfd A manifold object created by \code{\link{createM}}
#' @param X A matrix with each column being a point on the manifold.
#' @param ... Passed into specific methods
#' @export
#' @returns A matrix with 1 column containing the Fréchet mean
frechetMean <- function(mfd, X, ...) UseMethod("frechetMean", mfd)


#' Generate random variables on the manifold
#'
#' By default, random variables are generated by mapping isotropic Gaussian distributions on the tangent space back to the manifold using the exponential map
#' @param mfd A manifold object created by \code{\link{createM}}
#' @param n Sample size
#' @param dimIntrinsic The intrinsic dimension of the target manifold
#' @param ... Passed into specific methods
#' @export
#' @returns A matrix with n columns, each containing a random sample on the manifold
rmfd <- function(mfd, n, dimIntrinsic, ...) UseMethod("rmfd", mfd)


#' Returns the origin of the manifold
#' 
#' Each manifold defined in this package has been designated a more or less reasonable choice of origin.
#' @param mfd A manifold object created by \code{\link{createM}}
#' @param dimIntrinsic An integer for the intrinsic dimension
#' @param ... Passed into specific methods
#' @export
#' @returns A matrix with 1 column, a designation of the origin on the manifold
origin <- function(mfd, dimIntrinsic, ...) UseMethod('origin', mfd)


#' Transform the coordinates for the tangent space into a tangent vector
#'
#' Transform the coordinates w.r.t. a basis returned by [basisTan()] to tangent vectors
#' @param mfd A manifold object created by \code{\link{createM}}
#' @param p A vector for the base point 
#' @param V0 A matrix or a vector containing the coordinates
#' @param ... Other arguments passed into specific methods
#' @export
#' @returns A matrix in which each column correspond to a different tangent vector
coordToTanV <- function(mfd, p, V0, ...) {
  if (!missing(p)) {
    if (length(p) == 0) {
      stop('`p` cannot have length 0')
    }
    if (is.matrix(p) && ncol(p) > 1) {
      stop('`p` must be a single point')
    }
  }
  if (length(V0) == 0) {
    stop('`V` cannot have length 0')
  }
  UseMethod('coordToTanV', mfd)
}


#' Transform tangent vectors to their coordinates on the tangent space
#'
#' Transform the tangent vectors to coordinates w.r.t. the basis returned by [basisTan()]
#' @param mfd A manifold object created by \code{\link{createM}}
#' @param p A vector for the base point 
#' @param V A matrix or a vector containing the tangent vectors
#' @param ... Passed into specific methods
#' @export
#' @returns A matrix in which each column correspond to a coordinate
tanVToCoord <- function(mfd, p, V, ...) {

  if (!missing(p)) {
    if (length(p) == 0) {
      stop('`p` cannot have length 0')
    }
    if (is.matrix(p) && ncol(p) > 1) {
      stop('`p` must be a single point')
    }
  }

  if (length(V) == 0) {
    stop('`V` cannot have length 0')
  }
  UseMethod('tanVToCoord', mfd)
}


#' Obtain an orthonormal basis on the tangent space
#'
#' Parametrize the tangent space at location p, so that the parameterized version contains an open neighborhood around the origin. (The dimension of v is potentially reduced). 
#' @param mfd A manifold object created by \code{\link{createM}}
#' @param p A vector for a base point on the manifold
#' @returns An orthonormal basis matrix D, whose columns contain the basis vectors, so that `D^T v` give the coordinates `v0` for a tangent vector `v`, and `D%*%v0 = v`.
#' @export
basisTan <- function(mfd, p) UseMethod('basisTan')


# dotprod <- function(u, v) {
    # return(colSums(u*v))
# }


#' Dimensions in this package
#' 
#' `calcGeomPar` calculates geometric parameter.
#' `calcIntDim` calculates the intrinsic dimensions.
#' `calcAmbDim` calculates the ambient dimensions.
#' `calcTanDim` calculates the number of tuples used to represent a tangent vector.
#' @param mfd A manifold object created by \code{\link{createM}}
#' @param dimAmbient Dimension of the ambient space.
#' @param dimIntrinsic Intrinsic dimension of the manifold.
#' @param dimTangent The length of a tangent vector
#' @param geomPar The geometric parameter (e.g., SO(*))
#' @name dimensions
NULL
#> NULL


#' @export
#' @rdname dimensions
#' @returns A scalar value for the geometric parameter
calcGeomPar <- function(mfd, dimAmbient, dimIntrinsic, dimTangent) {
  if ((!missing(dimAmbient)) + 
      (!missing(dimIntrinsic)) + 
      (!missing(dimTangent)) > 1) {
    stop('Can only specify one of dimIntrinsic, dimAmbient, and dimTangent')
  }
  UseMethod('calcGeomPar', mfd)
}


#' @export
#' @rdname dimensions
#' @returns A scalar value for the intrinsic dimension
calcIntDim <- function(mfd, geomPar, dimAmbient, dimTangent) {
  if ((!missing(geomPar)) + 
      (!missing(dimAmbient)) + 
      (!missing(dimTangent)) > 1) {
    stop('Cannot specify two or more of the dimensions')
  }
  UseMethod('calcIntDim', mfd)
}


#' @export
#' @rdname dimensions
#' @returns A scalar value for the ambient dimension
calcAmbDim <- function(mfd, geomPar, dimIntrinsic, dimTangent) {
  if ((!missing(geomPar)) + 
      (!missing(dimIntrinsic)) + 
      (!missing(dimTangent)) > 1) {
    stop('Cannot specify two or more of the dimensions')
  }
  UseMethod('calcAmbDim', mfd)
}


#' @export
#' @rdname dimensions
#' @returns A scalar value for the number of components in the (implementation of the) tangent vector
calcTanDim <- function(mfd, geomPar, dimAmbient, dimIntrinsic) {
  if ((!missing(geomPar)) + 
      (!missing(dimAmbient)) + 
      (!missing(dimIntrinsic)) > 1) {
    stop('Cannot specify two or more of the dimensions')
  }
  UseMethod('calcTanDim', mfd)
}


#' @describeIn metric Method
#' @export
metric.default <- function(mfd, p=NULL, U, V) {

  U <- as.matrix(U)
  V <- as.matrix(V)

  if ((ncol(U) <= 1 && ncol(V) >= 1) || (ncol(U) >= 1 && ncol(V) <= 1)) {
    c(crossprod(U, V))
  } else {
    colSums(U * V)
  }
}


#' @describeIn norm Method
#' @export
norm.default <- function(mfd, p=NULL, U) {
  U <- as.matrix(U)
  if ((!missing(p)) && (!is.null(p)) && (is.matrix(p) && ncol(p) > 1) && (ncol(U) == 1)) {
    U <- matrix(U, nrow(U), ncol(p))
  }
  sqrt(colSums(U * U))
}


#' @param p Optionally, specify the base point of the tangent space, on which random tangent vectors will be generated. Default to the origin of `mfd`. 
#' @param dist Either a string or a function, describing the distributions. If it is character, then independent univariate r.v. following that distribution with total variance 1 is generated on the tangent space. If it is a function, it should specify the generation of the samples on the tangent space _coordinates_ (with dimension dimIntrinsic), and it must return a matrix for which the columns are the coordinates. Then it will be mapped to the tangent space at the origin and then to the manifold.
#' @param totalVar Total variance (sum of diagonal of the covariance matrix) of the tangent random vector
#' @export
#' @describeIn rmfd Default method
rmfd.default <- function(mfd, n, dimIntrinsic, p, dist=c('norm', 'unif', 'exp'), totalVar=1, ...) {

  if (!is.function(dist)) {
    dist <- match.arg(dist)
  }

  if (n == 0) {
    return(matrix(0, nrow=calcAmbDim(mfd, dimIntrinsic=dimIntrinsic), ncol=0))
  }

  if (dimIntrinsic == 0) {
    stop('`dimIntrinsic` must be at least one')
  }

  if (is.finiteDim(mfd)) {
    if (is.function(dist)) {
      V0 <- dist(n, dimIntrinsic, ...)
    } else {
      # Function for generating univariate r.v. with variance 1
      if (dist == 'norm') {
        runi <- function(x) stats::rnorm(x)
      } else if (dist == 'unif') {
        runi <- function(x) stats::runif(x, min=-sqrt(3), max=sqrt(3))
      } else if (dist == 'exp') {
        runi <- function(x) stats::rexp(x) - 1
      }

      V0 <- matrix(runi(n * dimIntrinsic), dimIntrinsic, n) * sqrt(totalVar / dimIntrinsic)

    } 
  
    if (missing(p)) {
      p <- origin(mfd, dimIntrinsic=dimIntrinsic)
    }
    V <- coordToTanV(mfd, p, V0)
  } else { # For functional manifolds
    stop('to be implemented')
  }

  rieExp(mfd, p, V)
}


#' @param ... Passed into specific methods
#' @export
#' @describeIn coordToTanV Method
coordToTanV.default <- function(mfd, p, V0, ...) {
  D <- basisTan(mfd, p)
  D %*% as.matrix(V0)
}


#' @param ... Passed into specific methods
#' @export
#' @describeIn tanVToCoord Method
tanVToCoord.default <- function(mfd, p, V, ...) {
  D <- basisTan(mfd, p)
  t(D) %*% as.matrix(V)
}


#' List all available manifold names
#' 
#' @returns A character vector of the names of available manifolds
#' @export
listAvailMfd <- function() {
  gsub('^.*\\.', '', as.character(utils::methods(distance)))
}


#' Tell whether a manifold is a finite-dimensional one
#' @param mfd A manifold object created by \code{\link{createM}}
#' @returns A logical scalar indicating whether the manifold is finite-dimensional
#' @export
is.finiteDim <- function(mfd) {
  if (!is.character(mfd)) {
    mfd <- class(mfd)
  }
  isFunctional <- mfd %in% c('HS', 'L2', 'Dens') 
  all(!isFunctional)
}

#' Create an object to symbolize the manifold. 
#'
#' Create an object to symbolize the manifold. Supports the Euclidean, the unit sphere, special orthogonal group, and symmetric positive definite (SPD) matrices with either the affine invariant or the log-Euclidean metric.
#'
#' Overall, most of the arguments in the manifold functions respect recycling rules. If the input arguments are matrices with n rows corresponding to n points, then the manifold function will be applied to each of the n pairs; the returned value will be a matrix. If one input corresponds to n points and the other just one point, then the one point will be recycled. If each component corresponds to one point, the returned value will correspond to just one point.
#' @param mfdName One of 'Euclidean', 'Sphere', 'SO', 'AffInv', or 'LogEu'. The name is case sensitive.
#' @returns A structure 1 with class being `mfdName`
#' @export
createM <- function(mfdName) {

  if (mfdName %in% c('AffInv', 'LogEu') ) {
    parent <- 'SPD'
  } else {
    parent <- character(0)
  }

  structure(1, class=c(mfdName, parent))
}


#' @param mu0 A matrix of starting points. Each column corresponds to a starting point. If there are multiple columns, then a multistart algorithm is used. 
#' @param weight A vector of weights for the observations
#' @param tol The threshold for determining convergence
#' @param maxit Maximum iteration
#' @param alpha The step size parameter. See Fletcher et al. (2008) CVPR
#' @export
#' @describeIn frechetMedian Default method
frechetMedian.default <- function(mfd, X, mu0, weight=NULL, tol=1e-9, maxit=1000, alpha=1, ...) {

  X <- as.matrix(X)
  dimAmbient <- nrow(X)
  dimTangent <- calcTanDim(mfd, dimAmbient=dimAmbient)
  
  if (is.null(weight)) {
    n <- dim(X)[2]
    weight <- rep(1 / n, n)
    if (missing(mu0)) {
      mu0 <- apply(X, 1, mean)
    }
  } else {
    weight <- weight / sum(weight)
    ind <- abs(weight) > 1e-14
    n <- sum(ind)
    X <- X[, ind, drop=FALSE]
    weight <- weight[ind]
    if (missing(mu0)) {
      mu0 <- apply(X*t(replicate(dimAmbient, weight)), 1, sum)
    }
  }

  mu0 <- as.matrix(mu0)
  nStart <- ncol(mu0)

  allRes <- lapply(seq_len(nStart), function(i) {

    mu00 <- mu0[, i]
    it <- 0
    dif <- Inf
    SS <- 0

    while(dif > tol && it < maxit) {
      mu00 <- project(mfd, mu00)
      it <- it + 1
      V <- rieLog(mfd, mu00, X)
      # The Weiszfeld Algorithm. c.f. Fletcher et al (2008)
      D <- norm(mfd, mu00, V)
      indPos <- D >= 1e-10
      nPos <- sum(indPos)

      if (nPos > 0) {
        W <- sum(weight[indPos] / D[indPos])
        sumV <- rowSums(matrix(weight[indPos] / D[indPos], 
                               dimTangent, 
                               sum(indPos), byrow=TRUE) * 
                        V[, indPos, drop=FALSE])
        vk <- sumV / W
      } else {
        vk <- rep(0, dimTangent)
      }

      muNew <- as.numeric(rieExp(mfd, mu00, vk * alpha))
      SSNew <- sum(D[indPos])
      dif <- abs(SS - SSNew) / (SS + tol)

      mu00 <- muNew
      SS <- SSNew
    }

    list(x = mu00, fx = SS, it = it)
  })

  minInd <- which.min(vapply(allRes, `[[`, 1, 'fx'))
  mu0 <- allRes[[minInd]]$x
  it <- allRes[[minInd]]$it

  if (it == maxit) {
    warning('Maximum iteration reached!')
  }

  matrix(mu0)
}


#' @param mu0 A matrix of starting points. Each column corresponds to a starting point. If there are multiple columns, then a multistart algorithm is used. 
#' @param weight A vector of weights for the observations
#' @param tol The threshold for determining convergence
#' @param maxit Maximum iteration
#' @param lam The step length
#' @export
#' @describeIn frechetMean Default method
frechetMean.default <- function(mfd, X, mu0, weight=NULL, tol=1e-9, maxit=1000, lam=0.99, ...) {

  X <- as.matrix(X)
  dimAmbient <- nrow(X)
  dimTangent <- calcTanDim(mfd, dimAmbient=dimAmbient)
  
  if (is.null(weight)) {
    n <- dim(X)[2]
    weight <- rep(1 / n, n)
    if (missing(mu0)) {
      mu0 <- apply(X, 1, mean)
    }
  } else {
    weight <- weight / sum(weight)
    ind <- abs(weight) > 1e-14
    n <- sum(ind)
    X <- X[, ind, drop=FALSE]
    weight <- weight[ind]
    if (missing(mu0)) {
      mu0 <- apply(X*t(replicate(dimAmbient, weight)), 1, sum)
    }
  }

  mu0 <- as.matrix(mu0)
  nStart <- ncol(mu0)

  allRes <- lapply(seq_len(nStart), function(i) {

    mu00 <- mu0[, i]
    it <- 0
    dif <- Inf
    SS <- 0

    while(dif > tol && it < maxit) {
      mu00 <- project(mfd, mu00)
      it <- it + 1
      V <- rieLog(mfd, mu00, X)
      vNew <- as.numeric(rowSums(V * matrix(weight, dimTangent, n, byrow=TRUE)))
      muNew <- as.numeric(rieExp(mfd, mu00, vNew * lam^it))
      SSNew <- sum(V^2)
      dif <- abs(SS - SSNew) / (SS + tol)

      mu00 <- muNew
      SS <- SSNew
    }

    list(x = mu00, fx = SS, it = it)
  })

  minInd <- which.min(vapply(allRes, `[[`, 1, 'fx'))
  mu0 <- allRes[[minInd]]$x
  it <- allRes[[minInd]]$it

  if (it == maxit) {
    warning('Maximum iteration reached!')
  }

  matrix(mu0)
}