R/plotreg.fd.R

In fda: Functional Data Analysis

plotreg.fd <- function(reglist) {
#  plots the output of registration function register.fd
#
#  Argument:
#  REGLIST ... The named list generated by function REGISTER.FD
#              The members of REGLIST that are required are:
#              REGFD  ... the registered functions
#              WFD    ... the functions W(t) defining the warping functions h(t)
#              YFD    ... the unregistered functions
#              Y0FD   ... the target functions
#  If required objects are missing, REGLIST was probably generated by
#  an older verson of REGISTERFD, and the registration should be redone.

#  Last modified  17 December 2010 by Jim Ramsay

#  check the argument REGLIST

if (!inherits(reglist, "list")) stop("REGLIST is not a list object.")

#  extract the required members of REGLIST

missinginfo <- FALSE
if (is.null(reglist$regfd)) {
    missinginfo <- TRUE
} else {
   yregfd  <- reglist$regfd
}
if (is.null(reglist$Wfd)) {
    missinginfo <- TRUE
} else {
    Wfd     <- reglist$Wfd
}
if (is.null(reglist$yfd)) {
    missinginfo <- TRUE
} else {
    yfd     <- reglist$yfd
}
if (is.null(reglist$y0fd)) {
    missinginfo <- TRUE
} else {
    y0fd    <- reglist$y0fd
}
if (missinginfo) stop("REGLIST does not containing required objects.")

#  generate a fine mesh of argument values

ybasis  <- yfd$basis
yrange  <- ybasis$rangeval
ynbasis <- ybasis$nbasis
nfine   <- 201
argfine <- seq(yrange[1],yrange[2],len=nfine)

#  evaluate the required functional on over this fine mesh

ymat    <- eval.fd(argfine, yfd)
y0mat   <- eval.fd(argfine, y0fd)
yregmat <- eval.fd(argfine, yregfd)
warpmat <- eval.monfd(argfine, Wfd)
warpmat <- yrange[1] + (yrange[2]-yrange[1])*warpmat/
                       (matrix(1,nfine,1)%*%warpmat[nfine,])

#  extract the number of functions NCURVE and the number of variables NVAR

ydim     <- dim(yfd$coef)
ncurve   <- ydim[2]
if (length(ydim) == 2) {
    nvar <- 1
} else {
    nvar <- ydim[3]
}

if (!is.null(names(yfd$fdname)[[3]])) {
    ylabel = names(yfd$fdname)[[3]]
}  else {
    ylabel = "Function value"
}
casename <- yfd$fdname[[2]]
if (length(casename) != ncurve) casename = as.character(1:ncurve)

if (nvar == 1) {

    #  Plot the results if the functions are univariate: NVAR == 1

    if (dim(y0mat)[2] == 1) y0mat = y0mat %*% matrix(1,1,ncurve)
    ylimit   <- c(min(ymat),max(ymat))
    par(mfrow=c(1,2),pty="s",ask=TRUE)
    for (i in 1:ncurve) {
        #  left panel:  plot registered, unregistered and target functions
        plot (argfine, ymat[,i], type="l", lty=2, col=4, ylim=ylimit,
          xlab="Argument value", ylab=ylabel[i], main=paste("Case",casename[i]))
        lines(argfine, y0mat[,i],   lty=2, col=2)
        lines(argfine, yregmat[,i], lty=1, col=4)
        #  right panel: plot warping functions
        plot (argfine, warpmat[,i], type="l",
              xlab="Unwarped argument value", ylab="Warped argument value")
        abline(0,1,lty=2)
    }

}  else {

    #  Plot the results if the functions are multivariate: NVAR > 1

    if (dim(y0mat)[2] == 1) {
        temp <- array(0,dim(ymat))
        for (j in 1:ncurve) temp[,j,] = y0mat
        y0mat <- temp
    }
    ylimit   <- matrix(0,nvar,2)
    varname  <- vector("list",nvar)
    for (j in 1:nvar) {
        ylimit[j,]    <- c(min(ymat[,j,]),max(ymat[,j,]))
        if (!is.null(yfd$fdname[[3]]))
        varname[[j]]  <- yfd$fdname[[3]][[j]]
    }
    for (i in 1:ncurve) {
        par(mfrow=c(nvar,1),pty="m",ask=TRUE)
        for (j in 1:nvar) {
            if (!is.null(varname)) {
                ylabel = varname[[j]]
            }  else {
                ylabel = "Function value"
            }
            #  left panel:  plot registered, unregistered and target functions
            plot (argfine, ymat[,i,j], type="l", lty=2, col=4, ylim=ylimit[j,],
                  xlab="Argument value", ylab=ylabel,
                  main=paste("Case",casename[i],"Variable",varname[[j]]))
            lines(argfine, y0mat[,i,j],   lty=2, col=2)
            lines(argfine, yregmat[,i,j], lty=1, col=4)
        }
        par(mfrow=c(1,1),pty="s",ask=TRUE)
        plot (argfine, warpmat[,i], type="l",
              xlab="Unwarped argument value", ylab="Warped argument value",
              main=paste("Warping function for",casename[i]))
        abline(0,1,lty=2)
    }

}

}