[R-genetics] SF.net SVN: r-genetics: [1198] trunk/GeneticsPed

Revision: 1198
          http://svn.sourceforge.net/r-genetics/?rev=1198&view=rev
Author:   dnadave
Date:     2007-04-09 11:28:31 -0700 (Mon, 09 Apr 2007)

Log Message:
-----------
Strange error, I think,

if you have a pointer to an int set within a class and you want to modify it via:

delete p;
p = new int;
*p = newp;

you get a double free error.  What does work is:

p = new int;
*p = newp;

Is this a bug in the compiler?

Also, I added in some new sorting routines that require only the pedigree information.

Modified Paths:
--------------
    trunk/GeneticsPed/R/inverseAdditive.R
    trunk/GeneticsPed/R/sort.pedigree.R
    trunk/GeneticsPed/inst/doc/pedigreeHandling.Rnw
    trunk/GeneticsPed/src/inverseAdditive.cc
    trunk/GeneticsPed/src/pedtemplate.cc
    trunk/GeneticsPed/src/register.cc
    trunk/GeneticsPed/src/sortped.cc

Added Paths:
-----------
    trunk/GeneticsPed/include/pedSort.h
    trunk/GeneticsPed/src/pedSort.cc

Modified: trunk/GeneticsPed/R/inverseAdditive.R
===================================================================
--- trunk/GeneticsPed/R/inverseAdditive.R	2007-04-06 12:03:00 UTC (rev 1197)
+++ trunk/GeneticsPed/R/inverseAdditive.R	2007-04-09 18:28:31 UTC (rev 1198)
@@ -2,7 +2,7 @@
 ###------------------------------------------------------------------------
 ### What: Inverse of Additive Relationship Matrix
 ### $Id$
-### Time-stamp: <2007-04-05 21:22:59 ggorjan>
+### Time-stamp: <2007-04-03 16:01:34 ggorjan>
 ###------------------------------------------------------------------------
 
 inverseAdditive <- function(x, sort=TRUE, names=TRUE, ...)
@@ -15,13 +15,17 @@
     stop("no method for pedigree with more than two ascendants")
   if(sort) idOrig <- as.character(x[[subject]])  # for sorting-back
 
-  ## FIXME: is all this needed! - check Dave's code
-  ## Pedigree must be sorted, extended
-  x <- GeneticsPed:::checkAttributes(x, sorted=TRUE, extended=TRUE)
+  ## FIXME: all this is not needed! - check Dave's code
+  ## Pedigree must be sorted, extended and unknown equal to NA
+  ##x <- GeneticsPed:::checkAttributes(x, sorted=TRUE, extended=TRUE, unknownNA=TRUE)
+  ## I've commented out the above as my code does sort and extend, but does require
+  ##  that unknowns be coded correctly.  This will likely save considerable time
+  ##  for large pedigrees.
+  x <- GeneticsPed:::checkAttributes(x, sorted=FALSE, extended=FALSE, unknownNA=TRUE)
 
   ## Find unknown value (NA, "NA", 0, ...)
   ## paste is used due to possibility of having NA and as.character returns
-  ## NA and not "NA" i.e. paste(NA) always returns "NA"
+  ## NA and not "NA" i.e. paste(NA) returns "NA"
   na.value <- paste(as.character(attr(x, ".unknown")$.id))
 
   ## --- Core ---

Modified: trunk/GeneticsPed/R/sort.pedigree.R
===================================================================
--- trunk/GeneticsPed/R/sort.pedigree.R	2007-04-06 12:03:00 UTC (rev 1197)
+++ trunk/GeneticsPed/R/sort.pedigree.R	2007-04-09 18:28:31 UTC (rev 1198)
@@ -29,8 +29,32 @@
 
   ## --- Sort on different information ---
 
-  if(by == "pedigree") {
-    warning("sorting by pedigree information only is not yet implemented")
+  if(by == "pedigree")
+  {
+    subject <- attr(x, ".subject")
+    ascendant <- attr(x, ".ascendant")
+    na.value <- paste(as.character(attr(x, ".unknown")$.id))
+    n.na <- length(na.value)
+    n <- nrow(x)
+    ret <- .C(name=R_pedSort,
+        ## 1 number of individuals
+        as.integer(n),
+        ## 2 individual ID
+        as.character(x[[subject]]),
+        ## 3 father ID
+        as.character(x[[ascendant[1]]]),
+        ## 4 mother ID
+        as.character(x[[ascendant[2]]]),
+        ## 5 matrix coefficients
+        ## 6 character string indicating missing values
+        as.character(na.value),
+        ## 7 length of character string indicating missing values
+        as.integer(n.na),
+        PACKAGE="GeneticsPed")
+    x[[subject]] <- ret[[2]]
+    x[[ascendant[1]]] <- ret[[3]]
+    x[[ascendant[2]]] <- ret[[4]]
+    attr(x, ".sorted") <- TRUE
     return(x)
   }
 

Added: trunk/GeneticsPed/include/pedSort.h
===================================================================
--- trunk/GeneticsPed/include/pedSort.h	                        (rev 0)
+++ trunk/GeneticsPed/include/pedSort.h	2007-04-09 18:28:31 UTC (rev 1198)
@@ -0,0 +1,30 @@
+/*                                                                 */
+/* Program to calculate directly the inverse of A                  */
+/* Written by DAH 06/28/2000                                       */
+/*                                                                 */
+/* The program works by reading a pedigree and then                */
+/* forming A^(-1) directly using established rules.                */
+/*                                                                 */
+/* Utility functions for forming A^(-1)                            */
+/*                                                                 */
+
+#ifndef PEDSORT_H
+#define PEDSORT_H
+
+using namespace std;
+
+#include <R.h>
+#include <Rdefines.h>
+#include <string>
+#include <map>
+#include <algorithm>
+#include "pedtemplate.h"
+#include "sortped.h"
+
+extern "C" {
+
+void pedSort(int *n, char **ind, char **father, char **mother, char **na_value, int *n_na);
+
+}
+
+#endif

Modified: trunk/GeneticsPed/inst/doc/pedigreeHandling.Rnw
===================================================================
--- trunk/GeneticsPed/inst/doc/pedigreeHandling.Rnw	2007-04-06 12:03:00 UTC (rev 1197)
+++ trunk/GeneticsPed/inst/doc/pedigreeHandling.Rnw	2007-04-09 18:28:31 UTC (rev 1198)
@@ -101,7 +101,13 @@
 
 \section*{Introduction}
 
-what is pedigree/genealogy
+Pedigrees are collections of related individuals.  Often we represent these as a
+linked list, a collection of trios that links (or almost everyone) everyone
+together: an individual and its two parents.  This simple representation allows
+the use of graph theory in analysis.  The GeneticsPed package provides utilities
+for managing pedigrees; inputing, sorting, and subsetting pedigrees; and
+computing on pedigrees by calculating relationship coefficients and other
+similar quantities.
 
 name some fields were pedigree is used
 \cite{Falconer:1996}
@@ -145,6 +151,110 @@
 1+1
 @
 
+\section{Check consistency of data in pedigree}
+
+check
+
+check.Pedigree
+checkId
+
+check performs a series of checks on
+  pedigree object to ensure consistency of data.
+
+check(x, $\ldots$)
+checkId(x)
+
+
+\begin{itemize}
+  \item[x] pedigree, object to be checked
+  \item[]$\ldots$] arguments to other methods, none for now
+\end{itemize}
+
+  checkId performs various checks on subjects and
+  their ascendants. These checks are:
+\begin{itemize}
+    \item idClass: all ids must have the same class
+    \item subjectIsNA: subject can not be NA
+    \item subjectNotUnique: subject must be unique
+    \item subjectEqualAscendant: subject can not be equal (in
+      identification) to its ascendant
+    \item ascendantEqualAscendant: ascendant can not be equal to another
+      ascendant
+    \item ascendantInAscendant: ascendant can not appear again as
+      asescendant of other sex i.e. father can not be a mother to someone
+      else
+    \item unusedLevels: in case factors are used for id presentation,
+      there might be unused levels for some ids - some functions rely on
+      number of levels and a check is provided for this
+\end{itemize}
+
+  checkAttributes is intended primarly for internal use and
+  performs a series of checks on attribute values needed in various
+  functions. It causes stop with error messages for all given attribute
+  checks.
+
+
+List of more or less self-explanatory errors and "pointers" to
+  these errors for ease of further work i.e. removing errors.
+
+\begin{verbatim}
+  ## EXAMPLES BELLOW ARE ONLY FOR TESTING PURPOSES AND ARE NOT INTENDED
+  ## FOR USERS, BUT IT CAN NOT DO ANY HARM.
+
+  ## --- checkAttributes ---
+  tmp <- generatePedigree(5)
+  attr(tmp, "sorted") <- FALSE
+  attr(tmp, "coded") <- FALSE
+  GeneticsPed:::checkAttributes(tmp)
+  try(GeneticsPed:::checkAttributes(tmp, sorted=TRUE, coded=TRUE))
+
+  ## --- idClass ---
+  tmp <- generatePedigree(5)
+  tmp$id <- factor(tmp$id)
+  class(tmp$id)
+  class(tmp$father)
+  try(GeneticsPed:::idClass(tmp))
+
+  ## --- subjectIsNA ---
+  tmp <- generatePedigree(2)
+  tmp[1, 1] <- NA
+  GeneticsPed:::subjectIsNA(tmp)
+
+  ## --- subjectNotUnique ---
+  tmp <- generatePedigree(2)
+  tmp[2, 1] <- 1
+  GeneticsPed:::subjectNotUnique(tmp)
+
+  ## --- subjectEqualAscendant ---
+  tmp <- generatePedigree(2)
+  tmp[3, 2] <- tmp[3, 1]
+  GeneticsPed:::subjectEqualAscendant(tmp)
+
+  ## --- ascendantEqualAscendant ---
+  tmp <- generatePedigree(2)
+  tmp[3, 2] <- tmp[3, 3]
+  GeneticsPed:::ascendantEqualAscendant(tmp)
+
+  ## --- ascendantInAscendant ---
+  tmp <- generatePedigree(2)
+  tmp[3, 2] <- tmp[5, 3]
+  GeneticsPed:::ascendantInAscendant(tmp)
+  ## Example with multiple parents
+  tmp <- data.frame(id=c("A", "B", "C", "D"),
+                    father1=c("E", NA, "F", "H"),
+                    father2=c("F", "E", "E", "I"),
+                    mother=c("G", NA, "H", "E"))
+  tmp <- Pedigree(tmp, ascendant=c("father1", "father2", "mother"),
+                  ascendantSex=c(1, 1, 2),
+                  ascendantLevel=c(1, 1, 1))
+  GeneticsPed:::ascendantInAscendant(tmp)
+
+  ## --- unusedLevels ---
+  tmp <- generatePedigree(2, colClass="factor")
+  tmp[3:4, 2] <- NA
+  GeneticsPed:::unusedLevels(tmp)
+\end{verbatim}
+
 \bibliographystyle{apalike}
 \bibliography{library}
 

Modified: trunk/GeneticsPed/src/inverseAdditive.cc
===================================================================
--- trunk/GeneticsPed/src/inverseAdditive.cc	2007-04-06 12:03:00 UTC (rev 1197)
+++ trunk/GeneticsPed/src/inverseAdditive.cc	2007-04-09 18:28:31 UTC (rev 1198)
@@ -15,7 +15,7 @@
   TPedVec pedtmp;
   Pedigree Ped;
   vector< string > missing_values;
-  for ( int i = 0; i <*n_na; i++ )
+  for ( int i = 0; i < *n_na; i++ )
   {
     ostringstream ms;
     ms << na_value[ i ];

Added: trunk/GeneticsPed/src/pedSort.cc
===================================================================
--- trunk/GeneticsPed/src/pedSort.cc	                        (rev 0)
+++ trunk/GeneticsPed/src/pedSort.cc	2007-04-09 18:28:31 UTC (rev 1198)
@@ -0,0 +1,81 @@
+#define R_NO_REMAP
+
+#include <R.h>
+#include <Rdefines.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <sstream>
+#include "../include/pedSort.h"
+#include <vector>
+
+extern "C" {
+
+void pedSort( int *n , char **ind , char **father , char **mother , char **na_value , int *n_na )
+{
+  TPedVec pedtmp;
+  Pedigree Ped;
+  vector< string > missing_values;
+  string miss = "";
+  for ( int i = 0; i < *n_na; i++ )
+  {
+    ostringstream ms;
+    ms << na_value[ i ];
+    missing_values.insert( missing_values.end() , ms.str()  );
+    if ( i == 0 )
+    {
+      miss = na_value[ i ];
+    } 
+  }
+  for ( unsigned int i=0 ; i < static_cast< unsigned int >( *n ) ; i++ )
+  {
+    ostringstream kss, pss, mss, tmss, tpss;
+    kss << ind[i];
+    vector< string >::const_iterator miss = missing_values.begin();
+    tpss << father[i];
+    if ( ( miss = find( missing_values.begin() , missing_values.end() , tpss.str() ) ) != missing_values.end() )
+    {
+      pss << "";
+    }
+    else
+    {
+      pss << father[i];
+    }
+    tmss << mother[i];
+    if ( ( miss = find( missing_values.begin() , missing_values.end() , tmss.str() ) ) != missing_values.end() )
+    {
+      mss << "";
+    }
+    else
+    {
+      mss << mother[i];
+    }
+    string kid = kss.str();
+    string pop = pss.str();
+    string mom = mss.str();
+    pedtmp.insert( pedtmp.end() , TPed( kid , pop , mom , static_cast< int >( i ) + 1 ) );
+  }
+  SortPed( Ped , pedtmp );
+  missing_values.erase( missing_values.begin() , missing_values.end() );
+  for ( unsigned int i=0 ; i < static_cast< unsigned int >( *n ) ; i++ )
+  {
+    ind[i] = const_cast< char* >( Ped.GetMember( i ).c_str() );
+    if ( Ped.GetParent( 0 , i ) )
+    {
+      father[i] = const_cast< char* >( Ped.GetMember( Ped.GetParentIndex( i , SIRE ) ).c_str() );
+    }
+    else
+    {
+      father[i] = const_cast< char* >( miss.c_str() );
+    }
+    if ( Ped.GetParent( 0 , i ) )
+    {
+      mother[i] = const_cast< char* >( Ped.GetMember( Ped.GetParentIndex( i , DAM ) ).c_str() );
+    }
+    else
+    {
+      mother[i] = const_cast< char* >( miss.c_str() );
+    }
+  }
+}
+
+}

Modified: trunk/GeneticsPed/src/pedtemplate.cc
===================================================================
--- trunk/GeneticsPed/src/pedtemplate.cc	2007-04-06 12:03:00 UTC (rev 1197)
+++ trunk/GeneticsPed/src/pedtemplate.cc	2007-04-09 18:28:31 UTC (rev 1198)
@@ -221,11 +221,17 @@
 int TPed::IsBase()
 {
   if ( !(hasparents) )
+  {
     return 1;
+  }
   else if ( sire != "" && dam != "" )
+  {
     return 0;
+  }
   else
+  {
     return -1;
+  }
 }
 
 void TPed::SetIndex( int &index , int par )
@@ -234,7 +240,7 @@
   {
     if ( index >= 0 )
     {
-      delete s_index;
+//      delete s_index;
       s_index = new int;
       *s_index = index;
     }
@@ -243,7 +249,7 @@
   {
     if ( index >= 0 )
     {
-      delete d_index;
+//      delete d_index;
       d_index = new int;
       *d_index = index;
     }
@@ -254,15 +260,31 @@
 {
   if ( par == SIRE )
   {
-    delete s_index;
-    s_index = new int;
-    *s_index = index;
+    if ( s_index )
+    {
+//      delete s_index;
+      s_index = new int;
+      *s_index = index;
+    }
+    else
+    {
+      s_index = new int;
+      *s_index = index;
+    }
   }
   else
   {
-    delete d_index;
-    d_index = new int;
-    *d_index = index;
+    if ( d_index )
+    {
+//      delete d_index;
+      d_index = new int;
+      *d_index = index;
+    }
+    else
+    {
+      d_index = new int;
+      *d_index = index;
+    }
   }
 }
 
@@ -277,16 +299,24 @@
   if ( par == SIRE )
   {
     if ( s_index )
+    {
       index = *s_index;
+    }
     else
+    {
       index = -1;
+    }
   }
   else
   {
     if ( d_index )
+    {
       index = *d_index;
+    }
     else
+    {
       index = -1;
+    }
   }
   return index;
 }
@@ -296,22 +326,33 @@
   if ( par == SIRE )
   {
     if ( sire != "" )
+    {
       return true;
+    }
     else
+    {
       return false;
+    }
   }
   else
   {
     if ( dam != "" )
+    {
       return true;
+    }
     else
+    {
       return false;
+    }
   }
 }
 
 void TPed::operator= ( const TPed &copy )
 {
-  if ( this == &copy ) return;
+  if ( this == &copy )
+  {
+    return;
+  }
   delete s_index;
   delete d_index;
   s_index = new int;
@@ -373,9 +414,13 @@
 {
   string tmp = arg;
   if ( animal == tmp )
+  {
     return true;
+  }
   else
+  {
     return false;
+  }
 }
 
 void TPed::ShowPed()
@@ -441,16 +486,24 @@
   if ( p == 0 )
   {
     if ( pedigree[i].GetIndex( SIRE ) >= 0 )
+    {
       return true;
+    }
     else
+    {
       return false;
+    }
   }
   else
   {
     if ( pedigree[i].GetIndex( DAM ) >= 0 )
+    {
       return true;
+    }
     else
+    {
       return false;
+    }
   }
 }
 
@@ -480,7 +533,11 @@
 {
   TPedVec::iterator p = find( pedigree.begin() , pedigree.end() , a );
   if ( p != pedigree.end() )
+  {
     return p - pedigree.begin();
+  }
   else
-    return 99;
+  {
+    return -1;
+  }
 }

Modified: trunk/GeneticsPed/src/register.cc
===================================================================
--- trunk/GeneticsPed/src/register.cc	2007-04-06 12:03:00 UTC (rev 1197)
+++ trunk/GeneticsPed/src/register.cc	2007-04-09 18:28:31 UTC (rev 1198)
@@ -7,6 +7,7 @@
 
 #include "../include/meuwissen.h"
 #include "../include/ainverse.h"
+#include "../include/pedSort.h"
 
 extern "C" {
 
@@ -46,11 +47,21 @@
   INTSXP   // 8 int    **n_na
 };
 
+static R_NativePrimitiveArgType pedSortType[] = {
+  INTSXP,  // 1 int    *n
+  STRSXP,  // 2 char   **ind
+  STRSXP,  // 3 char   **father
+  STRSXP,  // 4 char   **mother
+  STRSXP,  // 5 char   **na_value
+  INTSXP   // 6 int    **n_na
+};
+
 // Functions accessed via .C()
 
 static const R_CMethodDef cMethods[] = {
   {"meuwissen",       (DL_FUNC) &meuwissen,       8, meuwissenType},
   {"inverseAdditive", (DL_FUNC) &inverseAdditive, 8, inverseAdditiveType},
+  {"pedSort",         (DL_FUNC) &pedSort,         6, pedSortType},
   {NULL, NULL, 0}
 };
 

Modified: trunk/GeneticsPed/src/sortped.cc
===================================================================
--- trunk/GeneticsPed/src/sortped.cc	2007-04-06 12:03:00 UTC (rev 1197)
+++ trunk/GeneticsPed/src/sortped.cc	2007-04-09 18:28:31 UTC (rev 1198)
@@ -145,9 +145,13 @@
           all.insert( all.end() , *q );
           parents.erase( q );
           if ( q + 1 <= parents.end() )
+          {
             q++;
+          }
           else
+          {
             q = parents.begin();
+          }
         }
         else if ( q + 1 <= parents.end() )
         {


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.