write_expander_solution<-function(bics,name,data,minRows,minCols,sample_specific=F){ filename = paste(name,'_sol.txt',sep=""); numBics = length(bics[[1]]) write('[Bick]',file=filename) counter=1 for (i in 1:numBics){ conds = bics[[2]][[i]] genes = bics[[1]][[i]] if (length(conds)1 && C[i,j-1]==1){curr_pi = pi_1_to_1} if (j>1 && C[i,j-1]==0){curr_pi = pi_0_to_1} samp = rbinom(1,1,curr_pi) if (samp==1){ C[i,j] = 1 bic_cols = c(bic_cols,paste("subj",i,";TR",j,sep="")) } } } bic_rows = rownames(m)[1:n_v_bic] H = rep(0,n_voxels);names(H) = rownames(m) H[bic_rows] = 1 HS = matrix (nrow = nrow(m),ncol=n_subj);HS[,]=0 rownames(HS) = names(H); colnames(HS) = rownames(C) # set HS if (p_s == 1 && p_s_c == 0){ HS[bic_rows,] = 1 m[bic_rows,bic_cols] = 1 } else{ for (j in 1:ncol(HS)){ for (i in 1:nrow(HS)){ if (H[i] == 1 && rbinom(1,1,p_s)==1){HS[i,j]=1} if (H[i] == 0 && rbinom(1,1,p_s_c)==1){HS[i,j]=1} } curr_subj_trs = colnames(m)[grepl(pattern = paste("subj",j,";",sep=''),colnames(m))] curr_subj_trs = curr_subj_trs[C[j,]==1] m[HS[,j]==1,curr_subj_trs]=1 } } if (setAsGarbage){ m[,]=0 H[]=0 HS[,] = 0 C[,] = 0 } # add noise for (i in 1:nrow(m)){ for (j in 1:ncol(m)){ currp = out_bic_noise if (m[i,j]==1){ currp = in_bic_noise } samp = rbinom(1,1,currp) if (samp==1){m[i,j]=1-m[i,j]} } } # run the algorithm ptm <- proc.time() gibbs_solution = one_musbic_step(m,...) iters = gibbs_solution$iters gibbs_solution = gibbs_solution$gibbs_solution running_time = (proc.time() - ptm)[3] # get some Jaccard scores cols_j = length(intersect(gibbs_solution$bic_cols,bic_cols))/length(union(gibbs_solution$bic_cols,bic_cols)) rows_j = length(intersect(gibbs_solution$bic_rows,bic_rows))/length(union(gibbs_solution$bic_rows,bic_rows)) C_j = binary_matrix_jaccard(C,gibbs_solution$C) HS_j = binary_matrix_jaccard (HS,gibbs_solution$HS) H_j = binary_matrix_jaccard (H,gibbs_solution$H) H_fpr = binary_FPR(gibbs_solution$H,H) HS_fpr = binary_FPR(gibbs_solution$HS,HS) C_fpr = binary_FPR(gibbs_solution$C,C) H_tpr = binary_TPR(gibbs_solution$H,H) HS_tpr = binary_TPR(gibbs_solution$HS,HS) C_tpr = binary_TPR(gibbs_solution$C,C) H_pr = binary_precision(gibbs_solution$H,H) HS_pr = binary_precision(gibbs_solution$HS,HS) C_pr = binary_precision(gibbs_solution$C,C) return (list(H_j=H_j,C_j=C_j,HS_j=HS_j, H_fpr=H_fpr,C_fpr=C_fpr,HS_fpr=HS_fpr, H_tpr=H_tpr,C_tpr=C_tpr,HS_tpr=HS_tpr, H_pr=H_pr,C_pr=C_pr,HS_pr=HS_pr, running_time=running_time,gibbs_solution=gibbs_solution,iters=iters,I=m,C=C,HS=HS)) } run_simulation_normal<-function(n_voxels = 500, n_subj = 2, n_tr = 50, n_tr_bic = 10, n_v_bic = 20, pi_1_to_1 = 0.6, pi_0_to_1 = 0.1, pi_1_1 = 0.05, p_s_c = 0, p_s = 1, in_bic_noise = 0.1, out_bic_noise = 0.05,setAsGarbage = F,...){ m = matrix(nrow=n_voxels,ncol = n_subj*n_tr);m[,] = 0 rownames(m) = paste ("v",1:n_voxels,sep="_") colnames(m) = rep("",ncol(m)) C = matrix(nrow = n_subj,ncol = n_tr);C[,] = 0 colnames(C) = rep("",ncol(C));rownames(C) = rep("",nrow(C)) count = 1 for(i in 1:n_subj){ for (j in 1:n_tr){ colnames(m)[count] = paste("subj",i,";TR",j,sep="") colnames(C)[j] = paste("TR",j,sep="") rownames(C)[i] = paste("subj",i,sep="") count = count+1 } } # currently assume that w is 1 bic_cols = c() for(i in 1:n_subj){ for (j in 1:n_tr){ curr_pi = pi_1_1 if (j>1 && C[i,j-1]==1){curr_pi = pi_1_to_1} if (j>1 && C[i,j-1]==0){curr_pi = pi_0_to_1} samp = rbinom(1,1,curr_pi) if (samp==1){ C[i,j] = 1 bic_cols = c(bic_cols,paste("subj",i,";TR",j,sep="")) } } } bic_rows = rownames(m)[1:n_v_bic] H = rep(0,n_voxels);names(H) = rownames(m) H[bic_rows] = 1 HS = matrix (nrow = nrow(m),ncol=n_subj);HS[,]=0 rownames(HS) = names(H); colnames(HS) = rownames(C) # set HS if (p_s == 1 && p_s_c == 0){ HS[bic_rows,] = 1 m[bic_rows,bic_cols] = 1 } else{ for (j in 1:ncol(HS)){ for (i in 1:nrow(HS)){ if (H[i] == 1 && rbinom(1,1,p_s)==1){HS[i,j]=1} if (H[i] == 0 && rbinom(1,1,p_s_c)==1){HS[i,j]=1} } curr_subj_trs = colnames(m)[grepl(pattern = paste("subj",j,";",sep=''),colnames(m))] curr_subj_trs = curr_subj_trs[C[j,]==1] m[HS[,j]==1,curr_subj_trs]=1 } } if (setAsGarbage){ m[,]=0 H[]=0 HS[,] = 0 C[,] = 0 } # add noise for (i in 1:nrow(m)){ for (j in 1:ncol(m)){ currp = out_bic_noise if (m[i,j]==1){ currp = in_bic_noise } samp = rnorm(1,mean=0,sd=currp) m[i,j] = m[i,j] + samp } } # run the algorithm ptm <- proc.time() # tests # gibbs_solution = one_musbic_step(m,gibbs_reps=50,initial_method='bimax',method='norm') gibbs_solution = one_musbic_step(m,...) iters = gibbs_solution$iters gibbs_solution = gibbs_solution$gibbs_solution running_time = (proc.time() - ptm)[3] # get some Jaccard scores cols_j = length(intersect(gibbs_solution$bic_cols,bic_cols))/length(union(gibbs_solution$bic_cols,bic_cols)) rows_j = length(intersect(gibbs_solution$bic_rows,bic_rows))/length(union(gibbs_solution$bic_rows,bic_rows)) C_j = binary_matrix_jaccard(C,gibbs_solution$C) HS_j = binary_matrix_jaccard (HS,gibbs_solution$HS) H_j = binary_matrix_jaccard (H,gibbs_solution$H) H_fpr = binary_FPR(gibbs_solution$H,H) HS_fpr = binary_FPR(gibbs_solution$HS,HS) C_fpr = binary_FPR(gibbs_solution$C,C) H_tpr = binary_TPR(gibbs_solution$H,H) HS_tpr = binary_TPR(gibbs_solution$HS,HS) C_tpr = binary_TPR(gibbs_solution$C,C) H_pr = binary_precision(gibbs_solution$H,H) HS_pr = binary_precision(gibbs_solution$HS,HS) C_pr = binary_precision(gibbs_solution$C,C) return (list(H_j=H_j,C_j=C_j,HS_j=HS_j, H_fpr=H_fpr,C_fpr=C_fpr,HS_fpr=HS_fpr, H_tpr=H_tpr,C_tpr=C_tpr,HS_tpr=HS_tpr, H_pr=H_pr,C_pr=C_pr,HS_pr=HS_pr, running_time=running_time,gibbs_solution=gibbs_solution,iters=iters,I=m,C=C,HS=HS)) } run_simulation_multiple_bics<-function(n_bics = 2, wrapper = run_complete_musbic_algorithm, n_voxels = 500, n_subj = 2, n_tr = 50, n_tr_bic = 10, n_v_bic = 20, pi_1_to_1 = 0.6, pi_0_to_1 = 0.1, pi_1_1 = 0.05, p_s_c = 0, p_s = 1, in_bic_noise = 0.1, out_bic_noise = 0.05,setAsGarbage = F,...){ m = matrix(nrow=n_voxels,ncol = n_subj*n_tr);m[,] = 0 rownames(m) = paste ("v",1:n_voxels,sep="_") colnames(m) = rep("",ncol(m)) Cs = list();cols = list() for (k in 1:n_bics){ C = matrix(nrow = n_subj,ncol = n_tr);C[,] = 0 colnames(C) = rep("",ncol(C));rownames(C) = rep("",nrow(C)) count = 1 for(i in 1:n_subj){ for (j in 1:n_tr){ colnames(m)[count] = paste("subj",i,";TR",j,sep="") colnames(C)[j] = paste("TR",j,sep="") rownames(C)[i] = paste("subj",i,sep="") count = count+1 } } # currently assume that w is 1 bic_cols = c() for(i in 1:n_subj){ for (j in 1:n_tr){ curr_pi = pi_1_1 if (j>1 && C[i,j-1]==1){curr_pi = pi_1_to_1} if (j>1 && C[i,j-1]==0){curr_pi = pi_0_to_1} samp = rbinom(1,1,curr_pi) if (samp==1){ C[i,j] = 1 bic_cols = c(bic_cols,paste("subj",i,";TR",j,sep="")) } } } Cs[[k]] = C; cols[[k]] = bic_cols } Hs = list(); HSs = list(); rows = list() for (k in 1:n_bics){ C = Cs[[k]]; bic_cols = cols[[k]] bic_rows = sample(rownames(m))[1:n_v_bic] H = rep(0,n_voxels);names(H) = rownames(m) H[bic_rows] = 1 HS = matrix (nrow = nrow(m),ncol=n_subj);HS[,]=0 rownames(HS) = names(H); colnames(HS) = rownames(C) # set HS if (p_s == 1 && p_s_c == 0){ HS[bic_rows,] = 1 m[bic_rows,bic_cols] = 1 } else{ for (j in 1:ncol(HS)){ for (i in 1:nrow(HS)){ if (H[i] == 1 && rbinom(1,1,p_s)==1){HS[i,j]=1} if (H[i] == 0 && rbinom(1,1,p_s_c)==1){HS[i,j]=1} } curr_subj_trs = colnames(m)[grepl(pattern = paste("subj",j,";",sep=''),colnames(m))] curr_subj_trs = curr_subj_trs[C[j,]==1] m[HS[,j]==1,curr_subj_trs]=1 } } if (setAsGarbage){ m[,]=0;H[]=0;HS[,] = 0;C[,] = 0 } Hs[[k]] = H; HSs[[k]] = HS; rows[[k]] = bic_rows } # add noise to the data for (i in 1:nrow(m)){ for (j in 1:ncol(m)){ currp = out_bic_noise if (m[i,j]==1){currp = in_bic_noise} samp = rbinom(1,1,currp) if (samp==1){m[i,j]=1-m[i,j]} } } # run the algorithm ptm <- proc.time() solution = wrapper(m,...) running_time = (proc.time() - ptm)[3] # get some Jaccard scores max_j_matrix_H = matrix(nrow = n_bics,ncol = length(solution[[1]])) max_j_matrix_H[,] = 0 max_j_matrix_HS = matrix(nrow = n_bics,ncol = length(solution[[1]])) max_j_matrix_HS[,] = 0 max_j_matrix_C = matrix(nrow = n_bics,ncol = length(solution[[1]])) max_j_matrix_C[,] = 0 for (i in 1:nrow(max_j_matrix_H)){ if (ncol(max_j_matrix_H)==0){next} for (j in 1:ncol(max_j_matrix_H)){ gibbs_sol = extract_gibbs_solution(m,solution[[4]][[j]]) C_j = binary_matrix_jaccard(Cs[[i]],gibbs_sol$C) HS_j = binary_matrix_jaccard (HSs[[i]],gibbs_sol$HS) H_j = binary_matrix_jaccard (Hs[[i]],gibbs_sol$H) # update max_j_matrix_H[i,j] = max(max_j_matrix_H[i,j],H_j) max_j_matrix_HS[i,j] = max(max_j_matrix_HS[i,j],HS_j) max_j_matrix_C[i,j] = max(max_j_matrix_C[i,j],C_j) } } return (list(max_j_matrix_H = max_j_matrix_H, max_j_matrix_HS = max_j_matrix_HS, max_j_matrix_C = max_j_matrix_C, running_time = running_time)) } #n_voxels = 500; n_subj = 10; n_tr = 50; n_v_bic = 20; n_bics=5 #pi_1_to_1 = 0.6; pi_0_to_1 = 0.1; pi_1_1 = 0.05; p_s_c = 0.01; p_s = 0.9 #in_bic_noise = 1; out_bic_noise = 1 run_simulation_multiple_bics_normal<-function(n_bics = 2, wrapper = run_complete_musbic_algorithm_normal, n_voxels = 500, n_subj = 2, n_tr = 50, n_v_bic = 20, pi_1_to_1 = 0.6, pi_0_to_1 = 0.1, pi_1_1 = 0.05, p_s_c = 0, p_s = 1, in_bic_noise = 0.1, out_bic_noise = 0.05,setAsGarbage = F,...){ m = matrix(nrow=n_voxels,ncol = n_subj*n_tr);m[,] = 0 rownames(m) = paste ("v",1:n_voxels,sep="_") colnames(m) = rep("",ncol(m)) Cs = list();cols = list() for (k in 1:n_bics){ C = matrix(nrow = n_subj,ncol = n_tr);C[,] = 0 colnames(C) = rep("",ncol(C));rownames(C) = rep("",nrow(C)) count = 1 for(i in 1:n_subj){ for (j in 1:n_tr){ colnames(m)[count] = paste("subj",i,";TR",j,sep="") colnames(C)[j] = paste("TR",j,sep="") rownames(C)[i] = paste("subj",i,sep="") count = count+1 } } # currently assume that w is 1 bic_cols = c() for(i in 1:n_subj){ for (j in 1:n_tr){ curr_pi = pi_1_1 if (j>1 && C[i,j-1]==1){curr_pi = pi_1_to_1} if (j>1 && C[i,j-1]==0){curr_pi = pi_0_to_1} samp = rbinom(1,1,curr_pi) if (samp==1){ C[i,j] = 1 bic_cols = c(bic_cols,paste("subj",i,";TR",j,sep="")) } } } Cs[[k]] = C; cols[[k]] = bic_cols } Hs = list(); HSs = list(); rows = list() for (k in 1:n_bics){ C = Cs[[k]]; bic_cols = cols[[k]] bic_rows = sample(rownames(m))[1:n_v_bic] H = rep(0,n_voxels);names(H) = rownames(m) H[bic_rows] = 1 HS = matrix (nrow = nrow(m),ncol=n_subj);HS[,]=0 rownames(HS) = names(H); colnames(HS) = rownames(C) # set HS if (p_s == 1 && p_s_c == 0){ HS[bic_rows,] = 1 #m[bic_rows,bic_cols] = m[bic_rows,bic_cols]+1 m[bic_rows,bic_cols] = 1 } else{ for (j in 1:ncol(HS)){ for (i in 1:nrow(HS)){ if (H[i] == 1 && rbinom(1,1,p_s)==1){HS[i,j]=1} if (H[i] == 0 && rbinom(1,1,p_s_c)==1){HS[i,j]=1} } curr_subj_trs = colnames(m)[grepl(pattern = paste("subj",j,";",sep=''),colnames(m))] curr_subj_trs = curr_subj_trs[C[j,]==1] #m[HS[,j]==1,curr_subj_trs] = m[HS[,j]==1,curr_subj_trs]+1 m[HS[,j]==1,curr_subj_trs] = 1 } } if (setAsGarbage){ m[,]=0;H[]=0;HS[,] = 0;C[,] = 0 } Hs[[k]] = H; HSs[[k]] = HS; rows[[k]] = bic_rows } # add noise to the data for (i in 1:nrow(m)){ for (j in 1:ncol(m)){ currp = out_bic_noise if (m[i,j]==1){currp = in_bic_noise} samp = rnorm(1,mean = 0,sd = currp) #print (paste(m[i,j],currp,samp)) m[i,j] = m[i,j]+samp } } heatmap(t(m),scale='none',col = gray.colors(3)) # run the algorithm # for tests: # wrapper = wrapper2_algo # solution = wrapper(m,initial_algo='isa',gibbs_reps=0,model='norm') # solution = wrapper(m,initial_algo='bimax',gibbs_reps=0,model='norm') # wrapper = run_complete_musbic_algorithm_normal # solution = wrapper(m,initial_method='bimax',gibbs_reps=20) # solution = wrapper(m,initial_method='bimax',gibbs_reps=0) ptm <- proc.time() # test using the binary algorithm solution = wrapper(m,...) running_time = (proc.time() - ptm)[3] # get some Jaccard scores max_j_matrix_H = matrix(nrow = n_bics,ncol = length(solution[[1]])) max_j_matrix_H[,] = 0 max_j_matrix_HS = matrix(nrow = n_bics,ncol = length(solution[[1]])) max_j_matrix_HS[,] = 0 max_j_matrix_C = matrix(nrow = n_bics,ncol = length(solution[[1]])) max_j_matrix_C[,] = 0 for (i in 1:nrow(max_j_matrix_H)){ if (ncol(max_j_matrix_H)==0){next} for (j in 1:ncol(max_j_matrix_H)){ gibbs_sol = extract_gibbs_solution(m,solution[[4]][[j]]) C_j = binary_matrix_jaccard(Cs[[i]],gibbs_sol$C) HS_j = binary_matrix_jaccard (HSs[[i]],gibbs_sol$HS) H_j = binary_matrix_jaccard (Hs[[i]],gibbs_sol$H) # update max_j_matrix_H[i,j] = max(max_j_matrix_H[i,j],H_j) max_j_matrix_HS[i,j] = max(max_j_matrix_HS[i,j],HS_j) max_j_matrix_C[i,j] = max(max_j_matrix_C[i,j],C_j) } } return (list(max_j_matrix_H = max_j_matrix_H, max_j_matrix_HS = max_j_matrix_HS, max_j_matrix_C = max_j_matrix_C, running_time = running_time)) } run_simulation_multiple_bics_normal_binary_approx<-function(n_bics = 2, wrapper = run_complete_musbic_algorithm, n_voxels = 500, n_subj = 2, n_tr = 50, n_v_bic = 20, pi_1_to_1 = 0.6, pi_0_to_1 = 0.1, pi_1_1 = 0.05, p_s_c = 0, p_s = 1, in_bic_noise = 0.1, out_bic_noise = 0.05,setAsGarbage = F,bin_q_thr = 0.75,...){ m = matrix(nrow=n_voxels,ncol = n_subj*n_tr);m[,] = 0 rownames(m) = paste ("v",1:n_voxels,sep="_") colnames(m) = rep("",ncol(m)) Cs = list();cols = list() for (k in 1:n_bics){ C = matrix(nrow = n_subj,ncol = n_tr);C[,] = 0 colnames(C) = rep("",ncol(C));rownames(C) = rep("",nrow(C)) count = 1 for(i in 1:n_subj){ for (j in 1:n_tr){ colnames(m)[count] = paste("subj",i,";TR",j,sep="") colnames(C)[j] = paste("TR",j,sep="") rownames(C)[i] = paste("subj",i,sep="") count = count+1 } } # currently assume that w is 1 bic_cols = c() for(i in 1:n_subj){ for (j in 1:n_tr){ curr_pi = pi_1_1 if (j>1 && C[i,j-1]==1){curr_pi = pi_1_to_1} if (j>1 && C[i,j-1]==0){curr_pi = pi_0_to_1} samp = rbinom(1,1,curr_pi) if (samp==1){ C[i,j] = 1 bic_cols = c(bic_cols,paste("subj",i,";TR",j,sep="")) } } } Cs[[k]] = C; cols[[k]] = bic_cols } Hs = list(); HSs = list(); rows = list() for (k in 1:n_bics){ C = Cs[[k]]; bic_cols = cols[[k]] bic_rows = sample(rownames(m))[1:n_v_bic] H = rep(0,n_voxels);names(H) = rownames(m) H[bic_rows] = 1 HS = matrix (nrow = nrow(m),ncol=n_subj);HS[,]=0 rownames(HS) = names(H); colnames(HS) = rownames(C) # set HS if (p_s == 1 && p_s_c == 0){ HS[bic_rows,] = 1 #m[bic_rows,bic_cols] = m[bic_rows,bic_cols]+1 m[bic_rows,bic_cols] = 1 } else{ for (j in 1:ncol(HS)){ for (i in 1:nrow(HS)){ if (H[i] == 1 && rbinom(1,1,p_s)==1){HS[i,j]=1} if (H[i] == 0 && rbinom(1,1,p_s_c)==1){HS[i,j]=1} } curr_subj_trs = colnames(m)[grepl(pattern = paste("subj",j,";",sep=''),colnames(m))] curr_subj_trs = curr_subj_trs[C[j,]==1] #m[HS[,j]==1,curr_subj_trs] = m[HS[,j]==1,curr_subj_trs]+1 m[HS[,j]==1,curr_subj_trs] = 1 } } if (setAsGarbage){ m[,]=0;H[]=0;HS[,] = 0;C[,] = 0 } Hs[[k]] = H; HSs[[k]] = HS; rows[[k]] = bic_rows } # add noise to the data for (i in 1:nrow(m)){ for (j in 1:ncol(m)){ currp = out_bic_noise if (m[i,j]==1){currp = in_bic_noise} samp = rnorm(1,mean = 0,sd = currp) #print (paste(m[i,j],currp,samp)) m[i,j] = m[i,j]+samp } } #par(nfrow=c(1,2)) heatmap(t(m),scale='none',col = gray.colors(3)) #hist(m) ptm <- proc.time() # test using the binary algorithm bin_thr = quantile(c(m),bin_q_thr) print (bin_thr) binM = m>=bin_thr mode(binM) = 'numeric' heatmap(t(binM),scale='none',col = gray.colors(3)) solution = wrapper(binM,...) running_time = (proc.time() - ptm)[3] # get some Jaccard scores max_j_matrix_H = matrix(nrow = n_bics,ncol = length(solution[[1]])) max_j_matrix_H[,] = 0 max_j_matrix_HS = matrix(nrow = n_bics,ncol = length(solution[[1]])) max_j_matrix_HS[,] = 0 max_j_matrix_C = matrix(nrow = n_bics,ncol = length(solution[[1]])) max_j_matrix_C[,] = 0 for (i in 1:nrow(max_j_matrix_H)){ if (ncol(max_j_matrix_H)==0){next} for (j in 1:ncol(max_j_matrix_H)){ gibbs_sol = extract_gibbs_solution(m,solution[[4]][[j]]) C_j = binary_matrix_jaccard(Cs[[i]],gibbs_sol$C) HS_j = binary_matrix_jaccard (HSs[[i]],gibbs_sol$HS) H_j = binary_matrix_jaccard (Hs[[i]],gibbs_sol$H) # update max_j_matrix_H[i,j] = max(max_j_matrix_H[i,j],H_j) max_j_matrix_HS[i,j] = max(max_j_matrix_HS[i,j],HS_j) max_j_matrix_C[i,j] = max(max_j_matrix_C[i,j],C_j) } } return (list(max_j_matrix_H = max_j_matrix_H, max_j_matrix_HS = max_j_matrix_HS, max_j_matrix_C = max_j_matrix_C, running_time = running_time)) } binary_matrix_jaccard <-function(x,y){ x = c(x); y = c(y) i = sum(x&y);u = sum(x|y) return (i/u) } binary_FPR <-function(x,y){ x = c(x); y = c(y) fp = sum(x&!y);n = sum(!y) return (fp/n) } binary_TPR <-function(x,y){ x = c(x); y = c(y) tp = sum(x&y);n = sum(y) return (tp/n) } binary_precision <-function(x,y){ x = c(x); y = c(y) tp = sum(x&y);p = sum(x) return (tp/p) } ################ Evaluation of solutions for real data analysis ############### get_bicluster_mean <-function(x,rows,cols){ return (mean(x[rows,cols])) } get_bicluster_mean_cor <-function(x,rows,cols){ return (corr_means(cor(t(x[rows,cols])))) } get_bicluster_mean_cor_outside <-function(x,rows,cols){ col_inds = which(is.element(colnames(x),set=cols)) return (corr_means(cor(t(x[rows,-col_inds])))) } get_bicluster_rows_mean_cor <-function(x,rows){ return (corr_means(cor(t(x[rows,])))) } corr_means <-function(m){ return (mean(m[lower.tri(m)])) } run_max_avg<-function(x){ x = as.matrix(x) if (length(x)==0){return (0)} a1 = apply(x,1,max) a2 = apply(x,2,max) return (mean(c(a1,a2))) }