Nothing
##-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-##
## Compute and visualize the high dimensional function
##-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-##
plot.msc.svm <- function(x, drawStdDev=FALSE, span=0.5, nsamples=50, plot=TRUE, colorMap=0,...) {
plot.msc(x, drawStdDev, span, nsamples)
}
plot.msc.kd <- function(x, drawStdDev=FALSE, span=0.5, nsamples=50, plot=TRUE, colorMap=0,...) {
plot.msc(x, drawStdDev, span, nsamples)
}
plot.msc <- function(x, drawStdDev=FALSE, span=0.5, nsamples=50, plot=TRUE, colorMap=0,...) {
## Compute the geometry
geom <- geometry.mscPlot(x, span=span, nsamples=nsamples)
## Compute the vis scene
scene <- scene.mscPlot(geom, colorMap=colorMap)
## Create the mscPlot object
obj <- structure(list(geom=geom, scene=scene),class="mscPlot")
## Create the mscPlot
if(plot)
obj <- plot.mscPlot(obj, drawStdDev)
## Return the object
invisible(obj)
}
##-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-##
## Plot the rgl scene
##-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-##
plot.mscPlot <- function(x, drawStdDev=FALSE, axesOn=TRUE, ...) {
## Initialize the 3D scene
if(length(x$dev) == 0 || rgl.cur() == 0) {
x$dev <- open3d(mouseMode=c("none","trackball", "zoom"), windowRect=c(0, 0, 512, 512))
view3d(phi=0, theta=0, fov=15)
light3d(phi=-45, theta=-45)
light3d(phi=45, theta=45)
bg3d(sphere=FALSE, color="white")
}
else
rgl.clear()
## Draw the spheres
spheres3d(x$geom$pMins3d, radius=0.01, col=rgb(x$scene$colormap(x$geom$pMins3d[,3]), maxColorValue=256))
spheres3d(x$geom$pMaxs3d, radius=0.01, col=rgb(x$scene$colormap(x$geom$pMaxs3d[,3]), maxColorValue=256))
## Draw the cylinders
cIds = c()
for(i in 1:length(x$scene$cylinders))
cIds[i] <- shade3d(addNormals(x$scene$cylinders[[i]]), col=x$scene$cylColors1[[i]])
x$cIds <- cIds
## Draw the optional standard deviation tubes
if(drawStdDev) {
for(i in 1:length(x$scene$tubes)) {
shade3d(addNormals(x$scene$tubes[[i]]), col="white", alpha=0.25, back="fill", front="cull")
shade3d(addNormals(x$scene$tubes[[i]]), col="white", alpha=0.25, back="cull", front="fill")
}
}
## Create the orientation axes
if(axesOn) {
text3d(1, 0, 0, texts="x1")
text3d(0, 1, 0, texts="x2")
text3d(0, 0, 1, texts="f(x)")
lines3d(c(1, 0, 0),c(0, 0, 1), c(0,0, 0), color="black")
segments3d(c(0, 0),c(0, 0), c(0, 1), color="red")
}
## Create the mouse selection mode
x=mouseSelect(1, x)
## Return the object
invisible(x)
}
## @-@-@-@-@-@-@ Plot Curves @-@-@-@-@-@-@-@-@-@ ##
plotCurves <- function(obj) {
## Calculate the number of pages, columns and rows
if(length(obj$plotList) == 0)
obj$plotList <- c( 1:ncol(obj$geom$curvesXHD[[1]]))
nc <- length(obj$plotList)
numPages=nc%/%25
if(nc%%25 != 0)
numPages = numPages+1
col=(nc/numPages)%/%5+1
if(col > 1)
row=5
else
row=nc%%5
## Create the number of pages
if(length(dev.list()) == 0) {
for(d in 1:numPages)
dev.new()
}
for(c in 1:nc) {
if(c%%25 == 1) {
dev.set(which=c%/%25+2)
par(bg = "white", mfcol=c(row, col), mar=c(5, 4, 1, 1))
}
cindex <- obj$plotList[c]
title <- paste("Plot #", cindex, ": ", obj$geom$colNames[cindex], " vs ", "Y")
yLabs <- c(sprintf("%.2f", obj$geom$ylim[1]), sprintf("%.2f", (obj$geom$ylim[2]-obj$geom$ylim[1])/2),
sprintf("%.2f", obj$geom$ylim[2]))
plot(y=NA, x=NA, xlim=obj$geom$xlim[cindex,], ylim=range(0:1),
main=title, xlab=obj$geom$colNames[cindex], ylab="Y", yaxt="n", xaxt="n")
axis(1, labels=c(0.0,0.5,1.0), at=c(0.0,0.5,1.0))
axis(2, labels=yLabs, at=c(0.0,0.5,1.0))
for(i in 1:length(obj$scene$cylinders)){
if(obj$scene$state[[i]] == 2){
lightColor=hex(mixcolor(0.25, hex2RGB(obj$scene$cylColors2[[i]]), RGB(1,1,1)))
lines(x=obj$geom$curvesXHD[[i]][,cindex], y=obj$geom$crvs3d[[i]][,3],
col=obj$scene$cylColors2[[i]], lwd=3)
lines(x=obj$geom$curvesXHD[[i]][,cindex] +
obj$geom$residuals[[i]][,cindex],
y=obj$geom$crvs3d[[i]][,3], col=lightColor, lwd=1)
lines(x=obj$geom$curvesXHD[[i]][,cindex] -
obj$geom$residuals[[i]][,cindex],
y=obj$geom$crvs3d[[i]][,3], col=lightColor, lwd=1)
}
}
}
}
## - - - - - - - MOUSE SELECT - - - - - ##
mouseSelect <- function(button, obj) {
## (-)-(-)-(-)-(-) MOUSE BEGIN (-)-(-)-(-)-(-) ##
mouseBegin <- function(x, y) {
##-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-##
## Define the area of the window that is selected,
## return a function that tests if the cylinders is selected
createSelectFunction <- function(xa, ya, w=3, h=3) {
proj <- rgl.projection()
vp <- proj$view
llx <- (xa - vp[1]) / vp[3]
lly <- 1-(ya - vp[2]) / vp[4]
urx <- (xa -vp[1] + w) / vp[3]
ury <- 1-(ya - vp[2] +h) / vp[4]
if ( llx > urx ){
temp <- llx
llx <- urx
urx <- temp
}
if ( lly > ury ){
temp <- lly
lly <- ury
ury <- temp
}
function(x,y=NULL,z=NULL) {
pixel <- rgl.user2window(x,y,z,projection=proj)
x <- pixel[,1]
y <- pixel[,2]
z <- pixel[,3]
(llx <= x) & (x <= urx) & (lly <= y) & (y <= ury) &
(0 <= z) & (z <= 1)
}
}
## Create the picker function
picker <- createSelectFunction(x, y)
selected <- c(0,1)
## Cycle through all curves
for(i in 1:length(obj$scene$cylinders)){
## Test to see if this curve was selected
selector <- picker(t(obj$scene$cylinders[[i]]$vb[1:3,]))
## Check each point. If picked, save smallest Z
for(j in 1:length(selector)){
if(selector[j] == TRUE){
thisZ <- obj$scene$cylinders[[i]]$vb[3,j]
if(thisZ < selected[2])
selected <- c(i, thisZ)
}
}
}
## If we have selected more than is supported, throw an error
if(selected[1] != 0){
noSelect <- FALSE
if(length(obj$selColors)==1){
if(obj$scene$state[[selected[1]]] == 1)
{
noSelect <- TRUE
print("ERROR: Too Many Selected Curves!! Unselect a curve.")
}
}
## If in selected state, pop the selected object
if(!noSelect){
rgl.pop(id=obj$cIds[[selected[1]]])
}
## Re-draw the tubes if selected
addColor <- FALSE
takeColor <- FALSE
if(obj$scene$state[[selected[1]]] == 1 && !noSelect){
obj$scene$state[[selected[1]]] <<- 2
obj$scene$cylColors2[selected[1]] <<- obj$scene$selectColors[1]
obj$cIds[[selected[1]]] <<- shade3d(addNormals(obj$scene$cylinders[[selected[1]]]),
col=obj$scene$cylColors2[[selected[1]]])
takeColor <- TRUE
}
else if(obj$scene$state[[selected[1]]] == 2)
{
obj$scene$state[[selected[1]]] <<- 1
obj$cIds[[selected[1]]] <<- shade3d(addNormals(obj$scene$cylinders[[selected[1]]]),
col=obj$scene$cylColors1[[selected[1]]])
addColor <- TRUE
}
## Add or take colors from the color vector
if(addColor)
obj$scene$selectColors <<- c(obj$scene$cylColors2[[selected[1]]], obj$scene$selectColors)
else if(takeColor)
obj$scene$selectColors <<- obj$scene$selectColors[2:length(obj$scene$selectColors)]
## Plot the curves
plotCurves(obj)
}
invisible(obj)
}
## Set the rgl mouse call backs
rgl.setMouseCallbacks(button, mouseBegin)
invisible(obj)
}
##-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-##
## Generate the vis pieces
##-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-##
scene.mscPlot <- function(geom, cSides=8, colorMap=0) {
## Create a list of selected colors
selectColors <- c()
brewerSets=c("Set1", "Dark2", "Set3")
for(b in 1:length(brewerSets)){
cols <- brewer.pal(brewer.pal.info[brewerSets[b],]$maxcolors, brewerSets[b])
selectColors <- c(selectColors, cols)
}
selectColors <- c(selectColors, "#000000")
## The blue, green, red colormap
if(colorMap == 0)
colormap = colorRamp( c('#0066CC', '#CCCC00', '#D22905') , interpolate="linear", bias=0.5)
else if(colorMap == 1)
colormap = colorRamp(c('#0520B0', '#f7f7f7', '#CA0020'), interpolate="linear", bias=0.5)
else if(colorMap == 2)
colormap = colorRamp(c('#7b3294', '#f7f7f7', '#008837'), interpolate="linear", bias=0.5)
## Create the cylinders
cylinders <- c()
cylColors1 <- c()
cylColors2 <- c()
state <- c()
tubes <- c()
for(i in 1:geom$npartition){
color <- rgb(colormap(geom$crvs3d[[i]][,3]), maxColorValue=255)
fullColors <- c()
for(p in 1:length(geom$crvs3d[[i]][,3]))
fullColors[((p-1)*(cSides*4)+1) : (p*(cSides*4))] <- color[p]
## Create the tubes
cylColors1[[i]] <- fullColors
cylColors2[[i]] <- "black"
cylinders[[i]] <- cylinder3d(geom$crvs3d[[i]], twist = 0, sides = cSides,
e2=rbind(c(1,0,0),c(1,0,0)), radius = 0.008)
tubes[[i]] <- cylinder3d(geom$crvs3d[[i]], twist=0, sides=cSides,
e2=rbind(c(1,0,0),c(1,0,0)), radius = geom$totResiduals[[i]])
state[[i]] <- 1
}
scene <- structure(list(geom=geom,cylinders=cylinders, cylColors1=cylColors1, cylColors2=cylColors2,
state=state, tubes=tubes, colormap=colormap, selectColors=selectColors))
invisible(scene)
}
##-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-##
## Generate the geometry for the scene
##-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-##
geometry.mscPlot <- function(ms, span = 0.35, nsamples = 50)
{
msLevel <- ms$level[[ms$predictLevel]]
## Unique min and max indicies into ms$x/ms$y (domain/function value)
uMins <- unique(msLevel$mins)
uMaxs <- unique(msLevel$maxs)
nmins <- length(uMins)
nmaxs <- length(uMaxs)
## Do PCA on extrema
extrema <- ms$x[c(uMins, uMaxs), ]
pcaE <- prcomp(extrema)
## Location of mins and maxs
pMins <- pcaE$x[1:nmins, 1:2]
pMaxs <- pcaE$x[(nmins+1):(nmins+nmaxs), 1:2]
# Scale the function values to be within the bounding box (0-1)
miny <- min(ms$y)
maxy <- max(ms$y)
# ms$y <- (ms$y - miny) / (maxy-miny)
## Function values of mins and maxs
fMins <- ms$y[uMins]
fMaxs <- ms$y[uMaxs]
## The number of partition
npartition <- length(msLevel$partitionSize)
## For each crystal store curve points (curvesX) and associated function value (curvesY)
## Curves are of the from x_i = c_i(y), e.g. a n-dimensional curve
## curve projected into 2D PCA sapce
crvs3d <- c()
## curves in original space
curvesXHD <- c()
## residuals for each curve, for each dimension
residuals <- c()
for(i in 1:npartition){
crvs3d[[i]] <- matrix(nrow=nsamples, ncol=3)
residuals[[i]] <- matrix(nrow=nsamples, ncol=ncol(ms$x))
eMax <- msLevel$maxs[i]
eMin <- msLevel$mins[i]
miny <- ms$y[eMin]
maxy <- ms$y[eMax]
h <- (maxy-miny) / (nsamples-1)
sharedMax <- unique( which(msLevel$maxs == eMax) )
sharedMin <- unique( which(msLevel$mins == eMin) )
y <- ms$y[msLevel$partition == i]
x <- ms$x[msLevel$partition == i, ]
for(k in sharedMin){
if(k != i){
y2 <- ms$y[msLevel$partition == k]
y2 <- miny - (y2 - miny)
y <- c(y, y2)
x <- rbind(x, ms$x[msLevel$partition == k, ])
}
}
for(k in sharedMax){
if(k != i){
y2 <- ms$y[msLevel$partition == k]
y2 <- maxy + (maxy-y2)
y <- c(y, y2)
x <- rbind(x, ms$x[msLevel$partition == k, ] )
}
}
ye <- miny + (0:(nsamples-1))*h
ye[nsamples] <- maxy
xp <- matrix(nrow=nsamples, ncol=ncol(ms$x))
for(k in 1:ncol(xp)){
l <- loess(x[,k] ~ y, span=span, normalize=FALSE, degree=1)
rl <- loess(I(abs(l$residuals)) ~ y, span=span/2.5, normalize=FALSE, degree=1)
residuals[[i]][,k] <- predict(rl, ye)
xp[, k] <- predict(l, ye)
}
curvesXHD[[i]] <- xp
pcaC <- prcomp(xp)
iMin <- which(uMins == eMin)
iMax <- which(uMaxs == eMax)
xp <- pcaC$x[, 1:2]
if(nmins == 1)
anchor <- pMins
else
anchor <- pMins[iMin, ]
xp[, 1] <- xp[, 1] - xp[1, 1] + anchor[1]
xp[, 2] <- xp[, 2] - xp[1, 2] + anchor[2]
if(nmaxs == 1){
stretch <- pMaxs - xp[nsamples, ]
}
else
stretch <- pMaxs[iMax, ] - xp[nsamples, ]
for(k in 1:nsamples)
xp[k, ] <- xp[k, ] + (k-1)/(nsamples-1) * stretch
crvs3d[[i]][, 1:2] <- xp
crvs3d[[i]][,3] <- ye
}
## Extend the points & curves to 3D
if(nmins ==1){
pMins3d <- t(as.matrix(c(pMins, fMins)))
}
else{
pMins3d <- cbind(pMins, fMins)
}
if(nmaxs == 1){
pMaxs3d <- t(as.matrix(c(pMaxs, fMaxs)))
}
else{
pMaxs3d <- cbind(pMaxs, fMaxs)
}
pTmins <- c(min(c(pMins3d[,1], pMaxs3d[,1])), min(c(pMins3d[,2], pMaxs3d[,1])), min(ms$y))
pTmaxs <- c(max(c(pMins3d[,1], pMaxs3d[,1])), max(c(pMins3d[,2], pMaxs3d[,2])), max(ms$y))
for(i in 1:npartition){
crvs3d[[i]][,1] <- (crvs3d[[i]][,1] - pTmins[1])/(pTmaxs[1] - pTmins[1])
crvs3d[[i]][,2] <- (crvs3d[[i]][,2] - pTmins[2])/(pTmaxs[2] - pTmins[2])
crvs3d[[i]][,3] <- (crvs3d[[i]][,3] - pTmins[3])/(pTmaxs[3] - pTmins[3])
}
for(i in 1:nmins)
pMins3d[i,] <- (pMins3d[i,] - pTmins)/ (pTmaxs - pTmins)
for(i in 1:nmaxs)
pMaxs3d[i,] <- (pMaxs3d[i,] - pTmins)/ (pTmaxs - pTmins)
## Get the complete residual and scale appropriatly
totResiduals <- c()
totScale <- 1 / mean(dist(extrema))
for(i in 1:npartition)
totResiduals[[i]] <- sqrt(rowSums(residuals[[i]]^2))*totScale
# Set the limits of the x
xlim <- matrix(nrow=ncol(ms$x), ncol=2)
for( i in 1:ncol(ms$x))
xlim[i,] <- range(ms$x[, i])
## Return lots of info
mscGeom <- structure(list(npartition=npartition,
curvesXHD = curvesXHD,
crvs3d=crvs3d,
residuals = residuals,
totResiduals=totResiduals,
pMins3d=pMins3d,
pMaxs3d=pMaxs3d,
ylim= c(pTmins[3], pTmaxs[3]),
xlim = xlim,
colNames = colnames(ms$x)))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.