R/auxilliary.functions.R

In 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

  xanc <- ape::ace(xx, phy, method=method)$ace
  xall = c(xx,xanc)

  a0 = max(ages)-ages[phy$edge[,1]]
  a1 = max(ages)-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(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)
  }

  plot.args<-list(...)

  if(!"xlab"%in%names(plot.args)) plot.args$xlab<-""
  if(!"ylab"%in%names(plot.args)) plot.args$ylab<-""
  if(!"bty"%in%names(plot.args)) plot.args$bty<-"n"
  if(!"cex.axis"%in%names(plot.args)) plot.args$cex.axis<-0.8
  if(!"yaxt"%in%names(plot.args)) plot.args$yaxt<-"n"
  if(!"lwd"%in%names(plot.args)) plot.args$lwd<-2

  do.call(plot,c(list(x=NA,xlim=rev(range(c(a0,a1))),ylim=range(c(x0,x1))),plot.args))
  do.call(segments,c(list(x0=a0,y0=x0,x1=a1,y1=x1,col=colo),
                     plot.args[which(names(plot.args)%in%c("lwd","lty","lend","ljoin","lmitre"))]))


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

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

Plot_ConvexHull<-function(xcoord, ycoord,...){
  hpts <- chull(x = xcoord, y = ycoord)
  hpts <- c(hpts, hpts[1])

  poly.args<-list(...)
  do.call(polygon,c(list(x=xcoord[hpts],y=ycoord[hpts]),poly.args))
}

dosur <- function(scores,pcs,sel=NULL,mshape,radius=0){
  if(is.null(sel)==TRUE) {
    temp<-Morpho::restoreShapes(scores,pcs,mshape)
  } else {
    temp<-Morpho::restoreShapes(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))
}