R/denoise_1D.R

In shimlab/HMTree: Hidden Markov tree prior for analysis of multiple functional data

#' Simulate multiple curve data from a given model (either Normal or Poisson). It checks if parameters are valid and calls \code{sim.data} to simulate data. 
#'
#' @param data a matrix of type double, number of colum represent the length of data, number of rows represent number of groups.!!!!!!!!!(TODO)!!!!!!!!!!!
#' @param model.mode either Normal or Poisson; this function simulates \code{num.samples} number of curve data under a model provided as \code{model.mode}.
#' @return \code{simu.curves} returns a list of the following elements
#' \item{true.mean.curve}{a vector of length \code{curve.length}-vector; it contains true mean curve generated by the function \code{simulate_true_curve}}
#' \item{data}{a matrix of \code{num.samples} by \code{curve.length}; it contains \code{num.samples} curve data; each row contains simulated data under a given model \code{model.mode} with mean curve \code{true.mean.curve}}
#'
#' @examples
#' # set up parameters.
#' curve.length = 1024
#' model.mode = 'Poisson'
#' num.samples=10
#' set.seed(666)
#' res.pois = simu.curves(curve.length=curve.length, model.mode=model.mode, num.samples=num.samples)
#' simu.data<-res.pois$data
#' data<-HMTree.parameter<-function(simu.data,model.mode)
#' @export
#' 
HMTree.parameter<-function(data,model.mode){
  #check input
  #model.mode
  if(model.mode!="Poisson" && model.mode!="Normal"){
    stop("model.mode should be either 'Poisson' or 'Normal' ")
  }
  #data
  if(length(data)==0 && typeof(data)!="double"){
    stop("Please input data in correct format(details in description)")
  }
  #should I check data length for each group?????????????????
  x<-list(data=data,model.mode=model.mode)
}


#'Return alpha_vec,open to user for now
#' @export
#' 
HMTree.denoise_1D<-function(x,num.states = 2,tolerance = 0.0001){
  n<-ncol(x$data)
  model.mode<-x$model.mode
  if(model.mode=="Normal"){
    #just save the data(with noise)
    alpha_vec<-x$data
  }
  else if(model.mode=="Poisson"){
    MLE.dat<-MLE.se.logitp(x=x$data, g=NULL, reflect=FALSE, minobs=1, pseudocounts=0.5, all=FALSE, center=FALSE, repara=TRUE, forcebin=FALSE, lm.approx=TRUE, disp="add")
    # compute posterior mean and varince of log(overall intensity). The 1024th coefficient 
    overall.res<-estimate.logp.overallintensity(x = x$data, g = NULL, read.depth = NULL)
    alpha_vec <- t(cbind(MLE.dat$alpha, MLE.dat$alpha.se)) # Compute all 2^n-1 coefficients except root (mean)
    alpha_vec <- cbind(t(cbind(overall.res$lp.mean, overall.res$lp.var)), alpha_vec) # add last coefficient (mean)
    
    #    TEMP Correction: set all NA to zero mean and unit st.d.: 
    for(i in 1:n){
      if(is.na(alpha_vec[1,i])){ # set mean to zero if we have NA
        alpha_vec[1,i] <- 0.0; 
      }
      if(is.na(alpha_vec[2,i])){
        alpha_vec[2,i] <- 1.0;  # set st.d. to one if we have NA
      }
    }
  }
  r<-denoise_1D(alpha_vec,tolerance,n,model.mode,num.states)
  #Give C object names
  return(list(alpha_vec=r[1],wavelet_coef=r[2],denoised_wavelet_coef=r[3],denoised=r[4],param_signa_sqr=r[5],param_marginal_probs=r[6],param_transition_probs=r[7],post_probs=r[8],post_trans_probs=r[9]))
}