X <- "dataSim1fabia" p=10 alpha=0.03 cyc=50 spl=0 non_negative=1 write_file=0 norm=0 lap=100.0 samples=2000 iter=1 quant=0.01 lowerB=1.1 upperB=0.05*1000 eps=1e-5 spz = 0.5 random = 1 scale = 0 nL = 0 lL = 0 bL = 0 initL = 0, dorescale = FALSE, doini = FALSE eps1 = 1e-10 eps <- as.double(eps) eps1 <- as.double(eps1) init_lapla <- as.double(1) init_psi <- as.double(0.1) samples <- as.integer(sort.int(as.integer(unique(samples)))) initL <- as.integer(initL) iter <- as.integer(iter) norm <- as.integer(norm) cyc <- as.integer(cyc) nL <- as.integer(nL) lL <- as.integer(lL) bL <- as.integer(bL) quant <- as.double(quant) lowerB <- as.double(lowerB) upperB <- as.double(upperB) alpha <- as.double(alpha) non_negative <- as.integer(non_negative) write_file <- as.integer(write_file) p <- as.integer(p) spz <- as.double(spz) scale <- as.double(scale) lap <- as.double(lap) res <- .Call("spfabic", X, p, alpha, cyc, spl, spz, non_negative, random, write_file, init_psi, init_lapla, norm, scale, lap, nL, lL, bL, eps, eps1, samples, initL, iter, quant, lowerB, upperB, PACKAGE = "fabia") spfabia <- function (X, p = 5, alpha = 0.1, cyc = 500, spl = 0, spz = 0.5, non_negative = 0, random = 1, write_file = 1, norm = 1, scale = 0, lap = 1, nL = 0, lL = 0, bL = 0, samples = 0, initL = 0, iter = 1, quant = 0.001, lowerB = 0, upperB = 1000, dorescale = FALSE, doini = FALSE, eps = 0.001, eps1 = 1e-10) { if (missing(X)) { stop("Data file name missing. Stopped.") } message("Running sparse FABIA on a sparse matrix in file >", X, "<:") message(" Number of biclusters ---------------- p: ", p) message(" Sparseness factor --------------- alpha: ", alpha) message(" Number of iterations -------------- cyc: ", cyc) message(" Loading prior parameter ----------- spl: ", spl) message(" Factor prior parameter ------------ spz: ", spz) if (random > 0) { message(" Initialization loadings--------- random: ", random, " = positive") } else { message(" Initialization loadings--------- random: ", random, " = all values") } if (non_negative > 0) { message(" Nonnegative Loadings and Factors ------: ", non_negative, " = Yes") } else { message(" Nonnegative Loadings and Factors ------: ", non_negative, " = No") } if (write_file > 0) { message(" Write results files -------------------: ", write_file, " = Yes") } else { message(" Write results files -------------------: ", write_file, " = No") } if (norm > 0) { message(" Scaling to variance one: --------- norm: ", norm, " = Yes") } else { message(" Scaling -------------------------- norm: ", norm, " = No") } if (scale > 0) { message(" Scaling loadings per iteration -- scale: ", scale, " = to ", scale) } else { message(" Scaling loadings per iteration -- scale: ", scale, " = No") } message(" Constraint variational parameter -- lap: ", lap) if ((nL > 0) && (nL < p)) { message(" Max. number of biclusters per row -- nL: ", nL) if (bL > 0) { message(" starting at ------------------ bL: ", bL) } else { message(" starting at ------------------ bL: ", bL, " = from start") } } else { message(" Max. number of biclusters per row -- nL: ", nL, " = no limit") } if (lL > 0) { message(" Max. number of row elements / biclu. lL: ", lL) if (bL > 0) { message(" starting at ------------------ bL: ", bL) } else { message(" starting at ------------------ bL: ", bL, " = from start") } } else { message(" Max. number of row elements / biclu. lL: ", lL, " = no limit") } if (samples[1] == 0) { message(" Number of samples ------------- samples: ", samples, " = all samples") } else { message(" Number of samples ------------ samples: ", length(samples)) } if (initL[1] == 0) { message(" Random initialization of L ------ initL: ", initL) } else { message(" Biclusters init. by samples ----- initL: ", length(initL)) } message(" Iterations ------------------------iter: ", iter) if (iter > 1) { message(" Quantile of largest L removed --- quant: ", quant) } message(" Lower column sum bound --------- lowerB: ", lowerB) message(" Upper column sum bound --------- upperB: ", upperB) if (dorescale) { message(" Z and L are rescaled -------- dorescale: TRUE") } else { message(" Z and L are not rescaled ---- dorescale: FALSE") } if (doini && dorescale) { message(" Biclusters sorted (information) - doini: TRUE") } else { message(" Biclusters not sorted (infor.) -- doini: FALSE") } eps <- as.double(eps) eps1 <- as.double(eps1) init_lapla <- as.double(1) init_psi <- as.double(0.1) samples <- as.integer(sort.int(as.integer(unique(samples)))) initL <- as.integer(initL) iter <- as.integer(iter) norm <- as.integer(norm) cyc <- as.integer(cyc) nL <- as.integer(nL) lL <- as.integer(lL) bL <- as.integer(bL) quant <- as.double(quant) lowerB <- as.double(lowerB) upperB <- as.double(upperB) alpha <- as.double(alpha) non_negative <- as.integer(non_negative) write_file <- as.integer(write_file) p <- as.integer(p) spz <- as.double(spz) scale <- as.double(scale) lap <- as.double(lap) res <- .Call("spfabic", X, p, alpha, cyc, spl, spz, non_negative, random, write_file, init_psi, init_lapla, norm, scale, lap, nL, lL, bL, eps, eps1, samples, initL, iter, quant, lowerB, upperB, PACKAGE = "fabia") if (is.null(res)) { return(new("Factorization", parameters = list(), n = 1, p1 = 1, p2 = 1, l = 1, center = as.vector(1), scaleData = as.vector(1), X = as.matrix(1), L = noL, Z = nZ, M = as.matrix(1), LZ = as.matrix(1), U = as.matrix(1), avini = as.vector(1), xavini = as.vector(1), ini = as.matrix(1), Psi = as.vector(1), lapla = as.matrix(1))) } l = ncol(res$E_SX_n) n = nrow(res$L) if (dorescale) { pi = p * iter eps <- as.double(0.001) nvect <- as.vector(rep(1, n)) epsn <- eps * nvect iin <- 1/l vz <- iin * apply(res$E_SX_n, 1, function(x) sum(x^2)) vz <- sqrt(vz + 1e-10) ivz <- 1/vz if (length(ivz) == 1) { nZ <- ivz * res$E_SX_n noL <- vz * res$L } else { nZ <- ivz * res$E_SX_n noL <- t(vz * t(res$L)) } } else { nZ <- res$E_SX_n noL <- res$L } if (dorescale && doini) { ini <- matrix(0, l, (pi + 1)) avini <- as.vector(rep(0, (pi + 1))) xavini <- as.vector(rep(0, (l + 1))) idp <- diag(pi) ppL <- matrix(0, pi, pi) for (ite in 0:(iter - 1)) { ppL[((ite * p) + 1):((ite + 1) * p), ((ite * p) + 1):((ite + 1) * p)] <- crossprod(noL[, ((ite * p) + 1):((ite + 1) * p)], (1/(res$Psi[((ite * n) + 1):((ite + 1) * n)] + epsn)) * noL[, ((ite * p) + 1):((ite + 1) * p)]) } for (j in 1:l) { mat <- idp + ppL/res$lapla[j, ] ini[j, 1:pi] <- log(diag(mat)) s <- log(det(mat)) ini[j, pi + 1] <- s xavini[j] <- s } for (i in 1:pi) { avini[i] <- sum(ini[, i]) } ss <- sum(ini[, pi + 1]) xavini[l + 1] <- ss avini[pi + 1] <- ss if ((avini[pi + 1] > 1e-08) && (pi > 1)) { soo <- sort(avini[1:pi], decreasing = TRUE, index.return = TRUE) avini[1:pi] <- avini[soo$ix] noL <- noL[, soo$ix] nZ <- nZ[soo$ix, ] } } else { ini <- as.matrix(1) avini <- as.vector(1) xavini <- as.vector(1) } return(new("Factorization", parameters = list("spfabia", cyc, alpha, spl, spz, p, NULL, NULL, random, scale, norm, NULL, lap, nL, lL, bL, non_negative, write_file, init_lapla, init_psi, samples, initL, iter, quant, lowerB, upperB), n = n, p1 = pi, p2 = pi, l = l, center = as.vector(1), scaleData = as.vector(1), X = as.matrix(1), L = noL, Z = nZ, M = as.matrix(1), LZ = as.matrix(1), U = as.matrix(1), avini = avini, xavini = xavini, ini = ini, Psi = res$Psi, lapla = res$lapla)) }