/auteur/src/vmat.cpp

  1#include <Rcpp.h>
  2#include <vector>
  3/* 
  4 * FUNCTIONS TO GENERATE VARIANCE-COVARIANCE MATRIX from an S3 'phylo' object from the R-package ape (Paradis E, J Claude, and K Strimmer 2004 [BIOINFORMATICS 20:289-290])
  5 * author: Jonathan M Eastman 01.11.2011
  6 */
  7
  8/* INTERNAL C++ */
  9void sortedges(std::vector<double> const &unsortededges,
 10			   std::vector<double> &edges,
 11			   std::vector<int> const &des)
 12{
 13	std::vector<int>::size_type i,j;
 14	//int i, j;
 15	for(i=0; i<edges.size(); i++) {
 16		for(j=0; j<des.size(); j++){
 17			if(des.at(j)==(signed)i+1)
 18			{
 19				edges.at(i)=unsortededges.at(j);							 
 20			}
 21		}
 22	}
 23}
 24
 25/* INTERNAL C++ */
 26void gatherdescendants(int const &node,
 27					   int const &root,
 28					   int const &endofclade,
 29					   std::vector<int> &TIPS,
 30					   std::vector<int> const &anc, 
 31					   std::vector<int> const &des)
 32{
 33	/* 
 34	 * node: current internal node
 35	 * root: most internal node
 36	 * endofclade: number of rows in the phylogenetic edge matrix (of the 'phylo' object)
 37	 * TIPS: vector in which to store all (tip) descendants of node (by node label)
 38	 * anc: first column of 'phylo' edge matrix (e.g., phy$edge[,1])
 39	 * des: second column of 'phylo' edge matrix (e.g., phy$edge[,2])
 40	 */
 41	
 42	int i, eoc=endofclade, r=root, n=node;
 43
 44	for(i=0; i<eoc; i++) {
 45		if(anc.at(i) == n)
 46		{
 47			if(des.at(i)>r)
 48			{
 49				gatherdescendants(des.at(i), r, eoc, TIPS, anc, des);
 50			}
 51			else 
 52			{
 53				TIPS.push_back(des.at(i));
 54	
 55			}
 56		}
 57	}
 58}
 59	
 60					   
 61/* INTERNAL C++ */
 62void descend_vcv(int const &node,
 63				 double const &edge,
 64				 int const &root,
 65				 int const &endofclade,
 66				 std::vector<int> const &anc, 
 67				 std::vector<int> const &des, 
 68				 std::vector<double> &vcv)
 69{
 70	/* 
 71	 * node: internal node of which descendants are recorded
 72	 * edge: length to add to covariances of all descendants where 'node' is internal
 73	 * root: node to initiate vcv update
 74	 * endofclade: rows in the phylogenetic edge matrix (of the 'phylo' object)
 75	 * anc: first column of phylo edge matrix (e.g., phy$edge[,1]
 76	 * des: second column of phylo edge matrix (e.g., phy$edge[,2]
 77	 * vcv: the current variance-covariance matrix
 78	 * TIPS: an empty vector in which to store all (tip) descendants of node
 79	 */
 80	
 81	int n;
 82	std::vector<int>::size_type a,b,j;
 83	n=root-1;
 84	std::vector<int> TIPS;	
 85	TIPS.reserve(root-1);
 86	
 87	gatherdescendants(node, root, endofclade, TIPS, anc, des);
 88		
 89	/* add 'edge' length of 'node' to covariances V[a,b] where a and b are tips descended from 'node' and a!=b [OFFDIAGONAL ELEMENTS in V] */
 90	for(a=0; a<TIPS.size(); a++) {
 91		for(b=a; b<TIPS.size(); b++) {
 92			if(a!=b)
 93			{
 94				vcv[ TIPS.at(a)-1 + n*( TIPS.at(b)-1 ) ] = vcv[ TIPS.at(b)-1 + n*(TIPS.at(a)-1 ) ] += edge;
 95			}
 96		}
 97	}
 98	
 99	/* add 'edge' length of 'node' to variances of V[a,b] where a and b are tips descended from 'node' and a==b [DIAGONAL ELEMENTS in V]*/ 
100	for(j=0; j<TIPS.size(); j++) {
101		vcv[ TIPS.at(j)-1 + n*( TIPS.at(j)-1 ) ] += edge;
102	}
103	
104	std::vector<int>().swap(TIPS);
105}
106
107
108/* INTERNAL C++ */
109void vcv_internal (int const &maxnode,
110				   int const &root,
111				   int const &endofclade,
112				   std::vector<int> const &anc, 
113				   std::vector<int> const &des, 
114				   std::vector<double> const &edges, 
115				   std::vector<double> &vcv)
116{
117	/* 
118	 * maxnode: largest internal node
119	 * root: most internal node
120	 * endofclade: number of rows in the phylogenetic edge matrix (of the 'phylo' object)
121	 * anc: first column of phylo edge matrix (e.g., phy$edge[,1]
122	 * des: second column of phylo edge matrix (e.g., phy$edge[,2]
123	 * edges: sorted branch lengths (by node label), including the root (ntips+1)
124	 * vcv: the current state of the variance-covariance matrix 
125	 */
126	
127	int n, nn, r=root, mx=maxnode, eoc=endofclade;	
128	nn=r-1;
129	
130	/* cycle over internal nodes within the tree; tree stem does not contribute to VCV so is skipped */
131	for(n=r+1; n<=mx; n++) {
132		
133		/* find descendants of n; add edge length of n to variances for tips and covariances among tips */	
134		descend_vcv(n, edges.at(n-1), r, eoc, anc, des, vcv);
135		
136	}
137	
138	/* add leaves to lengths in V[n,n] where n is a tip [DIAGONAL ELEMENTS in V]*/
139	for(n=1; n<r; n++) {
140		vcv[ n-1 + nn*( n-1 ) ] += edges.at(n-1);
141	}
142}
143
144
145/* C++ | R INTERFACE; main function */
146RcppExport SEXP vmat (SEXP tree) 
147{
148	/* 
149	 * tree: a list of elements 
150		* ROOT: most internal node
151		* MAXNODE: least internal internal node (and largest valued in edge matrix)
152		* ENDOFCLADE: rows in edge matrix
153		* ANC: first column of 'phylo' edge matrix (e.g., phy$edge[,1])
154		* DES: second column of 'phylo' edge matrix (e.g., phy$edge[,2])
155		* EDGES: edge lengths, sorted by node label and including the root (at position Ntip(phy)+1)
156		* VCV: the current state of the variance-covariance matrix, initialized with 0s in R
157	 */
158		
159	try {
160		std::vector<int>::size_type i;
161		
162		/* call in parameters associated with 'phylo' object */
163		Rcpp::List phylo(tree);
164
165		int root = (int) Rcpp::as<int>(phylo["ROOT"]);
166		int maxnode =  Rcpp::as<int>(phylo["MAXNODE"]);
167		int endofclade =  Rcpp::as<int>(phylo["ENDOFCLADE"]);
168		std::vector<int> anc=phylo["ANC"];
169		std::vector<int> des=phylo["DES"];
170		std::vector<double> unsortededges=phylo["EDGES"];
171		std::vector<double> edges=phylo["EDGES"];
172		std::vector<double> V=phylo["VCV"];
173		
174		/* initialize edges */
175		for(i=0; i<edges.size(); i++) {
176			edges.at(i)=0;
177		}		
178		
179		/* sort edges by node label */
180		sortedges(unsortededges, edges, des);
181		
182		/* call to primary function that updates VCV matrix */
183		vcv_internal(maxnode, root, endofclade, anc, des, edges, V);
184		
185		/* PREPARE OUTPUT FOR R */
186		return Rcpp::List::create(Rcpp::Named("VCV",V));
187
188
189    } catch( std::exception &ex ) {		
190		forward_exception_to_r( ex );
191    } catch(...) { 
192		::Rf_error( "C++ exception: unknown reason" ); 
193    }
194    return R_NilValue; 
195}