1## Calculate square error projection from sscden objects
2project.sscden1 <- function(object,include,...)
3{
4    include <- unique(include)
5    term <- object$terms
6    mf <- object$mf
7    nobs <- dim(mf)[1]
8    xnames <- object$xnames
9    ynames <- object$ynames
10    xx.wt <- object$xx.wt
11    nx <- length(xx.wt)
12    xx <- mf[!object$x.dup.ind,xnames,drop=FALSE]
13    qd.pt <- object$rho$env$qd.pt
14    qd.wt <- object$rho$env$qd.wt
15    rho.d <- t(object$rho$fun(xx,qd.pt,object$rho$env,outer=TRUE))
16    rho.wt <- rho.d*qd.wt
17    rho.d <- t(t(log(rho.d))-apply(log(rho.d)*rho.wt,2,sum))
18    nmesh <- length(qd.wt)
19    ns <- length(object$id.s)
20    nr <- length(object$id.r)
21    nbasis <- length(object$id.basis)
22    s.rho <- ss <- 0
23    r.rho <- matrix(0,nbasis,nr)
24    sr <- array(0,c(ns,nbasis,nr))
25    rr <- array(0,c(nbasis,nbasis,nr,nr))
26    rho2 <- sum(xx.wt*apply(rho.d^2*rho.wt,2,sum))
27    for (k in 1:nx) {
28        id.x <- (1:nobs)[!object$x.dup.ind][k]
29        qd.s <- NULL
30        qd.r <- list(NULL)
31        iq <- 0
32        for (label in term$labels) {
33            vlist <- term[[label]]$vlist
34            x.list <- xnames[xnames%in%vlist]
35            y.list <- ynames[ynames%in%vlist]
36            xy.basis <- mf[object$id.basis,vlist]
37            qd.xy <- data.frame(matrix(0,nmesh,length(vlist)))
38            names(qd.xy) <- vlist
39            qd.xy[,y.list] <- qd.pt[,y.list]
40            if (length(x.list)) xx <- mf[rep(id.x,nmesh),x.list,drop=FALSE]
41            else xx <- NULL
42            nphi <- term[[label]]$nphi
43            nrk <- term[[label]]$nrk
44            if (nphi) {
45                phi <- term[[label]]$phi
46                for (i in 1:nphi) {
47                    if (is.null(xx)) {
48                        qd.s.wk <- phi$fun(qd.xy[,,drop=TRUE],nu=i,env=phi$env)
49                        qd.s <- cbind(qd.s,qd.s.wk)
50                    }
51                    else {
52                        if (length(y.list)>0) {
53                            qd.xy[,x.list] <- xx
54                            qd.s <- cbind(qd.s,phi$fun(qd.xy,i,phi$env))
55                        }
56                    }
57                }
58            }
59            if (nrk) {
60                rk <- term[[label]]$rk
61                for (i in 1:nrk) {
62                    if (is.null(xx)) {
63                        qd.r.wk <- rk$fun(qd.xy[,,drop=TRUE],xy.basis,nu=i,env=rk$env,TRUE)
64                        iq <- iq+1
65                        qd.r[[iq]] <- qd.r.wk
66                    }
67                    else {
68                        if (length(y.list)>0) {
69                            qd.xy[,x.list] <- xx
70                            iq <- iq+1
71                            qd.r[[iq]] <- rk$fun(qd.xy,xy.basis,i,rk$env,TRUE)
72                        }
73                    }
74                }
75            }
76        }
77        if (ns) {
78            qd.s <- sweep(qd.s,2,apply(qd.s*rho.wt[,k],2,sum))
79            s.rho <- s.rho + xx.wt[k]*apply(qd.s*rho.d[,k]*rho.wt[,k],2,sum)
80            ss <- ss + xx.wt[k]*t(rho.wt[,k]*qd.s)%*%qd.s
81        }
82        for (i in 1:iq) {
83            qd.r[[i]] <- sweep(qd.r[[i]],2,apply(qd.r[[i]]*rho.wt[,k],2,sum))
84            r.rho[,i] <- r.rho[,i] + xx.wt[k]*apply(qd.r[[i]]*rho.d[,k]*rho.wt[,k],2,sum)
85            if (ns) sr[,,i] <- sr[,,i] + xx.wt[k]*t(rho.wt[,k]*qd.s)%*%qd.r[[i]]
86            for (j in 1:i) {
87                rr.wk <- xx.wt[k]*t(rho.wt[,k]*qd.r[[i]])%*%qd.r[[j]]
88                rr[,,i,j] <- rr[,,i,j] + rr.wk
89                if (i-j) rr[,,j,i] <- rr[,,j,i] + t(rr.wk)
90            }
91        }
92    }
93    ## evaluate full model
94    if (ns) d <- object$d[object$id.s]
95    c <- object$c
96    theta <- object$theta[object$id.r]
97    nq <- length(theta)
98    if (ns) s.eta <- ss%*%d
99    r.eta <- tmp <- NULL
100    r.rho.wk <- sr.wk <- rr.wk <- 0
101    for (i in 1:nq) {
102        tmp <- c(tmp,10^(2*theta[i])*sum(diag(rr[,,i,i])))
103        if (ns) {
104            s.eta <- s.eta + 10^theta[i]*sr[,,i]%*%c
105            if (length(d)==1) r.eta.wk <- sr[,,i]*d
106            else r.eta.wk <- t(sr[,,i])%*%d
107            sr.wk <- sr.wk + 10^theta[i]*sr[,,i]
108        }
109        else r.eta.wk <- 0
110        r.rho.wk <- r.rho.wk + 10^theta[i]*r.rho[,i]
111        for (j in 1:nq) {
112            r.eta.wk <- r.eta.wk + 10^theta[j]*rr[,,i,j]%*%c
113            rr.wk <- rr.wk + 10^(theta[i]+theta[j])*rr[,,i,j]
114        }
115        r.eta <- cbind(r.eta,r.eta.wk)
116    }
117    rho.eta <- sum(r.rho.wk*c)
118    if (ns) rho.eta <- rho.eta + sum(r.rho.wk*c)
119    eta2 <- sum(c*(rr.wk%*%c))
120    if (ns) eta2 <- eta2 + sum(d*(ss%*%d)) + 2*sum(d*(sr.wk%*%c))
121    mse <- eta2 + rho2 + 2*rho.eta
122    ## extract terms in subspace
123    id.s <- id.r <- NULL
124    for (label in include) {
125        id.s <- c(id.s,object$id.s.list[[label]])
126        id.r <- c(id.r,object$id.r.list[[label]])
127    }
128    if (is.null(id.s)&is.null(id.r))
129        stop("gss error in project.sscden1: include some terms")
130    if (!all(id.s%in%object$id.s)|!all(id.r%in%object$id.r))
131        stop("gss error in project.sscden1: included terms are not in the model")
132    ## calculate projection
133    rkl <- function(theta1=NULL) {
134        theta.wk <- 1:nq0
135        theta.wk[fix] <- theta[fix]
136        if (nq0-1) theta.wk[-fix] <- theta1
137        ##
138        id.s0 <- (1:length(object$id.s))[object$id.s%in%id.s]
139        id.r0 <- (1:length(object$id.r))[object$id.r%in%id.r]
140        if (length(id.s0)) ss.wk <- ss[id.s0,id.s0,drop=FALSE]
141        if (length(id.r0)) {
142            r.eta.wk <- rr.wk <- 0
143            sr.wk <- matrix(0,length(id.s),nbasis)
144            for (i in 1:length(id.r0)) {
145                r.eta.wk <- r.eta.wk + 10^theta.wk[i]*r.eta[,id.r0[i]]
146                if (length(id.s0)) sr.wk <- sr.wk + 10^theta.wk[i]*sr[id.s0,,id.r0[i]]
147                for (j in 1:length(id.r0)) {
148                    rr.wk <- rr.wk + 10^(theta.wk[i]+theta.wk[j])*rr[,,id.r0[i],id.r0[j]]
149                }
150            }
151            if (length(id.s0)) {
152                v <- cbind(rbind(ss.wk,t(sr.wk)),rbind(sr.wk,rr.wk))
153                mu <- c(s.eta[id.s0],r.eta.wk)
154            }
155            else {
156                v <- rbind(sr.wk,rr.wk)
157                mu <- r.eta.wk
158            }
159        }
160        else {
161            v <- ss.wk
162            mu <- s.eta[id.s0]
163        }
164        nn <- length(mu)
165        suppressWarnings(z <- chol(v,pivot=TRUE))
166        v <- z
167        rkv <- attr(z,"rank")
168        m.eps <- .Machine$double.eps
169        while (v[rkv,rkv]<2*sqrt(m.eps)*v[1,1]) rkv <- rkv - 1
170        if (rkv<nn) v[(1:nn)>rkv,(1:nn)>rkv] <- diag(v[1,1],nn-rkv)
171        mu <- backsolve(v,mu[attr(z,"pivot")],transpose=TRUE)
172        eta2 - sum(mu[1:rkv]^2)
173    }
174    cv.wk <- function(theta) cv.scale*rkl(theta)+cv.shift
175    ## initialization
176    fix <- rev(order(tmp[id.r]))[1]
177    theta <- object$theta[id.r]
178    ## projection
179    nq0 <- length(id.r)
180    if (nq0>1) {
181        if (object$skip.iter) se <- rkl(theta[-fix])
182        else {
183            if (nq0-2) {
184                ## scale and shift cv
185                tmp <- abs(rkl(theta[-fix]))
186                cv.scale <- 1
187                cv.shift <- 0
188                if (tmp<1&tmp>10^(-4)) {
189                    cv.scale <- 10/tmp
190                    cv.shift <- 0
191                }
192                if (tmp<10^(-4)) {
193                    cv.scale <- 10^2
194                    cv.shift <- 10
195                }
196                zz <- nlm(cv.wk,theta[-fix],stepmax=.5,ndigit=7)
197            }
198            else {
199                the.wk <- theta[-fix]
200                repeat {
201                    mn <- the.wk-1
202                    mx <- the.wk+1
203                    zz <- nlm0(rkl,c(mn,mx))
204                    if (min(zz$est-mn,mx-zz$est)>=1e-3) break
205                    else the.wk <- zz$est
206                }
207            }
208            se <- rkl(zz$est)
209        }
210    }
211    else se <- rkl()
212    list(ratio=se/mse,se=se)
213}
214