Nothing
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)
}
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.