Rcpp code improvement

[Prof-ThiagoOliveira]

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;
}

[gregorgorjanc]

@RosCraddock let's work together on this issue/suggestion at some point - remind me to book one longer working slot and we can work on it together

[gregorgorjanc]

Now we are using for loop counters as integers instead of integers for better performance.

ChatGPT says: No — it isn’t faster. a good optimizer will emit essentially the same machine code. The only real difference is scope/lifetime (the second i dies at the end of the loop).

ped(i, j) to ped(i-1, j) because NumericMatrix indexing in Rcpp starts at 0.

Actually, we use this as a useful trick. Unknown individuals are coded as 0 and we use 0-th row for them where we have value 0 for them.

I see merit in other code. Let me do some testing to ensure all is correct.

[gregorgorjanc]

@Prof-ThiagoOliveira yeah, using // [[Rcpp::export]] makes everything much simpler!. Done now in 723c7be.