/working/auteur.extended/src/ouweights.cpp

  1#include <Rcpp.h>
  2#include <vector>
  3
  4/* internal C++: computes the contribution of each ancestral edge for a species to the weights matrix */
  5void process_epochs ( std::vector<double> const &epochs, 
  6					  double const &alphasq, 
  7					  std::vector<double> &y ) 
  8{
  9	int i, np = epochs.size();
 10	double t;
 11	
 12	for (i = 0; i < np; i++) {
 13		t = epochs[0]-epochs[i];
 14		y.push_back( exp(-alphasq*t) );
 15	}
 16	for (i = 0; i < np-1; i++) {
 17		y[i] -= y[i+1];
 18	}
 19}
 20
 21/* internal C++: builds a binary array for each species, signifying in which selective regime the ancestor of the species is bounded */
 22void process_beta ( std::vector<int> &beta,
 23				    std::vector<int> const &ancestors,
 24				    std::vector<double> const &regimes,
 25				    std::vector<double> const &optima )
 26{
 27	int np, r, k, nreg=optima.size();
 28
 29	np=ancestors.size();
 30	beta.assign(np*nreg, 0);
 31	
 32	for(r=0; r<np; r++) {
 33		for(k=0; k<nreg; k++){
 34			if(regimes.at( ancestors.at(r)-1 )==optima.at(k))
 35			{
 36				beta.at(r+k*np)=1;
 37			}
 38		}
 39	}
 40}
 41
 42
 43/* C++ | R INTERFACE; main function */
 44RcppExport SEXP ouweights ( SEXP PHYLO ) 
 45{
 46	/* 
 47	 * PHYLO: 
 48	 *	- epochs: list of branching times for each segment from tip to root [T=t1>t2>...>troot=0] for each species
 49	 *	- ancestors: list of nodes from tip to root for each species 
 50	 *	- regimes: array of optima associated with each edge; order from node 1 to node 2N-2 
 51	 *	- optima: array of unique selective regimes across the tree
 52	 *	- alpha: constraint factor for the Ornstein-Uhlenbeck process; herein, alpha is expected given as alpha^2
 53	 *  - tips: number of species in tree
 54	 */
 55			
 56	try {
 57		Rcpp::List phylo(PHYLO);
 58		
 59		Rcpp::List epochs = phylo["epochs"];
 60		Rcpp::List ancestors = phylo["ancestors"];
 61		std::vector<double> regimes =  phylo["regimes"];
 62		std::vector<double> optima = phylo["optima"];
 63		double alphasq = Rcpp::as<double> ( phylo["alpha"] );
 64		int ntip = Rcpp::as<int> ( phylo["tips"] );
 65		int nreg = optima.size();
 66		
 67		int t, r, q;
 68		int segments;
 69		std::vector<double> y;
 70		std::vector<double> W;
 71		W.assign(ntip*nreg,0);
 72		std::vector<int> tbeta;
 73
 74		for(t=0; t<ntip; t++) {
 75			std::vector<int> tancestors=ancestors[t];
 76			process_beta(tbeta, tancestors, regimes, optima); 
 77			std::vector<double> tepochs=epochs[t];
 78			segments=tepochs.size();
 79			y.reserve(segments);
 80			process_epochs(tepochs, alphasq, y);
 81			for(r=0; r<nreg; r++) {
 82				for(q=0; q<segments; q++) {
 83					W[t+r*ntip] += y[q]*tbeta[q+r*segments];
 84				}
 85			}
 86			std::vector<double>().swap(y);
 87			std::vector<double>().swap(tepochs);
 88			std::vector<int>().swap(tancestors);
 89			std::vector<int>().swap(tbeta);
 90		}
 91		
 92		/* PREPARE OUTPUT FOR R */
 93		return Rcpp::List::create(Rcpp::Named("weights", W));
 94
 95	} catch( std::exception &ex ) {		
 96		forward_exception_to_r( ex );
 97    } catch(...) { 
 98		::Rf_error( "c++ exception (unknown reason)" ); 
 99    }
100    return R_NilValue; 
101}
102