###################################################
### chunk number 1: getGolub
###################################################
library(edd)
library(golubEsets)
library(xtable)
data(Golub_Merge)


###################################################
### chunk number 2: filterGolub
###################################################
madvec <- apply(exprs(Golub_Merge),1,mad)
minvec <- apply(exprs(Golub_Merge),1,min)
keep <- (madvec > median(madvec)) & (minvec > 300)
gmfilt <- Golub_Merge[keep==TRUE,]


###################################################
### chunk number 3: splitGolub
###################################################
ALL <- gmfilt$ALL.AML=="ALL"
gall <- gmfilt[,ALL==TRUE]
gaml <- gmfilt[,ALL==FALSE]
show(gall)


###################################################
### chunk number 4: eddByType
###################################################
set.seed(12345)
alldists <- edd(gall, meth="nnet", size=10, decay=.2)
amldists <- edd(gaml, meth="nnet", size=10, decay=.2)


###################################################
### chunk number 5: 
###################################################
greo <- match(sapply(eddDistList,tag),names(table(alldists)))
if (any(ii <- is.na(greo))) greo = greo[-which(ii)]


###################################################
### chunk number 6: 
###################################################
gn <- row.names(exprs(gall))
alld <- as.character(alldists)[1:5]
names(alld) <- gn[1:5]
print(alld)


###################################################
### chunk number 7: lookKNN
###################################################
set.seed(123)
alldistsKNN <- edd(gall, meth="knn", k=1, l=0)
alldistsTEST <- edd(gall, meth="test", thresh=.3)


###################################################
### chunk number 8: compareKNNvNN
###################################################
cap <- "Comparison of distribution shape classification by nnet (rows) and by knn (columns) methods in edd."
print(try(xtable(
   latEDtable(table(alldists,alldistsKNN),reorder=greo),digits=rep(0,length(table(alldists))+1),caption=cap, label="conc1")))


###################################################
### chunk number 9: testTable
###################################################
print(table(alldistsTEST))


###################################################
### chunk number 10: ALLtable
###################################################
cap <- "Frequencies of distributional shapes in filtered ALL data."
print(xtable(latEDtable(
    table(alldists),reorder=greo),digits=rep(0,length(table(alldists))+1),caption=cap,label="marg1"))


###################################################
### chunk number 11: compare
###################################################
mymat <- matrix(c(1,2),nr=1)
layout(mymat)#,heights=c(lcm(6),lcm(6)),widths=c(lcm(6),lcm(6)),respect=TRUE)
barplot(table(alldists),las=2)
barplot(table(amldists),las=2)


###################################################
### chunk number 12: discordTable
###################################################
cap <- "Rows are gene-specific distribution shapes for ALL, columns for AML, and cell entries are counts of genes."
print(xtable(latEDtable(table(alldists,amldists),reord=greo),cap=cap,
 label="disco1"))


###################################################
### chunk number 13: findBiMod
###################################################
print((1:540)[alldists==".75N(0,1)+.25N(4,1)" &
  amldists=="N(0,1)"][1:5])


###################################################
### chunk number 14: lookD87593
###################################################
par(mfrow=c(2,2))
 plotED(MIXN1,data=exprs(gall)[65,])
 hist(exprs(gall)[65,],xlim=c(0,3500))
 plotED(N01,data=exprs(gaml)[65,])
 hist(exprs(gaml)[65,],xlim=c(0,3500))


###################################################
### chunk number 15: namesEDL
###################################################
names(eddDistList)


###################################################
### chunk number 16: plotRefs
###################################################
par(mfrow=c(4,2))
for (i in 1:8)
 plotED(eddDistList[[i]])


###################################################
### chunk number 17: addMIXN3
###################################################
#> median(rmixnorm(10000,.75,0,1,6,1))
#[1] 0.4337298
#> mad(rmixnorm(10000,.75,0,1,6,1))
#[1] 1.550768  -- these are resistant stats!
MIXN3 <- new("eddDist", stub="mixnorm", parms=c(p1=.75,m1=0,s1=1,m2=6,s2=1),
     median=.43, mad=1.55, tag=".75N(0,1)+.25N(6,1)", plotlim=c(-3,11),
     latexTag="$\\frac{3}{4}\\Phi+\\frac{1}{4}\\Phi_{6,1}$")
eddDistList[["MIXN3"]] <- MIXN3
set.seed(12345)
alldists2 <- edd(gall, meth="nnet", size=10, decay=.2)
print(alldists2[65])


###################################################
### chunk number 18: checkFit3
###################################################
plotED(MIXN3,data=exprs(gall)[65,])


###################################################
### chunk number 19: checkFit1
###################################################
plotED(MIXN1,data=exprs(gall)[65,])