The R Project SVN R-packages

                /* Methods for pedigrees */
#include "pedigree.h"

/**
 * Create the left Cholesky factor of the numerator relationship
 * matrix from a pedigree.
 *
 * @param x a pedigree object
 * @param tinv the inverse of T, a unit lower dtCMatrix
 * @param ttrans the transpose of T, a unit upper dtCMatrix
 * @param ans T stored as a non-unit, lower dtCMatrix
 * the pedigree
 *
 * @return ans with elements modified to incorporate D
 */
SEXP pedigree_chol(SEXP x, SEXP ans)
{
    SEXP Sire = GET_SLOT(x, install("sire"));
    int *ai = INTEGER(GET_SLOT(ans, Matrix_iSym)),
    *ap = INTEGER(GET_SLOT(ans, Matrix_pSym)),
    *dam = INTEGER(GET_SLOT(x, install("dam"))),
    *sire = INTEGER(Sire), 
    i, j, n = LENGTH(Sire);
    double *ax = REAL(GET_SLOT(ans, Matrix_xSym)), *F, Di, tmp;

    setAttrib(ans, install("F"), allocVector(REALSXP, n));
    F = REAL(getAttrib(ans, install("F")));
    for (i = 0; i < n; i++) {
    int k, p = sire[i] - 1, q = dam[i] - 1;
    if (sire[i] == NA_INTEGER) {
        F[i] = 0;
        Di = (dam[i] == NA_INTEGER) ? 1 : sqrt(0.75 - 0.25 * F[q]);
    } else {
        if (dam[i] == NA_INTEGER) { /* sire only */
        F[i] = 0;
        Di = sqrt(0.75 - 0.25 * F[p]);
        } else {        /* both parents in pedigree */
        Di = sqrt(0.5 - 0.25 * (F[p] + F[q]));
        F[i] = NA_REAL;
        if ((ap[i + 1] - ap[i]) > 1) {    /* skip if no progeny */
            if (p > q) {j = p; p = q; q = j;} /* ensure p <= q */
            for (j = 0, F[i] = 0; j <= p; j++) {
            for (k = ap[j], tmp = 0;
                 k < ap[j + 1] && ai[k] <= q; k++) {
                int ii = ai[k];
                if (ii == p) tmp = ax[k];
                if (ii == q) F[i] += tmp * ax[k]/2;
            }
            }
        }
        }
    }
    for (j = ap[i]; j < ap[i + 1]; j++) ax[j] *= Di;
    }
    return ans;
}

/* NOTE: This function requires that missing parents be coded as zero */
/**
 * Create the inbreeding coefficients according to the algorithm given
 * in "Comparison of four direct algorithms for computing inbreeding
 * coefficients" by Mehdi Sargolzaei and Hiroaki Iwaisaki, Animal
 * Science Journal (2005) 76, 401--406.  This function is a modified
 * version of the code published in an appendix to that paper.
 *
 * @param x a pedigree object
 *
 * @return a list of the inbreeding coefficients
 */

SEXP pedigree_inbreeding(SEXP x)
{
    SEXP ans, sp = GET_SLOT(x, install("sire"));
    int i, j, t, n = LENGTH(sp), S, D;
    int *sire = INTEGER(sp),
    *dam = INTEGER(GET_SLOT(x, install("dam"))),
    *Anc = Calloc(n + 1, int), /* ancestor */
    *SI, *MI;          /* start and minor */
    double *F = Calloc(n + 1, double), /* inbreeding coefficients */
      *L = Calloc(n + 1, double),
      *B = Calloc(n + 1, double); 
    int *LAP = Calloc(n + 1, int);  /* longest ancestoral path */
    
    F[0] =-1; LAP[0] =-1; /* set F and lap for unknown parents */
    for(i = 1, t = -1; i <= n; i++) {   /* evaluate LAP and its maximum */
    S = sire[i]; D = dam[i]; /* parents of animal i */
    LAP[i] = ((LAP[S] < LAP[D]) ? LAP[D] : LAP[S]) + 1;
    if (LAP[i] > t) t = LAP[i];
    }
    SI = Calloc(t + 1, int);
    MI = Calloc(t + 1, int);
    for(i = 0; i <= t ; ++i) SI[i] = MI[i] = 0; /* initialize start and minor */
    for(i = 1; i <= n; i++) {   /* evaluate F */
    S = sire[i]; D = dam[i]; /* parents of animal i */
    B[i] = 0.5 - 0.25 * (F[S] + F[D]); 
                /* adjust start and minor */
    for (j = 0; j < LAP[i]; j++) {++SI[j]; ++MI[j];} 
    if (S == 0 || D == 0) { /* both parents unknown */
        F[i] = L[i] = 0; continue;
    }
    if(S == sire[i-1] && D == dam[i-1]) { /* full-sib with last animal */
        F[i] = F[i-1]; L[i] = L[i-1]; continue;
    }
    
    F[i] = -1; L[i] = 1; 
    t = LAP[i]; /* largest lap group number in the animal's pedigree */
    Anc[MI[t]++] = i; /* initialize Anc and increment MI[t] */
    while(t > -1) { /* from the largest lap group to zero */
        j = Anc[--MI[t]]; /* next ancestor */
        S = sire[j]; D = dam[j]; /* parents of the ancestor */
        if (S) {
        if (!L[S]) Anc[MI[LAP[S]]++] = S; 
                /* add sire in its lap group in Anc
                 * array if it is not added yet and
                 * increment the minor index for the group */ 
        L[S] += 0.5 * L[j]; /* contribution to sire */
        }
        if (D) {
        if (!L[D]) Anc[MI[LAP[D]]++] = D;
        L[D] += 0.5 * L[j]; /* contribution to dam */
        }
        F[i] += L[j] * L[j] * B[j];
        L[j] = 0; /*clear L[j] for the evaluation of the next animal */
        if (MI[t] == SI[t]) --t; /* move to the next lap group when
                      * all ancestors in group t have been
                      * evaluated */
    } 
    }
    ans = PROTECT(allocVector(REALSXP, n));
    Memcpy(REAL(ans), F + 1, n);
    Free(Anc); Free(F); Free(L); Free(B); Free(LAP); Free(SI); Free(MI);
    UNPROTECT(1);
    return ans;
}