library(shiny)
library(mwaved)
# Check if ggplot2 is available and use it
ggplot2Avail <- 'ggplot2' %in% .packages(all.available = TRUE)
# Graphical settings
shrinkLabel <- 'Shrinkage Type:'
shrinkChoices <- list("Hard" = "Hard", "Soft" = "Soft", "Garrote" = "Garrote")
fill.col <- 'light blue'
outline.col <- 'blue'
density.col <- 'dark grey'
highlight.col <- 'red'
if (ggplot2Avail){
library(ggplot2)
blanklabs <- labs(x = '', y = '')
}
primenum <- c(83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199)
lsize=0.5;
hsize=1;
xcat <- function(obj){
paste('c(',paste(as.numeric(obj), collapse = ','),')',sep='')
}
numclean <- function(obj){
return(round(as.numeric(obj), 3))
}
shinyServer(function(input, output, session) {
###############################################
## Compute Signals and other bits
###############################################
currentShrinkage <- 'Hard'
sigList <- reactive({
n <- 2^input$J
m <- input$m;
if (m == 1){
SNR <- sample(input$SNR[1]:input$SNR[2], 1)
} else {
SNR <- sample(seq(from = input$SNR[1], to = input$SNR[2], length = m), m);
}
alpha <- sample(seq(from = input$alpha[1], to = input$alpha[2], length = m), m);
shape <- sample(seq(from = input$gshape[1], to = input$gshape[2], length = m), m);
scale <- sample(seq(from = input$gscale[1], to = input$gscale[2], length = m), m);
width <- 1/sqrt(sample(primenum, m))
G <- switch(input$blur,
'smooth'= gammaBlur(n, shape, scale),
'direct' = directBlur(n, m),
'box.car'= boxcarBlur(n, width))
signal <- switch(input$sig,
'lidar' = makeLIDAR(n),
'doppler' = makeDoppler(n),
'bumps' = makeBumps(n),
'blocks' = makeBlocks(n),
'cusp' = makeCusp(n),
'heavisine' = makeHeaviSine(n))
signalName <- switch(input$sig,'lidar' = "LIDAR Signal",
'doppler' = "Doppler Signal",
'bumps' = "Bumps Signal",
'blocks' = "Blocks Signal",
'cusp' = "Cusp Signal",
'heavisine' = "HeaviSine Signal")
signalBlur <- blurSignal(signal, G)
sigma <- sigmaSNR(signalBlur, SNR)
eps <- multiNoise(n, sigma, alpha)
Y <- signalBlur + eps;
x <- seq(from = 0, to = 1 - 1/n, length = n);
resolution <- input$resolution
list(m = m, n = n, signal = Y, G = G, blur = input$blur, resolution = input$resolution, SNR = SNR,
alpha = alpha, shape = shape, scale = scale , G = G, trueSignal = signal, eps = eps,
signalBlur = signalBlur, width = width, x = x, signalName = signalName, sigma = sigma)
})
mWaveDList <- reactive({
multiSig <- sigList()
if (input$etaChoose == "TRUE"){
eta <- input$eta
mlwvd <- multiWaveD(multiSig$signal, multiSig$G, alpha = multiSig$alpha, eta = eta,
resolution = multiSig$resolution, shrinkType = input$shrinkage1,
deg = as.integer(input$degree))
} else {
mlwvd <- multiWaveD(multiSig$signal, multiSig$G, alpha = multiSig$alpha,
resolution = multiSig$resolution, shrinkType = input$shrinkage1,
deg = as.integer(input$degree))
}
j0 <- mlwvd$j0
j1 <- mlwvd$j1
list(n = multiSig$n, m = multiSig$m, mWaveD = mlwvd, j0 = mlwvd$j0, j1 = mlwvd$j1, coef = mlwvd$coef, shrinkCoef = mlwvd$shrinkCoef)
})
###############################################
## Render WvdPlots
###############################################
output$reactiveSignal <- renderPlot({
sList <- sigList()
# Base signal selected
if (input$signalShow == 1){
plotTitle <- sList$signalName
if (ggplot2Avail){
plot.df <- data.frame(x = sList$x, signal = sList$trueSignal)
print(ggplot(data = plot.df, aes(x = x, y = signal)) + geom_line(size = lsize) + ggtitle(plotTitle) + blanklabs)
} else {
plot(sList$x, sList$trueSignal, type = 'l', main = plotTitle, ylab = '', xlab = '')
grid()
}
# + reduce.margins
}
if (input$signalShow == 2){
plotTitle <- paste('Blurred',sList$signalName)
if (ggplot2Avail){
plot.df <- data.frame(Y = as.vector(sList$signalBlur), x = rep(sList$x, sList$m),
Channel = rep(LETTERS[1:sList$m], each = sList$n))
print(ggplot(plot.df, aes(x = x, y = Y, colour=Channel)) + geom_line(size = lsize, alpha = 0.9) + ggtitle(paste('Blurred ', sList$signalName, sep = "")) + blanklabs )
} else {
par(oma = c(4, 1, 1, 1))
plot.out <- matplot(sList$x, sList$signalBlur, type = 'l', main = plotTitle, xlab = '', ylab = '')
grid()
par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
legend("bottom", paste('Channel ', LETTERS[1:sList$m]), xpd = TRUE, horiz = TRUE,
inset = c(0, 0), bty = "n", lty = 1:sList$m, col = 1:sList$m, cex = 1)
}
# + reduce.margins
# }
}
if (input$signalShow == 3){
plotTitle <- paste('Noisy blurred', sList$signalName)
if (ggplot2Avail){
plot.df <- data.frame(Y = as.vector(sList$signal), x = rep(sList$x, sList$m), Channel = rep(LETTERS[1:sList$m], each = sList$n))
print(ggplot(plot.df, aes(x = x, y = Y, colour = Channel)) + geom_line(size = lsize, alpha = 0.5) + ggtitle(plotTitle) + blanklabs)
} else {
par(oma = c(4, 1, 1, 1))
matplot(sList$x, sList$signal, type = 'l', main = plotTitle, xlab = '', ylab = '')
grid()
par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
legend("bottom", paste('Channel ', LETTERS[1:sList$m]), xpd = TRUE, horiz = TRUE,
inset = -0.1, bty = "n", lty = 1:sList$m, col = 1:sList$m, cex = 1)
}
}
}, height=function() { session$clientData$output_reactiveSignal_width * 3/4 })
output$multiPlot <- renderPlot({
mList <- mWaveDList()
mWaveD <- mList$mWaveD
plot(mWaveD, which = 4)
}, height=function() { session$clientData$output_multiPlot_width * 3/4 })
output$wvdPlot <- renderPlot({
sList <- sigList()
mList <- mWaveDList()
mlwvd <- mList$mWaveD
n <- mList$n
m <- mList$m
# Number of things to plot
k = 2
plotTitle <- paste(input$shrinkage1, ' thresholded estimate of ', sList$signalName,sep="")
if (ggplot2Avail){
out <- c(sList$trueSignal,mlwvd$estimate)
lty <- c(rep('solid',n),rep('dashed',n))
group <- rep(c("True Signal", "Estimate"), each = n)
if (input$wvdshow != 3){
k <- 3
extra <- rep('dashed', n)
if (m == 1){
out = c(out,multiEstimate(mlwvd$signal , G = mlwvd$G , resolution = mlwvd$resolutionMethod, alpha = mlwvd$alpha, shrinkType = mlwvd$shrinkType, deg = as.integer(input$degree)))
if (input$wvdshow == 2){
plot.lab <- 'Naive mWaveD average'
} else {
plot.lab <- 'Best channel'
}
group <- c(group, rep(plot.lab, n))
lty <- c(lty,extra)
} else {
if (input$wvdshow == 2){
wvdMean = apply(sapply(1:m, function(x)
multiEstimate(as.matrix(mlwvd$signal[, x]) ,as.matrix(mlwvd$G[, x]) , resolution = mlwvd$resolutionMethod, alpha = mlwvd$alpha[x], shrinkType = mlwvd$shrinkType, deg = as.integer(input$degree))), 1, mean)
out <- c(out,wvdMean)
group <- c(group,rep('Naive mWaveD average',n))
lty <- c(lty,extra)
} else {
resolution.ty <- mlwvd$resolution
if (resolution.ty == 'smooth'){
best <- which.max(mlwvd$blurInfo$freqCutoffs)
} else {
best <- which.min(mlwvd$sigmaEst)
}
wvdBest <- multiEstimate(as.matrix(mlwvd$signal[, best]), as.matrix(mlwvd$G[, best]), resolution = mlwvd$resolutionMethod, deg = as.integer(input$degree))
out <- c(out, wvdBest)
group <- c(group, rep("mWaveD best chan", n))
lty <- c(lty, extra)
}
}
}
est.df = data.frame(Y = out, x = rep(1:n/n, k), group = group, lty = lty)
est.plot <- ggplot(est.df, aes(x = x, y = Y, colour = group)) + geom_line(aes(linetype = lty, colour = group), size = lsize + 0.25) + scale_linetype(guide = 'none') + ggtitle(plotTitle)
print(est.plot)
} else {
if (input$wvdshow != 3){
k <- 3
if (input$wvdshow == 2){
xtra = apply(sapply(1:m, function(x)
multiEstimate(as.matrix(mlwvd$signal[, x]) ,as.matrix(mlwvd$G[, x]) , blur = 'smooth', alpha = mlwvd$alpha[x], shrinkType = mlwvd$shrinkType, deg = as.integer(input$degree))), 1, mean)
xtralab <- 'Naive mean'
} else {
resolution.ty <- mlwvd$resolution
if (resolution == 'smooth'){
best = mlwvd$blurInfo$bestChannel
} else {
best = which.min(mlwvd$sigmaEst)
}
xtra <- multiEstimate(as.matrix(mlwvd$signal[, best]),as.matrix(mlwvd$G[, best]), blur = 'smooth', deg = as.integer(input$degree))
xtralab <- 'Best channel'
}
} else {
xtra <- NULL
xtralab <- NULL
}
Y <- cbind(sList$trueSignal,mlwvd$estimate, xtra)
t = (0:(n-1))/n
labels <- c("True Signal",'mWaveD estimate',xtralab)
par(oma = c(4, 1, 1, 1))
matplot(t, Y, type = 'l', xlab = 'x', ylab = '', main = plotTitle)
grid()
par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
legend("bottom", labels, xpd = TRUE, horiz = TRUE,
inset = -0.1, bty = "n", lty = 1:3, col = 1:3, cex = 1)
}
}, height=function() { session$clientData$output_wvdPlot_width * 3/4 })
output$resolutionPlot <- renderPlot({
mList = mWaveDList()
mWaveD <- mList$mWaveD
n <- mList$n
m <- mList$m
iw = -(n/2 - 1):(n/2)
blurInfo = mWaveD$blurInfo
resType = mWaveD$resolutionMethod
blurDetected = mWaveD$blurDetected
plot(mWaveD, which = 3)
}, height=function() { session$clientData$output_resolutionPlot_width * 3/4 })
output$summaryOut <- renderPrint({
mList <- mWaveDList()
print(summary(mList$mWaveD))
})
output$summaryWVD <- renderPrint({
wList <- mWaveDList()
sList <- sigList()
cat("Degree of Meyer poly =", wList$mWaveD$degree, "\n")
cat("# obs per channel =", wList$n, "\n\n")
cat("Max possible resolution =", log2(wList$n) - 1, "\n")
cat("Coarse resolution, j0 =", wList$j0,"\n")
cat("mWaveD optimal res, j1 =", wList$j1)
cat("\n\n")
cat("Thresholding method:", wList$mWaveD$shrinkType, "\n")
cat("Tuning parameter: eta =",wList$mWaveD$eta,'\n\n')
threshMatrix <- cbind(round(wList$mWaveD$levelMax,4), round(wList$mWaveD$thresh,4), round(wList$mWaveD$percent,2))
rownames(threshMatrix) <- paste("Res", wList$j0:wList$j1,":" )
colnames(threshMatrix) <- c("max|coef|", "thresh", "%shrunk" )
print(threshMatrix)
cat('\n Measure of fit:\n')
cat('L_1 norm =', sum(abs(sList$trueSignal - wList$mWaveD$estimate))/sList$n, '\n')
cat('L_2 norm =', sqrt(sum((sList$trueSignal - wList$mWaveD$estimate)^2)/sList$n), '\n')
})
output$summaryResolution <- renderPrint({
wList <- mWaveDList()
cat("# obs per channel =", wList$n, "\n\n")
cat("# channels: m = ", wList$m, ".\n", sep ='')
cat("Max possible resolution =", log2(wList$n) - 1, "\n")
cat("Coarse resolution, j0 =", wList$j0,"\n")
cat("mWaveD optimal res, j1 =", wList$j1)
cat("\n\n")
cat("Estimated Channel information:\n\n")
if (wList$mWaveD$resolutionMethod == "smooth"){
if (wList$mWaveD$blurDetected == "direct") {
mat <- cbind(round(wList$mWaveD$sigmaEst, 3), round(wList$mWaveD$alpha, 3), wList$mWaveD$blurInfo$freq, rep(wList$mWaveD$j1, wList$m))
colnames(mat) <- c("sigma", "alpha", "cutoff", "best res")
rownames(mat) <- paste("Chan ", 1:wList$m,':', sep = '')
print(mat)
} else {
mat <- cbind(round(wList$mWaveD$sigmaEst, 3), round(wList$mWaveD$alpha, 3), wList$mWaveD$blurInfo$freq, wList$mWaveD$blurInfo$maxLevels)
colnames(mat) <- c("sigma", "alpha", "cutoff", "best res")
rownames(mat) <- paste("chan: ", 1:wList$m)
print(mat)
cat("\n")
cat("Estimated best channel = ", wList$mWaveD$blurInfo$bestChannel)
}
} else {
mat <- cbind(round(wList$mWaveD$sigma, 3), round(wList$mWaveD$alpha, 3))
colnames(mat) <- c("sigma", "alpha")
rownames(mat) <- paste("chan: ", 1:wList$m)
print(mat)
}
})
output$summaryMRA <- renderPrint({
wList <- mWaveDList()
cat("Degree of Meyer poly =", wList$mWaveD$degree, "\n")
cat("# obs per channel =", wList$n, "\n\n")
cat("Max possible resolution =", log2(wList$n) - 1, "\n")
cat("Coarse resolution, j0 =", wList$j0,"\n")
cat("mWaveD optimal res, j1 =", wList$j1)
cat("\n\n")
cat("Thresholding method:", wList$mWaveD$shrinkType, "\n")
cat("Tuning parameter: eta =",wList$mWaveD$eta,'\n\n')
threshMatrix <- cbind(round(wList$mWaveD$levelMax,4), round(wList$mWaveD$thresh,4), round(wList$mWaveD$percent,2))
rownames(threshMatrix) <- paste("Res", wList$j0:wList$j1,":" )
colnames(threshMatrix) <- c("max|coef|", "thresh", "%shrunk" )
print(threshMatrix)
})
output$summarySignal <- renderPrint({
sList <- sigList()
cat("# obs per channel = ", sList$n, ".\n", sep = '')
cat("# channels: m = ", sList$m, ".\n", sep ='')
cat('Blur type: ',sList$blur, '\n\n')
cat("Channel information:\n\n")
mat <- cbind(round(sList$sigma, 3), round(sList$alpha, 3), round(sList$SNR,3))
colnames(mat) <- c("sigma", "alpha", "SNR")
rownames(mat) <- paste("Chan ", 1:sList$m," :", sep ='')
if ( sList$blur == 'smooth'){
mat <- cbind(mat, round(sList$shape, 3), round(sList$scale, 3))
colnames(mat) <- c("sigma", "alpha","SNR", "shape", "scale")
} else {
if (sList$blur == 'box.car'){
mat <- cbind(mat, paste(round(sList$width, 3), " = 1/sqrt(", 1/sList$width^2 ,")", sep = '') )
colnames(mat) <- c("sigma", "alpha","SNR", "Box.car widths")
}
}
print(mat, quote = FALSE)
})
output$signalCalls <- renderPrint({
cat('Base R Function calls:\n\n')
sList <- sigList()
cat('J <-', log2(sList$n), '\n')
cat('n <- 2^J', '\n')
cat('m <-', sList$m, '\n')
switch(input$sig,
'lidar' = cat('signal <- makeLIDAR(n)', '\n'),
'doppler' = cat('signal <- makeDoppler(n)', '\n'),
'bumps' = cat('signal <- makeBumps(n)', '\n'),
'blocks' = cat('signal <- makeBlocks(n)', '\n'),
'cusp' = cat('signal <- makeCusp(n)\n'),
'heavisine' = cat('signal <- makeHeaviSine(n)\n'))
cat('t <- (0:(n-1))/n\n')
if (input$signalShow == '1'){
cat("plot(t, signal, type = 'l')\n")
}
if (input$signalShow == '2' || input$signalShow == '3'){
if (sList$blur == 'direct'){
cat('G <- directBlur(n, m)\n')
}
if (sList$blur == 'smooth'){
if (sList$m == 1){
cat('shape <-', numclean(sList$shape),'\n')
cat('scale <-', numclean(sList$scale),'\n')
} else {
cat('shape <-',xcat(numclean(sList$shape)),'\n', sep = '')
cat('scale <-',xcat(numclean(sList$scale)),'\n', sep = '')
}
cat('G <- gammaBlur(n, shape, scale)\n')
}
if (sList$blur == 'box.car'){
if (sList$m == 1){
cat('boxWindow <- 1/sqrt(', 1/as.numeric(sList$width)^2,')\n', sep = '')
} else {
cat('boxWindow <- 1/sqrt(', xcat(1/as.numeric(sList$width)^2),')\n', sep = '')
}
cat("G <- boxcarBlur(n, boxWindow)\n")
}
cat('X <- blurSignal(signal, G)\n')
if (input$signalShow == '2'){
cat("matplot(t, X, type = 'l')")
}
if (input$signalShow == '3'){
if (sList$m == 1){
if (any(input$alpha != 1)){
cat('alpha <-', numclean(sList$alpha),'\n')
}
cat('SNR <-', numclean(sList$SNR), '\n')
} else {
if (any(input$alpha != 1)){
cat('alpha <-', xcat(numclean(sList$alpha)),'\n')
}
cat('SNR <-', xcat(numclean(sList$SNR)), '\n')
}
cat('sigma <- sigmaSNR(X, SNR)\n')
cat('E <- multiNoise(n, sigma, ',ifelse(any(input$alpha != 1), paste0('alpha'),),')\n', sep ='')
cat('Y <- X + E\n')
cat("matplot(t, Y, type = 'l')")
}
}
})
output$resolutionCalls <- renderPrint({
cat('Base R Function calls:\n\n')
wList <- mWaveDList()
showG <- ifelse(input$blur != "direct", paste0(",G = G"), paste0(""))
showAlpha <- ifelse(any(input$alpha != 1), paste0(", alpha = alpha"), paste0(""))
showRes <- ifelse((wList$mWaveD$resolution == "block" && detectBlur(wList$mWaveD$G) != "box.car") || (wList$mWaveD$resolution == "smooth" && detectBlur(wList$mWaveD$G) == "box.car"), paste0(", resolution = '",wList$mWaveD$resolutionMethod,"'"),paste0(""))
showShrink <- ifelse(wList$mWaveD$shrinkType != "hard", paste0(", shrinkType = '",wList$mWaveD$shrinkType,"'"),paste0(""))
showDeg <- ifelse(wList$mWaveD$degree != 3, paste0(", deg = ",wList$mWaveD$degree), paste0(""))
cat("mWaveD.output <- multiWaveD(Y", showG, showAlpha, showRes, showShrink, showDeg,')\n', sep = '')
cat('plot(mWaveD.output, which = 3)')
})
output$mWaveDCalls <- renderPrint({
cat('Base R Function calls:\n\n')
wList <- mWaveDList()
showG <- ifelse(input$blur != "direct", paste0(",G = G"), paste0(""))
showAlpha <- ifelse(any(input$alpha != 1), paste0(", alpha = alpha"), paste0(""))
showRes <- ifelse((wList$mWaveD$resolution == "block" && detectBlur(wList$mWaveD$G) != "box.car") || (wList$mWaveD$resolution == "smooth" && detectBlur(wList$mWaveD$G) == "box.car"), paste0(", resolution = '",wList$mWaveD$resolutionMethod,"'"),paste0(""))
showShrink <- ifelse(wList$mWaveD$shrinkType != "hard", paste0(", shrinkType = '",wList$mWaveD$shrinkType,"'"),paste0(""))
showDeg <- ifelse(wList$mWaveD$degree != 3, paste0(", deg = ",wList$mWaveD$degree), paste0(""))
cat('# Note the following commands just give the mWaveD estimate only.\n')
cat("mWaveD.output <- multiWaveD(Y", showG, showAlpha, showRes, showShrink, showDeg,')\n', sep = '')
cat('plot(mWaveD.output, which = 2)\n')
cat('## Or alternatively\n')
cat('x = (0:(n-1))/n\n')
cat("estimate <- multiEstimate(Y", showG, showAlpha, showRes, showShrink, showDeg,')\n', sep = '')
cat("plot(x, estimate, type = 'l')")
})
output$mraCalls <- renderPrint({
cat('Base R Function calls:\n\n')
wList <- mWaveDList()
cat('j1 <- mWaveD.output$j1\n')
cat("rawOutputCoefs <- mWaveD.output$coef\n")
cat("shrinkCoefs <- mWaveD.output$shrinkCoef\n")
cat('plot(rawOutputCoefs, shrinkCoef = shrinkCoefs, highest = j1)\n')
cat('## Or alternatively\n')
cat('plot(mWaveD.output, which = 4)')
})
observe({
# This observer depends on shrinkage1 and updates shrinkage2 with any changes
if (currentShrinkage != input$shrinkage1){
# Then we assume text2 hasn't yet been updated
updateSelectInput(session, inputId = 'shrinkage2', label = NULL, choices = shrinkChoices,
selected = input$shrinkage1)
currentShrinkage <<- input$shrinkage1
}
})
observe({
# This observer depends on shrinkage2 and updates shrinkage1 with any changes
if (currentShrinkage != input$shrinkage2){
# Then we assume text2 hasn't yet been updated
updateSelectInput(session, inputId = 'shrinkage1', label = NULL, choices = shrinkChoices,
selected = input$shrinkage2)
currentShrinkage <<- input$shrinkage2
}
})
updateSelectInput(session, inputId = 'shrinkage1', label = shrinkLabel, choices = shrinkChoices,
selected = currentShrinkage)
updateSelectInput(session, inputId = 'shrinkage2', label = shrinkLabel, choices = shrinkChoices,
selected = currentShrinkage)
output$etaSlider <- renderUI({
sList <- sigList()
deta <- theoreticalEta(alpha = sList$alpha, blur = sList$blur, G = sList$G, sList$sigma)
sliderInput("eta", "Smoothing parameter range:",
min = 0.1, max = 4, value = deta)
})
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.