Kernels.hpp

#pragma once
/* (C) 2021 Roman Werpachowski. */
#include <memory>
#include <type_traits>
#include <Eigen/Core>
#include "dll.hpp"
 
namespace ml
{
    namespace Kernels
    {
        class Kernel
        {
        public:
            DLL_DECLSPEC virtual ~Kernel();
 
            DLL_DECLSPEC virtual double value(const Eigen::Ref<const Eigen::VectorXd> x1, const Eigen::Ref<const Eigen::VectorXd> x2) const = 0;
 
            virtual unsigned int dim() const = 0;
        protected:
            DLL_DECLSPEC void validate_arguments(const Eigen::Ref<const Eigen::VectorXd> x1, const Eigen::Ref<const Eigen::VectorXd> x2) const;
        };
 
        class DifferentiableKernel: public virtual Kernel
        {
            DLL_DECLSPEC virtual void gradient(const Eigen::Ref<const Eigen::VectorXd> x1, const Eigen::Ref<const Eigen::VectorXd> x2, Eigen::Ref<Eigen::VectorXd> dydx1) const = 0;
        };
 
        
        class DoubleDifferentiableKernel : public virtual DifferentiableKernel
        {
            DLL_DECLSPEC virtual void hessian(const Eigen::Ref<const Eigen::VectorXd> x1, const Eigen::Ref<const Eigen::VectorXd> x2, Eigen::Ref<Eigen::MatrixXd> H) const = 0;
        };
 
        class RadialBasisFunction
        {
        public:
            DLL_DECLSPEC virtual ~RadialBasisFunction();
 
            DLL_DECLSPEC virtual double value(double r2) const = 0;
        };
 
        class DifferentiableRadialBasisFunction : public RadialBasisFunction
        {
        public:
            DLL_DECLSPEC virtual double gradient(double r2) const = 0;
        };
 
        class DoubleDifferentiableRadialBasisFunction : public DifferentiableRadialBasisFunction
        {
        public:
            DLL_DECLSPEC virtual double second_derivative(double r2) const = 0;
        };
 
        class GaussianRBF: public DoubleDifferentiableRadialBasisFunction
        {
        public:
            DLL_DECLSPEC double value(double r2) const override;
            DLL_DECLSPEC double gradient(double r2) const override;
            DLL_DECLSPEC double second_derivative(double r2) const override;
        };
 
        template <class RBF = RadialBasisFunction> class RBFKernel : public virtual Kernel
        {
            static_assert(std::is_base_of_v<RadialBasisFunction, RBF>);
        public:
            RBFKernel(std::unique_ptr<const RBF>&& rbf, unsigned int dim)
                : rbf_(std::move(rbf)), dim_(dim)
            {
                if (rbf_ == nullptr) {
                    throw std::invalid_argument("Null RBF object");
                }
                if (!dim) {
                    throw std::domain_error("Kernel dimension must be positive");
                }
            }
 
            RBFKernel(const RBFKernel& other) = delete;
            RBFKernel& operator=(const RBFKernel& other) = delete;
 
            DLL_DECLSPEC RBFKernel(RBFKernel&& other) = default;
            DLL_DECLSPEC RBFKernel& operator=(RBFKernel&& other) = default;
 
            double value(const Eigen::Ref<const Eigen::VectorXd> x1, const Eigen::Ref<const Eigen::VectorXd> x2) const override
            {
                validate_arguments(x1, x2);
                return rbf_->value((x1 - x2).squaredNorm());
            }
 
            unsigned int dim() const override
            {
                return dim_;
            }
        protected:
            std::unique_ptr<const RBF> rbf_;
        private:
            unsigned int dim_;
        };
 
        template <class DiffRBF = DifferentiableRadialBasisFunction> class DifferentiableRBFKernel : public virtual DifferentiableKernel, public RBFKernel<DiffRBF>
        {
            static_assert(std::is_base_of_v<DifferentiableRadialBasisFunction, DiffRBF>);
        public:
            using RBFKernel<DiffRBF>::RBFKernel;
 
            void gradient(const Eigen::Ref<const Eigen::VectorXd> x1, const Eigen::Ref<const Eigen::VectorXd> x2, Eigen::Ref<Eigen::VectorXd> dydx1) const override
            {
                this->validate_arguments(x1, x2);
                if (dydx1.size() != this->dim()) {
                    throw std::invalid_argument("Wrong dimension of dydx1");
                }
                const double rbf1der = this->rbf_->gradient((x1 - x2).squaredNorm());
                for (Eigen::Index i = 0; i < x1.size(); ++i) {
                    dydx1[i] = 2 * (x1[i] - x2[i]) * rbf1der;
                }
            }
        };        
    }
}