R/CSTRfn.R In fda: Functional Data Analysis

Documented in CSTRfn

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

#%  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")

#[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]
n.Sres <- N*k12
Nval <- n.Sres + n.Lres
#
nbasis <- (dim(CSTR0\$Dres\$DSres)[2]/2)
Dres0 <- with(CSTR0\$Dres, rbind(DSres, DLres))
# 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.

# 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
r.D <- crossprod(res0, Dres0)
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)))
#
##
## 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 <- svd2(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)
}
coef0 <- coef1
res0 <- res1
Dres0 <- Dres1
F0 <- F1
# end while(( gradnorm1 > tolval) | (fundif > tolval)){
}
##
##% 4.  update coef
##
coef <- coef0
fitstruct\$coef0 <- coef
#
# 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 <- svd2(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
#
#  7.1.  set up basis and basis derivatve matrices
#
#
#    onesb = ones(1,nbasis)
onesb <- rep(1, nbasis)
#
#  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
#
#
#  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)
#
#
#
#    %  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
#
#
#      D2GDck <- zeros(2*nbasis,1)
#
D2GDck <- rep(0, 2*nbasis)
#        D2GDck(ind1,1) <- (lamC/Cwt).* ...
#            (lamT/Twt).* ...
D2GDck[ind1] <- ((lamC/Cwt)*
(lamT/Twt)*
#        D2GDck(ind2,1) <- (lamC/Cwt).* ...
#            (lamT/Twt).* ...
D2GDck[ind2] <- ((lamC/Cwt)*
(lamT/Twt)*
#
#        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
#
#        D2GDcE <- zeros(2*nbasis,1)
D2GDcE <- rep(0, 2*nbasis)
#        D2GDcE(ind1,1) <- (lamC/Cwt).* ...
#            (lamT/Twt).* ...
D2GDcE[ind1] <- ((lamC/Cwt)* DcLC.. + (lamT/Twt)* DcLT..)
#        D2GDcE(ind2,1) <- (lamC./Cwt).* ...
#            (lamT./Twt).* ...
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))
#
#
#        D2GDca <- zeros(2*nbasis,1)
D2GDca <- rep(0, 2*nbasis)
#
#        D2GDca(ind2,1) <- (lamT./Twt).* ...
D2GDca[ind2] <- ((lamT/Twt)*
#
#        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)
#
#
#        D2GDcb <- zeros(2*nbasis,1)
#        D2GDcb(ind2,1) <- (lamT./Twt).* ...
#
D2GDcb <- rep(0, 2*nbasis)
D2GDcb[ind2] <- ((lamT/Twt)*
#
#        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 <- 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) {
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(class(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
##