#include <iostream>
#include <complex>
#include <functional>


class Kernel {
public: 
    /*
    Covariance function : return the covariance between two R.V. for the entire kernel's domain definition. 
    */
    virtual double covarianceFunction(
        double   X,
        double   Y
    )const = 0 ;
    ~Kernel() = default;
};


template <typename Func>
class FooKernel : public Kernel {
public:

    FooKernel(Func&& fun) : fun_(std::forward<Func>(fun)) {}
    double covarianceFunction(
        double   X,
        double   Y
    ) const {
        return fun_(X, Y);
    }
    template<class T>
    auto operator+(const T b) const {
        return make_foo_kernel([b, this](double X, double Y) -> double {
            return this->covarianceFunction(X, Y) + b.covarianceFunction(X, Y);
        });
    }
    FooKernel operator=(const FooKernel other) const {
        return other;
    }
private:
   Func fun_;
};

template <typename Func>
auto make_foo_kernel(Func&& fun)
{
    return FooKernel<Func>(std::forward<Func>(fun));
}


class GaussianKernel : public Kernel {
public:
    GaussianKernel(double sigma, double scale) : m_sigma(sigma), m_scale(scale) {}
    GaussianKernel(double sigma) : m_sigma(sigma), m_scale(1) {}
    /*
    A well known covariance function that enforces smooth deformations
    Ref : Shape modeling using Gaussian process Morphable Models, Luethi et al.
    */
    double covarianceFunction(
        double   X,
        double   Y
    ) const 
    {
        //use diagonal matrix
    double result;
    result = m_scale  *  exp(-std::norm(X - Y) / (m_sigma*m_sigma));
    return result;      
    }
    template<class T>
    auto operator+(const T b) const {
        return make_foo_kernel([b, this](double X, double Y) -> double {
            auto debugBval = b.covarianceFunction(X, Y);
            auto debugAval = this->covarianceFunction(X, Y);
            auto test = debugBval + debugAval;
            return test;
        });
    }
private:
    double m_sigma;
    double m_scale;
};

int main()
{
    auto C = GaussianKernel(50,60) + GaussianKernel(100,200);
    auto result = C.covarianceFunction(30.0,40.0);

    return 0;
}