#pragma once

#ifndef clip_h
#define clip_h

#include "utils.h"

#include <Eigen/Dense>

class SigmaClip {
private:
    int maxiters = 5;
    double sigma = 3, minValue = INF, maxValue = -INF;
    double (*cenF)(std::vector<Eigen::Vector2d> data) = nullptr;
    double (*stdF)(std::vector<Eigen::Vector2d> data) = nullptr;

public:
    SigmaClip(double sigma = 3, int maxiters = 5,
              double (*cenF)(std::vector<Eigen::Vector2d> data) = nullptr,
              double (*stdF)(std::vector<Eigen::Vector2d> data) = nullptr);
    ~SigmaClip() {};
    std::vector<Eigen::Vector2d> clip(std::vector<Eigen::Vector2d> data);
    std::vector<Eigen::Vector2d> clipN(std::vector<Eigen::Vector2d> data);

private:
    void computeBound(std::vector<Eigen::Vector2d> data);
};

SigmaClip::SigmaClip(double sigma_, int maxiters_,
                     double (*cenF_)(std::vector<Eigen::Vector2d> data),
                     double (*stdF_)(std::vector<Eigen::Vector2d> data)) {
    sigma = sigma_;
    maxiters = maxiters_;
    cenF = cenF_ == nullptr ? dataAvg2D : cenF_;
    stdF = stdF_ == nullptr ? dataStd2D : stdF_;
}

void SigmaClip::computeBound(std::vector<Eigen::Vector2d> data) {
    double std = (*stdF)(data);
    double mean = (*cenF)(data);

    minValue = mean - (std * sigma);
    maxValue = mean + (std * sigma);
}

std::vector<Eigen::Vector2d> SigmaClip::clip(std::vector<Eigen::Vector2d> data) {
    std::vector<Eigen::Vector2d>::iterator itor;
    
    minValue = INF, maxValue = -INF;
    for (int k = 1; k <= maxiters; k++) {
        computeBound(data);
        for (itor = data.begin(); itor != data.end();) {
            if ((*itor)(0) < minValue || (*itor)(0) > maxValue)
                data.erase(itor);
            else
                itor++;
        }
    }

    return data;
}

#endif