Add interp and repmat function.

src/Function1D.cpp

@@ -1,4 +1,5 @@
 #include "Function1D.h"
+#include "Function.h"
 #include <complex>
 #include <cstring>
 #include <iostream>
@@ -10,6 +11,8 @@
 #include <Eigen/Eigen>
 #include <Eigen/Dense>
 
+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<double> resultData = std::shared_ptr<double>(Aurora::malloc(nx1), Aurora::free);
+    std::vector<int> 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<column; ++i)
+    {
+        for(int j=1; j<=aRowTimes; ++j)
+        {
+            std::copy(originalData, originalData+row, resultDataTemp);
+            resultDataTemp += row;
+        }
+        originalData += row;
+    }
+
+    resultDataTemp = resultData;
+    int step = originalDataSize * aRowTimes;
+    for(int i=1; i<aColumnTimes; ++i)
+    {
+        std::copy(resultData, resultData + step, resultData + i*step);
+    }
+
+    std::vector<int> resultInfo;
+    resultInfo.push_back(aRowTimes * row);
+    column = aColumnTimes*column;
+    if (column > 1)
+    {
+        resultInfo.push_back(column);
+    }
+
+    return Matrix(std::shared_ptr<double>(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<aSliceTimes; ++i)
+    {
+        std::copy(resultData, resultData + resultTempDataSize, resultData + i*resultTempDataSize);
+    }
+    std::vector<int> 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<double>(resultData, Aurora::free),resultInfo);
+}

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);
 };
 
 

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<double> resultData = std::shared_ptr<double>(Aurora::malloc(nx1), Aurora::free);
+    for (int i = 0; i < nx1; ++i)
+    {
+        std::shared_ptr<double> xResultData = std::shared_ptr<double>(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<int>{columnNum});
+        resultData.get()[i] = interp1(aX,xResult,aX1(i).toMatrix(),aMethod).getData()[0];
+    }
+    return Matrix(resultData,std::vector<int>{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);
+}

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);
 
 };
 

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<double> resultData = std::shared_ptr<double>(new double[nx1], std::default_delete<double[]>());
+    for (int i = 0; i < nx1; ++i)
+    {
+        std::shared_ptr<double> zResultData = std::shared_ptr<double>(new double[zAxisNum], std::default_delete<double[]>());
+        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<int>{columnNum});
+        resultData.get()[i] = interp1(aZ,xResult,aZ1(i).toMatrix(),aMethod).getData()[0];
+    }
+    return Matrix(resultData,std::vector<int>{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);
+}

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);
 
 };
 

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<double> resultData = std::shared_ptr<double>(new double[rows], std::default_delete<double[]>());
-        std::copy(mData.get() + (aColumn - 1) * rows, mData.get() + aColumn * rows, resultData.get());
-        std::vector<int> resultInfo = {rows};
-        Matrix result(resultData, resultInfo);
-        return result;
-    }
-
-    Matrix Matrix::getDataFromDims1(int aRow) {
-        if (1 != getDims() || aRow > mInfo.back()) {
-            return Matrix();
-        }
-        std::shared_ptr<double> resultData = std::shared_ptr<double>(new double[1], std::default_delete<double[]>());
-        resultData.get()[0] = mData.get()[aRow - 1];
-        std::vector<int> 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<int> vector;

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对象