Aurora/src/Function2D.cpp

#include <iostream>
#include "Function.h"
#include "Function2D.h"
#include "Function1D.h"
//必须在Eigen之前
#include "AuroraDefs.h"

#include <Eigen/Core>
#include <Eigen/Eigen>
#include <Eigen/Dense>

using namespace Aurora;

double Aurora::immse(const Aurora::Matrix &aImageA, const Aurora::Matrix &aImageB) {
    if (aImageA.getDims()!=2|| aImageB.getDims()!=2){
        std::cerr<<"Fail!immse args must all 2d matrix!";
        return 0.0;
    }
    if (!aImageB.compareShape(aImageA)){
        std::cerr<<"Fail!immse args must be same shape!";
        return 0.0;
    }
    if (aImageA.getValueType()!=Normal || aImageB.getValueType() != Normal) {
        std::cerr << "Fail!immse args must be normal value type!";
        return 0.0;
    }
    int size = aImageA.getDataSize();
    auto temp = malloc(size);
    vdSub(size, aImageA.getData(), aImageB.getData(), temp);
    vdSqr(size, temp, temp);
    double result = cblas_dasum(size, temp, 1) / (double) size;
    free(temp);
    return result;
}

Aurora::Matrix Aurora::inv(const Aurora::Matrix &aMatrix) {
    if (aMatrix.getDims() != 2) {
        std::cerr << "Fail!inv args must be 2d matrix!";
        return aMatrix;
    }
    if (aMatrix.getDimSize(0) != aMatrix.getDimSize(1)) {
        std::cerr << "Fail!inv args must be square matrix!";
        return aMatrix;
    }
    if (aMatrix.getValueType() != Normal) {
        std::cerr << "Fail!inv args must be normal value type!";
        return aMatrix;
    }
    int size = aMatrix.getDataSize();
    int *ipiv = new int[aMatrix.getDimSize(0)];
    auto result = malloc(size);
    cblas_dcopy(size,aMatrix.getData(), 1,result, 1);
    LAPACKE_dgetrf(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), aMatrix.getDimSize(0), result, aMatrix.getDimSize(0), ipiv);
    LAPACKE_dgetri(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), result, aMatrix.getDimSize(0), ipiv);
    delete[] ipiv;
    return Matrix::New(result,aMatrix);
}

Aurora::Matrix Aurora::inv(Aurora::Matrix&& aMatrix) {
    if (aMatrix.getDims() != 2) {
        std::cerr << "Fail!inv args must be 2d matrix!";
        return aMatrix;
    }
    if (aMatrix.getDimSize(0) != aMatrix.getDimSize(1)) {
        std::cerr << "Fail!inv args must be square matrix!";
        return aMatrix;
    }
    if (aMatrix.getValueType() != Normal) {
        std::cerr << "Fail!inv args must be normal value type!";
        return aMatrix;
    }
    int *ipiv = new int[aMatrix.getDimSize(0)];
    LAPACKE_dgetrf(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), aMatrix.getDimSize(0), aMatrix.getData(), aMatrix.getDimSize(0), ipiv);
    LAPACKE_dgetri(LAPACK_ROW_MAJOR, aMatrix.getDimSize(0), aMatrix.getData(), aMatrix.getDimSize(0), ipiv);
    delete[] ipiv;
    return aMatrix;
}

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

Matrix Aurora::std(const Matrix &aMatrix) {
    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.deepCopy();
    int calc_size = src.getDimSize(0) == 1 ? src.getDimSize(1) : src.getDimSize(0);
    int col = src.getDimSize(0) == 1?1:src.getDimSize(1) ;
    auto std = Aurora::malloc(col);
    for (int i = 0; i < col; ++i) {
        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(col, std, std);
    return Matrix::New(std, 1, col, aMatrix.getDimSize(2));
}

Matrix Aurora::min(const Matrix &aMatrix, FunctionDirection direction, long& rowIdx, long& colIdx) {
    if (aMatrix.getDimSize(2)>1 || aMatrix.isComplex()) {
        std::cerr
                << (aMatrix.getDimSize(2) > 1 ? "min() not support 3D data!" : "min() not support complex value type!")
                << std::endl;
        return Matrix();
    }
    //针对向量行等于列
    if (direction == Column && aMatrix.getDimSize(0)==1){
        direction = All;
    }
    switch (direction)
    {
        case All: {
            Eigen::Map<Eigen::MatrixXd> retV(aMatrix.getData(), aMatrix.getDimSize(0), aMatrix.getDimSize(1));
            double *ret = malloc(1);
            ret[0] = retV.array().minCoeff(&rowIdx, &colIdx);
            return Matrix::New(ret, 1);
        }
        case Row:
        {
            Eigen::Map<Eigen::MatrixXd> srcMatrix(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
            double * ret = malloc(aMatrix.getDimSize(0));

            if (aMatrix.getDimSize(0) == 1){
                ret[0] = srcMatrix.topRows(0).minCoeff(&rowIdx, &colIdx);
            }
            else{
                Eigen::Map<Eigen::MatrixXd> retMatrix(ret,aMatrix.getDimSize(0),1);
                retMatrix = srcMatrix.rowwise().minCoeff();
            }
            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(1));

            if (aMatrix.getDimSize(1) == 1){
                ret[0] = srcMatrix.col(0).minCoeff(&rowIdx, &colIdx);
            }
            else {
                Eigen::Map<Eigen::MatrixXd> retMatrix(ret,1,aMatrix.getDimSize(1));
                retMatrix = srcMatrix.colwise().minCoeff();
            }
            return Matrix::New(ret,1,aMatrix.getDimSize(1));
        }
    }
}

Matrix Aurora::min(const Matrix &aMatrix, FunctionDirection direction) {
    long a,b;
    return  min(aMatrix,direction,a,b);
}

Matrix Aurora::min(const Matrix &aMatrix, const Matrix &aOther) {
    if (aMatrix.getDimSize(2)>1 || aMatrix.isComplex()) {
        std::cerr
                << (aMatrix.getDimSize(2) > 1 ? "min() not support 3D data!" : "min() not support complex value type!")
                << std::endl;
        return Matrix();
    }
    if (aOther.getDimSize(2)>1 || aOther.isComplex()) {
        std::cerr
                << (aOther.getDimSize(2) > 1 ? "min() not support 3D data!" : "min() not support complex value type!")
                << std::endl;
        return Matrix();
    }
    //same shape
    if (aMatrix.compareShape(aOther)){
        double* output = malloc(aMatrix.getDataSize());
        vdFminI(aMatrix.getDataSize(),aMatrix.getData(),1,aOther.getData(),1,output,1);
        return Matrix::New(output,aMatrix);
    }
    // one is scalar
    else if (aMatrix.getDataSize() == 1 || aOther.getDataSize() == 1){
            double scalar = (aMatrix.getDataSize() == 1)?aMatrix.getData()[0]:aOther.getData()[0];
            auto matrix = (aMatrix.getDataSize() == 1)?aOther:aMatrix;
            double* output = malloc(matrix.getDataSize());
            vdFminI(matrix.getDataSize(),matrix.getData(),1,&scalar,0,output,1);
            return Matrix::New(output,matrix);
    }
    else if (aMatrix.getDimSize(1) == 1 || aOther.getDimSize(0) == 1) {
        if (aMatrix.getDimSize(1) == 1){
            double* output = malloc(aOther.getDataSize());
            for (int i = 0; i < aOther.getDimSize(1); ++i) {
                vdFminI(aMatrix.getDataSize(), aMatrix.getData(), 1, aOther.getData() + aOther.getDimSize(0) * i, 1,
                        output + aOther.getDimSize(0) * i, 1);
            }
            return Matrix::New(output,aOther);
        }
        else{
            double* output = malloc(aMatrix.getDataSize());
            for (int i = 0; i < aMatrix.getDimSize(0); ++i) {
                vdFminI(aOther.getDataSize(), aOther.getData(), 1, aMatrix.getData() +  i, aMatrix.getDimSize(0),
                        output +  i, aOther.getDimSize(0));
            }
            return Matrix::New(output,aMatrix);
        }
    }
    else{
        std::cerr
                << "min(A,B) with matrix must be like A[MxN] - B[1xN] or A[Mx1] - B[MxN]"
                << std::endl;
        return Matrix();
    }
}

Matrix Aurora::max(const Matrix &aMatrix, FunctionDirection direction) {
    long a,b;
    return  max(aMatrix,direction,a,b);
}

Matrix Aurora::max(const Matrix &aMatrix, FunctionDirection direction, long& rowIdx, long& colIdx) {
    if (aMatrix.getDimSize(2)>1 || aMatrix.isComplex()) {
        std::cerr
                << (aMatrix.getDimSize(2) > 1 ? "max() not support 3D data!" : "max() not support complex value type!")
                << std::endl;
        return Matrix();
    }
    //针对向量行等于列
    if (direction == Column && aMatrix.getDimSize(0)==1){
        direction = All;
    }
    switch (direction)
    {
        case All:
        {
            Eigen::Map<Eigen::MatrixXd> retV(aMatrix.getData(), aMatrix.getDimSize(0), aMatrix.getDimSize(1));
            double *ret = malloc(1);
            ret[0] = retV.array().maxCoeff(&rowIdx, &colIdx);
            return Matrix::New(ret,1);
        }
        case Row:
        {
            Eigen::Map<Eigen::MatrixXd> srcMatrix(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
            double * ret = malloc(aMatrix.getDimSize(0));
            if (aMatrix.getDimSize(0) == 1){
                ret[0] = srcMatrix.topRows(0).maxCoeff(&rowIdx, &colIdx);
            }
            else{
                Eigen::Map<Eigen::MatrixXd> retMatrix(ret,aMatrix.getDimSize(0),1);
                retMatrix = srcMatrix.rowwise().maxCoeff();
            }
            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));
            if (aMatrix.getDimSize(1) == 1){
                ret[0] = srcMatrix.col(0).maxCoeff(&rowIdx, &colIdx);
            }
            else {
                Eigen::Map<Eigen::MatrixXd> retMatrix(ret,1,aMatrix.getDimSize(1));
                retMatrix = srcMatrix.colwise().maxCoeff();
            }
            return Matrix::New(ret,1,aMatrix.getDimSize(1));
        }
    }
}

Matrix Aurora::sum(const Matrix &aMatrix, FunctionDirection direction) {
    if (aMatrix.getDimSize(2)>1 ) {
        std::cerr<<  "sum() not support 3D data!"
                << std::endl;
        return Matrix();
    }
    //针对向量行等于列
    if (direction == Column && aMatrix.getDimSize(0)==1){
        direction = Row;
    }
    if (aMatrix.isComplex()){
        switch (direction)
        {
            case All:
            {
                Eigen::Map<Eigen::VectorXcd> srcV((std::complex<double>*)aMatrix.getData(),aMatrix.getDataSize());
                std::complex<double>* ret = (std::complex<double>*)malloc(1,true);
                ret[0] = srcV.array().sum();
                return Matrix::New((double*)ret,1,1,1,Complex);
            }
            case Row:
            {
                Eigen::Map<Eigen::MatrixXcd> srcM((std::complex<double>*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
                std::complex<double> * ret = (std::complex<double>*)malloc(aMatrix.getDimSize(0),true);
                Eigen::Map<Eigen::VectorXcd> retV(ret,aMatrix.getDimSize(0));
                retV = srcM.rowwise().sum();
                return Matrix::New((double*)ret,aMatrix.getDimSize(0),1,1,Complex);
            }
            case Column:
            default:
            {
                Eigen::Map<Eigen::MatrixXcd> srcM((std::complex<double>*)aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
                std::complex<double>* ret = (std::complex<double>*)malloc(aMatrix.getDimSize(1),true);
                Eigen::Map<Eigen::VectorXcd> retV(ret,aMatrix.getDimSize(1));
                retV = srcM.colwise().sum();
                return Matrix::New((double*)ret,1,aMatrix.getDimSize(1),1,Complex);
            }
        }
    }
    else{
        switch (direction)
        {
            case All:
            {
                Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData(),aMatrix.getDataSize());
                double * ret = malloc(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));
                Eigen::Map<Eigen::VectorXd> retV(ret,aMatrix.getDimSize(0));
                retV = srcM.rowwise().sum();
                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(1));
                retV = srcM.colwise().sum();
                return Matrix::New(ret,1,aMatrix.getDimSize(1));
            }
        }
    }
}

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 (direction == Column && aMatrix.getDimSize(0)==1){
        direction = All;
    }
    if (aIncludeNan){
        switch (direction)
        {
            case All:
            {
                Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData(),aMatrix.getDataSize());
                double * ret = malloc(64);
                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(1));
                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));
            }
        }
    }
}

Matrix Aurora::sort(const Matrix &aMatrix) {
    if (aMatrix.getDimSize(2)>1 || aMatrix.isComplex()) {
        std::cerr
                << (aMatrix.getDimSize(2) > 1 ? "sort() not support 3D data!" : "sort() not support complex value type!")
                << std::endl;
        return Matrix();
    }
    return sort(std::forward<Matrix &&>(aMatrix.deepCopy()));
}

Matrix Aurora::sort(Matrix &&aMatrix) {
    if (aMatrix.getDimSize(2)>1 || aMatrix.isComplex()) {
        std::cerr
                << (aMatrix.getDimSize(2) > 1 ? "sort() not support 3D data!" : "sort() not support complex value type!")
                << std::endl;
        return Matrix();
    }
    //针对向量行等于列
    if (aMatrix.getDimSize(0)==1){
        return sortrows(aMatrix);
    }
    if (aMatrix.getDimSize(0)>=100000){
        #pragma omp parallel for
        for (int i = 0; i < aMatrix.getDimSize(1); ++i) {
            Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData()+i*aMatrix.getDimSize(0),aMatrix.getDimSize(0));
            std::sort(srcV.array().begin(),srcV.array().end());
        }
    }
    else
    {
        for (int i = 0; i < aMatrix.getDimSize(1); ++i) {
            Eigen::Map<Eigen::VectorXd> srcV(aMatrix.getData()+i*aMatrix.getDimSize(0),aMatrix.getDimSize(0));
            std::sort(srcV.array().begin(),srcV.array().end());
        }
    }

    return aMatrix;
}

Matrix Aurora::sortrows(const Matrix &aMatrix) {
    if (aMatrix.getDimSize(2)>1 || aMatrix.isComplex()) {
        std::cerr
                << (aMatrix.getDimSize(2) > 1 ? "sortrows() not support 3D data!" : "sortrows() not support complex value type!")
                << std::endl;
        return Matrix();
    }
    return sortrows(std::forward<Matrix &&>(aMatrix.deepCopy()));
}

Matrix Aurora::sortrows(Matrix &&aMatrix) {
    if (aMatrix.getDimSize(2)>1 || aMatrix.isComplex()) {
        std::cerr
                << (aMatrix.getDimSize(2) > 1 ? "sortrows() not support 3D data!" : "sortrows() not support complex value type!")
                << std::endl;
        return Matrix();
    }
    Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),aMatrix.getDimSize(0),aMatrix.getDimSize(1));
    if (aMatrix.getDimSize(1)>=100000){
        #pragma omp parallel for
        for (int i = 0; i < aMatrix.getDimSize(0); ++i) {
            std::sort(srcM.row(i).array().begin(),srcM.row(i).array().end());
        }
    }
    else
    {
        for (int i = 0; i < aMatrix.getDimSize(0); ++i) {
            std::sort(srcM.row(i).array().begin(),srcM.row(i).array().end());
        }
    }

    return aMatrix;
}

Matrix Aurora::median(const Matrix &aMatrix) {
    if (aMatrix.getDimSize(2) > 1 || aMatrix.isComplex()) {
        std::cerr
                << (aMatrix.getDimSize(2) > 1 ? "median() not support 3D data!"
                                              : "median() not support complex value type!")
                << std::endl;
        return Matrix();
    }
    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<Eigen::MatrixXd> srcM(sorted.getData(),rows,cols);
    bool flag = rows % 2 == 1;
    double* ret = malloc(cols);
    Eigen::Map<Eigen::VectorXd> retV(ret,cols);
    if (flag) {
        retV = srcM.row(rows/2);
        return Matrix::New(ret,1,cols);
    } else {
        retV = (srcM.row(rows/2-1).array()+srcM.row(rows/2).array())/2;
        return Matrix::New(ret,1,cols);
    }
}

Matrix Aurora::fft(const Matrix &aMatrix) {
    double *output = nullptr;
    mkl_free_buffers();
    output = malloc(aMatrix.getDataSize(), true);
    if (!aMatrix.isComplex()) {
        cblas_dcopy(aMatrix.getDataSize(), aMatrix.getData(), 1, output, 2);
        double zero = 0.0;
        cblas_dcopy(aMatrix.getDataSize(), &zero, 0, output+1, 2);
    } else {
        cblas_zcopy(aMatrix.getDataSize(), aMatrix.getData(), 1, output, 1);
    }

    DFTI_DESCRIPTOR_HANDLE my_desc_handle = NULL;
    MKL_LONG status;

    //创建 Descriptor, 精度 double , 输入类型实数, 维度1
    status = DftiCreateDescriptor(&my_desc_handle, DFTI_DOUBLE, DFTI_COMPLEX, 1, aMatrix.getDimSize(0));
    if (status != DFTI_NO_ERROR) goto error;

    //通过 setValue 配置Descriptor
    //使用单独的输出数据缓存
    status = DftiSetValue(my_desc_handle, DFTI_PLACEMENT, DFTI_INPLACE);
    if (status != DFTI_NO_ERROR) goto error;

    //每个傅里叶变换的输入距离
    status = DftiSetValue(my_desc_handle,DFTI_INPUT_DISTANCE,aMatrix.getDimSize(0));
    if (status != DFTI_NO_ERROR) goto error;

    //每个傅里叶变换的输出距离
    status = DftiSetValue(my_desc_handle,DFTI_OUTPUT_DISTANCE,aMatrix.getDimSize(0));
    if (status != DFTI_NO_ERROR) goto error;

    //傅里叶变换的数量
    status = DftiSetValue(my_desc_handle,DFTI_NUMBER_OF_TRANSFORMS,aMatrix.getDimSize(1));
    if (status != DFTI_NO_ERROR) goto error;

    //提交 修改配置后的Descriptor(实际上会进行FFT的计算初始化)
    status = DftiCommitDescriptor(my_desc_handle);
    if (status != DFTI_NO_ERROR) goto error;

    //执行计算
    status = DftiComputeForward(my_desc_handle, output, output);
    if (status != DFTI_NO_ERROR) goto error;

    //释放资源
    status = DftiFreeDescriptor(&my_desc_handle);
    if (status != DFTI_NO_ERROR) goto error;
    mkl_free_buffers();
    return Matrix::New(output, aMatrix.getDimSize(0), aMatrix.getDimSize(1), aMatrix.getDimSize(2), Complex);
    error:
        std::cerr<<"FFT fail, error message:"<<DftiErrorMessage(status)<<std::endl;
        return Matrix();
}

Matrix Aurora::ifft(const Matrix &aMatrix) {
    if (!aMatrix.isComplex()){
        std::cerr<<"ifft input must be complex value"<<std::endl;
        return Matrix();
    }
    DFTI_DESCRIPTOR_HANDLE my_desc_handle = NULL;
    mkl_free_buffers();
    auto output = malloc(aMatrix.getDataSize(),true);
    MKL_LONG status;

    //创建 Descriptor, 精度 double , 输入类型实数, 维度1
    int size = aMatrix.getDimSize(0);
    status = DftiCreateDescriptor(&my_desc_handle, DFTI_DOUBLE, DFTI_COMPLEX, 1, size);
    if (status != DFTI_NO_ERROR) goto error;
    //通过 setValue 配置Descriptor
    //使用单独的输出数据缓存
    status = DftiSetValue(my_desc_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
    if (status != DFTI_NO_ERROR) goto error;

    //设置DFTI_BACKWARD_SCALE ！！！很关键，不然值不对
    status = DftiSetValue(my_desc_handle, DFTI_BACKWARD_SCALE, 1.0f / size);
    if (status != DFTI_NO_ERROR) goto error;

    status = DftiSetValue(my_desc_handle,DFTI_INPUT_DISTANCE,size);
    if (status != DFTI_NO_ERROR) goto error;

    //每个傅里叶变换的输出距离
    status = DftiSetValue(my_desc_handle,DFTI_OUTPUT_DISTANCE,size);
    if (status != DFTI_NO_ERROR) goto error;

    //傅里叶变换的数量
    status = DftiSetValue(my_desc_handle,DFTI_NUMBER_OF_TRANSFORMS,aMatrix.getDimSize(1));
    if (status != DFTI_NO_ERROR) goto error;

    status = DftiCommitDescriptor(my_desc_handle);
    if (status != DFTI_NO_ERROR) goto error;

    //提交 修改配置后的Descriptor(实际上会进行FFT的计算初始化)
    status = DftiCommitDescriptor(my_desc_handle);
    if (status != DFTI_NO_ERROR) goto error;
    
    //执行计算
    status = DftiComputeBackward(my_desc_handle, aMatrix.getData(), output);
    if (status != DFTI_NO_ERROR) goto error;
    //释放资源
    status = DftiFreeDescriptor(&my_desc_handle);
    mkl_free_buffers();
    if (status != DFTI_NO_ERROR) goto error;
    {
        return Matrix::New(output, aMatrix.getDimSize(0), aMatrix.getDimSize(1), aMatrix.getDimSize(2), Complex);
    }
    error:
    std::cerr<<"FFT fail, error message:"<<DftiErrorMessage(status)<<std::endl;
    return Matrix();
}

Matrix Aurora::ifft_symmetric(const Matrix &aMatrix,long length)
{
    if(!aMatrix.isVector()){
        std::cerr<<"ifft_symmetric only support vector!"<<std::endl;
        return Matrix();
    }
    int matrixLength = aMatrix.getDataSize();
    int resultHalfLength = (int)std::ceil(((double)length*0.5));
    int copyLength = resultHalfLength<matrixLength?resultHalfLength:matrixLength;
    double* calcData = malloc(length,true);
    double zero = 0.0;
    //所有数据统一置0
    cblas_dcopy(length*2,&zero,0,calcData,1);
    //copy前半段数据
    cblas_zcopy(copyLength,aMatrix.getData(),1,calcData,1);
    //copy后半段数据,跳过index 0的值,并设置虚部共轭
    vdAddI(copyLength-1,&zero,0,(aMatrix.getData()+2),2,(calcData+(length-1)*2),-2);
    vdSubI(copyLength-1,&zero,0,(aMatrix.getData()+2+1),2,(calcData+(length-1)*2+1),-2);

    return real(ifft(Matrix::New(calcData,length,1,1,Complex)));
}

Matrix Aurora::hilbert(const Matrix &aMatrix) {
    auto x = fft(aMatrix);
    auto h = malloc(aMatrix.getDimSize(0));
    auto two = 2.0;
    auto zero = 0.0;
    cblas_dcopy(aMatrix.getDimSize(0), &zero, 0, h, 1);
    cblas_dcopy(aMatrix.getDimSize(0) / 2, &two, 0, h, 1);
    h[aMatrix.getDimSize(0) / 2] = ((aMatrix.getDimSize(0) << 31) >> 31) ? 2.0 : 1.0;
    h[0] = 1.0;
    for (int i = 0; i < aMatrix.getDimSize(1); ++i) {
        auto p = (double *)(x.getData() + aMatrix.getDimSize(0)* i*2);
        vdMulI(aMatrix.getDimSize(0), p, 2, h, 1, p, 2);
        vdMulI(aMatrix.getDimSize(0), p + 1, 2, h, 1, p + 1, 2);
    }
    auto result = ifft( x);
    free(h);
    return result;
}

Matrix Aurora::prod(const Matrix &aMatrix) {
    if (aMatrix.getDimSize(2) > 1 ) {
        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<Eigen::MatrixXcd> srcM((std::complex<double>*)aMatrix.getData(),row,col);
        auto ret = malloc(col,true);
        Eigen::Map<Eigen::VectorXcd> retV((std::complex<double>*)ret,col);
        retV = srcM.colwise().prod();
        return Matrix::New(ret,1,col,1,Complex);
    }
    else{
        Eigen::Map<Eigen::MatrixXd> srcM(aMatrix.getData(),row,col);
        auto ret = malloc(col);
        Eigen::Map<Eigen::VectorXd> retV(ret,col);
        retV = srcM.colwise().prod();
        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;
        return Matrix();
    }
    if (!aMatrix.compareShape(aOther)){
        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);
    }
    else{
        if (direction == Column)
        {
            auto ret = malloc(aMatrix.getDimSize(1));
            for (int i = 0; i < aMatrix.getDimSize(1); ++i) {
                ret[i]=cblas_ddot(aMatrix.getDimSize(0),aMatrix.getData()+i*aMatrix.getDimSize(0),1,
                                aOther.getData()+i*aMatrix.getDimSize(0),1);
            }
            return Matrix::New(ret,1,aMatrix.getDimSize(1),1);
        }
        else{
            auto ret = malloc(aMatrix.getDimSize(0));
            for (int i = 0; i < aMatrix.getDimSize(0); ++i) {
                ret[i] = cblas_ddot(aMatrix.getDimSize(1),aMatrix.getData()+i,aMatrix.getDimSize(0),
                                aOther.getData()+i,aMatrix.getDimSize(0));
            }
            return Matrix::New(ret,aMatrix.getDimSize(1),1,1);
        }
    }
}