finish adding flag to check removal of object

QUILT/R/RcppExports.R

@@ -182,8 +182,8 @@ rcpp_fly_weighter <- function(i, gCurrent, gAverage, ll_current, log_mult, ll_re
 }
 
 #' @export
-rcpp_forwardBackwardGibbsNIPT <- function(sampleReads, priorCurrent_m, alphaMatCurrent_tc, eHapsCurrent_tc, transMatRate_tc_H, ff, blocks_for_output, alphaHat_t1, betaHat_t1, alphaHat_t2, betaHat_t2, alphaHat_t3, betaHat_t3, hapSum_tc, wif0, grid_has_read, L_grid, smooth_cm, Jmax_local = 100L, maxDifferenceBetweenReads = 1000, maxEmissionMatrixDifference = 10000000000, run_fb_subset = FALSE, run_fb_grid_offset = 0L, return_genProbs = TRUE, return_hapProbs = TRUE, return_gamma = TRUE, return_alpha = FALSE, return_p_store = FALSE, return_extra = FALSE, grid = 0L, snp_start_1_based = -1L, snp_end_1_based = -1L, generate_fb_snp_offsets = FALSE, suppressOutput = 1L, n_gibbs_starts = 1L, n_gibbs_sample_its = 1L, n_gibbs_burn_in_its = 1L, use_starting_read_labels = FALSE, verbose = FALSE, double_list_of_starting_read_labels = NULL, seed_vector = -1L, prev_list_of_alphaBetaBlocks = NULL, i_snp_block_for_alpha_beta = 1L, haploid_gibbs_equal_weighting = TRUE, gibbs_initialize_iteratively = FALSE, gibbs_initialize_at_first_read = TRUE, do_block_resampling = FALSE, artificial_relabel = -1L, pass_in_alphaBeta = FALSE, update_in_place = FALSE, update_hapSum = FALSE, record_read_set = FALSE, class_sum_cutoff = 0.06, perform_block_gibbs = FALSE, shuffle_bin_radius = 5000L, block_gibbs_iterations = as.integer( c(0)), return_gibbs_block_output = FALSE, return_advanced_gibbs_block_output = FALSE, rescale_eMatRead_t = TRUE, use_smooth_cm_in_block_gibbs = FALSE, block_gibbs_quantile_prob = 0.9, do_shard_ff0_block_gibbs = TRUE) {
-    .Call('_QUILT_rcpp_forwardBackwardGibbsNIPT', PACKAGE = 'QUILT', sampleReads, priorCurrent_m, alphaMatCurrent_tc, eHapsCurrent_tc, transMatRate_tc_H, ff, blocks_for_output, alphaHat_t1, betaHat_t1, alphaHat_t2, betaHat_t2, alphaHat_t3, betaHat_t3, hapSum_tc, wif0, grid_has_read, L_grid, smooth_cm, Jmax_local, maxDifferenceBetweenReads, maxEmissionMatrixDifference, run_fb_subset, run_fb_grid_offset, return_genProbs, return_hapProbs, return_gamma, return_alpha, return_p_store, return_extra, grid, snp_start_1_based, snp_end_1_based, generate_fb_snp_offsets, suppressOutput, n_gibbs_starts, n_gibbs_sample_its, n_gibbs_burn_in_its, use_starting_read_labels, verbose, double_list_of_starting_read_labels, seed_vector, prev_list_of_alphaBetaBlocks, i_snp_block_for_alpha_beta, haploid_gibbs_equal_weighting, gibbs_initialize_iteratively, gibbs_initialize_at_first_read, do_block_resampling, artificial_relabel, pass_in_alphaBeta, update_in_place, update_hapSum, record_read_set, class_sum_cutoff, perform_block_gibbs, shuffle_bin_radius, block_gibbs_iterations, return_gibbs_block_output, return_advanced_gibbs_block_output, rescale_eMatRead_t, use_smooth_cm_in_block_gibbs, block_gibbs_quantile_prob, do_shard_ff0_block_gibbs)
+rcpp_forwardBackwardGibbsNIPT <- function(sampleReads, priorCurrent_m, alphaMatCurrent_tc, eHapsCurrent_tc, transMatRate_tc_H, ff, blocks_for_output, alphaHat_t1, betaHat_t1, alphaHat_t2, betaHat_t2, alphaHat_t3, betaHat_t3, hapSum_tc, hapMatcher, distinctHapsIE, rhb_t, ref_error, which_haps_to_use, wif0, grid_has_read, L_grid, smooth_cm, param_list, Jmax_local = 100L, maxDifferenceBetweenReads = 1000, maxEmissionMatrixDifference = 10000000000, run_fb_subset = FALSE, run_fb_grid_offset = 0L, grid = 0L, snp_start_1_based = -1L, snp_end_1_based = -1L, generate_fb_snp_offsets = FALSE, suppressOutput = 1L, n_gibbs_starts = 1L, n_gibbs_sample_its = 1L, n_gibbs_burn_in_its = 1L, use_starting_read_labels = FALSE, verbose = FALSE, double_list_of_starting_read_labels = NULL, seed_vector = -1L, prev_list_of_alphaBetaBlocks = NULL, i_snp_block_for_alpha_beta = 1L, haploid_gibbs_equal_weighting = TRUE, gibbs_initialize_iteratively = FALSE, gibbs_initialize_at_first_read = TRUE, do_block_resampling = FALSE, artificial_relabel = -1L, pass_in_alphaBeta = FALSE, update_in_place = FALSE, update_hapSum = FALSE, record_read_set = FALSE, class_sum_cutoff = 0.06, perform_block_gibbs = FALSE, shuffle_bin_radius = 5000L, block_gibbs_iterations = as.integer( c(0)), rescale_eMatRead_t = TRUE, use_smooth_cm_in_block_gibbs = FALSE, block_gibbs_quantile_prob = 0.9, do_shard_ff0_block_gibbs = TRUE, use_small_eHapsCurrent_tc = TRUE) {
+    .Call('_QUILT_rcpp_forwardBackwardGibbsNIPT', PACKAGE = 'QUILT', sampleReads, priorCurrent_m, alphaMatCurrent_tc, eHapsCurrent_tc, transMatRate_tc_H, ff, blocks_for_output, alphaHat_t1, betaHat_t1, alphaHat_t2, betaHat_t2, alphaHat_t3, betaHat_t3, hapSum_tc, hapMatcher, distinctHapsIE, rhb_t, ref_error, which_haps_to_use, wif0, grid_has_read, L_grid, smooth_cm, param_list, Jmax_local, maxDifferenceBetweenReads, maxEmissionMatrixDifference, run_fb_subset, run_fb_grid_offset, grid, snp_start_1_based, snp_end_1_based, generate_fb_snp_offsets, suppressOutput, n_gibbs_starts, n_gibbs_sample_its, n_gibbs_burn_in_its, use_starting_read_labels, verbose, double_list_of_starting_read_labels, seed_vector, prev_list_of_alphaBetaBlocks, i_snp_block_for_alpha_beta, haploid_gibbs_equal_weighting, gibbs_initialize_iteratively, gibbs_initialize_at_first_read, do_block_resampling, artificial_relabel, pass_in_alphaBeta, update_in_place, update_hapSum, record_read_set, class_sum_cutoff, perform_block_gibbs, shuffle_bin_radius, block_gibbs_iterations, rescale_eMatRead_t, use_smooth_cm_in_block_gibbs, block_gibbs_quantile_prob, do_shard_ff0_block_gibbs, use_small_eHapsCurrent_tc)
 }
 
 #' @export

QUILT/R/copied_from_stitch.R

@@ -9,7 +9,7 @@ print_message <- function(x, include_mem = FALSE) {
     } else {
         mem <- ""
     }
-    print(
+    message(
         paste0(
             "[", format(Sys.time(), "%Y-%m-%d %H:%M:%S"), "] ", mem, x
         )

QUILT/R/functions.R

@@ -528,7 +528,8 @@ get_and_impute_one_sample <- function(
     print_extra_timing_information,
     n_gibbs_burn_in_its,
     block_gibbs_iterations,
-    plot_per_sample_likelihoods
+    plot_per_sample_likelihoods,
+    use_small_eHapsCurrent_tc
 ) {
 
 
@@ -803,6 +804,7 @@ get_and_impute_one_sample <- function(
                 distinctHapsIE = distinctHapsIE,
                 hapMatcher = hapMatcher,
                 rhb_t = rhb_t,
+                ref_error = ref_error,
                 nSNPs = nSNPs,
                 sampleReads = sampleReads,
                 small_eHapsCurrent_tc = small_eHapsCurrent_tc,
@@ -824,7 +826,6 @@ get_and_impute_one_sample <- function(
                 block_gibbs_iterations = block_gibbs_iterations,
                 perform_block_gibbs = TRUE,
                 make_plots = make_plots,
-                ref_error = ref_error,
                 wif0 = wif0,
                 grid_has_read = grid_has_read,                              
                 plot_description = paste0("it", i_it, ".gibbs"),
@@ -858,7 +859,8 @@ get_and_impute_one_sample <- function(
                 sample_name = sample_name,
                 regionStart = regionStart,
                 regionEnd = regionEnd,
-                buffer = buffer
+                buffer = buffer,
+                use_small_eHapsCurrent_tc = use_small_eHapsCurrent_tc
             )
 
             if (hla_run) {            
@@ -1849,6 +1851,7 @@ impute_one_sample <- function(
     distinctHapsIE,
     hapMatcher,
     rhb_t,
+    ref_error,    
     nSNPs,
     sampleReads,
     small_eHapsCurrent_tc,
@@ -1871,7 +1874,6 @@ impute_one_sample <- function(
     perform_block_gibbs,
     make_plots,
     maxDifferenceBetweenReads,
-    ref_error,
     wif0,
     grid_has_read,
     verbose,
@@ -1908,7 +1910,8 @@ impute_one_sample <- function(
     suppressOutput = 1,
     use_smooth_cm_in_block_gibbs = TRUE,
     block_gibbs_quantile_prob = 0.95,
-    make_plots_block_gibbs = FALSE
+    make_plots_block_gibbs = FALSE,
+    use_small_eHapsCurrent_tc = TRUE
 ) {
     ##
     K <- length(which_haps_to_use)
@@ -1950,40 +1953,60 @@ impute_one_sample <- function(
 
     ## argh
     ## print(paste0("start = ", Sys.time()))
-    inflate_fhb_t_in_place(
-        rhb_t,
-        small_eHapsCurrent_tc,
-        haps_to_get = which_haps_to_use - 1,
-        nSNPs = nSNPs,
-        ref_error = ref_error
-    )
+    if (use_small_eHapsCurrent_tc) {
+        inflate_fhb_t_in_place(
+            rhb_t,
+            small_eHapsCurrent_tc,
+            haps_to_get = which_haps_to_use - 1,
+            nSNPs = nSNPs,
+            ref_error = ref_error
+        )
+    }
     if (make_plots_block_gibbs) {
         return_gibbs_block_output <- TRUE
         return_advanced_gibbs_block_output <- TRUE
     }
     ##
     ## ugh, alphaMatCurrent_tc is a CONSTANT
     ##
-
+    ## param_list <- list(
+    ##     return_alpha = FALSE,
+    ##     return_p_store = FALSE,
+    ##     return_extra = FALSE,        
+    ##     return_genProbs = TRUE,
+    ##     return_hapProbs = TRUE,
+    ##     return_gamma = TRUE,
+    ##     return_gibbs_block_output = FALSE,
+    ##     return_advanced_gibbs_block_output = FALSE
+    ## )
+    param_list <- list(
+        return_alpha = FALSE,
+        return_extra = FALSE,
+        return_genProbs = return_genProbs,
+        return_gamma = as.logical(return_gamma | make_plots),
+        return_hapProbs = return_hapProbs,
+        return_p_store = return_p_store,
+        return_gibbs_block_output = return_gibbs_block_output,
+        return_advanced_gibbs_block_output = return_advanced_gibbs_block_output
+    )
     out <- rcpp_forwardBackwardGibbsNIPT(
         sampleReads = sampleReads,
         priorCurrent_m = small_priorCurrent_m,        
         alphaMatCurrent_tc = small_alphaMatCurrent_tc,
         eHapsCurrent_tc = small_eHapsCurrent_tc,
         transMatRate_tc_H = small_transMatRate_tc_H,
+        hapMatcher = hapMatcher,
+        distinctHapsIE = distinctHapsIE,
+        rhb_t = rhb_t,
+        ref_error = ref_error,
+        which_haps_to_use = which_haps_to_use,
         ff = 0,
         maxDifferenceBetweenReads = maxDifferenceBetweenReads,
         Jmax = Jmax,
         maxEmissionMatrixDifference = maxEmissionMatrixDifference,
         run_fb_subset = FALSE,
         run_fb_grid_offset = FALSE,
         blocks_for_output = array(0, c(1, 1)),
-        return_alpha = FALSE,
-        return_extra = FALSE,
-        return_genProbs = return_genProbs,
-        return_gamma = as.logical(return_gamma | make_plots),
-        return_hapProbs = return_hapProbs,
-        return_p_store = return_p_store, ## return_p_store
         grid = grid,
         pass_in_alphaBeta = TRUE,
         alphaHat_t1 = alphaHat_t1,
@@ -2019,14 +2042,13 @@ impute_one_sample <- function(
         shuffle_bin_radius = shuffle_bin_radius,
         L_grid = L_grid,
         block_gibbs_iterations = block_gibbs_iterations,
-        return_gibbs_block_output = return_gibbs_block_output,
-        return_advanced_gibbs_block_output = return_advanced_gibbs_block_output,
         rescale_eMatRead_t = rescale_eMatRead_t,
         smooth_cm = smooth_cm,
+        param_list = param_list,
         use_smooth_cm_in_block_gibbs = use_smooth_cm_in_block_gibbs,
-        block_gibbs_quantile_prob = block_gibbs_quantile_prob
+        block_gibbs_quantile_prob = block_gibbs_quantile_prob,
+        use_small_eHapsCurrent_tc = use_small_eHapsCurrent_tc
     )
-
     ##
     genProbs_t <- out$genProbsM_t
     hapProbs_t <- out$happrobs_t

QUILT/R/quilt.R

@@ -64,6 +64,7 @@
 #' @param block_gibbs_iterations What iterations to perform block Gibbs sampling for the Gibbs sampler
 #' @param n_gibbs_burn_in_its How many iterations to run the Gibbs sampler for each time it is run
 #' @param plot_per_sample_likelihoods Plot per sample likelihoods i.e. the likelihood as the method progresses through the Gibbs sampling iterations
+#' @param use_small_eHapsCurrent_tc For testing purposes only
 #' @return Results in properly formatted version
 #' @author Robert Davies
 #' @export
@@ -131,7 +132,8 @@ QUILT <- function(
     print_extra_timing_information = FALSE,
     block_gibbs_iterations = c(3,6,9),
     n_gibbs_burn_in_its = 20,
-    plot_per_sample_likelihoods = FALSE
+    plot_per_sample_likelihoods = FALSE,
+    use_small_eHapsCurrent_tc = TRUE
 ) {
 
 
@@ -579,7 +581,11 @@ QUILT <- function(
         betaHat_t3 <- array(0, c(K, nGrids))
         small_priorCurrent_m <- array(1 / K, c(K, S))
         small_alphaMatCurrent_tc <- array(1 / K, c(K, nGrids - 1, S))
-        small_eHapsCurrent_tc <- array(0, c(K, nSNPs, S))
+        if (use_small_eHapsCurrent_tc) {
+            small_eHapsCurrent_tc <- array(0, c(K, nSNPs, S))
+        } else {
+            small_eHapsCurrent_tc <- array(0, c(1, 1, 1))
+        }
 
         results_across_samples <- as.list(sampleRange[2] - sampleRange[1] + 1)
 
@@ -661,7 +667,8 @@ QUILT <- function(
                 print_extra_timing_information = print_extra_timing_information,
                 n_gibbs_burn_in_its = n_gibbs_burn_in_its,
                 block_gibbs_iterations = block_gibbs_iterations,
-                plot_per_sample_likelihoods = plot_per_sample_likelihoods
+                plot_per_sample_likelihoods = plot_per_sample_likelihoods,
+                use_small_eHapsCurrent_tc = use_small_eHapsCurrent_tc
             )
 
             results_across_samples[[iSample - sampleRange[1] + 1]] <- out

QUILT/R/reference-single.R

@@ -366,7 +366,8 @@ make_rhb_t_equality <- function(
     rhb_t,
     nSNPs,
     ref_error,
-    nMaxDH = NA
+    nMaxDH = NA,
+    verbose = TRUE
 ) {
     if (is.na(nMaxDH)) {
         nMaxDH_default <- 2 ** 8 - 1
@@ -417,7 +418,9 @@ make_rhb_t_equality <- function(
             max(c(2 ** 4 - 1, suggested_value)),
             nMaxDH_default
         )
-        print_message(paste0("Using nMaxDH = ", nMaxDH))
+        if (verbose) {
+            print_message(paste0("Using nMaxDH = ", nMaxDH))
+        }
         distinctHapsB <- distinctHapsB[1:nMaxDH, ]
         hapMatcher[hapMatcher > (nMaxDH)] <- 0
     }