/auteur/src/oumat.cpp

 1#include <Rcpp.h>
 2#include <vector>
 3
 4/* internal C++ */
 5void process_OU_covariance ( std::vector<double> const &VCV, 
 6							 double const &alphasq, 
 7							 double const &sigmasq, 
 8							 int const &ntip,
 9							 std::vector<double> &V) 
10{	
11	unsigned int i, j, nt=ntip;
12	double sij, ti, tj, elti, eltj, W;
13	
14	for (i = 0; i < nt; i++) {
15		for (j = 0; j <= i; j++) {
16			ti = VCV[i+i*nt];
17			sij = VCV[i+j*nt];
18			tj = VCV[j+j*nt];
19			elti = exp(-alphasq*(ti-sij));
20			eltj = exp(-alphasq*(tj-sij));
21			W = elti*sigmasq*eltj/(alphasq+alphasq);
22			V[i+nt*j] += W;
23			if (j != i) 
24			{
25				V[j+nt*i] += W;
26			}
27		}
28	}
29}
30
31
32/* C++ | R INTERFACE; main function */
33RcppExport SEXP oumat (
34						   SEXP TIMES, 
35						   SEXP ALPHA, 
36						   SEXP SIGMA,
37						   SEXP TIPS) 
38{
39	/* 
40	 * TIMES: array form of VCV matrix (under Brownian motion): split times for pairs of species 
41	 * ALPHA: herein, alpha is expected given as alpha^2
42	 * SIGMA: herein, sigma is expected given as sigma^2
43	 * TIPS: number of species
44	 */
45	
46	try {
47		Rcpp::List phylo(TIMES);
48		std::vector<double> VCV =  phylo["vmat"];
49		double alphasq = Rcpp::as<double> (ALPHA);
50		double sigmasq = Rcpp::as<double> (SIGMA);
51		unsigned int ntip = Rcpp::as<int> (TIPS);
52		
53		std::vector<double> V;
54		V.assign(ntip*ntip,0);
55		
56		process_OU_covariance(VCV, alphasq, sigmasq, ntip, V);
57		
58		/* PREPARE OUTPUT FOR R */
59		return Rcpp::List::create(Rcpp::Named("oumat", V));
60		
61	} catch( std::exception &ex ) {		
62		forward_exception_to_r( ex );
63    } catch(...) { 
64		::Rf_error( "c++ exception (unknown reason)" ); 
65    }
66    return R_NilValue; 
67}
68