R/smse3misc.R
In WIAS-BERLIN/aws: Adaptive Weights Smoothing
Defines functions
suggestkappa
getkappas
matrm
betagamma
reduceparam
lkfullse3
lkfulls0
interpolatesphere1
getkappasmsh3
getkappas3
getkappas
getkappasmsh
lkfullse3msh
interpolatesphere
getnext3g
unifybvals
getsphwghts
sphtrarea1
sphtrarea1 <- function(g1,g2,g3){
## Compute area of sherical triangle spanned by vectors
## g1,g2,g3 on unit sphere
## use absolute values to identify opposite directions with each other
c12 <- abs(g1%*%g2)
c13 <- abs(g1%*%g3)
c23 <- abs(g2%*%g3)
s12 <- sqrt(1-c12^2)
s13 <- sqrt(1-c13^2)
s23 <- sqrt(1-c23^2)
a12 <- acos(c12)
a23 <- acos(c23)
a13 <- acos(c13)
## observed b1, b2, b3 outside [-1,1] due to numerics
if(s12*s13*s23>1e-8){
b1 <- min(max(1+(c23-cos(a12-a13))/s12/s13,-1),1)
b2 <- min(max(1+(c13-cos(a12-a23))/s12/s23,-1),1)
b3 <- min(max(1+(c12-cos(a13-a23))/s13/s23,-1),1)
acos(b1)+acos(b2)+acos(b3)-pi
} else {
0
}
}
getsphwghts <- function(g,g1,g2,g3){
##
## compute weights for linear interpolation in g using g1,g2,g3
##
ierr <- 0
w <- numeric(3)
w0 <- sphtrarea1(g1,g2,g3)
w[1] <- sphtrarea1(g,g2,g3)
w[2] <- sphtrarea1(g,g1,g3)
w[3] <- sphtrarea1(g,g1,g2)
if(w0 < sum(w) - 1e-6){
#cat("gradient does not belong to triangle")
#print(g)
#print(cbind(g1,g2,g3))
#cat("w0",w0,"w",w,"alt",w[1]+w[2]-w[3],"\n")
ierr <- 1
}
list(w=w/sum(w),ierr=ierr)
}
unifybvals <- function(bval,dbv=51){
nbv <- length(bval)
nbval <- bval
obval <- numeric(nbv)
while(any(nbval!=obval)){
obval <- nbval
sbv <- sort(obval)
dsbv <- (1:(nbv-1))[diff(sbv)<dbv]
sbv[dsbv+1] <- sbv[dsbv]
obv <- order(obval)
nbval[obv] <- sbv
}
for(bv in unique(nbval)) nbval[nbval==bv] <- trunc(mean(bval[nbval==bv]))
nbval
}
getnext3g <- function(grad,bv){
##
## calculate next neighbors on the sphere and weights for interpolation
##
binomial <- function(a,b) prod(1:a)/prod(1:(a-b))/prod(1:b)
##
grad <- grad[,bv>0]
ng <- dim(grad)[2]
ng2 <- max(4,ng/6)
perm <- matrix(1,(ng2-2)*(ng2-1)*ng2/6,3)
l <- 1
for(k in 1:(ng2-2)){
for(j in (k+1):(ng2-1)){
for(i in (j+1):ng2){
perm[l,] <- c(k,j,i)
l <- l+1
}
}
}
dperm <- l-1
perm <- perm[2:dperm,]
# thats all ordered triples from (1:ng2), without (1,2,3)
bv <- unifybvals(bv[bv>0]) # identifies b-values that differ by less than dbv=51
ubv <- sort(unique(bv[bv>max(bv)/50]))
nbv <- length(ubv)
ind <- array(0,c(3,nbv,ng))
w <- array(0,c(3,nbv,ng))
for(i in 1:nbv){
indb <- (1:ng)[bv==ubv[i]]
ind[1,i,indb] <- indb
ind[2,i,indb] <- indb
ind[3,i,indb] <- indb
w[1,i,indb] <- 1
w[2,i,indb] <- 0
w[3,i,indb] <- 0
for(j in (1:nbv)[-i]){
indbk <- (1:ng)[bv==ubv[j]]
perm0 <- perm
if(length(indbk)<ng2){
indp <- apply(perm < length(indbk),1,all)
perm0 <- perm0[indp,]
}
dperm <- dim(perm0)[1]+1
for(k in indb){
d <- abs(t(grad[,k])%*%grad[,indbk])
od <- order(d,decreasing = TRUE)
ijk0 <- (1:length(indbk))[od]
ijk <- indbk[od]
# cat("k",k,"d",d[od][1:8],"ijk",ijk[1:8],"\n")
ind[,j,k] <- ijk[1:3]
if(max(d)>1-1e-6){
w[,j,k] <- c(1,0,0)
} else {
z <- getsphwghts(grad[,k],grad[,ijk[1]],grad[,ijk[2]],grad[,ijk[3]])
l <- 1
if(z$ierr==1){
# order triplets in perm according
dd <- numeric(dperm-1)
d <- acos(pmin(1,d))
for(l in 1:(dperm-1)){
ijk1 <- ijk0[perm0[l,]]
dd[l] <- d[ijk1[1]]+d[ijk1[2]]+d[ijk1[3]]
}
odd <- order(dd)
l <- 1
}
while(z$ierr==1&&l<dperm){
ijk1 <- ijk[perm0[odd[l],]]
z <- getsphwghts(grad[,k],grad[,ijk1[1]],grad[,ijk1[2]],grad[,ijk1[3]])
ind[,j,k] <- ijk1[1:3]
l<-l+1
}
## if we were running out of options use ijk[1:3]
if(z$ierr==1){
ind[,j,k] <- ijk[1:3]
z <- getsphwghts(grad[,k],grad[,ijk[1]],grad[,ijk[2]],grad[,ijk[3]])
}
w[,j,k] <- z$w
}
}
}
}
# spheres identified by bvalues in ubv
# ind[j,k,] contains indices of gradients to be used in spherical interpolation
# on shell i for gradient k
# w[j,k,] contains the corresponding weights
list(ind=ind, w=w, ubv = ubv, nbv = nbv, bv=bv)
}
interpolatesphere <- function(theta,n3g){
## interpolate estimated thetas to get values on all spheres
## n3g generated by function getnext3g
## dim(theta) = c(n1,n2,n3,ngrad)
dtheta <- dim(theta)
mstheta <- array(0,c(n3g$nbv,dtheta))
for(i in 1:n3g$nbv){
for(j in 1:dim(theta)[2]){
mstheta[i,,j] <- theta[,n3g$ind[,i,j]]%*%n3g$w[,i,j]
}
}
mstheta
}
lkfullse3msh <- function(h,kappa,gradstats,vext,n){
nbv <- gradstats$nbv
dist <- gradstats$dist
ngrad <- gradstats$ngrad
bvind <- gradstats$bvind
ind <- matrix(0,5,n)
w <- numeric(n)
nn <- 0
dist <- 4
for(i in 1:nbv){
gshell <- list(k456=gradstats$k456[[i]],bghat=gradstats$bghat[[i]],dist=dist)
z <- lkfullse3(h[bvind[[i]]],kappa[bvind[[i]]],gshell,vext,n)
ind[1:3,nn+1:z$nind] <- z$ind[1:3,]
ind[4:5,nn+1:z$nind] <- matrix(bvind[[i]][z$ind[4:5,]],2, z$nind)
#
# convert indeces on selected shell to total indeces
#
w[nn+1:z$nind] <- z$w
nn <- nn + z$nind
}
list(h=h,kappa=kappa,ind=ind[,1:nn],w=w[1:nn],nind=nn)
}
gethseqfullse3msh <- function (kstar, gradstats, kappa, vext = c(1, 1))
{
#
# generate information on local bandwidths and variance reduction
# for smoothing on multiple shells
#
nbv <- gradstats$nbv
ngrad <- gradstats$ngrad
h <- vr <- matrix(1,ngrad,kstar+1)
n <- 0
dist <- 4
for(i in 1:nbv){
gshell <- list(k456=gradstats$k456[[i]],bghat=gradstats$bghat[[i]],dist=dist)
z <- gethseqfullse3(kstar, gshell, kappa=kappa, vext=vext)
h[gradstats$bvind[[i]],-1] <- z$h
vr[gradstats$bvind[[i]],-1] <- z$vred
vr[gradstats$bvind[[i]],1] <- vr[gradstats$bvind[[i]],2]/1.25
n <- n+z$n
}
cat("\n total number of positive weights:",n,"mean maximal bandwidth",signif(mean(h[,kstar]),3), "\n")
list(h=h,vred=vr,n=n)
}
gethseqfullse3 <- function (kstar, gradstats, kappa=NULL, vext = c(1, 1))
{
ngrad <- dim(gradstats$bghat)[2]
h <- vr <- matrix(0,ngrad,kstar)
if(is.null(kappa)) kappa <- pi/ngrad
# if(length(kappa)<kstar) kappa <- rep(kappa[1],kstar)
prt0 <- Sys.time()
cat("get sequences of bw, kappa up to kstar=", kstar, " ")
n <- 0
for(i in 1:ngrad){
z <- .Fortran(C_ghfse3i,
as.integer(i),#i4
as.integer(kstar),#kstar
as.double(gradstats$k456),
as.integer(ngrad),
as.double(kappa),#kappa
as.double(vext),#vext
h=double(kstar),
vr=double(kstar),#
n=integer(1),#
as.integer(gradstats$dist))[c("h","vr","n")]
h[i,] <- z$h
vr[i,] <- z$vr
n <- n+z$n
cat(".")
}
cat("\n number of positive weights:",n,"mean maximal bandwidth",signif(mean(h[,kstar]),3),"time elapsed:",
format(difftime(Sys.time(), prt0), digits = 3),
"\n")
list(h=h,kappa=kappa,vred=vr,n=n)
}
getkappasmsh <- function(grad, msstructure, trace = 0, dist = 1){
ngrad <- dim(grad)[2]
nbv <- msstructure$nbv
bv <- msstructure$bv
ubv <- msstructure$ubv
bvind <- k456 <- bghat <- list(NULL)
for(i in 1:nbv){
#
# collect information for spherical distances on each schell separately
#
ind <- (1:ngrad)[bv==ubv[i]]
z <- getkappas(grad[,ind], trace = trace, dist = dist)
bvind[[i]] <- ind
k456[[i]] <- z$k456
bghat[[i]] <- z$bghat
}
list(k456 = k456, bghat = bghat, bvind = bvind, dist=dist, nbv = nbv, ngrad = ngrad)
}
getkappas <- function(grad, trace = 0, dist = 1){
#
# dist = 1: k4^2+k5^2+|k6|
# dist = 2: k4^2+k5^2+k6^2
# dist = 3: acos(g_i%*%g_j)
# new version with analytic expm(par[1]*mx)
#
krit <- function(par, matm, beta){
## sum((matm-expm(par[1]*m4)%*%expm(par[2]*m5)%*%expm(par[3]*m6))^2)
.Fortran(C_k456krb,
as.double(par),
as.double(beta),
as.double(matm),
erg = double(1))$erg
}
krit5 <- function(x,p,pk4,matm,beta){
# for line search with respect to k5 to get second solution
p1 <- p+c(pk4/2,x,pi)
krit(p1,matm,beta)
}
ngrad <- dim(grad)[2]
if(dist<3){
prta <- Sys.time()
cat("Start computing spherical distances", format(Sys.time()), "\n")
kappa456 <- kappa456a <- array(0, c(3, ngrad, ngrad))
bghat <- betagamma(grad)
for (i in 1:ngrad) for (j in 1:ngrad) {
bg <- bghat[, i, j]
# fix for discontinuity
if(abs(cos(bg[1])) < 1.e-6) bg[1] = pi/2 - 1e-6*sign(cos(bg[1]))
matm <- matrm(bg[1], bg[2])
cbg1 <- cos(bg[1])
k456 <- runif(3, -.1, .1)
maxit <- 10000
fnscale <- 1
z <- optim(k456, krit, method = "BFGS", matm = matm, beta=bg[1],
control = list(trace = trace, fnscale=fnscale, maxit=maxit, reltol = 1e-12, abstol = 1e-16))
count <- 5
while (z$value > 1e-14&count>0) {
## cat("i",i,"j",j,"value",z$value,"par",z$par,"\n")
maxit <- maxit*1.5
fnscale <- fnscale/2
k456 <- runif(3, -1, 1)
z <- optim(k456, krit, method = "BFGS", matm = matm, beta=bg[1],
control = list(trace = trace, fnscale=fnscale, maxit=maxit, reltol = 1e-12, abstol = 1e-16))
## cat(" new value",z$value,"par",z$par,"\n")
count <- count - 1
if(count==0)
cat("failed for bg:",bg,"value",z$value,"par",z$par,"\n")
}
z$par[1] <- z$par[1]/cbg1
kappa456[, i, j] <- z$par
pk4 <- abs(2*pi*cbg1/(2-cbg1^2)^.5)
if(kappa456[1,i,j] < -pk4/2) kappa456[1,i,j] <- kappa456[1,i,j] - trunc(kappa456[1, i, j]/pk4-1)*pk4
if(kappa456[1,i,j] >= pk4/2) kappa456[1,i,j] <- kappa456[1,i,j] - trunc(kappa456[1, i, j]/pk4+1)*pk4
if(kappa456[2,i,j] < -pi) kappa456[2,i,j] <- kappa456[2,i,j] - trunc(kappa456[2, i, j]/2/pi-1)*2*pi
if(kappa456[2,i,j] > pi) kappa456[2,i,j] <- kappa456[2,i,j] - trunc(kappa456[2, i, j]/2/pi+1)*2*pi
if(kappa456[3,i,j] < -pi) kappa456[3,i,j] <- kappa456[3,i,j] - trunc(kappa456[3, i, j]/2/pi-1)*2*pi
if(kappa456[3,i,j] > pi) kappa456[3,i,j] <- kappa456[3,i,j] - trunc(kappa456[3, i, j]/2/pi+1)*2*pi
kappa456a[,i,j] <- kappa456[,i,j]+c(pk4/2,0,pi)
kpar <- kappa456[,i,j]*c(cbg1,1,1)
kappa456a[2,i,j] <- optimize(krit5,c(-pi,pi),p=kpar,pk4=pk4,
matm=matm,beta=bg[1],maximum=FALSE)$minimum
if(kappa456a[1,i,j] < -pk4/2) kappa456a[1,i,j] <- kappa456a[1,i,j] - trunc(kappa456a[1, i, j]/pk4-1)*pk4
if(kappa456a[1,i,j] >= pk4/2) kappa456a[1,i,j] <- kappa456a[1,i,j] - trunc(kappa456a[1, i, j]/pk4+1)*pk4
if(kappa456a[2,i,j] < -pi) kappa456a[2,i,j] <- kappa456a[2,i,j] - trunc(kappa456a[2, i, j]/2/pi-1)*2*pi
if(kappa456a[2,i,j] > pi) kappa456a[2,i,j] <- kappa456a[2,i,j] - trunc(kappa456a[2, i, j]/2/pi+1)*2*pi
if(kappa456a[3,i,j] < -pi) kappa456a[3,i,j] <- kappa456a[3,i,j] - trunc(kappa456a[3, i, j]/2/pi-1)*2*pi
if(kappa456a[3,i,j] > pi) kappa456a[3,i,j] <- kappa456a[3,i,j] - trunc(kappa456a[3, i, j]/2/pi+1)*2*pi
}
dka <- switch(dist,kappa456[1,,]^2+kappa456[2,,]^2+abs(kappa456[3,,]),
kappa456[1,,]^2+kappa456[2,,]^2+kappa456[3,,]^2)
dkb <- switch(dist,kappa456a[1,,]^2+kappa456a[2,,]^2+abs(kappa456a[3,,]),
kappa456a[1,,]^2+kappa456a[2,,]^2+kappa456a[3,,]^2)
dim(kappa456) <- dim(kappa456a) <- c(3,ngrad*ngrad)
kappa456[,dkb<dka] <- kappa456a[,dkb<dka]
dim(kappa456) <- c(3,ngrad,ngrad)
prtb <- Sys.time()
cat("End computing spherical distances", format(Sys.time()), "\n")
} else {
kappa456 <- array(0, c(3, ngrad, ngrad))
bghat <- betagamma(grad)
for(i in 1:ngrad) kappa456[1,i,] <- acos(pmin(1,abs(grad[,i]%*%grad)))
}
list(k456 = kappa456, bghat = bghat, dist=dist)
}
getkappas3 <- function(grad, trace = 0){
#
# dist: acos(g_i%*%g_j)
#
ngrad <- dim(grad)[2]
kappa456 <- array(0, c(3, ngrad, ngrad))
bghat <- betagamma(grad)
for(i in 1:ngrad) kappa456[1,i,] <- acos(pmin(1,abs(grad[,i]%*%grad)))
list(k456 = kappa456, bghat = bghat, dist=3)
}
getkappasmsh3 <- function(grad, msstructure, trace = 0){
ngrad <- dim(grad)[2]
nbv <- msstructure$nbv
bv <- msstructure$bv
ubv <- msstructure$ubv
bvind <- k456 <- bghat <- list(NULL)
for(i in 1:nbv){
#
# collect information for spherical distances on each schell separately
#
ind <- (1:ngrad)[bv==ubv[i]]
z <- getkappas3(grad[,ind], trace = trace)
bvind[[i]] <- ind
k456[[i]] <- z$k456
bghat[[i]] <- z$bghat
}
list(k456 = k456, bghat = bghat, bvind = bvind, dist=3, nbv = nbv, ngrad = ngrad)
}
interpolatesphere1 <- function(theta,th0,ni,ni0,n3g){
## interpolate estimated thetas to get values on all spheres
## n3g generated by function getnext3g
## dim(theta) = c(n1,n2,n3,ngrad)
nbv <- n3g$nbv
bv <- n3g$bv
ubv <- n3g$ubv
dtheta <- dim(theta)
dth0 <- dim(th0)
ng <- dtheta[2]
n <- dtheta[1]
#t1 <- Sys.time()
z <- .Fortran(C_ipolsp1,
as.double(theta),
as.double(th0),
as.double(ni),
as.double(ni0),
as.integer(n),
as.integer(ng),
as.integer(n3g$ind),
as.double(n3g$w),
as.integer(nbv),
as.integer(nbv+1),
msth=double((nbv+1)*n*ng),
msni=double((nbv+1)*n*ng))[c("msth","msni")]
#cat("time for sb-interpolation", format(difftime(Sys.time(),t1),digits=3),"\n")
# now fill vector for s0
dim(theta) <- dim(ni) <- c(n,ng)
msth0 <- msni0 <- array(0,c(nbv+1,n))
msth0[1,] <- th0
msni0[1,] <- ni0
for(i in 1:nbv){
indi <- (1:ng)[bv==ubv[i]]
lindi <- length(indi)
# msth0[i+1,mask] <- theta[mask,indi]%*%rep(1/lindi,lindi)
msth0[i+1,] <- .Fortran(C_getmsth0,
as.double(theta[,indi]),
as.integer(n),
as.integer(lindi),
msth0=double(n))$msth0
# correct value would be
# msni0[i+1,] <- 1/(1/(ni[,indi])%*%rep(1/lindi,lindi)^2)
# try to be less conservative by ignorin squares in w
# msni0[i+1,mask] <- 1/((1/ni[mask,indi])%*%rep(1/lindi,lindi))
msni0[i+1,] <- .Fortran(C_getmsni0,
as.double(ni[,indi]),
as.integer(n),
as.integer(lindi),
msni0=double(n))$msni0
}
dim(msth0) <- dim(msni0) <- c(nbv+1,n)
dim(z$msth) <- dim(z$msni) <- c(nbv+1,dtheta)
list(mstheta=z$msth,msni=z$msni,msth0=msth0,msni0=msni0)
}
lkfulls0 <- function(h,vext,n){
z <- .Fortran(C_lkfuls0,
as.double(h),
as.double(vext),
ind=integer(3*n),
w=double(n),
n=as.integer(n))[c("ind","w","n")]
dim(z$ind) <- c(3,n)
list(h=h,ind=z$ind[,1:z$n],w=z$w[1:z$n],nind=z$n)
}
lkfullse3 <- function(h,kappa,gradstats,vext,n){
ngrad <- dim(gradstats$bghat)[2]
if(length(h)<ngrad) h <- rep(h[1],ngrad)
z <- .Fortran(C_lkfulse3,
as.double(h),
as.double(kappa),
as.double(gradstats$k456),
as.integer(ngrad),
as.double(vext),
ind=integer(5*n),
w=double(n),
n=as.integer(n),
as.integer(gradstats$dist))[c("ind","w","n")]
dim(z$ind) <- c(5,n)
list(h=h,kappa=kappa,ind=z$ind[,1:z$n],w=z$w[1:z$n],nind=z$n)
}
reduceparam <- function(param){
ind <- param$ind[4,]==param$ind[5,]
param$ind <- param$ind[,ind]
param$w <- param$w[ind]
param$n <- sum(ind)
h <- max(param$h)
rind <- param$ind[1,]*h*2*h*2+param$ind[2,]*h*2+param$ind[3,]
oind <- order(rind)
param$ind <- param$ind[,oind]
param$w <- param$w[oind]
starts <- cumsum(rle(rind[oind])$lengths)
param$nstarts <- length(starts)
param$starts <- c(0,starts)
param
}
betagamma <- function(g){
dg <- dim(g)
ngrad <- if(!is.null(dg)) dg[2] else 1
bghat <- .Fortran(C_bgstats,
as.double(g),
as.integer(ngrad),
double(2*ngrad),
bghat = double(2*ngrad*ngrad))$bghat
dim(bghat) <- c(2, ngrad, ngrad)
## sphaerische Coordinaten fuer Gradienten-Paare
bghat
}
matrm <- function(b, g){
matrix(c(cos(b), 0, sin(b),
sin(b)*sin(g), cos(g), -cos(b)*sin(g),
-cos(g)*sin(b), sin(g), cos(b)*cos(g)),
3, 3)
}
#
#
#
getkappas <- function(grad, trace = 0, dist = 1){
#
# dist = 1: k4^2+k5^2+|k6|
# dist = 2: k4^2+k5^2+k6^2
# dist = 3: acos(g_i%*%g_j)
# new version with analytic expm(par[1]*mx)
#
krit <- function(par, matm, beta){
## sum((matm-expm(par[1]*m4)%*%expm(par[2]*m5)%*%expm(par[3]*m6))^2)
.Fortran(C_k456krb,
as.double(par),
as.double(beta),
as.double(matm),
erg = double(1))$erg
}
krit5 <- function(x,p,pk4,matm,beta){
# for line search with respect to k5 to get second solution
p1 <- p+c(pk4/2,x,pi)
krit(p1,matm,beta)
}
ngrad <- dim(grad)[2]
if(dist<3){
prta <- Sys.time()
cat("Start computing spherical distances", format(Sys.time()), "\n")
kappa456 <- kappa456a <- array(0, c(3, ngrad, ngrad))
bghat <- betagamma(grad)
for (i in 1:ngrad) for (j in 1:ngrad) {
bg <- bghat[, i, j]
# fix for discontinuity
if(abs(cos(bg[1])) < 1.e-6) bg[1] = pi/2 - 1e-6*sign(cos(bg[1]))
matm <- matrm(bg[1], bg[2])
cbg1 <- cos(bg[1])
k456 <- runif(3, -.1, .1)
maxit <- 10000
fnscale <- 1
z <- optim(k456, krit, method = "BFGS", matm = matm, beta=bg[1],
control = list(trace = trace, fnscale=fnscale, maxit=maxit, reltol = 1e-12, abstol = 1e-16))
count <- 5
while (z$value > 1e-14&count>0) {
## cat("i",i,"j",j,"value",z$value,"par",z$par,"\n")
maxit <- maxit*1.5
fnscale <- fnscale/2
k456 <- runif(3, -1, 1)
z <- optim(k456, krit, method = "BFGS", matm = matm, beta=bg[1],
control = list(trace = trace, fnscale=fnscale, maxit=maxit, reltol = 1e-12, abstol = 1e-16))
## cat(" new value",z$value,"par",z$par,"\n")
count <- count - 1
if(count==0)
cat("failed for bg:",bg,"value",z$value,"par",z$par,"\n")
}
z$par[1] <- z$par[1]/cbg1
kappa456[, i, j] <- z$par
pk4 <- abs(2*pi*cbg1/(2-cbg1^2)^.5)
if(kappa456[1,i,j] < -pk4/2) kappa456[1,i,j] <- kappa456[1,i,j] - trunc(kappa456[1, i, j]/pk4-1)*pk4
if(kappa456[1,i,j] >= pk4/2) kappa456[1,i,j] <- kappa456[1,i,j] - trunc(kappa456[1, i, j]/pk4+1)*pk4
if(kappa456[2,i,j] < -pi) kappa456[2,i,j] <- kappa456[2,i,j] - trunc(kappa456[2, i, j]/2/pi-1)*2*pi
if(kappa456[2,i,j] > pi) kappa456[2,i,j] <- kappa456[2,i,j] - trunc(kappa456[2, i, j]/2/pi+1)*2*pi
if(kappa456[3,i,j] < -pi) kappa456[3,i,j] <- kappa456[3,i,j] - trunc(kappa456[3, i, j]/2/pi-1)*2*pi
if(kappa456[3,i,j] > pi) kappa456[3,i,j] <- kappa456[3,i,j] - trunc(kappa456[3, i, j]/2/pi+1)*2*pi
kappa456a[,i,j] <- kappa456[,i,j]+c(pk4/2,0,pi)
kpar <- kappa456[,i,j]*c(cbg1,1,1)
kappa456a[2,i,j] <- optimize(krit5,c(-pi,pi),p=kpar,pk4=pk4,
matm=matm,beta=bg[1],maximum=FALSE)$minimum
if(kappa456a[1,i,j] < -pk4/2) kappa456a[1,i,j] <- kappa456a[1,i,j] - trunc(kappa456a[1, i, j]/pk4-1)*pk4
if(kappa456a[1,i,j] >= pk4/2) kappa456a[1,i,j] <- kappa456a[1,i,j] - trunc(kappa456a[1, i, j]/pk4+1)*pk4
if(kappa456a[2,i,j] < -pi) kappa456a[2,i,j] <- kappa456a[2,i,j] - trunc(kappa456a[2, i, j]/2/pi-1)*2*pi
if(kappa456a[2,i,j] > pi) kappa456a[2,i,j] <- kappa456a[2,i,j] - trunc(kappa456a[2, i, j]/2/pi+1)*2*pi
if(kappa456a[3,i,j] < -pi) kappa456a[3,i,j] <- kappa456a[3,i,j] - trunc(kappa456a[3, i, j]/2/pi-1)*2*pi
if(kappa456a[3,i,j] > pi) kappa456a[3,i,j] <- kappa456a[3,i,j] - trunc(kappa456a[3, i, j]/2/pi+1)*2*pi
}
dka <- switch(dist,kappa456[1,,]^2+kappa456[2,,]^2+abs(kappa456[3,,]),
kappa456[1,,]^2+kappa456[2,,]^2+kappa456[3,,]^2)
dkb <- switch(dist,kappa456a[1,,]^2+kappa456a[2,,]^2+abs(kappa456a[3,,]),
kappa456a[1,,]^2+kappa456a[2,,]^2+kappa456a[3,,]^2)
dim(kappa456) <- dim(kappa456a) <- c(3,ngrad*ngrad)
kappa456[,dkb<dka] <- kappa456a[,dkb<dka]
dim(kappa456) <- c(3,ngrad,ngrad)
prtb <- Sys.time()
cat("End computing spherical distances", format(Sys.time()), "\n")
} else {
kappa456 <- array(0, c(3, ngrad, ngrad))
bghat <- betagamma(grad)
for(i in 1:ngrad) kappa456[1,i,] <- acos(pmin(1,abs(grad[,i]%*%grad)))
}
list(k456 = kappa456, bghat = bghat, dist=dist)
}
suggestkappa <- function(grad,vred=1,dist=1){
#
# get a kappa value from variance reduction on the sphere
#
gstats <- getkappas(grad,dist=dist)
ngrad <- dim(grad)[2]
vred <- min(vred,ngrad-1)
d <- switch(dist,apply(gstats$k456[1:2,,]^2,2:3,sum)+abs(gstats$k456[3,,]),
apply(gstats$k456^2,2:3,sum),
apply(gstats$k456^2,2:3,sum))
kmin <- sqrt(min(d[d>1e-8]))# just to prevent from taking zero
kappa <- kmin
vredk <- 1
while(vredk < vred){
kappa <- kappa*1.005
w <- matrix(pmax(1-d/kappa^2,0),ngrad,ngrad)
vredk <- mean(apply(w,1,sum)^2/apply(w^2,1,sum))
}
list(kappa=kappa,vred=vredk)
}
WIAS-BERLIN/aws documentation built on Sept. 10, 2023, 6:20 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions suggestkappa getkappas matrm betagamma reduceparam lkfullse3 lkfulls0 interpolatesphere1 getkappasmsh3 getkappas3 getkappas getkappasmsh lkfullse3msh interpolatesphere getnext3g unifybvals getsphwghts sphtrarea1
sphtrarea1 <- function(g1,g2,g3){
## Compute area of sherical triangle spanned by vectors
## g1,g2,g3 on unit sphere
## use absolute values to identify opposite directions with each other
c12 <- abs(g1%*%g2)
c13 <- abs(g1%*%g3)
c23 <- abs(g2%*%g3)
s12 <- sqrt(1-c12^2)
s13 <- sqrt(1-c13^2)
s23 <- sqrt(1-c23^2)
a12 <- acos(c12)
a23 <- acos(c23)
a13 <- acos(c13)
## observed b1, b2, b3 outside [-1,1] due to numerics
if(s12*s13*s23>1e-8){
b1 <- min(max(1+(c23-cos(a12-a13))/s12/s13,-1),1)
b2 <- min(max(1+(c13-cos(a12-a23))/s12/s23,-1),1)
b3 <- min(max(1+(c12-cos(a13-a23))/s13/s23,-1),1)
acos(b1)+acos(b2)+acos(b3)-pi
} else {
0
}
}
getsphwghts <- function(g,g1,g2,g3){
##
## compute weights for linear interpolation in g using g1,g2,g3
##
ierr <- 0
w <- numeric(3)
w0 <- sphtrarea1(g1,g2,g3)
w[1] <- sphtrarea1(g,g2,g3)
w[2] <- sphtrarea1(g,g1,g3)
w[3] <- sphtrarea1(g,g1,g2)
if(w0 < sum(w) - 1e-6){
#cat("gradient does not belong to triangle")
#print(g)
#print(cbind(g1,g2,g3))
#cat("w0",w0,"w",w,"alt",w[1]+w[2]-w[3],"\n")
ierr <- 1
}
list(w=w/sum(w),ierr=ierr)
}
unifybvals <- function(bval,dbv=51){
nbv <- length(bval)
nbval <- bval
obval <- numeric(nbv)
while(any(nbval!=obval)){
obval <- nbval
sbv <- sort(obval)
dsbv <- (1:(nbv-1))[diff(sbv)<dbv]
sbv[dsbv+1] <- sbv[dsbv]
obv <- order(obval)
nbval[obv] <- sbv
}
for(bv in unique(nbval)) nbval[nbval==bv] <- trunc(mean(bval[nbval==bv]))
nbval
}
getnext3g <- function(grad,bv){
##
## calculate next neighbors on the sphere and weights for interpolation
##
binomial <- function(a,b) prod(1:a)/prod(1:(a-b))/prod(1:b)
##
grad <- grad[,bv>0]
ng <- dim(grad)[2]
ng2 <- max(4,ng/6)
perm <- matrix(1,(ng2-2)*(ng2-1)*ng2/6,3)
l <- 1
for(k in 1:(ng2-2)){
for(j in (k+1):(ng2-1)){
for(i in (j+1):ng2){
perm[l,] <- c(k,j,i)
l <- l+1
}
}
}
dperm <- l-1
perm <- perm[2:dperm,]
# thats all ordered triples from (1:ng2), without (1,2,3)
bv <- unifybvals(bv[bv>0]) # identifies b-values that differ by less than dbv=51
ubv <- sort(unique(bv[bv>max(bv)/50]))
nbv <- length(ubv)
ind <- array(0,c(3,nbv,ng))
w <- array(0,c(3,nbv,ng))
for(i in 1:nbv){
indb <- (1:ng)[bv==ubv[i]]
ind[1,i,indb] <- indb
ind[2,i,indb] <- indb
ind[3,i,indb] <- indb
w[1,i,indb] <- 1
w[2,i,indb] <- 0
w[3,i,indb] <- 0
for(j in (1:nbv)[-i]){
indbk <- (1:ng)[bv==ubv[j]]
perm0 <- perm
if(length(indbk)<ng2){
indp <- apply(perm < length(indbk),1,all)
perm0 <- perm0[indp,]
}
dperm <- dim(perm0)[1]+1
for(k in indb){
d <- abs(t(grad[,k])%*%grad[,indbk])
od <- order(d,decreasing = TRUE)
ijk0 <- (1:length(indbk))[od]
ijk <- indbk[od]
# cat("k",k,"d",d[od][1:8],"ijk",ijk[1:8],"\n")
ind[,j,k] <- ijk[1:3]
if(max(d)>1-1e-6){
w[,j,k] <- c(1,0,0)
} else {
z <- getsphwghts(grad[,k],grad[,ijk[1]],grad[,ijk[2]],grad[,ijk[3]])
l <- 1
if(z$ierr==1){
# order triplets in perm according
dd <- numeric(dperm-1)
d <- acos(pmin(1,d))
for(l in 1:(dperm-1)){
ijk1 <- ijk0[perm0[l,]]
dd[l] <- d[ijk1[1]]+d[ijk1[2]]+d[ijk1[3]]
}
odd <- order(dd)
l <- 1
}
while(z$ierr==1&&l<dperm){
ijk1 <- ijk[perm0[odd[l],]]
z <- getsphwghts(grad[,k],grad[,ijk1[1]],grad[,ijk1[2]],grad[,ijk1[3]])
ind[,j,k] <- ijk1[1:3]
l<-l+1
}
## if we were running out of options use ijk[1:3]
if(z$ierr==1){
ind[,j,k] <- ijk[1:3]
z <- getsphwghts(grad[,k],grad[,ijk[1]],grad[,ijk[2]],grad[,ijk[3]])
}
w[,j,k] <- z$w
}
}
}
}
# spheres identified by bvalues in ubv
# ind[j,k,] contains indices of gradients to be used in spherical interpolation
# on shell i for gradient k
# w[j,k,] contains the corresponding weights
list(ind=ind, w=w, ubv = ubv, nbv = nbv, bv=bv)
}
interpolatesphere <- function(theta,n3g){
## interpolate estimated thetas to get values on all spheres
## n3g generated by function getnext3g
## dim(theta) = c(n1,n2,n3,ngrad)
dtheta <- dim(theta)
mstheta <- array(0,c(n3g$nbv,dtheta))
for(i in 1:n3g$nbv){
for(j in 1:dim(theta)[2]){
mstheta[i,,j] <- theta[,n3g$ind[,i,j]]%*%n3g$w[,i,j]
}
}
mstheta
}
lkfullse3msh <- function(h,kappa,gradstats,vext,n){
nbv <- gradstats$nbv
dist <- gradstats$dist
ngrad <- gradstats$ngrad
bvind <- gradstats$bvind
ind <- matrix(0,5,n)
w <- numeric(n)
nn <- 0
dist <- 4
for(i in 1:nbv){
gshell <- list(k456=gradstats$k456[[i]],bghat=gradstats$bghat[[i]],dist=dist)
z <- lkfullse3(h[bvind[[i]]],kappa[bvind[[i]]],gshell,vext,n)
ind[1:3,nn+1:z$nind] <- z$ind[1:3,]
ind[4:5,nn+1:z$nind] <- matrix(bvind[[i]][z$ind[4:5,]],2, z$nind)
#
# convert indeces on selected shell to total indeces
#
w[nn+1:z$nind] <- z$w
nn <- nn + z$nind
}
list(h=h,kappa=kappa,ind=ind[,1:nn],w=w[1:nn],nind=nn)
}
gethseqfullse3msh <- function (kstar, gradstats, kappa, vext = c(1, 1))
{
#
# generate information on local bandwidths and variance reduction
# for smoothing on multiple shells
#
nbv <- gradstats$nbv
ngrad <- gradstats$ngrad
h <- vr <- matrix(1,ngrad,kstar+1)
n <- 0
dist <- 4
for(i in 1:nbv){
gshell <- list(k456=gradstats$k456[[i]],bghat=gradstats$bghat[[i]],dist=dist)
z <- gethseqfullse3(kstar, gshell, kappa=kappa, vext=vext)
h[gradstats$bvind[[i]],-1] <- z$h
vr[gradstats$bvind[[i]],-1] <- z$vred
vr[gradstats$bvind[[i]],1] <- vr[gradstats$bvind[[i]],2]/1.25
n <- n+z$n
}
cat("\n total number of positive weights:",n,"mean maximal bandwidth",signif(mean(h[,kstar]),3), "\n")
list(h=h,vred=vr,n=n)
}
gethseqfullse3 <- function (kstar, gradstats, kappa=NULL, vext = c(1, 1))
{
ngrad <- dim(gradstats$bghat)[2]
h <- vr <- matrix(0,ngrad,kstar)
if(is.null(kappa)) kappa <- pi/ngrad
# if(length(kappa)<kstar) kappa <- rep(kappa[1],kstar)
prt0 <- Sys.time()
cat("get sequences of bw, kappa up to kstar=", kstar, " ")
n <- 0
for(i in 1:ngrad){
z <- .Fortran(C_ghfse3i,
as.integer(i),#i4
as.integer(kstar),#kstar
as.double(gradstats$k456),
as.integer(ngrad),
as.double(kappa),#kappa
as.double(vext),#vext
h=double(kstar),
vr=double(kstar),#
n=integer(1),#
as.integer(gradstats$dist))[c("h","vr","n")]
h[i,] <- z$h
vr[i,] <- z$vr
n <- n+z$n
cat(".")
}
cat("\n number of positive weights:",n,"mean maximal bandwidth",signif(mean(h[,kstar]),3),"time elapsed:",
format(difftime(Sys.time(), prt0), digits = 3),
"\n")
list(h=h,kappa=kappa,vred=vr,n=n)
}
getkappasmsh <- function(grad, msstructure, trace = 0, dist = 1){
ngrad <- dim(grad)[2]
nbv <- msstructure$nbv
bv <- msstructure$bv
ubv <- msstructure$ubv
bvind <- k456 <- bghat <- list(NULL)
for(i in 1:nbv){
#
# collect information for spherical distances on each schell separately
#
ind <- (1:ngrad)[bv==ubv[i]]
z <- getkappas(grad[,ind], trace = trace, dist = dist)
bvind[[i]] <- ind
k456[[i]] <- z$k456
bghat[[i]] <- z$bghat
}
list(k456 = k456, bghat = bghat, bvind = bvind, dist=dist, nbv = nbv, ngrad = ngrad)
}
getkappas <- function(grad, trace = 0, dist = 1){
#
# dist = 1: k4^2+k5^2+|k6|
# dist = 2: k4^2+k5^2+k6^2
# dist = 3: acos(g_i%*%g_j)
# new version with analytic expm(par[1]*mx)
#
krit <- function(par, matm, beta){
## sum((matm-expm(par[1]*m4)%*%expm(par[2]*m5)%*%expm(par[3]*m6))^2)
.Fortran(C_k456krb,
as.double(par),
as.double(beta),
as.double(matm),
erg = double(1))$erg
}
krit5 <- function(x,p,pk4,matm,beta){
# for line search with respect to k5 to get second solution
p1 <- p+c(pk4/2,x,pi)
krit(p1,matm,beta)
}
ngrad <- dim(grad)[2]
if(dist<3){
prta <- Sys.time()
cat("Start computing spherical distances", format(Sys.time()), "\n")
kappa456 <- kappa456a <- array(0, c(3, ngrad, ngrad))
bghat <- betagamma(grad)
for (i in 1:ngrad) for (j in 1:ngrad) {
bg <- bghat[, i, j]
# fix for discontinuity
if(abs(cos(bg[1])) < 1.e-6) bg[1] = pi/2 - 1e-6*sign(cos(bg[1]))
matm <- matrm(bg[1], bg[2])
cbg1 <- cos(bg[1])
k456 <- runif(3, -.1, .1)
maxit <- 10000
fnscale <- 1
z <- optim(k456, krit, method = "BFGS", matm = matm, beta=bg[1],
control = list(trace = trace, fnscale=fnscale, maxit=maxit, reltol = 1e-12, abstol = 1e-16))
count <- 5
while (z$value > 1e-14&count>0) {
## cat("i",i,"j",j,"value",z$value,"par",z$par,"\n")
maxit <- maxit*1.5
fnscale <- fnscale/2
k456 <- runif(3, -1, 1)
z <- optim(k456, krit, method = "BFGS", matm = matm, beta=bg[1],
control = list(trace = trace, fnscale=fnscale, maxit=maxit, reltol = 1e-12, abstol = 1e-16))
## cat(" new value",z$value,"par",z$par,"\n")
count <- count - 1
if(count==0)
cat("failed for bg:",bg,"value",z$value,"par",z$par,"\n")
}
z$par[1] <- z$par[1]/cbg1
kappa456[, i, j] <- z$par
pk4 <- abs(2*pi*cbg1/(2-cbg1^2)^.5)
if(kappa456[1,i,j] < -pk4/2) kappa456[1,i,j] <- kappa456[1,i,j] - trunc(kappa456[1, i, j]/pk4-1)*pk4
if(kappa456[1,i,j] >= pk4/2) kappa456[1,i,j] <- kappa456[1,i,j] - trunc(kappa456[1, i, j]/pk4+1)*pk4
if(kappa456[2,i,j] < -pi) kappa456[2,i,j] <- kappa456[2,i,j] - trunc(kappa456[2, i, j]/2/pi-1)*2*pi
if(kappa456[2,i,j] > pi) kappa456[2,i,j] <- kappa456[2,i,j] - trunc(kappa456[2, i, j]/2/pi+1)*2*pi
if(kappa456[3,i,j] < -pi) kappa456[3,i,j] <- kappa456[3,i,j] - trunc(kappa456[3, i, j]/2/pi-1)*2*pi
if(kappa456[3,i,j] > pi) kappa456[3,i,j] <- kappa456[3,i,j] - trunc(kappa456[3, i, j]/2/pi+1)*2*pi
kappa456a[,i,j] <- kappa456[,i,j]+c(pk4/2,0,pi)
kpar <- kappa456[,i,j]*c(cbg1,1,1)
kappa456a[2,i,j] <- optimize(krit5,c(-pi,pi),p=kpar,pk4=pk4,
matm=matm,beta=bg[1],maximum=FALSE)$minimum
if(kappa456a[1,i,j] < -pk4/2) kappa456a[1,i,j] <- kappa456a[1,i,j] - trunc(kappa456a[1, i, j]/pk4-1)*pk4
if(kappa456a[1,i,j] >= pk4/2) kappa456a[1,i,j] <- kappa456a[1,i,j] - trunc(kappa456a[1, i, j]/pk4+1)*pk4
if(kappa456a[2,i,j] < -pi) kappa456a[2,i,j] <- kappa456a[2,i,j] - trunc(kappa456a[2, i, j]/2/pi-1)*2*pi
if(kappa456a[2,i,j] > pi) kappa456a[2,i,j] <- kappa456a[2,i,j] - trunc(kappa456a[2, i, j]/2/pi+1)*2*pi
if(kappa456a[3,i,j] < -pi) kappa456a[3,i,j] <- kappa456a[3,i,j] - trunc(kappa456a[3, i, j]/2/pi-1)*2*pi
if(kappa456a[3,i,j] > pi) kappa456a[3,i,j] <- kappa456a[3,i,j] - trunc(kappa456a[3, i, j]/2/pi+1)*2*pi
}
dka <- switch(dist,kappa456[1,,]^2+kappa456[2,,]^2+abs(kappa456[3,,]),
kappa456[1,,]^2+kappa456[2,,]^2+kappa456[3,,]^2)
dkb <- switch(dist,kappa456a[1,,]^2+kappa456a[2,,]^2+abs(kappa456a[3,,]),
kappa456a[1,,]^2+kappa456a[2,,]^2+kappa456a[3,,]^2)
dim(kappa456) <- dim(kappa456a) <- c(3,ngrad*ngrad)
kappa456[,dkb<dka] <- kappa456a[,dkb<dka]
dim(kappa456) <- c(3,ngrad,ngrad)
prtb <- Sys.time()
cat("End computing spherical distances", format(Sys.time()), "\n")
} else {
kappa456 <- array(0, c(3, ngrad, ngrad))
bghat <- betagamma(grad)
for(i in 1:ngrad) kappa456[1,i,] <- acos(pmin(1,abs(grad[,i]%*%grad)))
}
list(k456 = kappa456, bghat = bghat, dist=dist)
}
getkappas3 <- function(grad, trace = 0){
#
# dist: acos(g_i%*%g_j)
#
ngrad <- dim(grad)[2]
kappa456 <- array(0, c(3, ngrad, ngrad))
bghat <- betagamma(grad)
for(i in 1:ngrad) kappa456[1,i,] <- acos(pmin(1,abs(grad[,i]%*%grad)))
list(k456 = kappa456, bghat = bghat, dist=3)
}
getkappasmsh3 <- function(grad, msstructure, trace = 0){
ngrad <- dim(grad)[2]
nbv <- msstructure$nbv
bv <- msstructure$bv
ubv <- msstructure$ubv
bvind <- k456 <- bghat <- list(NULL)
for(i in 1:nbv){
#
# collect information for spherical distances on each schell separately
#
ind <- (1:ngrad)[bv==ubv[i]]
z <- getkappas3(grad[,ind], trace = trace)
bvind[[i]] <- ind
k456[[i]] <- z$k456
bghat[[i]] <- z$bghat
}
list(k456 = k456, bghat = bghat, bvind = bvind, dist=3, nbv = nbv, ngrad = ngrad)
}
interpolatesphere1 <- function(theta,th0,ni,ni0,n3g){
## interpolate estimated thetas to get values on all spheres
## n3g generated by function getnext3g
## dim(theta) = c(n1,n2,n3,ngrad)
nbv <- n3g$nbv
bv <- n3g$bv
ubv <- n3g$ubv
dtheta <- dim(theta)
dth0 <- dim(th0)
ng <- dtheta[2]
n <- dtheta[1]
#t1 <- Sys.time()
z <- .Fortran(C_ipolsp1,
as.double(theta),
as.double(th0),
as.double(ni),
as.double(ni0),
as.integer(n),
as.integer(ng),
as.integer(n3g$ind),
as.double(n3g$w),
as.integer(nbv),
as.integer(nbv+1),
msth=double((nbv+1)*n*ng),
msni=double((nbv+1)*n*ng))[c("msth","msni")]
#cat("time for sb-interpolation", format(difftime(Sys.time(),t1),digits=3),"\n")
# now fill vector for s0
dim(theta) <- dim(ni) <- c(n,ng)
msth0 <- msni0 <- array(0,c(nbv+1,n))
msth0[1,] <- th0
msni0[1,] <- ni0
for(i in 1:nbv){
indi <- (1:ng)[bv==ubv[i]]
lindi <- length(indi)
# msth0[i+1,mask] <- theta[mask,indi]%*%rep(1/lindi,lindi)
msth0[i+1,] <- .Fortran(C_getmsth0,
as.double(theta[,indi]),
as.integer(n),
as.integer(lindi),
msth0=double(n))$msth0
# correct value would be
# msni0[i+1,] <- 1/(1/(ni[,indi])%*%rep(1/lindi,lindi)^2)
# try to be less conservative by ignorin squares in w
# msni0[i+1,mask] <- 1/((1/ni[mask,indi])%*%rep(1/lindi,lindi))
msni0[i+1,] <- .Fortran(C_getmsni0,
as.double(ni[,indi]),
as.integer(n),
as.integer(lindi),
msni0=double(n))$msni0
}
dim(msth0) <- dim(msni0) <- c(nbv+1,n)
dim(z$msth) <- dim(z$msni) <- c(nbv+1,dtheta)
list(mstheta=z$msth,msni=z$msni,msth0=msth0,msni0=msni0)
}
lkfulls0 <- function(h,vext,n){
z <- .Fortran(C_lkfuls0,
as.double(h),
as.double(vext),
ind=integer(3*n),
w=double(n),
n=as.integer(n))[c("ind","w","n")]
dim(z$ind) <- c(3,n)
list(h=h,ind=z$ind[,1:z$n],w=z$w[1:z$n],nind=z$n)
}
lkfullse3 <- function(h,kappa,gradstats,vext,n){
ngrad <- dim(gradstats$bghat)[2]
if(length(h)<ngrad) h <- rep(h[1],ngrad)
z <- .Fortran(C_lkfulse3,
as.double(h),
as.double(kappa),
as.double(gradstats$k456),
as.integer(ngrad),
as.double(vext),
ind=integer(5*n),
w=double(n),
n=as.integer(n),
as.integer(gradstats$dist))[c("ind","w","n")]
dim(z$ind) <- c(5,n)
list(h=h,kappa=kappa,ind=z$ind[,1:z$n],w=z$w[1:z$n],nind=z$n)
}
reduceparam <- function(param){
ind <- param$ind[4,]==param$ind[5,]
param$ind <- param$ind[,ind]
param$w <- param$w[ind]
param$n <- sum(ind)
h <- max(param$h)
rind <- param$ind[1,]*h*2*h*2+param$ind[2,]*h*2+param$ind[3,]
oind <- order(rind)
param$ind <- param$ind[,oind]
param$w <- param$w[oind]
starts <- cumsum(rle(rind[oind])$lengths)
param$nstarts <- length(starts)
param$starts <- c(0,starts)
param
}
betagamma <- function(g){
dg <- dim(g)
ngrad <- if(!is.null(dg)) dg[2] else 1
bghat <- .Fortran(C_bgstats,
as.double(g),
as.integer(ngrad),
double(2*ngrad),
bghat = double(2*ngrad*ngrad))$bghat
dim(bghat) <- c(2, ngrad, ngrad)
## sphaerische Coordinaten fuer Gradienten-Paare
bghat
}
matrm <- function(b, g){
matrix(c(cos(b), 0, sin(b),
sin(b)*sin(g), cos(g), -cos(b)*sin(g),
-cos(g)*sin(b), sin(g), cos(b)*cos(g)),
3, 3)
}
#
#
#
getkappas <- function(grad, trace = 0, dist = 1){
#
# dist = 1: k4^2+k5^2+|k6|
# dist = 2: k4^2+k5^2+k6^2
# dist = 3: acos(g_i%*%g_j)
# new version with analytic expm(par[1]*mx)
#
krit <- function(par, matm, beta){
## sum((matm-expm(par[1]*m4)%*%expm(par[2]*m5)%*%expm(par[3]*m6))^2)
.Fortran(C_k456krb,
as.double(par),
as.double(beta),
as.double(matm),
erg = double(1))$erg
}
krit5 <- function(x,p,pk4,matm,beta){
# for line search with respect to k5 to get second solution
p1 <- p+c(pk4/2,x,pi)
krit(p1,matm,beta)
}
ngrad <- dim(grad)[2]
if(dist<3){
prta <- Sys.time()
cat("Start computing spherical distances", format(Sys.time()), "\n")
kappa456 <- kappa456a <- array(0, c(3, ngrad, ngrad))
bghat <- betagamma(grad)
for (i in 1:ngrad) for (j in 1:ngrad) {
bg <- bghat[, i, j]
# fix for discontinuity
if(abs(cos(bg[1])) < 1.e-6) bg[1] = pi/2 - 1e-6*sign(cos(bg[1]))
matm <- matrm(bg[1], bg[2])
cbg1 <- cos(bg[1])
k456 <- runif(3, -.1, .1)
maxit <- 10000
fnscale <- 1
z <- optim(k456, krit, method = "BFGS", matm = matm, beta=bg[1],
control = list(trace = trace, fnscale=fnscale, maxit=maxit, reltol = 1e-12, abstol = 1e-16))
count <- 5
while (z$value > 1e-14&count>0) {
## cat("i",i,"j",j,"value",z$value,"par",z$par,"\n")
maxit <- maxit*1.5
fnscale <- fnscale/2
k456 <- runif(3, -1, 1)
z <- optim(k456, krit, method = "BFGS", matm = matm, beta=bg[1],
control = list(trace = trace, fnscale=fnscale, maxit=maxit, reltol = 1e-12, abstol = 1e-16))
## cat(" new value",z$value,"par",z$par,"\n")
count <- count - 1
if(count==0)
cat("failed for bg:",bg,"value",z$value,"par",z$par,"\n")
}
z$par[1] <- z$par[1]/cbg1
kappa456[, i, j] <- z$par
pk4 <- abs(2*pi*cbg1/(2-cbg1^2)^.5)
if(kappa456[1,i,j] < -pk4/2) kappa456[1,i,j] <- kappa456[1,i,j] - trunc(kappa456[1, i, j]/pk4-1)*pk4
if(kappa456[1,i,j] >= pk4/2) kappa456[1,i,j] <- kappa456[1,i,j] - trunc(kappa456[1, i, j]/pk4+1)*pk4
if(kappa456[2,i,j] < -pi) kappa456[2,i,j] <- kappa456[2,i,j] - trunc(kappa456[2, i, j]/2/pi-1)*2*pi
if(kappa456[2,i,j] > pi) kappa456[2,i,j] <- kappa456[2,i,j] - trunc(kappa456[2, i, j]/2/pi+1)*2*pi
if(kappa456[3,i,j] < -pi) kappa456[3,i,j] <- kappa456[3,i,j] - trunc(kappa456[3, i, j]/2/pi-1)*2*pi
if(kappa456[3,i,j] > pi) kappa456[3,i,j] <- kappa456[3,i,j] - trunc(kappa456[3, i, j]/2/pi+1)*2*pi
kappa456a[,i,j] <- kappa456[,i,j]+c(pk4/2,0,pi)
kpar <- kappa456[,i,j]*c(cbg1,1,1)
kappa456a[2,i,j] <- optimize(krit5,c(-pi,pi),p=kpar,pk4=pk4,
matm=matm,beta=bg[1],maximum=FALSE)$minimum
if(kappa456a[1,i,j] < -pk4/2) kappa456a[1,i,j] <- kappa456a[1,i,j] - trunc(kappa456a[1, i, j]/pk4-1)*pk4
if(kappa456a[1,i,j] >= pk4/2) kappa456a[1,i,j] <- kappa456a[1,i,j] - trunc(kappa456a[1, i, j]/pk4+1)*pk4
if(kappa456a[2,i,j] < -pi) kappa456a[2,i,j] <- kappa456a[2,i,j] - trunc(kappa456a[2, i, j]/2/pi-1)*2*pi
if(kappa456a[2,i,j] > pi) kappa456a[2,i,j] <- kappa456a[2,i,j] - trunc(kappa456a[2, i, j]/2/pi+1)*2*pi
if(kappa456a[3,i,j] < -pi) kappa456a[3,i,j] <- kappa456a[3,i,j] - trunc(kappa456a[3, i, j]/2/pi-1)*2*pi
if(kappa456a[3,i,j] > pi) kappa456a[3,i,j] <- kappa456a[3,i,j] - trunc(kappa456a[3, i, j]/2/pi+1)*2*pi
}
dka <- switch(dist,kappa456[1,,]^2+kappa456[2,,]^2+abs(kappa456[3,,]),
kappa456[1,,]^2+kappa456[2,,]^2+kappa456[3,,]^2)
dkb <- switch(dist,kappa456a[1,,]^2+kappa456a[2,,]^2+abs(kappa456a[3,,]),
kappa456a[1,,]^2+kappa456a[2,,]^2+kappa456a[3,,]^2)
dim(kappa456) <- dim(kappa456a) <- c(3,ngrad*ngrad)
kappa456[,dkb<dka] <- kappa456a[,dkb<dka]
dim(kappa456) <- c(3,ngrad,ngrad)
prtb <- Sys.time()
cat("End computing spherical distances", format(Sys.time()), "\n")
} else {
kappa456 <- array(0, c(3, ngrad, ngrad))
bghat <- betagamma(grad)
for(i in 1:ngrad) kappa456[1,i,] <- acos(pmin(1,abs(grad[,i]%*%grad)))
}
list(k456 = kappa456, bghat = bghat, dist=dist)
}
suggestkappa <- function(grad,vred=1,dist=1){
#
# get a kappa value from variance reduction on the sphere
#
gstats <- getkappas(grad,dist=dist)
ngrad <- dim(grad)[2]
vred <- min(vred,ngrad-1)
d <- switch(dist,apply(gstats$k456[1:2,,]^2,2:3,sum)+abs(gstats$k456[3,,]),
apply(gstats$k456^2,2:3,sum),
apply(gstats$k456^2,2:3,sum))
kmin <- sqrt(min(d[d>1e-8]))# just to prevent from taking zero
kappa <- kmin
vredk <- 1
while(vredk < vred){
kappa <- kappa*1.005
w <- matrix(pmax(1-d/kappa^2,0),ngrad,ngrad)
vredk <- mean(apply(w,1,sum)^2/apply(w^2,1,sum))
}
list(kappa=kappa,vred=vredk)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.