From d73dc8be23d4bbd257cdab3acc479cf3dfcc6cf2 Mon Sep 17 00:00:00 2001 From: sunwen Date: Thu, 20 Apr 2023 15:34:38 +0800 Subject: [PATCH] Add interp and repmat function. --- src/Function1D.cpp | 130 +++++++++++++++++++++++++++++++++++++++++++++ src/Function1D.h | 12 +++++ src/Function2D.cpp | 44 +++++++++++++++ src/Function2D.h | 8 ++- src/Function3D.cpp | 45 ++++++++++++++++ src/Function3D.h | 8 ++- src/Matrix.cpp | 24 --------- src/Matrix.h | 4 -- 8 files changed, 245 insertions(+), 30 deletions(-) diff --git a/src/Function1D.cpp b/src/Function1D.cpp index e4d08dc..5a44f66 100644 --- a/src/Function1D.cpp +++ b/src/Function1D.cpp @@ -1,4 +1,5 @@ #include "Function1D.h" +#include "Function.h" #include #include #include @@ -10,6 +11,8 @@ #include #include +using namespace Aurora; + namespace { const int COMPLEX_STRIDE = 2; const int REAL_STRIDE = 1; @@ -179,3 +182,130 @@ Aurora::Matrix Aurora::sign(const Aurora::Matrix&& matrix) { return matrix; } } + +Matrix Aurora::interp1(const Matrix& aX, const Matrix& aV, const Matrix& aX1, InterpnMethod aMethod) +{ + const int nx = aX.getDimSize(0); + const int ny = 1; + int nx1 = aX1.getDimSize(0); + std::shared_ptr resultData = std::shared_ptr(Aurora::malloc(nx1), Aurora::free); + std::vector resultInfo = {nx1}; + Aurora::Matrix result(resultData, resultInfo); + DFTaskPtr task ; + int status = dfdNewTask1D(&task, nx, aX.getData(), DF_NO_HINT, ny, aV.getData(), DF_NO_HINT); + if (status != DF_STATUS_OK) + { + return Matrix(); + } + MKL_INT sorder; + MKL_INT stype; + if (aMethod == InterpnMethod::Spline) + { + sorder = DF_PP_CUBIC; + stype = DF_PP_BESSEL; + } + else + { + sorder = DF_PP_LINEAR; + stype = DF_PP_BESSEL; + } + double* scoeffs = Aurora::malloc(ny * (nx-1) * sorder); + status = dfdEditPPSpline1D(task, sorder,DF_PP_NATURAL , DF_BC_NOT_A_KNOT, 0, DF_NO_IC, 0, scoeffs, DF_NO_HINT); + if (status != DF_STATUS_OK) + { + return Matrix(); + } + status = dfdConstruct1D( task, DF_PP_SPLINE, DF_METHOD_STD ); + if (status != DF_STATUS_OK) + { + return Matrix(); + } + + int dorder = 1; + status = dfdInterpolate1D(task, DF_INTERP, DF_METHOD_PP, nx1, aX1.getData(), DF_NO_HINT, 1, &dorder, + DF_NO_APRIORI_INFO, resultData.get(), DF_MATRIX_STORAGE_ROWS, nullptr); + + status = dfDeleteTask(&task); + + Aurora::free(scoeffs); + return result; +} + +Matrix Aurora::repmat(const Matrix& aMatrix,int aRowTimes, int aColumnTimes) +{ + if(aRowTimes < 1 || aColumnTimes < 1 || aMatrix.getDims() > 2 || aMatrix.isNull()) + { + return Matrix(); + } + int originalDataSize = aMatrix.getDataSize(); + double* resultData = Aurora::malloc(originalDataSize * aRowTimes * aColumnTimes); + int row = aMatrix.getDimSize(0); + int column = 1; + if(aMatrix.getDims() > 1) + { + column = aMatrix.getDimSize(1); + } + + double* originalData = aMatrix.getData(); + double* resultDataTemp = resultData; + for(int i=0; i resultInfo; + resultInfo.push_back(aRowTimes * row); + column = aColumnTimes*column; + if (column > 1) + { + resultInfo.push_back(column); + } + + return Matrix(std::shared_ptr(resultData, Aurora::free),resultInfo); +} + +Matrix Aurora::repmat(const Matrix& aMatrix,int aRowTimes, int aColumnTimes, int aSliceTimes) +{ + if(aRowTimes < 1 || aColumnTimes < 1 || aSliceTimes < 1 || aMatrix.getDims() > 2 || aMatrix.isNull()) + { + return Matrix(); + } + Matrix resultTemp = Aurora::repmat(aMatrix, aRowTimes, aColumnTimes); + int resultTempDataSize = resultTemp.getDataSize(); + double* resultData = Aurora::malloc(resultTempDataSize * aSliceTimes); + std::copy(resultTemp.getData(), resultTemp.getData() + resultTempDataSize, resultData); + for(int i=1; i resultInfo; + int row = resultTemp.getDimSize(0); + int column = 1; + if(resultTemp.getDims() > 1) + { + column = resultTemp.getDimSize(1); + } + resultInfo.push_back(row); + if (column > 1 || aSliceTimes > 1) + { + resultInfo.push_back(column); + } + if (aSliceTimes > 1) + { + resultInfo.push_back(aSliceTimes); + } + + return Matrix(std::shared_ptr(resultData, Aurora::free),resultInfo); +} diff --git a/src/Function1D.h b/src/Function1D.h index 1872cba..755461a 100644 --- a/src/Function1D.h +++ b/src/Function1D.h @@ -4,6 +4,12 @@ #include "Matrix.h" namespace Aurora { + + enum InterpnMethod + { + Spline=0,Linear + }; + Matrix complex(const Matrix& matrix); Matrix real(const Matrix& matrix); @@ -34,6 +40,12 @@ namespace Aurora { Matrix sign(const Matrix& matrix); Matrix sign(const Matrix&& matrix); + + Matrix interp1(const Matrix& aX, const Matrix& aV, const Matrix& aX1, InterpnMethod aMethod); + + Matrix repmat(const Matrix& aMatrix,int aRowTimes, int aColumnTimes); + + Matrix repmat(const Matrix& aMatrix,int aRowTimes, int aColumnTimes, int aSliceTimes); }; diff --git a/src/Function2D.cpp b/src/Function2D.cpp index f08ca06..7241bea 100644 --- a/src/Function2D.cpp +++ b/src/Function2D.cpp @@ -1 +1,45 @@ #include "Function2D.h" +#include "Function1D.h" +#include "Function.h" + +using namespace Aurora; + +Matrix Aurora::interp2(const Matrix& aX, const Matrix& aY, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, InterpnMethod aMethod) +{ + if (aV.getDims() != 2) + { + return Matrix(); + } + + if (aX1.getDimSize(0) != aY1.getDimSize(0)) + { + return Matrix(); + } + + int columnNum = aV.getDimSize(1); + int rowNum = aV.getDimSize(0); + + if(aX.getDimSize(0) != columnNum || aY.getDimSize(0) != rowNum) + { + return Matrix(); + } + + int nx1 = aX1.getDimSize(0); + std::shared_ptr resultData = std::shared_ptr(Aurora::malloc(nx1), Aurora::free); + for (int i = 0; i < nx1; ++i) + { + std::shared_ptr xResultData = std::shared_ptr(Aurora::malloc(columnNum), Aurora::free); + for(int j =0; j < columnNum; ++j) + { + xResultData.get()[j] = interp1(aY,aV($,j).toMatrix(),aY1(i).toMatrix(),aMethod).getData()[0]; + } + Matrix xResult(xResultData,std::vector{columnNum}); + resultData.get()[i] = interp1(aX,xResult,aX1(i).toMatrix(),aMethod).getData()[0]; + } + return Matrix(resultData,std::vector{nx1}); +} + +Matrix Aurora::interpn(const Matrix& aX, const Matrix& aY, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, InterpnMethod aMethod) +{ + return Aurora::interp2(aY,aX,aV,aY1,aX1,aMethod); +} diff --git a/src/Function2D.h b/src/Function2D.h index 94f8e3e..1d42f81 100644 --- a/src/Function2D.h +++ b/src/Function2D.h @@ -1,8 +1,14 @@ #ifndef AURORA_FUNCTION2D_H #define AURORA_FUNCTION2D_H +#include "Matrix.h" +#include "Function1D.h" -class Function2D { + +namespace Aurora { + + Matrix interp2(const Matrix& aX, const Matrix& aY, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, InterpnMethod aMethod); + Matrix interpn(const Matrix& aX, const Matrix& aY, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, InterpnMethod aMethod); }; diff --git a/src/Function3D.cpp b/src/Function3D.cpp index c39d6c8..40b4a2c 100644 --- a/src/Function3D.cpp +++ b/src/Function3D.cpp @@ -1 +1,46 @@ #include "Function3D.h" +#include "Function2D.h" +#include "Function.h" + +using namespace Aurora; + +Matrix Aurora::interp3(const Matrix& aX, const Matrix& aY, const Matrix& aZ, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, const Matrix& aZ1,InterpnMethod aMethod) +{ + if (aV.getDims() != 3) + { + return Matrix(); + } + + if ( aX1.getDimSize(0) != aY1.getDimSize(0) || + aX1.getDimSize(0) != aZ1.getDimSize(0) || + aZ1.getDimSize(0) != aY1.getDimSize(0)) + { + return Matrix(); + } + int nx1 = aX1.getDimSize(0); + int zAxisNum = aV.getDimSize(2); + int columnNum = aV.getDimSize(1); + int rowNum = aV.getDimSize(0); + + if(aX.getDimSize(0) != columnNum || aY.getDimSize(0) != rowNum || aZ.getDimSize(0) != zAxisNum) + { + return Matrix(); + } + std::shared_ptr resultData = std::shared_ptr(new double[nx1], std::default_delete()); + for (int i = 0; i < nx1; ++i) + { + std::shared_ptr zResultData = std::shared_ptr(new double[zAxisNum], std::default_delete()); + for(int j =0; j < zAxisNum; ++j) + { + zResultData.get()[j] = interp2(aX,aY,aV($,$,j).toMatrix(),aX1(i).toMatrix(),aY1(i).toMatrix(),aMethod).getData()[0]; + } + Matrix xResult(zResultData,std::vector{columnNum}); + resultData.get()[i] = interp1(aZ,xResult,aZ1(i).toMatrix(),aMethod).getData()[0]; + } + return Matrix(resultData,std::vector{nx1}); +} + +Matrix Aurora::interpn(const Matrix& aX, const Matrix& aY, const Matrix& aZ, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, const Matrix& aZ1,InterpnMethod aMethod) +{ + return Aurora::interp3(aY,aX,aZ,aV,aY1,aX1,aZ1,aMethod); +} diff --git a/src/Function3D.h b/src/Function3D.h index f626bca..222816d 100644 --- a/src/Function3D.h +++ b/src/Function3D.h @@ -2,7 +2,13 @@ #define AURORA_FUNCTION3D_H -class Function3D { +#include "Matrix.h" +#include "Function1D.h" + +namespace Aurora { + + Matrix interp3(const Matrix& aX, const Matrix& aY, const Matrix& aZ, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, const Matrix& aZ1,InterpnMethod aMethod); + Matrix interpn(const Matrix& aX, const Matrix& aY, const Matrix& aZ, const Matrix& aV, const Matrix& aX1, const Matrix& aY1, const Matrix& aZ1,InterpnMethod aMethod); }; diff --git a/src/Matrix.cpp b/src/Matrix.cpp index 7f1cbae..6a85017 100644 --- a/src/Matrix.cpp +++ b/src/Matrix.cpp @@ -163,30 +163,6 @@ namespace Aurora { return true; } - Matrix Matrix::getDataFromDims2(int aColumn) { - if (2 != getDims() || aColumn > mInfo.back()) { - return Matrix(); - } - - int rows = mInfo.at(0); - std::shared_ptr resultData = std::shared_ptr(new double[rows], std::default_delete()); - std::copy(mData.get() + (aColumn - 1) * rows, mData.get() + aColumn * rows, resultData.get()); - std::vector resultInfo = {rows}; - Matrix result(resultData, resultInfo); - return result; - } - - Matrix Matrix::getDataFromDims1(int aRow) { - if (1 != getDims() || aRow > mInfo.back()) { - return Matrix(); - } - std::shared_ptr resultData = std::shared_ptr(new double[1], std::default_delete()); - resultData.get()[0] = mData.get()[aRow - 1]; - std::vector resultInfo{1}; - Matrix result(resultData, resultInfo); - return result; - } - Matrix Matrix::New(double *data, int rows, int cols, int slices, ValueType type) { if (!data) return Matrix(); std::vector vector; diff --git a/src/Matrix.h b/src/Matrix.h index 55763c6..b2a61d9 100644 --- a/src/Matrix.h +++ b/src/Matrix.h @@ -94,10 +94,6 @@ namespace Aurora { */ static Matrix New(const Matrix &shapeMatrix); - Matrix getDataFromDims2(int aColumn); - - Matrix getDataFromDims1(int aRow); - /** * 深拷贝操作 * @return 深拷贝的Matrix对象