Rcpp code improvement - Githubissues

Is your feature request related to a problem? Please describe. A possible improvement to the Rcpp functions.

Describe the solution you'd like

#include "AlphaPartDrop.h"
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List AlphaPartDrop(double c1, double c2, int nI, int nP, int nT, NumericMatrix ped, IntegerVector P, IntegerVector Px) {

  // --- Outputs ---
  NumericMatrix pa(nI+1, nT);    // parent average
  NumericMatrix w(nI+1, nT);     // Mendelian sampling
  NumericMatrix xa(nI+1, nP*nT); // Parts

  // --- Compute ---
  for(int i = 1; i < nI+1; i++) {
    for(int t = 0; t < nT; t++) {
      // Parent average (PA)
      pa(i, t) = c1 * ped(i-1, 3+t) + c2 * ped(i-1, 3+nT+t);

      // Mendelian sampling (MS)
      w(i, t) = ped(i-1, 3+t) - pa(i, t);

      // Parts

      // ... for the MS part
      int j = Px[t] + P[i];
      xa(i, j) = w(i, t);

      // ... for the PA parts
      for(int p = 0; p < nP; p++) {
        j = Px[t] + p;
        xa(i, j) += c1 * xa(ped(i-1, 3), j) + c2 * xa(ped(i-1, 3+nT), j);
      }
    }
  }

  // --- Return ---
  return List::create(Named("pa") = pa, Named("w") = w, Named("xa") = xa);
}

Changes:

Now we are using for loop counters as integers instead of integers for better performance.
ped(i, j) to ped(i-1, j) because NumericMatrix indexing in Rcpp starts at 0.
return statement by using the List::create function. (reason on page 6 - https://mran.microsoft.com/snapshot/2016-08-12/web/packages/Rcpp/vignettes/Rcpp-modules.pdf)

Same for AlphaPartDropGroup:

#include "AlphaPartDropGroup.h"
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericMatrix AlphaPartDropGroup(double c1, double c2, int nI, int nP, int nT, int nG,
                                 NumericMatrix ped, IntegerVector P, IntegerVector Px, IntegerVector g) {

  // --- Outputs ---

  NumericMatrix pa(nI+1, nT);    // parent average
  NumericMatrix w(nI+1, nT);     // Mendelian sampling
  NumericMatrix xa(nI+1, nP*nT); // Parts
  NumericMatrix xg(nG+1, nP*nT); // Parts for groups

  // --- Compute ---

  for (int i = 1; i < nI+1; i++) {
    for (int t = 0; t < nT; t++) {
      // Parent average (PA)
      pa(i, t) = c1 * ped(ped(i, 1), 3+t) + c2 * ped(ped(i, 2), 3+t);

      // Mendelian sampling (MS)
      w(i, t) = ped(i, 3+t) - pa(i, t);

      // Parts
      int j = Px[t] + P[i];

      // ... for the MS part
      xa(i, j) = w(i, t);

      // ... for the PA parts
      for (int p = 0; p < nP; p++) {
        j = Px[t] + p;
        xa(i, j) += c1 * xa(ped(i, 1), j) + c2 * xa(ped(i, 2), j);

        // Accumulate parts by group
        xg(g(i), j) += xa(i, j);
      }
    }
  }

  // --- Return ---

  return xg;
}

AlphaGenes / AlphaPart

Rcpp code improvement #3