From a895785319906b48780dbebf1a341ff234e5ca3d Mon Sep 17 00:00:00 2001 From: Krad Date: Mon, 24 Apr 2023 15:24:23 +0800 Subject: [PATCH] Add mean and mean's unit test. --- src/Function2D.cpp | 132 ++++++++++++++++++++++++++++++++++++--- src/Function2D.h | 9 +++ test/Function2D_Test.cpp | 96 +++++++++++++++++++++++++--- 3 files changed, 219 insertions(+), 18 deletions(-) diff --git a/src/Function2D.cpp b/src/Function2D.cpp index 91c2e8c..d4aec69 100644 --- a/src/Function2D.cpp +++ b/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 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 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 srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); double * ret = malloc(aMatrix.getDimSize(0)); - for (int i = 0; i < aMatrix.getDimSize(0); ++i) { - ret[i] = cblas_dasum(aMatrix.getDimSize(1), aMatrix.getData() + i, - aMatrix.getDimSize(0)); - } + Eigen::Map retV(ret,aMatrix.getDimSize(0)); + retV = srcM.rowwise().sum(); return Matrix::New(ret,aMatrix.getDimSize(0),1); } case Column: default: { - double * ret = malloc(aMatrix.getDimSize(0)); - for (int i = 0; i < aMatrix.getDimSize(1); ++i) { - ret[i] = cblas_dasum(aMatrix.getDimSize(0), aMatrix.getData()+aMatrix.getDimSize(0)*i, - 1); - } + Eigen::Map srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); + double * ret = malloc(aMatrix.getDimSize(1)); + Eigen::Map retV(ret,aMatrix.getDimSize(0)); + 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 srcV(aMatrix.getData(),aMatrix.getDataSize()); + double * ret = malloc(1); + ret[0] = srcV.array().mean(); + return Matrix::New(ret,1); + } + case Row: + { + Eigen::Map srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); + double * ret = malloc(aMatrix.getDimSize(0)); + Eigen::Map retV(ret,aMatrix.getDimSize(0)); + retV = srcM.rowwise().mean(); + return Matrix::New(ret,aMatrix.getDimSize(0),1); + } + case Column: + default: + { + Eigen::Map srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); + double * ret = malloc(aMatrix.getDimSize(1)); + Eigen::Map retV(ret,aMatrix.getDimSize(0)); + retV = srcM.colwise().mean(); + return Matrix::New(ret,1,aMatrix.getDimSize(1)); + } + } + } + else{ + switch (direction) + { + case All: + { + Eigen::Map srcV(aMatrix.getData(),aMatrix.getDataSize()); + double * retVd = malloc(aMatrix.getDataSize()); + Eigen::Map 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 srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); + double * retMd = malloc(aMatrix.getDataSize()); + Eigen::Map 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 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 srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); + double * retMd = malloc(aMatrix.getDataSize()); + Eigen::Map 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 retV(ret,aMatrix.getDimSize(1)); + retV = retM.colwise().sum(); + retV = retV.array()/countM.array(); + free(retMd); + return Matrix::New(ret,1,aMatrix.getDimSize(1)); + } + } + } +} \ No newline at end of file diff --git a/src/Function2D.h b/src/Function2D.h index 77b697a..9830463 100644 --- a/src/Function2D.h +++ b/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); + }; diff --git a/test/Function2D_Test.cpp b/test/Function2D_Test.cpp index 566db96..66647ec 100644 --- a/test/Function2D_Test.cpp +++ b/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){