#include <iostream>
#include "Function.h"
#include "Function2D.h"
#include "mkl.h"

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