setwd("C:/sepp/work/fabia/experiments")


library(fabia)
library(biclust)
library(truecluster)
library(BicARE)


readBicatResults <- function(filename,n,l,method="isa",pre="standardization",tc="2.0",tg="2.0") {

    require(biclust)

    r1 <- readLines(filename)
    a1 <- strsplit(r1," ")

    la <- length(a1)
    p <- la/3

     if ((p==0)||(la==0)) {
        L <- matrix(0,n,1)
        Z <- matrix(0,1,l)
        p <- 0
    } else {

        L <- matrix(0,n,p)
        Z <- matrix(0,p,l)

        for (i in 1:p) {

            nn <- as.numeric(a1[[3*i-2]][1])
            nl <- as.numeric(a1[[3*i-2]][2])
            for (k in 1:nn) {
                L[as.numeric(a1[[3*i-1]][k]),i] <- 1
            }
            for (k in 1:nl) {
                Z[i,as.numeric(a1[[3*i]][k])] <- 1
            }
        }
    }

    pp<- list()
    pp[[1]]=method # method
    pp[[2]]=pre    # preprocessing: standardization
    pp[[3]]=tg     # gene score threshold
    pp[[4]]=tc     # condition score threshold

    bicB <- BiclustResult(pp,L,Z,p)

    return(bicB)
}



indices<- function (bicA, bicB)
{
    require(truecluster)

    pA <- bicA@Number
    pB <- bicB@Number
    if ((pA>0)&&(pB>0)) {
        n <- length(bicA@RowxNumber[,1])
        l <- length(bicA@NumberxCol[1,])
        u <- n*l

        jamat <- matrix(0,pA,pB)
        kumat <- matrix(0,pA,pB)
        ocmat <- matrix(0,pA,pB)
        somat <- matrix(0,pA,pB)
        for (i in 1:pA) {
            bcA <- tcrossprod(bicA@RowxNumber[,i],bicA@NumberxCol[i,])
            apos <- bcA > 0
            sa <- sum(apos)
            if (sa > 0.5*u) {
                 bcA[apos] <- 0
                 sa <- 0
            }
            for (j in 1:pB) {
                bcB <- tcrossprod(bicB@RowxNumber[,j],bicB@NumberxCol[j,])
                bpos <- bcB > 0
                sb <- sum(bpos)
                if (sb > 0.5*u) {
                    bcB[bpos] <- 0
                    sb <- 0
                }
                bcAB <- bcA + bcB
                abpos <- bcAB > 0
                sab <- sum(abpos)

                if (sab>0) {
                    jamat[i,j] <- (sa + sb)/sab - 1
                    somat[i,j] <- 2.0-2.0*sab/(sa+sb)
                } else {
                    jamat[i,j] <- 0
                    somat[i,j] <- 0
                }
                if ((sa>0)&&(sb>0))
                {
                    kumat[i,j] <- 1.0+0.5*( (sa-sab)/sb + (sb -sab)/sa )
                    ocmat[i,j] <- (sa+sb-sab)/sqrt(sb*sa)
                }else {
                    kumat[i,j] <- 0
                    ocmat[i,j] <- 0
                }
            }
        }

        fac <- sqrt(pB*pA)

        mm <- min(pA,pB)
        ma <- max(pA,pB)

        sja <- sum(jamat)/fac
        sku <- sum(kumat)/fac
        soc <- sum(ocmat)/fac
        sso <- sum(somat)/fac

        indja <- munkres(-jamat,tieorder = FALSE)
        indku <- munkres(-kumat,tieorder = FALSE)
        indoc <- munkres(-ocmat,tieorder = FALSE)
        indso <- munkres(-somat,tieorder = FALSE)

        rjat <-  sum(diag(jamat[indja$row, indja$col]))
        rkut <-  sum(diag(kumat[indku$row, indku$col]))
        roct <-  sum(diag(ocmat[indoc$row, indoc$col]))
        rsot <-  sum(diag(somat[indso$row, indso$col]))

        rja <-  rjat/mm
        rku <-  rkut/mm
        roc <-  roct/mm
        rso <-  rsot/mm

        rja1 <-  rjat/ma
        rku1 <-  rkut/ma
        roc1 <-  roct/ma
        rso1 <-  rsot/ma
    } else {
        rja <- 0
        rku <- 0
        roc <- 0
        rso <- 0

        rja1 <- 0
        rku1 <- 0
        roc1 <- 0
        rso1 <- 0

        sja <- 0
        sku <- 0
        soc <- 0
        sso <- 0

    }



    return(as.vector(c(rja,rku,roc,rso,rja1,rku1,roc1,rso1,sja,sku,soc,sso)))
}




indrisa_2_s<-list()
indrspec_1_s<-list()
indrspec_2_s<-list()

    sindrisa_2_s<- as.vector(rep(0,12))
    sindrspec_1_s<-  as.vector(rep(0,12))
    sindrspec_2_s<-  as.vector(rep(0,12))

#=====================================================================
#=====================================================================

m <- 100

n=1000
l=100
p=13

#=====================================================================
#=====================================================================

for (i in 1:m) {


    load(file=paste("exp_",as.character(i),".RData",sep = ""))

    load(file=paste("exp_",as.character(i),"_Biclust_orig.RData",sep = ""))

    X <- dat[[1]]
    X <- X- rowMeans(X)
    XX <- (1/ncol(X))*tcrossprod(X)
    dXX <- 1/sqrt(diag(XX)+0.001*as.vector(rep(1,nrow(X))))
    X <- dXX*X


    risa_2 <- readBicatResults(paste("exp_",as.character(i),"_isa_2.res",sep = ""),1000,100,method="isa",pre="standardization",tc="1.1",tg="0.7")

    save(risa_2,file=paste("exp_",as.character(i),"_Biclust_isa_2.RData",sep = ""))


    rspec_1<- biclust(exp(X), method=BCSpectral(), normalization="log", numberOfEigenvalues=1, minr=30, minc=5, withinVar=100)
    save(rspec_1,file=paste("exp_",as.character(i),"_Biclust_spec_1.RData",sep = ""))


    rspec_2<- biclust(exp(X), method=BCSpectral(), normalization="log", numberOfEigenvalues=3, minr=30, minc=5, withinVar=20)
    save(rspec_2,file=paste("exp_",as.character(i),"_Biclust_spec_2.RData",sep = ""))

    indrisa_2<-indices(rorig,risa_2)
    indrspec_1<-indices(rorig,rspec_1)
    indrspec_2<-indices(rorig,rspec_2)


    save(indrisa_2,indrspec_1,indrspec_2,file=paste("exp_",as.character(i),"_ind_res_1.RData",sep = ""))

     write.table(t(as.vector(indrisa_2)), file = paste("exp_",as.character(i),"_in_res_1.txt",sep = ""), quote = FALSE, sep = "\t",row.names = FALSE,col.names = FALSE,append =FALSE)
     write.table(t(as.vector(indrspec_1)), file = paste("exp_",as.character(i),"_in_res_1.txt",sep = ""), quote = FALSE, sep = "\t",row.names = FALSE,col.names = FALSE,append =TRUE)
     write.table(t(as.vector(indrspec_2)), file = paste("exp_",as.character(i),"_in_res_1.txt",sep = ""), quote = FALSE, sep = "\t",row.names = FALSE,col.names = FALSE,append =TRUE)

    indrisa_2_s[[i]]<-indrisa_2
    indrspec_1_s[[i]]<-indrspec_1
    indrspec_2_s[[i]]<-indrspec_2

    sindrisa_2_s<-sindrisa_2_s  + indrisa_2
    sindrspec_1_s<- sindrspec_1_s + indrspec_1
    sindrspec_2_s<- sindrspec_2_s + indrspec_2

    rm(dat)
    rm(X)
    rm(XX)



}

    sindrisa_2_s<-sindrisa_2_s/m
    sindrspec_1_s<- sindrspec_1_s/m
    sindrspec_2_s<- sindrspec_2_s/m


    sindrisa_2_v<- as.vector(rep(0,12))
    sindrspec_1_v<-  as.vector(rep(0,12))
    sindrspec_2_v<-  as.vector(rep(0,12))

    for (i in 1:m) {


    sindrisa_2_v<-(indrisa_2_s[[i]]-sindrisa_2_s)^2
    sindrspec_1_v<-(indrspec_1_s[[i]]-sindrspec_1_s)^2
    sindrspec_2_v<-(indrspec_2_s[[i]]-sindrspec_2_s)^2


}

    sindrisa_2_v<-sindrisa_2_v/(m-1)
    sindrspec_1_v<- sindrspec_1_v/(m-1)
    sindrspec_2_v<- sindrspec_2_v/(m-1)

     write.table(t(as.vector(sindrisa_2_s)), file ="final_result_1.txt", quote = FALSE, sep = "\t",row.names = FALSE,col.names = FALSE,append =FALSE)
     write.table(t(as.vector(sindrspec_1_s)), file ="final_result_1.txt", quote = FALSE, sep = "\t",row.names = FALSE,col.names = FALSE,append =TRUE)
     write.table(t(as.vector(sindrspec_2_s)), file ="final_result_1.txt", quote = FALSE, sep = "\t",row.names = FALSE,col.names = FALSE,append =TRUE)



     write.table(t(as.vector(sindrisa_2_v)), file ="final_result_1.txt", quote = FALSE, sep = "\t",row.names = FALSE,col.names = FALSE,append =TRUE)
     write.table(t(as.vector(sindrspec_1_v)), file ="final_result_1.txt", quote = FALSE, sep = "\t",row.names = FALSE,col.names = FALSE,append =TRUE)
     write.table(t(as.vector(sindrspec_2_v)), file ="final_result_1.txt", quote = FALSE, sep = "\t",row.names = FALSE,col.names = FALSE,append =TRUE)