/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#ifndef _GCTL_KERNEL_DENSITY_ESTIMATION_H
#define _GCTL_KERNEL_DENSITY_ESTIMATION_H

#include "../core.h"
#include "../maths.h"
#include "../io.h"

namespace gctl
{
    enum kde_kernel_e
    {
        KDE_Gaussian,
        KDE_Epanechnikov,
        KDE_Rectangular,
        KDE_Triangular,
    };

    enum kde_norm_e
    {
        KDE_MAX2ONE,
        KDE_SUM2ONE,
        KDE_CUMSTOM,
    };

    class kde
    {
    public:
        kde();
        kde(double h, const array<double> &x);
        virtual ~kde();

        void init(double h, const array<double> &x);
        double get_density_at(double x, kde_kernel_e k_type = KDE_Gaussian);

        /**
         * @brief Get the probability density of a single kernel. Note the value is not normalized by the kernel number.
         * 
         * @param k_id kernel index
         * @param x inquiring location
         * @param k_type kernel type
         * @return kernel value
         */
        double get_kernel_density_at(size_t k_id, double x, kde_kernel_e k_type = KDE_Gaussian);
        double get_gradient_at(double x, kde_kernel_e k_type = KDE_Gaussian);
        double get_kernel_gradient_at(size_t k_id, double x, kde_kernel_e k_type = KDE_Gaussian);
        void get_distribution(const array<double> x, array<double> &d, array<double> &dx, 
            kde_kernel_e k_type = KDE_Gaussian, kde_norm_e n_type = KDE_MAX2ONE, double norm = 1.0);

    private:
        double h_;
        array<double> x_;

        double gaussian_kernel(double x);
        double epanechnikov_kernel(double x, bool gradient = false);
        double rectangular_kernel(double x, bool gradient = false);
        double triangular_kernel(double x, bool gradient = false);
    };

    class kde2d
    {
    public:
        kde2d();
        kde2d(double h, const array<double> &x, const array<double> &y);
        kde2d(double h, const std::vector<double> &x, const std::vector<double> &y);
        virtual ~kde2d();

        void init(double h, const array<double> &x, const array<double> &y);
        void init(double h, const std::vector<double> &x, const std::vector<double> &y);
        double get_density_at(double x, double y, kde_kernel_e k_type = KDE_Gaussian);

        /**
         * @brief Get the probability density of a single kernel. Note the value is not normalized by the kernel number.
         * 
         * @param k_id kernel index
         * @param x inquiring location x
         * @param y inquiring location y
         * @param k_type kernel type
         * @return kernel value
         */
        double get_kernel_density_at(size_t k_id, double x, double y, kde_kernel_e k_type = KDE_Gaussian);
        void get_gradient_at(double x, double y, double &gx, double &gy, kde_kernel_e k_type = KDE_Gaussian);
        void get_kernel_gradient_at(size_t k_id, double x, double y, double &gx, double &gy, kde_kernel_e k_type = KDE_Gaussian);
        void get_distribution(const array<double> x, const array<double> y, array<double> &d, array<double> &dx, 
            array<double> &dy, kde_kernel_e k_type = KDE_Gaussian, kde_norm_e n_type = KDE_MAX2ONE, double norm = 1.0);
        void save(double xmin, double xmax, double ymin, double ymax, int xnum, int ynum, std::string file);
    
    private:
        double h_;
        array<double> x_, y_;

        double gaussian_kernel(double x, double y);
    };
}

#endif // _GCTL_KERNEL_DENSITY_ESTIMATION_H