-
Notifications
You must be signed in to change notification settings - Fork 1
/
plotABC_3D.R
196 lines (171 loc) · 8.74 KB
/
plotABC_3D.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
#' 3D ABCplots
#'
#' Plot posterior density distribution for each generation in 3d plot window
#'
#' This opens a new interactive 3d plotting window and plots the posterior
#' density distribution of accepted particles from each generation. Several
#' options are available to add to the plot: plotting particles by weight or
#' distance, plotting particle parantage, and plotting the real parameter
#' values (if known).
#'
#' As of version 0.6.0, rejected particles are not saved for
#' outputting by the parallelized algorithm,
#' and thus they are no longer displayed by this function,
#' unlike previous versions.
#' @note
#' This function requires access to the function \code{triangulate}
#' and the \code{as} method for class \code{gpc.poly} from package \code{gpclib}.
#' As of 08-22-19, this package was not available from CRAN as a
#' Windows binary, and thus this function is likely
#' unavailable to many (if not all) Windows users.
#'
#' This function also requires the package \code{rgl}, which is usually easier to
#' obtain than package \code{gpclib}, but may not be buildable on some UNIX workstations.
#'
#' @param particleDataFrame A \code{particleDataFrame} object, as found among the output from \code{\link{doRun}} functions.
#' @param parameter Column number of parameter of interest from
#' \code{particleDataFrame}.
#' @param show.particles Option to show particles on 3d plot as "none" or as a
#' function of "weights" or "distance".
#' @param plot.parent Option to plot lines on the floor of the 3d plot to show
#' particle parantage.
#' @param realParam Option to display real parameter value as a solid line,
#' also must give actual value for this (realParamValues). Note: this should
#' only be done with simulated data where real param values are recorded.
#' @param realParamValues Value for \code{realParam}.
#' @author Barb Banbury
# @references O'Meara and Banbury, unpublished
# @keywords plotABC_3D
#' @examples
#'
#' # need to check for required suggested packages
#' if(requireNamespace("gpclib", quietly = TRUE) & requireNamespace("rgl", quietly = TRUE)){
#'
#' data(simRunExample)
#' plotABC_3D(
#' particleDataFrame = resultsBMExample[[1]]$particleDataFrame,
#' parameter = 7,
#' show.particles = "none",
#' plot.parent = FALSE,
#' realParam = FALSE,
#' realParamValues = NA
#' )
#'
#' }
#' @name plotABC_3D
#' @rdname plotABC_3D
#' @export
plotABC_3D <- function(particleDataFrame, parameter, show.particles = "none",
plot.parent = FALSE, realParam = FALSE, realParamValues = NA) {
# check if gpclib exists, if not - FAIL
has_gpclib <- requireNamespace("gpclib", quietly = TRUE)
if(!has_gpclib){stop(
"This function cannot be run without package gpclib available (Note: Windows binaries of gpclib were not unavailable as of 08-22-19).")
}
# check if rgl exists, if not - FAIL
has_rgl <- requireNamespace("rgl", quietly = TRUE)
if(!has_rgl){stop(
"This function cannot be run without package rgl available.")
}
#ugh
#library(gpclib)
#library(rgl)
#library(geometry)
# ugh ugh
#generation <- NULL #to appease R CMD CHECK
# yes??? I think this is right, not sure
generation <- particleDataFrame$generation
param.position <- parameter
pdfMaxGen <- particleDataFrame[particleDataFrame$weight>0 & particleDataFrame$generation == max(particleDataFrame$generation), ]
nParticles <- dim(pdfMaxGen)[1]
nparams <- dim(particleDataFrame)[2]-6
#
q <- vector() #vector of x vals
r <- vector() #vector of y vals
s <- vector() #generation each x-y coord is found
#
rangeParam <- (max(particleDataFrame[particleDataFrame$weight>0, param.position])
-min(particleDataFrame[particleDataFrame$weight>0, param.position]))
#
if (rangeParam != 0) {
v <- vector("list", max(particleDataFrame$generation))
for (i in 1:max(particleDataFrame$generation)){
which.gen <- (i+1)-1 # the hell is this?
pdfParam <- particleDataFrame[(particleDataFrame$weight>0 & particleDataFrame$generation == i), param.position]
pdfSix <- particleDataFrame[particleDataFrame$weight>0 & particleDataFrame$generation == i, 6]
v[[i]] <- density(pdfParam, weights = nParticles*pdfSix/sum(nParticles*pdfSix))
q <- c(q, v[[i]]$x)
r <- c(r, v[[i]]$y)
#s <- c(s, v[[i]]$x) # return a$generation which v[[i]]
#
s <- c(s, rep(which.gen, times = length(v[[i]]$x)))
T <- as.matrix(cbind(q, r, s))
}
x <- T[, 1]
y <- T[, 2]
z <- T[, 3]
#
# plotting with rgl
rgl::open3d() #make bigger window
#bg3d("color) #gives background color for plot
rgl::plot3d(x, y, z, col = "black", box = FALSE, type = "n", xlab = "", ylab = "", zlab = "",
zlim = c(0, max(particleDataFrame$generation)), ylim = c(0, max(y)))
#message(paste("HERE"))
rgl::rgl.viewpoint(35, 1, 90) #sets viewpoint for initial plot
rgl::title3d(colnames(x)[param.position], col = 'red', pos = c(NA, -2, max(z)))
#text3d(x = min(x), y = mean(y), z = max(z), text = "Density" col = 'blue')
#title3d("Starting States", col = 'purple', pos = c(NA, 0, max(which.gen)))
#
for (i in 1:max(s)){
ngen <- (i+1)-1
triangles <- gpclib::triangulate(as(cbind(x[which(z == i)], y[which(z == i)]), "gpc.poly"))
xyzGen <- data.frame(zG = z[which(z == i)], xG = x[which(z == i)], yG = y[which(z == i)])
zfit <- predict(lm(zG ~ xG + yG, xyzGen), newdata = data.frame(xG = triangles[, 1], yG = triangles[, 2]))
opacity <- 0.8*(ngen/length(v))
rgl::rgl.material(color = "black", alpha = opacity, lit = FALSE)
#message(dim(triangles))
#message(dim(zfit))
if(length(zfit) != (dim(triangles)[1])){
stop("Error, predict() not properly returning vector for same number of input")
}else{
rgl::triangles3d(cbind(triangles, zfit), col = "red")
}
#readline(prompt = "hit enter ")
if (realParam) {
rgl::rgl.material(color = "blue", lwd = 2)
rgl::lines3d(x = c(realParamValues[1], realParamValues[1]), y = c(0, 0), z = c(min(s), max(s)))
}
}
}
else {
warning(paste("You are attempting to plot", colnames(x)[param.position], ", which is a fixed value")) #return fixed val
}
show.particles <- match.arg(arg = show.particles, choices = c("none", "weights", "distance"), several.ok = FALSE)
#kept <- subset(particleDataFrame[which(particleDataFrame$id>0), ])[, ]
#reject <- subset(particleDataFrame[which(particleDataFrame$id<0), ])[, ]
#short.kept <- subset(kept[which(kept$generation>1), ])[, ]
#
kept <- particleDataFrame[particleDataFrame$id>0, ]
short.kept <- kept[kept$generation>1, ]
if (show.particles == "none"){
message("Note: Currently not plotting particles. To plot particles modify the show.particles = argument.")
}
else if (show.particles == "weights") {
for (j in 1:(dim(kept)[1])) {
circle.size <- (kept[j, 6]/max(kept[, 6]))*(.1*(max(kept[, param.position])-min(kept[, param.position])))*10
rgl::spheres3d(x = kept[j, param.position], y = -1*(.05*(max(y)-min(y))), z = kept[j, 1], radius = circle.size, type = "s", col = "black", add = TRUE, aspect = TRUE)
}
}
else if (show.particles == "distance") {
for (j in 1:(dim(kept)[1])) {
circle.size <- (kept[j, 5]/max(kept[, 5]))*(.1*(max(kept[, param.position])-min(kept[, param.position])))*15
rgl::spheres3d(x = kept[j, param.position], y = -1*(.05*(max(y)-min(y))), z = kept[j, 1], radius = circle.size, type = "s", col = "black", add = TRUE, aspect = TRUE)
}
}
if (plot.parent){
for (k in 1:(dim(short.kept)[1])) {
prev.gen <- subset(kept[which(kept$generation == short.kept[k, 1]-1), ])[, ] #works to retreive prev gen
rgl::rgl.material(color = "black", alpha = 0.5)
rgl::lines3d(x = c(short.kept[k, param.position], prev.gen[short.kept[k, ]$parentid, param.position]), y = c(0, 0), z = c(short.kept[k, 1], short.kept[k, 1]-1)) }
}
}