apply flashr with OCV

wrapper function for OCV

CVPMD_softImpute=function(Y,c_s,K,fold = 10, method = "PMD"){
  N = dim(Y)[1]
  P = dim(Y)[2]
  colindex = matrix(sample(P,P),ncol = fold)
  rowindex = matrix(sample(N,N),ncol = fold)

  missing= array(0,dim = c(fold,N,P))
  foldindex = array(0,dim = c(fold,fold,2))
  for(i in 1:fold){
    for(j in 1:fold){
      foldindex[i,j,1] = i
      foldindex[i,j,2] = (i+j) %% fold
    }
  }
  foldindex[which(foldindex == 0)] = fold
  for(i in 1:fold){
    missing[i, , ] = Y
    for(j in 1:fold){
      missing[i,rowindex[,foldindex[j,i,1]],colindex[,foldindex[j,i,2]]] = NA
    }
    missing[i,,which(colSums(missing[i,,],na.rm = T) ==0)] = Y[,which(colSums(missing[i,,],na.rm = T) ==0)]
  }
  # c_s is the candicate of shrinkage parameter
  n_s = length(c_s)
  # rmse for each grids
  CVRMSE = rep(0,n_s)
  minrmse = Inf
  opt_s = 0
  # for each candidate, we run it N_sim times
  for(t_s in 1:n_s){
    # for each grid
    # each time we set the rmse to zeros
    rmse = rep(0,fold)
    for(i in 1:fold){
      if(method == "PMD"){
        res_log = capture.output({out = PMD(missing[i,,], sumabs = c_s[t_s], sumabsv = NULL, sumabsu = NULL,K = K)})
      }else{
        out = softImpute(missing[i,,], rank.max = K,lambda = c_s[t_s])
      }
      if(length(out$d)==1){
        misshat = (out$d) * out$u %*% t(out$v)
      }else{
        misshat = out$u %*%  diag(out$d) %*% t(out$v)
      }
      for(j in 1:fold){
        # for each fold j
        rmse[i] = rmse[i] + sum((Y[rowindex[,foldindex[j,i,1]],colindex[,foldindex[j,i,2]]] -
                                   misshat[rowindex[,foldindex[j,i,1]],colindex[,foldindex[j,i,2]]])^2,na.rm = TRUE)
      }
    } #get the result for one run
    CVRMSE[t_s] = CVRMSE[t_s] + sqrt(sum(rmse)/(N*P))
    if(CVRMSE[t_s] < minrmse){
      minrmse = CVRMSE[t_s]
      opt_s = c_s[t_s]
    }
  }
  return(list(opt_s = opt_s, output = CVRMSE))
}


PMA.wrapper = function(Y,ngrids = 10,K,fold = 10){
  library(PMA)
  N = dim(Y)[1]
  P = dim(Y)[2]
  c_s = seq(0.3,0.8,len=ngrids)
  cvout = CVPMD_softImpute(Y,c_s,K ,fold , method = "PMD")
  res_log = capture.output({out = PMD(Y,sumabsu = NULL, sumabsv = NULL, sumabs = cvout$opt_s ,K = K)})
  return(list(d = out$d, u = out$u, v = out$v))
}

softImpute.wrapper = function(Y,ngrids = 10,K,fold = 10){
  library(softImpute)
  N = dim(Y)[1]
  P = dim(Y)[2]
  c_s = seq(0,100,len=ngrids)
  cvout = CVPMD_softImpute(Y,c_s,K ,fold , method = "softImpute")
  out = softImpute(Y, rank.max = K,lambda = cvout$opt_s)
  return(list(d = out$d, u = out$u, v = out$v))
}


OCV_index=function(Y,k_fold = 5){
  N = dim(Y)[1]
  P = dim(Y)[2]
  colindex = matrix(sample(P,P),ncol = k_fold)
  rowindex = matrix(sample(N,N),ncol = k_fold)

  foldindex = array(0,dim = c(k_fold,k_fold,2))
  for(i in 1:k_fold){
    for(j in 1:k_fold){
      foldindex[i,j,1] = i
      foldindex[i,j,2] = (i+j) %% k_fold
    }
  }
  foldindex[which(foldindex == 0)] = k_fold

  return(list(foldindex = foldindex, rowindex = rowindex, colindex = colindex))
}
# OCVindex = OCV_index(Y,k_fold = 5)

OCV_data = function(Y,OCVindex,k_fold = 5){
  N = dim(Y)[1]
  P = dim(Y)[2]
  colindex = OCVindex$colindex
  rowindex = OCVindex$rowindex
  foldindex = OCVindex$foldindex
  missing= array(0,dim = c(k_fold,N,P))
  for(i in 1:k_fold){
    missing[i, , ] = Y
    for(j in 1:k_fold){
      missing[i,rowindex[,foldindex[j,i,1]],colindex[,foldindex[j,i,2]]] = NA
    }
    missing[i,,which(colSums(missing[i,,],na.rm = T) ==0)] = Y[,which(colSums(missing[i,,],na.rm = T) ==0)]
  }
  return(missing)
}
# OCVdata = OCV_data(Y,OCVindex,k_fold = 5)


OCV_SSE = function(Y,OCVindex,OCVdata,k_fold = 5,method = "flash",Kmax = 50){
  colindex = OCVindex$colindex
  rowindex = OCVindex$rowindex
  foldindex = OCVindex$foldindex
  missing = OCVdata
  SSE = rep(0,k_fold)
  for(i in 1:k_fold){
    miss_hat = call_method(missing[i,,], method = method, Kmax = Kmax)
    for(j in 1:k_fold){
      SSE[i] = SSE[i] + sum((Y[rowindex[,foldindex[j,i,1]],colindex[,foldindex[j,i,2]]] -
                               miss_hat[rowindex[,foldindex[j,i,1]],colindex[,foldindex[j,i,2]]])^2,na.rm = TRUE)
    }
  }
  RMSE = sqrt(sum(SSE)/(dim(Y)[1] * dim(Y)[2]))
  return(RMSE)
}

call_method = function(Y_data,method,Kmax = 50){
  if(method == "flash"){
    data = flashr::flash_set_data(Y_data)
    f_greedy = flashr::flash_add_greedy(data,K=Kmax)
    Y_hat = f_greedy$EL %*% t(f_greedy$EF)
  }else if(method == "flash_wn"){
    data = flashr::flash_set_data(Y_data)
    f_greedy = flashr::flash_add_greedy(data,K=Kmax,nullcheck=FALSE)
    Y_hat = f_greedy$EL %*% t(f_greedy$EF)
  }else if(method == "flash_gb"){
    data = flashr::flash_set_data(Y_data)
    f_greedy = flashr::flash_add_greedy(data,K=Kmax)
    f = flashr::flash_backfit(data,f_greedy)
    Y_hat = f$EL %*% t(f$EF)
  }else if(method == "flash_bf"){
    data = flashr::flash_set_data(Y_data)
    f_data = flashr::flash_add_factors_from_data(data,K = Kmax)
    f = flashr::flash_backfit(data,f_data)
    Y_hat = f$EL %*% t(f$EF)
  }else if(method == "PMD"){
    res_log = capture.output({out = PMA.wrapper(Y_data,ngrids = 10,K = Kmax,fold = 10)})
    if(length(out$d)==1){
      Y_hat =  (out$d) * out$u %*% t(out$v)
    }else{
      Y_hat =  out$u %*% diag(out$d) %*% t(out$v)
    }
  }else if(method == 'PN'){
    library(ebnm)
    library(flashr)
    data = flashr::flash_set_data(Y_data)
    f_greedy = flashr::flash_add_greedy(data,K=Kmax,ebnm_fn = ebnm_pn)
    Y_hat = f_greedy$EL %*% t(f_greedy$EF)
  }else if(method == 'PMD_d'){
    out = PMA::PMD(Y_data, K = Kmax)
    Y_hat =  out$u %*% diag(out$d) %*% t(out$v)
  }else if(method == 'softImpute_d'){
    gsoft = softImpute::softImpute(Y_data, rank.max = Kmax)
    Y_hat =  gsoft$u %*% diag(gsoft$d) %*% t(gsoft$v)
  }else if(method == "softImpute"){
    gsoft = try(softImpute.wrapper(Y_data,ngrids = 10,K = Kmax,fold = 10))
    if(length(gsoft$d)==1){
      Y_hat =  try((gsoft$d) * gsoft$u %*% t(gsoft$v))
    }else{
      Y_hat =  try(gsoft$u %*% diag(gsoft$d) %*% t(gsoft$v))
    }
  }else if(method == "SF_flash"){
    gsoft = try(softImpute.wrapper(Y_data,ngrids = 10,K = Kmax,fold = 10))
    LL = gsoft$u %*% diag(gsoft$d)
    FF = gsoft$v
    data = flashr::flash_set_data(Y_data)
    f_data = flashr::flash_add_lf(data, LL=LL, FF=FF)
    f = flashr::flash_backfit(data,f_data)
    Y_hat = f$EL %*% t(f$EF)
  }else{
    # stop("the method is not implemented yet, please check it out")
    Y_hat = matrix(0,ncol = dim(Y_data)[2],nrow = dim(Y_data)[1])
  }
  return(Y_hat)
}

NSF_OCV = function(OCVindex){
  if(file.exists("NBSFAout")){
    unlink("NBSFAout", recursive= T)
  }
  system("mkdir NBSFAout")
  # write the missing index

  writeMat("./NBSFAout/OCVindex.mat", OCVindex = OCVindex)

  system("matlab -nosplash -nodesktop -r \"addpath(\'../\'); run(\'missvalue.m\');exit;\" ")
  g_nsfa = readMat("./NBSFAout/NSFAresult.mat")
  nsfa_rmse = as.vector(g_nsfa$RMSE)
  return(nsfa_rmse)
}

This file name is wrapper.R

matlab code

In this folder

ddpath('~/HG/flash/data/missingvalue/methods/NBSF/nsfa-master/');
addpath('~/HG/flash/data/missingvalue/methods/NBSF/nsfa-master/utils/');

% you need to change the path and the data file (not centered)
load ../Ydata.mat;
Ycentered=Y;
settings=defaultsettings();
[settings.D,settings.N]=size(Ycentered);
settings.iterations=100;
% now the missing is decide by the this is the missing index.(or you can upload the missing matrix)

file_path = fullfile(pwd,'NBSFAout/OCVindex.mat');
load(file_path);

% missindex = ones(settings.D,settings.N); this is not missing index this
% is observed index.
missindex = (~isnan(Y));
k_fold = 10;
foldindex = zeros(k_fold,k_fold,2);
for k = 1:2
    for j = 1:k_fold
        for i = 1:k_fold
            foldindex(i,j,k) = OCVindex.foldindex(i + k_fold*(j - 1) + k_fold * k_fold * (k - 1));
        end
    end
end

SSE = zeros(1,k_fold);
for i = 1:k_fold
    mvmask = missindex;
    for j = 1:k_fold
        mvmask(OCVindex.rowindex(:,foldindex(j,i,1)),OCVindex.colindex(:,foldindex(j,i,2))) = 0;
    end
    mvmask(:,sum(mvmask) == 0) = (missindex(:,sum(mvmask) == 0));
    initialsample=init_nsfa(settings);
    [finalsample,resultstable]=nsfa(Ycentered,mvmask,initialsample,settings);
    Y_hat = (finalsample.G *  finalsample.X);
    for j = 1:k_fold
        SSE(i) = SSE(i) + ...
            sum(nansum((Ycentered(OCVindex.rowindex(:,foldindex(j,i,1)),OCVindex.colindex(:,foldindex(j,i,2)))- ...
            Y_hat(OCVindex.rowindex(:,foldindex(j,i,1)),OCVindex.colindex(:,foldindex(j,i,2)))).^2));
    end

end

RMSE = sqrt(sum(SSE)/(settings.D*settings.N));
RMSE;

save_path = fullfile(pwd,'NBSFAout/NSFAresult');
save(save_path,'RMSE');

you need change the path ~/HG/flash/data/missingvalue/methods/NBSF/nsfa-master/ and ~/HG/flash/data/missingvalue/methods/NBSF/nsfa-master/utils/ which have the NBSFA matlab code.

creat the run.R file to run the code

In this case, we have problem in running flashr on RCC. so we just try other methods first.

source("../wrapper.R")
library(R.matlab)
## run the code
Y_centered = readMat("../example.mat")
Y = Y_centered$Y
# in the matlab package of NSF, the use the centered data by rows
N = dim(Y)[1]
P = dim(Y)[2]
Y = Y - rowMeans(Y) %*% t(rep(1,P))
OCVindex = OCV_index(Y,k_fold = 10)
OCVdata = OCV_data(Y,OCVindex,k_fold = 10)

PN_rmse = OCV_SSE(Y,OCVindex,OCVdata,k_fold = 10,method = "PN",Kmax = 40)
flashG_rmse = OCV_SSE(Y,OCVindex,OCVdata,k_fold = 10,method = "flash",Kmax = 40)
flashGwn_rmse = OCV_SSE(Y,OCVindex,OCVdata,k_fold = 10,method = "flash_wn",Kmax = 40)
flashB_rmse = OCV_SSE(Y,OCVindex,OCVdata,k_fold = 10,method = "flash_bf",Kmax = 40)
pmd_rmse = OCV_SSE(Y,OCVindex,OCVdata,k_fold = 10,method = "PMD",Kmax = 40)
soft_rmse = OCV_SSE(Y,OCVindex,OCVdata,k_fold = 10,method = "softImpute",Kmax = 40)
nsfa_rmse = NSF_OCV(OCVindex)

result = c(PN_rmse,flashG_rmse,flashGwn_rmse,flashB_rmse,nsfa_rmse,pmd_rmse,soft_rmse)
saveRDS(result, "./NBSFAout/output.rds")
print(sessionInfo())

this file name is runOCV.R

To read different data sets, please click here for more details

creat folder to track the result

mkdir outlog

creat .sbatch file

#!/bin/bash

#SBATCH --job-name=arrayJob
#SBATCH --output=./outlog/arrayJob_%A_%a.out
#SBATCH --error=./outlog/arrayJob_%A_%a.err
#SBATCH --array=1-20
#SBATCH --time=26:00:00
#SBATCH --partition=mstephens
#SBATCH --ntasks=1
#SBATCH --mem-per-cpu=8000

######################
# Begin work section #
######################

# Print this sub-job's task ID
module load R
module load matlab

mkdir test${SLURM_ARRAY_TASK_ID}
cd test${SLURM_ARRAY_TASK_ID}
Rscript --verbose ../runOCV.R

this file is name as OCV.sbatch