#Dmel
# dmel_1 <-read.table('/media/sdc/D.melanogaster/rpm_rpkm/exon/Dmel_rep1.txt', header=F,sep='\t');
# dmel_2 <-read.table('/media/sdc/D.melanogaster/rpm_rpkm/exon/Dmel_rep2.txt', header=F,sep='\t');
# dmel_3 <-read.table('/media/sdc/D.melanogaster/rpm_rpkm/exon/Dmel_rep3.txt', header=F,sep='\t');
# dmel_4 <-read.table('/media/sdc/D.melanogaster/rpm_rpkm/exon/Dmel_rep4.txt', header=F,sep='\t');
# dmel_5 <-read.table('/media/sdc/D.melanogaster/rpm_rpkm/exon/Dmel_rep5.txt', header=F,sep='\t');
# cmat <- cbind( as.numeric(dmel_1$V2), as.numeric(dmel_2$V2), as.numeric(dmel_3$V2), as.numeric(dmel_4$V2), as.numeric(dmel_5$V2));
# 
# 
# 
# exon<-dmel_1[,3]
# ecount <- rowSums(cmat)  #total for each exon
# keep <- which(ecount>0)
# ymat <- cmat[keep,]  #exons with count of 1 have no information
# ecount <- ecount[keep]
# lambda <- ecount/sum(ymat)
# 
# lcount <- colSums(ymat)     #total for each lane
# expect <- outer(lambda,lcount)
# 
# rpkm.scale <- 1e9* outer(1/ratExons[keep], 1/lcount)
# rpkm <- ymat*rpkm.scale

agreeFun <- function(lane1, lane2, cuts) {
    bin1 <- cut(lane1, cuts, include.lowest=T)
    bin2 <- cut(lane2, cuts, include.lowest=T)
    table(bin1, bin2)
}

kappaFun <- function(tab) {
    agree <- sum(diag(tab))
    chance <- sum(rowSums(tab)*colSums(tab))/sum(tab)
    kappa <-  (agree-chance) / (sum(tab) - chance)
    c(agree=agree/sum(tab),kappa=kappa)
}

# test12 <- agree(1,2)


# What is the expected table, if the data were truely Poisson?
# etable <- function(lane1, lane2, cuts=c(10,25)){
#     newcut1 <- outer(rpkm.scale[,lane1], cuts)
#     p1 <- cbind(ppois(newcut1, expect[,lane1]),1)
#     for (i in ncol(p1):2) p1[,i] <- p1[,i] - p1[,i-1]
#    
#     newcut2 <- outer(rpkm.scale[,lane2], cuts)
#     p2 <- cbind(ppois(newcut2, expect[,lane2]),1)
#     for (i in ncol(p1):2) p2[,i] <- p2[,i] - p2[,i-1]
# 
#     t(p1) %*% p2
# }
# 
# #correlation coeeficient
# for(i in 1:ncol(cmat)) {
#   for(j in i:ncol(cmat)) {
#     if(i == j){}
#     else{
#       cat("lane ",i," with lane ",j,": ", cor(log10(rpkm[,i]+1), log10(rpkm[,j]+1),method="pearson"),cor(rpkm[,i], rpkm[,j],method="spearman"),"\n"); 
#     }
#   }
# }