@@ -26,21 +26,20 @@ namespace Rcpp{
26
26
namespace sugar {
27
27
28
28
template <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> > {
30
30
public:
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;
33
34
34
- Mean ( const VEC_TYPE& object_ ) : object(object_){}
35
+ Mean (const VEC_TYPE& object_) : object(object_) {}
35
36
36
37
STORAGE get () const {
37
- // return sum(object).get() / object.size() ;
38
- NumericVector input = object;
39
-
38
+ VECTOR input = object;
40
39
int n = input.size (); // double pass (as in summary.c)
41
40
long double s = std::accumulate (input.begin (), input.end (), 0 .0L );
42
41
s /= n;
43
- if (R_FINITE ((double )s)) {
42
+ if (R_FINITE ((double )s)) {
44
43
long double t = 0.0 ;
45
44
for (int i = 0 ; i < n; i++) {
46
45
t += input[i] - s;
@@ -51,16 +50,63 @@ class Mean : public Lazy< typename Rcpp::traits::storage_type<RTYPE>::type , Mea
51
50
}
52
51
private:
53
52
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
+ };
55
91
56
92
} // sugar
57
93
58
94
template <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);
61
102
}
62
103
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
+ }
63
108
64
109
} // Rcpp
65
110
#endif
66
111
112
+
0 commit comments