0.1 Make plots

f <-  "collman14v2_tightAnnotationCubes20171101T2000.csv.h5"
#h5ls(f)
loc <- h5read(f, name = "Locations")[, c(3,2,1)]
w <- h5read(f, name = "w")
colnames(loc) <- c("x", "y", "z")

## First cleaning the data and removing rows with all zeros
## it seems that NaN's in python get read by R as 0
#tmp <- read.csv("collman14v2_tightAnnotationCubes20171101T1630.csv")
tmp <- read.csv("collman14v2_tightAnnotationCubes20171101T2000.csv")
ze <- apply(tmp, 1, function(x) all(x == 0))
#write.csv(loc[ze,c(3,2,1)], file = "collman14v2_ZeroLocations.csv", row.names = FALSE)
tight <- tmp[!ze,]
loc   <- loc[!ze,]

w <- w[!ze] ## Weight given by number of annotated pixels in cube

for(i in 1:length(w)){
  tight[i,] <- tight[i,]/w[i]
}
  
stight <- scale(tight, center = TRUE, scale = TRUE)

ccol <- c('blue', 'blue', 'blue', 'red', 'red', 
          'red', 'red', 'darkgreen', 'darkgreen', 'darkgreen',  
          'darkgreen')

if(FALSE){
  set.seed(317)
  Lt <- runAll(stight, ccol = ccol)
  Lt[[1]] <- mlocation(tight, ccol = ccol)
  saveRDS(Lt, file = "collman14v2_medaRun_seed317_tightAnnotations.rds")
} else {
  Lt <- readRDS("collman14v2_medaRun_seed317_tightAnnotations.rds")
}


w = 720
h = 720 

png("meda_plots_tight/d1heat.png", width = w, height = h)
p2 <- plot(Lt[[2]]) 
show(p2)
dev.off()

png("meda_plots_tight/mlocation.png", width = w, height = 0.5*h)
p1 <- plot(Lt[[1]]) 
show(p1)
dev.off()

png("meda_plots_tight/cumulativeVariance.png", width = w, height = h)
p3 <- plot(Lt[[3]]) 
show(p3)
dev.off()

png("meda_plots_tight/outliers.png", width = w, height = h)
p4 <- plot(Lt[[4]]) 
show(p4)
dev.off()

png("meda_plots_tight/cor.png", width = w, height = h)
p5 <- plot(Lt[[6]]) 
show(p5)
dev.off()

png("meda_plots_tight/pairhexTight.png", width = 2*w, height = 2*h)
pairhex(stight)
dev.off()


png("meda_plots_tight/hmcClassificationsTight.png", width = 2*w, height = 2*h)
cr <- viridis(max(Lt[[7]]$dat$labels$col))
pairs(Lt[[7]]$dat$data, pch = 19, cex = 0.5, col = cr[Lt[[7]]$dat$labels$col])
dev.off()


png("meda_plots_tight/dendrograms.png", width = w, height = h)
p8 <- plotDend(Lt[[7]])
show(p8)
dev.off()

png("meda_plots_tight/stackMeans.png", width = w, height = 2*h)
p9 <- stackM(Lt[[7]], ccol = ccol, depth = 3, centered = TRUE)
show(p9)
dev.off()

png("meda_plots_tight/clusterMeans.png", width = w, height = h)
p10 <- clusterMeans(Lt[[7]], ccol = ccol)
show(p10)
dev.off()
K <- 5
write.csv(L[[7]]$dat$labels$col, file = "meda_plots_tight/hmcLabels.csv", row.names = FALSE)
c5 <- meda::closestK(Lt[[7]], K = K, locs = loc)
saveRDS(c5, file = "meda_plots_tight/hmcC5.rds")
oc5 <- Reduce(rbind, c5)

row.names(oc5) <- rep(names(c5), each = K)
write.csv(oc5, file ="meda_plots_tight/hmcC5.csv") 
write.csv(oc5, file ="meda_plots_tight/hmcC5synaptogramLocations.csv", row.names=FALSE) 
#f <- 'collman14v2_annotationCubes.csv.h5'
f <- 'collman14v2_fullCubes_20171101T1630.csv.h5'
#h5ls(f)
F0 <- h5read(f, name = "F0")
loc <- h5read(f, name = "Locations")
chan <- h5read(f, name = "Channels")
colnames(F0) <- chan
sF0 <- scale(F0, center = TRUE, scale = TRUE)

ccol <- c('blue', 'blue', 'blue', 'red', 'red', 
          'red', 'red', 'darkgreen', 'darkgreen', 'darkgreen', 
          'darkgreen')


if(FALSE){
  set.seed(317)
  Lc <- runAll(sF0, ccol = ccol)
  Lc[[1]] <- mlocation(F0, ccol = ccol)
  saveRDS(Lc, file = "collman14v2_medaRun_seed317_fullCubesAnnotations.rds")
} else {
  Lc <- readRDS("collman14v2_medaRun_seed317_fullCubesAnnotations.rds")
}


w = 720
h = 720 

png("meda_plots_cubes/d1heat.png", width = w, height = h)
p2 <- plot(Lc[[2]]) 
show(p2)
dev.off()

png("meda_plots_cubes/mlocation.png", width = w, height = 0.5*h)
p1 <- plot(Lc[[1]]) 
show(p1)
dev.off()

png("meda_plots_cubes/cumulativeVariance.png", width = w, height = h)
p3 <- plot(Lc[[3]]) 
show(p3)
dev.off()

png("meda_plots_cubes/outliers.png", width = w, height = h)
p4 <- plot(Lc[[4]]) 
show(p4)
dev.off()

png("meda_plots_cubes/cor.png", width = w, height = h)
plot(Lc[[6]]) 
dev.off()

png("meda_plots_cubes/pairhex.png", width = 2*w, height = 2*h)
pairhex(sF0)
dev.off()


png("meda_plots_cubes/hmcClassifications.png", width = 2*w, height = 2*h)
cr <- viridis(max(Lc[[7]]$dat$labels$col))
pairs(Lc[[7]]$dat$data, pch = 19, cex = 0.5, col = cr[Lc[[7]]$dat$labels$col])
dev.off()


png("meda_plots_cubes/dendrograms.png", width = w, height = h)
plotDend(Lc[[7]])
dev.off()

png("meda_plots_cubes/stackMeans.png", width = w, height = 2*h)
p7 <- stackM(Lc[[7]], ccol = ccol, depth = 3, centered = TRUE)
show(p7)
dev.off()

png("meda_plots_cubes/clusterMeans.png", width = w, height = h)
p8 <- clusterMeans(Lc[[7]], ccol = ccol)
show(p8)
dev.off()

1 Tight F0 plots

1.1 1-d Heatmap

1.2 Location meda_plots

1.3 Outliers as given by randomForest

1.4 Correlation Matrix

1.5 Cumulative Variance with Elbows

1.6 Paired Hex-binned plot

1.7 Hierarchical GMM Classifications

1.8 Hierarchical GMM Dendrogram

1.9 Stacked Means

1.10 Cluster Means

2 Half micron cube F0 plots

2.1 1-d Heatmap

2.2 Location meda_plots

2.3 Outliers as given by randomForest

2.4 Correlation Matrix

2.5 Cumulative Variance with Elbows

2.6 Paired Hex-binned plot

2.7 Hierarchical GMM Classifications

2.8 Hierarchical GMM Dendrogram

2.9 Stacked Means

2.10 Cluster Means