alphagenes / alphapart Goto Github PK

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

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