#ifndef kernels_hpp
#define kernels_hpp
#include <vector>
#include <functional>
#include <type_traits>
#include <cmath>
#include "NestedMatrix.hpp"
// K:R(n) X R(n) -> R (or complex)
template <typename T>
using VectorType = std::vector<T>;
template <typename T>
using MatrixType = NestedMatrix<T>;
// Kernels with vector arguments
template <typename T>
using KernelType = std::function<T(const VectorType<T>& a, const VectorType<T>& b)>;
template <typename T1, typename T2, template <typename T2> class S >
using KernelTypeGeneral = std::function<T1 (const S<T2>& a, const S<T2>& b)>;
template <typename T1, typename T2>
using KernelVectorType = KernelTypeGeneral<T1, T2, VectorType>;
template <typename T1, typename T2>
using KernelMatrixType = KernelTypeGeneral<T1, T2, MatrixType>;
// Evaluate a kernel function k(x,y)
// Uses concepts + template template parametsrs
template <typename T, template <typename T> class K >
concept KernelClient = requires (K<T> k, const VectorType<T>&a, const VectorType<T>&b)
{ // Evaluate kernel at two vector points a and b
k(a, b);
};
// Evaluate a similarity function for kernels
template <typename T, template <typename T> class K >
concept KernelClientII = requires (K<T> k, MatrixType<T>& p, MatrixType<T>& q)
{ // Evaluate similarity function for two point sets p and q
k.similarity(p, q);
};
// Evaluate a distance function for kernels
template <typename T, template <typename T> class K >
concept KernelDistance = requires (K<T> k, MatrixType<T>& p, MatrixType<T>& q)
{ // Evaluate similarity function for two point sets p and q
k.distance(p, q);
};
template <typename T, template <typename T> class K>
concept KernelProtocol = KernelClient<T, K> && KernelClientII<T, K> && KernelDistance<T, K>;
// Kernel classes
template <typename T, typename D>
class BaseKernel
{ // Using CRTP and Template Method Pattern
private:
public:
BaseKernel() {}
// K(x,y) for derived classes
T operator() (const VectorType<T>& x, const VectorType<T>& y)
{
return static_cast<D*>(this)->kernel(x,y);
}
T kernel(const VectorType<T>& x, const VectorType<T>& y)
{
return (*this)(x, y);
}
NestedMatrix<T> kernelMatrix(const NestedMatrix<T>& x, const NestedMatrix<T>& y)
{ // x[0], x[1], ...,x[n] and y[0], y[1], ...,y[m] are vectors of size D
std::size_t n = x.size1();
std::size_t m = y.size1();
NestedMatrix<T> mat(n, m);
for (std::size_t i = 0; i < n; ++i)
{
for (std::size_t j = 0; j < m; ++j)
{
mat(i, j) = kernel(x.row(i), y.row(j));
}
}
return mat;
}
NestedMatrix<T> distanceMatrix(const NestedMatrix<T>& x, const NestedMatrix<T>& y)
{ // x[0], x[1], ...,x[n] and y[0], y[1], ...,y[m] are vectors of size D
std::size_t n = x.size1();
std::size_t m = y.size1();
NestedMatrix<T> mat(n, m);
for (std::size_t i = 0; i < n; ++i)
{
for (std::size_t j = 0; j < m; ++j)
{
mat(i, j) = kernel(x.row(i), y.row(j));
}
}
return mat;
}
T similarity (MatrixType<T>& p, MatrixType<T>& q)
{
T result = 0.0;
for (std::size_t i = 0; i < p.size1(); ++i)
{
for (std::size_t j = 0; j < q.size1(); ++j)
{
// Kick down to the derived class
result += static_cast<D*>(this)->kernel(p.row(i), q.row(j));
}
}
return result;
}
T distance (MatrixType<T>& p, MatrixType<T>& q)
{ // aka sqrt(discrepancy error)
return std::sqrt(similarity(p, p) + similarity(q, q) - 2.0* similarity(p, q));
}
T distance(const VectorType<T>& p, const VectorType<T>& q)
{
return std::sqrt(kernel(p, p) + kernel(q, q) - 2.0 * kernel(p, q));
}
NestedMatrix<T> distaneMatrix(const NestedMatrix<T>& x, const NestedMatrix<T>& y)
{ // x[0], x[1], ...,x[n] and y[0], y[1], ...,y[m] are vectors of size D
std::size_t n = x.size1();
std::size_t m = y.size1();
NestedMatrix<T> mat(n, m);
for (std::size_t i = 0; i < n; ++i)
{
for (std::size_t j = 0; j < m; ++j)
{
mat(i, j) = distance(x.row(i), y.row(j));
}
}
return mat;
}
};
template <typename T>
class LinearKernel : public BaseKernel<T, LinearKernel<T>>
{
private:
T c;
public:
LinearKernel(T offset = 0.0) : c(offset) {}
// K(x,y)
T operator() (const VectorType<T>& x, const VectorType<T>& y)
{
T result = x[0] * y[0] + c;
for (std::size_t i = 1; i < x.size(); ++i)
{
result += x[i] * y[i];
}
return result;
}
T kernel(const VectorType<T>& x, const VectorType<T>& y)
{
return (*this)(x, y);
}
};
template <typename T>
class RELUKernel : public BaseKernel<T, LinearKernel<T>>
{
private:
T c;
public:
RELUKernel(T offset = 0.0) : c(offset) {}
// K(x,y)
T operator() (const VectorType<T>& x, const VectorType<T>& y)
{
T result = x[0] * y[0] + c;
for (std::size_t i = 1; i < x.size(); ++i)
{
result += std::abs(x[i] - y[i]);
}
return 1.0 - result;
}
T kernel(const VectorType<T>& x, const VectorType<T>& y)
{
return (*this)(x, y);
}
};
template <typename T>
class MaternKernel : public BaseKernel<T, LinearKernel<T>>
{
private:
T c;
public:
MaternKernel(T offset = 0.0) : c(offset) {}
// K(x,y)
T operator() (const VectorType<T>& x, const VectorType<T>& y)
{
T result = x[0] * y[0] + c;
for (std::size_t i = 1; i < x.size(); ++i)
{
result += std::abs(x[i] - y[i]);
}
return 1.0 - result;
}
T kernel(const VectorType<T>& x, const VectorType<T>& y)
{
return (*this)(x, y);
}
};
template <typename T>
class GaussianKernel : public BaseKernel<T, LinearKernel<T>>
{
private:
T sig;
public:
GaussianKernel(T sigma) : sig(sigma) {}
// K(x,y)
T kernel (const VectorType<T>& x, const VectorType<T>& y) const
{
T diff = 0.0;
T tmp;
for (std::size_t i = 1; i < x.size(); ++i)
{
tmp = x[i] - y[i];
diff += tmp*tmp;
}
return std::exp(-diff / (2.0 * sig * sig));
}
};
template <typename T>
class LaplacianKernel : public BaseKernel<T, LinearKernel<T>>
{
private:
T sig;
public:
LaplacianKernel(T sigma) : sig(sigma) {}
// K(x,y)
T kernel(const VectorType<T>& x, const VectorType<T>& y) const
{
T tmp = 0.0;
for (std::size_t i = 1; i < x.size(); ++i)
{
tmp += std::abs(x[i] - y[i]);
}
return std::exp(-tmp / sig);
}
};
template <typename T>
class MLPKernel : public BaseKernel<T, MLPKernel<T>>
{ //A multilayer perceptron(MLP) is a class of feedforward artificial neural network(ANN).
private:
// K(x,y) = tanh(a x*y + c)
T c;
T a;
public:
MLPKernel(T intercept, T slope) : a(slope), c (intercept) {}
// K(x,y) = tanh(
T kernel (const VectorType<T>& x, const VectorType<T>& y) const
{
T sum = 0.0;
for (std::size_t i = 0; i < x.size(); ++i)
{
sum += a*x[i]*y[i] + c;
}
return std::tanh(sum + c);
}
};
template <typename T>
class PolynomialKernel : public BaseKernel<T, MLPKernel<T>>
{ //A multilayer perceptron(MLP) is a class of feedforward artificial neural network(ANN).
private:
// K(x,y) = (a x*y + c)^d
T c;
T a;
T d;
public:
PolynomialKernel(T slope, T intercept, T degree) : a(slope), c(intercept), d(degree) {}
// K(x,y) =
T kernel (const VectorType<T>& x, const VectorType<T>& y) const
{
T sum = 0.0;
for (std::size_t i = 0; i < x.size(); ++i)
{
sum += a * x[i] * y[i] + c;
}
return std::pow(sum + c, d);
}
};
template <typename T>
class ANOVAKernel : public BaseKernel<T, LinearKernel<T>>
{ // radial basis function, useful in multidimensional regression prob
private:
T sig;
T d;
public:
ANOVAKernel(T sigma, T degree) : sig(sigma), d(degree) {}
// K(x,y)
T kernel(const VectorType<T>& x, const VectorType<T>& y) const
{
T sum = 0.0;
T tmp;
for (std::size_t i = 1; i < x.size(); ++i)
{
tmp = x[i] - y[i];
sum += std::exp(std::pow(- sig*tmp*tmp,d));
}
return sum;
}
};
#endif
