#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 "spdlog/spdlog.h"
#include "spdlog/sinks/stdout_color_sinks.h"
#include "spdlog/sinks/basic_file_sink.h"
#include "MatlabReader.h"

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

    m.read(&i,&o,"/mnt/d/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 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: %f4\r\n",medianNoise.getScalar());
    long minCol,row;
    min(abs(highNoiseScore - median(highNoiseScore)), Aurora::Column, row, minCol);
    //
    minCol = 998;
    auto sumDiff = Aurora::zeros(1, matchedFilter.getDimSize(1));
    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];
    }

    auto indexe = sumDiff > (mean(sumDiff)+2.596 * Aurora::std(sumDiff));
    for (int l = 0; l < indexe.getDataSize(); ++l) {
        if ((bool)indexe.getData()[l]){
            matchedFilter(Aurora::$,l) = matchedFilter(Aurora::$,minCol);
        }
    }
    if (measuredCEused){
        auto mFTime2 = ifft(-matchedFilter);
        max(Aurora::abs(mFTime2)).printfShape();
        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);
    }
    return 0;
}