## \subsection{Jacobian and transformation functions for S4}
## \label{sec:jaco}

## \begin{verbatim}

".jacobian" <- function(formulae, variables = NULL) {
  if(is.null(variables))
    variables = all.vars(formulae)
  switch(mode(variables),
         numeric = {
           values = as.list(variables)
           variables = names(values)
         },
         list = {
           values = variables
           variables = names(variables)
         },
         character = {
           values <- NULL
         },
         stop("invalid `variables'"))
  nf = length(formulae)
  nv = length(variables)
  fnames = names(formulae)
  if(is.null(fnames))
    fnames = unlist(lapply(formulae, deparse))
  if(any(i <- fnames == ""))
    fnames[i] = sapply(formulae[i], deparse)
  if(is.null(values)) {
    J = list()
    for(i in 1:nf) {
      tmp = list()
      for(j in 1:nv)
        tmp[[variables[j]]] = D(formulae[[i]],
                                variables[j])
      J[[fnames[i]]] = tmp
    }
    oldClass(J) = "jacobian"
    return(J)
  } else {
    M = matrix(0, nf, nv,
      dimnames = list(fnames, variables))
    for(i in 1:nf)
      for(j in 1:nv)
        M[i, j] = eval(D(formulae[[i]], variables[j]),
                       local = values)
    M
  }
}

jacobian <- function(...) {
  m = list(...)
  l = unlist(lapply(m, is.list))
  f = unlist(lapply(m, inherits, "formula"))
  formulae = c(unlist(m[l], recursive = F), m[f])
  formulae = as.expression(lapply(formulae,
      function(x) x[[length(x)]]))
  variables = unlist(m[!(l | f)])
  .jacobian(formulae, variables)
}

print.jacobian = function(x, ...) {
  rjustify = function(str) {
    m = max(n <- nchar(str))
    blanks = paste(rep(" ", m),
      collapse = "")
    paste(substring(blanks, 0, m - n), str,
          sep = "")
  }
  j = lapply(x, function(y)
             lapply(y, deparse))
  M = matrix(rjustify(as.vector(unlist(j))),
    length(j), length(j[[1]]),
    dimnames = list(paste(names(j[[1]]), ":"),
      paste(names(j), ":")))
  print(t(M), quote = F)
  invisible(x)
}

setClass("transf",
         representation(beta = "named",
                        SE = "named",
                        Vcov = "matrix"))

print.transf = function(x, ...) {
  print(data.frame(beta = x@beta, SE = x@SE))
  invisible(x)
}

setMethod("show", "transf",
          function(object)
          print.transf(object))

setGeneric("transf", function(theta, Vcov, ...)
           standardGeneric("transf"))

setDefaultMethod("transf",
                 function(theta, Vcov, ...) {
  m = as.expression(substitute(f(...))[-1])
  J = .jacobian(m, theta)
  theta = as.list(theta)
  beta = unlist(lapply(m, eval, local = theta))
  Vb = J %*% Vcov %*% t(J)
  SE = sqrt(diag(Vb))
  names(SE) = names(beta)
  new("transf", beta = beta, SE = SE, Vcov = Vb)
})

setMethod("transf", signature(theta = "transf"),
         function(theta, Vcov, ...) {
           m = match.call()
           V = theta@Vcov
           beta = theta@beta
           m$theta = Quote(beta)
           m$Vcov = Quote(V)
           eval(m)
         })
## \end{verbatim}


## \section{A preliminary version of an S4 polynomial class}
## \label{sec:poly}

## \begin{verbatim}
setClass("polynomial", representation(.Data = "numeric"),
         validity = function(object) {
           if(!is(object@.Data, "numeric"))
             return("bad coefficient vector")
           if(any(is.na(object@.Data) | is.inf(object@.Data)))
             return("NAs in coefficient vector")
           ld = length(object@.Data)
           if(ld == 1) return(TRUE)
           if(object@.Data[ld] == 0)
             return("zero leading coefficient")
           TRUE
         }
)

clip = function(p) {
  p = as(p, "numeric")
  while((l = length(p)) > 1 &&
        p[l] == 0) p = p[-l]
  p
}

polynomial = function(p = 0:1) 
  new("polynomial", .Data = clip(p))

setMethod("Math", "polynomial",
          function(x) {
            switch(.Generic,
                   ceiling =,
                   floor =,
                   trunc = polynomial(callGeneric(x@.Data)),
                   stop(paste(.Generic, 
                              "not allowed on polynomials"))
                   )}
          )

setMethod("Math2", "polynomial",
          function(x, digits) {
            x = x@.Data
            polynomial(callGeneric())
          })

.dividePolynomial = function(p, q) {
  lq = length(q)
  if(lq == 1)
    return(p/q)
  p = rev(p)
  q = rev(q)
  r = numeric(0)
  while(length(p) >= length(q)) {
    d = p[1]/q[1]
    r = c(d, r)
    p[1:lq] = p[1:lq] - d * q
    p = p[-1]
  }
  if(length(r) == 0)
    r = 0
  r
}

.multiplyPolynomial = function(p, q) {
  if(length(q) == 1 || length(p) == 1)
    return(p * q)
  m = outer(p, q)
  m = unlist(lapply(split(m, row(m) + col(m)), sum))
  as.vector(m)
}

.remainderPolynomial = function(p, q) {
  lq = length(q)
  if(lq == 1)
    return(0)
  p = rev(p)
  q = rev(q)
  r = numeric(0)
  while(length(p) >= length(q)) {
    d = p[1]/q[1]
    r = c(d, r)
    p[1:lq] = p[1:lq] - d * q
    p = p[-1]
  }
  rev(p)
}

doArithPolynomial = function(e1, e2) {
  e1 = unlist(unclass(e1))
  e2 = unlist(unclass(e2))
  e1.op.e2 = switch(.Generic,
    "+" = ,
    "-" = {
      l = length(e1) - length(e2)
      e1 = c(e1, rep(0, max(0,  - l)))
      e2 = c(e2, rep(0, max(0, l)))
      callGeneric(e1, e2)
    },
    "*" = .multiplyPolynomial(e1, e2),
    "/" = ,
    "%/%" = .dividePolynomial(e1, e2),
    "%%" = .remainderPolynomial(e1, e2),
    "^" = if(length(e2) > 1 || e2 < 0 ||
      (e2 %% 1) != 0) {
      stop("positive integer powers only")
    } else {
      switch(as.character(e2),
             "0" = 1, "1" = e1,
             {
               m = e1
               for(i in 2:e2)
                 m = .multiplyPolynomial(m, e1)
               m
             })
    })
  polynomial(e1.op.e2)
}

setMethod("Arith", signature(e1 = "polynomial", e2 = "missing"),
          function(e1, e2) {
            e1 = e1@.Data
            polynomial(callGeneric())
          })

setMethod("Arith", signature(e1 = "polynomial", 
                             e2 = "polynomial"),
          doArithPolynomial)

setMethod("Arith", signature(e1 = "numeric", e2 =  "polynomial"),
          doArithPolynomial)

setMethod("Arith", signature(e1 = "integer", e2 =  "polynomial"),
          doArithPolynomial)

setMethod("Arith", signature(e1 = "polynomial", e2 =  "numeric"),
          doArithPolynomial)

setMethod("Arith", signature(e1 = "polynomial", e2 =  "integer"),
          doArithPolynomial)

makeBrace <- function(p) {
  p <- rev(p)
  m <- length(p)
  ex <- vector("expression", m + 2)
  ex[[1]] <- structure(expression(p, 0), 
                       mode = "<-")
  for(i in 1:length(p)) {
    ex[[i+1]] <- substitute(p <- v + x*p, 
                            list(v = p[i]))
  }
  ex[[m + 2]] <- as("p", "name")
  mode(ex) <- "{" #}
  ex
}

Polynomial2Function <- function(object) 
          substitute(function(x) body, 
                     list(body = makeBrace(object@.Data)))

setIs("polynomial", "function", coerce = Polynomial2Function)

"poly.calc.coef" = function(x, y) {
  if(missing(y)) {
    ## Polynomial with specifice zeros
    p = 1
    for(xi in x)
      p = c(0, p) - c(xi * p, 0)
    return(p)
  }
  ## Lagrange interpolation formula
  if(any(toss = duplicated(x))) {
    y = y[!toss]
    x = x[!toss]
  }
  if((m = length(x)) != length(y))
    stop("length(x) != length(y)")
  if(m <= 1)
    return(y)
  r = 0
  for(i in 1:m)
    r = r + (y[i] * Recall(x[ - i]))/prod(x[i] - x[ - i])
  r
}

ZeroPoly = function(x) polynomial(poly.calc.coef(x))

Lagrange = function(x, y) polynomial(poly.calc.coef(x, y))          

make.poly.char = function(p) {
  lp = length(p) - 1
  names(p) = 0:lp
  tap = as(abs(p), "character")
  if(any(toss <- match(tap, c("0", "-0"),
                       nomatch = F))) {
    p = p[!toss]
    if(length(p) == 0)
      return("0")
    tap = tap[!toss]
  }
  sp = c("- ", "+ ")[(p > 0) + 1]
  np = names(p)
  p = tap
  p[p == "1" & np != "0"] = ""
  pow = paste("x", "^", np, sep = "")
  pow[np == "0"] = ""
  pow[np == "1"] = "x"
  stars = rep("*", length(p))
  stars[p == "" | pow == ""] = ""
  p = paste(sp, p, stars, pow, sep = "", collapse = " ")
  while((ch = substring(p, 1, 1)) == "+" || ch == " ") 
    p = substring(p, 2)
  p
}

setIs("polynomial", "character", coerce =
      function(object) make.poly.char(object@.Data))

print.polynomial = function(x, ...) {
  p = make.poly.char(zapsmall(x@.Data))
  ow = max(35, options("width")$width - 3)
  m2 = 0
  np = nchar(p)
  while(m2 < np) {
    m1 = m2 + 1
    m2 = m2 + ow
    if(m2 < np)
      while(m2 > m1 && substring(p, m2, m2) != " ") m2 = m2 - 1
    cat("# ", substring(p, m1, m2), "\n")
  }
  invisible(x)
}

setMethod("show", "polynomial",
          function(object) print.polynomial(object))          

polyZeros <- function(z) {
        r = polyroot(z@.Data)
        zapsmall(Re(r)) + zapsmall(Im(r)) * (1i)
}

solve.polynomial <- function(a, b, ...) 
  if(missing(b)) polyZeros(a) else polyZeros(a-b)

## \end{verbatim}
## \end{document}