# R/Allfuncs.R In wgmao/DataRemix: DataRemix: a universal data transformation for optimal inference from gene expression datasets

#### Documented in DataRemixDataRemix_display

```#library(MASS)

tr <- function(matrix){
return(sum(diag(matrix)))
}#tr

normF <- function(x){
return(sqrt(sum(x^2)))
}#fastcor

kernel <- function(x,y){
#x <- as.matrix(x)
#y <- as.matrix(y)
#if (ncol(x)==1){
#  x <- t(x)
#}#if
#if (ncol(y)==1){
#  y<-t(y)
#}#if
res <- matrix(0, nrow = nrow(y), ncol = nrow(x))
for (i in 1:nrow(y)){
for (j in 1:nrow(x)){
res[i,j] <- exp(-0.5*normF(y[i,]-x[j,]))
}#for j
}#for i
return(res)
}#kernel

marginal_likelihood <- function(sigma, K, y){
delta_matrix <- diag(sigma^2, nrow = nrow(K))
ml <- -0.5*log(det(K+delta_matrix))-0.5*t(y)%*%solve(K+delta_matrix)%*%y
return(-ml)
}#marginal likelihood

feature_map <- function(x, basis, b){
x_row <- nrow(x) #num of points
x_col <- ncol(x) #dimension
mt <- nrow(basis)
res <- matrix(1, nrow = mt, ncol = x_row) #mt: number of bases

for (j in 1:x_row){
for (i in 1:nrow(res)){
res[i,j] <- sqrt(2/mt)*cos(sum(x[j,]*basis[i,])+b[i])
}#for i
}#for j

return(res)
}#feature_map

bayesian <- function(x, theta, basis, b){
dim_x <- length(x)
x <- matrix(x, nrow = 1)
phi <- feature_map(x, basis, b)

approx <- t(phi)%*%theta
return(-approx)
}#bayesian

bayesian_A <- function(record, y, sigma, basis, b){
phi_rec <- feature_map(record, basis, b)

#posterior infer
A <- 1/sigma^2*phi_rec%*%t(phi_rec)+diag(1, nrow = nrow(basis))
A_inv <- solve(A)
mu <- 1/sigma^2*A_inv%*%phi_rec%*%y
variance <- A_inv
#sample from posterior
theta <- mvrnorm(1, mu = mu, Sigma = variance)
theta <- matrix(theta,ncol = 1)

return(theta)
}#bayesian

thompson_sampling <- function (lower, upper, record, K, y, basis, b){
#infer sigma
set.seed(1)
#sigma_res <- optim(0,marginal_likelihood,method = "Brent", lower=-10, upper = 10, K=K,y=y)
#sigma <- sigma_res\$par
sigma <- 0.1

#infer the bayesian
theta <- bayesian_A(record, y, sigma, basis, b)

par <- c()
minimum <- 0
seed_num <- 20
sp_list1 <- runif(seed_num, min = lower[1], max = upper[1])
sp_list2 <- runif(seed_num, min = lower[2], max = upper[2])
sp_list3 <- runif(seed_num, min = lower[3], max = upper[3])
sp_list <- cbind(sp_list1,sp_list2,sp_list3)
for (i in 1:nrow(sp_list)){
sp <- sp_list[i,]
bayesian_res <- optim(sp, bayesian, method = "L-BFGS-B", lower = lower, upper= upper, theta = theta, basis = basis, b = b)
if (bayesian_res\$value < minimum){
minimum <- bayesian_res\$value
par <- bayesian_res\$par
}#if
}#for i

return(list(par=par))
}#thompson_sampling

main <- function(record, lower_limit, upper_limit, basis, b){
#log scale on mu
record[,length(lower_limit)] <- log10(record[,length(lower_limit)])
lower_limit[length(lower_limit)]<- log10(lower_limit[length(lower_limit)])
upper_limit[length(upper_limit)] <- log10(upper_limit[length(upper_limit)])

#scale
mean_val <- apply(record,2,mean)
sd_val <- apply(record,2,sd)
for (j in 1:length(sd_val)){
if (sd_val[j]==0){
sd_val[j] <-1
mean_val[j] <- 0
}#if
}#for j
record_scale <- sweep(record,2,mean_val,"-")
record_scale <- sweep(record_scale,2,sd_val,"/")

#pre-computation
y <- matrix(record_scale[,ncol(record_scale)],ncol=1)
record <- record_scale[,-ncol(record_scale)]
K <- kernel(record, record)

#scale the limit
lower_scale <- c()
upper_scale <- c()
for (j in 1:length(lower_limit)){
lower_scale <- c(lower_scale, (lower_limit[j]-mean_val[j])/sd_val[j] )
upper_scale <- c(upper_scale, (upper_limit[j]-mean_val[j])/sd_val[j] )
}

res <- thompson_sampling(lower = lower_scale, upper = upper_scale, record, K, y, basis, b)

#evaluate the new x
new_comb <- c()
for (j in 1:length(lower_limit)){
if (j==1){
new_comb <- c(new_comb, min(max(round(res\$par[j]*sd_val[j]+mean_val[j]),lower_limit[j]), upper_limit[j])  )
}else if (j==3){
new_comb <- c(new_comb, 10^min(max(res\$par[j]*sd_val[j]+mean_val[j],lower_limit[j]), upper_limit[j])   )
}else{
new_comb <- c(new_comb, min(max(res\$par[j]*sd_val[j]+mean_val[j],lower_limit[j]), upper_limit[j])   )
}#else
}#for j
return(new_comb)
}#main

random_sample <- function(lower, upper){
para1 <- sample(lower[1]:upper[1],1)
para2 <- runif(1, min = lower[2], max = upper[2])
para3 <- runif(1, min = lower[3], max = upper[3])
return(c(para1,para2,para3))
}#random_sample

coarsed_grid <- function(lower, upper){
loc <- c(0.25, 0.5, 0.75)
para <- c()
for (i in 1:length(loc)){
for (j in 1:length(loc)){
for (k in 1:length(loc)){
para1 <- floor(lower[1]+loc[i]*(upper[1]-lower[1]))
para2 <- lower[2]+loc[j]*(upper[2]-lower[2])
para3 <- 10^(log10(lower[3])+loc[k]*(log10(upper[3])-log10(lower[3])))
para <- rbind(para, c(para1, para2, para3))
}#for k
}#for j
}#for i
return(para)
}#coarsed_grid

SVDcombine<-function(svdData, k, p=1, lambda=0){
if (k >1){
SVDcombine=svdData\$u[, 1:k] %*% diag(svdData\$d[1:k]^(p)) %*% t(svdData\$v[, 1:k])
}else{
SVDcombine= (svdData\$d[1:k]^(p)) *matrix(svdData\$u[, 1:k],ncol = 1) %*% matrix(svdData\$v[, 1:k], nrow = 1)
}

if(lambda>0){
nc=ncol(svdData\$u)
if (k ==nc-1){
SVDcombine=SVDcombine+svdData\$d[(k+1):nc]* lambda* matrix(svdData\$u[, (k+1):nc],ncol=1) %*% matrix(t(svdData\$v[, (k+1):nc]),nrow=1)
}else if (k < nc-1){
SVDcombine=SVDcombine+lambda*svdData\$u[, (k+1):nc] %*% diag(svdData\$d[(k+1):nc]) %*% t(svdData\$v[, (k+1):nc])
}#else if
}#if
rownames(SVDcombine)=rownames(svdData\$u)
colnames(SVDcombine)=rownames(svdData\$v)
SVDcombine
}#SVDcombine

SVDcombineFast <- function(svdData, matrix, k, p=1, lambda=0){
if (k >1){
SVDcombine=svdData\$u[, 1:k] %*% diag(svdData\$d[1:k]^(p)) %*% t(svdData\$v[, 1:k])
}else{
SVDcombine= (svdData\$d[1:k]^(p)) *matrix(svdData\$u[, 1:k],ncol = 1) %*% matrix(svdData\$v[, 1:k], nrow = 1)
}

if(lambda>0){
if (k == 1){
principal <- (svdData\$d[1:k]) *matrix(svdData\$u[, 1:k],ncol = 1) %*% matrix(svdData\$v[, 1:k], nrow = 1)
}else{
principal <- svdData\$u[, 1:k] %*% diag(svdData\$d[1:k]) %*% t(svdData\$v[, 1:k])
}#else
SVDcombine <- SVDcombine+lambda*(matrix-principal)
}#if
rownames(SVDcombine)=rownames(svdData\$u)
colnames(SVDcombine)=rownames(svdData\$v)
SVDcombine
}#SVDcombine

DataRemix <- function(svdres, matrix=NULL, fn, k_limits = c(1, ceiling(length(svdres\$d)/2)), p_limits = c(-1,1), mu_limits = c(1e-12,1), num_of_initialization = 5, num_of_thompson = 600, basis = "omega", basis_size = 2000,xi = 0.1, full = T, verbose = T, ...){

if (basis == "omega_Gaussian"){
basis <- DataRemix::omega_Gaussian
}else if (basis == "omega_Laplacian"){
basis <- DataRemix::omega_Laplacian
}else{
basis <- DataRemix::omega
}#else

mt <- nrow(basis)
set.seed(1)
b <- runif(mt)*2*pi
fast <- !is.null(matrix)

lower_limit <- c(k_limits[1], p_limits[1], mu_limits[1])
upper_limit <- c(k_limits[2], p_limits[2], mu_limits[2])

#Initialization
history <- c()
if (num_of_initialization==0){
para <- coarsed_grid(lower_limit, upper_limit)
for (i in 1:nrow(para)){
if (fast){
rec.term <- fn(SVDcombineFast(svdres, matrix, k= para[i,1], p=para[i,2], lambda = para[i,3]), ...)
}else{
rec.term <- fn(SVDcombine(svdres, k= para[i,1], p=para[i,2], lambda = para[i,3]), ...)
}#else
history <- rbind(history,c(para[i,], rec.term))
}#for i
}else{
for (i in 1:num_of_initialization){
para <- random_sample(lower_limit, upper_limit)
if (fast){
rec.term <- fn(SVDcombineFast(svdres, matrix, k= para[1], p=para[2], lambda = para[3]), ...)
}else{
rec.term <- fn(SVDcombine(svdres, k= para[1], p=para[2], lambda = para[3]), ...)
}#else
history <- rbind(history,c(para, rec.term))
}#for i
}#Random Search

#Thompson
record <- history[,c(1:length(lower_limit),ncol(history))]

if (num_of_thompson > 0){
for (i in 1:num_of_thompson){
uniform <- runif(1)
if (uniform < xi){
para <- random_sample(lower_limit, upper_limit)
}else{
para <- main(record,lower_limit, upper_limit, basis, b)
}#else

#duplicated para
while (length(which(record[,1]==para[1] & record[,2]==para[2] & record[,3] == para[3])) > 0){
para <- random_sample(lower_limit, upper_limit)
}#while

if (fast){
rec.term <- fn(SVDcombineFast(svdres, matrix, k= para[1], p=para[2], lambda = para[3]), ...)
}else{
rec.term <- fn(SVDcombine(svdres, k= para[1], p=para[2], lambda = para[3]), ...)
}#else
history <- rbind(history,c(para, rec.term))
record <- rbind(record, c(para, rec.term[length(rec.term)]))
if (verbose){
print(c(i, para, rec.term[length(rec.term)]))
}#if
}#for i
}#if

if (class(record)=="numeric"){
record <- matrix(record, nrow = 1)
history <- matrix(history, nrow = 1)
}#if

#output record or only the best
if (full){
return(list(para = record, full = history))
}else{
index = which.max(record[,ncol(record)])
return(list(para = matrix(record[index,], nrow=1), full = matrix(history[index,], nrow = 1)))
}#else
}#DataRemix

DataRemix_display <- function(DataRemix.res, col.names = c("Rank", "k", "p", "mu", "mean AUPR", "mean AUC"), top.rank = 15){
if (nrow(DataRemix.res\$full) < top.rank){
if (nrow(DataRemix.res\$full)==1){
knitr::kable(matrix(c(1, DataRemix.res\$full), nrow=1), align = "l", col.names = col.names)
}else{
top.rank <- nrow(DataRemix.res\$full)
knitr::kable(cbind(1:top.rank,DataRemix.res\$full[order(DataRemix.res\$para[,4],
decreasing = T)[1:top.rank],]), align = "l", col.names = col.names)
}#else
}else{
knitr::kable(cbind(1:top.rank,DataRemix.res\$full[order(DataRemix.res\$para[,4],
decreasing = T)[1:top.rank],]), align = "l", col.names = col.names)
}#else
}#DataRemix.display
```
wgmao/DataRemix documentation built on Aug. 6, 2020, 4:49 p.m.