do_outline_alpha <- function(rp, alpha)
{
ah = alphahull::ashape(rp,alpha=alpha)
return(ah)
}
#do_outline_concave <- function(rp, concavity)
#{
# cm = concaveman::concaveman(rp,concavity=concavity, length_threshold=0)
# return(cm)
#}
do_outline_ball <- function(rp, radius)
{
gb = terra::buffer(sp::SpatialPoints(rp), quadsegs=2, width=radius)
return(gb)
}
do_outline_raster <- function(pts,res)
{
pts <- as.matrix(pts)
pr <- padded_range(pts,multiply.interval.amount=0.25)
e <- terra::ext(pr[1,1], pr[2,1], pr[1,2], pr[2,2])
r <- terra::rast(e, ncols=res, nrows=res)
x <- terra::rasterize(pts, r, rep(1, nrow(pts)), fun=mean,background=NA)
w <- terra::as.polygons(x,dissolve=TRUE)
return(w)
}
plot.Hypervolume <- function(x, ...)
{
templist = new("HypervolumeList")
templist@HVList=list(x)
plot.HypervolumeList(templist, ...)
}
extendrange <- function(x,factor=0.5)
{
xmin <- min(x,na.rm=TRUE)
xmax <- max(x,na.rm=TRUE)
xminf <- xmin - (xmax - xmin)*factor
xmaxf <- xmax + (xmax - xmin)*factor
result <- c(xminf, xmaxf)
return(result)
}
plot.HypervolumeList <- function(x,
show.3d=FALSE,plot.3d.axes.id=NULL,
show.axes=TRUE, show.frame=TRUE,
show.random=TRUE, show.density=TRUE,show.data=TRUE,
names=NULL, show.legend=TRUE, limits=NULL,
show.contour=TRUE, contour.lwd=1.5,
contour.type='kde',
contour.alphahull.alpha=0.25,
contour.ball.radius.factor=1,
contour.kde.level=1e-4,
contour.raster.resolution=20,
show.centroid=TRUE, cex.centroid=2,
colors=rainbow(floor(length(x@HVList)*1.5),alpha=0.8),
point.alpha.min=0.2, point.dark.factor=0.5,
cex.random=0.5,cex.data=0.75,cex.axis=0.75,cex.names=1.0,cex.legend=0.75,
num.points.max.data = 1000, num.points.max.random = 2000, reshuffle=TRUE,
plot.function.additional=NULL,
verbose=FALSE,
...)
{
#### ALEX !!!!!! check the method, if n_occupancy remove any 0
method_is_occupancy <- FALSE
if(inherits(x, "Hypervolume"))
{
# with the new Method n_overlap and n_overlap_test it makes sense to calculate
# weighted mean
if(identical(x@Method, "n_occupancy") | identical(x@Method, "n_occupancy_test") | identical(x@Method, "n_occupancy_permute") |
identical(x@Method, "occupancy_to_union") | identical(x@Method, "occupancy_to_unshared") | identical(x@Method, "occupancy_to_intersection")){
method_is_occupancy <- TRUE
}
if(nrow(x@RandomPoints) == 0){
stop("The hypervolume has no random points.")
}
}
if(inherits(x, "HypervolumeList")){
method_list <- unique(unlist(lapply(x@HVList, function(x) x@Method)))
if(identical(method_list, "n_occupancy") | identical(method_list, "n_occupancy_test") | identical(method_list, "n_occupancy_permute") |
identical(method_list, "occupancy_to_union") | identical(method_list, "occupancy_to_unshared") | identical(method_list, "occupancy_to_intersection")){
method_is_occupancy <- TRUE
}
rp_list <- sapply(x@HVList, function(x) nrow(x@RandomPoints))
if(all(rp_list == 0)){
stop("All the hypervolumes have no random points.")
} else if(any(rp_list == 0)){
x <- x[[which(rp_list > 0)]]
names_list <- sapply(x@HVList, function(x) x@Name)
names_list <- names_list[rp_list > 0]
names_list <- paste("Hypervolumes", paste(names_list, collapse = ", "), "were removed because had no random points")
message(names_list)
}
}
#### ALEX !!!!!!
# remove 0 values from ValueAtRandomPoints
# remove the corresponding coordinates
if(method_is_occupancy){
if(inherits(x, "HypervolumeList")){
for(i in 1:length(x@HVList)){
hv_temp <- x@HVList[[i]]
x@HVList[[i]]@RandomPoints <- hv_temp@RandomPoints[! is.na(hv_temp@ValueAtRandomPoints), ]
x@HVList[[i]]@ValueAtRandomPoints <- hv_temp@ValueAtRandomPoints[! is.na(hv_temp@ValueAtRandomPoints)]
hv_temp <- x@HVList[[i]]
x@HVList[[i]]@RandomPoints <- hv_temp@RandomPoints[hv_temp@ValueAtRandomPoints != 0, ]
x@HVList[[i]]@ValueAtRandomPoints <- hv_temp@ValueAtRandomPoints[hv_temp@ValueAtRandomPoints != 0]
}
}
}
#### ALEX !!!!!!
# remove columns that are not to be plotted, put here but it will be useful later
if(method_is_occupancy){
columns_to_remove <- 3
} else {
columns_to_remove <- 2
}
sapply(x@HVList, function(z)
{
if (verbose==TRUE)
{
cat(sprintf("Showing %d random points of %d for %s\n",min(nrow(z@RandomPoints), num.points.max.random), nrow(z@RandomPoints), z@Name))
}
if (show.data && length(z@Data) > 0)
{
npd <- ifelse(all(is.nan(z@Data)), 0, nrow(z@Data))
if (verbose==TRUE)
{
cat(sprintf("Showing %d data points of %d for %s\n",min(num.points.max.data, npd), npd, z@Name))
}
}
})
if (!requireNamespace("alphahull", quietly = TRUE)) {
warning("The package 'alphahull' is needed for contour plotting with contour.type='alphahull'. Please install it to continue.\n\n *** Temporarily setting contour.type='kde'.", call. = FALSE)
contour.type <- 'kde'
}
alldims = sapply(x@HVList, function(z) { z@Dimensionality })
allnames = sapply(x@HVList, function(z) { z@Name })
stopifnot(all(alldims[1] == alldims))
######################################################################################################
### ALEX if-else, if the method is occupancy or occupancy test adds the occupancy column
# otherwise keep it normal
if(method_is_occupancy){
all <- NULL
alldata <- NULL
for (i in 1:length(x@HVList))
{
ivals = sample(nrow(x@HVList[[i]]@RandomPoints), min(c(num.points.max.random, nrow(x@HVList[[i]]@RandomPoints))))
#### ALEX !!!!!!
# include the column "occupancy", it will be needed later for setting cex
subsampledpoints = data.frame(x@HVList[[i]]@RandomPoints[ivals,,drop=FALSE])
densityvals = x@HVList[[i]]@ValueAtRandomPoints[ivals]
if (nrow(subsampledpoints) > 0)
{
subsampledpoints = cbind(subsampledpoints, ID=rep(i, nrow(subsampledpoints)), Density=(densityvals-min(densityvals,na.rm=TRUE))/(max(densityvals,na.rm=TRUE)-min(densityvals,na.rm=TRUE)), Occupancy = abs(x@HVList[[i]]@ValueAtRandomPoints[ivals]))
subsampledpoints[is.nan(subsampledpoints[,"Density"]),"Density"] <- 1
all <- rbind(all, subsampledpoints)
}
thisdata=x@HVList[[i]]@Data
alldata <- rbind(alldata, cbind(thisdata, ID=rep(i,nrow(thisdata))))
}
} else {
all <- NULL
alldata <- NULL
for (i in 1:length(x@HVList))
{
ivals = sample(nrow(x@HVList[[i]]@RandomPoints), min(c(num.points.max.random, nrow(x@HVList[[i]]@RandomPoints))))
subsampledpoints = data.frame(x@HVList[[i]]@RandomPoints[ivals,,drop=FALSE])
densityvals = x@HVList[[i]]@ValueAtRandomPoints[ivals]
if (nrow(subsampledpoints) > 0)
{
subsampledpoints = cbind(subsampledpoints, ID=rep(i, nrow(subsampledpoints)), Density=(densityvals-min(densityvals,na.rm=TRUE))/(max(densityvals,na.rm=TRUE)-min(densityvals,na.rm=TRUE)))
subsampledpoints[is.nan(subsampledpoints[,"Density"]),"Density"] <- 1
all <- rbind(all, subsampledpoints)
}
thisdata=x@HVList[[i]]@Data
alldata <- rbind(alldata, cbind(thisdata, ID=rep(i,nrow(thisdata))))
}
}
######################################################################################################
alldata <- as.data.frame(alldata)
if (num.points.max.data < nrow(alldata) && !is.null(num.points.max.data))
{
alldata <- alldata[sample(nrow(alldata), min(c(num.points.max.data, nrow(alldata)))),]
}
if (is.null(all))
{
warning('No random points to plot.')
if (is.null(dimnames(x@HVList[[1]]@RandomPoints)[[2]]))
{
all <- matrix(0,ncol=2+alldims,nrow=1,dimnames=list(NULL,c(paste("X",1:alldims,sep=""),"ID","Density")))
}
else
{
all <- matrix(0,ncol=2+alldims,nrow=1,dimnames=list(NULL,c(dimnames(x@HVList[[1]]@RandomPoints)[[2]],"ID","Density")))
}
all <- as.data.frame(all)
}
if (reshuffle==TRUE)
{
all <- all[sample(nrow(all),replace=FALSE),,drop=FALSE] # reorder to shuffle colors
alldata <- alldata[sample(nrow(alldata),replace=FALSE),,drop=FALSE]
}
no_names_supplied = FALSE
if (is.null(names))
{
dn = dimnames(all)[[2]]
names = dn[1:(ncol(all)-columns_to_remove)]
no_names_supplied = TRUE
}
if (!is.null(limits) & !is.list(limits))
{
varlimlist = vector('list',ncol(all)-2)
for (i in 1:length(varlimlist))
{
varlimlist[[i]] <- limits
}
limits = varlimlist
}
colorlist <- colors[all$ID]
alphavals <- (all$Density - quantile(all$Density, 0.025, na.rm=T)) / (quantile(all$Density, 0.975, na.rm=T) - quantile(all$Density,0.025, na.rm=T))
alphavals[is.nan(alphavals)] <- 0.5 # in case the quantile is un-informative
alphavals[alphavals < 0] <- 0
alphavals[alphavals > 1] <- 1
alphavals <- point.alpha.min + (1 - point.alpha.min)*alphavals
if (show.density==FALSE)
{
alphavals <- rep(1, length(colorlist))
}
for (i in 1:length(colorlist))
{
colorlist[i] <- rgb_2_rgba(colorlist[i], alphavals[i])
}
colorlistdata = colors[alldata$ID]
for (i in 1:length(colorlistdata))
{
colorlistdata[i] <- rgb_2_set_hsv(colorlistdata[i], v=1-point.dark.factor)
}
if (ncol(all) < 2)
{
stop('Plotting only available in n>=2 dimensions.')
}
if (show.3d==FALSE)
{
op = par(no.readonly = T)
par(mfrow=c(ncol(all)-columns_to_remove, ncol(all)-columns_to_remove))
par(mar=c(0,0,0,0))
par(oma=c(0.5,0.5,0.5,0.5))
for (i in 1:(ncol(all)-columns_to_remove))
{
for (j in 1:(ncol(all)-columns_to_remove))
{
if (j > i)
{
# set up axes with right limits
plot(all[,j], all[,i],type="n",axes=FALSE,xlim=limits[[j]], ylim=limits[[i]],bty='n')
# draw random points
if(show.random==TRUE)
{
if(method_is_occupancy){
#### ALEX !!!!!!
# here I set the cex
cex.occupancy <- all[, "Occupancy"]
points(all[,j], all[,i], col=colorlist, cex= cex.occupancy / max(cex.occupancy) * cex.random, pch = 16)
} else {
points(all[,j], all[,i], col=colorlist,cex=cex.random,pch=16)
}
}
# show data
if (show.data & nrow(alldata) > 0)
{
points(alldata[,j], alldata[,i], col=colorlistdata,cex=cex.data,pch=16)
}
if (show.centroid == TRUE)
{
for (whichid in 1:length(unique(all$ID)))
{
allss <- subset(all, all$ID==whichid)
if(method_is_occupancy){
#### ALEX !!!!!!
# weighted centroids
centroid_x <- weighted.mean(allss[,j], cex.occupancy[all$ID==whichid], na.rm=TRUE)
centroid_y <- weighted.mean(allss[,i], cex.occupancy[all$ID==whichid], na.rm=TRUE)
} else{
centroid_x <- mean(allss[,j],na.rm=TRUE)
centroid_y <- mean(allss[,i],na.rm=TRUE)
}
# draw point
points(centroid_x, centroid_y, col=colors[whichid],cex=cex.centroid,pch=16)
# add a white boundary for clarity
points(centroid_x, centroid_y, col='white',cex=cex.centroid,pch=1,lwd=1.5)
}
}
# calculate contours
if (show.contour==TRUE)
{
# draw shaded centers
for (whichid in 1:length(unique(all$ID)))
{
allss <- subset(all, all$ID==whichid) # remove oversampling of some points
if (nrow(allss) > 0)
{
contourx <- allss[,j]
contoury <- allss[,i]
rp = cbind(contourx, contoury)
vol_this = x@HVList[[whichid]]@Volume
density_this = nrow(rp) / vol_this
dim_this = x@HVList[[whichid]]@Dimensionality
radius_critical <- density_this^(-1/dim_this) * contour.ball.radius.factor
if (contour.type=='alphahull')
{
poly_outline = do_outline_alpha(rp=rp, alpha=contour.alphahull.alpha)
plot(poly_outline,add=TRUE,wpoints=FALSE,wlines='none',lwd=contour.lwd,col=colors[whichid])
}
else if (contour.type=='ball')
{
poly_outline <- do_outline_ball(rp=rp, radius=radius_critical)
sp::plot(poly_outline, add=TRUE,lwd=contour.lwd,col=colors[whichid])
}
else if (contour.type=='kde')
{
if (nrow(rp) > 1)
{
m_kde = kde2d(rp[,1], rp[,2], n=50, h=radius_critical)
contour(m_kde, add=TRUE, levels=contour.kde.level,drawlabels=FALSE,lwd=contour.lwd,col=colors[whichid])
}
}
else if (contour.type=='raster')
{
poly_raster <- do_outline_raster(as.matrix(rp),res=contour.raster.resolution)
sp::plot(poly_raster, add=TRUE, lwd=contour.lwd,col=colors[whichid])
}
}
}
}
if (!is.null(plot.function.additional))
{
plot.function.additional(j,i)
}
if (show.frame==TRUE)
{
box()
}
}
else if (j == i)
{
plot(0,0,type="n",xlim=c(0,1),ylim=c(0,1),axes=FALSE)
text(0.5, 0.5, names[j],cex=cex.names)
}
else if (j==1 & i == (ncol(all) - columns_to_remove))
{
plot(0,0,type="n",xlim=c(0,1),ylim=c(0,1),axes=FALSE)
if (show.legend == TRUE)
{
legend('topleft',legend=allnames,text.col=colors,bty='n',cex=cex.legend)
}
}
else
{
plot(0,0,type="n",axes=FALSE)
}
if (j==i+1)
{
if (show.axes==TRUE)
{
axis(side=1,cex.axis=cex.axis)
axis(side=2,cex.axis=cex.axis)
}
}
}
}
par(op)
}
else
{
if (is.null(plot.3d.axes.id))
{
plot.3d.axes.id=1:3
}
if (no_names_supplied==TRUE)
{
axesnames <- names[plot.3d.axes.id]
}
else
{
axesnames <- names
}
if(length(plot.3d.axes.id)!=3) { stop('Must specify three axes') }
if (show.density==TRUE)
{
for (i in 1:length(colorlist))
{
colorlist[i] <- rgb_2_set_hsv(colorlist[i], s=(alphavals[i]^2)) # should do this with alpha, but a workaround for non-transparent OpenGL implementations...
}
}
rgl::plot3d(all[,plot.3d.axes.id],col=colorlist,expand=1.05, xlab=axesnames[1], ylab=axesnames[2], zlab=axesnames[3], xlim=limits[[1]],ylim=limits[[2]],zlim=limits[[3]],size=cex.random,type='p',box=show.frame,axes=show.axes)
if (show.legend==TRUE)
{
for (i in 1:length(allnames))
{
rgl::mtext3d(allnames[i],edge='x-+',line=1+i*cex.legend*1.25,color=colors[i],cex=cex.legend)
}
}
if (show.data)
{
if (!any(is.nan(as.matrix(alldata[,plot.3d.axes.id]))))
{
rgl::points3d(x=alldata[,plot.3d.axes.id[1]], y=alldata[,plot.3d.axes.id[2]], z=alldata[,plot.3d.axes.id[3]], col=colorlistdata,cex=cex.data,pch=16)
}
}
if (show.centroid == TRUE)
{
for (whichid in 1:length(unique(all$ID)))
{
allss <- subset(all, all$ID==whichid)
centroid_1 <- mean(allss[,plot.3d.axes.id[1]],na.rm=TRUE)
centroid_2 <- mean(allss[,plot.3d.axes.id[2]],na.rm=TRUE)
centroid_3 <- mean(allss[,plot.3d.axes.id[3]],na.rm=TRUE)
# draw point
rgl::points3d(x=centroid_1, y=centroid_2, z=centroid_3, col=colors[whichid],cex=cex.centroid,pch=16)
}
}
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.