@@ -26,21 +26,20 @@ namespace Rcpp{
2626namespace sugar {
2727
2828template <int RTYPE, bool NA, typename T>
29- class Mean : public Lazy < typename Rcpp::traits::storage_type<RTYPE>::type , Mean<RTYPE,NA,T> > {
29+ class Mean : public Lazy <typename Rcpp::traits::storage_type<RTYPE>::type, Mean<RTYPE,NA,T> > {
3030public:
31- typedef typename Rcpp::VectorBase<RTYPE,NA,T> VEC_TYPE ;
32- typedef typename Rcpp::traits::storage_type<RTYPE>::type STORAGE ;
31+ typedef typename Rcpp::VectorBase<RTYPE,NA,T> VEC_TYPE;
32+ typedef typename Rcpp::traits::storage_type<RTYPE>::type STORAGE;
33+ typedef Rcpp::Vector<RTYPE> VECTOR;
3334
34- Mean ( const VEC_TYPE& object_ ) : object(object_){}
35+ Mean (const VEC_TYPE& object_) : object(object_) {}
3536
3637 STORAGE get () const {
37- // return sum(object).get() / object.size() ;
38- NumericVector input = object;
39-
38+ VECTOR input = object;
4039 int n = input.size (); // double pass (as in summary.c)
4140 long double s = std::accumulate (input.begin (), input.end (), 0 .0L );
4241 s /= n;
43- if (R_FINITE ((double )s)) {
42+ if (R_FINITE ((double )s)) {
4443 long double t = 0.0 ;
4544 for (int i = 0 ; i < n; i++) {
4645 t += input[i] - s;
@@ -51,16 +50,63 @@ class Mean : public Lazy< typename Rcpp::traits::storage_type<RTYPE>::type , Mea
5150 }
5251private:
5352 const VEC_TYPE& object ;
54- } ;
53+ };
54+
55+ template <bool NA, typename T>
56+ class Mean <CPLXSXP,NA,T> : public Lazy<Rcomplex, Mean<CPLXSXP,NA,T> > {
57+ public:
58+ typedef typename Rcpp::VectorBase<CPLXSXP,NA,T> VEC_TYPE;
59+
60+ Mean (const VEC_TYPE& object_) : object(object_) {}
61+
62+ Rcomplex get () const {
63+ ComplexVector input = object;
64+ int n = input.size (); // double pass (as in summary.c)
65+ long double s = 0.0 , si = 0.0 ;
66+ for (int i=0 ; i<n; i++) {
67+ Rcomplex z = input[i];
68+ s += z.r ;
69+ si += z.i ;
70+ }
71+ s /= n;
72+ si /= n;
73+ if (R_FINITE ((double )s) && R_FINITE ((double )si)) {
74+ long double t = 0.0 , ti = 0.0 ;
75+ for (int i = 0 ; i < n; i++) {
76+ Rcomplex z = input[i];
77+ t += z.r - s;
78+ ti += z.i - si;
79+ }
80+ s += t/n;
81+ si += ti/n;
82+ }
83+ Rcomplex z;
84+ z.r = s;
85+ z.i = si;
86+ return z;
87+ }
88+ private:
89+ const VEC_TYPE& object ;
90+ };
5591
5692} // sugar
5793
5894template <bool NA, typename T>
59- inline sugar::Mean<REALSXP,NA,T> mean ( const VectorBase<REALSXP,NA,T>& t){
60- return sugar::Mean<REALSXP,NA,T>( t ) ;
95+ inline sugar::Mean<REALSXP,NA,T> mean (const VectorBase<REALSXP,NA,T>& t) {
96+ return sugar::Mean<REALSXP,NA,T>(t);
97+ }
98+
99+ template <bool NA, typename T>
100+ inline sugar::Mean<INTSXP,NA,T> mean (const VectorBase<INTSXP,NA,T>& t) {
101+ return sugar::Mean<INTSXP,NA,T>(t);
61102}
62103
104+ template <bool NA, typename T>
105+ inline sugar::Mean<CPLXSXP,NA,T> mean (const VectorBase<CPLXSXP,NA,T>& t) {
106+ return sugar::Mean<CPLXSXP,NA,T>(t);
107+ }
63108
64109} // Rcpp
65110#endif
66111
112+
0 commit comments