basisfd.R

#  Generator function of class basisfd

basisfd <- function(type, rangeval, nbasis, params, dropind=vector("list",0),
                    quadvals=vector("list",0), values=vector("list",0),
                    basisvalues=vector("list",0))
{
  #  BASISFD  generator function of "basisfd" class.
  #  Arguments:
  #  TYPE    ...a string indicating the type of basisobj.
  #             This may be one of:
  #             "Bspline", "bspline", "Bsp", "bsp",
  #             "con", "const", "constant"
  #             "exp", "exponen", "exponential"
  #             "fdVariance"
  #             "FEM"
  #             "Fourier", "fourier", "Fou", "fou",
  #             "mon", "monom", "monomial",
  #             "polyg", "polygon", "polygonal"
  #             "power" "pow"
  #  RANGEVAL...an array of length 2 containing the lower and upper
  #             boundaries for (the rangeval of argument values
  #             If basis is of FEM type, rangeval is not used
  #  NBASIS ... the number of basis functions
  #  PARAMS ... If the basis is "fourier", this is a single number indicating
  #               the period.  That is, the basis functions are periodic on
  #               the interval (0,PARAMS) or any translation of it.
  #             If the basis is "bspline", the values are interior points at
  #               which the piecewise polynomials join.
  #               Note that the number of basis functions NBASIS is equal
  #               to the order of the Bspline functions plus the number of
  #               interior knots, that is the length of PARAMS.
  #             This means that NBASIS must be at least 1 larger than the
  #               length of PARAMS.
  #  DROPIND...A set of indices in 1:NBASIS of basis functions to drop when
  #              basis objects are arguments.  Default is vector("list",0)
  #              Note that argument NBASIS is reduced by the number of
  #              indices, and the derivative matrices in VALUES are also clipped.
  #  QUADVALS...A NQUAD by 2 matrix.  The firs t column contains quadrature
  #              points to be used in a fixed point quadrature.  The second
  #              contains quadrature weights.  For example, for (Simpson"s
  #              rule for (NQUAD = 7, the points are equally spaced and the
  #              weights are delta.*[1, 4, 2, 4, 2, 4, 1]/3.  DELTA is the
  #              spacing between quadrature points.  The default is
  #              matrix("numeric",0,0).
  #  VALUES ...A list, with entries containing the values of
  #              the basis function derivatives starting with 0 and
  #              going up to the highest derivative needed.  The values
  #              correspond to quadrature points in QUADVALS and it is
  #              up to the user to decide whether or not to multiply
  #              the derivative values by the square roots of the
  #              quadrature weights so as to make numerical integration
  #              a simple matrix multiplication.
  #              Values are checked against QUADVALS to ensure the correct
  #              number of rows, and against NBASIS to ensure the correct
  #              number of columns.
  #              The default value of is VALUES is vector("list",0).
  #              VALUES contains values of basis functions and derivatives at
  #              quadrature points weighted by square root of quadrature weights.
  #              These values are only generated as required, and only if slot
  #              QUADVALS is not matrix("numeric",0,0).
  #  BASISVALUES...A vector of lists, allocated by code such as
  #              vector("list",1).
  #              This field is designed to avoid evaluation of a
  #              basis system repeatedly at a set of argument values.
  #              Each list within the vector corresponds to a specific set
  #              of argument values, and must have at least two components,
  #              which may be tagged as you wish.
  #              The first component in an element of the list vector contains the
  #              argument values.
  #              The second component in an element of the list vector
  #              contains a matrix of values of the basis functions evaluated
  #              at the arguments in the first component.
  #              The third and subsequent components, if present, contain
  #              matrices of values their derivatives up to a maximum
  #              derivative order.
  #              Whenever function getbasismatrix is called, it checks
  #              the first list in each row to see, first, if the number of
  #              argument values corresponds to the size of the first dimension,
  #              and if this test succeeds, checks that all of the argument
  #              values match.  This takes time, of course, but is much
  #              faster than re-evaluation of the basis system.  Even this
  #              time can be avoided by direct retrieval of the desired
  #              array.
  #              For example, you might set up a vector of argument values
  #              called "evalargs" along with a matrix of basis function
  #              values for these argument values called "basismat".
  #              You might want too use tags like "args" and "values",
  #              respectively for these.  You would then assign them
  #              to BASISVALUES with code such as
  #                basisobj$basisvalues <- vector("list",1)
  #                basisobj$basisvalues[[1]] <-
  #                             list(args=evalargs, values=basismat)
  #
  #  Returns
  #  BASISOBJ  ... a basisfd object with slots
  #         type
  #         rangeval
  #         nbasis
  #         params
  #         dropind
  #         quadvals
  #         values
  #         basisvalues
  #  Slot VALUES contains values of basis functions and derivatives at
  #   quadrature points weighted by square root of quadrature weights.
  #   These values are only generated as required, and only if slot
  #   quadvals is not empty.
  #
  #  An alternative name for (this function is CREATE.BASIS, but PARAMS argument
  #     must be supplied.
  #  Specific types of bases may be set up more conveniently using functions
  #  CREATE.BSPLINE.BASIS     ...  creates a b-spline basis
  #  CREATE.CONSTANT.BASIS    ...  creates a constant basis
  #  CREATE.EXPONENTIAL.BASIS ...  creates an exponential basis
  #  CREATE.FDVARIANCE.BASIS  ...  creates an fdVariance basis
  #  CREATE.FEM.BASIS         ...  creates an FEM basis  
  #  CREATE.FOURIER.BASIS     ...  creates a fourier basis
  #  CREATE.MONOMIAL.BASIS    ...  creates a monomial basis
  #  CREATE.POLYGON.BASIS     ...  creates a polygonal basis
  #  CREATE.POLYNOMIAL.BASIS  ...  creates a polynomial basis
  #  CREATE.POWER.BASIS       ...  creates a monomial basis
  
  #  Last modified 3 January 2020 by Jim Ramsay
  # value -> values 2012.12.27 by spencer graves
  
  #  Set up default basis if there are no arguments:
  #     order 2 monomial basis over [0,1]
  
  if (nargs()==0) {
    type        <- "bspline"
    rangeval    <- c(0,1)
    nbasis      <- 2
    params      <- vector("list",0)
    dropind     <- vector("list",0)
    quadvals    <- vector("list",0)
    values      <- vector("list",0)
    basisvalues <- vector("list",0)
    
    basisobj  <- list(type=type,     rangeval=rangeval, nbasis=nbasis,
                      params=params, dropind=dropind,   quadvals=quadvals,
                      values=values, basisvalues=basisvalues)
    oldClass(basisobj) <- "basisfd"
    return(basisobj)
  }
  
  #  if first argument is a basis object, return
  
  if (inherits(type,"basisfd")) {
    basisobj <- type
    return(basisobj)
  }
  
  #  check basistype
  
  # type <- moreNames(type)
  
  #  recognize type of basis by use of several variant spellings
  
  if(type == "bspline" ||
     type == "Bspline" ||
     type == "spline"  ||
     type == "Bsp"     ||
     type == "bsp") {
    type = "bspline"
  }
  else if(type == "con"      ||
          type == "const"    ||
          type == "constant") {
    type = "const"
  }
  else if(type == "exp"    ||
          type == "expon"  ||
          type == "exponential") {
    type = "expon"
  }
  else if(type == "fdVariance" ||
          type == "fdVar") {
    type = "fdVariance"
  }
  else if(type == "FEM") {
    type = "FEM"
  }
  else if(type == "Fourier" ||
          type == "fourier" ||
          type == "Fou"     ||
          type == "fou") {
    type = "fourier"
  }
  else if(type == "mon" ||
          type == "monom"  ||
          type == "monomial") {
    type = "monom"
  }
  else if(type == "polyg"    ||
          type == "polygon"  ||
          type == "polygonal") {
    type = "polygonal"
  }
  else if(type == "polynomial"    ||
          type == "polynom") {
    type = "polynomial"
  }
  else if(type == "pow"    ||
          type == "power") {
    type = "power"
  }
  else {
    type = "unknown"
  }
  
  if (type=="unknown"){
    stop("'type' unrecognizable.")
  }
  
  #  check rangeval if the object is not of type FEM
  
  if (!type == "FEM") {
    rangeval = as.vector(rangeval)
    if (!is.numeric(rangeval)) stop("Argument rangeval is not numeric.")
    if (length(rangeval) != 2) stop("Argument rangeval is not of length 2.")
    if (!(rangeval[2] > rangeval[1]))
      stop("Argument rangeval is not strictly increasing.")
  }
  
  #  check nbasis
  
  if (nbasis <= 0)             stop("Argument nbasis is not positive.")
  if (round(nbasis) != nbasis) stop("Argument nbasis is not an integer.")
  
  #  check if QUADVALS is present, and set to default if not
  
  if (missing(quadvals)) quadvals <- vector("list",0)
  else if(!(length(quadvals) == 0 || is.null(quadvals))){
    nquad <- dim(quadvals)[1]
    ncol  <- dim(quadvals)[2]
    if ((nquad == 2) && (ncol > 2)){
      quadvals <- t(quadvals)
      nquad    <- dim(quadvals)[1]
      ncol     <-dim(quadvals)[2]
    }
    if (nquad < 2) stop("Less than two quadrature points are supplied.")
    if (ncol != 2) stop("'quadvals' does not have two columns.")
  }
  
  #  check VALUES is present, and set to a single empty list if not.
  if(!(length(values) == 0 || missing(values) || is.null(values))) {
    n <- dim(values)[1]
    k <- dim(values)[2]
    if (n != nquad)
      stop(paste("Number of rows in 'values' not equal to number of",
                 "quadrature points."))
    if (k != nbasis)
      stop(paste("Number of columns in 'values' not equal to number of",
                 "basis functions."))
  }
  else values <- vector("list",0)
  
  #  check BASISVALUES is present, and set to vector("list",0) if not.
  #  If present, it must be a two-dimensional list created by a command like
  #  listobj <- matrix("list", 2, 3)
  
  if(!(length(basisvalues) == 0 || missing(basisvalues) || !is.null(basisvalues))) {
    if (!is.list(basisvalues)) stop("BASISVALUES is not a list object.")
    sizevec <- dim(basisvalues)
    if (length(sizevec) != 2) stop("BASISVALUES is not 2-dimensional.")
    for (i in 1:sizevec[1]) {
      if (length(basisvalues[[i,1]]) != dim(basisvalues[[i,2]])[1]) stop(
        paste("Number of argument values not equal number",
              "of values."))
    }
  }
  else basisvalues <- vector("list",0)
  
  #  check if DROPIND is present, and set to default if not
  
  if(missing(dropind)) dropind <- vector("list",0)
  
  if (length(dropind) > 0) {
    #  check DROPIND
    ndrop = length(dropind)
    if (ndrop >= nbasis) stop('Too many index values in DROPIND.')
    dropind = sort(dropind)
    if (ndrop > 1 && any(diff(dropind)) == 0)
      stop('Multiple index values in DROPIND.')
    for (i in 1:ndrop) {
      if (dropind[i] < 1 || dropind[i] > nbasis)
        stop('A DROPIND index value is out of range.')
    }
    #  drop columns from VALUES cells if present
    nvalues = length(values)
    if (nvalues > 0 && length(values[[1]] > 0)) {
      for (ivalue in 1:nvalues) {
        derivvals = values[[ivalue]]
        derivvals = derivvals[,-dropind]
        values[[ivalue]] = derivvals
      }
    }
  }
  
  #  select the appropriate type and process
  
  if (type=="fourier"){
    paramvec   <- rangeval[2] - rangeval[1]
    period     <- params[1]
    if (period <= 0)  stop("Period must be positive for (a Fourier basis")
    params <- period
    if ((2*floor(nbasis/2)) == nbasis)  nbasis <- nbasis + 1
  } else if(type=="bspline"){
    if (!missing(params)){
      nparams  <- length(params)
      if(nparams>0){
        if (params[1] <= rangeval[1])
          stop("Smallest value in BREAKS not within RANGEVAL")
        if (params[nparams] >= rangeval[2])
          stop("Largest value in BREAKS not within RANGEVAL")
      }
    }
  } else if(type=="expon") {
    if (length(params) != nbasis)
      stop("No. of parameters not equal to no. of basis fns for (exponential basisobj$")
  }   else if(type=="fdVariance") {
    if (length(params) != 2)
      stop("No. of parameters not equal to 8 for (FEM basisobj$")
  } else if(type=="FEM") {
    # if (length(params) != 8)
    #   stop("No. of parameters not equal to 8 for (FEM basisobj$")
  } else if(type=="polynomial") {
    if (length(params) != nbasis)
      stop("No. of parameters not equal to no. of basis fns for (polynomialal basisobj$")
  } else if(type=="power") {
    if (length(params) != nbasis)
      stop("No. of parameters not equal to no. of basis fns for (power basisobj$")
  } else if(type=="const") {
    params <- 0
  } else if(type=="monom") {
    if (length(params) != nbasis)
      stop("No. of parameters not equal to no. of basis fns for (monomial basisobj$")
  } else if(type=="polygonal") {
    if (length(params) != nbasis)
      stop("No. of parameters not equal to no. of basis fns for polygonal basisobj$")
  } else stop("Unrecognizable basis")
  
  #  Save call
  
  obj.call <- match.call()
  
  #  S4 definition
  
  # basisobj <- new("basisfd", call=obj.call, type=type, rangeval=rangeval,
  #                 nbasis=nbasis,  params=params, dropind=dropind,
  #                 quadvals=quadvals, values=values, basisvalues=basisvalues)
  
  #  S3 definition
  
  basisobj <- list(call=obj.call, type=type, rangeval=rangeval, nbasis=nbasis,
                   params=params, dropind=dropind, quadvals=quadvals,
                   values=values, basisvalues=basisvalues)
  oldClass(basisobj) <- "basisfd"
  
  basisobj
  
}

#  --------------------------------------------------------------------------
#                  print for basisfd class
#  --------------------------------------------------------------------------

print.basisfd <- function(x, ...)
{
  
  #  Last modified 3 January 2008 by Jim Ramsay
  
  basisobj <- x
  cat("\nBasis object:\n")
  if (!inherits(basisobj, "basisfd"))
    stop("Argument not a functional data object")
  
  #  print type
  
  cat(paste("\n  Type:  ", basisobj$type,"\n"))
  
  #  print range
  
  cat(paste("\n  Range: ", basisobj$rangeval[1],
            " to ",        basisobj$rangeval[2],"\n"))
  
  #  return if a constant basis
  
  if (basisobj$type == "const") {
    return()
  }
  
  #  print number of basis functions
  
  cat(paste("\n  Number of basis functions: ",
            basisobj$nbasis,     "\n"))
  
  #  print parameters according to type of basis
  
  if (basisobj$type == "fourier")
    cat(paste("\n  Period: ",basisobj$params,"\n"))
  if (basisobj$type == "bspline") {
    norder <- basisobj$nbasis - length(basisobj$params)
    cat(paste("\n  Order of spline: ", norder, "\n"))
    if (length(basisobj$params) > 0) {
      print("  Interior knots")
      print(basisobj$params)
    } else {
      print("  There are no interior knots.")
    }
  }
  if (basisobj$type == "polyg") {
    print("  Argument values")
    print(basisobj$params)
  }
  if (basisobj$type == "expon") {
    print("  Rate coefficients")
    print(basisobj$params)
  }
  if (basisobj$type == "monom") {
    print("  Exponents")
    print(basisobj$params)
  }
  if (basisobj$type == "power") {
    print("  Exponents")
    print(basisobj$params)
  }
  
  
  #  display indices of basis functions to be dropped
  
  if (length(basisobj$dropind) > 0) {
    print("  Indices of basis functions to be dropped")
    print(basisobj$dropind)
  }
  
}

#  --------------------------------------------------------------------------
#                  summary for basisfd class
#  --------------------------------------------------------------------------

summary.basisfd <- function(object, ...)
{
  basisobj <- object
  cat("\nBasis object:\n")
  if (!inherits(basisobj, "basisfd"))
    stop("Argument not a functional data object")
  cat(paste("\n  Type:  ", basisobj$type,"\n"))
  cat(paste("\n  Range: ", basisobj$rangeval[1],
            " to ",        basisobj$rangeval[2],"\n"))
  if (basisobj$type == "const") {
    return()
  }
  cat(paste("\n  Number of basis functions: ",
            basisobj$nbasis,     "\n"))
  if (basisobj$type == "fourier")
    cat(paste("\n  Period: ",basisobj$params,"\n"))
  if (length(basisobj$dropind) > 0) {
    print(paste(length(basisobj$dropind),
                "indices of basis functions to be dropped"))
  }
}

#  --------------------------------------------------------------------------
#  equality for basisfd class
#  --------------------------------------------------------------------------

"==.basisfd" <- function(basis1, basis2)
{
  
  # EQ assesses whether two bases are equivalent.
  
  #  Last modified 1 January 2007
  
  type1   <- basis1$type
  range1  <- basis1$rangeval
  nbasis1 <- basis1$nbasis
  pars1   <- basis1$params
  drop1   <- basis1$dropind
  
  type2   <- basis2$type
  range2  <- basis2$rangeval
  nbasis2 <- basis2$nbasis
  pars2   <- basis2$params
  drop2   <- basis2$dropind
  
  basisequal <- TRUE
  
  #  check types
  
  if (!(type1 == type2)) {
    basisequal <- FALSE
    return(basisequal)
  }
  
  #  check ranges
  
  if (range1[1] != range2[1] || range1[2] != range2[2]) {
    basisequal <- FALSE
    return(basisequal)
  }
  
  #  check numbers of basis functions
  
  if (nbasis1 != nbasis2) {
    basisequal <- FALSE
    return(basisequal)
  }
  
  #  check parameter vectors
  
  if (!(all(pars1 == pars2))) {
    basisequal <- FALSE
    return(basisequal)
  }
  
  #  check indices of basis function to drop
  
  if (!(all(drop1 == drop2))) {
    basisequal <- FALSE
    return(basisequal)
  }
  
  return(basisequal)
  
}

#  --------------------------------------------------------------------------
#  pointwise multiplication method for basisfd class
#  --------------------------------------------------------------------------

"*.basisfd" <- function (basisobj1, basisobj2)
{
  # TIMES for (two basis objects sets up a basis suitable for (
  #  expanding the pointwise product of two functional data
  #  objects with these respective bases.
  # In the absence of a true product basis system in this code,
  #  the rules followed are inevitably a compromise:
  #  (1) if both bases are B-splines, the norder is the sum of the
  #      two orders - 1, and the breaks are the union of the
  #      two knot sequences, each knot multiplicity being the maximum
  #      of the multiplicities of the value in the two break sequences.
  #      Order, however, is not allowed to exceed 20.
  #      That is, no knot in the product knot sequence will have a
  #      multiplicity greater than the multiplicities of this value
  #      in the two knot sequences.
  #      The rationale this rule is that order of differentiability
  #      of the product at each value will be controlled  by
  #      whichever knot sequence has the greater multiplicity.
  #      In the case where one of the splines is order 1, or a step
  #      function, the problem is dealt with by replacing the
  #      original knot values by multiple values at that location
  #      to give a discontinuous derivative.
  #  (2) if both bases are Fourier bases, AND the periods are the
  #      the same, the product is a Fourier basis with number of
  #      basis functions the sum of the two numbers of basis fns.
  #  (3) if only one of the bases is B-spline, the product basis
  #      is B-spline with the same knot sequence and order two
  #      higher.
  #  (4) in all other cases, the product is a B-spline basis with
  #      number of basis functions equal to the sum of the two
  #      numbers of bases and equally spaced knots.
  
  #  Of course the ranges must also match.
  
  #  Last modified 2022.04.19 by Jim Ramsay
  
  #  check the ranges
  
  range1 <- basisobj1$rangeval
  range2 <- basisobj2$rangeval
  if (range1[1] != range2[1] || range1[2] != range2[2])
    stop("Ranges are not equal.")
  
  #  get the types
  
  type1 <- basisobj1$type
  type2 <- basisobj2$type
  
  #  deal with constant bases
  
  if (type1 == "const" && type2 == "const") {
    prodbasisobj <- create.constant.basis(range1)
    return(prodbasisobj)
  }
  
  if (type1 == "const") {
    prodbasisobj <- basisobj2
    return(prodbasisobj)
  }
  
  if (type2 == "const") {
    prodbasisobj <- basisobj1
    return(prodbasisobj)
  }
  
  #  get the numbers of basis functions
  
  nbasis1 <- basisobj1$nbasis
  nbasis2 <- basisobj2$nbasis
  
  #  work through the cases
  
  if (type1 == "bspline" && type2 == "bspline") {
    #  both are bases B-splines
    #  get orders
    interiorknots1 <- basisobj1$params
    interiorknots2 <- basisobj2$params
    #    uniqueknots    <- sort(union(interiorknots1, interiorknots2))
    interiorknots1.2 <- union(interiorknots1, interiorknots2)
    uniqueknots <- {
      if(is.null(interiorknots1.2)) NULL else sort(interiorknots1.2)
    }
    nunique <- length(uniqueknots)
    multunique <- rep(0,nunique)
    for (i in seq(length=nunique)) {
      mult1 <- {
        if(length(interiorknots1)>0)
          length(interiorknots1[interiorknots1==uniqueknots[i]])
        else 0
      }
      mult2 <- {
        if(length(interiorknots2)>0)
          length(interiorknots2[interiorknots2==uniqueknots[i]])
        else 0
      }
      multunique[i] <- max(mult1,mult2)
    }
    #
    allknots <- rep(0,sum(multunique))
    m2 <- 0
    for (i in seq(length=nunique)) {
      m1 <- m2 + 1
      m2 <- m2 + multunique[i]
      allknots[m1:m2] <- uniqueknots[i]
    }
    norder1 <- nbasis1 - length(interiorknots1)
    norder2 <- nbasis2 - length(interiorknots2)
    #  norder is not allowed to exceed 20
    norder  <- min(c(norder1 + norder2 - 1,20))
    allbreaks  <- c(range1[1], allknots, range1[2])
    nbasis <- length(allbreaks) + norder - 2
    prodbasisobj <-
      create.bspline.basis(range1, nbasis, norder, allbreaks)
    return(prodbasisobj)
  }
  
  if (type1 == "fourier" && type2 == "fourier") {
    #  both bases Fourier
    #  check whether periods match
    #  if they do not, default to the basis below.
    period1 <- basisobj1$params
    period2 <- basisobj2$params
    nbasis  <- nbasis1 + nbasis2-1
    if (period1 == period2) {
      prodbasisobj <- create.fourier.basis(range1, nbasis, period1)
      return(prodbasisobj)
    }
  }
  
  #  Default case when all else fails: the product basis is B-spline
  #  When neither basis is a B-spline basis, the order
  #  is the sum of numbers of bases, but no more than 8.
  #  When one of the bases if B-spline and the other isn"t,
  #  the order is the smaller of 8 or the order of the spline
  #  plus 2.  Under no circumstances can the order exceed 20, however.
  #  See BsplineS where this restriction is tested.
  
  if (type1 == "bspline" || type2 == "bspline") {
    norder <- 8
    if (type1 == "bspline") {
      interiorknots1 <- basisobj1$params
      norder1        <- nbasis1 - length(interiorknots1)
      norder         <- min(c(norder1+2, norder))
    }
    if (type2 == "bspline") {
      interiorknots2 <- basisobj2$params
      norder2        <- nbasis2 - length(interiorknots2)
      norder         <- min(c(norder2+2, norder))
    }
  } else {
    #  neither basis is B-spline
    norder <- min(c(8, nbasis1+nbasis2))
  }
  #  set up the default B-spline product basis
  nbasis <- max(c(nbasis1+nbasis2, norder+1))
  prodbasisobj <- create.bspline.basis(range1, nbasis, norder)
  return(prodbasisobj)
}

#  ---------------------------------------------------------
#        Subscripted reference to a basis object
#  ---------------------------------------------------------

#  Last modified 22 December 2007

"[.basisfd" <- function(basisobj, subs=TRUE)
{
  #  select subsets of basis functions in a basis object
  
  dropind = vector("numeric", 0)
  nbasis <- basisobj$nbasis
  for (i in 1:nbasis) {
    if (!any(subs==i)) dropind = c(dropind, i)
  }
  basisobj$dropind <- dropind
  return(basisobj)
}