#define EIGEN_USE_MKL_ALL

#include <iostream>
#include <algorithm>


#include <complex>

#include "Matrix.h"
#include "Function.h"
#include "Function1D.h"
#include "Function2D.h"
#include "Function3D.h"
#include "MatlabReader.h"

int main()
{
    MatlabReader m;
    Input i;
    MatchedFilter o;

    m.read(&i, &o, "/home/krad/TestData/testCreateMatchedFilter.mat");
    bool measuredCEused = true, findDefects = i.mParams->mFindDefects, removeOutliers = i.mRemoveOutliersFromCEMeasured;

    Aurora::Matrix mFTime = Aurora::Matrix::fromRawData(i.mCe, 4000, 2304);
    if (removeOutliers)
    {
        auto normSTD = Aurora::std(Aurora::abs(mFTime));
        Aurora::nantoval(normSTD, 0.0);
        auto sortSTD = Aurora::sort(normSTD);
        int t = (int)std::round(0.4 * mFTime.getDimSize(1)) - 1;
        t = t <= 0 ? 1.0 : t;
        auto absFTime = abs(mFTime);
        auto maxAbsFTime = max(absFTime);
        auto maxFlag = maxAbsFTime < (0.1 * max(maxAbsFTime));
        auto lessFlag = normSTD < sortSTD(0, t).toMatrix();
        long maxCol, maxRow;
        max(normSTD, Aurora::Column, maxRow, maxCol);
        for (int j = 0; j < mFTime.getDimSize(1); ++j)
        {
            if ((bool)(lessFlag.getData()[j]) || (bool)(maxFlag.getData()[j]))
            {
                mFTime(Aurora::$, j) = mFTime(Aurora::$, maxCol);
            }
        }
    }
    auto matchedFilter = fft(mFTime);
    auto sumDiff = Aurora::zeros(1, matchedFilter.getDimSize(1));
    long minCol = 998, row;

    auto absMatchedFilter = abs(matchedFilter);
    auto highNoiseScore = mean(absMatchedFilter) * Aurora::std(absMatchedFilter);
    // auto highNoiseScore = mean( abs(matchedFilter)) * Aurora::std(absMatchedFilter);
    printf("highNoiseScore[998]: %f4, highNoiseScore[2053]:%f4\r\n", highNoiseScore.getData()[998], highNoiseScore.getData()[2053]);
    auto medianNoise = median(highNoiseScore);
    printf("median: %f , should be: 1724817516.8468074\r\n", medianNoise.getScalar());

    min(abs(highNoiseScore - median(highNoiseScore)), Aurora::Column, row, minCol);
    //
    // minCol = 998;

    auto maxMatchFilter = matchedFilter(Aurora::$, minCol).toMatrix();
    for (int k = 0; k < matchedFilter.getDimSize(1); ++k)
    {
        sumDiff.getData()[k] = sum(abs(matchedFilter(Aurora::$, k).toMatrix() - maxMatchFilter)).getData()[0];
    }

    // printf("meanSumDiff\r\n");
    // auto meanSumDiff = mean(sumDiff);
    // printf("meanSumDiff finish\r\n");
    {
        double meansumDiff = mean(sumDiff).getScalar() ;
        double stdSumDiff =  Aurora::std(sumDiff).getScalar();
        double sumDiffJ =meansumDiff + 2.596 * stdSumDiff;
        // auto indexe = sumDiff > sumDiffJ;
        printf("indexe finish\r\n");

        for (int l = 0; l < sumDiff.getDataSize(); ++l)
        {
            
            if (sumDiff.getData()[l]> sumDiffJ)
            {
                matchedFilter(Aurora::$, l) = matchedFilter(Aurora::$, minCol);
            }
        }
        printf("\r\n");
    }
    

    if (measuredCEused)
    {
        auto mFTime2 = real(ifft(-matchedFilter));
        mFTime2 = mFTime2 / repmat(max(Aurora::abs(mFTime2)), mFTime2.getDimSize(0), 1);
        mFTime2 = mFTime2 / repmat(sum(Aurora::abs(mFTime2)), mFTime2.getDimSize(0), 1);
        matchedFilter = fft(mFTime2);
    }
    for (size_t i = 0; i < 10000; i += 500)
    {
        printf("index :%d,origin output:%f4,%f4, output:%f4,%f4\r\n",
               i, o.mReal[i], o.mImag[i], matchedFilter.getData()[i * 2], matchedFilter.getData()[i * 2 + 1]);
    }

    return 0;
}