Nothing
## code for soap film smoothing to deal with difficult boundary regions
## Copyright Simon Wood 2006-2012
unconstrain <- function(object,beta) {
## function to produce full version of constrained coefficients of
## smooth object. Returned vector may have an attribute "constant"
## to be subtraccted from results.
## NOTE: sum to zero on some parameters only branch is not fully
## tested (also unused at present)!
del.index <- attr(object,"del.index")
if (!is.null(del.index)) {
beta.full <- rep(0,length(beta)+length(del.index))
k <- 1;j <- 1
for (i in 1:length(beta.full)) {
if (j <= length(del.index) && i==del.index[j]) {
beta.full[i] <- 0;j <- j + 1
} else {
beta.full[i] <- beta[k];k <- k + 1
}
}
beta <- beta.full
} ## end of del.index handling
qrc <- attr(object,"qrc")
if (!is.null(qrc)) { ## then smoothCon absorbed constraints
j <- attr(object,"nCons")
if (j>0) { ## there were constraints to absorb - need to untransform
k <- length(beta) + j
if (inherits(qrc,"qr")) {
indi <- attr(object,"indi") ## index of constrained parameters
if (is.null(indi)) {
## X <- t(qr.qty(qrc,t(X))[(j+1):k,,drop=FALSE]) ## XZ
beta <- qr.qy(qrc,c(rep(0,j),beta))
} else { ## only some parameters are subject to constraint
## NOTE: this branch not fully tested
##nx <- length(indi)
##nc <- j;nz <- nx - nc
Xbeta <- qr.qy(qrc,c(rep(0,j),beta[indi]))
beta.full <- rep(0,length(beta)+j)
ib <- 1;ii <- 1
for (i in 1:length(beta.full)) {
if (i==indi[ii]) {
beta.full[i] <- Xbeta[ii]; ii <- ii + 1
} else {
beta.full[i] <- beta[ib]; ib <- ib + 1
}
}
##X[,indi[1:nz]]<-t(qr.qty(qrc,t(X[,indi,drop=FALSE]))[(nc+1):nx,,drop=FALSE])
## X <- X[,-indi[(nz+1):nx]]
beta <- beta.full
}
} else if (inherits(qrc,"sweepDrop")) {
## Sweep and drop constraints. First element is index to drop.
## Remainder are constants to be swept out of remaining columns
## X <- sweep(X[,-qrc[1],drop=FALSE],2,qrc[-1])
#X <- X[,-qrc[1],drop=FALSE] - matrix(qrc[-1],nrow(X),ncol(X)-1,byrow=TRUE)
cnst <- sum(beta*qrc[-1])
if (qrc[1]==1) beta <- c(0,beta) else
if (qrc[1]==length(beta)+1) beta <- c(beta,0) else
beta <- c(beta[1:(qrc[1]-1)],0,beta[qrc[1]:length(beta)])
attr(beta,"constant") <- cnst
} else if (qrc>0) { ## simple set to zero constraint
##X <- X[,-qrc]
if (qrc==1) beta <- c(0,beta) else
if (qrc==length(beta)+1) beta <- c(beta,0) else
beta <- c(beta[1:(qrc-1)],0,beta[qrc:length(beta)])
} else if (qrc<0) { ## params sum to zero
# X <- t(diff(t(X)))
beta <- t(diff(diag(length(beta)+1)))%*%beta
}
} ## end if (j>0)
} ## end if qrc exists
beta
} ## end of unconstrain
bnd2C <- function(bnd) {
## converts boundary loop list to form required in C code.
n.loop <- 1
if (is.null(bnd$x)) { ## translate into form that C routine needs
bn <- list(x=bnd[[1]]$x,y=bnd[[1]]$y)
n.loop <- length(bnd)
if (length(bnd)>1) for (i in 2:n.loop) {
bn$x <- c(bn$x,NA,bnd[[i]]$x);bn$y <- c(bn$y,NA,bnd[[i]]$y)
}
bnd <- bn
}
## replace NA segment separators with a numeric code
lowLim <- min(c(bnd$x,bnd$y),na.rm=TRUE)-1
ind <- is.na(bnd$x)|is.na(bnd$y)
bnd$x[ind] <- bnd$y[ind] <- lowLim - 1
bnd$n <- length(bnd$x)
if (bnd$n != length(bnd$y)) stop("x and y must be same length")
bnd$breakCode <-lowLim
bnd$n.loop <- n.loop
bnd
} ## end bnd2C
inSide <- function(bnd,x,y)
## tests whether each point x[i],y[i] is inside the boundary defined
## by bnd$x, bnd$y, or by multiple boundary loops in bnd[[1]]$x,
## bnd[[1]]$y, bnd[[2]]$x, ... etc.
## names in bnd must match those of x and y, but do not need to be "x" and "y"
{ ## match the names up first...
xname <- deparse(substitute(x))
yname <- deparse(substitute(y))
bnd.name <- names(bnd)
if (is.null(bnd.name)) for (i in 1:length(bnd)) {
bnd.name <- names(bnd[[i]])
if (xname%in%bnd.name==FALSE||yname%in%bnd.name==FALSE) stop("variable names don't match boundary names")
bnd.name[xname==bnd.name] <- "x"
bnd.name[yname==bnd.name] <- "y"
names(bnd[[i]]) <- bnd.name
} else {
if (xname%in%bnd.name==FALSE||yname%in%bnd.name==FALSE) stop("variable names don't match boundary names")
bnd.name[xname==bnd.name] <- "x"
bnd.name[yname==bnd.name] <- "y"
names(bnd) <- bnd.name
}
## now do the real stuff...
bnd <- bnd2C(bnd)
um <-.C(C_in_out,bx=as.double(bnd$x),by=as.double(bnd$y),break.code=as.double(bnd$breakCode),
x=as.double(x),y=as.double(y),inside=as.integer(y*0),nb=as.integer(bnd$n),
n=as.integer(length(x)))
as.logical(um$inside)
} ## end inSide
process.boundary <- function(bnd)
## takes a list of boundary loops, makes sure that they join up
## and add a distance along loop array, d to each list element.
{ for (i in 1:length(bnd)) {
x <- bnd[[i]]$x;y<-bnd[[i]]$y;n <- length(x)
if (length(y)!=n) stop("x and y not same length")
if (x[1]!=x[n]||y[1]!=y[n]) { ## then loop not closed, so close it
n<-n+1;x[n] <- x[1];y[n] <- y[1]
if (inherits(bnd[[i]],"data.frame")) bnd[[i]][n,] <-bnd[[i]][1,]
else { ## hopefully a list!
bnd[[i]]$x[n] <- x[1];bnd[[i]]$y[n] <- y[1]
if (!is.null(bnd[[i]]$f)) bnd[[i]]$f[n] <- bnd[[i]]$f[1]
}
}
len <- c(0,sqrt((x[1:(n-1)]-x[2:n])^2+(y[1:(n-1)]-y[2:n])^2)) ## seg lengths
bnd[[i]]$d<-cumsum(len) ## distance along boundary
}
bnd
} ## end process.boundary
crunch.knots <- function(G,knots,x0,y0,dx,dy)
## finds indices of knot locations in solution grid
## the knot x,y locations are given in the `knots' argument.
{ nk <- length(knots$x)
nx <- ncol(G);ny <- nrow(G)
ki <- rep(0,nk)
if (nk==0) return(ki)
for (k in 1:nk) {
i <- round((knots$x[k]-x0)/dx)+1
j <- round((knots$y[k]-y0)/dy)+1
if (i>1&&i<=nx&&j>1&&j<=ny) {
ki[k] <- G[j,i]
if (ki[k] <= 0) {
str <- paste("knot",k,"is on or outside boundary")
stop(str)
}
}
} ## all knots done
ki ## ki[k] indexes kth knot in solution grid
} ## end crunch.knots
setup.soap <- function(bnd,knots,nmax=100,k=10,bndSpec=NULL) {
## setup soap film smooth - nmax is number of grid cells for longest side
## it's important that grid cells are square!
## check boundary...
if (!inherits(bnd,"list")) stop("bnd must be a list.")
n.loops <- length(bnd)
if (n.loops!=length(k)) {
if (length(k)==1) k <- rep(k,n.loops)
else stop("lengths of k and bnd are not compatible.")
}
bnd <- process.boundary(bnd) ## add distances and close any open loops
## create grid on which to solve Laplace equation
## Obtain grid limits from boundary 'bnd'....
x0 <- min(bnd[[1]]$x);x1 <- max(bnd[[1]]$x)
y0 <- min(bnd[[1]]$y);y1 <- max(bnd[[1]]$y)
if (length(bnd)>1) for (i in 2:length(bnd)) {
x0 <- min(c(x0,bnd[[i]]$x)); x1 <- max(c(x1,bnd[[i]]$x))
y0 <- min(c(y0,bnd[[i]]$y)); y1 <- max(c(y1,bnd[[i]]$y))
} ## now got the grid limits, can set it up
if (x1-x0>y1-y0) { ## x is longest side
dy <- dx <- (x1-x0) /(nmax-1)
nx <- nmax
ny <- ceiling((y1-y0)/dy)+1
} else { ## y is longest side
dy <- dx <- (y1-y0) /(nmax-1)
ny <- nmax
nx <- ceiling((x1-x0)/dy)+1
}
## so grid is now nx by ny, cell size is dx by dy (but dx=dy)
## x0, y0 is "lower left" cell centre
## Create grid index G
bnc <- bnd2C(bnd) ## convert boundary to form required in C code
G <- matrix(0,ny,nx)
nb <- rep(0,bnc$n.loop)
oo <- .C(C_boundary,G=as.integer(G), d=as.double(G), dto=as.double(G), x0=as.double(x0),
y0 = as.double(y0), dx=as.double(dx), dy = as.double(dy),
nx=as.integer(nx),as.integer(ny), x=as.double(bnc$x),y=as.double(bnc$y),
breakCode=as.double(bnc$breakCode),n=as.integer(bnc$n),nb=as.integer(nb))
ret <- list(G=matrix(oo$G,ny,nx),nb=oo$nb,d=oo$d[oo$d >= 0],x0=x0,y0=y0,dx=dx,dy=dy,bnd=bnd)
rm(oo)
## Now create the PDE coefficient matrix
n.inside <- sum(ret$G > - nx*ny)
xx <- rep(0,5*n.inside)
o1 <- .C(C_pde_coeffs,as.integer(ret$G),xx=as.double(xx),ii=as.integer(xx),jj=as.integer(xx),
n=as.integer(0),as.integer(nx),as.integer(ny),as.double(dx),as.double(dy))
ind <- 1:o1$n
X <- sparseMatrix(i=o1$ii[ind]+1,j=o1$jj[ind]+1,x=o1$xx[ind])
er <- expand(lu(X))
ret$Q <- er$Q;ret$U <- er$U;ret$L <- er$L;ret$P <- er$P
ret$ng <- n.inside ## the number of cells to solve for
rm(er);rm(X)
## ... so the sparse LU decomposition of X can be used to solve PDE.
## X = PLUQ where P and Q are permuation matrices.
## now obtain location of knots in solution ...
ret$ki <- crunch.knots(ret$G,knots,x0,y0,dx,dy)
## setup the boundary conditions/boundary splines
bc <- list() ## to hold boundary conditions
start <- 1
for (i in 1:length(bnd)) {
stop <- start - 1 + ret$nb[i] ## ret$d[start:stop] = dist along boundary loop i
if (is.null(bnd[[i]]$f)) { ## this boundary is free
d <- c(ret$d[start:stop],0) # boundary gridpoint distances along smooth
if (is.null(bndSpec)) {
bsm <- smooth.construct(s(d,bs="cc",k=k[i]),data=data.frame(d=d),knots=NULL)
} else if (bndSpec$bs=="cc"){
if (bndSpec$knot.space=="even")
knots <- seq(min(d),max(d),length=k[i])
else
knots <- quantile(unique(d),seq(0,1,length=k[i]))
bsm <- smooth.construct(s(d,bs="cc",k=k[i]),data=data.frame(d=d),knots=NULL)
} else { ## use "cp" P-spline
bsm <- smooth.construct(s(d,bs="cp",k=k[i],m=bndSpec$m),data=data.frame(d=d),knots=NULL)
}
bc[[i]] <- list(bsm=bsm,X=bsm$X[1:ret$nb[i],],S=bsm$S[[1]],free.bound=TRUE)
} else { ## boundary is fixed
## pmax/pmin needed to avoid rounding error induced NA's
d <- pmax(pmin(ret$d[start:stop],max(bnd[[i]]$d)),min(bnd[[i]]$d))
ui <- !duplicated(bnd[[i]]$d)
ff <- approx(bnd[[i]]$d[ui],bnd[[i]]$f[ui],d)$y ## fixed values for BC
bc[[i]] <- list(f=ff,free.bound=FALSE)
}
start <- stop + 1
}
ret$bc <- bc
ret
} ## end of setup.soap
soap.basis <- function(sd,x=NA,y=NA,film=TRUE,wiggly=TRUE,penalty=TRUE,plot=FALSE,beta=1) {
## function to evaluate soap basis using soap definition object 'sd'
## returned by setup.soap. x and y are values at which to evaluate.
## If plot==TRUE then then data suitable for plotting are returned at the resolution
## of the solution grid. Then beta contains either the coefficients, or a single number
## representing the single basis function to return (0 for the offset).
if (!plot) {
indout <- inSide(sd$bnd,x,y); n <- length(x)
} else {
penalty <- FALSE
## Some constraints result in the need to add a constant
## to the field (e.g. sweep and drop)
cnst <- attr(beta,"constant")
if (is.null(cnst)) cnst <- 0 else cnst <- -cnst
}
offset.needed <- FALSE;
nc <- length(sd$ki)*as.numeric(wiggly) ## number of interior knots
nb <- 0 ## boundary basis dimension
offset <- NULL
if (film) {
stop <- 0
for (i in 1:length(sd$bc)) { ## work through boundary loops
start <- stop + 1;stop <- start - 1 + sd$nb[i]
if (sd$bc[[i]]$free.bound) nb <- nb + ncol(sd$bc[[i]]$X)
else { ## fixed boundary, so offset required
if (!offset.needed) {
bndOff <- rep(0,sd$ng) ## array for fixed boundary conditions
offset.needed <- TRUE
}
bndOff[start:stop] <- sd$bc[[i]]$f
} ## fixed boundary done
} ## finished first pass through loops
} ## finished first if film
if (plot) { ## preliminaries for plotting info
if (length(beta)==1) { ## just one basis function to be returned
if (beta<0||beta>nc+nb||(beta==0&&!offset.needed)) stop("attempt to select non existent basis function")
select.basis <- TRUE
} else { ## whole smooth to be returned
if (length(beta)!=nc+nb) stop("coefficient vector wrong length")
select.basis <- FALSE
}
G <- sd$G ## solution grid
G[G < - length(G)] <- NA ## exterior marked as NA
ind <- !is.na(G)
gind <- G[ind] <- abs(G[ind])+1
## need to create the indices such that G[gind] <- g is correct...
gind[G[ind]] <- (1:length(G))[ind]
G[ind] <- cnst ## now clear interior of G
} ## finished preliminary if (plot)
if (film) {
if (offset.needed) { ## solve for offset soap film
bndOff <- solve(sd$Q,solve(sd$U,solve(sd$L,solve(t(sd$P),bndOff))))
if (plot) { ## grid is all that's needed
if (select.basis&&beta==0||!select.basis) {
G[gind] <- bndOff
}
} else { ## need full interpolation
NAcode <- max(bndOff)*2
offset <- .C(C_gridder,z=as.double(x),as.double(x),as.double(y),as.integer(length(x)),as.double(bndOff),
as.integer(sd$G),nx=as.integer(ncol(sd$G)),ny=as.integer(nrow(sd$G)),as.double(sd$x0),
as.double(sd$y0),as.double(sd$dx),as.double(sd$dy),as.double(NAcode*2))$z
offset[offset>NAcode] <- NA
offset[!indout] <- NA
}
}
} ## finished preliminary if (film)
if (!plot) {
X <- matrix(0,n,nb+nc) ## model matrix
if (penalty) { S <- list();off <- 1;nS=0} else {off <- S <- NULL}
}
k <- 1 ## model matrix column
if (film&&nb>0) {
## now work through boundary bases
stop <- 0
for (i in 1:length(sd$bc)) { ## work through boundary loops
start <- stop + 1;stop <- start - 1 + sd$nb[i]
ind <- start:stop ## index of this loop in solution grid
if (sd$bc[[i]]$free.bound) {
if (penalty) {
nS <- nS + 1
off[nS] <- k
S[[nS]] <- sd$bc[[i]]$S
} ## penalty done
for (j in 1:ncol(sd$bc[[i]]$X)) { ## loop over loop basis cols
z <- rep(0,sd$ng)
z[ind] <- sd$bc[[i]]$X[,j] ## PDE rhs
z <- solve(sd$Q,solve(sd$U,solve(sd$L,solve(t(sd$P),z))))
if (plot) {
if (select.basis) {
if (beta==k) G[gind] <- z
} else G[gind] <- G[gind] + beta[k]*z
} else {
NAcode <- max(z)*2
Xj <- .C(C_gridder,z=as.double(x),as.double(x),as.double(y),as.integer(length(x)),as.double(z),
as.integer(sd$G),nx=as.integer(ncol(sd$G)),ny=as.integer(nrow(sd$G)),as.double(sd$x0),
as.double(sd$y0),as.double(sd$dx),as.double(sd$dy),as.double(NAcode*2))$z
Xj[Xj>NAcode] <- NA;X[,k] <- Xj;
}
k <- k + 1
} ## basis done
} ## end of free boundary
} ## end of boundary loops
} ## end of film processing
if (wiggly) { ## interior basis functions required
g <- matrix(0,sd$ng,nc)
for (i in 1:nc) g[sd$ki[i],i] <- 1
g <- as(solve(sd$Q,solve(sd$U,solve(sd$L,solve(t(sd$P),g)))),"matrix")
g <- sweep(g,2,apply(g,2,max),"/") ## normalize - not really needed
if (penalty) { ## get soap penalty
nS <- nS + 1;off[nS] <- k
S[[nS]] <- crossprod(g) * sd$dx * sd$dy
}
g <- solve(sd$Q,solve(sd$U,solve(sd$L,solve(t(sd$P),g))))
NAcode <- max(g)*2
for (i in 1:nc) {
if (plot) {
if (select.basis) {
if (k==beta) G[gind] <- g[,i]
} else G[gind] <- G[gind] + beta[k]*g[,i]
} else {
Xj <- .C(C_gridder,z=as.double(x),as.double(x),as.double(y),as.integer(length(x)),as.double(g[,i]),
as.integer(sd$G),nx=as.integer(ncol(sd$G)),ny=as.integer(nrow(sd$G)),as.double(sd$x0),
as.double(sd$y0),as.double(sd$dx),as.double(sd$dy),as.double(NAcode*2))$z
Xj[Xj>NAcode] <- NA;X[,k] <- Xj
}
k <- k + 1
}
}
if (plot) {
return(t(G))
} else {
X[!indout,] <- NA
return(list(X=X,S=S,off=off,offset=offset))
}
} ## end soap.basis
smooth.construct.so.smooth.spec<-function(object,data,knots)
## a full soap film smooth constructor method function for
## integration with mgcv::gam
{ if (is.null(knots)) stop("knots must be specified for soap")
if (object$dim!=2) stop("soap films are bivariate only")
x <- data[[object$term[1]]]
y <- data[[object$term[2]]]
knt <- list(x=knots[[object$term[1]]],y=knots[[object$term[2]]])
if (length(knt$x)<1) stop("need at least one interior knot")
bnd <- object$xt$bnd
if (is.null(bnd)) stop("can't soap smooth without a boundary")
if (!inherits(bnd,"list")) stop("bnd must be a list of boundary loops")
## check knots within boundary...
kin <- in.out(bnd,cbind(knt[[1]],knt[[2]]))
if (any(!kin)) warning("dropping soap knots not inside boundary - use 'in.out' to investigate.")
knt[[1]] <- knt[[1]][kin];knt[[2]] <- knt[[2]][kin]
for (i in 1:length(bnd)) { ## re-lable boundary
nm <- names(bnd[[i]])
ind <- nm==object$term[1]
if (sum(ind)!=1) stop("faulty bnd")
names(bnd[[i]])[ind] <- "x"
ind <- nm==object$term[2]
if (sum(ind)!=1) stop("faulty bnd")
names(bnd[[i]])[ind] <- "y"
}
if (length(object$bs.dim)==1) k <- rep(object$bs.dim,length(bnd))
else {
if (length(object$bs.dim)==length(bnd)) k <- object$bs.dim else
stop("k and bnd lengths are inconsistent")
}
if (is.null(object$xt$nmax)) nmax <- 200 else nmax <- object$xt$nmax
## setup the soap defining structures
sd <- setup.soap(bnd,knots=knt,nmax=nmax,k=k,bndSpec=object$xt$bndSpec)
b <- soap.basis(sd,x,y,film=TRUE,wiggly=TRUE,penalty=TRUE)
if (sum(is.na(b$X))>0) stop("data outside soap boundary")
# b <- soap.construct(x,y,bnd,knots=knt,k=k,n.grid=n.grid,basis.type=2,
# depth=depth,rel.eps=rel.eps,abs.eps=abs.eps,
# bndSpec=object$xt$bndSpec)
## get penalty null space for the term...
ns.dim <- 0;n <- length(sd$bc)
if (n>0) for (i in 1:n) if (sd$bc[[i]]$free.bound)
ns.dim <- ns.dim + sd$bc[[i]]$bsm$null.space.dim
object$null.space.dim <- ns.dim
need.con <- TRUE
for (i in 1:length(sd$bc))
if (!sd$bc[[i]]$free.bound) need.con <- FALSE
## rescale basis for nice conditioning....
irng <- 1/as.numeric(apply(b$X,2,max)-apply(b$X,2,min))
b$X <- t(t(b$X)*irng)
## now apply rescaling
for (i in 1:length(b$S)) {
a <- irng[b$off[i]:(b$off[i]+ncol(b$S[[i]])-1)]
b$S[[i]] <- diag(a)%*%b$S[[i]]%*%diag(a)
}
object$irng <- irng ## the column scaling factor
if (any(!is.finite(irng))) stop("soap basis ill-conditioned - changing 'xt$nmax' may help")
object$X <- b$X ## model matrix
attr(object$X,"offset") <- b$offset
if (!object$fixed) { ## have to unpack a bit...
S <- list();n <- ncol(object$X)
for (i in 1:length(b$S)) {
S[[i]] <- matrix(0,n,n)
m <- ncol(b$S[[i]])
ind <- b$off[i]:(b$off[i]+m-1)
S[[i]][ind,ind] <- b$S[[i]]
}
object$S <- S ## penalties
}
rr <- ncol(b$S[[1]])-1
if (length(b$S)>1) for (i in 2:length(b$S)) rr <- c(rr,ncol(b$S[[i]])-1)
rr[length(rr)] <- rr[length(rr)]+1
object$rank <- rr # penalty ranks
if (!need.con) object$C <- matrix(0,0,ncol(object$X)) ## no con
object$df <- ncol(object$X) # -nrow(object$C)
for (i in 1:length(sd$bc)) {
sd$bc[[i]]$bsm <- sd$bc[[i]]$S <- NULL
}
object$sd <- sd
class(object)<-"soap.film" # Give object a class
object
} ## end of full soap constructor
smooth.construct.sf.smooth.spec<-function(object,data,knots)
## a soap film smooth boundary interpolating film only constructor
## method function for integration with mgcv::gam
{ if (is.null(knots)) stop("knots must be specified for soap")
if (object$dim!=2) stop("soap films are bivariate only")
x <- data[[object$term[1]]]
y <- data[[object$term[2]]]
knt <- list(x=knots[[object$term[1]]],y=knots[[object$term[2]]])
## if (length(knt$x)<1) stop("need at least one interior knot")
bnd <- object$xt$bnd
if (is.null(bnd)) stop("can't soap smooth without a boundary")
if (!inherits(bnd,"list")) stop("bnd must be a list of boundary loops")
for (i in 1:length(bnd)) { ## re-lable boundary
nm <- names(bnd[[i]])
ind <- nm==object$term[1]
if (sum(ind)!=1) stop("faulty bnd")
names(bnd[[i]])[ind] <- "x"
ind <- nm==object$term[2]
if (sum(ind)!=1) stop("faulty bnd")
names(bnd[[i]])[ind] <- "y"
}
if (length(object$bs.dim)==1) k <- rep(object$bs.dim,length(bnd))
else {
if (length(object$bs.dim)==length(bnd)) k <- object$bs.dim else
stop("k and bnd lengths are inconsistent")
}
if (is.null(object$xt$nmax)) nmax <- 200 else nmax <- object$xt$nmax
## setup the soap defining structures
sd <- setup.soap(bnd,knots=knt,nmax=nmax,k=k,bndSpec=object$xt$bndSpec)
b <- soap.basis(sd,x,y,film=TRUE,wiggly=FALSE,penalty=TRUE)
if (sum(is.na(b$X))>0) stop("data outside soap boundary")
if (ncol(b$X)==0) stop("no free coefs in sf smooth")
# b <- soap.construct(x,y,bnd,knots=knt,k=k,n.grid=n.grid,basis.type=2,
# depth=depth,rel.eps=rel.eps,abs.eps=abs.eps,film=TRUE,
# wiggly=FALSE,bndSpec=object$xt$bndSpec)
## get penalty null space for term
ns.dim <- 0;n <- length(sd$bc)
k <- 0 ## counter for b$S
rr <- rep(0,length(b$S))
if (n>0) for (i in 1:n) if (sd$bc[[i]]$free.bound) {
nsd <- sd$bc[[i]]$bsm$null.space.dim
ns.dim <- ns.dim + nsd
k <- k + 1
rr[k] <- ncol(b$S[[k]]) - nsd ## rank of b$S[[k]]
}
object$null.space.dim <- ns.dim
object$rank <- rr # penalty ranks
need.con <- TRUE
for (i in 1:length(sd$bc))
if (!sd$bc[[i]]$free.bound) need.con <- FALSE
## rescale basis for nice conditioning....
irng <- 1/as.numeric(apply(b$X,2,max)-apply(b$X,2,min))
b$X <- t(t(b$X)*irng)
## now apply rescaling
if (length(b$S)>0) for (i in 1:length(b$S)) {
a <- irng[b$off[i]:(b$off[i]+ncol(b$S[[i]])-1)]
b$S[[i]] <- diag(a)%*%b$S[[i]]%*%diag(a)
}
object$irng <- irng ## the column scaling factor
object$X <- b$X ## model matrix
attr(object$X,"offset") <- b$offset
if (!object$fixed) { ## have to unpack a bit...
S <- list();n <- ncol(object$X)
if (length(b$S)>0) for (i in 1:length(b$S)) {
S[[i]] <- matrix(0,n,n)
m <- ncol(b$S[[i]])
ind <- b$off[i]:(b$off[i]+m-1)
S[[i]][ind,ind] <- b$S[[i]]
}
object$S <- S ## penalties
}
if (!need.con) object$C <- matrix(0,0,ncol(object$X)) ## no con
object$df <- ncol(object$X) # -nrow(object$C)
for (i in 1:length(sd$bc)) {
sd$bc[[i]]$bsm <- sd$bc[[i]]$S <- NULL
}
object$sd <- sd
class(object)<-c("sf","soap.film") # Give object a class
object
} ## end of boundary film component soap constructor
smooth.construct.sw.smooth.spec<-function(object,data,knots)
## a soap film smooth wiggly component only constructor method function for
## integration with mgcv::gam
{ if (is.null(knots)) stop("knots must be specified for soap")
if (object$dim!=2) stop("soap films are bivariate only")
x <- data[[object$term[1]]]
y <- data[[object$term[2]]]
knt <- list(x=knots[[object$term[1]]],y=knots[[object$term[2]]])
if (length(knt$x)<1) stop("need at least one interior knot")
bnd <- object$xt$bnd
if (is.null(bnd)) stop("can't soap smooth without a boundary")
if (!inherits(bnd,"list")) stop("bnd must be a list of boundary loops")
for (i in 1:length(bnd)) { ## re-lable boundary
nm <- names(bnd[[i]])
ind <- nm==object$term[1]
if (sum(ind)!=1) stop("faulty bnd")
names(bnd[[i]])[ind] <- "x"
ind <- nm==object$term[2]
if (sum(ind)!=1) stop("faulty bnd")
names(bnd[[i]])[ind] <- "y"
}
if (length(object$bs.dim)==1) k <- rep(object$bs.dim,length(bnd))
else {
if (length(object$bs.dim)==length(bnd)) k <- object$bs.dim else
stop("k and bnd lengths are inconsistent")
}
if (is.null(object$xt$nmax)) nmax <- 200 else nmax <- object$xt$nmax
## setup the soap defining structures
sd <- setup.soap(bnd,knots=knt,nmax=nmax,k=k,bndSpec=object$xt$bndSpec)
b <- soap.basis(sd,x,y,film=FALSE,wiggly=TRUE,penalty=TRUE)
if (sum(is.na(b$X))>0) stop("data outside soap boundary")
object$null.space.dim <- 0 ## penalty is full rank, for this case
## rescale basis for nice conditioning....
irng <- 1/as.numeric(apply(b$X,2,max)-apply(b$X,2,min))
b$X <- t(t(b$X)*irng)
## now apply rescaling
for (i in 1:length(b$S)) {
a <- irng[b$off[i]:(b$off[i]+ncol(b$S[[i]])-1)]
b$S[[i]] <- diag(a)%*%b$S[[i]]%*%diag(a)
}
object$irng <- irng ## the column scaling factor
object$X <- b$X ## model matrix
if (!object$fixed) { ## have to unpack a bit...
S <- list();n <- ncol(object$X)
for (i in 1:length(b$S)) {
S[[i]] <- matrix(0,n,n)
m <- ncol(b$S[[i]])
ind <- b$off[i]:(b$off[i]+m-1)
S[[i]][ind,ind] <- b$S[[i]]
}
object$S <- S ## penalties
}
rr <- ncol(b$S[[1]])-1
if (length(b$S)>1) for (i in 2:length(b$S)) rr <- c(rr,ncol(b$S[[i]])-1)
rr[length(rr)] <- rr[length(rr)]+1
object$rank <- rr # penalty ranks
object$df <- ncol(object$X) # -nrow(object$C)
for (i in 1:length(sd$bc)) {
sd$bc[[i]]$bsm <- sd$bc[[i]]$S <- NULL
}
object$sd <- sd
object$C <- matrix(0,0,ncol(object$X)) ## this is tied to zero
class(object)<-c("sw","soap.film") # Give object a class
object
} ## end of wiggly component of soap constructor
Predict.matrix.soap.film<-function(object,data)
# prediction method function for the soap.film smooth class
{ x <- get.var(object$term[1],data)
y <- get.var(object$term[2],data)
b <- soap.basis(object$sd,x,y,film=TRUE,wiggly=TRUE,penalty=FALSE)
X <- t(object$irng*t(b$X))
attr(X,"offset") <- b$offset
X
}
Predict.matrix.sf <- function(object,data)
# prediction method function for the sf smooth class --- the boundary interpolating film
# component of a soap film smooth
{ x <- get.var(object$term[1],data)
y <- get.var(object$term[2],data)
b <- soap.basis(object$sd,x,y,film=TRUE,wiggly=FALSE,penalty=FALSE)
X <- t(object$irng*t(b$X))
attr(X,"offset") <- b$offset
X
}
Predict.matrix.sw <- function(object,data)
# prediction method function for the sw smooth class --- the wiggly
# component of a soap film smooth
{ x <- get.var(object$term[1],data)
y <- get.var(object$term[2],data)
X <- soap.basis(object$sd,x,y,film=FALSE,wiggly=TRUE,penalty=FALSE)$X
X <- t(object$irng*t(X))
X
}
plot.soap.film <- function(x,P=NULL,data=NULL,label="",se1.mult=1,se2.mult=2,
partial.resids=FALSE,rug=TRUE,se=TRUE,scale=-1,n=100,n2=40,n3=3,
theta=30,phi=30,jit=FALSE,xlab=NULL,ylab=NULL,main=NULL,
ylim=NULL,xlim=NULL,too.far=0.1,shade=FALSE,shade.col="gray80",
shift=0,trans=I,by.resids=FALSE,scheme=0,hcolors=heat.colors(100),
contour.col=1,...) {
## plot method function for soap.smooth terms
if (scheme==3) {
if (is.null(P)) outline <- FALSE else outline <- TRUE
if (is.null(xlim)) xlim <- c(x$sd$x0,x$sd$x0+ncol(x$sd$G)*x$sd$dx)
if (is.null(ylim)) ylim <- c(x$sd$y0,x$sd$y0+nrow(x$sd$G)*x$sd$dy)
P0 <- plot.mgcv.smooth(x=x,P=P,data=data,label=label,se1.mult=se1.mult,se2.mult=se2.mult,
partial.resids=partial.resids,rug=rug,se=se,scale=scale,n=n,n2=n2,
theta=theta,phi=phi,jit=jit,xlab=xlab,ylab=ylab,main=main,
ylim=ylim,xlim=xlim,too.far=too.far,shade=shade,shade.col=shade.col,
shift=shift,trans=trans,by.resids=by.resids,hcolors=hcolors, ## don't pass scheme!!
contour.col=contour.col,...)
if (outline) { if (is.null(names(P$bnd))) {
for (i in 1:length(P$bnd)) lines(P$bnd[[i]],lwd=2)
} else lines(P$bnd,lwd=2)
} else { P0$bnd <- x$sd$bnd}
return(P0)
}
if (is.null(P)) { ## get plotting information...
if (!x$plot.me) return(NULL) ## shouldn't or can't plot
## get basic plot data
beta <- unconstrain(x,attr(x,"coefficients"))*x$irng ## coefs
raw <- data[x$term]
film <- wiggly <- TRUE
if (inherits(x,"sw")) film <- FALSE else if (inherits(x,"sf")) wiggly <- FALSE
soap.basis(x$sd,film=film,wiggly=wiggly,plot=TRUE,beta=beta) -> G
if (is.null(xlab)) xlabel<- x$term[1] else xlabel <- xlab
if (is.null(ylab)) ylabel <- x$term[2] else ylabel <- ylab
xscale <- x$sd$x0 + 0:(nrow(G)-1) * x$sd$dx
yscale <- x$sd$y0 + 0:(ncol(G)-1) * x$sd$dy
main <- if (is.null(main)) label
return(list(fit=G,scale=FALSE,se=FALSE,raw=raw,xlab=xlabel,ylab=ylabel,
xscale=xscale,yscale=yscale,main=main,bnd=x$sd$bnd))
} else { ## do plot
if (scheme==0) {
xlim <- range(P$xscale);dx = xlim[2] - xlim[1]
ylim <- range(P$yscale);dy = ylim[2] - ylim[1]
plot(P$xscale[1],P$yscale[1],xlab=P$xlab,ylab=P$ylab,main=P$main,xlim=xlim,ylim=ylim,...)
rect(xlim[1]-dx,ylim[1]-dy,xlim[2]+dx,ylim[2]+dy,col="lightgrey")
image(P$xscale,P$yscale,P$fit,add=TRUE,col=hcolors,...)
contour(P$xscale,P$yscale,P$fit,add=TRUE,col=contour.col,...)
} else if (scheme==1) {
image(P$xscale,P$yscale,P$fit,col=grey(0:50/50),xlab=P$xlab,
ylab=P$ylab,main=P$main,...)
contour(P$xscale,P$yscale,P$fit,add=TRUE,col=contour.col,...)
} else if (scheme==2) {
contour(P$xscale,P$yscale,P$fit,xlab=P$xlab,
ylab=P$ylab,main=P$main,col=contour.col,...)
if (is.null(names(P$bnd))) {
for (i in 1:length(P$bnd)) lines(P$bnd[[i]],lwd=2)
} else lines(P$bnd,lwd=2)
}
}
} ## end plot.soap.smooth
fs.test <- function(x,y,r0=.1,r=.5,l=3,b=1,exclude=TRUE)
## test function based on Tim Ramsay (2002) J.R.Statist. Soc. B
## 64(2):307-319 "Spline smoothing over difficult regions"
{ q <- pi*r/2 ## 1/2 length of semi-circle part of centre curve
a <- d <- x*0 ## along and distance to arrays
## convert x,y to along curve and distance to curve (a,d)
## co-ordinates. 0 distance along is at (x=-r,y=0)
ind <- x>=0 & y>0
a[ind] <- q + x[ind]
d[ind] <- y[ind]-r
ind <- x>=0 & y<=0
a[ind] <- -q - x[ind]
d[ind] <- -r - y[ind]
ind <- x < 0
a[ind] <- -atan(y[ind]/x[ind])*r
d[ind] <- sqrt(x[ind]^2+y[ind]^2) - r
## create exclusion index
ind <- abs(d)>r-r0 | (x>l & (x-l)^2+d^2 > (r-r0)^2)
# f <- a*b # the original
f <- a*b+d^2
if (exclude) f[ind] <- NA
attr(f,"exclude") <- ind
f
}
fs.boundary <- function(r0=.1,r=.5,l=3,n.theta=20)
## produce boundary file for fs.test
{ rr <- r+(r-r0)
theta <- seq(pi,pi/2,length=n.theta)
x <- rr*cos(theta); y <- rr*sin(theta)
theta <- seq(pi/2,-pi/2,length=2*n.theta)
x <- c(x,(r-r0)*cos(theta)+l); y <- c(y,(r-r0)*sin(theta)+r)
theta <- seq(pi/2,pi,length=n.theta)
x <- c(x,r0*cos(theta)); y <- c(y,r0*sin(theta))
n<-length(x)
x <- c(x,x[n:1]);y <- c(y,-y[n:1])
return(list(x=x,y=y))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.