1## Make RK for thin-plate splines
2mkrk.tp <- function(dm,order,mesh,weight=1)
3{
4    ## Check inputs
5    if (!((2*order>dm)&(dm>=1))) {
6        stop("gss error: thin-plate spline undefined for the parameters")
7    }
8    if (xor(is.vector(mesh),dm==1)
9        |xor(is.matrix(mesh),dm>=2)) {
10        stop("gss error in mkrk.tp: mismatched inputs")
11    }
12    if ((min(weight)<0)|(max(weight)<=0)) {
13        stop("gss error in mkrk.tp: negative weights")
14    }
15    ## Set weights
16    if (is.vector(mesh)) N <- length(mesh)
17    else N <- dim(mesh)[1]
18    weight <- rep(weight,len=N)
19    weight <- sqrt(weight/sum(weight))
20    ## Obtain orthonormal basis
21    phi.p <- mkphi.tp.p(dm,order)
22    nnull <- choose(dm+order-1,dm)
23    s <- NULL
24    for (nu in 1:nnull) s <- cbind(s,phi.p$fun(mesh,nu,phi.p$env))
25    s <- qr(weight*s)
26    if (s$rank<nnull) {
27        stop("gss error in mkrk.tp: insufficient normalizing mesh for thin-plate spline")
28    }
29    q <- qr.Q(s)
30    r <- qr.R(s)
31    ## Set Q^{T}E(|u_{i}-u_{j}|)Q
32    rk.p <- mkrk.tp.p(dm,order)
33    pep <- weight*t(weight*rk.p$fun(mesh,mesh,rk.p$env,out=TRUE))
34    pep <- t(q)%*%pep%*%q
35    ## Create the environment
36    env <- list(dim=dm,order=order,weight=weight,
37                phi.p=phi.p,rk.p=rk.p,q=q,r=r,mesh=mesh,pep=pep)
38    ## Create the rk function
39    fun <- function(x,y,env,outer.prod=FALSE) {
40        ## Check inputs
41        if (env$dim==1) {
42            if (!(is.vector(x)&is.vector(y))) {
43                stop("gss error in rk: inputs are of wrong types")
44            }
45            nx <- length(x)
46            ny <- length(y)
47        }
48        else {
49            if (is.vector(x)) x <- t(as.matrix(x))
50            if (env$dim!=dim(x)[2]) {
51                stop("gss error in rk: inputs are of wrong dimensions")
52            }
53            nx <- dim(x)[1]
54            if (is.vector(y)) y <- t(as.matrix(y))
55            if (env$dim!=dim(y)[2]) {
56                stop("gss error in rk: inputs are of wrong dimensions")
57            }
58            ny <- dim(y)[1]
59        }
60        ## Return the results
61        nnull <- choose(env$dim+env$order-1,env$dim)
62        if (outer.prod) {
63            phix <- phiy <- NULL
64            for (nu in 1:nnull) {
65                phix <- rbind(phix,env$phi.p$fun(x,nu,env$phi.p$env))
66                phiy <- rbind(phiy,env$phi.p$fun(y,nu,env$phi.p$env))
67            }
68            phix <- backsolve(env$r,phix,transpose=TRUE)
69            phiy <- backsolve(env$r,phiy,transpose=TRUE)
70            ex <- env$rk.p$fun(env$mesh,x,env$rk.p$env,out=TRUE)
71            ex <- env$weight*ex
72            ex <- t(env$q)%*%ex
73            ey <- env$rk.p$fun(env$mesh,y,env$rk.p$env,out=TRUE)
74            ey <- env$weight*ey
75            ey <- t(env$q)%*%ey
76            env$rk.p$fun(x,y,env$rk.p$env,out=TRUE)-t(phix)%*%ey-
77                t(ex)%*%phiy+t(phix)%*%env$pep%*%phiy
78        }
79        else {
80            N <- max(nx,ny)
81            phix <- phiy <- NULL
82            for (nu in 1:nnull) {
83                phix <- rbind(phix,env$phi.p$fun(x,nu,env$phi.p$env))
84                phiy <- rbind(phiy,env$phi.p$fun(y,nu,env$phi.p$env))
85            }
86            phix <- backsolve(env$r,phix,transpose=TRUE)
87            phix <- matrix(phix,nnull,N)
88            phiy <- backsolve(env$r,phiy,transpose=TRUE)
89            phiy <- matrix(phiy,nnull,N)
90            ex <- env$rk.p$fun(env$mesh,x,env$rk.p$env,out=TRUE)
91            ex <- env$weight*ex
92            ex <- t(env$q)%*%ex
93            ex <- matrix(ex,nnull,N)
94            ey <- env$rk.p$fun(env$mesh,y,env$rk.p$env,out=TRUE)
95            ey <- env$weight*ey
96            ey <- t(env$q)%*%ey
97            ey <- matrix(ey,nnull,N)
98            fn1 <- function(x,n) x[1:n]%*%x[n+(1:n)]
99            fn2 <- function(x,pep,n) t(x[1:n])%*%pep%*%x[n+(1:n)]
100            env$rk.p$fun(x,y,env$rk.p$env)-apply(rbind(phix,ey),2,fn1,nnull)-
101                apply(rbind(phiy,ex),2,fn1,nnull)+
102                    apply(rbind(phix,phiy),2,fn2,env$pep,nnull)
103        }
104    }
105    ## Return the function and the environment
106    list(fun=fun,env=env)
107}
108
109## Make phi function for thin-plate splines
110mkphi.tp <-  function(dm,order,mesh,weight)
111{
112    ## Check inputs
113    if (!((2*order>dm)&(dm>=1))) {
114        stop("gss error: thin-plate spline undefined for the parameters")
115    }
116    if (xor(is.vector(mesh),dm==1)
117        |xor(is.matrix(mesh),dm>=2)) {
118        stop("gss error in mkphi.tp: mismatched inputs")
119    }
120    if ((min(weight)<0)|(max(weight)<=0)) {
121        stop("gss error in mkphi.tp: negative weights")
122    }
123    ## Set weights
124    if (is.vector(mesh)) N <- length(mesh)
125    else N <- dim(mesh)[1]
126    weight <- rep(weight,len=N)
127    weight <- sqrt(weight/sum(weight))
128    ## Create the environment
129    phi.p <- mkphi.tp.p(dm,order)
130    nnull <- choose(dm+order-1,dm)
131    s <- NULL
132    for (nu in 1:nnull) s <- cbind(s,phi.p$fun(mesh,nu,phi.p$env))
133    s <- qr(weight*s)
134    if (s$rank<nnull) {
135        stop("gss error in mkphi: insufficient normalizing mesh for thin-plate spline")
136    }
137    r <- qr.R(s)
138    env <- list(dim=dm,order=order,phi.p=phi.p,r=r)
139    ## Create the phi function
140    fun <- function(x,nu,env) {
141        nnull <- choose(env$dim+env$order-1,env$dim)
142        phix <- NULL
143        for(i in 1:nnull)
144            phix <- rbind(phix,env$phi.p$fun(x,i,env$phi.p$env))
145        t(backsolve(env$r,phix,transpose=TRUE))[,nu+1]
146    }
147    ## Return the function and the environment
148    list(fun=fun,env=env)
149}
150
151## Make pseudo RK for thin-plate splines
152mkrk.tp.p <- function(dm,order)
153{
154    ## Check inputs
155    if (!((2*order>dm)&(dm>=1))) {
156        stop("gss error: thin-plate spline undefined for the parameters")
157    }
158    ## Create the environment
159    if (dm%%2) {
160        theta <- gamma(dm/2-order)/2^(2*order)/pi^(dm/2)/gamma(order)
161    }
162    else {
163        theta <- (-1)^(dm/2+order+1)/2^(2*order-1)/pi^(dm/2)/
164            gamma(order)/gamma(order-dm/2+1)
165    }
166    env <- list(dim=dm,order=order,theta=theta)
167    ## Create the rk.p function
168    fun <- function(x,y,env,outer.prod=FALSE) {
169        ## Check inputs
170        if (env$dim==1) {
171            if (!(is.vector(x)&is.vector(y))) {
172                stop("gss error in rk: inputs are of wrong types")
173            }
174        }
175        else {
176            if (is.vector(x)) x <- t(as.matrix(x))
177            if (env$dim!=dim(x)[2]) {
178                stop("gss error in rk: inputs are of wrong dimensions")
179            }
180            if (is.vector(y)) y <- t(as.matrix(y))
181            if (env$dim!=dim(y)[2]) {
182                stop("gss error in rk: inputs are of wrong dimensions")
183            }
184        }
185        ## Return the results
186        if (outer.prod) {
187            if (env$dim==1) {
188                fn1 <- function(x,y) abs(x-y)
189                d <- outer(x,y,fn1)
190            }
191            else {
192                fn2 <- function(x,y) sqrt(sum((x-y)^2))
193                d <- NULL
194                for (i in 1:dim(y)[1]) d <- cbind(d,apply(x,1,fn2,y[i,]))
195            }
196        }
197        else {
198            if (env$dim==1) d <- abs(x-y)
199            else {
200                N <- max(dim(x)[1],dim(y)[1])
201                x <- t(matrix(t(x),env$dim,N))
202                y <- t(matrix(t(y),env$dim,N))
203                fn <- function(x) sqrt(sum(x^2))
204                d <- apply(x-y,1,fn)
205            }
206        }
207        power <- 2*env$order-env$dim
208        switch(1+env$dim%%2,
209               env$theta*d^power*log(ifelse(d>0,d,1)),
210               env$theta*d^power)
211    }
212    ## Return the function and the environment
213    list(fun=fun,env=env)
214}
215
216## Make pseudo phi function for thin-plate splines
217mkphi.tp.p <- function(dm,order)
218{
219    ## Check inputs
220    if (!((2*order>dm)&(dm>=1))) {
221        stop("gss error: thin-plate spline undefined for the parameters")
222    }
223    ## Create the environment
224    pol.code <- NULL
225    for (i in 0:(order^dm-1)) {
226        ind <- i; code <- NULL
227        for (j in 1:dm) {
228            code <- c(code,ind%%order)
229            ind <- ind%/%order
230        }
231        if (sum(code)<order) pol.code <- cbind(pol.code,code)
232    }
233    env <- list(dim=dm,pol.code=pol.code)
234    ## Create the phi function
235    fun <- function(x,nu,env) {
236        if (env$dim==1) x <- as.matrix(x)
237        if (env$dim!=dim(x)[2]) {
238            stop("gss error in phi: inputs are of wrong dimensions")
239        }
240        apply(t(x)^env$pol.code[,nu],2,prod)
241    }
242    ## Return the function and the environment
243    list(fun=fun,env=env)
244}
245
246## Make RK for spherical splines
247mkrk.sphere <- function(order)
248{
249    ## Create the environment
250    env <- list(order=order)
251    ## Create the rk function
252    fun <- function(x,y,env,outer.prod=FALSE) {
253        ##% Check the inputs
254        if (is.vector(x)) x <- t(as.matrix(x))
255        if (is.vector(y)) y <- t(as.matrix(y))
256        if (!(is.matrix(x)&is.matrix(y))) {
257            stop("gss error in rk: inputs are of wrong types")
258        }
259        if ((dim(x)[2]!=2)|(dim(y)[2]!=2)) {
260            stop("gss error in rk: inputs are of wrong dimensions")
261        }
262        if ((max(abs(x[,1]),abs(y[,1]))>90)|(max(abs(x[,2]),abs(y[,2]))>180)) {
263            stop("gss error in rk: inputs are out of range")
264        }
265        ##% Convert to radian
266        lat.x <- x[,1]/180*pi; lon.x <- x[,2]/180*pi
267        lat.y <- y[,1]/180*pi; lon.y <- y[,2]/180*pi
268        ##% Return the result
269        rk <- function(lat.x,lon.x,lat.y,lon.y,order) {
270            z <- cos(lat.x)*cos(lat.y)*cos(lon.x-lon.y)+sin(lat.x)*sin(lat.y)
271            W <- ifelse(z<1-10^(-10),(1-z)/2,0)
272            A <- ifelse(W>0,log(1+1/sqrt(W)),0)
273            C <- ifelse(W>0,2*sqrt(W),0)
274            switch(order-1,
275                   (A*4*W*(3*W-1)+6*W*(1-C)+1)/2,
276                   (W*W*(A*((840*W-720)*W+72)+420*W*(1-C)+220*C-150)-4*W+3)/12,
277                   (W*W*W*(A*(((27720*W-37800)*W+12600)*W-600)+
278                           (13860*(1-C)*W+14280*C-11970)*W-2772*C+1470)+
279                    15*W*W-3*W+5)/30) - 1/(2*order-1)
280        }
281        if (outer.prod) {
282            zz <- NULL
283            for (i in 1:length(lat.y))
284                zz <- cbind(zz,rk(lat.x,lon.x,lat.y[i],lon.y[i],env$order))
285        }
286        else zz <- rk(lat.x,lon.x,lat.y,lon.y,env$order)
287        zz
288    }
289    ## Return the function and the environment
290    list(fun=fun,env=env)
291}
292