/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */

/* 
 * File:   Spectrum.hh
 * Author: hacene
 *
 * Created on December 31, 2020, 10:14 AM
 */

#ifndef SPECTRUM_HH
#define SPECTRUM_HH

#include "Parameter.hh"
#include "ParamData.hh"
#include "DataTypeMath.hh"
#include "Operation.hh"
#include "AMDA_exception.hh"


#include <iterator>
namespace AMDA {
    namespace Parameters {
        namespace Spectrum {

#define AVERAGE_TIME 1200 // (seconds)
#define MAX_GAP_SIZE 3600 // (seconds)

            using namespace std;

            typedef enum {
                CENTER,
                FORWARD,
                BACKWARD,
            } samplingEnum;

            static std::map<std::string, samplingEnum> strToSamplingEnum{

                {"center", samplingEnum::CENTER},
                {"forward", samplingEnum::FORWARD},
                {"backward", samplingEnum::BACKWARD},};

            template<typename DataType, typename OutputType>
            class SpectrumBase : public Operation {
            public:

                /**
                 * @brief Constructor.
                 * @details Create the ParamData type of the input ParamData.
                 */
                SpectrumBase(Process& pProcess, TimeIntervalListSPtr pTimeIntervalList, ParamDataSpec<DataType >&paramInput, double windowtime, std::string sampling = "center") :
                Operation(pProcess),
                _timeIntervalList(pTimeIntervalList),
                _currentTimeInterval(pTimeIntervalList->begin()),
                _paramInput(paramInput),
                _paramOutput(new ParamDataSpec<std::vector<OutputType>>()),
                _windowtime(windowtime),
                _sampling(sampling),
                _needInit(true) {
                    _paramDataOutput = _paramOutput;
                };

                virtual ~SpectrumBase() {
                }

                double getWindowTime() {
                    return _windowtime;
                }

                bool inWindow(double time) {
                    return (time >= _minWindow && time <= _maxWindow);
                }

                double getIntStartTime() {
                    return _currentTimeInterval->_startTime;
                }

                double getIntStopTime() {
                    return _currentTimeInterval->_stopTime;
                }

                bool inInt(double time) {
                    return (time >= _currentTimeInterval->_startTime && time <= _currentTimeInterval->_stopTime);
                }

                double getTarget() {
                    return _target;
                }

                bool needInit() {
                    return _needInit;
                }

                void setNeedInit(bool needInit) {
                    _needInit = needInit;
                }

                bool setTarget(double time) {

                    if (!inInt(time)) {
                        return false;
                    }

                    _target = time;

                    switch (strToSamplingEnum[_sampling]) {
                        case samplingEnum::FORWARD:
                            _minWindow = _target;
                            if (_minWindow < _currentTimeInterval->_startTime) {
                                _minWindow = _currentTimeInterval->_startTime;
                            }
                            _maxWindow = _target + _windowtime;
                            if (_maxWindow > _currentTimeInterval->_stopTime) {
                                _maxWindow = _currentTimeInterval->_stopTime;
                            }
                        case samplingEnum::BACKWARD:
                            _minWindow = _target - _windowtime;
                            if (_minWindow < _currentTimeInterval->_startTime) {
                                _minWindow = _currentTimeInterval->_startTime;
                            }
                            _maxWindow = _target;
                            if (_maxWindow > _currentTimeInterval->_stopTime) {
                                _maxWindow = _currentTimeInterval->_stopTime;
                            }
                            break;
                        default:
                            _minWindow = _target - _windowtime / 2.;
                            if (_minWindow < _currentTimeInterval->_startTime) {
                                _minWindow = _currentTimeInterval->_startTime;
                            }
                            _maxWindow = _target + _windowtime / 2.;
                            if (_maxWindow > _currentTimeInterval->_stopTime) {
                                _maxWindow = _currentTimeInterval->_stopTime;
                            }
                            break;
                    }
                    return true;
                }

                /**
                 * @overload Operation::reset(double pStartTime, double pTimeInt)
                 * @brief reset static data to process another TimeInterval
                 */
                virtual void reset() {
                    Operation::reset();
                    if (_currentTimeInterval == _timeIntervalList->end()) {
                        return;
                    }
                    ++_currentTimeInterval;
                    _needInit = true;
                    resetFunc();
                }

                virtual void init() {
                    setTarget(getIntStartTime());
                    setNeedInit(false);
                }

                virtual bool nextTarget() {
                    double target = getTarget() + getWindowTime();
                    bool res = setTarget(target);
                    while (!_mem.empty() && !inWindow(_mem.front().first)) {
                        _mem.pop_front();
                    }
                    return res;
                }

                virtual bool needToChangeTarget(double crtTime) {
                    return !needInit() && !inWindow(crtTime);
                }

                virtual double getSampling() {
                    return getWindowTime();
                }

                virtual void pushData(double time, DataType& elem) {
                    _mem.push_back(std::make_pair(time, elem));
                }

                virtual void resetFunc() {
                    _mem.clear();
                }

                double getInputParamSampling() {
                    return _paramInput.getMinSampling();
                }

                virtual std::vector<OutputType> compute() = 0;

                /**
                 * @overload Operation::write(ParamDataIndexInfo &pParamDataIndexInfo)
                 */
                void write(ParamDataIndexInfo &pParamDataIndexInfo) {
                    if ((pParamDataIndexInfo._nbDataToProcess > 0)) {
                        for (unsigned int _index = pParamDataIndexInfo._startIndex;
                                _index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess;
                                ++_index) {
                            double crtTime = _paramInput.getTime(_index);
                            DataType crtVal = _paramInput.get(_index);
                            if (needToChangeTarget(crtTime)) {
                                _paramOutput->pushTime(getTarget());
                                _paramOutput->push(compute());
                                pushData(crtTime, crtVal);
                                nextTarget();
                                bool skip = false;
                                while (!skip && needToChangeTarget(crtTime)) {
                                    _paramOutput->pushTime(getTarget());
                                    _paramOutput->push(compute());
                                    skip = nextTarget();
                                }
                            } else {
                                pushData(crtTime, crtVal);
                                if (needInit()) {
                                    init();
                                }
                            }
                        }
                    }
                    if (pParamDataIndexInfo._timeIntToProcessChanged || pParamDataIndexInfo._noMoreTimeInt) {
                        if (!needInit()) {
                            do {
                                if (inInt(getTarget())) {
                                    _paramOutput->pushTime(getTarget());
                                    _paramOutput->push(compute());
                                }
                            } while (nextTarget());
                        }
                    }
                }
            protected:

                TimeIntervalListSPtr _timeIntervalList;

                TimeIntervalList::iterator _currentTimeInterval;

                ParamDataSpec<DataType>& _paramInput;

                ParamDataSpec<std::vector<OutputType>>*_paramOutput;

                double _windowtime;

                std::string _sampling;

                double _minWindow;

                double _maxWindow;

                double _target;

                bool _needInit;

                std::list<std::pair<double, DataType> > _mem;

            };

            template<typename DataType, typename OutputType>
            class SpectrumFourier : public SpectrumBase<DataType, OutputType> {
            public:

                SpectrumFourier(Process& pProcess, TimeIntervalListSPtr pTimeIntervalList, ParamDataSpec<DataType >&paramInput, double windowtime, std::string sampling = "center") :
                SpectrumBase<DataType, OutputType>(pProcess, pTimeIntervalList, paramInput, windowtime, sampling = "center") {
                    _spectroSize = std::ceil(windowtime / paramInput.getMinSampling()) - 1;
                }

                virtual ~SpectrumFourier() {
                }

                virtual void init() {
                    if (!_targets.empty()) {
                        SpectrumBase<DataType, OutputType>::setTarget(_targets.front());
                        SpectrumBase<DataType, OutputType>::setNeedInit(false);
                        _targets.pop_front();
                    }
                }

                virtual bool nextTarget() {
                    if (!_targets.empty()) {
                        bool res = SpectrumBase<DataType, OutputType>::setTarget(_targets.front());
                        _targets.pop_front();
                        while (!_mem.empty() && !SpectrumBase<DataType, OutputType>::inWindow(_mem.front().first)) {
                            _mem.pop_front();
                        }
                        return res;
                    }
                    return false;
                }

                virtual bool needToChangeTarget(double crtTime) {
                    return !SpectrumBase<DataType, OutputType>::needInit() && !SpectrumBase<DataType, OutputType>::inWindow(crtTime) && !_targets.empty();
                }

                virtual double getSampling() {
                    return SpectrumBase<DataType, OutputType>::getInputParamSampling();
                }

                virtual void pushData(double time, DataType& elem) {
                    _mem.push_back(std::make_pair(time, elem));
                    _targets.push_back(time);
                }

                virtual void resetFunc() {
                    _mem.clear();
                    _targets.clear();
                }

                bool isUniform() {
                    if (_mem.empty())
                        return false;

                    auto it_1 = _mem.begin();
                    it_1++;
                    for (auto it = _mem.begin(); it != _mem.end(); it++) {
                        if (*it != _mem.back()) {
                            if (it_1->first - it->first != getSampling())
                                return false;
                            it_1++;
                        }
                    }
                    return true;
                }

                std::vector<OutputType> fft(std::vector<DataType> &xn, int N) {
                    /**
                     * 
                     * @param N number of points in the dft
                     * @param xn date vector 
                     */

                    int len = xn.size();
                    std::vector<OutputType> output;
                    output.resize(N);
                    output << NotANumber();

                    double Xr;
                    double Xi;
                    int counter;

                    int k, n = 0;
                    for (k = 0; k < N; k++) {
                        Xr = 0;
                        Xi = 0;
                        counter = 0;
                        for (n = 0; n < len; n++) {
                            if (!isNAN(xn[n])) {
                                Xr = Xr + (OutputType) xn[n] * cos(2 * 3.141592 * k * n / N);
                                Xi = Xi + (OutputType) xn[n] * sin(2 * 3.141592 * k * n / N);
                                counter += 1;
                            }
                        }
                        if (counter == 0)
                            break;
                        output[k] = std::sqrt(Xr * Xr + Xi * Xi);
                    }
                    return output;
                }

                DataType computeInterpolation(std::list<std::pair<double, DataType> > & input, double time) {
                    double min_t = time - AVERAGE_TIME / 2.;
                    double max_t = time + AVERAGE_TIME / 2.;
                    std::vector<std::pair<double, DataType> > values_for_mean;
                    DataType nanVal;
                    nanVal << NotANumber();
                    std::pair<double, DataType> prev_value(NAN, nanVal);
                    std::pair<double, DataType> next_value(NAN, nanVal);
                    DataType value = nanVal;
                    for (auto it = input.begin(); it != input.end(); ++it) {
                        if (it->first == time) {
                            value = it->second;
                            return value;
                            break;
                        } else if (isNAN(it->second))
                            continue;
                        else if (it->first > max_t) {
                            next_value = *it;
                            break;
                        } else if (it->first < min_t) {
                            prev_value = *it;
                        } else {
                            values_for_mean.push_back(*it);
                        }
                    }
                    if (!values_for_mean.empty()) {
                        //Compute mean
                        DataType sum = 0;
                        for (auto it = values_for_mean.begin(); it != values_for_mean.end(); ++it) {
                            sum += it->second;
                        }
                        value = sum / (DataType) values_for_mean.size();
                    } else {
                        if (!isNAN(prev_value.first) && !isNAN(next_value.first) && (next_value.first - prev_value.first <= MAX_GAP_SIZE)) {
                            //Compute interpolated value
                            value = prev_value.second + (time - prev_value.first) / (next_value.first - prev_value.first) * (next_value.second - prev_value.second);
                        }
                    }

                    return value;

                }

                void standarize(std::list<std::pair<double, DataType> > & mem, int size, double target, double inputSampling, std::string sampling) {
                    auto it = mem.begin();
                    double time_;
                    std::list<std::pair<double, DataType>> valuesToAdd;
                    for (int i = 0; i < size; i++) {
                        switch (strToSamplingEnum[sampling]) {
                            case samplingEnum::FORWARD:
                                time_ = target + i*inputSampling;
                                break;
                            case samplingEnum::BACKWARD:
                                time_ = target + (i - size + 1) * inputSampling;
                                break;
                            default:
                                time_ = target + (i - size / 2 + 1) * inputSampling;
                                break;
                        }
                        if (i < mem.size()) {
                            if (it->first != time_ || isNAN(it->second)) {
                                std::pair<double, DataType> valueToAdd(time_, computeInterpolation(mem, time_));
                                valuesToAdd.push_back(valueToAdd);
                            }
                        } else {
                            std::pair<double, DataType> valueToAdd(time_, computeInterpolation(mem, time_));
                            valuesToAdd.push_back(valueToAdd);
                        }
                        it++;
                    }
                    if (valuesToAdd.empty() || valuesToAdd.size() < size - mem.size())
                        BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_PROCESS_ERR) << AMDA::ex_msg(std::string("Spectrum:: ERROR at data interpolation")));

                    mem.merge(valuesToAdd);
                    return;
                }
                
                int numberOfNanValues(std::list<std::pair<double, DataType> > & mem){
                    int res = 0;
                     for (auto elt : mem){
                         if(isNAN(elt.second))
                             res ++;
                     }
                    return res;
                }

                std::vector<OutputType> compute() {
                    std::vector<OutputType> res;
                    res.resize(_spectroSize);
                    res << NotANumber();
                    bool is_uniform = isUniform();
                    int nans =  numberOfNanValues(_mem);
                    if (!is_uniform || _mem.empty() || _mem.size() < 2*_spectroSize/3 || nans < 2*_spectroSize/3)
                        return res;

                    if (_mem.size() < _spectroSize)
                        standarize(_mem, _spectroSize, SpectrumBase<DataType, OutputType>::_target, SpectrumBase<DataType, OutputType>::getInputParamSampling(), SpectrumBase<DataType, OutputType>::_sampling);

                    std::vector<DataType> inputData;
                    
                    for (auto elt : _mem)
                        inputData.push_back(elt.second);
                    
                    res = fft(inputData, _spectroSize);
                    
                    return res;
                }

            protected:
                std::list<double> _targets;

                std::list<std::pair<double, DataType> > _mem;

                int _spectroSize;
            };




        } // end Spectrum
    } // end Parameters
} // end AMDA


#endif /* SPECTRUM_HH */