R/susiF_workhorse.R

In stephenslab/susiF.alpha: Sum of Single Functions

#' @title Sum of Single Function
#'
#' @description Implementation of the SuSiF method
#'
#' @details Implementation of the SuSiF method
#'
#'
#' @param obj an object of class susiF
#' @param W a list in which element D contains matrix of wavelet d coefficients and
#'  element C contains the vector of scaling coefficients
#' @param X matrix of size n by p contains the covariates
#'

#' @param verbose If \code{verbose = TRUE}, the algorithm's progress,
#' and a summary of the optimization settings are printed to the
#' console.
#' @param lowc_wc  vector of index of which wavelet coefficients display pathologic distribution
#'
#' @param tol a small, non-negative number specifying the convergence
#' tolerance for the IBSS fitting procedure. The fitting procedure
#' will halt when the difference in the variational lower bound, or
#' \dQuote{ELBO} (the objective function to be maximized), is less
#' than \code{tol}.
#'
#' @param maxit Maximum number of IBSS iterations.
#'
#' @param cov_lev numeric between 0 and 1, corresponding to the
#' expected level of coverage of the cs if not specified set to 0.95
#'
#' @param min_purity minimum purity for estimated credible sets
#' @param init_pi0_w starting value of weight on null compoenent in mixsqp
#'  (between 0 and 1)
#' @param control_mixsqp list of parameter for mixsqp function see  mixsqp package
#' @param  cal_obj logical if set as TRUE compute ELBO for convergence monitoring
#' @param nullweight numeric value for penalizing likelihood at point mass 0 (should be between 0 and 1)
#' (usefull in small sample size)
#' @param lowc_wc list of wavelet coefficients that exhibit too little variance
#' @param indx_lst list generated by gen_wavelet_indx for the given level of resolution
#'
#' @param tt output of the cal_Bhat_Shat function
#'
#' @param max_step_EM max_step_EM
#' @param parallel if true use parallel computation
#' @param max_SNP_EM check susiF  description
#' @param cor_small check susiF  description
#' @param is.pois check susiF  description
#' @param e threshold value to avoid computing posterior that have low alpha value. Set it to 0 to compute the entire posterio. default value is 0.001
#' @export
susiF.workhorse <- function(obj,
                            W,
                            X,
                            tol,
                            init_pi0_w ,
                            control_mixsqp ,
                            indx_lst,
                            lowc_wc,
                            nullweight,
                            cal_obj,
                            verbose,
                            cov_lev,
                            min_purity,
                            maxit,
                            tt,
                            parallel=FALSE,
                            max_SNP_EM=500,
                            max_step_EM=1,
                            cor_small=FALSE,
                            is.pois=FALSE,
                            e = 0.001){
#browser()
  G_prior  <- get_G_prior(obj )
  Y_f      <-  cbind( W$D,W$C)
  update_Y <- Y_f
  # numerical value to check breaking condition of while
  # numerical value to check breaking condition of while
  check <- 3*tol
  init  <- TRUE
  v1    <-  rep(1, dim(X)[1])


  if (cor_small){
    df = nrow(X)-1
  }else{
    df =NULL
  }

  if( obj$L_max==1)
  {
    tt   <-cal_Bhat_Shat(Y       = update_Y,
                         X       = X,
                         v1      = v1 ,
                         lowc_wc = lowc_wc )
    Bhat <- tt$Bhat
    Shat <- tt$Shat #UPDATE. could be nicer
    if(is.pois){
      Shat <- update_Shat_pois(Shat     = Shat,
                               indx_lst = indx_lst,
                               lowc_wc  = lowc_wc
                               )
    }
    tpi  <- get_pi(obj,1)
    G_prior <- update_prior(G_prior = G_prior,
                            tpi     = tpi ) #allow EM to start close to previous solution (to double check)

    EM_out  <- EM_pi(G_prior        = G_prior,
                     Bhat           = Bhat,
                     Shat           = Shat,
                     indx_lst       = indx_lst,
                     init_pi0_w     = init_pi0_w,
                     control_mixsqp = control_mixsqp,
                     lowc_wc        = lowc_wc,
                     nullweight     = nullweight,
                     max_SNP_EM     = max_SNP_EM,
                     max_step       = max_step_EM,
                     df             = df

    )

    obj <-  update_susiF_obj(obj = obj ,
                                   l         = 1,
                                   EM_pi     = EM_out,
                                   Bhat      = Bhat,
                                   Shat      = Shat,
                                   indx_lst  = indx_lst,
                                   lowc_wc   = lowc_wc,
                                   cov_lev   =  cov_lev,
                                   e         = e
    )
    obj <- update_ELBO(obj  = obj,
                             get_objective( obj = obj,
                                            Y         = Y_f,
                                            X         = X,
                                            D         = W$D,
                                            C         = W$C,
                                            indx_lst  = indx_lst
                             )
    )

  }else{
    ##### Start While -----
    iter <- 1


    while( (check >tol & iter <maxit))
    {

      for( l in 1:obj$L)
      {

        #print(obj$alpha[[l]])
        update_Y  <-  cal_partial_resid(
          obj = obj,
          l         = (l-1)  ,
          X         = X,
          D         = W$D,
          C         = W$C,
          indx_lst  = indx_lst
        )

        if(verbose){
          print(paste("Fitting effect ", l,", iter" ,  iter ))
        }

        if(init){#recycle operation used to fit the prior


          EM_out <- gen_EM_out (tpi_k= get_pi_G_prior(G_prior),
                                lBF  = log_BF  (G_prior  = G_prior,
                                                Bhat     = tt$Bhat,
                                                Shat     = tt$Shat,
                                                lowc_wc  = lowc_wc,
                                                indx_lst = indx_lst,
                                                df       = df
                                )
          )

          init <- FALSE
        }else{


          tt   <- cal_Bhat_Shat(Y = update_Y,
                                X = X,
                                v1 =v1 ,
                                lowc_wc =lowc_wc )

          tpi <-  get_pi(obj,l)
          G_prior <- update_prior(G_prior, tpi= tpi ) #allow EM to start close to previous solution (to double check)
          if(is.pois){
            tt$Shat <- update_Shat_pois(Shat     = tt$Shat,
                                        indx_lst = indx_lst,
                                        lowc_wc  = lowc_wc
            )
         #   print( tt$Shat)
          }

          EM_out  <- EM_pi(G_prior        = G_prior,
                           Bhat           = tt$Bhat,
                           Shat           = tt$Shat,
                           indx_lst       = indx_lst,
                           init_pi0_w     = init_pi0_w,
                           control_mixsqp = control_mixsqp,
                           lowc_wc        = lowc_wc,
                           nullweight     = nullweight,
                           max_SNP_EM     = max_SNP_EM,
                           max_step       = max_step_EM,
                           df             = df,
                           tol_null_prior = obj$tol_null_prior
          )
           #plot(EM_out$lBF/(sum( EM_out$lBFlBF)))
        }

        #print(h)
        # print(EM_out$lBF[1:10])
        obj <-  update_susiF_obj(obj   = obj ,
                                       l           = l,
                                       EM_pi       = EM_out,
                                       Bhat        = tt$Bhat,
                                       Shat        = tt$Shat,
                                       indx_lst    = indx_lst,
                                       lowc_wc     = lowc_wc,
                                       df          = df
        )


      }#end for l in 1:L  -----

      # plot(obj$lBF[[1]])
      #plot(obj$lBF[[2]])

      #plot(obj$lBF[[3]])
      # save(obj, file ="D:/Document/Serieux/Travail/Package/susiF.alpha/pb_object.RData")
      # break
      ####Check greedy/backfit and stopping condition -----
      obj <- greedy_backfit (obj,
                                   verbose    = verbose,
                                   cov_lev    = cov_lev,
                                   X          = X,
                                   min_purity = min_purity
      )
      sigma2    <- estimate_residual_variance(obj,
                                              Y         = Y_f,
                                              X         = X)
      #print(sigma2)
      obj <- update_residual_variance(obj     = obj,
                                            sigma2    = sigma2 )
      obj <- test_stop_cond(obj      = obj,
                                  check     = check,
                                  cal_obj   = cal_obj,
                                  Y         = Y_f,
                                  X         = X,
                                  D         = W$D,
                                  C         = W$C,
                                  indx_lst  = indx_lst)
      #  print(obj$alpha)
      #print(obj$ELBO)
      check <- obj$check


      iter <- iter +1


    }#end while
  }#end else in if(L==1)

 return(obj)
}