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



#########################
##Blocks
make_fabi_additive_data_blocks_pos<- function(n,l,p,f1,f2,of1,of2,sd_noise,
sd_off,
mean_off,
sd_off_noise,
sd_row,
mean_row,
sd_row_noise,
sd_col,
mean_col,
sd_col_noise
){
	
    if(!is.numeric(n)) {
        stop("n must be numeric")
    }
    if(!is.numeric(l)) {
        stop("l must be numeric")
    }
	
    if(!is.numeric(p)) {
		stop("p must be numeric")
    }
    if(!is.numeric(f1)) {
        stop("f1 must be numeric")
    }
	
	
    za <- floor(runif(p)*(l/f1)+of1)
    la <- floor(runif(p)*(n/f2)+of2)
    ZC <- list()
    LC <- list()
	
	
    X <- matrix(rnorm(l*n,mean = 0, sd = sd_noise),n,l)
    Y <- matrix(0,n,l)
	
    for (i in 1:p){
    	
    	offset <- matrix(rnorm(l*n,mean = 0, sd = sd_off_noise),n,l) + rnorm(1, mean = mean_off, sd = sd_off)
    	
    	row_effect <- matrix(rnorm(l*n,mean = 0, sd = sd_row_noise),n,l) + matrix(rep(rnorm(n, mean = mean_row, sd = sd_row),l), n,l)
    	
    	col_effect <- matrix(rnorm(l*n,mean = 0, sd = sd_col_noise),n,l) + matrix(rep(rnorm(l, mean = mean_col, sd = sd_col),n), n,l, byrow=T)
		
		tf <- floor(runif(1)*(l-za[i]))
        zf <- tf:(tf+za[i]-1)
        ZC[[i]] <- zf
        
        uf <- floor(runif(1)*(n-la[i]))
        lf <- uf:(uf+la[i]-1)
        LC[[i]] <- lf
		
		
        X[lf,zf] <- X[lf,zf] + offset[lf,zf] + row_effect[lf,zf] + col_effect[lf,zf]
        Y[lf,zf] <- Y[lf,zf] + offset[lf,zf] + row_effect[lf,zf] + col_effect[lf,zf]
		
		
    }
    return(list(X=X,Y=Y,ZC=ZC,LC=LC))
}



#########################
## NO Blocks
make_fabi_additive_data <- function(n,l,p,f1,f2,of1,of2,sd_noise,sd_off,mean_off,sd_off_noise,sd_row,mean_row,sd_row_noise,sd_col,mean_col,sd_col_noise){
	
    if(!is.numeric(n)) {
        stop("n must be numeric")
    }
    if(!is.numeric(l)) {
        stop("l must be numeric")
    }
	
    if(!is.numeric(p)) {
        stop("p must be numeric")
    }
    if(!is.numeric(f1)) {
        stop("f1 must be numeric")
    }
	
	
	
    za <- floor(runif(p)*(l/f1)+of1)
    la <- floor(runif(p)*(n/f2)+of2)
    ZC <- list()
    LC <- list()
	
	X <- matrix(rnorm(l*n,mean = 0, sd = sd_noise),n,l)
    Y <- matrix(0,n,l)
	
    for (i in 1:p){
		
		neg <- sign(rnorm(1,mean=1,sd=2))
		
    	offset <- matrix(rnorm(l*n,mean = 0, sd = sd_off_noise),n,l) + rnorm(1, mean = neg * mean_off, sd = sd_off)
    	
    	row_effect <- matrix(rnorm(l*n,mean = 0, sd = sd_row_noise),n,l) + matrix(rep(rnorm(n, mean = mean_row, sd = sd_row),l), n,l)
    	
    	col_effect <- matrix(rnorm(l*n,mean = 0, sd = sd_col_noise),n,l) + matrix(rep(rnorm(l, mean = mean_col, sd = sd_col),n), n,l, byrow=T)
		
#tf <- floor(runif(1)*(l-za[i]))
        zf <- as.integer(runif(za[i], min=1,max=l))
		
        ZC[[i]] <- zf
        
#      uf <- floor(runif(1)*(n-la[i]))
        lf <- as.integer(runif(la[i], min=1,max=n))
		
        LC[[i]] <- lf
		
		
        X[lf,zf] <- X[lf,zf] + offset[lf,zf] + row_effect[lf,zf] + col_effect[lf,zf]
        Y[lf,zf] <- Y[lf,zf] + offset[lf,zf] + row_effect[lf,zf] + col_effect[lf,zf]
		
		
    }
    return(list(X=X,Y=Y,ZC=ZC,LC=LC))
}







convertDat <- function(dat,n=1000,l=100)
{
	
    require(biclust)
	
	
    p <- length(dat$LC)
	
    L <- matrix(0,n,p)
    Z <- matrix(0,p,l)
	
	
    for (i in 1:p) {
		
        for (k in 1:length(dat$LC[[i]])) {
            L[dat$LC[[i]][k],i] <- 1
        }
        for (k in 1:length(dat$ZC[[i]])) {
            Z[i,dat$ZC[[i]][k]] <- 1
        }
    }
	
	
    lp<- list()
    lp[[1]]="make_fabi_data"
    lp[[2]]=n
    lp[[3]]=l
	
    bicB <- BiclustResult(lp,L,Z,p)
	
    return(bicB)
}



convertFabia <- function(rFab,n=1000,l=100,minc=5,minr=30,method="fabia",cyc=200,alpha=0.3,spl=1.0,spz=1.0,p=13,sL=0.0,sZ=0.0,La=NULL,Za=NULL,lapla=NULL,Psi=NULL)
{
	
    require(biclust)
	
	
    p <- length(rFab$bic[,1])
	
    pp <- 0
    for (i in 1:p) {
        if (rFab$bic[i,1]$binp[1]>=minr)
        {
            if (rFab$bic[i,1]$binp[2]>=minc) {
                pp <- pp + 1
            }
        }
    }
	
	if ((p==0)||(pp==0)) {
        L <- matrix(0,n,1)
        Z <- matrix(0,1,l)
        pp <- 0
    } else {
		
		
        L <- matrix(0,n,pp)
        Z <- matrix(0,pp,l)
		
        j <- 0
		
        for (i in 1:p) {
			
            if (rFab$bic[i,1]$binp[1]>=minr)
            {
                if (rFab$bic[i,1]$binp[2]>=minc) {
                    j <- j + 1
                    for (k in 1:rFab$bic[i,1]$binp[2]) {
                        Z[j,rFab$numn[i,2]$numnp[k]] <- 1
                    }
                    for (k in 1:rFab$bic[i,1]$binp[1]) {
                        L[rFab$numn[i,1]$numng[k],j] <- 1
                    }
                }
            }
        }
    }
	
    lp<- list()
    lp[[1]]=method
    lp[[2]]=cyc
    lp[[3]]=alpha
    lp[[4]]=spl
    lp[[5]]=spz
    lp[[6]]=p
    lp[[7]]=sL
    lp[[8]]=sZ
    lp[[9]]=La
    lp[[10]]=Za
    lp[[11]]=lapla
    lp[[12]]=Psi
	
	
	
    bicB <- BiclustResult(lp,L,Z,pp)
	
    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)
# print(c(sa, u))
            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)
                
#print(c(sb, u))
                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)))
}


readPlaidResults <- function(filename,n,l,ab="ab",ss="ss",iter="default") {
	
    require(biclust)
	
    r1 <- readLines(filename)
    a1 <- strsplit(r1," ")
	
    la <- length(a1)
    p <- as.numeric(a1[[la-4]][1])
	
	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)
		
        j <- 0
        r <- 0
		
        for (i in 1:la) {
			
            if (!is.na(a1[[i]][1]))
            {
				
				
                if (a1[[i]][1]=="row")
                {
                    j <- j+1
                    r <- 1
                }
				
                if (a1[[i]][1]=="col")
                {
                    r <- 2
                }
				
                if (a1[[i]][1]==as.character(j))
                {
                    if(r==1) {
                        L[as.numeric(a1[[i]][2]),j] <- 1
                    }
                    if(r==2) {
                        Z[j,as.numeric(a1[[i]][2])] <- 1
                    }
                }
            }
        }
    }
	
    pp<- list()
    pp[[1]]="plaid" # method
    pp[[2]]=ab    # a=alpha(probe) effects; b= beta(conditions) effects
    pp[[3]]=ss    # seek=ss prefer large, possibly diffuse layers
# seek=ms prefer intense, possibly small layers
    pp[[4]]=iter # how many seeking iterations
	
    bicB <- BiclustResult(pp,L,Z,p)
	
    return(bicB)
}


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


readSambaResults <- function(filename,n,l,pre="standardization",opt="valsp_3ap",over="0.5") {
	
    require(biclust)
	
    r1 <- readLines(filename)
    a1 <- strsplit(r1,"\t")
	
    la <- length(a1)
	
    s <- 1
    while ((a1[[s]][1]!="[Bicd]")&&(s<=la)) { s <- s+1 }
	
    p <- as.numeric(a1[[s-1]][1])+1
	
	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)
		
        s <- s+1
        for (i in s:la) {
			
            if (a1[[i]][2]=="1") {
                L[as.numeric(a1[[i]][3]),(as.numeric(a1[[i]][1])+1)] <- 1
            }
			
            if (a1[[i]][2]=="0") {
                Z[(as.numeric(a1[[i]][1])+1),as.numeric(a1[[i]][3])] <- 1
            }
        }
    }
	
    pp<- list()
    pp[[1]]="samba" # method
    pp[[2]]=pre     # preprocessing: standardization
    pp[[3]]=opt     # options: valsp_3ap
    pp[[4]]=over    # overlap genes 0.1 or 0.5
	
    bicB <- BiclustResult(pp,L,Z,p)
	
    return(bicB)
}