#' Internal_functions
#'
#' Here are reported a collection of internal functions
#' @author Paolo Piras, Antonio Profico
#' @export
energies2d<-function(source,target,mag=1,mu1=1,mu2=1,doopa=T,links=NULL,PB=NA,PA=NA,S=NA,SB=NA,H=NA,V=3,a=NULL,q=NULL,Y=FALSE,j=FALSE,D=FALSE,St=Inf,Q=TRUE,graph=T,body=T){
library(shapes)
library(Morpho)
library(matrixcalc)
library(RTriangle)
library(ggm)
library(geometry)
library(tripack)
library(geoR)
library(gstat)
library(sp)
library(fields)
mate<-source
mate2<-target
mate2<-mate+(mate2-mate)*mag
if(doopa==T){
theopa<-rotonto(mate,mate2,scale=F,reflection=F)
mate<-mate
mate2<-theopa$yrot}
matr<-mate
if(!is.null(links)&is.na(S)==T){S<-list2matrix(links)}
po<-areasip(matr,links=links,PB=PB,PA=PA,S=S,SB=SB,H=H,V=V,a=a,q=q,Y=Y,j=j,D=D,St=St,Q=Q)
ip<-po$ip
areas<-po$areas
M<-po$centros
myl<-tpsdry2(mate,mate2,doopa=F)
m<-dim(mate)[2]
WtY<-t(myl$W)%*%mate2
trace<-c(0)
for (i in 1:m) {
trace <- trace + WtY[i, i]
}
be <- (16 * pi) * trace
ge<-gen1(mate,mate2,mu1=mu1,mu2=mu2)
veclist<-NULL
for(j in 1:nrow(M)){
vecj<-NULL
for(i in 1:nrow(matr)){
vec1ji<-2*(M[j,1]-matr[i,1])+2*(M[j,1]-matr[i,1])*log((M[j,1]-matr[i,1])^2+(M[j,2]-matr[i,2])^2)
vec2ji<-2*(M[j,2]-matr[i,2])+2*(M[j,2]-matr[i,2])*log((M[j,1]-matr[i,1])^2+(M[j,2]-matr[i,2])^2)
vecji<-c(vec1ji,vec2ji)
vecj<-rbind(vecj,vecji)
}
veclist<-c(veclist,list(vecj))
}
jac<-NULL
for(i in 1: length(veclist)){
jaci<-myl$at+t(myl$W)%*%veclist[[i]]
jac<-c(jac,list(as.matrix(jaci)))
}
deltus<-NULL
dpl<-NULL
sens<-NULL
for(i in 1:length(jac)){
deltusi<-jac[[i]]-diag(nrow(jac[[i]]))
dpli<-0.5*(deltusi+t(deltusi))
seni<-0.5*matrix.trace(t(as.matrix(dpli))%*%as.matrix(dpli))
deltus<-c(deltus,list(as.matrix(deltusi)))
dpl<-c(dpl,list(as.matrix(dpli)))
sens<-c(sens,seni)
}
sens2<-sens*areas
jac2<-list2array(jac)
mean11<-mean(jac2[1,1,])
mean12<-mean(jac2[1,2,])
mean21<-mean(jac2[2,1,])
mean22<-mean(jac2[2,2,])
detmean<-det(matrix(c(mean11,mean21,mean12,mean22),ncol=2))
myj<-unlist(lapply(jac,det))
if(body==T){
veclist2<-NULL
for(j in 1:nrow(M)){
vecj<-NULL
for(i in 1:nrow(matr)){
vec1ji<-2+(4*(M[j,1]-matr[i,1])^2/((M[j,1]-matr[i,1])^2+(M[j,2]-matr[i,2])^2))+(2*log((M[j,1]-matr[i,1])^2+(M[j,2]-matr[i,2])^2))
vec2ji<-2+(4*(M[j,2]-matr[i,2])^2/((M[j,1]-matr[i,1])^2+(M[j,2]-matr[i,2])^2))+(2*log((M[j,1]-matr[i,1])^2+(M[j,2]-matr[i,2])^2))
vec12ji<-4*(M[j,1]-matr[i,1])*(M[j,2]-matr[i,2])/((M[j,1]-matr[i,1])^2+(M[j,2]-matr[i,2])^2)
#
# vec1ji<-4*(M[j,1]-matr[i,1])^2/((M[j,1]-matr[i,1])^2+(M[j,2]-matr[i,2])^2)
# vec2ji<-4*(M[j,2]-matr[i,2])^2/((M[j,1]-matr[i,1])^2+(M[j,2]-matr[i,2])^2)
# vec12ji<-4*(M[j,1]-matr[i,1])*(M[j,2]-matr[i,2])/((M[j,1]-matr[i,1])^2+(M[j,2]-matr[i,2])^2)
arrji<-array(matrix(c(vec1ji,vec12ji,vec12ji,vec2ji),ncol=2),dim=c(2,2,1))
vecj<-abind::abind(vecj,arrji)
}
veclist2<-c(veclist2,list(vecj))
}
veclist2reshaped<-NULL
for(z in 1:length(veclist2)){
veclist2reshapedz<-NULL
for(k in 1:dim(veclist2[[z]])[[3]]){
veclist2reshapedzk<-c(veclist2[[z]][,,k])
veclist2reshapedz<-rbind(veclist2reshapedz,veclist2reshapedzk)
}
veclist2reshaped<-c(veclist2reshaped,list(veclist2reshapedz))
}
jac2<-NULL
for(i in 1: length(veclist2reshaped)){
jac2i<-t(myl$W)%*%veclist2reshaped[[i]]
jac2<-c(jac2,list(as.matrix(jac2i)))
}
beb<-NULL
for(i in 1:length(areas)){
bebi<-norm(jac2[[i]],type="F")^2*areas[[i]]
beb<-c(beb,bebi)
}
beb<-sum(beb)
}else{beb<-"No beb"}
out<-list(be=be,ge=ge,se=sum(sens2),area=sum(areas),ip=ip,beb=beb)
out
}
#' @export
matrix2arrayasis<-function(matrix,nland){
require(abind)
fin<-NULL
div<-newtypelist(nland,nrow(matrix)/nland)
for(i in 1:length(div)){
piece<-array(matrix[div[[i]],],dim=c(nland,ncol(matrix),1))
fin<-abind::abind(fin,piece)
}
return(fin)
}
#' export
newtypelist<-function(nlevelinteger,ncolinteger){####internal function
mynewlist<-NULL
myouter<-as.vector(outer(1:ncolinteger,nlevelinteger))
myouter2<-(myouter-nlevelinteger)+1
for(i in 1:length(myouter)){
myveci<-list(c(myouter2[i]:myouter[i]))
mynewlist<-append(mynewlist,myveci)
}
mynewlist
}
#' export
plot2dhull<-function(matrix,group,scalevariable=1,asp=1,extl=F,legend=T,xl=NULL,yl=NULL,posl=c("topright"),labels=NULL,clabel=0,pch=19,lwd=0.5,colhull=NULL,col=as.numeric(group),xlab=NULL,ylab=NULL,grey=F,xlimi=range(matrix[,1])*1.3,ylimi=range(matrix[,2])*1.3,reorder=T,alpha=rep(0,nlevels(group))){
library(MASS)
library(compositions)
library(spatstat)
library(gdata)
library(TeachingDemos)
library(ade4)
library(calibrate)
if (!is.factor(group)) stop("'group' must be a factor")
makeTransparent = function(..., alpha=0.5) {
if(alpha<0 | alpha>1) stop("alpha must be between 0 and 1")
alpha = floor(255*alpha)
newColor = col2rgb(col=unlist(list(...)), alpha=FALSE)
.makeTransparent = function(col, alpha) {
rgb(red=col[1], green=col[2], blue=col[3], alpha=alpha, maxColorValue=255)
}
newColor = apply(newColor, 2, .makeTransparent, alpha=alpha)
return(newColor)
}
x<-matrix[,1]
y<-matrix[,2]
if(reorder==T){group<-factor(group,levels=unique(group))}
species<-as.numeric(group)
col=col
if(grey==T){col=col2grey(col)}
coli<-data.frame(group,as.numeric(group),col,pch)
if(is.null(colhull)){colhull=aggregate(coli[,-1],by=list(coli[,1]),mean)[,-3][,2]}else{colhull=colhull}
colim<-cbind(aggregate(coli[,-1],by=list(coli[,1]),mean)[,-c(3,4)],colhull)
plot(x,y,type="p",cex=scalevariable,col=as.character(coli[,3]),pch=pch,xlim=xlimi,ylim=ylimi,xlab=xlab,ylab=ylab,asp=asp)
if(!is.null(labels)){textxy(x,y,labels)}else{NULL}
for(i in 1:(max(species)))
{
abline(h=0,v=0)
if(length(x[species==i])<3){NULL}else{
hulli<-convexhull.xy(subset(cbind(x,y),species==i))
par(new=T)
plot(hulli,col=makeTransparent(colim[,3][i], 1,alpha=alpha[i]),lwd=lwd,lty=i,add=T,asp=asp)
par(new=T)
daplotx<-c(hulli$bdry[[1]][[1]][length(hulli$bdry[[1]][[1]])],hulli$bdry[[1]][[1]])
daploty<-c(hulli$bdry[[1]][[2]][length(hulli$bdry[[1]][[2]])],hulli$bdry[[1]][[2]])
par(new=T)
plot(daplotx,daploty,lwd=lwd,type="l",lty=i,col=colim[,3][i],xlim=xlimi,ylim=ylimi,axes=F,xlab="",ylab="",asp=asp)
}
}
opar <- par(mar = par("mar"))
par(mar = c(0.1, 0.1, 0.1, 0.1))
on.exit(par(opar))
meansx<-aggregate(matrix[,1],by=list(group),FUN=mean)
meansy<-aggregate(matrix[,2],by=list(group),FUN=mean)
if (clabel > 0)
for(i in 1:nlevels(group)){ scatterutil.eti(meansx[i,-1],meansy[i,-1], meansx[,1][i], clabel,coul=1)}
if(legend==T){
if(extl==T){
x11()
plot(x,y,col="white")
legend(min(x),max(y),unique(coli[,1]), cex=1, col=unique(coli[,3]), pch=unique(coli[,4]),box.col="white")
}else{
if(is.null(xl)==T&is.null(yl)==T){
legend(posl,legend=unique(coli[,1]), col=unique(coli[,3]), pch=unique(coli[,4]),bty='n')}else{
legend(xl,yl,unique(coli[,1]),col=unique(coli[,3]),pch=unique(coli[,4]),bty='n')}
}
}
}
#' export
ptau6<-function(array,factor,CSinit=T,sepure=F,polyn=1,CR=NULL,locs=NULL,perm=999){
library(Morpho)
library(shapes)
library(vegan)
warning("WARNING: this function reorders data (if they are not) in increasing size order within each level")
k<-dim(array)[1]
m<-dim(array)[2]
n<-dim(array)[3]
newarray<-array[,,order(factor,centroid.size(array))]
factor<-factor[order(factor)]
cbind(dimnames(newarray)[[3]],as.character(factor))
depepure<-array2mat(procSym(newarray,pcAlign=F,CSinit=CSinit,scale=F)$orpdata,n,k*m)
if(sepure==T){depepure<-array2mat(sepgpa(newarray,factor,CSinit=CSinit,scale=F)$mountedorpdata,n,k*m)}
indepepure<-apply(newarray,3,cSize)
print("Individual multivariate (linear) regression between shape and size" )
print(manymultgr(depepure,indepepure,factor,steps=perm))
print("Pairwise *linear* mancova p-values on original data")
print(pwpermancova(depepure,indepepure,factor,nperm=perm)$p_adonis_pred1_pred2)
thedatapure<-data.frame(indepure=indepepure,depure=depepure)
thelmlistpure<-NULL
mypredictpure<-NULL
myresidpure<-NULL
for(i in 1:nlevels(factor)){
thelmpurei<-lm(as.matrix(thedatapure[,-1][as.numeric(factor)==i,])~poly(indepure[as.numeric(factor)==i],degree=polyn,raw=TRUE),data=thedatapure)
mypredictpurei<-predict(thelmpurei)
myresidpurei<-resid(thelmpurei)
thelmlistpure<-c(thelmlistpure,list(thelmpurei))
mypredictpure<-rbind(mypredictpure,mypredictpurei)
myresidpure<-rbind(myresidpure,myresidpurei)
}
mypredictpure<-read.inn(mypredictpure,k,m)
myresidpure<-read.inn(myresidpure,k,m)
if(is.null(CR)==T){CR<-procSym(mypredictpure[,,firstsfac(factor)],scale=F,pcAlign=F,reflect=F,CSinit=F)$mshape}else{CR<-CR}
if(is.null(locs)==T){locs<-mypredictpure[,,firstsfac(factor)]}else{locs<-locs}
prls<-lshift2(mypredictpure,factor,CSinit=CSinit,CR=CR,locs=locs)
common<-array2mat(procSym(newarray,pcAlign=,scale=F,CSinit=CSinit)$orpdata,n,k*m)
print("Pairwise multivariate (linear) regression between shape and size" )
print(pwpermancova(common,indepepure,factor,nperm=perm)$p_adonis_pred1_pred2)
space1<-prcomp(array2mat(prls$transported,n,k*m))
space1mshape<-procSym(prls$transported,pcAlign=F)$mshape
origtrasp<-prls$transported+myresidpure
origproj<-predict(space1,array2mat(origtrasp,n,k*m))
print("Pairwise multivariate (linear) regression between shape of transported data and size" )
print(pwpermancova(array2mat(origtrasp,n,k*m),indepepure,factor,nperm=perm)$p_adonis_pred1_pred2)
depepure2<-origtrasp
print("Individual multivariate (linear) regression between shape of transported data and size" )
print(manymultgr(array2mat(depepure2,n,k*m),indepepure,factor,steps=perm))
thedatapure2<-data.frame(indepure2=indepepure,depure2=array2mat(depepure2,n,k*m))
thelmlistpure2<-NULL
mypredictpure2<-NULL
myresidpure2<-NULL
eqpredtrasp<-NULL
eqsizes<-NULL
for(i in 1:nlevels(factor)){
thelmpure2i<-lm(as.matrix(thedatapure2[,-1][as.numeric(factor)==i,])~poly(indepure2[as.numeric(factor)==i],degree=polyn,raw=TRUE),data=thedatapure2)
mypredictpure2i<-predict(thelmpure2i)
myresidpure2i<-resid(thelmpure2i)
thelmlistpure2<-c(thelmlistpure2,list(thelmpure2i))
mypredictpure2<-rbind(mypredictpure2,mypredictpure2i)
myresidpure2<-rbind(myresidpure2,myresidpure2i)
eqpredtraspi<-serpred(array2mat(depepure2,n,k*m)[as.numeric(factor)==i,],indepepure[as.numeric(factor)==i],polyn=polyn)
eqpredtrasp<-rbind(eqpredtrasp,eqpredtraspi)
eqsizesi<-seq(min(indepepure[as.numeric(factor)==i]),max(indepepure[as.numeric(factor)==i]),length.out=10)
eqsizes<-c(eqsizes,eqsizesi)
}
eqpredtrasp<-read.inn(eqpredtrasp,k,m)
out<-list(k=k,m=m,n=n,arrayord=newarray,factorord=factor,CR=CR,locs=locs,depepure=depepure,indepepure=indepepure,thelmlistpure=thelmlistpure,predictpure=mypredictpure,residpure=myresidpure,shifted=array2mat(prls$transported,n,k*m),thelmlistpure2=thelmlistpure2,predictpure2=array2mat(prls$transported,n,k*m),predictpure3=mypredictpure2,residpure2=myresidpure2,space1=space1,origtrasp=array2mat(origtrasp,n,k*m),origproj=origproj,space1mshape=space1mshape,eqpredtrasp=eqpredtrasp,eqsizes=eqsizes,sizerange=aggregate(indepepure,by=list(factor),range))
}
#' export
plotptau5<-function(objptau,links=NULL,projorig=T,whichaxes=c(1,2),cexscale=round(max(centroid.size(objptau$arrayord)),digits=0),shapescale=10,mag=1,subplotdim=2,shiftnegy=1,shiftposy=1,shiftnegx=1,shiftposx=1,col=as.numeric(objptau$factorord),pch=as.numeric(objptau$factorord),triang=NULL,heat=T,from=NULL,scaleramp=F,topcs=NULL,tore=NULL,plotsource=T){
library(TeachingDemos)
k<-dim(objptau$arrayord)[1]
m<-dim(objptau$arrayord)[2]
n<-dim(objptau$arrayord)[3]
transp<-read.inn(objptau$shifted,k,m)
if(m>2){
if(heat==T&is.null(triang)==T){stop("with heat=T in 3D I need triangulation")}
}
origsize<-objptau$indepepure
if(projorig==T){
ordiwithshapes(objptau$space1mshape,objptau$space1$x,objptau$space1$rotation,procSym=F,whichaxes=whichaxes,addata=objptau$origproj[,whichaxes],asp=1,factraj=objptau$factorord,cex=origsize/cexscale,triang=triang,from=from,to=topcs,links=links,mag=mag,shiftnegy=1.5,shiftposy=2,col=col,pch=pch,subplotdim=subplotdim,plotsource=plotsource)
title("LS on predictions in PCA space and original data re-projected")}else{
ordiwithshapes(objptau$space1mshape,objptau$space1$x,objptau$space1$rotation,procSym=F,whichaxes=whichaxes,asp=1,factraj=objptau$factorord,cex=origsize/cexscale,triang=triang,from=from,to=topcs,links=links,mag=mag,shiftnegy=1.5,shiftposy=2,col=col,pch=pch,subplotdim=subplotdim,plotsource=plotsource)
title("LS on predictions in PCA space")}
ranges<-apply(rbind(objptau$space1$x,objptau$origproj),2,range)
ratesorig<-ratesbygroup(objptau$predictpure,objptau$factorord,objptau$indepepure,plot=F)
ratestrasp<-ratesbygroup(read.inn(objptau$predictpure2,k,m),objptau$factorord,objptau$indepepure,plot=F)
plot(objptau$indepepure,ratesorig,pch=pch,col=col)
for(i in 1:nlevels(objptau$factorord)){
lines(objptau$indepepure[as.numeric(objptau$factorord)==i][-firstsfac(objptau$factorord)],ratesorig[as.numeric(objptau$factorord)==i][-firstsfac(objptau$factorord)],col=col[firstsfac(objptau$factorord)][i])
}
title("Rates of shape change among consecutive predictions per unit size in original data")
if(m>2){
if(is.null(triang)){open3d()}
if(is.null(triang)){heat=F}
rgl.close()
if(is.null(triang)){triang2<-NULL}else{triang2<-t(triang)}
mydef<-plotdefo3d(list(mshape=objptau$space1mshape,PCs=objptau$space1$rotation,PCscores=objptau$space1$x),procSym=F,triang=triang2,from=from,to=topcs,links=links,heat=heat,out=T)
}else{
mydef<-plotdefo(list(mshape=objptau$space1mshape,PCs=objptau$space1$rotation,PCscores=objptau$space1$x),procSym=F,linkss=links,linkst=links,heat=heat)
}
reset.par()
plot(objptau$indepepure,ratestrasp,pch=pch,col=col)
for(i in 1:nlevels(objptau$factorord)){
lines(objptau$indepepure[as.numeric(objptau$factorord)==i][-firstsfac(objptau$factorord)],ratestrasp[as.numeric(objptau$factorord)==i][-firstsfac(objptau$factorord)],col=col[firstsfac(objptau$factorord)][i])
}
title("Rates of shape change among consecutive predictions per unit size in transported data")
plot(ratesorig,ratestrasp,col=col,pch=pch)
title("Original rates of shape change per unit size vs those of transported data")
adults<-NULL
adultsdaplot<-NULL
for(i in 1:nlevels(objptau$factorord)){
adulti<-array(showPC2(objptau$space1$x[lastsfac(objptau$factorord),1:3][i,],objptau$space1$rotation[,1:3],objptau$space1mshape),dim=c(k,m,1))
adultdaploti<-array(showPC2(objptau$space1$x[lastsfac(objptau$factorord),1:3][i,]*mag,objptau$space1$rotation[,1:3],objptau$space1mshape),dim=c(k,m,1))
adults<-abind::abind(adults,adulti)
adultsdaplot<-abind::abind(adultsdaplot,adultdaploti)
}
if(m<3){adultsdaplot<-abind::abind(adultsdaplot,array(rep(0,k*n),dim=c(k,1,nlevels(objptau$factorord))),along=2)}
init<-array(showPC2(objptau$space1$x[firstsfac(objptau$factorord),1:3][1,],objptau$space1$rotation[,1:3],objptau$space1mshape),dim=c(k,m,1))
initdaplot<-array(showPC2(objptau$space1$x[firstsfac(objptau$factorord),1:3][1,]*mag,objptau$space1$rotation[,1:3],objptau$space1mshape),dim=c(k,m,1))
if(m<3){initdaplot<-abind::abind(initdaplot,array(rep(0,k*n),dim=c(k,1,1)),along=2)}
if(m>2){
if(heat==T){
adultsdaplotvis<-NULL
alldists<-NULL
for(i in 1:nlevels(objptau$factorord)){
adultsdaplotvisi<-heat3d(initdaplot[,,1],adultsdaplot[,,i],t(triang),graphics=F,to=tore,from=from,scaleramp=scaleramp)
alldists<-c(alldists,list(adultsdaplotvisi$obs))
daagg<-objptau$space1$x[1:n,1:3][lastsfac(objptau$factorord),][i,]
adultsdaplotvisi<-translate3d(scalemesh(adultsdaplotvisi$obm$colMesh,1/shapescale,center="none"),daagg[1],daagg[2],daagg[3])
adultsdaplotvis<-c(adultsdaplotvis,list(adultsdaplotvisi))
}
}else{adultsdaplotvis<-NULL;alldists<-NULL}
}
open3d()
r3dDefaults$windowRect <- c(0,50, 8000, 8000)
mfrow3d(1,2,byrow=T,sharedMouse=T)
for(i in 1:nlevels(objptau$factorord)){
lines3d(objptau$space1$x[1:n,1:3][as.numeric(objptau$factorord)==i,],col=i,add=T)
}
title3d("LS data in the PCA; first three PC scores")
decorate3d(xlab="PC1",ylab="PC2",zlab="PC3")
if(projorig==T){
plot3D(objptau$origproj[,1:3],col=as.numeric(objptau$factorord),add=T,bbox=F,cex=objptau$indepepure/cexscale)
}
if(heat==F){
babedaagg<-rep.row(objptau$space1$x[firstsfac(objptau$factorord),1:3][1,],k)
for(i in 1:nlevels(objptau$factorord)){
daplottarei<-adultsdaplot[,,i]
daagg<-rep.row(objptau$space1$x[1:n,1:3][lastsfac(objptau$factorord),][i,],k)
plot3D((daplottarei/shapescale)+daagg,col=i,add=T,size=1,bbox=F)
if(!is.null(links)){lineplot((daplottarei/shapescale)+daagg,links,col=i)}
if(!is.null(triang)){plotsurf((daplottarei/shapescale)+daagg,triang,col=i)}
}
plot3D((initdaplot[,,1]/shapescale)+babedaagg,bbox=F,add=T)
if(!is.null(links)){lineplot((initdaplot[,,1]/shapescale)+babedaagg,links)}
if(!is.null(triang)){plotsurf((initdaplot[,,1]/shapescale)+babedaagg,triang,col=1)}
}else{
if(m>2){
for(i in 1:nlevels(objptau$factorord)){
meshDist(adultsdaplotvis[[i]],distvec=alldists[[i]],from=range(alldists)[1],to=range(alldists)[2],add=T,scaleramp=scaleramp)
}
babedaagg<-rep.row(objptau$space1$x[firstsfac(objptau$factorord),1:3][1,],k)
shade3d(plotsurf((initdaplot[,,1]/shapescale)+babedaagg,triang,plot=F),add=T,alpha=0.5)
}
}
next3d()
babedaagg<-objptau$space1$x[firstsfac(objptau$factorord),1:3][1,]
for(i in 1:nlevels(objptau$factorord)){
daagg<-rep.row(objptau$space1$x[1:n,1:3][lastsfac(objptau$factorord),][i,],k)
text3d(daagg,texts=levels(objptau$factorord)[i],col=i,add=T,size=1)
}
text3d(babedaagg,texts="CR",bbox=F,add=T)
decorate3d(xlab="PC1",ylab="PC2",zlab="PC3")
for(i in 1:nlevels(objptau$factorord)){
lines3d(objptau$space1$x[1:n,1:3][as.numeric(objptau$factorord)==i,],col=i,add=T)
}
title3d("LS data in the PCA; first three PC scores")
out<-list(heat=heat,ratesorig=ratesorig,ratestrasp=ratestrasp,alldists=alldists,init=init,initdaplot=initdaplot,adults=adults,adultsdaplot=adultsdaplot,adultsdaplotvis=adultsdaplotvis,links=links,triang=triang,plotdefout=mydef)
}
#' export
benfrombygroup<-function(locs,array,factor,doopa=T){
resu<-NULL
for(j in 1:nlevels(factor)){
for(k in 1:table(factor)[j]){
resui<-jbe(locs[,,j],array[,,as.numeric(factor)==j][,,k],doopa=doopa)
resu<-append(resu,resui)
print(paste(j,k,sep="_"))
}
}
resu
}
#' @export
plotancova<-function(y,x,group=NULL,polyn=1,pch=NULL,col=1,shade=T,alpha=0.1,confint=T,cex=0.5,asp=NULL,legend=T,labels=NULL,plot=T,xlab=NULL,ylab=NULL,xlim=range(x),ylim=range(y)){
library(phia)
if(is.null(group)==T){group<-factor(c(rep(1,length(y))))}else{pch=pch}
if(is.null(pch)==T){pch=19}else{pch=pch}
dati<-data.frame(group,as.numeric(group),col,pch)
if(is.null(xlab)==T){xlab<-c("x")}else{xlab<-xlab}
if(is.null(ylab)==T){ylab<-c("y")}else{ylab<-ylab}
if(plot==T){
par(mar=c(par('mar')[1:3], 0)) # optional, removes extraneous right inner margin space
plot.new()
l <- legend(0, 0, bty='n', xpd=NA,unique(group), pch=unique(pch), col=unique(col),plot=FALSE)
# calculate right margin width in ndc
w <- grconvertX(l$rect$w, to='ndc') - grconvertX(0, to='ndc')
par(omd=c(0, 1-w, 0, 1))
plot(x,y,xlim=xlim,ylim=ylim,col=col,pch=pch,cex=cex,xlab=xlab,ylab=ylab,asp=asp)
if(legend==T){legend(par('usr')[2], par('usr')[4], bty='n', xpd=NA,unique(group), pch=unique(pch), col=unique(col))}
if(!is.null(labels)){textxy(x,y,labels)}
}
lmlist<-NULL
for(i in 1:nlevels(group)){
localy<-y[as.numeric(group)==i]
localx<-x[as.numeric(group)==i]
lmi<-lm(localy~poly(localx,degree=polyn,raw=TRUE))
lmlist<-c(lmlist,list(lmi))
}
names(lmlist)<-levels(group)
datapredinflist<-NULL
for(i in 1:nlevels(group)){
datapredinfi<-cbind(x[as.numeric(group)==i], predict(lmlist[[i]], interval="confidence")[,2])
datapredinfi<-datapredinfi[order(datapredinfi[,2]),]
datapredinflist<-c(datapredinflist,list(datapredinfi))
}
names(datapredinflist)<-levels(group)
datapredsuplist<-NULL
for(i in 1:nlevels(group)){
datapredsupi<-cbind(x[as.numeric(group)==i], predict(lmlist[[i]], interval="confidence")[,3])
datapredsupi<-datapredsupi[order(datapredsupi[,2]),]
datapredsuplist<-c(datapredsuplist,list(datapredsupi))
}
names(datapredsuplist)<-levels(group)
if(plot==T){
if(polyn<2){
for(i in 1: length(lmlist)){
abline(lmlist[[i]],col=dati[as.numeric(group)==i,][1,3],ylim=ylim,xlim=xlim)
} }else{
for(i in 1: length(lmlist)){
newx <- seq(min(x[as.numeric(group)==i]), max(x[as.numeric(group)==i]), length.out=500)
preds <- predict(lmlist[[i]],newdata = data.frame(localx=newx))
matord<-cbind(newx,preds)
matord<-matord[order(matord[,1]),]
lines(matord,cex=0.1,col=dati[as.numeric(group)==i,][1,3],ylim=ylim,xlim=xlim)
}
}
}
if(plot==T){
if(confint==T){
for(i in 1: length(lmlist)){
if(anova(lmlist[[i]])$Pr[1]<0.05){
matord<-datapredinflist[[i]][order(datapredinflist[[i]][,1]),]
lines(matord,cex=0.1,lty = 'dashed',col=dati[as.numeric(group)==i,][1,3],ylim=ylim,xlim=xlim)
}}
for(i in 1: length(lmlist)){
if(anova(lmlist[[i]])$Pr[1]<0.05){
matord<-datapredsuplist[[i]][order(datapredsuplist[[i]][,1]),]
lines(matord,cex=0.1,lty = 'dashed',col=dati[as.numeric(group)==i,][1,3],ylim=ylim,xlim=xlim)
}}}
if(shade==T){
for(i in 1: length(lmlist)){
if(anova(lmlist[[i]])$Pr[1]<0.05){
newx <- seq(min(x[as.numeric(group)==i]), max(x[as.numeric(group)==i]), length.out=1000)
preds <- predict(lmlist[[i]],newdata = data.frame(localx=newx), interval = 'confidence')
polygon(c(rev(newx), newx), c(rev(preds[ ,3]), preds[ ,2]), col = makeTransparent(dati[as.numeric(group)==i,][1,3],alpha=alpha), border = NA)
}
}
}
}
summarylist<-NULL
for(i in 1:length(lmlist)){
summarylisti<-summary(lmlist[[i]])
summarylist<-c(summarylist,list(summarylisti))
}
names(summarylist)<-levels(group)
if(polyn<2){
sint_results<-NULL
for(i in 1:length(summarylist)){
inti<-summarylist[[i]]$coefficients[1,1]
p_inti<-summarylist[[i]]$coefficients[1,4]
betai<-summarylist[[i]]$coefficients[2,1]
p_betai<-summarylist[[i]]$coefficients[2,4]
r_sqi<-summarylist[[i]]$r.squared
sinti<-c(inti,p_inti,betai,p_betai,r_sqi)
sint_results<-rbind(sint_results,sinti)
}
rownames(sint_results)<-levels(group)
colnames(sint_results)<-c("Intercept","p-value Intercept","Beta","p-value Beta","R squared")
print(sint_results)}else{sint_results<-NULL}
if(polyn<2){
if(nlevels(group)>1){
df<-data.frame(y=y,x=x,group=group)
model<-lm(y~x*group)
slopint<-testInteractions(model, pairwise="group", slope="x")
slopintmat<- matrix(NA, nlevels(group), nlevels(group))
slopintmat[lower.tri(slopintmat) ] <-round( slopint[1:(nrow(slopint)-1),5], 5)
slopintmat<-t(slopintmat)
colnames(slopintmat)<-levels(group)
rownames(slopintmat)<-levels(group)
slopintmat[lower.tri(slopintmat)]<-slopintmat[upper.tri(slopintmat)]
elevint<-testInteractions(lm(y~x+group))
elevintmat<- matrix(NA, nlevels(group), nlevels(group))
elevintmat[lower.tri(elevintmat) ] <-round( elevint[1:(nrow(elevint)-1),5], 5)
elevintmat<-t(elevintmat)
colnames(elevintmat)<-levels(group)
rownames(elevintmat)<-levels(group)
elevintmat[lower.tri(elevintmat)]<-elevintmat[upper.tri(elevintmat)]
}else{
slopintmat<-NULL
elevintmat<-NULL
}}else{
slopintmat<-NULL
elevintmat<-NULL
}
out<-list(lmlist=lmlist,datapredsuplist=datapredsuplist,datapredinflist=datapredinflist,summarylist=summarylist,sint_results=sint_results,slopepval=slopintmat,elevpval=elevintmat)
return(out)
}
#' @export
heat3d<-function(source,target,triang,iter=0,dointerp=F,linkss=NULL,linkst=NULL,plotlands=F,legend=T,cols=1,colt=2,plotsource=T,plottarget=T,collinkss=1,lwds=2,collinkst=2,lwdt=2,cexs=0.5,cext=0.5,colors=c("blue","cyan","yellow","red"),alpha=1,ngrid=0,mag=1,graphics=T,to=NULL,from=NULL,scaleramp=F,lines=T){
mate<-source
mate2<-target
mate2<-mate+(mate2-mate)*mag
library(fields)
mmate<-plotsurf(mate,t(triang),plot=F)
mmate2<-plotsurf(mate2,t(triang),plot=F)
mesh2listri<-function(mat,tri){
if(ncol(tri)>3){tri<-t(tri)}
res <- apply(tri, 1, function(x) mat[x,])
res <- lapply(as.data.frame(res), function(x) matrix(x, nrow = 3, ncol = ncol(tri)))
res
}
ar01<-vcgArea(mmate,perface=T)$pertriangle
ar02<-vcgArea(mmate2,perface=T)$pertriangle
obs0<-log(ar02/ar01)
M0<-vcgBary(mmate)
M02<-vcgBary(mmate2)
tes<-tessell3d(triang,mate,iter)
class(tes) <- "mesh3d"
matr<-mate
M<-vcgBary(tes)
tes2<-tessell3d(triang,mate2,iter)
class(tes2) <- "mesh3d"
M2<-vcgBary(tes2)
triangtes<-t(tes$it)
if(iter>0){
ar1<-vcgArea(tes,perface=T)$pertriangle
ar2<-vcgArea(tes2,perface=T)$pertriangle
}
if(iter>0){obsref<-log(ar2/ar1)}else{obsref<-obs0}
if(dointerp==T){fit<- fastTps(M02,obs0,theta=0.009)}else{fit<-NULL}
if(dointerp==T){obs2<-predict(fit,rbind(mate2,M2))}else{obs2<-NULL}
if(dointerp==T){obs3<-c(obs2)}else{obs3<-NULL}
if(dointerp==T){obs3[is.na(obs3)]<-mean(obs3,na.rm=T)}
if(dointerp==T){obm<-meshDist(plotsurf(rbind(mate2,M2),mmate$it,plot=F),distv=obs3,add=T,rampcolors =colors,to=to,from=from,scaleramp=scaleramp,plot=F,alpha=alpha)}else{obm<-NULL}
if(graphics==T){
if(plotlands==T){deformGrid3d(mate,mate2,ngrid=ngrid,lines=lines,col1=cols,col2=colt)}
if(dointerp==T){shade3d(obm$colMesh,alpha=alpha)}
}
if(dointerp==F){
if(is.null(from)&is.null(to)==T){
colo<-plotrix::color.scale(obs0,c(0,0,1,1),c(0,0.93,1,0),c(0.78,0.93,0,0),color.spec="rgb")
if(graphics==T){shade3d(mmate2,col=rep(colo,each=3))}
}else{
colo1<-plotrix::color.scale(c(seq(from,to,length.out=100),obs0),c(0,0,1,1),c(0,0.93,1,0),c(0.78,0.93,0,0),color.spec="rgb")
colo<-colo1[-c(1:100)]
if(graphics==T){shade3d(mmate2,col=rep(colo,each=3),alpha=alpha)}
}
obm<-list(colMesh=list(vb=mmate$vb,it=mmate$it,material=list(color=as.matrix(rep(colo,each=3)))))
class(obm$colMesh)<-"mesh3d"
}
if(is.null(from)&is.null(to)==T){
colo1<-plotrix::color.scale(seq(min(obs0),max(obs0),length.out=100),c(0,0,1,1),c(0,0.93,1,0),c(0.78,0.93,0,0),color.spec="rgb")
if(graphics==T){
plot(cbind(seq(min(obs0),max(obs0),length.out=100),0),col=colo1,cex=10,pch=15,xaxt="n",yaxt="n",axes=F,xlab="",ylab="")
axis(1)
}
}else{
colo1<-plotrix::color.scale(c(seq(from,to,length.out=100),obs0),c(0,0,1,1),c(0,0.93,1,0),c(0.78,0.93,0,0),color.spec="rgb")
if(graphics==T){
plot(cbind(c(seq(from,to,length.out=100)),0),col=colo1[1:100],cex=10,pch=15,xaxt="n",yaxt="n",axes=F,xlab="",ylab="")
axis(1)
}
}
if(graphics==T){
if(!is.null(linkst)){lineplot(mate2,linkst,col=collinkst,lwd=lwds)}
if(plotsource==T){
shade3d(plotsurf(source,t(triang),plot=F),alpha=0.5)
if(!is.null(linkss)){lineplot(mate,linkss,col=collinkss,lwd=lwdt)}}
}
if(iter>0){ar1<-ar1}else{ar1<-c("no iter")}
if(iter>0){ar2<-ar2}else{ar2<-c("no iter")}
out<-list(mate=mate,mate2=mate2,mmate=mmate,mmate2=mmate2,ar01=ar01,ar02=ar02,ar1=ar1,ar2=ar2,M0=M0,M02=M02,M=M,M2=M2,obsref=obsref,obs2=obs2,obs=obs0,obs3=obs3,fit=fit,cols=makeTransparent(colorRampPalette(colors)(n = length(obs3)),alpha=alpha),tes=tes,tes2=tes2,obm=obm)
out
}
#' @export
plotsurf<-function(lmatrix,triang,col=1,alpha=0.5,plot=T){
thesurf_1=t(cbind(lmatrix,1))
triangolazioni=triang
thesurf=list(vb=thesurf_1,it=triangolazioni)
class(thesurf) <- "mesh3d"
if(plot==T){shade3d(thesurf,col=col,alpha=alpha)}
return(thesurf)
}
#' @export
tessell3d=function(tri,ver,iter){
tri_t=tri
ver_t=ver
for(j in 1:iter){
nrows=nrow(tri_t)
numb=range(tri_t)[2]
for(i in 1:nrows){
tri_i=ver_t[tri_t[i,],]
cen_i=colMeans(tri_i)
tri_1=c(tri_t[i,c(1,2)],numb+i)
tri_2=c(tri_t[i,1],numb+i,tri_t[i,3])
tri_3=c(tri_t[i,c(2,3)],numb+i)
tri_t=rbind(tri_t,tri_1,tri_2,tri_3)
ver_t=rbind(ver_t,cen_i)
}}
vertici_out=cbind(ver_t,1)
rownames(vertici_out)=NULL
triangoli_out=tri_t[(dim(tri)[1]+1):dim(tri_t)[1],]
rownames(triangoli_out)=NULL
mesh.out=list(vb=t(vertici_out),it=t(triangoli_out))
}
#' @export
plotdefo3d<-function(procsymobject,procSym=T,triang=NULL,links=NULL,collinks=1,lwd=1,axtit=NULL,mags=rep(1,3),heat=F,scaleramp=F,heateuc=F,sign=T,colors=c("blue4","cyan2","yellow","red4"),alpha=1,from=NULL,to=NULL,out=F,plot=T){
if(is.null(axtit)==T){axtit=c("PC1+","PC2+","PC3+","PC1-","PC2-","PC3-")}else{axtit=axtit}
texts<-axtit
if(!is.null(triang)==T&&ncol(triang)>3){stop("I want triangles matrix in the form nx3")}
mshape<-procsymobject$mshape
if(procSym==T){
pc1pos<-showPC(max(procsymobject$PCscores[,1])*mags[1],procsymobject$PCs[,1],procsymobject$mshape)
pc1neg<-showPC(min(procsymobject$PCscores[,1])*mags[1],procsymobject$PCs[,1],procsymobject$mshape)
pc2pos<-showPC(max(procsymobject$PCscores[,2])*mags[2],procsymobject$PCs[,2],procsymobject$mshape)
pc2neg<-showPC(min(procsymobject$PCscores[,2])*mags[2],procsymobject$PCs[,2],procsymobject$mshape)
pc3pos<-showPC(max(procsymobject$PCscores[,3])*mags[3],procsymobject$PCs[,3],procsymobject$mshape)
pc3neg<-showPC(min(procsymobject$PCscores[,3])*mags[3],procsymobject$PCs[,3],procsymobject$mshape)
}else{
pc1pos<-showPC2(max(procsymobject$PCscores[,1])*mags[1],procsymobject$PCs[,1],procsymobject$mshape)
pc1neg<-showPC2(min(procsymobject$PCscores[,1])*mags[1],procsymobject$PCs[,1],procsymobject$mshape)
pc2pos<-showPC2(max(procsymobject$PCscores[,2])*mags[2],procsymobject$PCs[,2],procsymobject$mshape)
pc2neg<-showPC2(min(procsymobject$PCscores[,2])*mags[2],procsymobject$PCs[,2],procsymobject$mshape)
pc3pos<-showPC2(max(procsymobject$PCscores[,3])*mags[3],procsymobject$PCs[,3],procsymobject$mshape)
pc3neg<-showPC2(min(procsymobject$PCscores[,3])*mags[3],procsymobject$PCs[,3],procsymobject$mshape)
}
if(heat==F&heateuc==F){
allshapes<-rbind(pc1pos,pc1neg,pc2pos,pc2neg,pc3pos,pc3neg)
allshapes<-matrix2arrayasis(allshapes,nrow(pc1pos))
css<-apply(allshapes,3,cSize)
themax<-allshapes[,,which.max(css)]
if(plot==T){
open3d(windowRect=c(100,100,1000,1000))
mat <- matrix(1:12, ncol=2)
layout3d(mat, height = rep(c(3,1), 3), model = "inherit")
plot3d(themax*1.2,box=F,axes=F,col="white",xlab="",ylab="",zlab="",type="s",size=0,aspect=F)
plot3d(pc1pos,bbox=F,type="s",asp=F,axes=F,box=F,size=0.6,xlab = "", ylab = "", zlab = "",add=T)
if(is.null(links)==F){lineplot(pc1pos,links,col=collinks,lwd=lwd)}
if(is.null(triang)==F){plotsurf(pc1pos,t(triang))}
next3d()
text3d(0,0,0, texts[1])
next3d()
plot3d(themax*1.2,box=F,axes=F,col="white",xlab="",ylab="",zlab="",type="s",size=0,aspect=F)
plot3d(pc2pos,bbox=F,type="s",asp=F,axes=F,box=F,size=0.6,xlab = "", ylab = "", zlab = "",add=T)
if(is.null(links)==F){lineplot(pc2pos,links,col=collinks,lwd=lwd)}
if(is.null(triang)==F){plotsurf(pc2pos,t(triang))}
next3d()
text3d(0,0,0, texts[2])
next3d()
plot3d(themax*1.2,box=F,axes=F,col="white",xlab="",ylab="",zlab="",type="s",size=0,aspect=F)
plot3d(pc3pos,bbox=F,type="s",asp=F,axes=F,box=F,size=0.6,xlab = "", ylab = "", zlab = "",add=T)
if(is.null(links)==F){lineplot(pc3pos,links,col=collinks,lwd=lwd)}
if(is.null(triang)==F){plotsurf(pc3pos,t(triang))}
next3d()
text3d(0,0,0, texts[3])
next3d()
plot3d(themax*1.2,box=F,axes=F,col="white",xlab="",ylab="",zlab="",type="s",size=0,aspect=F)
plot3d(pc1neg,bbox=F,type="s",asp=F,axes=F,box=F,size=0.6,xlab = "", ylab = "", zlab = "",add=T)
if(is.null(links)==F){lineplot(pc1neg,links,col=collinks,lwd=lwd)}
if(is.null(triang)==F){plotsurf(pc1neg,t(triang))}
next3d()
text3d(0,0,0, texts[4])
next3d()
plot3d(themax*1.2,box=F,axes=F,col="white",xlab="",ylab="",zlab="",type="s",size=0,aspect=F)
plot3d(pc2neg,bbox=F,type="s",asp=F,axes=F,box=F,size=0.6,xlab = "", ylab = "", zlab = "",add=T)
if(is.null(links)==F){lineplot(pc2neg,links,col=collinks,lwd=lwd)}
if(is.null(triang)==F){plotsurf(pc2neg,t(triang))}
next3d()
text3d(0,0,0, texts[5])
next3d()
plot3d(themax*1.2,box=F,axes=F,col="white",xlab="",ylab="",zlab="",type="s",size=0,aspect=F)
plot3d(pc3neg,bbox=F,type="s",asp=F,axes=F,box=F,size=0.6,xlab = "", ylab = "", zlab = "",add=T)
if(is.null(links)==F){lineplot(pc3neg,links,col=collinks,lwd=lwd)}
if(is.null(triang)==F){plotsurf(pc3neg,t(triang))}
next3d()
text3d(0,0,0, texts[6])}
}
if (heat==T&heateuc==F){
myh1pos<-heat3d(mshape,pc1pos,triang,lines=F,linkss=links,graphics=F,from=from,to=to,colors=colors,alpha=alpha,scaleramp=scaleramp)
myh1neg<-heat3d(mshape,pc1neg,triang,lines=F,linkss=links,graphics=F,from=from,to=to,colors=colors,alpha=alpha,scaleramp=scaleramp)
myh2pos<-heat3d(mshape,pc2pos,triang,lines=F,linkss=links,graphics=F,from=from,to=to,colors=colors,alpha=alpha,scaleramp=scaleramp)
myh2neg<-heat3d(mshape,pc2neg,triang,lines=F,linkss=links,graphics=F,from=from,to=to,colors=colors,alpha=alpha,scaleramp=scaleramp)
myh3pos<-heat3d(mshape,pc3pos,triang,lines=F,linkss=links,graphics=F,from=from,to=to,colors=colors,alpha=alpha,scaleramp=scaleramp)
myh3neg<-heat3d(mshape,pc3neg,triang,lines=F,linkss=links,graphics=F,from=from,to=to,colors=colors,alpha=alpha,scaleramp=scaleramp)
distrange<-range(c(myh1pos$obs3,myh1neg$obs3,myh2pos$obs3,myh2neg$obs3,myh3pos$obs3,myh3neg$obs3))
}
if (heat==T&heateuc==T){
myh1pos<-diffonmesh(mshape,pc1pos,t(triang),graph=F,from=from,to=to,rampcolors=colors,alphas=c(alpha,0.7),sign=sign)
myh1neg<-diffonmesh(mshape,pc1neg,t(triang),graph=F,from=from,to=to,rampcolors=colors,alphas=c(alpha,0.7),sign=sign)
myh2pos<-diffonmesh(mshape,pc2pos,t(triang),graph=F,from=from,to=to,rampcolors=colors,alphas=c(alpha,0.7),sign=sign)
myh2neg<-diffonmesh(mshape,pc2neg,t(triang),graph=F,from=from,to=to,rampcolors=colors,alphas=c(alpha,0.7),sign=sign)
myh3pos<-diffonmesh(mshape,pc3pos,t(triang),graph=F,from=from,to=to,rampcolors=colors,alphas=c(alpha,0.7),sign=sign)
myh3neg<-diffonmesh(mshape,pc3neg,t(triang),graph=F,from=from,to=to,rampcolors=colors,alphas=c(alpha,0.7),sign=sign)
distrange<-range(c(myh1pos$obm$dists,myh1neg$obm$dists,myh2pos$obm$dists,myh2neg$obm$dists,myh3pos$obm$dists,myh3neg$obm$dists))
}
if(heat==T){
allshapes<-rbind(t(myh1pos$obm$colMesh$vb[-4,]),t(myh1neg$obm$colMesh$vb[-4,]),t(myh2pos$obm$colMesh$vb[-4,]),t(myh2neg$obm$colMesh$vb[-4,]),t(myh3pos$obm$colMesh$vb[-4,]),t(myh3neg$obm$colMesh$vb[-4,]))
allshapes<-matrix2arrayasis(allshapes,nrow(t(myh1pos$obm$colMesh$vb[-4,])))
css<-apply(allshapes,3,cSize)
themax<-allshapes[,,which.max(css)]}
if(heat==T){
if(plot==T){
open3d(windowRect=c(100,100,1000,1000))
mat <- matrix(1:12, ncol=2)
layout3d(mat, height = rep(c(3,1), 3), model = "inherit")
plot3d(themax*1.2,box=F,axes=F,col="white",xlab="",ylab="",zlab="",type="s",size=0,aspect=F)
if(heateuc==F){meshDist(myh1pos$obm$colMesh,alpha=alpha,distvec=myh1pos$obs3,from=distrange[1],to=distrange[2],scaleramp=scaleramp,add=T)}else{
meshDist(myh1pos$obm$colMesh,alpha=alpha,distvec=myh1pos$obm$dists,from=distrange[1],to=distrange[2],scaleramp=scaleramp,add=T)
}
if(is.null(links)==F){lineplot(t(myh1pos$obm$colMesh$vb[-4,]),links,col=collinks,lwd=lwd)}
next3d()
text3d(0,0,0, texts[1])
next3d()
plot3d(themax*1.2,box=F,axes=F,col="white",xlab="",ylab="",zlab="",type="s",size=0,aspect=F)
if(heateuc==F){meshDist(myh2pos$obm$colMesh,alpha=alpha,distvec=myh2pos$obs3,from=distrange[1],to=distrange[2],scaleramp=scaleramp,add=T)}else{
meshDist(myh2pos$obm$colMesh,alpha=alpha,distvec=myh2pos$obm$dists,from=distrange[1],to=distrange[2],scaleramp=scaleramp,add=T)
}
if(is.null(links)==F){lineplot(t(myh2pos$obm$colMesh$vb[-4,]),links,col=collinks,lwd=lwd)}
next3d()
text3d(0,0,0, texts[2])
next3d()
plot3d(themax*1.2,box=F,axes=F,col="white",xlab="",ylab="",zlab="",type="s",size=0,aspect=F)
if(heateuc==F){meshDist(myh3pos$obm$colMesh,alpha=alpha,distvec=myh3pos$obs3,from=distrange[1],to=distrange[2],scaleramp=scaleramp,add=T)}else{
meshDist(myh3pos$obm$colMesh,alpha=alpha,distvec=myh3pos$obm$dists,from=distrange[1],to=distrange[2],scaleramp=scaleramp,add=T)
}
if(is.null(links)==F){lineplot(t(myh3pos$obm$colMesh$vb[-4,]),links,col=collinks,lwd=lwd)}
next3d()
text3d(0,0,0, texts[3])
next3d()
plot3d(themax*1.2,box=F,axes=F,col="white",xlab="",ylab="",zlab="",type="s",size=0,aspect=F)
if(heateuc==F){meshDist(myh1neg$obm$colMesh,alpha=alpha,distvec=myh1neg$obs3,from=distrange[1],to=distrange[2],scaleramp=scaleramp,add=T)}else{
meshDist(myh1neg$obm$colMesh,alpha=alpha,distvec=myh1neg$obm$dists,from=distrange[1],to=distrange[2],scaleramp=scaleramp,add=T)
}
if(is.null(links)==F){lineplot(t(myh1neg$obm$colMesh$vb[-4,]),links,col=collinks,lwd=lwd)}
next3d()
text3d(0,0,0, texts[4])
next3d()
plot3d(themax*1.2,box=F,axes=F,col="white",xlab="",ylab="",zlab="",type="s",size=0,aspect=F)
if(heateuc==F){meshDist(myh2neg$obm$colMesh,alpha=alpha,distvec=myh2neg$obs3,from=distrange[1],to=distrange[2],scaleramp=scaleramp,add=T)}else{
meshDist(myh2neg$obm$colMesh,alpha=alpha,distvec=myh2neg$obm$dists,from=distrange[1],to=distrange[2],scaleramp=scaleramp,add=T)
}
if(is.null(links)==F){lineplot(t(myh2neg$obm$colMesh$vb[-4,]),links,col=collinks,lwd=lwd)}
next3d()
text3d(0,0,0, texts[5])
next3d()
plot3d(themax*1.2,box=F,axes=F,col="white",xlab="",ylab="",zlab="",type="s",size=0,aspect=F)
if(heateuc==F){meshDist(myh3neg$obm$colMesh,alpha=alpha,distvec=myh3neg$obs3,from=distrange[1],to=distrange[2],scaleramp=scaleramp,add=T)}else{
meshDist(myh3neg$obm$colMesh,alpha=alpha,distvec=myh3neg$obm$dists,from=distrange[1],to=distrange[2],scaleramp=scaleramp,add=T)
}
if(is.null(links)==F){lineplot(t(myh3neg$obm$colMesh$vb[-4,]),links,col=collinks,lwd=lwd)}
next3d()
text3d(0,0,0, texts[6])
}
}
dimnames(allshapes)[[3]]<-c("pc1pos","pc1neg","pc2pos","pc2neg","pc3pos","pc3neg")
if(out==T){
if(heat==T){outp<-list(pc1pos=myh1pos$obm$colMesh,pc1neg=myh1neg$obm$colMesh,pc2pos=myh2pos$obm$colMesh,pc2neg=myh2neg$obm$colMesh,pc3pos=myh3pos$obm$colMesh,pc3neg=myh3neg$obm$colMesh,distrange=distrange)}else{outp<-allshapes}
}else{outp<-NULL}
outp
}
#' @export
conslinks<-function(number,open=T){
k=seq(1:(number-1))
aw=NULL
for(i in k){
b=list(c(i,i+1))
aw<-c(aw,b)
}
if(open==T){return(aw)}else{
aw<-c(aw,list(c(1,number)))
}
aw
}
#' @export
list2array<-function(mylist){
require(abind)
final<-NULL
for(i in 1:length(mylist)){
temp<-array(mylist[[i]],dim=c(nrow(mylist[[i]]),ncol(mylist[[i]]),1))
final<-abind::abind(final,temp)
}
return(final)
}
#' @export
serpred<-function(dep,indep,polyn=1,length.out=10,predat=NULL,plot=F){
if(is.null(predat)){sizes<-seq(min(indep),max(indep),length.out=length.out)}else{sizes<-predat}
sizesvetts<-NULL
thelm<-lm(dep~poly(indep,degree=polyn,raw=TRUE))
for(i in 1:length(sizes)){
if(polyn==1){sizesvettsi<-c(1,sizes[i])}else{sizesvettsi<-c(1,poly(c(sizes[i]),degree=polyn,raw=T))}
sizesvetts<-rbind(sizesvetts,sizesvettsi)}
thepredicts<-NULL
for(j in 1:nrow(sizesvetts)){
thepredictsj<-c(crossprod(coef(thelm),sizesvetts[j,]))
thepredicts<-rbind(thepredicts,thepredictsj)
}
if(plot==T){
plotancova(dep,indep,legend=F,xlim=range(c(indep,predat)),ylim=range(c(dep,thepredicts)),polyn=polyn)
points(cbind(sizes,thepredicts),col=2,pch=2,cex=2)
}
thepredicts
}
#' @export
plot2dhull<-function(matrix,group,scalevariable=1,asp=1,extl=F,legend=T,xl=NULL,yl=NULL,posl=c("topright"),labels=NULL,clabel=0,pch=19,lwd=0.5,colhull=NULL,col=as.numeric(group),xlab=NULL,ylab=NULL,grey=F,xlimi=range(matrix[,1])*1.3,ylimi=range(matrix[,2])*1.3,reorder=T,alpha=rep(0,nlevels(group))){
library(MASS)
library(compositions)
library(spatstat)
library(gdata)
library(TeachingDemos)
library(ade4)
library(calibrate)
if (!is.factor(group)) stop("'group' must be a factor")
makeTransparent = function(..., alpha=0.5) {
if(alpha<0 | alpha>1) stop("alpha must be between 0 and 1")
alpha = floor(255*alpha)
newColor = col2rgb(col=unlist(list(...)), alpha=FALSE)
.makeTransparent = function(col, alpha) {
rgb(red=col[1], green=col[2], blue=col[3], alpha=alpha, maxColorValue=255)
}
newColor = apply(newColor, 2, .makeTransparent, alpha=alpha)
return(newColor)
}
x<-matrix[,1]
y<-matrix[,2]
if(reorder==T){group<-factor(group,levels=unique(group))}
species<-as.numeric(group)
col=col
if(grey==T){col=col2grey(col)}
coli<-data.frame(group,as.numeric(group),col,pch)
if(is.null(colhull)){colhull=aggregate(coli[,-1],by=list(coli[,1]),mean)[,-3][,2]}else{colhull=colhull}
colim<-cbind(aggregate(coli[,-1],by=list(coli[,1]),mean)[,-c(3,4)],colhull)
plot(x,y,type="p",cex=scalevariable,col=as.character(coli[,3]),pch=pch,xlim=xlimi,ylim=ylimi,xlab=xlab,ylab=ylab,asp=asp)
if(!is.null(labels)){textxy(x,y,labels)}else{NULL}
for(i in 1:(max(species)))
{
abline(h=0,v=0)
if(length(x[species==i])<3){NULL}else{
hulli<-convexhull.xy(subset(cbind(x,y),species==i))
par(new=T)
plot(hulli,col=makeTransparent(colim[,3][i], 1,alpha=alpha[i]),lwd=lwd,lty=i,add=T,asp=asp)
par(new=T)
daplotx<-c(hulli$bdry[[1]][[1]][length(hulli$bdry[[1]][[1]])],hulli$bdry[[1]][[1]])
daploty<-c(hulli$bdry[[1]][[2]][length(hulli$bdry[[1]][[2]])],hulli$bdry[[1]][[2]])
par(new=T)
plot(daplotx,daploty,lwd=lwd,type="l",lty=i,col=colim[,3][i],xlim=xlimi,ylim=ylimi,axes=F,xlab="",ylab="",asp=asp)
}
}
opar <- par(mar = par("mar"))
par(mar = c(0.1, 0.1, 0.1, 0.1))
on.exit(par(opar))
meansx<-aggregate(matrix[,1],by=list(group),FUN=mean)
meansy<-aggregate(matrix[,2],by=list(group),FUN=mean)
if (clabel > 0)
for(i in 1:nlevels(group)){ scatterutil.eti(meansx[i,-1],meansy[i,-1], meansx[,1][i], clabel,coul=1)}
if(legend==T){
if(extl==T){
x11()
plot(x,y,col="white")
legend(min(x),max(y),unique(coli[,1]), cex=1, col=unique(coli[,3]), pch=unique(coli[,4]),box.col="white")
}else{
if(is.null(xl)==T&is.null(yl)==T){
legend(posl,legend=unique(coli[,1]), col=unique(coli[,3]), pch=unique(coli[,4]),bty='n')}else{
legend(xl,yl,unique(coli[,1]),col=unique(coli[,3]),pch=unique(coli[,4]),bty='n')}
}
}
}
#' @export
reparray<-function(array,n,type=c("sequ","block")){
if(is.null(dimnames(array)[3])){dimnames(array)[3]<-list(c(1:dim(array)[3]))}else{dimnames(array)[3]<-dimnames(array)[3]}
require(abind)
myarray2<-NULL
steps<-n
if(type=="sequ"){
for(i in 1:dim(array)[3]){
temp1<-rep(list(array(array[,,i],dim=c(dim(array)[1],dim(array)[2],1))),steps)
temp2<-list2array(temp1)
myarray2<-abind::abind(myarray2,temp2)
}
return(myarray2)}else{NULL}
if(type=="block"){
temp1<-list2matrix(rep(array2list(array),n))
temp2<-matrix2arrayasis(temp1,dim(array)[1])
return(temp2)
}else{NULL}
}
#' @export
pt2dvaruno<-function(mua,mub,va,doopa=T,sss=T,tol=0.001){
library(shapes)
library(matrixcalc)
if(doopa==T){
theopamuamub<-procOPA(mub,mua,scale=F,reflect=F)
mua<-theopamuamub$Bhat
rotmua<-theopamuamub$R
if(rotmua[1,1]+1<tol){mua<--mua}
if(sss==F){
va<-va-((matrix.trace(t(mua)%*%va)/matrix.trace(t(mua)%*%mua))*mua)
}
theopavamua<-procOPA(mua,va,scale=F)
va<-theopavamua$Bhat
rotva<-theopavamua$R
if(rotva[1,1]+1<tol){va<--va}
}else{
if(sss==F){
va<-va-((matrix.trace(t(mua)%*%va)/matrix.trace(t(mua)%*%mua))*mua)
}
rotmua<-NULL
rotva<-NULL
}
muaz <- complex(real = mua[,1], imaginary = mua[,2])
mubz <- complex(real = mub[,1], imaginary = mub[,2])
vaz <- complex(real = va[,1], imaginary = va[,2])
muazbarra<-muaz/(sqrt(muaz%*%(Conj(muaz))))
mubzbarra<-mubz/(sqrt(mubz%*%(Conj(mubz))))
if(sss==T){
vbz<-vaz-1i*((Im(Conj(mubzbarra)%*%vaz))/(1+Conj(muazbarra)%*%mubzbarra))*(muazbarra+mubzbarra)}else{
vbz<-vaz-1i*((Im(Conj(mubzbarra)%*%vaz))/(1+Conj(muazbarra)%*%mubzbarra))*(muazbarra+mubzbarra)-((Re(Conj(mubzbarra)%*%vaz))/(1+Conj(muazbarra)%*%mubzbarra))*(muazbarra+mubzbarra)}
vazbarra<-vaz/(sqrt(vaz%*%(Conj(vaz))))
eps<-Re(sqrt(2*Im(Conj(mubzbarra)%*%vazbarra)^2/(1+Conj(muazbarra)%*%mubzbarra)))
eps2<-Re(sqrt(2*Mod(Conj(mubzbarra)%*%vazbarra)^2/(1+Conj(muazbarra)%*%mubzbarra)))
vb<-cbind(Re(vbz),Im(vbz))
out<-list(vb=vb,rotmua=rotmua,rotva=rotva,eps=eps,eps2=eps2)
return(out)
}
#' @export
opaloop<-function(GM,array,scale=F,reflect=F){
library(abind)
k<-dim(array)[1]
m<-dim(array)[2]
n<-dim(array)[3]
looped<-NULL
rots<-NULL
for(i in 1:n){
opai<-procOPA(GM,array[,,i],scale=scale,reflect=reflect)
loopedi<-array(opai$Bhat,dim=c(k,m,1))
looped<-abind::abind(looped,loopedi)
rotsi<-opai$R
rots<-c(rots,list(rotsi))
}
out<-list(looped=looped,rots=rots)
out
}
#' @export
diffonmesh<-function(lmsource,lmtarget,triang,colsource=1,alphas=c(0.7,0.3),grid=F,aff=F,nonaff=F,distscale=1, scaleramp=F,from = NULL,to = NULL,tol=NULL,sign=F,title=T,displace = FALSE,plotsource=T, steps = 20,rampcolors = colorRamps::blue2green2red(steps - 1),graph=T){
for ( i in c("Morpho","Rvcg","rgl")) {
if (!require(i,character.only = TRUE))
stop(paste0("please install package ",i))
}
lmsource <- rotonto(lmtarget,lmsource)$yrot
thetarget <- t(cbind(lmtarget,1))
thesource <- t(cbind(lmsource,1))
getPlotNr <- length(which(c(aff,nonaff,grid) == TRUE))+1
layout(matrix(1:getPlotNr,1,getPlotNr,byrow = T))
triangolazioni <- triang
thetarget_mesh <- list(vb=thetarget,it=triangolazioni)
class(thetarget_mesh) <- "mesh3d"
class(thetarget_mesh) <- "mesh3d"
thesource_mesh<-list(vb=thesource,it=triangolazioni)
class(thesource_mesh) <- "mesh3d"
class(thesource_mesh) <- "mesh3d"
if(grid==T){
shade3d(thetarget_mesh)
###between pointclouds
distvec <- sqrt(rowSums((lmsource-lmtarget)^2))/distscale
if(graph==T){
them<-meshDist(lmtarget,distvec = distvec,from=from,to=to,tol=tol,sign=sign,steps = 20,rampcolors =rampcolors,scaleramp=scaleramp)
}else{
them<-meshDist(lmtarget,distvec = distvec,from=from,to=to,tol=tol,sign=sign,steps = 20,rampcolors =rampcolors,shade=F,scaleramp=scaleramp)
}
##if you want add a deformed cube
if(graph==T){
deformGrid3d(lmsource,lmtarget,type="p",size=10,ngrid = 5,add=T)
title("Distance between point-clouds")}
}
### between surfaces
thetarget_warped <- tps3d(thesource_mesh,lmsource,lmtarget,lambda = 1e-8)
distvec2 <- sqrt(rowSums((vert2points(thesource_mesh)-vert2points(thetarget_warped))^2))/distscale
if(graph==T){
if(plotsource==T){shade3d(thesource_mesh,col=colsource,alpha=alphas[2])}
them<-meshDist(thetarget_mesh,thesource_mesh,distvec=distvec2,alpha=alphas[1],add=T,from=from,to=to,tol=tol,scaleramp=scaleramp,sign=sign,displace=displace,steps = 20,rampcolors =rampcolors)
if(title==T){title3d("total")
title("total")}
}else{
them<-meshDist(thetarget_mesh,distvec=distvec2,alpha=alphas[1],add=T,from=from,to=to,tol=tol,scaleramp=scaleramp,sign=sign,displace=displace,steps = 20,rampcolors =rampcolors,shade=F)
}
if (nonaff || aff) {
affinetrafo <- computeTransform(vert2points(thetarget_mesh),vert2points(thesource_mesh),type="a")
affineshort <- applyTransform(thesource_mesh,affinetrafo)
if(nonaff) {### visualize non-affine deform between surfaces
##create affine transfrom to longnose.mesh to remove affine differences and calculate distance from affine transform to target
distvec3 <- sqrt(rowSums((vert2points(affineshort)-vert2points(thetarget_warped))^2))/distscale
if(graph==T){
open3d()
shade3d(thetarget_mesh,col=1,alpha=alphas[2])
themna<-meshDist(thesource_mesh,distvec=distvec3,add=T,from=from,to=to,tol=tol,scaleramp=scaleramp,sign=sign,displace=displace,steps = 20,rampcolors =rampcolors)
if(title==T){title3d("non affine")
title("non-affine")}
}else{
themna<-meshDist(thesource_mesh,distvec=distvec3,add=T,from=from,to=to,tol=tol,scaleramp=scaleramp,sign=sign,displace=displace,steps = 20,rampcolors =rampcolors,shade=F)
}
}
if(aff) {### the affine transform looks like that:
distaffnonaffon4<-sqrt(rowSums((vert2points(affineshort)-vert2points(thesource_mesh))^2))/distscale
if(graph==T){
open3d()
shade3d(thetarget_mesh,col=1,alpha=alphas[2])
thema<-meshDist(thesource_mesh,distvec=distaffnonaffon4,add=T,from=from,to=to,tol=tol,scaleramp=scaleramp,sign=sign,displace=displace,steps = 20,rampcolors =rampcolors)
if(title==T){title3d("affine")
title("affine")}
}else{
thema<-meshDist(thesource_mesh,distvec=distaffnonaffon4,add=T,from=from,to=to,tol=tol,scaleramp=scaleramp,sign=sign,displace=displace,steps = 20,rampcolors =rampcolors,shade=F)
}
}
}else{themna<-NULL
thema<-NULL
}
out<-list(obm=them,themna=themna,thema=thema)
}
#' @export
mycca<-function(x,y,pch=19,col=1,group=NULL,labels=NULL,extl=F,legend=T,xl=NULL,yl=NULL,posl=c("topright"),cex=1,xlab=NULL,ylab=NULL,asp=NULL){
require(CCA)
require(vegan)
require(calibrate)
myrda<-vegan::rda(y~x)
rsqadj<-RsquareAdj(myrda)
print(rsqadj)
anovarda<-anova(myrda,step=2000)
print(anovarda)
if(is.null(xlab)==T){xlab="Independent"}else{xlab=xlab}
if(is.null(ylab)==T){ylab="Dependent"}else{ylab=ylab}
if (!is.null(group)){
col=col
species<-as.numeric(group)
coli<-data.frame(group,as.numeric(group),col,pch)
colim<-aggregate(coli[,-1],by=list(coli[,1]),mean)}else{
coli<-data.frame(1,1,col,pch)
colim<-aggregate(coli[,-1],by=list(coli[,1]),mean)
}
thecca<-rcc(as.matrix(x),as.matrix(y),0.1,0.1)
if(is.null(group)==T){plot(x,rcc(as.matrix(x),as.matrix(y),0.1,0.1)$scores$yscores[,1],col=col,pch=pch,cex=cex,xlab=xlab,ylab=ylab)}
if (!is.null(group)){
plot2dhull(cbind(x,thecca$scores$yscores[,1]),group,scalevariable=cex,pch=pch,col=col,colhull=colim[,3],clabel=0,legend=F,reorder=F,xlimi=range(x),ylimi=range(thecca$scores$yscores[,1]),asp=asp,xlab=xlab,ylab=ylab)
}
if (!is.null(labels)){textxy(x,thecca$scores$yscores,labels)}
if(legend==T){
if(extl==T){
x11()
plot(x,y,col="white")
legend(min(x),max(y),colim[,1], cex=1, col=colim[,3], pch=colim[,4],box.col="white")
}else{
if(is.null(xl)==T&is.null(yl)==T){
legend(posl,legend=colim[,1], cex=1, col=colim[,3], pch=colim[,4],bty='n')}else{
legend(xl,yl,colim[,1], cex=1, col=colim[,3], pch=colim[,4],bty='n')}
}
}
if (!is.null(group)){
if(dim(as.matrix(y))[2]>1){bygroup<-manymultgr(y,x,group)}else{
bygroup<-plotancova(y,x,group,legend=F,plot=F)$sint_results}
}else{bygroup=c("no group structure")}
out<-list(xscores=thecca$scores$xscores,yscores=thecca$scores$yscores,AdjRsq=rsqadj,anovarda=anovarda,bygroup=bygroup)
return(out)
}
#' @export
posfac<-function(factor){
positions<-NULL
for(k in 1:max(as.numeric(factor))){
factor<-factor(factor,levels=unique(factor))
positioni<-which(as.numeric(factor)==k)
positions<-c(positions,list(positioni))
}
names(positions)<-levels(factor)
sequ<-NULL
for(i in 1:max(length(positions))){
seqi<-c(1:length(positions[[i]]))
sequ<-c(sequ,list(seqi))
}
pure<-rep(NA,length(unlist(positions)))
for(j in 1:max(length(positions))){
pure<-replace(pure,positions[[j]],c(sequ[[j]]))
}
pure}
#' @export
areasip<-function(matrix,links=NULL,PB=NA,PA=NA,S=NA,SB=NA,H=NA,V=3,a=NULL,q=NULL,Y=FALSE,j=FALSE,D=FALSE,St=Inf,Q=TRUE,graph=T,extpol=F){
if(!is.null(links)){warning("Links **must** identify, among other structures, a closed contour")}
library(geometry)
library(tripack)
library(geoR)
library(gstat)
library(sp)
library(fields)
library(RTriangle)
mate<-matrix
if(!is.null(links)&is.na(S)==T){S<-list2matrix(links)}
if(extpol==T){
step1<-pslg(mate,S=S,PB=PA,PA=PA,SB=SB,H=NA)
posimpnoh<-RTriangle::triangulate(step1,V=V,a=NULL,S=St,q=q,Y=Y,j=j,D=D,Q=Q)
step2<-matrix[unique(posimpnoh$E[ posimpnoh$EB%in%c(1)]),]
maulinks<-posimpnoh$E[posimpnoh$EB%in%c(1),]
newmaulinks<-maulinks
newmaulinks[,1]<-c(1:nrow(maulinks))
newmaulinks[,2]<-match(maulinks[,2],maulinks[,1])
cycles2 <- function(links) {
require(ggm)
if (any(is.na(links))) stop("missing value/s in links")
mat <- matrix(0L, max(links), max(links))
mat[links] <- 1
lapply(ggm::fundCycles(mat), function(xa) rev(xa[ , 1L]))
}
mypol<-step2[cycles2(newmaulinks)[[1]],]}else{mypol<-c("you do not want external polygon")}
p<-pslg(mate,PB=PB,PA=PA,S=S,SB=SB,H=H)
au<-RTriangle::triangulate(p,V=1,a=a,q=q,Y=Y,j=j,D=D,S=St,Q=Q)
if(graph==T){
plot(au,asp=1)
lineplot(mate,links,col=3)
points(au$P,pch=19)
points(mate,pch=21,cex=2)
}
centros<-NULL
areas<-NULL
for(i in 1:nrow(au$T)){
centrosi<-apply(au$P[au$T[i,],],2,mean)
centros<-rbind(centros,centrosi)
areasi<-convhulln(au$P[au$T[i,],],option="FA")$vol
areas<-c(areas,areasi)
}
if(graph==T){points(centros,col=2)}
M<-centros
thecentr<-centroids(mate)
dasomm<-NULL
for(i in 1:nrow(M)){
dasommi<-dist(rbind(thecentr,M[i,]))^2*areas[i]
dasomm<-c(dasomm,dasommi)
}
ip<-sum(dasomm)
are<-sum(areas)
deltri<-au$T
ptri<-au$P
origcentros<-rbind(mate,M)
out<-list(area=are,areas=areas,ip=ip,centros=centros,deltri=deltri,ptri=ptri,origcentros=origcentros,triangob=au,ext=mypol)
out
}
#' @export
array2list<-function(array){
thelist<-NULL
for(i in 1:dim(array)[3]){
eli<-array[,,i]
thelist<-c(thelist,list(eli))
}
if(is.null(dimnames(array)[[3]])==F){names(thelist)<-dimnames(array)[[3]]}
return(thelist)
}
#' @export
array2mat<-function(array,inds,vars){
if(class(array)=="matrix"){array<-array(array,dim=c(nrow(array),ncol(array),1))}
X1 <-aperm(array,c(3,2,1))
dim(X1)<- c(inds, vars)
if(!is.null(dimnames(array)[3])){rownames(X1)<-unlist(dimnames(array)[3])}else{rownames(X1)<-c(1:nrow(X1))}
return(X1)
}
#' @export
biharm.new<-function(vec1,vec2){
dim<-length(vec1)
n <- sqrt(sum((vec1-vec2)^2))
if(dim<3){
n^2*log(n^2)}else{
-n
}
}
#' @export
centershapes<-function(array){
if(is.matrix(array)==T){array<-array(array,dim=c(nrow(array),ncol(array),1))}
centros<-centroids(array)
k<-dim(array)[1]
m<-dim(array)[2]
n<-dim(array)[3]
prov<-array(rep(t(centros),each=k),dim=c(k,m,n))
array2<-array-prov
return(array2)
}
#' @export
centroids<-function(array){
meanmat<-function(mat){
mean<-apply(mat,2,mean)
mean
}
if(class(array)=="matrix"){array<-array(array,dim=c(dim(array)[1],dim(array)[2],1))}
centroidi<-t(apply(array,3,meanmat))
return(centroidi)
}
#' @export
firstsfac<-function(group){
firstsres<-NULL
for (i in levels(group)){
firstresi<-min(which(group==i))
firstsres<-rbind(firstsres,firstresi)
}
return(firstsres)
}
#' @export
heat2d<-function(init,fin,invc=F,constr=T,sen=F,logsen=T,nadd=5000,linkss=NULL,zlim=NULL,legend=T,plottarget=T,ext=0.1,collinkss=1,lwds=2,collinkst=2,lwdt=2,pchs=19,pcht=19,cexs=0.5,cext=0.5,colors=c("blue4","cyan2","yellow","red4"),alpha=1,ngrid=30,oma=c(5,5,5,5),mai=c(3,3,3,3),mar=c(3,3,3,3),mag=1,tol=0.1,graphics=T,PB=NA,PA=NA,S=NA,SB=NA,H=NA,V=3,a=NULL,q=NULL,Y=FALSE,j=FALSE,D=FALSE,St=Inf,Q=TRUE,pholes=NA){
library(shapes)
library(Morpho)
library(matrixcalc)
library(geometry)
library(tripack)
library(geoR)
library(gstat)
library(sp)
library(fields)
library(RTriangle)
library(sp)
library(alphahull)
library(ggm)
if(!is.na(H)&is.na(pholes)==T){stop("If you specify 'H' you need to specify 'pholes' also")}
if(!is.na(pholes)&is.null(linkss)==T){stop("If you specify 'H' and 'pholes' you need to specify 'linkss' also")}
mate<-fin
mate2<-init
fin<-init+(fin-init)*mag
if(!is.null(linkss)&is.na(S)==T){S<-list2matrix(linkss)}
posimpfin<-areasip(fin,S=S,PB=PA,PA=PA,SB=SB,H=H,V=V,a=NULL,q=q,Y=Y,j=j,D=D,St=St,Q=Q,graph=F,extpol=F)
posimp<-areasip(init,S=S,PB=PA,PA=PA,SB=SB,H=H,V=V,a=NULL,q=q,Y=Y,j=j,D=D,St=St,Q=Q,graph=F,extpol=F)
pofin<-areasip(fin,S=S,PB=PA,PA=PA,SB=SB,H=H,V=V,a=a,q=q,Y=Y,j=j,D=D,St=St,Q=Q,graph=F,extpol=F)
po<-areasip(init,S=S,PB=PA,PA=PA,SB=SB,H=H,V=V,a=a,q=q,Y=Y,j=j,D=D,St=St,Q=Q,graph=F,extpol=F)
posimpnohfin<-areasip(fin,S=S,PB=PA,PA=PA,SB=SB,H=NA,V=V,a=NULL,q=q,Y=Y,j=j,D=D,St=St,Q=Q,graph=F,extpol=F)
posimpnoh<-areasip(init,S=S,PB=PA,PA=PA,SB=SB,H=NA,V=V,a=NULL,q=q,Y=Y,j=j,D=D,St=St,Q=Q,graph=F,extpol=F)
matr<-init
M<-po$centros
areas<-po$areas
tpsgrid<-tpsgridpaolo(init,fin,linksTT=linkss,linksYY=linkss,axes2d=T,ext=ext,graphics=F,ngrid=ngrid)
veclist<-NULL
for(j in 1:nrow(M)){
vecj<-NULL
for(i in 1:nrow(matr)){
vec1ji<-2*(M[j,1]-matr[i,1])+2*(M[j,1]-matr[i,1])*log((M[j,1]-matr[i,1])^2+(M[j,2]-matr[i,2])^2)
vec2ji<-2*(M[j,2]-matr[i,2])+2*(M[j,2]-matr[i,2])*log((M[j,1]-matr[i,1])^2+(M[j,2]-matr[i,2])^2)
vecji<-c(vec1ji,vec2ji)
vecj<-rbind(vecj,vecji)
}
veclist<-c(veclist,list(vecj))
}
jac<-NULL
for(i in 1: length(veclist)){
jaci<-t(tpsgrid$B)+t(tpsgrid$W)%*%veclist[[i]]
jac<-c(jac,list(jaci))
}
deltus<-NULL
dpl<-NULL
sens<-NULL
for(i in 1:length(jac)){
deltusi<-jac[[i]]-diag(nrow(jac[[i]]))
dpli<-0.5*(deltusi+t(deltusi))
seni<-0.5*matrix.trace(t(dpli)%*%dpli)
deltus<-c(deltus,list(deltusi))
dpl<-c(dpl,list(dpli))
sens<-c(sens,seni)
}
lsens<-log2(sens)
sens2<-sens*areas
jac2<-list2array(jac)
mean11<-mean(jac2[1,1,])
mean12<-mean(jac2[1,2,])
mean21<-mean(jac2[2,1,])
mean22<-mean(jac2[2,2,])
detmean<-det(matrix(c(mean11,mean21,mean12,mean22),ncol=2))
myj<-unlist(lapply(jac,det))
if(sen==F){obs<-log2(myj/detmean)}else{
if(logsen==T){obs=lsens}else{obs=sens}
}
fit<- Tps(M, obs)
summary(fit)
if(constr==F){
hull <- chull(fin)
indx=hull
if(invc==F){points <- fin[indx,]}else{points <- init[indx,]}
points <- rbind(points,points[1,])
sfe1<-spsample(Polygon(points),nadd,type="regular")
sr1<-sfe1@coords}else{
if(invc==F){mau<-fin[unique(posimpnoh$triangob$E[posimpnoh$triangob$EB%in%c(1)]),]}else{mau<-init[unique(posimpnoh$triangob$E[posimpnoh$triangob$EB%in%c(1)]),]}
matepro<-fin
matepro[1:nrow(fin)%in%unique(posimpnoh$triangob$E[posimpnoh$triangob$EB%in%c(1)])==F,]<-NA
faclinks<-as.factor(is.na(matepro[,1]))
newindlinks<-posfac(faclinks)
maulinks<-posimpnoh$triangob$E[posimpnoh$triangob$EB%in%c(1),]
#plotmyarrays(mate,links=array2list(matrix2arrayasis(maulinks,1)))
newmaulinks<-maulinks
newmaulinks[,1]<-c(1:nrow(maulinks))
newmaulinks[,2]<-match(maulinks[,2],maulinks[,1])
#plotmyarrays(mau,links=array2list(matrix2arrayasis(newmaulinks,1)),xlim=c(25,70))
cycles2 <- function(links) {
require(ggm)
if (any(is.na(links))) stop("missing value/s in links")
mat <- matrix(0L, max(links), max(links))
mat[links] <- 1
lapply(ggm::fundCycles(mat), function(xa) rev(xa[ , 1L]))
}
#plotmyarrays(mau[cycles2(newmaulinks)[[1]],],links=conslinks(nrow(mau),open=F))
sfe1<-spsample(Polygon(rbind(mau[cycles2(newmaulinks)[[1]],],mau[cycles2(newmaulinks)[[1]],][1,])),nadd,type="regular")
sr1<-sfe1@coords
}
#points(sr1)
pred<-predict(fit,sr1)
if(is.null(zlim)==T){zlim<-range(pred)}else{zlim<-zlim}
if(graphics==T){
par(oma=oma,mai=mai,mar=mar)
if(legend==F){ima<-image(xyz2img(cbind(sr1,pred),tolerance=tol),asp=1,xlim=tpsgrid$grid$xlim,ylim=tpsgrid$grid$ylim,col=makeTransparent(colorRampPalette(colors)(n = length(obs)),alpha=alpha),xaxs="r",yaxs="r",xaxt="n",yaxt="n",bty="n",xlab="",ylab="",zlim=zlim)}else{
ima<-image.plot(xyz2img(cbind(sr1,pred),tolerance=tol),asp=1,xlim=tpsgrid$grid$xlim,ylim=tpsgrid$grid$ylim,col=makeTransparent(colorRampPalette(colors)(n = length(obs)),alpha=alpha),xaxs="r",yaxs="r",xaxt="n",yaxt="n",bty="n",xlab="",ylab="",zlim=zlim)
}
par(new=T)
plot(fin,xlim=tpsgrid$grid$xlim,ylim=tpsgrid$grid$ylim,asp=1,pch=pchs,cex=cexs,xaxt="n",yaxt="n",bty="n",xlab="",ylab="")
if(is.null(linkss)==F){lineplot(fin,linkss,col=collinkss,lwd=lwds)}
lines(tpsgrid$grid$ngrid,xlim=tpsgrid$grid$xlim,ylim=tpsgrid$grid$ylim)
lines(tpsgrid$grid$ngrid2,xlim=tpsgrid$grid$xlim,ylim=tpsgrid$grid$ylim)
if(plottarget==T){
par(new=T)
plot(init,xlim=tpsgrid$grid$xlim,ylim=tpsgrid$grid$ylim,asp=1,pch=pcht,cex=cext,xaxt="n",yaxt="n",bty="n",xlab="",ylab="")
lineplot(init,linkss,col=collinkst,lwd=lwdt)
}
}
pols<-NULL
pols2<-NULL
if(!is.na(pholes)){
for(i in 1:length(pholes)){
myindex<-which(apply(matrix(list2matrix(linkss)%in%pholes[[i]],ncol=2),1,sum)>1)
matili<-list2matrix(linkss)[myindex,]
holi<-init[unique(sort(c(unique(matili[,1]),unique(matili[,2])))),]
hol2i<-fin[unique(sort(c(unique(matili[,1]),unique(matili[,2])))),]
si<-matrix(as.numeric(ordered(matili)),ncol=2)
holtriangi<-RTriangle::triangulate(pslg(holi,S=si))
holtriang2i<-RTriangle::triangulate(pslg(hol2i,S=si))
holui<-holi[unique(holtriangi$E[ holtriangi$EB%in%c(1)]),]
holu2i<-hol2i[unique( holtriang2i$E[ holtriang2i$EB%in%c(1)]),]
holiproi<-holi
holiproi[1:nrow(holi)%in%unique(holtriangi$E[holtriangi$EB%in%c(1)])==F,]<-NA
faclinkshi<-as.factor(is.na(holiproi[,1]))
newindlinkshi<-posfac(faclinkshi)
hlinks<-holtriangi$E[holtriangi$EB%in%c(1),]
newhlinks<-hlinks
newhlinks[,1]<-c(1:nrow(hlinks))
newhlinks[,2]<-match(hlinks[,2],hlinks[,1])
pols<-c(pols,list(holui[cycles2(newhlinks)[[1]],]))
pols2<-c(pols2,list(holu2i[cycles2(newhlinks)[[1]],]))
if(graphics==T){
polygon(holui[cycles2(newhlinks)[[1]],],col="white",border=collinkss)
polygon(holu2i[cycles2(newhlinks)[[1]],],col="white",border=collinkst)}
}}
out<-list(mate=fin,mate2=init,centros=M,interpcoords=sr1,jacs=jac,detjac=myj,detmean=detmean,obs=obs,fit=fit,pred=pred,xlim=tpsgrid$grid$xlim,ylim=tpsgrid$grid$ylim,tol=tol,cols=makeTransparent(colorRampPalette(colors)(n = length(obs)),alpha=alpha),tpsgrid=tpsgrid,sumse=sum(sens2),pols=pols,pols2=pols2,areasipob=po)
out
}
#' @export
helmert<-function(p)
{H<-matrix(0, p, p)
diag(H)<--(0:(p-1)) * (-((0:(p-1))*((0:(p-1))+1))^(-0.5))
for (i in 2:p){H[i,1:(i-1)]<- -((i-1)*(i))^(-0.5)}
H[1,]<-1/sqrt(p)
H}
#' @export
lastsfac<-function(group){
lastsres<-NULL
for (i in levels(group)){
lastresi<-max(which(group==i))
lastsres<-rbind(lastsres,lastresi)
}
return(lastsres)
}
#' @export
list2matrix<-function(mylist){
final<-NULL
for(i in 1:length(mylist)){
temp<-mylist[[i]]
final<-rbind(final,temp)
}
#if(is.null(names(mylist))==F){rownames(final)<-names(mylist)}
return(final)
}
#' @export
makeTransparent = function(..., alpha=0.5) {
if(alpha<0 | alpha>1) stop("alpha must be between 0 and 1")
alpha = floor(255*alpha)
newColor = col2rgb(col=unlist(list(...)), alpha=FALSE)
.makeTransparent = function(col, alpha) {
rgb(red=col[1], green=col[2], blue=col[3], alpha=alpha, maxColorValue=255)
}
newColor = apply(newColor, 2, .makeTransparent, alpha=alpha)
return(newColor)
}
#' @export
manymultgr<-function(y,x,group,steps=5000){
res<-NULL
rsq<-NULL
adjrsq<-NULL
pval<-NULL
for(i in 1:max(as.numeric((group)))){
rsqi<-RsquareAdj(vegan::rda(y[as.numeric(group)==i,],x[as.numeric(group)==i]))$r.squared
adjrsqi<-RsquareAdj(vegan::rda(y[as.numeric(group)==i,],x[as.numeric(group)==i]))$adj.r.squared
pvali<-anova(vegan::rda(y[as.numeric(group)==i,],x[as.numeric(group)==i]),step=steps)$Pr[1]
rsq<-append(rsq,rsqi)
adjrsq<-append(adjrsq,adjrsqi)
pval<-append(pval,pvali)
}
Sig<-ifelse(pval<0.05,c("*"),c(""))
res<-data.frame(round(cbind(rsq,adjrsq,pval),digits=5),Sig)
rownames(res)<-levels(group)
return(res)
}
#' @export
mopa<-function(tt,yy,rot=c("mopa","opa"),CSinit = FALSE,volinit=FALSE,center=TRUE){
warning("the SECOND input matrix will be rotated on the FIRST one")
if(is.matrix(tt)==T){
k<-nrow(tt)
m<-ncol(tt)
ttar<-array(tt,dim=c(k,m,1))
yyar<-array(yy,dim=c(k,m,1))}else {
ttar=tt
yyar=yy
k<-dim(tt)[1]
m<-dim(tt)[2]
}
if(center==T){ttcs<-centershapes(ttar)[,,1]
yycs<-centershapes(yyar)[,,1]}
if(CSinit==T){
yycs<-scaleshapes(yycs)
ttcs<-scaleshapes(ttcs)
}else{yycs<-yycs;ttcs<-ttcs}
if(volinit==T){
at<-tpsdry2(ttcs,yycs,meth="mor",g11=F)$at
detat<-det(at)
if(ncol(yycs)>2){yycs<-yycs/(detat^(1/3))}else{yycs<-yycs/(detat^(1/2))}
}
if(rot=="opa"){rotprocstep1<-svd(t(ttcs)%*%yycs)}
stef<-CreateL(ttcs)
W<-(stef$Lsubk%*%yycs)
appo<-stef$Linv[(k+1):(k+(m+1)),1:k]%*%yycs
cT<-appo[1,]
At<-appo[-1,]
if(rot=="mopa"){
rotprocstep1<-svd(At)}
rotprocstep2<-rotprocstep1$v%*%t(rotprocstep1$u)
opizzata1<-yycs%*%rotprocstep2
opizzata2<-array(opizzata1,dim=c(k,m,1))[,,1,drop=F]####centershapes(array(opizzata1,dim=c(k,m,1)))[,,1,drop=F]
Wafter<-stef$Lsubk%*%opizzata2[,,1]
appoafter<-stef$Linv[(k+1):(k+(m+1)),1:k]%*%opizzata2[,,1]
cTafter<-appoafter[1,]
Atafter<-appoafter[-1,]
out=list(opizzata=opizzata2,W=W,cT=cT,At=At,Wafter=Wafter,cTafter=cTafter,Atafter=Atafter)
return(out)
}
#' @export
newmb<-function(source,target){
k<-nrow(source)
m<-ncol(source)
h<-1*(helmert(k)[-1,])
cm<-t(h)%*%h
source<-cm%*%source
target<-cm%*%target
s1<-sigm.new(source)
source<-h%*%source
target<-h%*%target
s<-h%*%s1%*%t(h)
S<-s
Q<-source
Gam_11=solve(S)-solve(S)%*%Q%*%solve(t(Q)%*%solve(S)%*%Q)%*%t(Q)%*%solve(S)
Gam_21=solve(t(Q)%*%solve(S)%*%Q)%*%t(Q)%*%solve(S)
Gam_22=-solve(t(Q)%*%solve(S)%*%Q)
W=Gam_11%*%target
At=Gam_21%*%target
out<-list(gamma11=Gam_11,gamma21=Gam_21,gamma22=Gam_22,W=W,At=At,sold=s1,snew=s,h=h)
return(out)
}
#' @export
opaloop2<-function(GM,array,scale=F,reflect=F){
library(abind)
library(Morpho)
k<-dim(array)[1]
m<-dim(array)[2]
n<-dim(array)[3]
looped<-NULL
rots<-NULL
for(i in 1:n){
print(i)
opai<-rotonto(GM,array[,,i],scale=scale,reflection=reflect,signref=T)
loopedi<-array(opai$yrot,dim=c(k,m,1))
looped<-abind::abind(looped,loopedi)
rotsi<-opai$gamm
rots<-c(rots,list(rotsi))
}
out<-list(looped=looped,rots=rots)
out
}
#' @export
ordiwithshapes<-function(mshape,scores,rotation,whichaxes=c(1,2),addata=NULL,asp=1,xlab="PC1",ylab="PC2",triang=NULL,factraj=NULL,coltraj=NULL,distscale=1,scaleramp=F,from=0,to=NULL,mag=1,procSym=T,subplotdim=2,legendgroup=NULL,shiftposx=1.5,shiftnegx=1.5,shiftposy=1.5,shiftnegy=1.5,links=NULL,collinks=1,grouphull=NULL,colhull=NULL,labels=NULL,labelscex=1,pch=19,col=1,cex=1,mult=0.5,tit3d=T,plotsource=T,xlim=extendrange(extendrange(range(scores[,whichaxes[1]]),f=mult)),ylim=extendrange(extendrange(range(scores[,whichaxes[2]]),f=mult))){
library(Morpho)
library(spatstat)
library(calibrate)
if(procSym==T){
pc1pos<-showPC(max(scores[,whichaxes[1]])*mag,rotation[,whichaxes[1]],mshape)
pc1neg<-showPC(min(scores[,whichaxes[1]])*mag,rotation[,whichaxes[1]],mshape)
pc2pos<-showPC(max(scores[,whichaxes[2]])*mag,rotation[,whichaxes[2]],mshape)
pc2neg<-showPC(min(scores[,whichaxes[2]])*mag,rotation[,whichaxes[2]],mshape)}else{
pc1pos<-showPC2(max(scores[,whichaxes[1]])*mag,rotation[,whichaxes[1]],mshape)
pc1neg<-showPC2(min(scores[,whichaxes[1]])*mag,rotation[,whichaxes[1]],mshape)
pc2pos<-showPC2(max(scores[,whichaxes[2]])*mag,rotation[,whichaxes[2]],mshape)
pc2neg<-showPC2(min(scores[,whichaxes[2]])*mag,rotation[,whichaxes[2]],mshape)
}
if(dim(mshape)[2]<3){
library(TeachingDemos)
library(Morpho)
library(calibrate)
library(spatstat)
library(shapes)
plot(scores[,whichaxes[1]],scores[,whichaxes[2]],axes=F,xlim=xlim,ylim=ylim,pch=pch,col=col,cex=cex,asp=asp,xlab=xlab,ylab=ylab)
axis(1,pos=0,at=round(seq(from=min(scores[,whichaxes[1]]),to=max(scores[,whichaxes[1]]),length.out=10),2))
axis(2,pos=0,at=round(seq(from=min(scores[,whichaxes[2]]),to=max(scores[,whichaxes[2]]),length.out=10),2))
if(!is.null(addata)){
par(new=T)
plot(addata,asp=asp,xlim=xlim,ylim=ylim,pch=pch,col=col,cex=cex,axes=F,xlab="",ylab="")}
if(!is.null(factraj)){
if(is.null(coltraj)){coltraj=rep(1,nlevels(factraj))}else{coltraj=coltraj}
for (i in 1 :nlevels(factraj)){
lines(scores[,whichaxes[1]][as.numeric(factraj)==i],scores[,whichaxes[2]][as.numeric(factraj)==i],xlim=xlim,ylim=ylim,col=coltraj[i])
par(new=T)
}
}
if(!is.null(labels)==T){textxy(scores[,whichaxes[1]],scores[,whichaxes[2]],labs=labels,cex=labelscex)}
if(!is.null(grouphull)==T){
if(is.null(colhull)==T){colhull<-c(1:nlevels(grouphull))}else{colhull<-colhull}
for(i in 1:(max(as.numeric(grouphull))))
{if(length(scores[,1][as.numeric(grouphull)==i])<3){NULL}
else{
hulli<-convexhull.xy(subset(cbind(scores[,whichaxes[1]],scores[,whichaxes[2]]),as.numeric(grouphull)==i))
par(new=T)
daplotx<-c(hulli$bdry[[1]][[1]][length(hulli$bdry[[1]][[1]])],hulli$bdry[[1]][[1]])
daploty<-c(hulli$bdry[[1]][[2]][length(hulli$bdry[[1]][[2]])],hulli$bdry[[1]][[2]])
plot(daplotx,daploty,lwd=2,type="l",lty=i,col=colhull[i],xlim=xlim,ylim=ylim,axes=F,asp=asp,xlab="",ylab="")
}
}}
subplot(tpsgridpaolo(mshape,pc1pos,linksYY=links,collinksYY=collinks,axes2d=F,displ=F,collinksTT=makeTransparent("white",alpha=0)),max(scores[,whichaxes[1]])*shiftposx,0,size=c(subplotdim,subplotdim))
subplot(tpsgridpaolo(mshape,pc1neg,linksYY=links,collinksYY=collinks,axes2d=F,displ=F,collinksTT=makeTransparent("white",alpha=0)),min(scores[,whichaxes[1]])*shiftnegx,0,size=c(subplotdim,subplotdim))
subplot(tpsgridpaolo(mshape,pc2pos,linksYY=links,collinksYY=collinks,axes2d=F,displ=F,collinksTT=makeTransparent("white",alpha=0)),0,max(scores[,whichaxes[2]])*shiftposy,size=c(subplotdim,subplotdim))
subplot(tpsgridpaolo(mshape,pc2neg,linksYY=links,collinksYY=collinks,axes2d=F,displ=F,collinksTT=makeTransparent("white",alpha=0)),0,min(scores[,whichaxes[2]])*shiftnegy,size=c(subplotdim,subplotdim))
}else{
if(is.null(triang)==F&is.null(links)==T){
open3d(windowRect=c(100,100,1000,1000))
mat <- matrix(1:8, ncol=2)
layout3d(mat, height = rep(c(2,1), 2), model = "inherit")
diffonmesh(mshape,pc1pos,triang,title=F,distscale=distscale,from=from,to=to,plotsource=plotsource,scaleramp=scaleramp)
next3d()
text3d(0,0,0,paste("Scores",whichaxes[1],"pos"))
next3d()
diffonmesh(mshape,pc2pos,triang,title=F,distscale=distscale,from=from,to=to,plotsource=plotsource,scaleramp=scaleramp)
next3d()
text3d(0,0,0,paste("Scores",whichaxes[2],"pos"))
next3d()
diffonmesh(mshape,pc1neg,triang,title=F,distscale=distscale,from=from,to=to,plotsource=plotsource,scaleramp=scaleramp)
next3d();text3d(0,0,0,paste("Scores",whichaxes[1],"neg"))
next3d()
diffonmesh(mshape,pc2neg,triang,title=F,distscale=distscale,from=from,to=to,plotsource=plotsource,scaleramp=scaleramp)
next3d();text3d(0,0,0,paste("Scores",whichaxes[2],"neg"))
}
if(is.null(links)==F){
open3d(windowRect=c(100,100,1000,1000))
mat <- matrix(1:8, ncol=2)
layout3d(mat, height = rep(c(2,1), 2), model = "inherit")
plot3D(pc1pos,add=T,bbox=F)
lineplot(pc1pos,links)
if(tit3d==T){next3d();text3d(0,0,0,paste("Scores",whichaxes[1],"pos"))}
next3d()
plot3D(pc2pos,add=T,bbox=F)
lineplot(pc2pos,links)
if(tit3d==T){next3d();text3d(0,0,0,paste("Scores",whichaxes[2],"pos"))}
next3d()
plot3D(pc1neg,add=T,bbox=F)
lineplot(pc1neg,links)
if(tit3d==T){next3d();text3d(0,0,0,paste("Scores",whichaxes[1],"neg"))}
next3d()
plot3D(pc2neg,add=T,bbox=F)
lineplot(pc2neg,links)
if(tit3d==T){next3d();text3d(0,0,0,paste("Scores",whichaxes[2],"neg"))}}
plot(scores[,whichaxes[1]],scores[,whichaxes[2]],axes=F,xlim=xlim,ylim=ylim,pch=pch,col=col,cex=cex,asp=asp,xlab=xlab,ylab=ylab)
axis(1,pos=0,at=round(seq(from=min(scores[,whichaxes[1]]),to=max(scores[,whichaxes[1]]),length.out=10),2))
axis(2,pos=0,at=round(seq(from=min(scores[,whichaxes[2]]),to=max(scores[,whichaxes[1]]),length.out=10),2))
if(!is.null(addata)){
par(new=T)
plot(addata,xlim=xlim,ylim=ylim,pch=pch,col=col,cex=cex,axes=F,asp=asp,xlab="",ylab="")}
if(!is.null(factraj)){
if(is.null(coltraj)){coltraj=rep(1,nlevels(factraj))}else{coltraj=coltraj}
for (i in 1 :nlevels(factraj)){
lines(scores[,whichaxes[1]][as.numeric(factraj)==i],scores[,whichaxes[2]][as.numeric(factraj)==i],xlim=xlim,ylim=ylim,col=coltraj[i])
par(new=T)
}
}
if(!is.null(labels)==T){textxy(scores[,whichaxes[1]],scores[,whichaxes[2]],labs=labels)}
if(!is.null(grouphull)==T){
if(is.null(colhull)==T){colhull<-c(1:nlevels(grouphull))}else{colhull<-colhull}
for(i in 1:(max(as.numeric(grouphull))))
{if(length(scores[,1][as.numeric(grouphull)==i])<3){NULL}
else{
hulli<-convexhull.xy(subset(cbind(scores[,whichaxes[1]],scores[,whichaxes[2]]),as.numeric(grouphull)==i))
par(new=T)
daplotx<-c(hulli$bdry[[1]][[1]][length(hulli$bdry[[1]][[1]])],hulli$bdry[[1]][[1]])
daploty<-c(hulli$bdry[[1]][[2]][length(hulli$bdry[[1]][[2]])],hulli$bdry[[1]][[2]])
plot(daplotx,daploty,lwd=2,type="l",lty=i,col=colhull[i],xlim=xlim,ylim=ylim,axes=F,asp=asp)
}
}}}
if(is.null(legendgroup)==F){
x11()
plot(scores[,1],scores[,2],col="white")
legend(min(scores[,1]),max(scores[,2]),unique(legendgroup), cex=1, col=unique(col), pch=unique(pch),box.col="white")}
}
#' @export
plotptau5<-function(objptau,links=NULL,projorig=T,whichaxes=c(1,2),cexscale=round(max(centroid.size(objptau$arrayord)),digits=0),shapescale=10,mag=1,subplotdim=2,shiftnegy=1,shiftposy=1,shiftnegx=1,shiftposx=1,col=as.numeric(objptau$factorord),pch=19,triang=NULL,from=0,topcs=NULL,tore=NULL,plotsource=T){
library(TeachingDemos)
k<-dim(objptau$arrayord)[1]
m<-dim(objptau$arrayord)[2]
n<-dim(objptau$arrayord)[3]
transp<-read.inn(objptau$shifted,k,m)
origsize<-objptau$indepepure
if(projorig==T){
ordiwithshapes(objptau$space1mshape,objptau$space1$x,objptau$space1$rotation,procSym=F,whichaxes=whichaxes,addata=objptau$origproj[,whichaxes],asp=1,factraj=objptau$factorord,cex=origsize/cexscale,triang=triang,from=from,to=topcs,links=links,mag=mag,shiftnegy=1.5,shiftposy=2,col=col,pch=pch,subplotdim=subplotdim,plotsource=plotsource)
title("LS on predictions on Shape Space and original data re-projected")}else{
ordiwithshapes(objptau$space1mshape,objptau$space1$x,objptau$space1$rotation,procSym=F,whichaxes=whichaxes,asp=1,factraj=objptau$factorord,cex=origsize/cexscale,triang=triang,from=from,to=topcs,links=links,mag=mag,shiftnegy=1.5,shiftposy=2,col=col,pch=pch,subplotdim=subplotdim,plotsource=plotsource)
title("LS on predictions on Shape Space")}
ranges<-apply(rbind(objptau$space1$x,objptau$origproj),2,range)
ratesorig<-ratesbygroup(read.inn(objptau$predictpure2,k,m),objptau$factorord,objptau$indepepure)
plot(objptau$indepepure,ratesorig,pch=pch,col=col)
for(i in 1:nlevels(objptau$factorord)){
lines(objptau$indepepure[as.numeric(objptau$factorord)==i][-firstsfac(objptau$factorord)],ratesorig[as.numeric(objptau$factorord)==i][-firstsfac(objptau$factorord)],col=col[firstsfac(objptau$factorord)][i])
}
title("Rates of shape change among consecutive predictions per unit size in original data")
ratestrasp<-ratesbygroup(objptau$predictpure,objptau$factorord,objptau$indepepure)
plot(objptau$indepepure,ratestrasp,pch=pch,col=col)
for(i in 1:nlevels(objptau$factorord)){
lines(objptau$indepepure[as.numeric(objptau$factorord)==i][-firstsfac(objptau$factorord)],ratestrasp[as.numeric(objptau$factorord)==i][-firstsfac(objptau$factorord)],col=col[firstsfac(objptau$factorord)][i])
}
title("Rates of shape change among consecutive predictions per unit size in transported data")
plot(ratesorig,ratestrasp,col=col,pch=pch)
title("Original rates of shape change per unit size vs those of transported data")
adults<-NULL
adultsdaplot<-NULL
for(i in 1:nlevels(objptau$factorord)){
adulti<-array(showPC2(objptau$space1$x[lastsfac(objptau$factorord),1:3][i,],objptau$space1$rotation[,1:3],objptau$space1mshape),dim=c(k,m,1))
adultdaploti<-array(showPC2(objptau$space1$x[lastsfac(objptau$factorord),1:3][i,]*mag,objptau$space1$rotation[,1:3],objptau$space1mshape),dim=c(k,m,1))
adults<-abind::abind(adults,adulti)
adultsdaplot<-abind::abind(adultsdaplot,adultdaploti)
}
if(m<3){adultsdaplot<-abind::abind(adultsdaplot,array(rep(0,k*n),dim=c(k,1,nlevels(objptau$factorord))),along=2)}
init<-array(showPC2(objptau$space1$x[firstsfac(objptau$factorord),1:3][1,],objptau$space1$rotation[,1:3],objptau$space1mshape),dim=c(k,m,1))
initdaplot<-array(showPC2(objptau$space1$x[firstsfac(objptau$factorord),1:3][1,]*mag,objptau$space1$rotation[,1:3],objptau$space1mshape),dim=c(k,m,1))
if(m<3){initdaplot<-abind::abind(initdaplot,array(rep(0,k*n),dim=c(k,1,1)),along=2)}
if(is.null(triang)==F){
adultsdaplotvis<-NULL
for(i in 1:nlevels(objptau$factorord)){
adultsdaplotvisi<-meshDist(plotsurf(adultsdaplot[,,i],triang,plot=F),plotsurf(initdaplot[,,1],triang,plot=F),to=tore,from=from,plot=F)
daagg<-objptau$space1$x[1:n,1:3][lastsfac(objptau$factorord),][i,]
adultsdaplotvisi<-translate3d(scalemesh(adultsdaplotvisi$colMesh,1/shapescale,center="none"),daagg[1],daagg[2],daagg[3])
adultsdaplotvis<-c(adultsdaplotvis,list(adultsdaplotvisi))
}
}
open3d()
for(i in 1:nlevels(objptau$factorord)){
lines3d(objptau$space1$x[1:n,1:3][as.numeric(objptau$factorord)==i,],col=i,add=T)
}
if(is.null(triang)==T){
babedaagg<-rep.row(objptau$space1$x[firstsfac(objptau$factorord),1:3][1,],k)
for(i in 1:nlevels(objptau$factorord)){
daplottarei<-adultsdaplot[,,i]
daagg<-rep.row(objptau$space1$x[1:n,1:3][lastsfac(objptau$factorord),][i,],k)
plot3D((daplottarei/shapescale)+daagg,col=i,add=T,size=1,bbox=F)
lineplot((daplottarei/shapescale)+daagg,links,col=i)
}
plot3D((initdaplot[,,1]/shapescale)+babedaagg,bbox=F,add=T)
lineplot((initdaplot[,,1]/shapescale)+babedaagg,links)
}else{
for(i in 1:nlevels(objptau$factorord)){
shade3d(adultsdaplotvis[[i]],add=T)
}
babedaagg<-rep.row(objptau$space1$x[firstsfac(objptau$factorord),1:3][1,],k)
shade3d(plotsurf((initdaplot[,,1]/shapescale)+babedaagg,triang,plot=F),add=T,alpha=0.5)
}
title3d("LS data in the Shape Space; first three PC scores")
decorate3d()
out<-list(ratesorig=ratesorig,ratestrasp=ratestrasp)
}
#' @export
pwpermancova<-function(y, pred1,pred2,nperm=999){
library(vegan)
if(is.factor(pred1)==TRUE){pred1<-factor(pred1,levels=unique(pred1))}else{NULL}
if(is.factor(pred2)==TRUE){pred2<-factor(pred2,levels=unique(pred2))}else{NULL}
if(is.factor(pred1)==TRUE){species<-as.numeric(pred1)}else{NULL}
if(is.factor(pred2)==TRUE){species<-as.numeric(pred2)}else{NULL}
if(is.factor(pred1)==TRUE){group<-pred1}else{NULL}
if(is.factor(pred2)==TRUE){group<-pred2}else{NULL}
fat_species<-group
r_adonis_pred1<-matrix(0,nrow=max(species),ncol=max(species))
p_adonis_pred1<-matrix(0,nrow=max(species),ncol=max(species))
r_adonis_pred2<-matrix(0,nrow=max(species),ncol=max(species))
p_adonis_pred2<-matrix(0,nrow=max(species),ncol=max(species))
r_adonis_pred1_pred2<-matrix(0,nrow=max(species),ncol=max(species))
p_adonis_pred1_pred2<-matrix(0,nrow=max(species),ncol=max(species))
for(i in 1:(max(species)-1))
{
for(j in (i+1):max(species)){
if(ncol(as.matrix(pred1))>1){
ado<-adonis(y[species==i|species==j,]~pred1[species==i|species==j,]*pred2[species==i|species==j],permutations=nperm,method="euclidean")
}else{NULL}
if(ncol(as.matrix(pred2))>1){
ado<-adonis(y[species==i|species==j,]~pred1[species==i|species==j]*pred2[species==i|species==j,],permutations=nperm,method="euclidean")
}else{NULL}
if(ncol(as.matrix(pred1))<2&ncol(as.matrix(pred2))<2){
ado<-adonis(y[species==i|species==j,]~pred1[species==i|species==j]*pred2[species==i|species==j],permutations=nperm,method="euclidean")
}else{NULL}
r_adonis_pred1[i,j]<-ado$aov.tab$R2[1]
p_adonis_pred1[i,j]<-ado$aov.tab$Pr[1]
r_adonis_pred2[i,j]<-ado$aov.tab$R2[2]
p_adonis_pred2[i,j]<-ado$aov.tab$Pr[2]
r_adonis_pred1_pred2[i,j]<-ado$aov.tab$R2[3]
p_adonis_pred1_pred2[i,j]<-ado$aov.tab$Pr[3]
}
}
rownames(r_adonis_pred1)=colnames(r_adonis_pred1)<-levels(group)
rownames(p_adonis_pred1)=colnames(p_adonis_pred1)<-levels(group)
rownames(r_adonis_pred2)=colnames(r_adonis_pred2)<-levels(group)
rownames(p_adonis_pred2)=colnames(p_adonis_pred2)<-levels(group)
rownames(r_adonis_pred1_pred2)=colnames(r_adonis_pred1_pred2)<-levels(group)
rownames(p_adonis_pred1_pred2)=colnames(p_adonis_pred1_pred2)<-levels(group)
out<-list(r_adonis_pred1=r_adonis_pred1,p_adonis_pred1=p_adonis_pred1,r_adonis_pred2=r_adonis_pred2,p_adonis_pred2=p_adonis_pred2,r_adonis_pred1_pred2=r_adonis_pred1_pred2,p_adonis_pred1_pred2=p_adonis_pred1_pred2)
return(out)
}
#' @export
ratesbygroup<-function(predictions,factor,indep,sss=F,plot=T){
prova3<-seqrate(predictions,factor,vector=indep,sss=sss)
prova4<-NULL
for(i in 1:nlevels(factor)){
prova4i<-c(rep(0,nlevels(factor))[i],prova3$rate[as.numeric(factor[-firstsfac(factor)])==i])
prova4<-append(prova4,prova4i)}
prova5<-prova4
if(plot==T){
plot(indep,prova5,col=as.numeric(factor),ylim=c(0,max(prova5)),xlim=range(indep))
for(i in 1:nlevels(factor)){
lines(indep[as.numeric(factor)==i][-firstsfac(factor)],prova4[as.numeric(factor)==i][-firstsfac(factor)],col=i)
}}
return(prova5)
}
#' @export
rep.row<-function(x,n){
matrix(rep(x,each=n),nrow=n)
}
#' @export
scaleshapes<-function(array){
k<-dim(array)[1]
m<-dim(array)[2]
n<-dim(array)[3]
library(abind)
library(shapes)
library(Morpho)
divmat<-function(mat){
mat2<-mat/cSize(mat)
mat2
}
if(is.matrix(array)==T){array<-array(array,dim=c(nrow(array),ncol(array),1))}
scaled<-vecx(t(apply(array,3,divmat)),revert=T,lmdim=m)
return(scaled)
}
#' @export
seqrate<-function(preds,factor,vector=centroid.size(preds),sss=F){
mytable<-table(factor)
chrate<-NULL
sizes<-NULL
seqdist<-NULL
for(i in 1: nlevels(factor)){
for(j in 2:mytable[i]){
sizei<-(vector[as.numeric(factor)==i][j])-(vector[as.numeric(factor)==i][j-1])
if(sss==T){
seqdisti<-ssriemdist2(preds[,,as.numeric(factor)==i][,,j],preds[,,as.numeric(factor)==i][,,j-1])}else{
seqdisti<-kendalldist(preds[,,as.numeric(factor)==i][,,j],preds[,,as.numeric(factor)==i][,,j-1])
}
chratei<-seqdisti/(sizei)
sizes<-append(sizes,sizei)
chrate<-append(chrate,chratei)
seqdist<-append(seqdist,seqdisti)
}
}
out<-list(rate=chrate,sizes=sizes,seqdist=seqdist)
return(out)
}
#' @export
showPC2<-function(x,rotation,mshape) {
k<-dim(mshape)[1]
m<-dim(mshape)[2]
mshape2<-array2mat(array(mshape,dim=c(k,m,1)),1,k*m)
dims <- dim(mshape)
if (length(x) > 1){predPC<-c(rotation %*% x)}else{predPC<-rotation*x}
modell <- read.inn(mshape2+predPC,k,m)[,,1]
return(modell)
}
#' @export
sigm.new<-function(mat1,mat2=NULL){
tt<-mat1
if(is.null(mat2)){yy<-mat1}else{yy<-mat2}
SGMr <- apply(tt,1,function(t)apply(yy,1,biharm.new,t))
replace(SGMr,which(SGMr=="NaN",arr.ind=T),0)
}
#' @export
tpsgridpaolo<-function (TT, YY, xbegin = -999, ybegin = -999, xlim=NULL,ylim=NULL,xwidth = -999, opt = 1, ext = 0.0, asp=1,ngrid = 22, cex = 1, pch = 20,colshift=1,zslice = 0, graphics=T,mag = 1, axes3 = FALSE,linksTT=NULL,linksYY=NULL,collinksTT=1,collinksYY=2,lwdtt=2,lwdyy=2,colgrid=1,axes2d=T,collandsTT=collinksTT,collandsYY=collinksYY,displ=T,coldispl=4){
######### some SMALL changes from Ian Dryden's function from package "shapes"
k <- dim(TT)[1]
m <- dim(TT)[2]
YY <- TT + (YY - TT) * mag
bb <- array(TT, c(dim(TT), 1))
aa <- defplotsize2(bb)
if (xwidth == -999) {
xwidth <- aa$width
}
if (xbegin == -999) {
xbegin <- aa$xl
}
if (ybegin == -999) {
ybegin <- aa$yl
}
if (m == 3) {
zup <- max(TT[, 3])
zlo <- min(TT[, 3])
zpos <- zslice
for (ii in 1:length(zslice)) {
zpos[ii] <- (zup + zlo)/2 + (zup - zlo)/2 * zslice[ii]
}
}
xstart <- xbegin
ystart <- ybegin
ngrid <- trunc(ngrid/2) * 2
kx <- ngrid
ky <- ngrid - 1
l <- kx * ky
step <- xwidth/(kx - 1)
r <- 0
X <- rep(0, times = kx)
Y2 <- rep(0, times = ky)
for (p in 1:kx) {
ystart <- ybegin
xstart <- xstart + step
for (q in 1:ky) {
ystart <- ystart + step
r <- r + 1
X[r] <- xstart
Y2[r] <- ystart
}
}
TPS <- bendingenergy(TT)
gamma11 <- TPS$gamma11
gamma21 <- TPS$gamma21
gamma31 <- TPS$gamma31
W <- gamma11 %*% YY
ta <- t(gamma21 %*% YY)
B <- gamma31 %*% YY
WtY <- t(W) %*% YY
trace <- c(0)
for (i in 1:m) {
trace <- trace + WtY[i, i]
}
if(m==2){benergy <- (16 * pi) * trace}else{benergy <- (8 * pi) * trace}
l <- kx * ky
phi <- matrix(0, l, m)
s <- matrix(0, k, 1)
for (islice in 1:length(zslice)) {
if (m == 3) {
refc <- matrix(c(X, Y2, rep(zpos[islice], times = kx *
ky)), kx * ky, m)
}
if (m == 2) {
refc <- matrix(c(X, Y2), kx * ky, m)
}
for (i in 1:l) {
s <- matrix(0, k, 1)
for (im in 1:k) {
s[im, ] <- shapes::sigmacov(refc[i, ] - TT[im, ])
}
phi[i, ] <- ta + t(B) %*% refc[i, ] + t(W) %*% s
}
if(graphics==T){
if (m == 3) {
if (opt == 2) {
shapes3d(TT, color = 2, axes3 = axes3, rglopen = FALSE)
shapes3d(YY, color = 4, rglopen = FALSE)
for (i in 1:k) {
lines3d(rbind(TT[i, ], YY[i, ]), col = 1)
}
for (j in 1:kx) {
lines3d(refc[((j - 1) * ky + 1):(ky * j), ],
color = 6)
}
for (j in 1:ky) {
lines3d(refc[(0:(kx - 1) * ky) + j, ], color = 6)
}
}
shapes3d(TT, color = collandsTT, axes3 = axes3, rglopen = FALSE)
shapes3d(YY, color = collandsYY, rglopen = FALSE)
for (i in 1:k) {
lines3d(rbind(TT[i, ], YY[i, ]), col = colshift)
}
for (j in 1:kx) {
lines3d(phi[((j - 1) * ky + 1):(ky * j), ], color = colgrid)
}
for (j in 1:ky) {
lines3d(phi[(0:(kx - 1) * ky) + j, ], color = colgrid)
}
}
}
}
if (m == 2) {
par(pty = "s")
if (opt == 2) {
par(mfrow = c(1, 2))
order <- linegrid(refc, kx, ky)
if(is.null(xlim)==T){xlim = c(xbegin -xwidth * ext, xbegin + xwidth * (1 + ext))}else{xlim=xlim}
if(is.null(ylim)==T){ylim = c(ybegin - (xwidth * ky)/kx * ext, ybegin + (xwidth * ky)/kx * (1 + ext))}else{ylim=ylim}
if(graphics==T){
plot(order[1:l, 1], order[1:l, 2], type = "l", xlim = xlim, ylim = ylim, xlab = " ", ylab = " ",asp=asp,,axes=axes2d)
lines(order[(l + 1):(2 * l), 1], order[(l + 1):(2 * l), 2], type = "l",col=colgrid,xlim = xlim, ylim = ylim,asp=asp)
points(TT, cex = cex, pch = pch, col = collandsTT)
}}
if(is.null(xlim)==T){xlim = c(xbegin -xwidth * ext, xbegin + xwidth * (1 + ext))}else{xlim=xlim}
if(is.null(ylim)==T){ylim = c(ybegin - (xwidth * ky)/kx * ext, ybegin + (xwidth * ky)/kx * (1 + ext))}else{ylim=ylim}
order <- linegrid(phi, kx, ky)
if(graphics==T){plot(order[1:l, 1], order[1:l, 2], type = "l", xlim = xlim, ylim = ylim, xlab = " ", ylab = " ",col=colgrid,asp=asp,axes=axes2d)
lines(order[(l + 1):(2 * l), 1], order[(l + 1):(2 * l), 2], type = "l",col=colgrid,xlim = xlim, ylim = ylim,asp=asp)}
if(graphics==T){points(YY, cex = cex, pch = pch, col = collandsYY)}
if(graphics==T){points(TT, cex = cex, pch = pch, col = collandsTT)}
if(graphics==T){
if(displ==T){
for (i in 1:(k)) {
arrows(TT[i, 1], TT[i, 2], YY[i, 1], YY[i, 2], col = coldispl, length = 0.1, angle = 20)
}
}}
firstcol<-order[1:l,][1:kx,][-kx,]
firstrow<-order[(l + 1):(2 * l),][((kx*ky-1)-ky):(kx*ky-1),][-c(1:2),]
lastcol<-order[1:l,][((kx*ky)-ky):(kx*ky),][-1,]
lastrow<-order[(l + 1):(2 * l),][2:ky,][order((ky:2)),]
bound = rbind(firstcol,firstrow,lastcol,lastrow)
them<-nrow(order[1:l,])/ngrid
}
if(m==2){grid<-list(ngrid=order[1:l,],ngrid2=order,m=them,n=ngrid,bound=bound,l=l,xlim=xlim,ylim=ylim,TT=TT,YY=YY)}else{
grid<-list(n=ngrid,l=l,TT=TT,YY=YY)
}
out<-list(YY=YY,gamma11=gamma11,gamma21=gamma21,gamma31=gamma31,W=W,ta=ta,B=B,WtY=WtY,trace=trace,benergy=benergy,grid=grid,kx=kx,ky=ky)
if(graphics==T){if(!is.null(linksTT)){lineplot(TT,linksTT,lwd=lwdtt,col=collinksTT)}}
if(graphics==T){if(!is.null(linksYY)){lineplot(YY,linksYY,lwd=lwdyy,col=collinksYY)}}
return(out)
}
#' @export
replotptau5<-function(objptau,objplotptau5,projorig=T,whichaxes=c(1,2),cexscale=round(max(centroid.size(objptau$arrayord)),digits=0),shapescale=10,mag=1,subplotdim=2,shiftnegy=1,shiftposy=1,shiftnegx=1,shiftposx=1,col=as.numeric(objptau$factorord),pch=as.numeric(objptau$factorord),triang=NULL,from=NULL,scaleramp=F,topcs=NULL,tore=NULL,plotsource=T){
library(TeachingDemos)
links=objplotptau5$links
triang<-objplotptau5$triang
k<-dim(objptau$arrayord)[1]
m<-dim(objptau$arrayord)[2]
n<-dim(objptau$arrayord)[3]
transp<-read.inn(objptau$shifted,k,m)
origsize<-objptau$indepepure
if(projorig==T){
ordiwithshapes(objptau$space1mshape,objptau$space1$x,objptau$space1$rotation,procSym=F,whichaxes=whichaxes,addata=objptau$origproj[,whichaxes],asp=1,factraj=objptau$factorord,cex=origsize/cexscale,triang=triang,from=from,to=topcs,links=links,mag=mag,shiftnegy=1.5,shiftposy=2,col=col,pch=pch,subplotdim=subplotdim,plotsource=plotsource)
title("LS on predictions in PCA space and original data re-projected")}else{
ordiwithshapes(objptau$space1mshape,objptau$space1$x,objptau$space1$rotation,procSym=F,whichaxes=whichaxes,asp=1,factraj=objptau$factorord,cex=origsize/cexscale,triang=triang,from=from,to=topcs,links=links,mag=mag,shiftnegy=1.5,shiftposy=2,col=col,pch=pch,subplotdim=subplotdim,plotsource=plotsource)
title("LS on predictions in PCA space")}
ranges<-apply(rbind(objptau$space1$x,objptau$origproj),2,range)
ratesorig<-ratesbygroup(objptau$predictpure,objptau$factorord,objptau$indepepure,plot=F)
ratestrasp<-ratesbygroup(read.inn(objptau$predictpure2,k,m),objptau$factorord,objptau$indepepure,plot=F)
plot(objptau$indepepure,ratesorig,pch=pch,col=col)
for(i in 1:nlevels(objptau$factorord)){
lines(objptau$indepepure[as.numeric(objptau$factorord)==i][-firstsfac(objptau$factorord)],ratesorig[as.numeric(objptau$factorord)==i][-firstsfac(objptau$factorord)],col=col[firstsfac(objptau$factorord)][i])
}
title("Rates of shape change among consecutive predictions per unit size in original data")
if(is.null(triang)){open3d()}
if(is.null(triang)){heat=F}else{heat=T}
rgl.close()
if(is.null(triang)){triang2<-NULL}else{triang2<-t(triang)}
mydef<-plotdefo3d(list(mshape=objptau$space1mshape,PCs=objptau$space1$rotation,PCscores=objptau$space1$x),procSym=F,triang=triang2,from=from,to=topcs,links=links,heat=heat)
plot(objptau$indepepure,ratestrasp,pch=pch,col=col)
for(i in 1:nlevels(objptau$factorord)){
lines(objptau$indepepure[as.numeric(objptau$factorord)==i][-firstsfac(objptau$factorord)],ratestrasp[as.numeric(objptau$factorord)==i][-firstsfac(objptau$factorord)],col=col[firstsfac(objptau$factorord)][i])
}
title("Rates of shape change among consecutive predictions per unit size in transported data")
plot(ratesorig,ratestrasp,col=col,pch=pch)
title("Original rates of shape change per unit size vs those of transported data")
adults<-NULL
adultsdaplot<-objplotptau5$adultsdaplot
init<-objplotptau5$init
initdaplot<-objplotptau5$initdaplot
adultsdaplotvis<-objplotptau5$adultsdaplotvis
alldists<-objplotptau5$alldists
open3d()
r3dDefaults$windowRect <- c(0,50, 8000, 8000)
mfrow3d(1,2,byrow=T,sharedMouse=T)
for(i in 1:nlevels(objptau$factorord)){
lines3d(objptau$space1$x[1:n,1:3][as.numeric(objptau$factorord)==i,],col=i,add=T)
}
title3d("LS data in the PCA; first three PC scores")
decorate3d(xlab="PC1",ylab="PC2",zlab="PC3")
if(projorig==T){
plot3D(objptau$origproj[,1:3],col=as.numeric(objptau$factorord),add=T,bbox=F,cex=objptau$indepepure/cexscale)
}
if(heat==F){
babedaagg<-rep.row(objptau$space1$x[firstsfac(objptau$factorord),1:3][1,],k)
for(i in 1:nlevels(objptau$factorord)){
daplottarei<-adultsdaplot[,,i]
daagg<-rep.row(objptau$space1$x[1:n,1:3][lastsfac(objptau$factorord),][i,],k)
plot3D((daplottarei/shapescale)+daagg,col=i,add=T,size=1,bbox=F)
if(!is.null(links)){lineplot((daplottarei/shapescale)+daagg,links,col=i)}
if(!is.null(triang)){plotsurf((daplottarei/shapescale)+daagg,triang,col=i)}
}
plot3D((initdaplot[,,1]/shapescale)+babedaagg,bbox=F,add=T)
if(!is.null(links)){lineplot((initdaplot[,,1]/shapescale)+babedaagg,links)}
if(!is.null(triang)){plotsurf((initdaplot[,,1]/shapescale)+babedaagg,triang,col=1)}
}else{
for(i in 1:nlevels(objptau$factorord)){
meshDist(adultsdaplotvis[[i]],distvec=alldists[[i]],from=range(alldists)[1],to=range(alldists)[2],add=T,scaleramp=scaleramp)
}
babedaagg<-rep.row(objptau$space1$x[firstsfac(objptau$factorord),1:3][1,],k)
shade3d(plotsurf((initdaplot[,,1]/shapescale)+babedaagg,triang,plot=F),add=T,alpha=0.5)
}
next3d()
babedaagg<-objptau$space1$x[firstsfac(objptau$factorord),1:3][1,]
for(i in 1:nlevels(objptau$factorord)){
daagg<-rep.row(objptau$space1$x[1:n,1:3][lastsfac(objptau$factorord),][i,],k)
text3d(daagg,texts=levels(objptau$factorord)[i],col=i,add=T,size=1)
}
text3d(babedaagg,texts="CR",bbox=F,add=T)
decorate3d(xlab="PC1",ylab="PC2",zlab="PC3")
for(i in 1:nlevels(objptau$factorord)){
lines3d(objptau$space1$x[1:n,1:3][as.numeric(objptau$factorord)==i,],col=i,add=T)
}
title3d("LS data in the PCA; first three PC scores")
}
#' @export
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.