stufield · May 20, 2019 21:49
diff --git a/reprex_ssDAMN.md b/reprex_ssDAMN.md
diff --git a/ssDAMN.cpp b/ssDAMN.cpp
 #include <numeric>
 #include <Rcpp.h>
 using namespace Rcpp;

 /*
 intended to be used in conjunction with apply(x, 1, ssDAMN, ...)
 to loop over samples/rows of the data.frame passed 'x'.
 */

 // [[Rcpp::export]]
 List ssDAMN(const NumericVector& x,              // named vector length ~ 5500
            NumericVector& ref_median,     // named vector length ~ 5500
            const NumericVector& ref_median_log, // named vector ~ 5500
            const NumericVector& ref_mad,        // named vector ~ 5500
            List dil_split,                      // list of characters length = 3
            List dil_split_protein,              // list of characters length = 3
            double effect_cut) {                 // effect size default = 2.0

  CharacterVector apts = x.names();
  int L = x.length();
  int max = 100;
  CharacterVector dil_names = dil_split.names();
  int dil_L = dil_names.length();
  NumericMatrix SF_mat(max, dil_L);
  SF_mat.attr("dimnames") = List::create(R_NilValue, dil_names);
  NumericVector apt_DAMN(L);
  CharacterVector apts_old(L);
  NumericVector sample_vec = clone(x);

  /* Initialize empty objects for use inside loop */
  NumericVector dat_effect(L);
  LogicalVector cut(L);
  CharacterVector apts_cut(L);
  CharacterVector apts_use(L);
  bool lgl1 = NA_LOGICAL;
  bool lgl2 = NA_LOGICAL;
  CharacterVector mix_prot(L);
  CharacterVector inter(L);
  IntegerVector match_id;
  IntegerVector sample_id(L);
  NumericVector new_vals(L);
  double SF;

  /* list of indices for each mix */
  List mix_idx(dil_L);
  CharacterVector mix_apts;
  for (int i = 0; i < dil_L; ++i) {
    mix_apts = dil_split[i];
    mix_idx[i] = match(mix_apts, apts) - 1;
  }

  /* The main loop; convergence */
  int counter = 0;
  while( counter <= max ) {     // typically converge in < 15
    //Rcout << "iter: " << counter << std::endl;
    dat_effect = (log10(sample_vec) - ref_median_log) / ref_mad;   // vectorized effect size ~ 5500
    cut        = abs(dat_effect) > effect_cut;
    apts_cut   = apts[cut];
    apts_use   = setdiff(apts, apts_cut);
    /*
    Rcout << "cut: " << apts_cut << std::endl;
    Rcout << "use: " << apts_use << std::endl;
    Rcout << "old: " << apts_old << std::endl;
    */
    lgl1 = apts_cut.size() == apts_old.size(); // check lengths equal
    lgl2 = is_true(all(apts_cut == apts_old)); // check all equal
    if (counter == max) {
      apt_DAMN = dat_effect;
      break;
    } else if (lgl1 & lgl2) {
      apt_DAMN = dat_effect;
      break;
    } else {
      apts_old = apts_cut;
      counter += 1;          // iterate the counter
    }

    for (int k = 0; k < dil_L; ++k) {   // k = 0:2
      mix_prot = dil_split_protein[k];
      inter    = intersect(apts_use, mix_prot);
      match_id = match(inter, apts) - 1;
      /* calculate scale factor for current dilution mix */
      SF = median(ref_median[ match_id ] / sample_vec[ match_id ], true);
      //Rcout << "SF: " << SF << std::endl;
      sample_id = mix_idx[k];                                      // get indices from list to scale
      new_vals = as<NumericVector>(sample_vec[sample_id]) * SF;    // scale values
      sample_vec[sample_id] = new_vals;  // update sample_vec
      /*
      IntegerVector::iterator j;
      IntegerVector::iterator j_end = sample_id.end();
      for (j = sample_id.begin(); j != j_end; ++j)
        sample_vec[*j] = sample_vec[*j] * SF;
      */
      SF_mat(counter - 1, k) = SF;       // store scale factor in SF_mat
    }
  }
  /* end loops */

  /* Clean up and store relevant objects */
  /* Set up SF_vec and its names */
  CharacterVector sf_names(dil_L, "NormScale.");
  for (int i = 0; i < dil_L; ++i) {
    sf_names[i] += dil_names[i];
  }

  /* replaces apply(SF_mat, 2, prod); better option? */
  NumericVector SF_vec(dil_L);
  NumericVector v(max);
  for (int p = 0; p < dil_L; ++p) {    // length dil_L = 3
    v = SF_mat.column( p );
    v = v[ v > 0 ];
    SF_vec[p] = std::accumulate(v.begin(), v.end(), 1.0, std::multiplies<double>());
  }

  SF_vec.names() = sf_names;
  apt_DAMN.names() = apts;

  return List::create(Named("MedDat", sample_vec),
                      Named("DAMN", apt_DAMN),
                      Named("SF", SF_vec));
 }



 /***R
  # set.seed(100)
  # n   <- 25
  # nms <- LETTERS[1:n]
  # mat <- matrix(rnorm(5 * n, mean = 100), ncol = n, dimnames = list(NULL, nms))
  # a   <- apply(mat, 2, median)
  # b   <- n:1 %>% purrr::set_names(nms)
  # c   <- rexp(n) %>% purrr::set_names(nms)
  # dil <- split(nms, sample(1:3, n, replace = TRUE))
  # t0  <- Sys.time()
  # ret <- apply(mat, 1, ssDAMN, ref_median = a,
  #              ref_median_log = c, ref_mad = b,
  #              dil_split = dil, dil_split_protein = dil,
  #              effect_cut = 0.5) %>%
  #   purrr::transpose() %>%
  #   purrr::map(~ do.call(rbind, .x) %>% data.frame())
  # t1 <- Sys.time()
  # all.equal(ret, readRDS("ret.rds")); t1 - t0
 */

 /***R
  df <- system.file("medianNorm_references",
                    "Population_Reference_V4_Covance_plasma.rds",
                    package = "SomaNormalization", mustWork = TRUE) %>%
    readRDS() %>%
    filter(SITE_ID == "Boise" & BMI < 21) %>%
    dplyr::select(SampleType, SampleMatrix, SlideId, Subarray, SampleId, SampleGroup,
           getNormNames(.), sample(getAptamers(.), 500))
  mark <- bench::mark(
    orig_ssDAMN = suppressMessages(singleSampleDAMN(df)),   # original
    cpp_ssDAMN  = suppressMessages(singleSampleDAMN2(df)),  # uses Rcpp
    filter_gc = FALSE, relative = F, iterations = 5
  )
  mark
  #summary(bnch, relative = TRUE)
 */
	#include <numeric>
	#include <Rcpp.h>
	using namespace Rcpp;

	/*
	intended to be used in conjunction with apply(x, 1, ssDAMN, ...)
	to loop over samples/rows of the data.frame passed 'x'.
	*/

	// [[Rcpp::export]]
	List ssDAMN(const NumericVector& x, // named vector length ~ 5500
	NumericVector& ref_median, // named vector length ~ 5500
	const NumericVector& ref_median_log, // named vector ~ 5500
	const NumericVector& ref_mad, // named vector ~ 5500
	List dil_split, // list of characters length = 3
	List dil_split_protein, // list of characters length = 3
	double effect_cut) { // effect size default = 2.0

	CharacterVector apts = x.names();
	int L = x.length();
	int max = 100;
	CharacterVector dil_names = dil_split.names();
	int dil_L = dil_names.length();
	NumericMatrix SF_mat(max, dil_L);
	SF_mat.attr("dimnames") = List::create(R_NilValue, dil_names);
	NumericVector apt_DAMN(L);
	CharacterVector apts_old(L);
	NumericVector sample_vec = clone(x);

	/* Initialize empty objects for use inside loop */
	NumericVector dat_effect(L);
	LogicalVector cut(L);
	CharacterVector apts_cut(L);
	CharacterVector apts_use(L);
	bool lgl1 = NA_LOGICAL;
	bool lgl2 = NA_LOGICAL;
	CharacterVector mix_prot(L);
	CharacterVector inter(L);
	IntegerVector match_id;
	IntegerVector sample_id(L);
	NumericVector new_vals(L);
	double SF;

	/* list of indices for each mix */
	List mix_idx(dil_L);
	CharacterVector mix_apts;
	for (int i = 0; i < dil_L; ++i) {
	mix_apts = dil_split[i];
	mix_idx[i] = match(mix_apts, apts) - 1;
	}

	/* The main loop; convergence */
	int counter = 0;
	while( counter <= max ) { // typically converge in < 15
	//Rcout << "iter: " << counter << std::endl;
	dat_effect = (log10(sample_vec) - ref_median_log) / ref_mad; // vectorized effect size ~ 5500
	cut = abs(dat_effect) > effect_cut;
	apts_cut = apts[cut];
	apts_use = setdiff(apts, apts_cut);
	/*
	Rcout << "cut: " << apts_cut << std::endl;
	Rcout << "use: " << apts_use << std::endl;
	Rcout << "old: " << apts_old << std::endl;
	*/
	lgl1 = apts_cut.size() == apts_old.size(); // check lengths equal
	lgl2 = is_true(all(apts_cut == apts_old)); // check all equal
	if (counter == max) {
	apt_DAMN = dat_effect;
	break;
	} else if (lgl1 & lgl2) {
	apt_DAMN = dat_effect;
	break;
	} else {
	apts_old = apts_cut;
	counter += 1; // iterate the counter
	}

	for (int k = 0; k < dil_L; ++k) { // k = 0:2
	mix_prot = dil_split_protein[k];
	inter = intersect(apts_use, mix_prot);
	match_id = match(inter, apts) - 1;
	/* calculate scale factor for current dilution mix */
	SF = median(ref_median[ match_id ] / sample_vec[ match_id ], true);
	//Rcout << "SF: " << SF << std::endl;
	sample_id = mix_idx[k]; // get indices from list to scale
	new_vals = as<NumericVector>(sample_vec[sample_id]) * SF; // scale values
	sample_vec[sample_id] = new_vals; // update sample_vec
	/*
	IntegerVector::iterator j;
	IntegerVector::iterator j_end = sample_id.end();
	for (j = sample_id.begin(); j != j_end; ++j)
	sample_vec[j] = sample_vec[j] * SF;
	*/
	SF_mat(counter - 1, k) = SF; // store scale factor in SF_mat
	}
	}
	/* end loops */

	/* Clean up and store relevant objects */
	/* Set up SF_vec and its names */
	CharacterVector sf_names(dil_L, "NormScale.");
	for (int i = 0; i < dil_L; ++i) {
	sf_names[i] += dil_names[i];
	}

	/* replaces apply(SF_mat, 2, prod); better option? */
	NumericVector SF_vec(dil_L);
	NumericVector v(max);
	for (int p = 0; p < dil_L; ++p) { // length dil_L = 3
	v = SF_mat.column( p );
	v = v[ v > 0 ];
	SF_vec[p] = std::accumulate(v.begin(), v.end(), 1.0, std::multiplies<double>());
	}

	SF_vec.names() = sf_names;
	apt_DAMN.names() = apts;

	return List::create(Named("MedDat", sample_vec),
	Named("DAMN", apt_DAMN),
	Named("SF", SF_vec));
	}



	/***R
	# set.seed(100)
	# n <- 25
	# nms <- LETTERS[1:n]
	# mat <- matrix(rnorm(5 * n, mean = 100), ncol = n, dimnames = list(NULL, nms))
	# a <- apply(mat, 2, median)
	# b <- n:1 %>% purrr::set_names(nms)
	# c <- rexp(n) %>% purrr::set_names(nms)
	# dil <- split(nms, sample(1:3, n, replace = TRUE))
	# t0 <- Sys.time()
	# ret <- apply(mat, 1, ssDAMN, ref_median = a,
	# ref_median_log = c, ref_mad = b,
	# dil_split = dil, dil_split_protein = dil,
	# effect_cut = 0.5) %>%
	# purrr::transpose() %>%
	# purrr::map(~ do.call(rbind, .x) %>% data.frame())
	# t1 <- Sys.time()
	# all.equal(ret, readRDS("ret.rds")); t1 - t0
	*/

	/***R
	df <- system.file("medianNorm_references",
	"Population_Reference_V4_Covance_plasma.rds",
	package = "SomaNormalization", mustWork = TRUE) %>%
	readRDS() %>%
	filter(SITE_ID == "Boise" & BMI < 21) %>%
	dplyr::select(SampleType, SampleMatrix, SlideId, Subarray, SampleId, SampleGroup,
	getNormNames(.), sample(getAptamers(.), 500))
	mark <- bench::mark(
	orig_ssDAMN = suppressMessages(singleSampleDAMN(df)), # original
	cpp_ssDAMN = suppressMessages(singleSampleDAMN2(df)), # uses Rcpp
	filter_gc = FALSE, relative = F, iterations = 5
	)
	mark
	#summary(bnch, relative = TRUE)
	*/