###################################################
### chunk number 1: stage0
###################################################
library(flowMerge)


###################################################
### chunk number 2: stage1
###################################################
data(rituximab)
summary(rituximab)
flowClust.res <- flowClust(rituximab, varNames=c("FSC.H", "SSC.H"), K=1:5, B=1000,B.init=100,tol=1e-5,tol.init=1e-2,randomStart=50,nu=4,nu.est=1,trans=1);


###################################################
### chunk number 3: stage2
###################################################
plot(flowMerge:::BIC(flowClust.res),
main="BIC for 1 through 10 cluster flowClust solutions",xlab="K",ylab="BIC",type="o");
flowClust.maxBIC<-flowClust.res[[which.max(flowMerge:::BIC(flowClust.res))]];


###################################################
### chunk number 4: stage3
###################################################
flowClust.flowobj<-flowObj(flowClust.maxBIC,rituximab);
flowClust.merge<-merge(flowClust.flowobj);
i<-fitPiecewiseLinreg(flowClust.merge,plot=T);


###################################################
### chunk number 5: stage4
###################################################
par(mfrow=c(2,2));
flowClust.mergeopt<-flowClust.merge[[i]];
plot(flowClust.res[[4]],data=rituximab,main="Max BIC solution");
plot(flowClust.res[[which.max(flowMerge:::ICL(flowClust.res))]],data=rituximab,main="Max ICL solution");
plot(flowClust.mergeopt,level=0.9,pch=20,main="Merged Solution");


###################################################
### chunk number 6: stage5
###################################################
pop<-which(apply(apply(getEstimates(flowClust.mergeopt)$locations,2,function(x)order(x)==2),1,all));
lymphocytes<-split(getData(flowClust.mergeopt),flowClust.mergeopt,population=list("lymphocytes"=pop))$lymphocytes;
lymphocytes<-lymphocytes[,c(3,5)];
l.flowC<-flowClust(lymphocytes,varNames=c("FL1.H","FL3.H"),K=1:8,B=1000,B.init=100,tol=1e-5,tol.init=1e-2,randomStart=50,nu=4,nu.est=1,trans=1);


###################################################
### chunk number 7: stage6
###################################################
par(mfrow=c(2,2));
l.flowO<-flowObj(l.flowC[[which.max(flowMerge:::BIC(l.flowC))]],lymphocytes);
plot(l.flowO,main="max BIC solution",new.window=F,pch=20,level=0.9);
plot(flowObj(l.flowC[[which.max(flowMerge:::ICL(l.flowC))]],lymphocytes),main="max ICL solution",new.window=F,pch=20,level=0.9);
l.flowM<-merge(l.flowO);
i<-fitPiecewiseLinreg(l.flowM,plot=T);
plot(l.flowM[[i]],new.window=F,main="Best Merged Solution",pch=20,level=0.9);