Nothing
vpaws <- function(y,
kstar = 16,
sigma2 = 1,
invcov = NULL,
mask = NULL,
scorr = 0,
spmin = 0.25,
ladjust = 1,
wghts = NULL,
u = NULL,
patchsize = 1) {
args <- match.call()
dy <- dim(y)
nvec <- dy[1]
n <- prod(dy[-1])
if(is.null(mask)){
dy <- dy[-1]
} else {
dy <- dim(mask)
if(is.null(dy)) dy <- length(mask)
}
d <- length(dy)
if (d > 3)
stop("Vector AWS for more than 3 dimensional grids is not implemented")
#
# set appropriate defaults
#
if(is.null(mask)) mask <- array(TRUE,dy)
nvoxel <- sum(mask)
position <- array(0,dy)
position[mask] <- 1:nvoxel
condensed <- n==nvoxel & n<prod(dy)
if(!condensed & prod(dy)!=n) stop("incompatible mask and data")
if(!condensed & !is.null(u)){
if(length(dim(u))==d){
u <- u[mask]
} else if(length(dim(u))==(d+1)&dim(u)[1]==nvec){
dim(u) <- c(nvec,prod(dy))
u <- u[,mask]
}
}
if(!is.null(invcov)) sigma2 <- 1
# set sigma2 to an uninformative value if invcov is given
lambda <-
2 * sigma2 * ladjust * switch(d,
qchisq(pchisq(14.6, 1), nvec),
## 1D
qchisq(pchisq(9.72, 1), nvec),
## 2D
qchisq(pchisq(8.82, 1), nvec))## 3D
if (is.null(wghts))
wghts <- c(1, 1, 1)
wghts <-
switch(length(dy), c(0, 0), c(wghts[1] / wghts[2], 0), wghts[1] / wghts[2:3])
n1 <- switch(d, dy, dy[1], dy[1])
n2 <- switch(d, 1, dy[2], dy[2])
n3 <- switch(d, 1, 1, dy[3])
n <- n1 * n2 * n3
if(!condensed){
dim(y) <- c(nvec,prod(dy))
y <- y[,mask]
}
if(!is.null(invcov)){
dinvcov <- dim(invcov)
nvd <- dinvcov[1]
if(nvd!=(nvec+1)*nvec/2) stop(paste("First dimension of invcov should be",
nvd,"(dense storage)"))
if(condensed){
if(dinvcov[2]!=nvoxel) stop("invcov should contain information for all voxel in mask")
} else {
if(prod(dinvcov[-1])!=n) stop("invcov should contain information for all voxel")
dim(invcov) <- c(nvd,n)
invcov <- invcov[,mask]
}
}
h0 <- 0
if (any(scorr > 0)) {
h0 <- numeric(length(scorr))
for (i in 1:length(h0))
h0[i] <- geth.gauss(scorr[i])
if (length(h0) < d)
h0 <- rep(h0[1], d)
cat("Corresponding bandwiths for specified correlation:",
h0,
"\n")
}
hseq <- 1
zobj <- list(bi = rep(1, nvoxel), theta = y)
bi <- zobj$bi
cat("Progress:")
total <- cumsum(1.25 ^ (1:kstar)) / sum(1.25 ^ (1:kstar))
mc.cores <- setCores(, reprt = FALSE)
np1 <- 2 * patchsize + 1
np2 <- if (n2 > 1)
2 * patchsize + 1
else
1
np3 <- if (n3 > 1)
2 * patchsize + 1
else
1
k <- 1
hmax <- 1.25 ^ (kstar / d)
lambda0 <- lambda
mae <- NULL
while (k <= kstar) {
hakt0 <- gethani(1, 1.25 * hmax, 2, 1.25 ^ (k - 1), wghts, 1e-4)
hakt <- gethani(1, 1.25 * hmax, 2, 1.25 ^ k, wghts, 1e-4)
cat("step", k, "hakt", hakt, "time", format(Sys.time()), "\n")
hseq <- c(hseq, hakt)
dlw <- (2 * trunc(hakt / c(1, wghts)) + 1)[1:d]
if (scorr[1] >= 0.1)
lambda0 <-
lambda0 * Spatialvar.gauss(hakt0 / 0.42445 / 4, h0, d) /
Spatialvar.gauss(hakt0 / 0.42445 / 4, 1e-5, d)
if(is.null(invcov)){
zobj <- .Fortran(C_pvaws,
as.double(y),
as.integer(position),
as.integer(nvec),
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(lambda0),
as.double(zobj$theta),
as.double(zobj$bi),
bi = double(nvoxel),
theta = double(nvec * nvoxel),
as.integer(mc.cores),
as.double(spmin),
double(prod(dlw)),
as.double(wghts),
double(nvec * mc.cores),
as.integer(np1),
as.integer(np2),
as.integer(np3))[c("bi", "theta", "hakt")]
} else {
zobj <- .Fortran(C_pvaws2,
as.double(y),
as.integer(position),
as.integer(nvec),
as.integer(nvd),
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(lambda0),
as.double(zobj$theta),
as.double(zobj$bi),
bi = double(nvoxel), #binn
theta = double(nvec * nvoxel),
as.double(invcov),
as.integer(mc.cores),
as.double(spmin),
double(prod(dlw)),
as.double(wghts),
double(nvec * mc.cores),
as.integer(np1),
as.integer(np2),
as.integer(np3))[c("bi", "theta", "hakt")]
}
if (!is.null(u)) {
cat(
"bandwidth: ",
signif(hakt, 3),
" MSE: ",
signif(mean((zobj$theta - u) ^ 2), 3),
" MAE: ",
signif(mean(abs(zobj$theta - u)), 3),
" mean(bi)=",
signif(mean(zobj$bi), 3),
"\n"
)
mae <- c(mae, signif(mean(abs(zobj$theta - u)), 3))
}
x <- 1.25 ^ k
scorrfactor <- x / (3 ^ d * prod(scorr) * prod(h0) + x)
lambda0 <- lambda * scorrfactor
if (max(total) > 0) {
cat(signif(total[k], 2) * 100, "% . ", sep = "")
}
k <- k + 1
gc()
}
cat("\n")
# list(y=y,theta=zobj$theta,hakt=hakt,sigma2=sigma2,lambda=lambda,
# ladjust=ladjust,args=args,wghts=wghts,mae=mae,ni=zobj$bi)
sigma2new <- if(is.null(invcov)) zobj$bi/sigma2 else
sweep(invcov, 2, zobj$bi, "*")
# this is the new inverse covariance
if(!condensed){
# expand results
y0 <- y
y <- array(0,c(nvec,prod(dy)))
y[,mask] <- y0
rm(y0)
dim(y) <- c(nvec,dy)
theta <- array(0,c(nvec,prod(dy)))
theta[,mask] <- zobj$theta
dim(theta) <- c(nvec,dy)
if(!is.null(invcov)){
sigma20 <- sigma2new
sigma2new <- array(0,c(nvd,prod(dy)))
sigma2new[,mask] <- sigma20
rm(sigma20)
dim(sigma2new) <- c(nvd,dy)
}
ni <- array(0,dy)
ni[mask] <- zobj$bi
} else {
# set dimensions
dim(y) <- c(nvec,nvoxel)
theta <- array(zobj$theta,c(nvec,nvoxel))
ni <- array(zobj$bi,nvoxel)
}
awsobj(
y,
theta,
sigma2new,
hakt,
sigma2,
lkern = 1L,
lambda,
ladjust,
aws = TRUE,
nvec = as.integer(nvec),
memory = FALSE,
args,
hseq = hseq,
homogen = FALSE,
mask = mask,
earlystop = FALSE,
family = "Gaussian",
wghts = wghts,
mae = mae,
psnr = NULL,
ni = ni
)
}
vpawscov <- function(y,
kstar = 16,
invcov = NULL,
mask = NULL,
scorr = 0,
spmin = 0.25,
ladjust = 1,
wghts = NULL,
maxni = FALSE,
patchsize = 1) {
##
## this is the version with full size invcov (triangular storage)
## needed for MPM
##
args <- match.call()
dy <- dim(y)
nvec <- dy[1]
dy <- dy[-1]
d <- length(dy)
if (length(dy) > 3)
stop("Vector AWS for more than 3 dimensional grids is not implemented")
lambda <- 2 * ladjust * switch(d,
qchisq(pchisq(14.6, 1), nvec),
## 1D
qchisq(pchisq(9.72, 1), nvec),
## 2D
qchisq(pchisq(8.82, 1), nvec))## 3D
lambda <- lambda * switch(patchsize+1,1,1.3,1.6)
if (is.null(wghts))
wghts <- c(1, 1, 1)
wghts <-
switch(length(dy), c(0, 0), c(wghts[1] / wghts[2], 0), wghts[1] / wghts[2:3])
n1 <- switch(d, dy, dy[1], dy[1])
n2 <- switch(d, 1, dy[2], dy[2])
n3 <- switch(d, 1, 1, dy[3])
n <- n1 * n2 * n3
if (is.null(mask))
mask <- rep(TRUE, n)
h0 <- 0
if (any(scorr > 0)) {
h0 <- numeric(length(scorr))
for (i in 1:length(h0))
h0[i] <- geth.gauss(scorr[i])
if (length(h0) < d)
h0 <- rep(h0[1], d)
cat("Corresponding bandwiths for specified correlation:",
h0,
"\n")
}
## create index information for voxel in mask
nvoxel <- sum(mask)
position <- array(0,dy)
position[mask] <- 1:nvoxel
dim(mask) <- NULL
dim(y) <- c(nvec,n)
dim(invcov) <- c(nvec*(nvec+1)/2,n)
hseq <- 1
zobj <- list(bi = rep(1, nvoxel), theta = y[,mask])
bi <- zobj$bi
cat("Progress:")
total <- cumsum(1.25 ^ (1:kstar)) / sum(1.25 ^ (1:kstar))
mc.cores <- setCores(, reprt = FALSE)
np1 <- 2 * patchsize + 1
np2 <- if (n2 > 1) 2 * patchsize + 1 else 1
np3 <- if (n3 > 1) 2 * patchsize + 1 else 1
k <- 1
hmax <- 1.25 ^ (kstar / d)
lambda0 <- lambda
while (k <= kstar) {
hakt0 <- gethani(1, 1.25 * hmax, 2, 1.25 ^ (k - 1), wghts, 1e-4)
hakt <- gethani(1, 1.25 * hmax, 2, 1.25 ^ k, wghts, 1e-4)
cat("step", k, "hakt", hakt, "time", format(Sys.time()), "\n")
hseq <- c(hseq, hakt)
dlw <- (2 * trunc(hakt / c(1, wghts)) + 1)[1:d]
if (scorr[1] >= 0.1)
lambda0 <-
lambda0 * Spatialvar.gauss(hakt0 / 0.42445 / 4, h0, d) / Spatialvar.gauss(hakt0 /
0.42445 / 4, 1e-5, d)
zobj <- .Fortran(C_pvaws2,
as.double(y[,mask]),
as.integer(position),
as.integer(nvec),
as.integer(nvec * (nvec + 1) / 2),
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(lambda0),
as.double(zobj$theta),
as.double(zobj$bi),
bi = double(nvoxel), #binn
theta = double(nvec * nvoxel),
as.double(invcov[,mask]),# compact storage
as.integer(mc.cores),
as.double(spmin),
double(prod(dlw)),
as.double(wghts),
double(nvec * mc.cores),
as.integer(np1),
as.integer(np2),
as.integer(np3))[c("bi", "theta", "hakt")]
if (maxni)
bi <- zobj$bi <- pmax(bi, zobj$bi)
cat(
"bandwidth: ",
signif(hakt, 3),
" mean(bi)=",
signif(mean(zobj$bi), 3),
"\n"
)
x <- 1.25 ^ k
scorrfactor <- x / (3 ^ d * prod(scorr) * prod(h0) + x)
lambda0 <- lambda * scorrfactor
if (max(total) > 0) {
cat(signif(total[k], 2) * 100, "% . ", sep = "")
}
k <- k + 1
gc()
}
cat("\n")
theta <- array(0,c(nvec,n))
theta[,mask] <- zobj$theta
dim(theta) <- c(nvec,dy)
dim(y) <- c(nvec,dy)
dim(invcov) <- c(nvec*(nvec+1)/2,dy)
bi <- array(0,dy)
bi[mask] <- zobj$bi
rm(zobj)
awsobj(
y,
theta,
sweep(invcov, 2:(d + 1), bi, "*"),
hakt,
invcov,
lkern = 1L,
lambda,
ladjust,
aws = TRUE,
memory = FALSE,
args,
hseq = hseq,
homogen = FALSE,
earlystop = FALSE,
family = "Gaussian",
wghts = wghts,
mae = NULL,
psnr = NULL,
ni = bi
)
}
vpawscov2 <- function(y,
kstar = 16,
invcov = NULL,
mask = NULL,
scorr = 0,
spmin = 0.25,
lambda = NULL,
ladjust = 1,
wghts = NULL,
patchsize = 1,
data = NULL,
verbose = TRUE) {#1
##
## this is the version with full size invcov (triangular storage)
## and optional smoothing of vector-valued images supplied in data
## for internal use in package qMRI
## Uses condensed data (voxel within mask only)
## returns a list
##
dy <- dim(y)
nvec <- dy[1]
if(nvec>5) stop("limited to 5 parameters")
indcov <- switch(nvec,1,
c(1,2,4),
c(1,2,5,3,6,9),
c(1,2,6,3,7,11,4,8,12,16),
c(1,2,7,3,8,13,4,9,14,19,5,10,15,20,25))
if(!is.null(data)) nsample <- dim(data)[1]
dy <- dim(mask)
d <- length(dy)
if (d != 3)
stop("need 3D mask")
if(is.null(lambda)){#2
lambda <- 2 * ladjust * qchisq(pchisq(8.82, 1), nvec)
lambda <- lambda * switch(patchsize+1,1,1.3,1.6)
}#2
if (is.null(wghts)) wghts <- c(1, 1, 1)
wghts <- wghts[1] / wghts[2:3]
n1 <- dy[1]
n2 <- dy[2]
n3 <- dy[3]
h0 <- 0
if (any(scorr > 0)) {#3
h0 <- numeric(length(scorr))
for (i in 1:length(h0))
h0[i] <- geth.gauss(scorr[i])
if (length(h0) < d)
h0 <- rep(h0[1], d)
if(verbose) cat("Corresponding bandwiths for specified correlation:",
h0,
"\n")
}#3
## create index information for voxel in mask
nvoxel <- sum(mask)
position <- array(0,dy)
position[mask] <- 1:nvoxel
dim(mask) <- NULL
dim(y) <- c(nvec,nvoxel)
dim(invcov) <- c(nvec * nvec,nvoxel)
hseq <- 1
zobj <- list(bi = rep(1, nvoxel), theta = y)
bi <- zobj$bi
if(verbose) cat("Progress:")
total <- cumsum(1.25 ^ (1:kstar)) / sum(1.25 ^ (1:kstar))
mc.cores <- setCores(, reprt = FALSE)
np1 <- 2 * patchsize + 1
np2 <- if (n2 > 1) 2 * patchsize + 1 else 1
np3 <- if (n3 > 1) 2 * patchsize + 1 else 1
k <- 1
hmax <- 1.25 ^ (kstar / d)
lambda0 <- 1e32
mae <- NULL
while (k <= kstar) {#4
hakt0 <- gethani(1, 1.25 * hmax, 2, 1.25 ^ (k - 1), wghts, 1e-4)
hakt <- gethani(1, 1.25 * hmax, 2, 1.25 ^ k, wghts, 1e-4)
if(verbose) cat("step", k, "hakt", hakt, "time", format(Sys.time()), "\n")
hseq <- c(hseq, hakt)
dlw <- (2 * trunc(hakt / c(1, wghts)) + 1)
if(k==kstar & !is.null(data)){#5
dim(data) <- c(nsample,nvoxel)
zobj <- .Fortran(C_pvawsme,
as.double(y),
as.double(data), ## data to smooth additionally
as.integer(position),
as.integer(nvec),
as.integer(nvec * (nvec + 1) / 2),
as.integer(nsample), ## leading dimension of data
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(lambda0),
as.double(zobj$theta),
as.double(zobj$bi),
bi = double(nvoxel), #binn
theta = double(nvec * nvoxel),
data = double(nsample*nvoxel),
as.double(invcov[indcov,]),#
as.integer(mc.cores),
as.double(spmin),
double(prod(dlw)),
as.double(wghts),
double(nvec * mc.cores),
double(nsample * mc.cores),
as.integer(np1),
as.integer(np2),
as.integer(np3))[c("bi", "theta", "hakt","data")]
dim(zobj$data) <- c(nsample, nvoxel)
} else {#6
zobj <- .Fortran(C_pvaws2,
as.double(y),
as.integer(position),
as.integer(nvec),
as.integer(nvec * (nvec + 1) / 2),
as.integer(n1),
as.integer(n2),
as.integer(n3),
hakt = as.double(hakt),
as.double(lambda0),
as.double(zobj$theta),
as.double(zobj$bi),
bi = double(nvoxel), #binn
theta = double(nvec * nvoxel),
as.double(invcov[indcov,]),# compact storage
as.integer(mc.cores),
as.double(spmin),
double(prod(dlw)),
as.double(wghts),
double(nvec * mc.cores),
as.integer(np1),
as.integer(np2),
as.integer(np3))[c("bi", "theta", "hakt")]
}#6
x <- 1.25 ^ k
scorrfactor <- x / (3 ^ d * prod(scorr) * prod(h0) + x)
lambda0 <- lambda * scorrfactor
if (verbose & max(total) > 0) {#7
cat(signif(total[k], 2) * 100, "% ", sep = "")
cat("mean(bi)", signif(mean(zobj$bi),3)," ")
}#7
k <- k + 1
gc()
}
dim(zobj$theta) <- c(nvec, nvoxel)
if(verbose) cat("\n")
list(
theta=zobj$theta,
hakt=hakt,
lambda=lambda,
hseq = hseq,
bi = zobj$bi,
data= if(!is.null(zobj$data)) zobj$data else NULL
)
}
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.