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)))
}

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

    pA <- bicA@Number
    pB <- bicB@Number
    if ((pA>0)&&(pB>0)) {
        l <- length(bicA@NumberxCol[1,])
        u <- 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 <- bicA@NumberxCol[i,]
            apos <- bcA > 0
            sa <- sum(apos)
            if (sa > 0.8*u) {
                 bcA[apos] <- 0
                 sa <- 0
            }
            for (j in 1:pB) {
                bcB <- bicB@NumberxCol[j,]
                bpos <- bcB > 0
                sb <- sum(bpos)
                if (sb > 0.8*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)))
}



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

alpha_fabia = 0.1
alpha_fabias = 0.4
spl=1.0
spz=1.0
cyc=200
p=5


#=====================================================================
#=====================================================================
    data(Breast_A)

    l <- ncol(XBreast)
    n <- nrow(XBreast)

    porg <- 3


    load("exp_breast_Biclust_orig.RData")

    X <- as.matrix(XBreast)
    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("exp_breast_isa_2.res",n,l,method="isa",pre="standardization",tc="1.1",tg="0.7")

    save(risa_2,file="exp_breast_Biclust_isa_2.RData")

    rspec_1 <- biclust(exp(X), method=BCSpectral(), normalization="log", numberOfEigenvalues=1, minr=30, minc=5, withinVar=100)
    save(rspec_1,file="exp_breast_Biclust_spec_1.RData")

    rspec_2 <- biclust(exp(X), method=BCSpectral(), normalization="log", numberOfEigenvalues=3, minr=30, minc=5, withinVar=20)
    save(rspec_2,file="exp_breast_Biclust_spec_2.RData")

    indrisa_2<-col_indices(breast_biclust,risa_2)

    indrspec_1<-col_indices(breast_biclust,rspec_1)
    indrspec_2<-col_indices(breast_biclust,rspec_2)

    save(indrisa_2,indrspec_1,indrspec_2,file="exp_breast_ind_res_1.RData")

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


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

    data(Multi_A)

    l <- ncol(XMulti)
    n <- nrow(XMulti)

    porg <- 4

    load("exp_multi_Biclust_orig.RData")

    X <- as.matrix(XMulti)
    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






    rspec_1 <- biclust(exp(X), method=BCSpectral(), normalization="log", numberOfEigenvalues=1, minr=30, minc=5, withinVar=100)
    save(rspec_1,file="exp_multi_Biclust_spec_1.RData")

    rspec_2 <- biclust(exp(X), method=BCSpectral(), normalization="log", numberOfEigenvalues=3, minr=30, minc=5, withinVar=20)
    save(rspec_2,file="exp_multi_Biclust_spec_2.RData")


    indrspec_1<-col_indices(multi_biclust,rspec_1)
    indrspec_2<-col_indices(multi_biclust,rspec_2)


    save(indrspec_1,indrspec_1,file="exp_multi_ind_res_1.RData")

     write.table(t(as.vector(indrspec_1)), file = "exp_multi_in_res_1.txt", quote = FALSE, sep = "\t",row.names = FALSE,col.names = FALSE,append =FALSE)
     write.table(t(as.vector(indrspec_2)), file = "exp_multi_in_res_1.txt", quote = FALSE, sep = "\t",row.names = FALSE,col.names = FALSE,append =TRUE)

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

    data(DLBCL_B)

    l <- ncol(XDLBCL)
    n <- nrow(XDLBCL)

    porg <- 3

    load("exp_dlbcl_Biclust_orig.RData")

    X <- as.matrix(XDLBCL)
    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("exp_dlbcl_isa_2.res",n,l,method="isa",pre="standardization",tc="1.1",tg="0.7")

    save(risa_2,file="exp_dlbcl_Biclust_isa_2.RData")

    rspec_1 <- biclust(exp(X), method=BCSpectral(), normalization="log", numberOfEigenvalues=1, minr=30, minc=5, withinVar=100)
    save(rspec_1,file="exp_dlbcl_Biclust_spec_1.RData")
    rspec_2 <- biclust(exp(X), method=BCSpectral(), normalization="log", numberOfEigenvalues=3, minr=30, minc=5, withinVar=20)
    save(rspec_2,file="exp_dlbcl_Biclust_spec_2.RData")

    indrisa_2<-col_indices(dlbcl_biclust,risa_2)
    indrspec_1<-col_indices(dlbcl_biclust,rspec_1)
    indrspec_2<-col_indices(dlbcl_biclust,rspec_2)

    save(indrisa_2,indrspec_1,indrspec_2,file="exp_dlbcl_ind_res_1.RData")

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