R/auxilliary.functions.R

In pasraia/RRphylo: Phylogenetic Ridge Regression Methods for Comparative Studies

range01 <- function(x, ...){(x - min(x, ...)) / (max(x, ...) - min(x, ...))}

unitV <- function(x) sum(x^2)^0.5

deg2rad <- function(deg) (deg * pi)/(180)

rad2deg <- function(rad)  (rad * 180)/(pi)

traitgram = function(
  x, phy,
  xaxt='s',
  underscore = FALSE,
  show.names = TRUE,
  show.xaxis.values = TRUE,
  method = c('ML','pic'),
  col=NULL,
  lwd=NULL,
  mgp=NULL,...)
{

  method <- match.arg(method)
  Ntaxa = length(phy$tip.label)
  Ntot = Ntaxa + phy$Nnode
  phy = picante::node.age(phy)
  ages = phy$ages[match(1:Ntot,phy$edge[,2])]
  ages[Ntaxa+1]=0


  if (class(x) %in% c('matrix','array')) {
    xx = as.numeric(x)
    names(xx) = row.names(x)
  } else xx = x

  if (!is.null(names(xx))) {
    umar = 0.1
    if (!all(names(xx) %in% phy$tip.label)) {
      print('trait and phy names do not match')
      return()
    }
    xx = xx[match(phy$tip.label,names(xx))]
  } else umar = 0.1

  lmar = 0.2
  if (xaxt=='s') if (show.xaxis.values) lmar = 1 else lmar = 0.5
  xanc <- ape::ace(xx, phy, method=method)$ace
  xall = c(xx,xanc)

  a0 = ages[phy$edge[,1]]
  a1 = ages[phy$edge[,2]]
  x0 = xall[phy$edge[,1]]
  x1 = xall[phy$edge[,2]]


  if (show.names) {
    maxNameLength = max(nchar(names(xx)))
    ylim = c(min(ages),max(ages)*(1+maxNameLength/50))
    if (!underscore) names(xx) = gsub('_',' ',names(xx))
  } else ylim = range(ages)

  if(is.null(mgp)) magp<-NULL else{
    mgp->magp
  }

  if(is.null(col)) colo<-par("fg") else{
    col->colo
    colo[match(names(x1),names(colo))]->colo
    data.frame(x0,x1,a0,a1,colo)->dato
    dato[order(dato[,5],decreasing=TRUE),]->dato
    dato[,1]->x0
    dato[,2]->x1
    dato[,3]->a0
    dato[,4]->a1
    as.character(dato[,5])->colo
    rownames(dato)->names(x1)->names(x0)->names(a1)->names(a0)->names(colo)
  }

  par(mar = c(3, 2.5, 2, 1))
  plot(range(c(a0,a1)),range(c(x0,x1)),

       type='n',xaxt='n',yaxt='n',
       xlab='',ylab='',bty='n',cex.axis=0.8)
  if (xaxt=='s') if (show.xaxis.values) axis(1,labels=TRUE,mgp=magp) else axis(1,labels=FALSE,mgp=magp)

  if(is.null(lwd)) linwd<-1 else{
    lwd->linwd
    linwd[match(names(x1),names(linwd))]->linwd
  }

  segments(a0,x0,a1,x1,col=colo,lwd=linwd)

  if (show.names) {
    text(max(ages),sort(xx),
         labels = names(xx)[order(xx)],
         adj = -0,
         srt=90,
         cex=.3)
  }

  return(data.frame(a1,x1))
}

Plot_ConvexHull<-function(xcoord, ycoord, lcolor,lwd=NULL, lty=NULL,col.p=NULL){
  hpts <- chull(x = xcoord, y = ycoord)
  hpts <- c(hpts, hpts[1])
  lines(xcoord[hpts], ycoord[hpts], col = lcolor,lwd=lwd, lty=lty)
  polygon(xcoord[hpts], ycoord[hpts], col=col.p, border=NA)
}

dosur <- function(scores,pcs,sel=NULL,mshape,radius=0){
  if(is.null(sel)==TRUE) {
    temp<-Morpho::showPC(scores,pcs,mshape)
  } else {
    temp<-Morpho::showPC(scores[sel],pcs[,sel],mshape)
  }
  mshape<-Rvcg::vcgBallPivoting(mshape, radius = radius)
  sur<-mshape
  sur$vb[1:3,]<-t(temp)
  return(sur)

}

areadiff<-function(mesh1,mesh2,out.rem=FALSE,scale01=TRUE,fact=1.5){

  area_shape1<-Rvcg::vcgArea(mesh1,perface=T)$pertriangle
  area_shape2<-Rvcg::vcgArea(mesh2,perface=T)$pertriangle
  diff_areas<-(area_shape1-area_shape2)/area_shape1
  sel<-which(is.na(diff_areas))

  if(length(sel)>0){
    mesh1$it<-mesh1$it[,-sel]
    mesh2$it<-mesh2$it[,-sel]
    mesh1<-Morpho::rmUnrefVertex(mesh1)
    mesh2<-Morpho::rmUnrefVertex(mesh2)
    area_shape1<-Rvcg::vcgArea(mesh1,perface=T)$pertriangle
    area_shape2<-Rvcg::vcgArea(mesh2,perface=T)$pertriangle
    diff_areas<-(area_shape1-area_shape2)/area_shape1
  }

  if(out.rem==TRUE){
    x=diff_areas
    qq <- quantile(x, c(1,3)/4, names=FALSE)
    r <- diff(qq) * fact
    tst <- x < qq[1] - r | x > qq[2] + r
    tstp<-qq[2] + r
    tstn<-qq[1] - r
    diff_areas[x>tstp]<-tstp
    diff_areas[x<tstn]<-tstn
  }else diff_areas=diff_areas

  if(scale01==TRUE) diff_areas<-range01(diff_areas)

  return(list("ash1"=area_shape1,"ash2"=area_shape2,"dareas"=diff_areas))
}

localmeshdiff <- function(mesh1, mesh2, ploton, paltot = rainbow(200),
                          from = NULL, to = NULL, n.int = 200, out.rem = TRUE, fact = 1.5,
                          visual = 1, scale01 = TRUE, vec = NULL, plot = FALSE,
                          densityplot = TRUE) {


  if(is.null(vec)) {
    areadiff(mesh1,mesh2,out.rem=out.rem,scale01=scale01,fact=fact)->adiff
    area_shape1<-adiff$ash1
    area_shape2<-adiff$ash2
    diff_areas<-adiff$dareas
    } else diff_areas<-vec

  cat("the range of diff_areas is ", range(diff_areas),
      sep = "\n")
  if (is.null(to) == TRUE) {
    to <- max(diff_areas) * 1.01
  }
  if (is.null(from) == TRUE) {
    from <- min(diff_areas) * 1.01
  }
  selfromto <- which(diff_areas < to & diff_areas >= from)

  colmap_tot <- colorRampPalette(paltot)
  breaks_tot <- cut(c(from, diff_areas, to), n.int)
  cols_tot <- colmap_tot(n.int)[breaks_tot]
  cols_tot <- cols_tot[-c(1, length(cols_tot))]
  selfromto <- which(diff_areas < to & diff_areas >= from)
  cols_tot[-selfromto] <- "#FFFFFF"

  if (isTRUE(densityplot)) {
    plot(density(c(from, diff_areas, to)), main = "",
         xlab = "", ylab = "")
    abline(v = seq(from, to, length.out = n.int), col = colmap_tot(n.int),
           lwd = 5)
    points(density(diff_areas), type = "l", lwd = 2)
  }
  if (ploton == 1) {
    meshtobeplotted <- mesh1
  }
  if (ploton == 2) {
    meshtobeplotted <- mesh2
  }
  if (plot == TRUE) {
    if (visual == 1) {
      rgl::triangles3d(t(meshtobeplotted$vb[, meshtobeplotted$it]),
                       col = rep(cols_tot, each = 3), alpha = 1, lit = TRUE,
                       specular = "black")
    }
    if (visual == 2) {
      rgl::triangles3d(t(meshtobeplotted$vb[, meshtobeplotted$it]),
                       col = rep(cols_tot, each = 3), alpha = 1, lit = TRUE,
                       specular = "black")
      rgl::wire3d(meshtobeplotted, col = "blue",
                  lit = FALSE)
    }
  }
  diffvert <- NULL
  track_v <- Rvcg::vcgVFadj(meshtobeplotted)
  for (i in 1:length(track_v)) {
    diffvert[i] <- mean(diff_areas[track_v[[i]]])
  }
  selfromto <- which(diffvert < to & diffvert >= from)
  colmap_tot <- colorRampPalette(paltot)
  breaks_tot <- cut(c(from, diffvert, to), n.int)
  cols_tot <- colmap_tot(n.int)[breaks_tot]
  cols_tot <- cols_tot[-c(1, length(cols_tot))]
  selfromto <- which(diffvert < to & diffvert >= from)
  cols_tot[-selfromto] <- "#FFFFFF"
  mesh <- meshtobeplotted
  mesh$material$color <- cols_tot
  if (!is.null(vec)) list(mesh = mesh) else
    return(list(ash1 = area_shape1, ash2 = area_shape2,
                dareas = diff_areas, mesh = mesh))
}