Add mean and mean's unit test.

src/Function2D.cpp

@@ -156,6 +156,7 @@ Matrix Aurora::min(const Matrix &aMatrix, FunctionDirection direction) {
             return Matrix::New(ret,aMatrix.getDimSize(0),1);
         }
         case Column:
+        default:
         {
             Eigen::Map<Eigen::MatrixXd> srcMatrix(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
             double * ret = malloc(aMatrix.getDimSize(0));
@@ -266,28 +267,139 @@ Matrix Aurora::sum(const Matrix &aMatrix, FunctionDirection direction) {
     {
         case All:
         {
+            Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData(),aMatrix.getDataSize());
             double * ret = malloc(1);
-            ret[0] = cblas_dasum(aMatrix.getDataSize(),aMatrix.getData(),1);
+            ret[0] = srcV.array().sum();
             return Matrix::New(ret,1);
         }
         case Row:
         {
+            Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
             double * ret = malloc(aMatrix.getDimSize(0));
-            for (int i = 0; i < aMatrix.getDimSize(0); ++i) {
+            Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(0));
-                ret[i] = cblas_dasum(aMatrix.getDimSize(1), aMatrix.getData() + i,
+            retV = srcM.rowwise().sum();
-                                     aMatrix.getDimSize(0));
-            }
             return Matrix::New(ret,aMatrix.getDimSize(0),1);
         }
         case Column:
         default:
         {
-            double * ret = malloc(aMatrix.getDimSize(0));
+            Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
-            for (int i = 0; i < aMatrix.getDimSize(1); ++i) {
+            double * ret = malloc(aMatrix.getDimSize(1));
-                ret[i] = cblas_dasum(aMatrix.getDimSize(0), aMatrix.getData()+aMatrix.getDimSize(0)*i,
+            Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(0));
-                                     1);
+            retV = srcM.colwise().sum();
-            }
             return Matrix::New(ret,1,aMatrix.getDimSize(1));
         }
     }
 }
+
+void excludeNan(double * aInput, int aInputSize,double* aOutput,int& aOutputSize){
+    aOutputSize = 0;
+    for (int i = 0; i < aInputSize; ++i) {
+        if (std::isnan(aInput[i])) continue;
+        aOutput[aOutputSize]=aInput[i];
+        aOutputSize++;
+    }
+}
+
+Matrix Aurora::mean(const Matrix &aMatrix, FunctionDirection direction, bool aIncludeNan) {
+    if (aMatrix.getDimSize(2)>1 || aMatrix.isComplex()) {
+        std::cerr
+                << (aMatrix.getDimSize(2) > 1 ? "sum() not support 3D data!" : "sum() not support complex value type!")
+                << std::endl;
+        return Matrix();
+    }
+    if (aIncludeNan){
+        switch (direction)
+        {
+            case All:
+            {
+                Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData(),aMatrix.getDataSize());
+                double * ret = malloc(1);
+                ret[0] = srcV.array().mean();
+                return Matrix::New(ret,1);
+            }
+            case Row:
+            {
+                Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                double * ret = malloc(aMatrix.getDimSize(0));
+                Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(0));
+                retV = srcM.rowwise().mean();
+                return Matrix::New(ret,aMatrix.getDimSize(0),1);
+            }
+            case Column:
+            default:
+            {
+                Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                double * ret = malloc(aMatrix.getDimSize(1));
+                Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(0));
+                retV = srcM.colwise().mean();
+                return Matrix::New(ret,1,aMatrix.getDimSize(1));
+            }
+        }
+    }
+    else{
+        switch (direction)
+        {
+            case All:
+            {
+                Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData(),aMatrix.getDataSize());
+                double * retVd = malloc(aMatrix.getDataSize());
+                Eigen::Map<Eigen::VectorXd> retV(retVd,aMatrix.getDataSize());
+                Eigen::VectorXd ones = Eigen::VectorXd(aMatrix.getDataSize());
+                ones.setConstant(1.0);
+                retV = srcV.array().isNaN().select(0.0,ones);
+                int count = retV.sum();
+                if (count == 0){
+                    free(retVd);
+                    double *ret = malloc(1);
+                    ret[0]=0;
+                    return Matrix::New(ret,1);
+                }
+                else {
+                    double *ret = malloc(1);
+                    retV = srcV.array().isNaN().select(0.0,srcV);
+                    ret[0] = retV.sum() / count;
+                    free(retVd);
+                    return Matrix::New(ret, 1);
+                }
+            }
+            case Row:
+            {
+                Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                double * retMd = malloc(aMatrix.getDataSize());
+                Eigen::Map<Eigen::MatrixXd> retM(retMd,aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                Eigen::MatrixXd zeros = Eigen::MatrixXd(aMatrix.getDimSize(0), aMatrix.getDimSize(1));
+                zeros.setConstant(0.0);
+                retM = srcM.array().isNaN().select(zeros,1.0);
+                Eigen::VectorXd countM = retM.rowwise().sum();
+                countM = (countM.array()==0.0).select(1.0,countM);
+                retM = srcM.array().isNaN().select(0.0,srcM);
+                double * ret = malloc(aMatrix.getDimSize(0));
+                Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(0));
+                retV = retM.rowwise().sum();
+                retV =retV.array()/countM.array();
+                free(retMd);
+                return Matrix::New(ret,aMatrix.getDimSize(0),1);
+            }
+            case Column:
+            default:
+            {
+                Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                double * retMd = malloc(aMatrix.getDataSize());
+                Eigen::Map<Eigen::MatrixXd> retM(retMd,aMatrix.getDimSize(0),aMatrix.getDimSize(1));
+                Eigen::MatrixXd zeros = Eigen::MatrixXd(aMatrix.getDimSize(0), aMatrix.getDimSize(1));
+                zeros.setConstant(0.0);
+                retM = srcM.array().isNaN().select(zeros,1.0);
+                Eigen::VectorXd countM = retM.colwise().sum();
+                countM = (countM.array()==0).select(1.0,countM);
+                retM = srcM.array().isNaN().select(0.0,srcM);
+                double * ret = malloc(aMatrix.getDimSize(1));
+                Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(1));
+                retV = retM.colwise().sum();
+                retV = retV.array()/countM.array();
+                free(retMd);
+                return Matrix::New(ret,1,aMatrix.getDimSize(1));
+            }
+        }
+    }
+}

src/Function2D.h

@@ -51,6 +51,15 @@ namespace Aurora {
      */
     Matrix sum(const Matrix& aMatrix,FunctionDirection direction = Column);
 
+    /**
+     * 求矩阵平均值，可按行、列、单元, 目前不支持三维，不支持复数
+     * @param aMatrix 矩阵
+     * @param direction 方向，Column, Row, All
+     * @param aIncludeNan 是否包含nan
+     * @return
+     */
+    Matrix mean(const Matrix& aMatrix,FunctionDirection direction = Column, bool aIncludeNan = true);
+
 };
 
 

test/Function2D_Test.cpp

@@ -74,7 +74,7 @@ TEST_F(Function2D_Test, std){
 
 TEST_F(Function2D_Test, min) {
     double *dataA = new double[3]{1, 2, 3};
-    double *dataB = new double[9]{2, 3, 3, 2, 2, 1, 3, 3, 3};
+    double *dataB = new double[9]{2, 3, 3, 2, 2, -1, 3, 3, 3};
     double *dataC = new double[1]{1.5};
     auto A = Aurora::Matrix::fromRawData(dataA, 3, 1);
     auto B = Aurora::Matrix::fromRawData(dataB, 3, 3);
@@ -84,18 +84,18 @@ TEST_F(Function2D_Test, min) {
     EXPECT_EQ(1, ret.getDimSize(0));
     EXPECT_EQ(3, ret.getDimSize(1));
     EXPECT_DOUBLE_EQ(2, ret.getData()[0]);
-    EXPECT_DOUBLE_EQ(1, ret.getData()[1]);
+    EXPECT_DOUBLE_EQ(-1, ret.getData()[1]);
     EXPECT_DOUBLE_EQ(3, ret.getData()[2]);
     ret = Aurora::min(B, Aurora::All);
     EXPECT_DOUBLE_EQ(1, ret.getDataSize());
-    EXPECT_DOUBLE_EQ(1, ret.getData()[0]);
+    EXPECT_DOUBLE_EQ(-1, ret.getData()[0]);
     ret = Aurora::min(B, Aurora::Row);
     EXPECT_DOUBLE_EQ(3, ret.getDataSize());
     EXPECT_EQ(3, ret.getDimSize(0));
     EXPECT_EQ(1, ret.getDimSize(1));
     EXPECT_DOUBLE_EQ(2, ret.getData()[0]);
     EXPECT_DOUBLE_EQ(2, ret.getData()[1]);
-    EXPECT_DOUBLE_EQ(1, ret.getData()[2]);
+    EXPECT_DOUBLE_EQ(-1, ret.getData()[2]);
     ret = Aurora::min(A, C);
     EXPECT_DOUBLE_EQ(3, ret.getDataSize());
     EXPECT_DOUBLE_EQ(1, ret.getData()[0]);
@@ -132,24 +132,104 @@ TEST_F(Function2D_Test, max) {
 }
 
 TEST_F(Function2D_Test, sum) {
-    double *dataB = new double[9]{2, 3, 3, 2, 2, 1, 3, 3, 3};
+    double *dataB = new double[9]{2, 3, 3, 2, 2, 1, 3, 3, -3};
     auto B = Aurora::Matrix::fromRawData(dataB, 3, 3);
     Aurora::Matrix ret = Aurora::sum(B);
     EXPECT_EQ(1, ret.getDimSize(0));
     EXPECT_EQ(3, ret.getDimSize(1));
     EXPECT_DOUBLE_EQ(8, ret.getData()[0]);
     EXPECT_DOUBLE_EQ(5, ret.getData()[1]);
-    EXPECT_DOUBLE_EQ(9, ret.getData()[2]);
+    EXPECT_DOUBLE_EQ(3, ret.getData()[2]);
     ret = Aurora::sum(B, Aurora::All);
     EXPECT_DOUBLE_EQ(1, ret.getDataSize());
-    EXPECT_DOUBLE_EQ(22, ret.getData()[0]);
+    EXPECT_DOUBLE_EQ(16, ret.getData()[0]);
     ret = Aurora::sum(B, Aurora::Row);
     EXPECT_DOUBLE_EQ(3, ret.getDataSize());
     EXPECT_EQ(3, ret.getDimSize(0));
     EXPECT_EQ(1, ret.getDimSize(1));
     EXPECT_DOUBLE_EQ(7, ret.getData()[0]);
     EXPECT_DOUBLE_EQ(8, ret.getData()[1]);
-    EXPECT_DOUBLE_EQ(7, ret.getData()[2]);
+    EXPECT_DOUBLE_EQ(1, ret.getData()[2]);
+}
+
+TEST_F(Function2D_Test, mean) {
+    {
+        double *dataB = new double[16]{1.1, 2.6, 3.8, 6.2,
+                                       4.3, 5.7, 6.9, 10.6,
+                                       7.1, 8.3, 9.7, 11.2,
+                                       17.8, 13.3, 26.5, -7.7};
+        auto B = Aurora::Matrix::fromRawData(dataB, 4, 4);
+        auto r = Aurora::mean(B, Aurora::All);
+        EXPECT_EQ(1, r.getDimSize(0));
+        EXPECT_EQ(1, r.getDimSize(1));
+        EXPECT_DOUBLE_EQ(7.9625, fourDecimalRound(r.getData()[0]));
+
+        r = Aurora::mean(B);
+        EXPECT_EQ(1, r.getDimSize(0));
+        EXPECT_EQ(4, r.getDimSize(1));
+        EXPECT_DOUBLE_EQ(3.4250, fourDecimalRound(r.getData()[0]));
+        EXPECT_DOUBLE_EQ(6.8750, fourDecimalRound(r.getData()[1]));
+        EXPECT_DOUBLE_EQ(9.0750, fourDecimalRound(r.getData()[2]));
+        EXPECT_DOUBLE_EQ(12.4750, fourDecimalRound(r.getData()[3]));
+
+        r = Aurora::mean(B, Aurora::Row);
+        EXPECT_EQ(4, r.getDimSize(0));
+        EXPECT_EQ(1, r.getDimSize(1));
+        EXPECT_DOUBLE_EQ(7.5750, fourDecimalRound(r.getData()[0]));
+        EXPECT_DOUBLE_EQ(7.4750, fourDecimalRound(r.getData()[1]));
+        EXPECT_DOUBLE_EQ(11.7250, fourDecimalRound(r.getData()[2]));
+        EXPECT_DOUBLE_EQ(5.0750, fourDecimalRound(r.getData()[3]));
+    }
+    //with nan
+    {
+        double tnan = std::nan("");
+        double *dataB = new double[16]{1.1, 2.6, 3.8, 6.2,
+                                       4.3, 5.7, 6.9, 10.6,
+                                       7.1, 8.3, 9.7, 11.2,
+                                       17.8, 13.3,tnan , -7.7};
+        auto B = Aurora::Matrix::fromRawData(dataB, 4, 4);
+        auto r = Aurora::mean(B, Aurora::All);
+        EXPECT_EQ(1, r.getDimSize(0));
+        EXPECT_EQ(1, r.getDimSize(1));
+        EXPECT_TRUE(std::isnan(r.getData()[0]));
+
+        r = Aurora::mean(B);
+        EXPECT_EQ(1, r.getDimSize(0));
+        EXPECT_EQ(4, r.getDimSize(1));
+        EXPECT_DOUBLE_EQ(3.4250, fourDecimalRound(r.getData()[0]));
+        EXPECT_DOUBLE_EQ(6.8750, fourDecimalRound(r.getData()[1]));
+        EXPECT_DOUBLE_EQ(9.0750, fourDecimalRound(r.getData()[2]));
+        EXPECT_TRUE(std::isnan(r.getData()[3]));
+
+        r = Aurora::mean(B, Aurora::Row);
+        EXPECT_EQ(4, r.getDimSize(0));
+        EXPECT_EQ(1, r.getDimSize(1));
+        EXPECT_DOUBLE_EQ(7.5750, fourDecimalRound(r.getData()[0]));
+        EXPECT_DOUBLE_EQ(7.4750, fourDecimalRound(r.getData()[1]));
+        EXPECT_TRUE(std::isnan(r.getData()[2]));
+        EXPECT_DOUBLE_EQ(5.0750, fourDecimalRound(r.getData()[3]));
+
+        r = Aurora::mean(B, Aurora::All,false);
+        EXPECT_EQ(1, r.getDimSize(0));
+        EXPECT_EQ(1, r.getDimSize(1));
+        EXPECT_DOUBLE_EQ(6.7267, fourDecimalRound(r.getData()[0]));
+
+        r = Aurora::mean(B, Aurora::Column, false);
+        EXPECT_EQ(1, r.getDimSize(0));
+        EXPECT_EQ(4, r.getDimSize(1));
+        EXPECT_DOUBLE_EQ(3.4250, fourDecimalRound(r.getData()[0]));
+        EXPECT_DOUBLE_EQ(6.8750, fourDecimalRound(r.getData()[1]));
+        EXPECT_DOUBLE_EQ(9.0750, fourDecimalRound(r.getData()[2]));
+        EXPECT_DOUBLE_EQ(7.8, fourDecimalRound(r.getData()[3]));
+
+        r = Aurora::mean(B, Aurora::Row, false);
+        EXPECT_EQ(4, r.getDimSize(0));
+        EXPECT_EQ(1, r.getDimSize(1));
+        EXPECT_DOUBLE_EQ(7.5750, fourDecimalRound(r.getData()[0]));
+        EXPECT_DOUBLE_EQ(7.4750, fourDecimalRound(r.getData()[1]));
+        EXPECT_DOUBLE_EQ(6.8, fourDecimalRound(r.getData()[2]));
+        EXPECT_DOUBLE_EQ(5.0750, fourDecimalRound(r.getData()[3]));
+    }
 }
 
 TEST_F(Function2D_Test, fftAndComplexAndIfft){