# R/backtransformAffine.R In HenrikBengtsson/aroma.light: Light-Weight Methods for Normalization and Visualization of Microarray Data using Only Basic R Data Types

#########################################################################/**
# @RdocGeneric backtransformAffine
# @alias backtransformAffine.matrix
#
# @title "Reverse affine transformation"
#
# \description{
#   @get "title".
# }
#
# \usage{
# @usage backtransformAffine,matrix
# }
#
# \arguments{
#  \item{X}{An NxK @matrix containing data to be backtransformed.}
#  \item{a}{A scalar, @vector, a @matrix, or a @list.
#    First, if a @list, it is assumed to contained the elements \code{a}
#    and \code{b}, which are the used as if they were passed as separate
#    arguments.
#    If a @vector, a matrix of size NxK is created which is then filled
#    \emph{row by row} with the values in the vector. Commonly, the
#    vector is of length K, which means that the matrix will consist of
#    copies of this vector stacked on top of each other.
#    If a @matrix, a matrix of size NxK is created which is then filled
#    \emph{column by column} with the values in the matrix (collected
#    column by column. Commonly, the matrix is of size NxK, or NxL with
#    L < K and then the resulting matrix consists of copies sitting
#    next to each other.
#    The resulting NxK matrix is subtracted from the NxK matrix \code{X}.
#  }
#  \item{b}{A scalar, @vector, a @matrix.
#    A NxK matrix is created from this argument. For details see
#    argument \code{a}.
#    The NxK matrix \code{X-a} is divided by the resulting NxK matrix.
#  }
#  \item{project}{
#    returned (K values per data point are returned).
#    If @TRUE, the backtransformed values "\code{(X-a)/b}" are projected
#    onto the line L(a,b) so that all columns
#    will be identical.
#  }
#  \item{...}{Not used.}
# }
#
# \value{
#   The "\code{(X-a)/b}" backtransformed NxK @matrix is returned.
#   If \code{project} is @TRUE, an Nx1 @matrix is returned, because
#   all columns are identical anyway.
# }
#
# \section{Missing values}{
#   Missing values remain missing values. If projected, data points that
#   contain missing values are projected without these.
# }
#
# @examples "../incl/backtransformAffine.matrix.Rex"
#
#*/#########################################################################
setMethodS3("backtransformAffine", "matrix", function(X, a=NULL, b=NULL, project=FALSE, ...) {

# Dimensions of 'X'
nobs  <- nrow(X);
ndims <- ncol(X);
if (ndims == 1L) {
stop("Can not fit affine multiscan model. Matrix must contain at least two columns (scans): ", ndims);
}

# If argument 'a' is a list assume it contains the elements 'a' and 'b'.
if (is.list(a)) {
b <- a$b; a <- a$a;
}

# If 'a' and/or 'b' are vector convert them to row matrices.
if (is.vector(a)) {
# Create a full matrix and filled row by row with 'a'
a <- matrix(a, nrow=nobs, ncol=ndims, byrow=TRUE);
} else if (is.matrix(a)) {
# Create a full matrix and filled column by column by the columns in 'a'
t <- a;
a <- matrix(NA_real_, nrow=nobs, ncol=ndims);
for (cc in 1:ndims) {
# Loop over the columns in a0 too.
col <- ((cc-1) %% ncol(t)) + 1L;
value <- rep(t[,col], length.out=nobs);
a[,cc] <- value;
}
# Not needed anymore
t <- NULL;
} else if (!is.null(a)) {
stop(paste("Unknown data type of argument 'a':", class(a)[1]));
}

if (!project) {
if (is.vector(b)) {
# Create a full matrix and filled row by row with 'b'
b <- matrix(b, nrow=nobs, ncol=ndims, byrow=TRUE);
} else if (is.matrix(b)) {
# Create a full matrix and filled column by column by the columns in 'b'
t <- b;
b <- matrix(NA_real_, nrow=nobs, ncol=ndims);
for (cc in 1:ndims) {
# Loop over the columns in a0 too.
col <- ((cc-1) %% ncol(t)) + 1L;
value <- rep(t[,col], length.out=nobs);
b[,cc] <- value;
}
} else if (!is.null(b)) {
stop(paste("Unknown data type of argument 'b':", class(b)[1]));
}
}

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# 2. Subtract the bias and rescale
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if (!is.null(a))
X <- X - a;

########################################################################
# Alternative 2
#
#   i) Translate the fitted line L to L' such that L' goes through
#      (a,a,...,a) and project the data y onto L' to obtain ytilde
#  ii) from which we can calculate xtilde = (ytilde - a) / b
# iii) Since all xtilde are the same take the first component to be our
#      estimate xhat in y = a + b*xhat.
########################################################################
if (project) {
# In theory:
#  ytilde <- projectUontoV(y-a,b) + a;
#  xtilde <- (ytilde-a)/b;
# In practice:
X <- t(X);

# 'b' standardized to a unit vector such that <v,v> == 1.
v <- b / sqrt(sum(b^2));

# projectUontoV(): U should be an NxK matrix and v an N vector.
X <- projectUontoV(X,v, na.rm=TRUE);

X <- X[1,] / b[1];                     # Note that here 'b' is a vector!
} else {
if (!is.null(b))
X <- X / b;                          # Note that here 'b' is a matrix!
}

as.matrix(X);
}) # backtransformAffine()

############################################################################
# HISTORY:
# 2011-04-12
# o Now using as.double(NA) instead of NA.
# 2006-06-03
# o Minor to merge two different threads of this code.
# o Method passes the tests in the example code (again).
# 2005-01-24
# o Now missing values are excluded before projection.
# 2005-01-08
# o Now, if project is TRUE, only one column is returned.
# o Now the method is guaranteed to return a matrix by calling as.matrix().
# o The average of projected data is now no the same scale as the average on
#   non-projected data. Before the data was max(b) times too large.