CSTRfn.R

CSTRfn <- function(parvec, datstruct, fitstruct,
                   CSTRbasis, lambda, gradwrd=TRUE){
#function [res, Dres, fitstruct, df, gcv] =
#    CSTRfn(parvec, datstruct, fitstruct, CSTRbasis, lambda, gradwrd)

# Last modified with R port 2007.07.21 by Spencer Graves
#%  Previously modified 29 July 2005

  max.log.betaCC <- (log(.Machine$double.xmax)/3)
# For certain values of 'coef',
# naive computation of betaCC will return +/-Inf,
# which generates NAs in Dres.
# Avoid this by clipping betaCC
#
# log(.Machine$double.xmax)/2 is too big,
# because a multiple of it is squared in CSTRfn ...
#

#if nargin < 6, gradwrd = 1;  end
##
## 1.  Modify fitstruct for use with CSTRfitLS
##
#  cat("CSTRfn: parvec =", parvec, "\n")
#
  eps <- .Machine$double.eps
  eps1 <- (1+2*eps)
  fit <- fitstruct$fit
  if(is.null(fit))
    stop("fitstruct has no 'fit' component")

#%  load parameters
#[kref, EoverR, a, b] = par2vals(parvec, fitstruct)
#%  set up fitstruct
#fitstruct.kref   = kref
#fitstruct.EoverR = EoverR
#fitstruct.a      = a
#fitstruct.b      = b

  estimate <- fitstruct$estimate
  if(is.null(estimate))
    stop("fitstruct has no 'estimate' component")
#  Compare estimate with parvec
  if(sum(estimate) != length(parvec)){
    cat("ERROR in CSTRfn:  sum(estimate) != length(parvec)\n")
    cat("parvec = ", parvec, "; \n")
    stop("fitstruct$estimate = ",
         paste(estimate, collapse=" "))
  }
  m <- 0
  {
#  1.1.  Estimate kref starting from parvec or use fitstruct?
    if(estimate[1]){
      m <- m+1
      kref <- parvec[m]
      fitstruct$kref <- kref
    }
    else
      kref <- fitstruct$kref
  }
  {
#  1.2.  Estimate EoverR starting from parvec or use fitstruct?
    if(estimate[2]){
      m <- m+1
      EoverR <- parvec[m]
      fitstruct$EoverR <- EoverR
    }
    else
      EoverR<- fitstruct$EoverR
  }
  {
#  1.3.  Estimate 'a' starting from parvec or use fitstruct?
    if(estimate[3]){
      m <- m+1
      a <- parvec[m]
      fitstruct$a <- a
    }
    else
      a <- fitstruct$a
  }
  {
#  1.4.  Estimate b starting from parvec or use fitstruct?
    if(estimate[4]){
      m <- m+1
      b <- parvec[m]
      fitstruct$b <- b
    }
    else
      b<- fitstruct$b
  }
##
## 2.  Set up inner optimization:  optimize fit with respect to coef
##
#  2.1.  Set up
  tolval <- 1e-10
  itermax <- 10
  coef0 <- as.matrix(fitstruct$coef0)
  dim.coef0 <- dim(coef0)
  if(length(dim.coef0)>1){
    coef0 <- as.vector(coef0)
    if(length(dim.coef0)==2 && (dim.coef0[2]==2)){
      coefNames <- outer(c("Conc", "Temp"), 1:dim.coef0[1], paste,
                         sep="")
      names(coef0) <- t(coefNames)
    }
  }
  ncoef <- length(coef0)
#  2.2.  initial value
#[res0, Dres] = CSTRfitLS(coef0, datstruct, fitstruct, lambda, 1)
  CSTR0 <- CSTRfitLS(coef0, datstruct, fitstruct, lambda, 1)

#  2.3.  Reshape for easier manipulation
  res0 <- with(CSTR0$res, c(Sres, Lres))
#
  N <- dim(CSTR0$res$Sres)[1]
  k12 <- dim(CSTR0$res$Sres)[2]
  nquad <- dim(CSTR0$res$Lres)[1]
  n.Sres <- N*k12
  n.Lres <- nquad*2
  Nval <- n.Sres + n.Lres
#
  nbasis <- (dim(CSTR0$Dres$DSres)[2]/2)
  Dres0 <- with(CSTR0$Dres, rbind(DSres, DLres))
# res0.mat <- readMat("CSTRfitLSres0.mat")# OK
# d.res0 <- (res0-res0.mat$res0)
# quantile(d.res0)
#      0%      25%      50%      75%     100%
#-7.1e-06 -1.4e-07  1.3e-08  2.1e-07  6.0e-06
# res0 good.

# Dres.mat <- read.csv("CSTRfitLSDres.csv", header=FALSE)
# d.Dres <- (Dres0-as.matrix(Dres.mat))
# quantile(d.Dres)
#           0%           25%           50%           75%          100%
#-0.0006900152  0.0000000000  0.0000000000  0.0000000000  0.0006814540
# sqrt(mean(Dres0^2)) # 0.747
# sqrt(mean(d.Dres^2)) # 9.8e-6
# Good.

# F0 = mean(res0.^2)
  F0 <- mean(res0^2)
  F1 <- 0
  iter <- 0
  fundif <- 1
# gradnorm0 = mean(res0'*Dres)
  r.D <- crossprod(res0, Dres0)
  gradnorm0 <- mean(r.D)
  if(is.na(gradnorm0))
    stop("Initial call to CSTRfitLS returned NAs with parvec = ",
         paste(parvec, collapse=", "), ";  sum(is.na(res0)) = ",
         sum(is.na(res0)), "; sum(is.na(Dres0)) = ", sum(is.na(Dres0)))
#
  gradnorm1 <- 1
##
## 3.  Gauss-Newton optimization loop:
##     optimize fit with respect to coef
##
  while(( gradnorm1 > tolval) | (fundif > tolval)){
#
    iter <- iter + 1
    if( iter > itermax) break
#
#    Dcoef = Dres\res0
    nNA.res0 <- sum(is.na(res0))
    nNA.Dres0 <- sum(is.na(Dres0))
    if(nNA.res0 | nNA.Dres0) {
      dump("parvec", "parvecError.R")
      cat("Error: ", nNA.res0, "and", nNA.Dres0,
          "NAs found in res0 and Dres0 with parvec =",
          parvec, "; parvec dumped to parvecError.R\n")
    }
#
#    Dcoef <- (lm.fit(Dres0, res0)$coefficients)
#
    Dres.svd <- svd(Dres0)
    ikeep <- with(Dres.svd, which(d > eps*max(d)))
    Dres.rank <- length(ikeep)
    if(Dres.rank < min(dim(Dres0)))
      warning("Dres0 has rank ", Dres.rank, " < dim(Dres0) = ",
              paste(dim(Dres0), collapse=", "),
              " in iteration ", iter,
              ".  Ignoring singular subspace.  ",
              "First (rank+1) singular values = ",
              paste(Dres.svd$d[1:(Dres.rank+1)], collapse=", "))
    Dcoef <- with(Dres.svd, v %*% ((t(u) %*% res0) / d))
#
#    %  initial step:  alpha = 1
    coef1 <- coef0 - Dcoef
#
#[res1, Dres] = CSTRfitLS(coef1, datstruct, fitstruct, lambda, 1)
    CSTR1 <- CSTRfitLS(coef1, datstruct, fitstruct, lambda, 1)
    res1 <- with(CSTR1$res, c(Sres, Lres))
    Dres1 <- with(CSTR1$Dres, rbind(DSres, DLres))
#
    F1 <- mean(res1^2)
    alpha <- 1
    fundif <- abs(F0-F1)/abs(F0)

#%     %  smaller steps as required, halving alpha each time
#%     while F1 >= F0*(1+2*eps)
    while(is.na(F1) || F1>= (eps1*F0) || any(is.na(Dres1))){
      alpha <- alpha/2
      if(is.na(F1)){
        n.na <- sum(is.na(res1))
        attr(coef1, "n.na.in.res1") <- n.na
#        ..CSTRfn.coef1.gen.NA <<- list(parvec=parvec, coef1=coef1)
        warning(n.na, " NAs in res1.")
#        warning(n.na, " NAs in res1;  coef1 saved in ",
#                "'..CSTRfn.coef1.gen.NA'")
      }
      if(alpha< 1e-4){
        pv <- paste(signif(parvec, 3), collapse=", ")
#        ..CSTRfn.coef1.gen.5 <<- list(parvec=parvec, coef1=coef1)
#        warning("Stepsize reduced below the minimum with parvec = ",
#             pv, " on iteration ", iter,
#             " in trying to optimize ", length(coef0),
#             " coefficients;  using suboptimal coefficients;  ",
#                "saved in '..CSTRfn.coef1.gen.5'")
        warning("Stepsize reduced below the minimum with parvec = ",
             pv, " on iteration ", iter,
             " in trying to optimize ", length(coef0),
             " coefficients;  using suboptimal coefficients. ")
        break
      }
#
      coef1 <- coef0 - alpha*Dcoef
#%    [res1, Dres] = CSTRfitLS(coef1, datstruct, fitstruct, lambda, 1)
      CSTR1 <- CSTRfitLS(coef1, datstruct, fitstruct, lambda, 1)
      res1 <- with(CSTR1$res, c(Sres, Lres))
      Dres1 <- with(CSTR1$Dres, rbind(DSres, DLres))
      F1 <- mean(res1^2)
      fundif <- abs(F0-F1)/abs(F0)
    }
#    gradnorm1 = mean(res1'*Dres)
    gradnorm1 <- mean(crossprod(res1, Dres1))
#%     disp(num2str([iter, F1, fundif, gradnorm1]))
    coef0 <- coef1
    res0 <- res1
    Dres0 <- Dres1
    F0 <- F1
    gradnorm0 <- gradnorm1
# end while(( gradnorm1 > tolval) | (fundif > tolval)){
  }
##
##% 4.  update coef
##
  coef <- coef0
  fitstruct$coef0 <- coef
#
# DresMat <- read.csv("CSTRfnDres.csv", header=FALSE)
# d.Dres <- Dres - as.matrix(DresMat)
# quantile(d.Dres)
#           0%           25%           50%           75%          100%
#-0.0006207360  0.0000000000  0.0000000000  0.0000000000  0.0006317204
# sqrt(mean(Dres^2)) # 0.747
# sqrt(mean(d.Dres^2))# 1.0e-5
#
##
##% 5.  compute df and gcv
##
  Zmat <- Dres0[1:ncoef,]
#  Nval = length(res1)
  Rfac <- Dres0[(ncoef+1):Nval,]
#  Smat = Zmat*inv(Zmat'*Zmat + Rfac'*Rfac)*Zmat'
# Use singular value decomposition so we never have to worry about
# ill conditioning.
  Zsvd <- svd(Zmat)
# Zmat = with(Zsvd, u %*% diag(d) %*% t(v))
# so Z'Z+R'R = v d^2 v' + R'R
#          = v %*% (d^2 + (R%*%v)'(R%*%v))
  Rv <- (Rfac %*% Zsvd$v)
  d2.vR.Rv <- (diag(Zsvd$d^2)+crossprod(Rv))
  d2.eig <- eigen(d2.vR.Rv, symmetric=TRUE)
# Check for ill conditioning
  d2.ev <- d2.eig$values
  ZR.rank <- sum(d2.ev>0)
  {
    if(ZR.rank<1){
      warning("Z'Z+R'R = 0")
      Smat <- array(0, dim=c(ncoef, ncoef))
    }
    else{
      if(ZR.rank<ncoef)
        warning("Z'Z+R'R is not positive definite.  rank = ",
                ZR.rank, ";  dim = ", ncoef, "; increasing the ",
                ncoef-ZR.rank, " smallest eigenvalues ",
                "to improve numeric stability.")
      d2.evMin <- eps*d2.ev[1]
      ZR.rank1 <- sum(d2.ev >= d2.evMin)
      if(ZR.rank1 < ZR.rank)
        warning("Z'Z+R'R is ill conditioned.  Increasing the ",
                ncoef-ZR.rank1, " smallest eigenvalues ",
                "to improve numeric stability.")
#  Z'(inv(Z'Z+R'R)Z
#     = (u%*%d%*%w) solve(LAM) (udw)'
      udw <- ((Zsvd$u * rep(Zsvd$d, each=ncoef)) %*% d2.eig$vectors)
# or
#udw1<-((Zsvd$u[,1:ZR.rank1, drop=FALSE]*rep(Zsvd$d[1:ZR.rank1],e=ncoef))%*%d2.eig$vectors[1:ZR.rank1,,drop=FALSE])
      d2.ev2 <- pmax(d2.ev, d2.evMin)
      Smat <- ((udw / rep(d2.ev2, each=ncoef)) %*% t(udw))
# or
#Smat1<-((udw1/rep(d2.ev2, each=ncoef)) %*% t(udw1))
    }
  }
  df. <- sum(diag(Smat))
  dfe <- ncoef - df.
  gcv <- (ncoef/dfe)*sum(res1[1:ncoef]^2)/dfe
##
##  6.  compute fits and residuals
##
# [N, nbasis] = size(datstruct.basismat)
  ind1 <- 1:nbasis
  ind2 <- (nbasis+1):(2*nbasis)
  Ccoef <- coef[ind1]
  Tcoef <- coef[ind2]
#
  phimat <- datstruct$basismat
#
  Chat0 <- phimat %*% Ccoef
  That0 <- phimat %*% Tcoef
#
  yobs <- datstruct$y
#
  Cwt <- as.vector(datstruct$Cwt)
  Twt <- as.vector(datstruct$Twt)
#
#  res = []
#  if fit(1) res = [res; (yobs(:,1) - Chat0)./sqrt(Cwt)]
#  if fit(2)res = [res; (yobs(:,2) - That0)./sqrt(Twt)]
#
  yNames <- c("Conc", "Temp")
  fitNames <- yNames[as.logical(fit)]
  fit12 <- length(fitNames)
  basisNames <- dimnames(phimat)
  res <- array(NA, dim=c(N, fit12), dimnames=
                list(basisNames[[1]], fitNames))
  if(fit[1])res[, 1] <- ((yobs[,1] - Chat0)/sqrt(Cwt))
#  if(fit[2])res[, fit12] <- ((yobs[,fit12] - That0)/sqrt(Twt))
  if(fit[2])res[, fit12] <- ((yobs[,2] - That0)/sqrt(Twt))
#  res[1:5] matches Matlab 2007.05.29
##
## 7.  Derivatives?
##
  Dres <- Dres0
  if( gradwrd){
#
#  7.1.  set up basis and basis derivatve matrices
#
    quadmat <- datstruct$quadbasismat
    Dquadmat <- datstruct$Dquadbasismat
#
#    [nquad, nbasis] = size(quadmat)
#    onesb = ones(1,nbasis)
#    onesq = ones(nquad, 1)
    onesb <- rep(1, nbasis)
    onesq <- rep(1, nquad)
#
#  7.2.  set up some constants that are required
#
    V      <- fitstruct$V
    rho    <- fitstruct$rho
    rhoc   <- fitstruct$rhoc
    delH   <- fitstruct$delH
    Cp     <- fitstruct$Cp
    Cpc    <- fitstruct$Cpc
    Tref   <- fitstruct$Tref
#
#  7.3.  Set up input arrays
#
    F.  <- datstruct$F.
    CA0 <- datstruct$CA0
    T0  <- datstruct$T0
    Tc  <- datstruct$Tcin
    Fc  <- datstruct$Fc
#
#  7.4.  C and T values at fine grid
#
    Chat  <- as.vector(quadmat%*%Ccoef)
    That  <- as.vector(quadmat%*%Tcoef)
    DChat <- as.vector(Dquadmat%*%Ccoef)
    DThat <- as.vector(Dquadmat%*%Tcoef)
#
#  7.5.  betaCC and betaTC depend on kref and Eover R
    Tdif   <- 1/That - 1/Tref
#    temp   = exp(-1e4*EoverR*Tdif)
    log.temp <- (-1e4*EoverR*Tdif)
    oops <- (log.temp > max.log.betaCC)
    if(any(oops)){
      warning(sum(oops), " of ", length(log.temp),
              " values of (-1e4*EoverR*Tdif) exceed the max = ",
              max.log.betaCC, ";  thresholding.")
      log.temp[oops] <- max.log.betaCC
    }
    temp <- exp(log.temp)
#
    betaCC <- kref*temp
    TCfac  <- (-delH/(rho*Cp))
    betaTC <- TCfac*betaCC
#
#  7.6.  betaTT depends on a and b
    Fc2b    <- Fc^b
    aFc2b   <- a*Fc2b
    K1      <- V*rho*Cp
    K2      <- 1./(2.*rhoc*Cpc)
    betaTT  <- Fc*aFc2b/(K1*(Fc + K2*aFc2b))
    betaTT0 <- F./V
#
#  7.7.  compute derivatives of residuals
#
#    %  L values
#
    LC <- DChat + (betaTT0 + betaCC)*Chat - betaTT0*CA0*onesq
    LT <- DThat + ((betaTT0 + betaTT)*That - betaTC*Chat -
                 (betaTT0*T0 + betaTT*Tc))
#
#    %  first order derivatives of L values
#    %  derivatives of L values with respect to
#    %  coefficient vectors c and t
#
    DtbetaCC <- (1e4*EoverR/That^2)*betaCC
    DtbetaTC <- TCfac*DtbetaCC
#
    DcDChat  <- (-(betaCC + betaTT0))
    DtDChat  <- (-DtbetaCC*Chat)
    DcDThat  <-  betaTC
    DtDThat  <- (-(betaTT+betaTT0) + DtbetaTC*Chat)
#
    DcLC   <- (Dquadmat - outer(DcDChat, onesb)*quadmat)
    DtLC   <-          - outer(DtDChat, onesb)*quadmat
    DcLT   <-          - outer(DcDThat, onesb)*quadmat
    DtLT   <- (Dquadmat - outer(DtDThat, onesb)*quadmat)
#
    quadwts    <- datstruct$quadwts
    rootwts    <- sqrt(quadwts)
    quadwtsmat <- outer(quadwts, onesb)
#
#    %  k derivatives
#
    lamC <- lambda[1]
    lamT <- lambda[2]

#  7.8.  assemble the Jacobian matrix

#    DLC <- sqrt(lamC/Cwt).*[DcLC, DtLC]
#    DLT <- sqrt(lamT/Twt).*[DcLT, DtLT]
    DLC <- sqrt(lamC/Cwt)*cbind(DcLC, DtLC)
    DLT <- sqrt(lamT/Twt)*cbind(DcLT, DtLT)
#
#    Jacobian <- [DLC; DLT]
    Jacobian <- rbind(DLC, DLT)
#
#  7.9.  compute derivatives with respect to parameters
#
#    %  set up right hand side of equation D2GDc
#
#    D2GDc <- []
    D2GDc. <- vector('list', 4)
    names(D2GDc.) <- c("kref", "EoverR", "a", "b")
#
#    %  kref
#
    if( estimate[1]) {
#
#        %  first derivative of L values
#
      DkbetaCC <- temp
      DkbetaTC <- TCfac*DkbetaCC
#
      DkLC <-  DkbetaCC*Chat
      DkLT <- -DkbetaTC*Chat
#
#        %  second derivative of L values
#
      DktbetaCC <- (1e4*EoverR/That^2)*temp
      DktbetaTC <- TCfac*DktbetaCC
#
      DkcLC <- outer(  DkbetaCC,   onesb)*quadmat
      DkcLT <- outer( -DkbetaTC, onesb)*quadmat
      DktLC <- outer( DktbetaCC*Chat, onesb)*quadmat
      DktLT <- outer(-DktbetaTC*Chat, onesb)*quadmat
#
#      D2GDck <- zeros(2*nbasis,1)
#
      D2GDck <- rep(0, 2*nbasis)
#        D2GDck(ind1,1) <- (lamC/Cwt).* ...
#            (DcLC'*(DkLC.*quadwts) + DkcLC'*(LC.*quadwts)) + ...
#            (lamT/Twt).* ...
#            (DcLT'*(DkLT.*quadwts) + DkcLT'*(LT.*quadwts))
      D2GDck[ind1] <- ((lamC/Cwt)*
            (t(DcLC)%*%(DkLC*quadwts) + t(DkcLC)%*%(LC*quadwts)) +
            (lamT/Twt)*
            (t(DcLT)%*%(DkLT*quadwts) + t(DkcLT)%*%(LT*quadwts)) )
#        D2GDck(ind2,1) <- (lamC/Cwt).* ...
#            (DtLC'*(DkLC.*quadwts) + DktLC'*(LC.*quadwts)) + ...
#            (lamT/Twt).* ...
#            (DtLT'*(DkLT.*quadwts) + DktLT'*(LT.*quadwts))
      D2GDck[ind2] <- ((lamC/Cwt)*
            (t(DtLC)%*%(DkLC*quadwts) + t(DktLC)%*%(LC*quadwts)) +
            (lamT/Twt)*
            (t(DtLT)%*%(DkLT*quadwts) + t(DktLT)%*%(LT*quadwts)) )
#
#        D2GDc <- [D2GDc, D2GDck]
      D2GDc.$kref <- D2GDck
#
#    end kref
    }

#    %  EoverR
    if(estimate[2]){
#
#        %  first derivative of L values
#
      Dtemp    <- (-1e4*kref*Tdif*temp)
      DEbetaCC <- Dtemp
      DEbetaTC <- TCfac*DEbetaCC
#
      DELC  <-  DEbetaCC*Chat
      DELT  <- (-DEbetaTC*Chat)
#
#      DEtbetaCC <- (1e4.*kref  ./That.^2).* ...
#            (1 - 1e4.*EoverR.*Tdif).*temp
      DEtbetaCC <- ((1e4*kref/That^2)*
            (1 - 1e4*EoverR*Tdif)*temp)
      DEtbetaTC <- TCfac*DEtbetaCC
#
#        DEcLC <- (  DEbetaCC        *onesb).*quadmat
#        DEcLT <- ( -DEbetaTC        *onesb).*quadmat
#        DEtLC <- (( DEtbetaCC.*Chat)*onesb).*quadmat
#        DEtLT <- ((-DEtbetaTC.*Chat)*onesb).*quadmat
      DEcLC <- outer(  DEbetaCC, onesb)*quadmat
      DEcLT <- outer( -DEbetaTC, onesb)*quadmat
      DEtLC <- outer(( DEtbetaCC*Chat), onesb)*quadmat
      DEtLT <- outer((-DEtbetaTC*Chat), onesb)*quadmat
#
#        D2GDcE <- zeros(2*nbasis,1)
      D2GDcE <- rep(0, 2*nbasis)
#        D2GDcE(ind1,1) <- (lamC/Cwt).* ...
#            (DcLC'*(DELC.*quadwts) + DEcLC'*(LC.*quadwts)) + ...
#            (lamT/Twt).* ...
#            (DcLT'*(DELT.*quadwts) + DEcLT'*(LT.*quadwts))
      DcLC.. <- (crossprod(DcLC, DELC*quadwts) +
                  crossprod(DEcLC, LC*quadwts))
      DcLT.. <- (crossprod(DcLT, DELT*quadwts) +
                  crossprod(DEcLT, LT*quadwts))
      D2GDcE[ind1] <- ((lamC/Cwt)* DcLC.. + (lamT/Twt)* DcLT..)
#        D2GDcE(ind2,1) <- (lamC./Cwt).* ...
#            (DtLC'*(DELC.*quadwts) + DEtLC'*(LC.*quadwts)) + ...
#            (lamT./Twt).* ...
#            (DtLT'*(DELT.*quadwts) + DEtLT'*(LT.*quadwts))
      DtLC.. <- (crossprod(DtLC, DELC*quadwts) +
                  crossprod(DEtLC, LC*quadwts))
      DtLT.. <- (crossprod(DtLT, DELT*quadwts) +
                  crossprod(DEtLT, LT*quadwts))
      D2GDcE[ind2] <- ((lamC/Cwt)* DtLC.. + (lamT/Twt)* DtLT..)
#        D2GDc <- [D2GDc, D2GDcE]
      D2GDc.$EoverR <- D2GDcE
#
#    end EoverR
    }

#    %  a
    if(estimate[3]){
#
#        %  first derivative of L values
#
#        DhbetaTT <- (betaTT./aFc2b).*(1 - K1.*K2.*betaTT./Fc)
#        DabetaTT <- DhbetaTT.*Fc2b
      DhbetaTT <- ((betaTT/aFc2b)*(1 - K1*K2*betaTT/Fc))
      DabetaTT <- DhbetaTT*Fc2b
#
      DaLT <- (DabetaTT*(That - Tc))
#
      DatLT <- outer(DabetaTT, onesb)*quadmat
#
#        D2GDca <- zeros(2*nbasis,1)
      D2GDca <- rep(0, 2*nbasis)
#
#        D2GDca(ind1,1) <- (lamT/Twt).*(DcLT'*(DaLT.*quadwts))
#        D2GDca(ind2,1) <- (lamT./Twt).* ...
#            (DtLT'*(DaLT.*quadwts) + DatLT'*(LT.*quadwts))
      D2GDca[ind1] <- ((lamT/Twt)*(t(DcLT) %*%(DaLT*quadwts)))
      D2GDca[ind2] <- ((lamT/Twt)*
            (t(DtLT)%*%(DaLT*quadwts) + t(DatLT)%*%(LT*quadwts)) )
#
#        D2GDc <- [D2GDc, D2GDca]
      D2GDc.$a <- D2GDca
#
#    end 'a'
    }

    if( estimate[4]){
#
#        %  b derivative of L values
#
#        DhbetaTT <- (betaTT./aFc2b).*(1 - K1.*K2.*betaTT./Fc)
#        DbbetaTT <- DhbetaTT.*b.*aFc2b./Fc
      DhbetaTT <- (betaTT/aFc2b)*(1 - K1*K2*betaTT/Fc)
      DbbetaTT <- DhbetaTT*b*aFc2b/Fc
#
      DbLT <- DbbetaTT*(That - Tc)
#
#        DbtLT <- (DbbetaTT*onesb).*quadmat
      DbtLT <- outer(DbbetaTT, onesb)*quadmat
#
#        D2GDcb <- zeros(2*nbasis,1)
#        D2GDcb(ind1,1) <- (lamT/Twt).*(DcLT'*(DbLT.*quadwts))
#        D2GDcb(ind2,1) <- (lamT./Twt).* ...
#            (DtLT'*(DbLT.*quadwts) + DbtLT'*(LT.*quadwts))
#
      D2GDcb <- rep(0, 2*nbasis)
      D2GDcb[ind1] <- ((lamT/Twt)*(t(DcLT)%*%(DbLT*quadwts)) )
      D2GDcb[ind2] <- ((lamT/Twt)*
            (t(DtLT)%*%(DbLT*quadwts) + t(DbtLT)%*%(LT*quadwts)) )
#
#        D2GDc <- [D2GDc, D2GDcb]
      D2GDc.$b <- D2GDcb
#
#    end 'b'
    }
#   Convert from a list to a matrix,
#   dropping columns not in estimate
    D2GDc <- do.call(cbind, D2GDc.)
##
## 8.  Construct D2GDc2
##
#  8.1.  First part
#    Wmat <- [quadwtsmat, quadwtsmat; quadwtsmat, quadwtsmat]
    W.5 <- cbind(quadwtsmat, quadwtsmat)
    Wmat <- rbind(W.5, W.5)
#
#    D2GDc2  <- (Jacobian.*Wmat)'*Jacobian
    D2GDc2  <- crossprod(Jacobian*Wmat, Jacobian)
#
#    ZtZmat <- phimat'*phimat
    ZtZmat <- crossprod(phimat)

    if( fit[1])
        D2GDc2[ind1,ind1] <- D2GDc2[ind1,ind1] + ZtZmat/Cwt
#    end
    if( fit[2])
        D2GDc2[ind2,ind2] <- D2GDc2[ind2,ind2] + ZtZmat/Twt
#    end

#  8.2.  Add second derivative information

    DttbetaCC <- ((1e4*kref*EoverR/That^2)*
                 (1e4*EoverR/That^2 - 2/That)*temp)
    DttbetaTC <- TCfac*DttbetaCC
#
#    DctLC <- sparse(zeros(nbasis,nbasis))
#    DttLC <- sparse(zeros(nbasis,nbasis))
#    DctLT <- sparse(zeros(nbasis,nbasis))
#    DttLT <- sparse(zeros(nbasis,nbasis))
    DctLC <- array(0, dim=c(nbasis, nbasis))
    DttLT <- DctLT <- DttLC <- DctLC

#    norder <- nbasis - length(getbasispar(CSTRbasis))
## *** 'getbasispar' <- interior knots
    norder <- nbasis - length(CSTRbasis$params)
#
    for( i in 1:nbasis){
      jstart <- max(c(1,i-norder+1))
      for( j in jstart:i) {
        qijvec <- quadmat[,i]*quadmat[,j]*quadwts
        DctLC[i,j] <- sum(qijvec*LC*DtbetaCC)
        DttLC[i,j] <- sum(qijvec*LC*DttbetaCC*Chat)
        DctLT[i,j] <- sum(qijvec*LT*DtbetaTC)
        DttLT[i,j] <- sum(qijvec*LT*DttbetaTC*Chat)
        if( i != j){
          DctLC[j,i] <- DctLC[i,j]
          DttLC[j,i] <- DttLC[i,j]
          DctLT[j,i] <- DctLT[i,j]
          DttLT[j,i] <- DttLT[i,j]
#       end if(i!=j)
        }
#    end for(j in jstart:i)
      }
#   end for(i in 1:nbasis)
    }
#    DctL <- lamC.*DctLC./Cwt + lamT.*DctLT./Twt
    DctL <- lamC*DctLC/Cwt + lamT*DctLT/Twt
#      DttL <- lamC.*DttLC./Cwt + lamT.*DttLT./Twt
    DttL <- lamC*DttLC/Cwt + lamT*DttLT/Twt

#  8.3.  modify D2GDc2

    D2GDc2[ind1,ind2] <- D2GDc2[ind1,ind2] + DctL
    D2GDc2[ind2,ind1] <- D2GDc2[ind2,ind1] + t(DctL)
    D2GDc2[ind2,ind2] <- D2GDc2[ind2,ind2] + DttL

#  8.4.  compute (D2GDc2)^{-1} D2GDc

#    DcDtheta <- D2GDc2\D2GDc
    DcDtheta <- try(solve(D2GDc2, D2GDc))
    if(inherits(DcDtheta,"try-error")) {
#
      D2GDc2.eig <- eigen(D2GDc2, symmetric=TRUE)
      Dc.ev <- D2GDc2.eig$values
      Dc.rank <- sum(Dc.ev>0)
#
      if(Dc.rank<ncoef)
        warning("D2GDc2 has reduced rank ", Dc.rank, "; using ginverse.")
      Dc.rank1 <- sum(Dc.ev > (eps*Dc.ev[1]))
      if(Dc.rank1 < Dc.rank)
        warning("D2GDc2 is ill conditioned.  Reducing rank to ",
                Dc.rank1, " from ", Dc.rank)
      jrank <- 1:Dc.rank1
      DcDtheta <- with(D2GDc2.eig, (vectors[, jrank] /
          rep(Dc.ev[jrank], each=ncoef)) %*% crossprod(vectors[, jrank], D2GDc))
    }
#
#  8.5.  set up Dres
#
#    Dres <- []
    Dres <- NULL
#
    if( fit[1])
      Dres <- phimat%*%DcDtheta[ind1,]/sqrt(Cwt)
#       end
    if( fit[2])
      Dres <- rbind(Dres, phimat%*%DcDtheta[ind2,]/sqrt(Twt))
  }
##
## 9.  Done
##
#  list(res=res, Dres=Dres, fitstruct=fitstruct, df=df, lambda=lambda, gradwrd=gradwrd)
  list(res=res, Dres=Dres, fitstruct=fitstruct, df=df., gcv=gcv)
}