1## Calculate Kullback-Leibler projection from ssden objects
2project.ssden <- function(object,include,mesh=FALSE,...)
3{
4    qd.pt <- object$quad$pt
5    qd.wt <- object$quad$wt
6    bias <- object$bias
7    ## evaluate full model
8    mesh0 <- dssden(object,qd.pt)
9    qd.wt <- qd.wt*bias$qd.wt
10    qd.wt <- t(t(qd.wt)/apply(qd.wt*mesh0,2,sum))
11    ## extract terms in subspace
12    nqd <- dim(qd.wt)[1]
13    nxi <- length(object$id.basis)
14    qd.s <- qd.r <- q <- NULL
15    theta <- d <- NULL
16    n0.wk <- nq.wk <- nq <- 0
17    for (label in object$terms$labels) {
18        x.basis <- object$mf[object$id.basis,object$term[[label]]$vlist]
19        qd.x <- qd.pt[,object$term[[label]]$vlist]
20        nphi <- object$term[[label]]$nphi
21        nrk <- object$term[[label]]$nrk
22        if (nphi) {
23            phi <- object$term[[label]]$phi
24            for (i in 1:nphi) {
25                n0.wk <- n0.wk + 1
26                if (!any(label==include)) next
27                d <- c(d,object$d[n0.wk])
28                qd.s <- cbind(qd.s,phi$fun(qd.x,nu=i,env=phi$env))
29            }
30        }
31        if (nrk) {
32            rk <- object$term[[label]]$rk
33            for (i in 1:nrk) {
34                nq.wk <- nq.wk + 1
35                if (!any(label==include)) next
36                nq <- nq + 1
37                theta <- c(theta,object$theta[nq.wk])
38                qd.r <- array(c(qd.r,rk$fun(x.basis,qd.x,nu=i,env=rk$env,out=TRUE)),
39                              c(nxi,nqd,nq))
40                q <- cbind(q,rk$fun(x.basis,x.basis,nu=i,env=rk$env,out=FALSE))
41            }
42        }
43    }
44    if (is.null(qd.s)&is.null(qd.r))
45        stop("gss error in project.ssden: include some terms")
46    if (!is.null(qd.s)) {
47        nn <- nxi + ncol(qd.s)
48        qd.s <- t(qd.s)
49    }
50    else nn <- nxi
51    ## calculate projection
52    rkl <- function(theta1=NULL) {
53        theta.wk <- 1:nq
54        theta.wk[fix] <- theta[fix]
55        if (nq-1) theta.wk[-fix] <- theta1
56        qd.rs <- 0
57        for (i in 1:nq) qd.rs <- qd.rs + 10^theta.wk[i]*qd.r[,,i]
58        qd.rs <- rbind(qd.rs,qd.s)
59        z <- .Fortran("drkl",
60                      cd=as.double(cd), as.integer(nn),
61                      as.double(t(qd.rs)), as.integer(nqd), as.integer(bias$nt),
62                      as.double(bias$wt), as.double(t(qd.wt)),
63                      mesh=as.double(mesh0), as.double(.Machine$double.eps),
64                      as.double(1e-6), as.integer(30), integer(nn),
65                      double(2*bias$nt*(nqd+1)+nn*(2*nn+4)), info=integer(1),
66                      PACKAGE="gss")
67        if (z$info==1)
68            stop("gss error in project.ssden: Newton iteration diverges")
69        if (z$info==2)
70            warning("gss warning in project.ssden: Newton iteration fails to converge")
71        assign("cd",z$cd,inherits=TRUE)
72        assign("mesh1",z$mesh,inherits=TRUE)
73        sum(bias$wt*(apply(qd.wt*log(mesh0/mesh1)*mesh0,2,sum)+
74                     log(apply(qd.wt*mesh1,2,sum))))
75    }
76    cv.wk <- function(theta) cv.scale*rkl(theta)+cv.shift
77    if (nq) {
78        ## initialization
79        if (is.null(qd.s)) theta.wk <- 0
80        else {
81            qd.r.wk <- 0
82            for (i in 1:nq) qd.r.wk <- qd.r.wk + 10^theta[i]*qd.r[,,i]
83            vv.s <- vv.r <- 0
84            for (i in 1:bias$nt) {
85                mu.s <- apply(qd.wt[,i]*qd.s,2,sum)/sum(qd.wt[,i])
86                v.s <- apply(qd.wt[,i]*qd.s^2,2,sum)/sum(qd.wt[,i])
87                v.s <- v.s - mu.s^2
88                mu.r <- apply(qd.wt[,i]*qd.r.wk,2,sum)/sum(qd.wt[,i])
89                v.r <- apply(qd.wt[,i]*qd.r.wk^2,2,sum)/sum(qd.wt[,i])
90                v.r <- v.r - mu.r^2
91                vv.s <- vv.s + bias$wt[i]*v.s
92                vv.r <- vv.r + bias$wt[i]*v.r
93            }
94            theta.wk <- log10(sum(vv.s)/(nn-nxi)/sum(vv.r)*nxi) / 2
95        }
96        theta <- theta + theta.wk
97        tmp <- NULL
98        for (i in 1:nq) tmp <- c(tmp,10^theta[i]*sum(q[,i]))
99        fix <- rev(order(tmp))[1]
100        ## projection
101        cd <- c(10^(-theta.wk)*object$c,d)
102        mesh1 <- NULL
103        if (nq-1) {
104            if (object$skip.iter) kl <- rkl(theta[-fix])
105            else {
106                if (nq-2) {
107                    ## scale and shift cv
108                    tmp <- abs(rkl(theta[-fix]))
109                    cv.scale <- 1
110                    cv.shift <- 0
111                    if (tmp<1&tmp>10^(-4)) {
112                        cv.scale <- 10/tmp
113                        cv.shift <- 0
114                    }
115                    if (tmp<10^(-4)) {
116                        cv.scale <- 10^2
117                        cv.shift <- 10
118                    }
119                    zz <- nlm(cv.wk,theta[-fix],stepmax=.5,ndigit=7)
120                }
121                else {
122                    the.wk <- theta[-fix]
123                    repeat {
124                        mn <- the.wk-1
125                        mx <- the.wk+1
126                        zz <- nlm0(rkl,c(mn,mx))
127                        if (min(zz$est-mn,mx-zz$est)>=1e-3) break
128                        else the.wk <- zz$est
129                    }
130                }
131                kl <- rkl(zz$est)
132            }
133        }
134        else kl <- rkl()
135    }
136    else {
137        nn <- nrow(qd.s)
138        z <- .Fortran("drkl",
139                      cd=as.double(d), as.integer(nn),
140                      as.double(qd.s), as.integer(nqd), as.integer(bias$nt),
141                      as.double(bias$wt), as.double(t(qd.wt)),
142                      mesh=as.double(mesh0), as.double(.Machine$double.eps),
143                      as.double(1e-6), as.integer(30), integer(nn),
144                      double(2*bias$nt*(nqd+1)+nn*(2*nn+4)), info=integer(1),
145                      PACKAGE="gss")
146        if (z$info==1)
147            stop("gss error in project.ssden: Newton iteration diverges")
148        if (z$info==2)
149            warning("gss warning in project.ssden: Newton iteration fails to converge")
150        mesh1 <- z$mesh
151        kl <- sum(bias$wt*(apply(qd.wt*log(mesh0/mesh1)*mesh0,2,sum)+
152                           log(apply(qd.wt*mesh1,2,sum))))
153    }
154    kl0 <- sum(bias$wt*(apply(qd.wt*log(mesh0)*mesh0,2,sum)+
155                        log(apply(qd.wt,2,sum))))
156    kl <- sum(bias$wt*(apply(qd.wt*log(mesh0/mesh1)*mesh0,2,sum)+
157                       log(apply(qd.wt*mesh1,2,sum))))
158    wt.wk <- t(t(qd.wt)/apply(qd.wt*mesh1,2,sum))
159    kl1 <- sum(bias$wt*(apply(wt.wk*log(mesh1)*mesh1,2,sum)+
160                        log(apply(wt.wk,2,sum))))
161    obj <- list(ratio=kl/kl0,kl=kl,check=(kl+kl1)/kl0)
162    if (mesh) obj$mesh <- mesh1
163    obj
164}
165