From b4423b756e8d52bd8eb3a50f750965c13c1faf08 Mon Sep 17 00:00:00 2001 From: Krad Date: Wed, 26 Apr 2023 14:58:18 +0800 Subject: [PATCH] Add column vector support to min, max, sum, prod. --- src/Function1D.h | 2 + src/Function2D.cpp | 91 +++++++++++++++++++++++++++------------- test/Function2D_Test.cpp | 27 ++++++++++++ 3 files changed, 90 insertions(+), 30 deletions(-) diff --git a/src/Function1D.h b/src/Function1D.h index 15138ea..238b570 100644 --- a/src/Function1D.h +++ b/src/Function1D.h @@ -83,6 +83,8 @@ namespace Aurora { * @return 查询结果 */ Matrix polyval(const Matrix& aP, const Matrix& aX); + + void nantoval(Matrix& aMatrix,double val); }; diff --git a/src/Function2D.cpp b/src/Function2D.cpp index a42570e..d43868b 100644 --- a/src/Function2D.cpp +++ b/src/Function2D.cpp @@ -117,18 +117,26 @@ Matrix Aurora::interpn(const Matrix& aX, const Matrix& aY, const Matrix& aV, con } Matrix Aurora::std(const Matrix &aMatrix) { - auto std = Aurora::malloc(aMatrix.getDimSize(1) * aMatrix.getDimSize(2)); + if (aMatrix.getDimSize(2) > 1 || aMatrix.isComplex()) { + std::cerr + << (aMatrix.getDimSize(2) > 1 ? "std() not support 3D data!" : "std() not support complex value type!") + << std::endl; + return Matrix(); + } + Matrix src = aMatrix.isComplex() ? Aurora::abs(aMatrix) : aMatrix; + int calc_size = src.getDimSize(0) == 1 ? src.getDimSize(1) : src.getDimSize(0); + auto std = Aurora::malloc(aMatrix.getDimSize(1)); for (int i = 0; i < src.getDimSize(1); ++i) { - double *p = src.getData() + i * src.getDimSize(0); - double mean = cblas_dasum(src.getDimSize(0), p, 1) / src.getDimSize(0); - vdSubI(src.getDimSize(0), p, 1, &mean, 0, p, 1); - vdSqr(src.getDimSize(0), p, p); - std[i] = cblas_dasum(src.getDimSize(0), p, 1) / (src.getDimSize(0) - 1); + double *p = src.getData() + i * calc_size; + double mean = cblas_dasum(calc_size, p, 1) / calc_size; + vdSubI(calc_size, p, 1, &mean, 0, p, 1); + vdSqr(calc_size, p, p); + std[i] = cblas_dasum(calc_size, p, 1) / (calc_size - 1); } vdSqrt(src.getDimSize(1), std, std); + return Matrix::New(std, 1, aMatrix.getDimSize(1), aMatrix.getDimSize(2)); - return Matrix::New(std,1,aMatrix.getDimSize(1), aMatrix.getDimSize(2)); } Matrix Aurora::min(const Matrix &aMatrix, FunctionDirection direction, long& rowIdx, long& colIdx) { @@ -138,6 +146,10 @@ Matrix Aurora::min(const Matrix &aMatrix, FunctionDirection direction, long& row << std::endl; return Matrix(); } + //针对向量行等于列 + if (direction == Column && aMatrix.getDimSize(0)==1){ + direction = All; + } switch (direction) { case All: { @@ -248,6 +260,10 @@ Matrix Aurora::max(const Matrix &aMatrix, FunctionDirection direction, long& row << std::endl; return Matrix(); } + //针对向量行等于列 + if (direction == Column && aMatrix.getDimSize(0)==1){ + direction = All; + } switch (direction) { case All: @@ -293,6 +309,10 @@ Matrix Aurora::sum(const Matrix &aMatrix, FunctionDirection direction) { << std::endl; return Matrix(); } + //针对向量行等于列 + if (direction == Column && aMatrix.getDimSize(0)==1){ + direction = Row; + } if (aMatrix.isComplex()){ switch (direction) { @@ -360,6 +380,10 @@ Matrix Aurora::mean(const Matrix &aMatrix, FunctionDirection direction, bool aIn << std::endl; return Matrix(); } + //针对向量行等于列 + if (direction == Column && aMatrix.getDimSize(0)==1){ + direction = All; + } if (aIncludeNan){ switch (direction) { @@ -465,6 +489,7 @@ Matrix Aurora::sort(const Matrix &aMatrix) { } return sort(std::forward(aMatrix.deepCopy())); } + Matrix Aurora::sort(Matrix &&aMatrix) { if (aMatrix.getDimSize(2)>1 || aMatrix.isComplex()) { std::cerr @@ -532,17 +557,20 @@ Matrix Aurora::median(const Matrix &aMatrix) { << std::endl; return Matrix(); } - Matrix sorted = sort(aMatrix); - Eigen::Map srcM(sorted.getData(),sorted.getDimSize(0),sorted.getDimSize(1)); + bool horVector = aMatrix.getDimSize(0)==1; + Matrix sorted = horVector?sortrows(aMatrix):sort(aMatrix); + int rows = horVector?sorted.getDimSize(1):sorted.getDimSize(0); + int cols = horVector?sorted.getDimSize(0):sorted.getDimSize(1); + Eigen::Map srcM(sorted.getData(),rows,cols); bool flag = aMatrix.getDimSize(0) % 2 == 1; - double* ret = malloc(aMatrix.getDimSize(1)); - Eigen::Map retV(ret,aMatrix.getDimSize(1)); + double* ret = malloc(cols); + Eigen::Map retV(ret,cols); if (flag) { - retV = srcM.row(aMatrix.getDimSize(0)/2); - return Matrix::New(ret,1,aMatrix.getDimSize(1)); + retV = srcM.row(rows/2); + return Matrix::New(ret,1,cols); } else { - retV = (srcM.row(aMatrix.getDimSize(0)/2-1).array()+srcM.row(aMatrix.getDimSize(0)/2).array())/2; - return Matrix::New(ret,1,aMatrix.getDimSize(1)); + retV = (srcM.row(rows/2-1).array()+srcM.row(rows/2).array())/2; + return Matrix::New(ret,1,cols); } } @@ -667,43 +695,46 @@ Matrix Aurora::hilbert(const Matrix &aMatrix) { vdMulI(aMatrix.getDimSize(0), p + 1, 2, h, 1, p + 1, 2); } auto result = ifft( x); - delete[] h; + free(h); return result; } Matrix Aurora::prod(const Matrix &aMatrix) { if (aMatrix.getDimSize(2) > 1 ) { - std::cerr<< "prod() not support 3D data!" - << std::endl; + std::cerr<< "prod() not support 3D data!"<< std::endl; return Matrix(); } + int row = aMatrix.getDimSize(0)==1?aMatrix.getDimSize(1):aMatrix.getDimSize(0); + int col = aMatrix.getDimSize(0)==1?1:aMatrix.getDimSize(1); if (aMatrix.isComplex()){ - Eigen::Map srcM((std::complex*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); - auto ret = malloc(aMatrix.getDimSize(1),true); - Eigen::Map retV((std::complex*)ret,aMatrix.getDimSize(1)); + Eigen::Map srcM((std::complex*)aMatrix.getData(),row,col); + auto ret = malloc(col,true); + Eigen::Map retV((std::complex*)ret,col); retV = srcM.colwise().prod(); - return Matrix::New(ret,1,aMatrix.getDimSize(1),1,Complex); + return Matrix::New(ret,1,col,1,Complex); } else{ - Eigen::Map srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1)); - auto ret = malloc(aMatrix.getDimSize(1)); - Eigen::Map retV(ret,aMatrix.getDimSize(1)); + Eigen::Map srcM(aMatrix.getData(),row,col); + auto ret = malloc(col); + Eigen::Map retV(ret,col); retV = srcM.colwise().prod(); - return Matrix::New(ret,1,aMatrix.getDimSize(1)); + return Matrix::New(ret,1,col); } } Matrix Aurora::dot(const Matrix &aMatrix,const Matrix& aOther,FunctionDirection direction ) { if ( direction == All){ - std::cerr<< "dot() not support 3D data!" - << std::endl; + std::cerr<< "dot() not support 3D data!"<< std::endl; return Matrix(); } if (!aMatrix.compareShape(aOther)){ - std::cerr<< "dot() matrix must be same shape!" - << std::endl; + std::cerr<< "dot() matrix must be same shape!"<< std::endl; return Matrix(); } + //针对向量行等于列 + if (direction == Column && aMatrix.getDimSize(0)==1){ + direction = Row; + } if (aMatrix.isComplex()){ return sum(conj(aMatrix)*aOther,direction); } diff --git a/test/Function2D_Test.cpp b/test/Function2D_Test.cpp index 0c2ff6e..9f62724 100644 --- a/test/Function2D_Test.cpp +++ b/test/Function2D_Test.cpp @@ -102,6 +102,10 @@ TEST_F(Function2D_Test, min) { EXPECT_DOUBLE_EQ(1, ret.getData()[0]); EXPECT_EQ(0, r); EXPECT_EQ(0, c); + ret = Aurora::min(D); + EXPECT_EQ(1, ret.getDimSize(0)); + EXPECT_EQ(1, ret.getDimSize(1)); + EXPECT_DOUBLE_EQ(1, ret.getData()[0]); ret = Aurora::min(A, C); EXPECT_DOUBLE_EQ(3, ret.getDataSize()); EXPECT_DOUBLE_EQ(1, ret.getData()[0]); @@ -141,6 +145,11 @@ TEST_F(Function2D_Test, max) { EXPECT_DOUBLE_EQ(3, ret.getData()[0]); EXPECT_EQ(2, r); EXPECT_EQ(0, c); + auto D = Aurora::Matrix::copyFromRawData(dataA, 1, 3); + ret = Aurora::max(D); + EXPECT_EQ(1, ret.getDimSize(0)); + EXPECT_EQ(1, ret.getDimSize(1)); + EXPECT_DOUBLE_EQ(3, ret.getData()[0]); } TEST_F(Function2D_Test, sum) { @@ -162,6 +171,12 @@ TEST_F(Function2D_Test, sum) { EXPECT_DOUBLE_EQ(7, ret.getData()[0]); EXPECT_DOUBLE_EQ(8, ret.getData()[1]); EXPECT_DOUBLE_EQ(1, ret.getData()[2]); + double *dataA = new double[3]{1, 2, 3}; + auto D = Aurora::Matrix::copyFromRawData(dataA, 1, 3); + ret = Aurora::sum(D); + EXPECT_EQ(1, ret.getDimSize(0)); + EXPECT_EQ(1, ret.getDimSize(1)); + EXPECT_DOUBLE_EQ(6, ret.getData()[0]); } TEST_F(Function2D_Test, mean) { @@ -308,6 +323,12 @@ TEST_F(Function2D_Test, median) { EXPECT_DOUBLE_EQ(9, ret.getData()[2]); EXPECT_DOUBLE_EQ(15.55, ret.getData()[3]); EXPECT_DOUBLE_EQ(7.25, ret.getData()[4]); + double *dataA = new double[3]{1, 2, 3}; + auto D = Aurora::Matrix::copyFromRawData(dataA, 1, 3); + ret = Aurora::median(D); + EXPECT_EQ(1, ret.getDimSize(0)); + EXPECT_EQ(1, ret.getDimSize(1)); + EXPECT_DOUBLE_EQ(2, ret.getData()[0]); } TEST_F(Function2D_Test, fftAndComplexAndIfft){ @@ -357,6 +378,12 @@ TEST_F(Function2D_Test, prod) { EXPECT_DOUBLE_EQ(0.6402, fourDecimalRound(ret.getData()[2]/10000)); EXPECT_DOUBLE_EQ(-4.8307, fourDecimalRound(ret.getData()[3]/10000)); EXPECT_DOUBLE_EQ(0.2608, fourDecimalRound(ret.getData()[4]/10000)); + double *dataA = new double[3]{1, 2, 3}; + auto D = Aurora::Matrix::copyFromRawData(dataA, 1, 3); + ret = Aurora::prod(D); + EXPECT_EQ(1, ret.getDimSize(0)); + EXPECT_EQ(1, ret.getDimSize(1)); + EXPECT_DOUBLE_EQ(6, ret.getData()[0]); } TEST_F(Function2D_Test, dot) {