Merge pull request #233 from wush978/master

Resolve #233

Author: eddelbuettel

inst/include/Rcpp/sugar/functions/var.h

@@ -3,6 +3,7 @@
 // var.h: Rcpp R/C++ interface class library -- var
 //
 // Copyright (C) 2011 Dirk Eddelbuettel and Romain Francois
+// Copyright (C) 2015 Wush Wu
 //
 // This file is part of Rcpp.
 //
@@ -26,31 +27,65 @@ namespace Rcpp{
 namespace sugar{
 
 template <int RTYPE, bool NA, typename T>
-class Var : public Lazy< typename Rcpp::traits::storage_type<RTYPE>::type , Var<RTYPE,NA,T> > {
+class Var : public Lazy< double , Var<RTYPE,NA,T> > {
 public:
-	typedef typename Rcpp::VectorBase<RTYPE,NA,T> VEC_TYPE ;
-	typedef typename Rcpp::traits::storage_type<RTYPE>::type STORAGE ;
+    typedef typename Rcpp::VectorBase<RTYPE,NA,T> VEC_TYPE ;
 
-	Var( const VEC_TYPE& object_ ) : object(object_){}
+    Var( const VEC_TYPE& object_ ) : object(object_){}
 
-	STORAGE get() const{
-	    STORAGE m = mean(object).get() ;
-		Minus_Vector_Primitive<RTYPE,NA,T> mm( object, m) ;
-		STORAGE ssq  = sum( pow(mm,2.0) ).get() ;
-	    return ssq / (object.size() - 1 ) ;
-	}
+    double get() const{
+        double
+            m = mean(object).get(),
+            ssq = sum( pow(object, 2.0) ).get();
+        return (ssq - m * m * object.size()) / (object.size() - 1);
+    }
 
 private:
-	const VEC_TYPE& object ;
+    const VEC_TYPE& object ;
+} ;
+
+template <bool NA, typename T>
+class Var<CPLXSXP,NA,T> : public Lazy< double , Var<CPLXSXP,NA,T> > {
+public:
+    typedef typename Rcpp::VectorBase<CPLXSXP,NA,T> VEC_TYPE ;
+
+    Var( const VEC_TYPE& object_ ) : object(object_){}
+
+    double get() const{
+        double sq = 0, ssq = 0;
+        for(int i = 0;i < object.size();i++) {
+            Rcomplex z = object[i];
+            sq += z.r;
+            ssq += z.r * z.r;
+        }
+        return (ssq - sq * sq / object.size()) / (object.size() - 1);
+    }
+
+private:
+    const VEC_TYPE& object ;
 } ;
 
 } // sugar
 
 template <bool NA, typename T>
 inline sugar::Var<REALSXP,NA,T> var( const VectorBase<REALSXP,NA,T>& t){
-	return sugar::Var<REALSXP,NA,T>( t ) ;
+    return sugar::Var<REALSXP,NA,T>( t ) ;
 }
 
+template <bool NA, typename T>
+inline sugar::Var<INTSXP,NA,T> var( const VectorBase<INTSXP,NA,T>& t){
+    return sugar::Var<INTSXP,NA,T>( t ) ;
+}
+
+template <bool NA, typename T>
+inline sugar::Var<LGLSXP,NA,T> var( const VectorBase<LGLSXP,NA,T>& t){
+    return sugar::Var<LGLSXP,NA,T>( t ) ;
+}
+
+template <bool NA, typename T>
+inline sugar::Var<CPLXSXP,NA,T> var( const VectorBase<CPLXSXP,NA,T>& t){
+    return sugar::Var<CPLXSXP,NA,T>( t ) ;
+}
 
 } // Rcpp
 #endif

inst/unitTests/runit.sugar.var.R

@@ -0,0 +1,36 @@
+#!/usr/bin/r -t
+#                     -*- mode: R; ess-indent-level: 4; indent-tabs-mode: nil; -*-
+#
+# Copyright (C) 2015 Wush Wu
+#
+# This file is part of Rcpp.
+#
+# Rcpp is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 2 of the License, or
+# (at your option) any later version.
+#
+# Rcpp is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Rcpp.  If not, see <http://www.gnu.org/licenses/>.
+
+.runThisTest <- Sys.getenv("RunAllRcppTests") == "yes"
+
+if (.runThisTest) {
+
+    test.Sugar.var <- function() {
+        f1 <- Rcpp::cppFunction('double myVar(NumericVector x) { return(var(x)); }')
+        f2 <- Rcpp::cppFunction('double myVar(IntegerVector x) { return(var(x)); }')
+        f3 <- Rcpp::cppFunction('double myVar(ComplexVector x) { return(var(x)); }')
+        f4 <- Rcpp::cppFunction('double myVar(LogicalVector x) { return(var(x)); }')
+        checkEquals(f1((1:10) * 1.1), var((1:10) * 1.1))
+        checkEquals(f2(1:10), var(1:10))
+        checkEquals(f3(1:10 + (1 + 1i)), var(1:10 + (1 + 1i)))
+        checkEquals(f4(c(T, F, T, F, T)), var(c(T, F, T, F, T)))
+    }
+
+}