diff --git a/docs/src/kernels.md b/docs/src/kernels.md
index a0fde4d..e6d75e1 100644
--- a/docs/src/kernels.md
+++ b/docs/src/kernels.md
@@ -1,26 +1,198 @@
 # Kernels
 
-## Mercer Kernels
+| Kernel | Mercer | Negative Definite | Stationary | Isotropic |
+| --- | :-: | :-: | :-: | :-: |
+| [Exponential Kernel](#Exponential-Kernel-1) | ✓ | | ✓ | ✓ |
+| [Rational Quadratic Kernel](#Rational-Quadratic-Kernel-1) | ✓ | | ✓ | ✓ |
+| [Exponentiated Kernel](#Exponentiated-Kernel-1) | ✓ | | | |
+
+## Exponential Kernel
+
+**Exponential Kernel**
+
+The exponential kernel (see [`ExponentialKernel`](@ref)) is an isotropic Mercer kernel of
+the form:
+
+```math
+\kappa(\mathbf{x},\mathbf{y})
+= \exp\left(-\alpha ||\mathbf{x} - \mathbf{y}||\right)
+\qquad \alpha > 0
+```
+where ``\alpha`` is a positive scaling parameter of the Euclidean distance. This kernel may
+also be referred to as the Laplacian kernel (see [`LaplacianKernel`](@ref)).
+
+**Squared-Exponential Kernel**
+
+A similar form of the exponential kernel squares the Euclidean distance:
+
+```math
+\kappa(\mathbf{x},\mathbf{y})
+= \exp\left(-\alpha ||\mathbf{x} - \mathbf{y}||^2\right)
+\qquad \alpha > 0
+```
+In this case, the kernel is often referred to as the squared exponential kernel (see
+[`SquaredExponentialKernel`](@ref)) or the Gaussian kernel (see [`GaussianKernel`](@ref)).
+
+**``\gamma``-Exponential Kernel**
+
+Both the exponential and the squared exponential kernels are specific cases of the more
+general ``\gamma``-exponential kernel:
+
+```math
+\kappa(\mathbf{x},\mathbf{y})
+= \exp\left(-\alpha ||\mathbf{x} - \mathbf{y}||^{2\gamma}\right)
+\qquad \alpha > 0, \;\; 0 < \gamma \leq 1
+```
+where ``\gamma`` is an additional shape parameter of the Euclidean distance.
+
+### Interface
+
 ```@docs
 ExponentialKernel
 SquaredExponentialKernel
 GammaExponentialKernel
+LaplacianKernel
+GaussianKernel
+RadialBasisKernel
+```
+
+## Rational-Quadratic Kernel
+
+**Rational-Quadratic Kernel**
+
+The rational-quadratic kernel (see [`RationalQuadraticKernel`](@ref)) is an isotropic
+Mercer kernel given by the formula:
+
+```math
+\kappa(\mathbf{x},\mathbf{y})
+= \left(1 +\alpha ||\mathbf{x} - \mathbf{y}||^{2}\right)^{-\beta}
+\qquad \alpha > 0, \;\; \beta > 0
+```
+where ``\alpha`` is a positive scaling parameter and ``\beta`` is a shape parameter of the
+Euclidean distance.
+
+**``\gamma``-Rational-Quadratic Kernel**
+
+The rational-quadratic kernel is a special case with ``\gamma = 1`` of the more general
+``\gamma``-rational-quadratic kernel (see [`GammaRationalQuadraticKernel`](@ref)):
+
+```math
+\kappa(\mathbf{x},\mathbf{y})
+= \left(1 +\alpha ||\mathbf{x} - \mathbf{y}||^{2\gamma}\right)^{-\beta}
+\qquad \alpha > 0, \; \beta > 0, \; 0 < \gamma \leq 1
+```
+where ``\alpha`` is a positive scaling parameter, ``\beta`` is a positive shape parameter
+and ``\gamma`` is a shape parameter of the Euclidean distance.
+
+### Interface
+```@docs
 RationalQuadraticKernel
-GammaRationalKernel
+GammaRationalQuadraticKernel
+```
+
+## Exponentiated Kernel
+
+The exponentiated kernel (see [`ExponentiatedKernel`](@ref)) is a Mercer kernel given by:
+
+```math
+\kappa(\mathbf{x},\mathbf{y}) = \exp\left(a \mathbf{x}^\intercal \mathbf{y} \right)
+\qquad a > 0
+```
+
+where ``\alpha`` is a positive shape parameter.
+
+### Interface
+```@docs
+ExponentiatedKernel
+```
+
+## Matern Kernel
+
+The Matern kernel is a Mercer kernel given by:
+
+```math
+\kappa(\mathbf{x},\mathbf{y}) =
+\frac{1}{2^{\nu-1}\Gamma(\nu)}
+\left(\frac{\sqrt{2\nu}||\mathbf{x}-\mathbf{y}||}{\theta}\right)^{\nu}
+K_{\nu}\left(\frac{\sqrt{2\nu}||\mathbf{x}-\mathbf{y}||}{\theta}\right)
+```
+where ``\nu`` and ``\rho`` are positive shape parameters.
+
+### Interface
+```@docs
 MaternKernel
-LinearKernel
+```
+
+## Polynomial Kernel
+The polynomial kernel is a Mercer kernel given by:
+
+```math
+\kappa(\mathbf{x},\mathbf{y}) =
+(a \mathbf{x}^\intercal \mathbf{y} + c)^d
+\qquad \alpha > 0, \; c \geq 0, \; d \in \mathbb{Z}_{+}
+```
+where ``a`` is a positive scale parameter, ``c`` is a non-negative shape parameter and ``d``
+is a shape parameter that determines the degree of the resulting polynomial.
+
+### Interface
+```@docs
 PolynomialKernel
-ExponentiatedKernel
+```
+
+## Periodic Kernel
+The periodic kernel is given by:
+
+```math
+\kappa(\mathbf{x},\mathbf{y}) =
+\exp\left(-\alpha \sum_{i=1}^n \sin(x_i - y_i)^2\right)
+\qquad \alpha > 0
+```
+where ``a`` is a positive scale parameter.
+
+### Interface
+```@docs
 PeriodicKernel
 ```
 
-## Negative Definite Kernels
+## Power Kernel
+The power kernel is given by:
+
+```math
+\kappa(\mathbf{x},\mathbf{y}) =
+\|\mathbf{x} - \mathbf{y} \|^{2\gamma}
+\qquad \gamma \in (0,1]
+```
+where ``\gamma`` is a shape parameter of the Euclidean distance.
+
+### Interface
 ```@docs
 PowerKernel
+```
+
+## Log Kernel
+The log kernel is a negative definite kernel given by:
+
+```math
+\kappa(\mathbf{x},\mathbf{y}) =
+\log \left(1 + \alpha\|\mathbf{x} - \mathbf{y} \|^{2\gamma}\right)
+\qquad \alpha > 0, \; \gamma \in (0,1]
+```
+where ``\alpha`` is a positive scaling parameter and ``\gamma`` is a shape parameter.
+
+### Interface
+```@docs
 LogKernel
 ```
 
-## Other Kernels
+## Sigmoid Kernel
+The Sigmoid Kernel is given by:
+
+```math
+\kappa(\mathbf{x},\mathbf{y}) =
+\tanh(a \mathbf{x}^\intercal \mathbf{y} + c)
+\qquad \alpha > 0, \; c \geq 0
+```
+The sigmoid kernel is a not a true kernel, although it has been used in application.
 
 ```@docs
 SigmoidKernel
diff --git a/src/HyperParameters.jl b/src/HyperParameters.jl
deleted file mode 100644
index 593e228..0000000
--- a/src/HyperParameters.jl
+++ /dev/null
@@ -1,232 +0,0 @@
-module HyperParameters
-
-import Base: convert, eltype, promote_type, show, string, ==, *, /, +, -, ^, isless, depwarn
-
-export
-    Bound,
-        LeftBound,
-        RightBound,
-        NullBound,
-
-    Interval,
-    interval,
-    checkbounds,
-
-    HyperParameter,
-    getvalue,
-    setvalue!,
-    checkvalue,
-    gettheta,
-    settheta!,
-    checktheta,
-    upperboundtheta,
-    lowerboundtheta
-
-
-#= bound.jl =#
-abstract type Bound{T} end
-
-eltype(::Bound{T}) where {T} = T
-
-struct OpenBound{T<:Real} <: Bound{T}
-    value::T
-    function OpenBound{T}(x::Real) where {T<:Real}
-        !(T <: Integer) || error("Bounds must be closed for integers")
-        if T <: AbstractFloat
-            !isnan(x) || error("Bound value must not be NaN")
-            !isinf(x) || error("Bound value must not be Inf/-Inf")
-        end
-        new{T}(x)
-    end
-end
-OpenBound(x::T) where {T<:Real} = OpenBound{T}(x)
-
-convert(::Type{OpenBound{T}}, b::OpenBound{T}) where {T<:Real} = b
-convert(::Type{OpenBound{T}}, b::OpenBound) where {T<:Real} = OpenBound{T}(b.value)
-string(b::OpenBound) = string("OpenBound(", b.value, ")")
-
-struct ClosedBound{T<:Real} <: Bound{T}
-    value::T
-    function ClosedBound{T}(x::Real) where {T<:Real}
-        if T <: AbstractFloat
-            !isnan(x) || error("Bound value must not be NaN")
-            !isinf(x) || error("Bound value must not be Inf/-Inf")
-        end
-        new{T}(x)
-    end
-end
-ClosedBound(x::T) where {T<:Real} = ClosedBound{T}(x)
-
-convert(::Type{ClosedBound{T}}, b::ClosedBound{T}) where {T<:Real} = b
-convert(::Type{ClosedBound{T}}, b::ClosedBound) where {T<:Real} = ClosedBound{T}(b.value)
-string(b::ClosedBound) = string("ClosedBound(", b.value, ")")
-
-
-struct NullBound{T<:Real} <: Bound{T} end
-NullBound(::Type{T}) where {T<:Real} = NullBound{T}()
-
-convert(::Type{NullBound{T}}, b::NullBound{T}) where {T<:Real} = b
-convert(::Type{NullBound{T}}, b::NullBound) where {T<:Real} = NullBound{T}()
-string(b::NullBound{T}) where {T} = string("NullBound(", T, ")")
-
-
-checkvalue(a::NullBound,   x::Real) = true
-checkvalue(a::OpenBound,   x::Real) = a.value <  x
-checkvalue(a::ClosedBound, x::Real) = a.value <= x
-
-checkvalue(x::Real, a::NullBound)   = true
-checkvalue(x::Real, a::OpenBound)   = x <  a.value
-checkvalue(x::Real, a::ClosedBound) = x <= a.value
-
-promote_rule(::Type{Bound{T1}},::Type{Bound{T2}}) where {T1,T2} = Bound{promote_rule(T1,T2)}
-
-function show(io::IO, b::T) where {T<:Bound}
-    print(io, string(b))
-end
-
-#= interval.jl =#
-
-struct Interval{T<:Real,A<:Bound{T},B<:Bound{T}}
-    a::A
-    b::B
-    function Interval{T}(a::A, b::B) where {T<:Real,A<:Bound{T},B<:Bound{T}}
-        if !(A <: NullBound || B <: NullBound)
-            va = a.value
-            vb = b.value
-            if A <: ClosedBound && B <: ClosedBound
-                va <= vb || error("Invalid bounds: a=$va must be less than or equal to b=$vb")
-            else
-                va < vb || error("Invalid bounds: a=$va must be less than b=$vb")
-            end
-        end
-        new{T,A,B}(a,b)
-    end
-end
-Interval(a::Bound{T}, b::Bound{T}) where {T<:Real} = Interval{T}(a,b)
-
-eltype(::Interval{T}) where {T} = T
-
-interval(a::Nothing, b::Nothing) = Interval(NullBound{Float64}(), NullBound{Float64}())
-interval(a::Bound{T}, b::Nothing) where {T<:Real} = Interval(a, NullBound{T}())
-interval(a::Nothing, b::Bound{T}) where {T<:Real} = Interval(NullBound{T}(), b)
-interval(a::Bound{T}, b::Bound{T}) where {T<:Real} = Interval(a,b)
-interval(::Type{T}) where {T<:Real} = Interval(NullBound{T}(), NullBound{T}())
-
-
-checkvalue(I::Interval, x::Real) = checkvalue(I.a, x) && checkvalue(x, I.b)
-
-function theta(I::Interval{T,A,B}, x::T) where {T<:AbstractFloat,A,B}
-    depwarn("theta will be removed entirely in a future release", :theta)
-    checkvalue(I,x) || throw(DomainError(x, "Not in $I"))
-    if A <: OpenBound
-        return B <: OpenBound ? log(x-I.a.value) - log(I.b.value-x) : log(x-I.a.value)
-    else
-        return B <: OpenBound ? log(I.b.value-x) : x
-    end
-end
-
-function upperboundtheta(I::Interval{T,A,B}) where {T<:AbstractFloat,A,B}
-    if B <: ClosedBound
-        return A <: OpenBound ? log(I.b.value - I.a.value) : I.b.value
-    elseif B <: OpenBound
-        return A <: ClosedBound ? log(I.b.value - I.a.value) : convert(T,Inf)
-    else
-        return convert(T,Inf)
-    end
-end
-
-function lowerboundtheta(I::Interval{T,A,B}) where {T<:AbstractFloat,A,B}
-    A <: ClosedBound && !(B <: OpenBound) ? I.a.value : convert(T,-Inf)
-end
-
-function checktheta(I::Interval{T}, x::T) where {T<:AbstractFloat}
-    lowerboundtheta(I) <= x <= upperboundtheta(I)
-end
-
-function eta(I::Interval{T,A,B}, x::T) where {T<:AbstractFloat,A,B}
-    depwarn("eta will be removed entirely in a future release", :eta)
-    checktheta(I,x) || throw(DomainError(x, "Not in $I"))
-    if A <: OpenBound
-        if B <: OpenBound
-            return (I.b.value*exp(x) + I.a.value)/(one(T) + exp(x))
-        else
-            return exp(x) + I.a.value
-        end
-    else
-        return B <: OpenBound ? I.b.value - exp(x) : x
-    end
-end
-
-function string(I::Interval{T1,T2,T3}) where {T1,T2,T3}
-    if T2 <: NullBound
-        if T3 <: NullBound
-            string("interval(", T1, ")")
-        else
-            string("interval(nothing,", string(I.b), ")")
-        end
-    else
-        string("interval(", string(I.a), ",", T3 <: NullBound ? "nothing" : string(I.b), ")")
-    end
-end
-
-function show(io::IO, I::Interval)
-    print(io, string(I))
-end
-
-#= hyperparameter.jl =#
-
-struct HyperParameter{T<:Real}
-    value::Base.RefValue{T}
-    interval::Interval{T}
-    function HyperParameter{T}(x::T, I::Interval{T}) where {T<:Real}
-        checkvalue(I, x) || error("Value $(x) must be in range " * string(I))
-        new{T}(Ref(x), I)
-    end
-end
-HyperParameter(x::T, I::Interval{T} = interval(T)) where {T<:Real} = HyperParameter{T}(x, I)
-
-eltype(::HyperParameter{T}) where {T} = T
-
-@inline getvalue(θ::HyperParameter{T}) where {T} = getindex(θ.value)
-
-function setvalue!(θ::HyperParameter{T}, x::T) where {T}
-    checkvalue(θ.interval, x) || error("Value $(x) must be in range " * string(θ.interval))
-    setindex!(θ.value, x)
-    return θ
-end
-
-checkvalue(θ::HyperParameter{T}, x::T) where {T} = checkvalue(θ.interval, x)
-
-convert(::Type{HyperParameter{T}}, θ::HyperParameter{T}) where {T<:Real} = θ
-function convert(::Type{HyperParameter{T}}, θ::HyperParameter) where {T<:Real}
-    HyperParameter{T}(convert(T, getvalue(θ)), convert(Interval{T}, θ.bounds))
-end
-
-function show(io::IO, θ::HyperParameter{T}) where {T}
-    print(io, string("HyperParameter(", getvalue(θ), ",", string(θ.interval), ")"))
-end
-
-function gettheta(θ::HyperParameter)
-    depwarn("gettheta will be removed entirely in a future release", :gettheta)
-    theta(θ.interval, getvalue(θ))
-end
-
-function settheta!(θ::HyperParameter, x::T) where {T}
-    depwarn("settheta! will be removed entirely in a future release", :(settheta!))
-    setvalue!(θ, eta(θ.interval,x))
-end
-
-function checktheta(θ::HyperParameter, x::T) where {T}
-    depwarn("checktheta will be removed entirely in a future release", :checktheta)
-    checktheta(θ.interval, x)
-end
-
-for op in (:isless, :(==), :+, :-, :*, :/)
-    @eval begin
-        $op(θ1::HyperParameter, θ2::HyperParameter) = $op(getvalue(θ1), getvalue(θ2))
-        $op(a::Number, θ::HyperParameter) = $op(a, getvalue(θ))
-        $op(θ::HyperParameter, a::Number) = $op(getvalue(θ), a)
-    end
-end
-
-end # End HyperParameter
\ No newline at end of file
diff --git a/src/MLKernels.jl b/src/MLKernels.jl
index 9bf2a23..854203d 100644
--- a/src/MLKernels.jl
+++ b/src/MLKernels.jl
@@ -8,46 +8,31 @@ import Base: convert, eltype, print, show, string, ==, *, /, +, -, ^, exp, tanh
 
 export
 
-    # Hyper Parameters
-    Bound,
-        OpenBound,
-        ClosedBound,
-        NullBound,
-
-    Interval,
-    interval,
-
-    HyperParameter,
-    getvalue,
-    setvalue!,
-    checkvalue,
-    gettheta,
-    settheta!,
-    checktheta,
-
     # Memory
     Orientation,
 
     # Kernel Functions
     Kernel,
-    MercerKernel,
-    ExponentialKernel,
-    LaplacianKernel,
-    SquaredExponentialKernel,
-    GaussianKernel,
-    RadialBasisKernel,
-    GammaExponentialKernel,
-    RationalQuadraticKernel,
-    GammaRationalKernel,
-    MaternKernel,
-    LinearKernel,
-    PolynomialKernel,
-    ExponentiatedKernel,
-    PeriodicKernel,
-    NegativeDefiniteKernel,
-    PowerKernel,
-    LogKernel,
-    SigmoidKernel,
+        MercerKernel,
+            AbstractExponentialKernel,
+                ExponentialKernel,
+                LaplacianKernel,
+                SquaredExponentialKernel,
+                GaussianKernel,
+                RadialBasisKernel,
+                GammaExponentialKernel,
+            AbstractRationalQuadraticKernel,
+                RationalQuadraticKernel,
+                GammaRationalQuadraticKernel,
+            MaternKernel,
+            LinearKernel,
+            PolynomialKernel,
+            ExponentiatedKernel,
+            PeriodicKernel,
+        NegativeDefiniteKernel,
+            PowerKernel,
+            LogKernel,
+        SigmoidKernel,
 
     # Kernel Function Properties
     ismercer,
@@ -71,28 +56,6 @@ using SpecialFunctions: besselk, gamma
 import LinearAlgebra
 import Statistics
 
-include("HyperParameters.jl")
-using MLKernels.HyperParameters:
-    Bound,
-        OpenBound,
-        ClosedBound,
-        NullBound,
-
-    Interval,
-    interval,
-
-    HyperParameter,
-    getvalue,
-    setvalue!,
-    checkvalue,
-    gettheta,
-    checktheta,
-    settheta!,
-    lowerboundtheta,
-    upperboundtheta
-
-import MLKernels.HyperParameters: gettheta, checktheta, settheta!
-
 @doc raw"""
     Orientation
 
diff --git a/src/basefunctions.jl b/src/basefunctions.jl
index 0c7fb92..5d5032d 100644
--- a/src/basefunctions.jl
+++ b/src/basefunctions.jl
@@ -23,15 +23,6 @@ end
 
 abstract type PreMetric <: BaseFunction end
 
-const pre_metrics = [
-    "chisquared",
-    "sinesquared"
-]
-
-for fname in pre_metrics
-    include(joinpath("basefunctions", "$(fname).jl"))
-end
-
 
 # Metrics ==================================================================================
 
@@ -43,4 +34,4 @@ const metrics = [
 
 for fname in metrics
     include(joinpath("basefunctions", "$(fname).jl"))
-end
+end
\ No newline at end of file
diff --git a/src/basefunctions/chisquared.jl b/src/basefunctions/chisquared.jl
deleted file mode 100644
index dfdbfa1..0000000
--- a/src/basefunctions/chisquared.jl
+++ /dev/null
@@ -1,14 +0,0 @@
-@doc raw"""
-    ChiSquared()
-
-The Chi-Squared base function is given by:
-
-```math
-f(\mathbf{x}, \mathbf{y}) = \sum_i \frac{(x_i - y_i)^2}{x_i + y_i}
-```
-"""
-struct ChiSquared <: PreMetric end
-
-@inline function base_aggregate(::ChiSquared, s::T, x::T, y::T) where {T}
-    x == y == zero(T) ? s : s + (x-y)^2/(x+y)
-end
diff --git a/src/basefunctions/sinesquared.jl b/src/basefunctions/sinesquared.jl
deleted file mode 100644
index e976c17..0000000
--- a/src/basefunctions/sinesquared.jl
+++ /dev/null
@@ -1,12 +0,0 @@
-@doc raw"""
-    SineSquared
-
-The Sine-Squared base function is given by:
-
-```math
-f(\mathbf{x}, \mathbf{y}) = \sine^2 \left(x_i - y_i\right)
-``` 
-"""
-struct SineSquared <: PreMetric end
-@inline base_aggregate(::SineSquared, s::T, x::T, y::T) where {T} = s + sin(x-y)^2
-@inline isstationary(::SineSquared) = true
\ No newline at end of file
diff --git a/src/basefunctions/squaredeuclidean.jl b/src/basefunctions/squaredeuclidean.jl
index 9ef1eb9..dc1cea4 100644
--- a/src/basefunctions/squaredeuclidean.jl
+++ b/src/basefunctions/squaredeuclidean.jl
@@ -7,7 +7,7 @@ The squared Euclidean function is defined by:
 f(\mathbf{x}, \mathbf{y}) = (\mathbf{x} - \mathbf{y})^{\intercal}(\mathbf{x} - \mathbf{y})
 ```
 """
-struct SquaredEuclidean <: PreMetric end
+struct SquaredEuclidean <: Metric end
 
 @inline base_aggregate(::SquaredEuclidean, s::T, x::T, y::T) where {T} = s + (x-y)^2
 
diff --git a/src/basematrix.jl b/src/basematrix.jl
index 42fbaa3..adc852a 100644
--- a/src/basematrix.jl
+++ b/src/basematrix.jl
@@ -145,41 +145,6 @@ function basematrix!(
 end
 
 
-function basematrix(
-        σ::Orientation,
-        f::BaseFunction,
-        X::AbstractMatrix{T},
-        symmetrize::Bool = true
-    ) where {T<:AbstractFloat}
-    basematrix!(σ, allocate_basematrix(σ, X), f, X, symmetrize)
-end
-
-function basematrix(
-        f::BaseFunction,
-        X::AbstractMatrix,
-        symmetrize::Bool = true
-    )
-    basematrix(Val(:row), f, X, symmetrize)
-end
-
-function basematrix(
-        σ::Orientation,
-        f::BaseFunction,
-        X::AbstractMatrix{T},
-        Y::AbstractMatrix{T}
-    ) where {T<:AbstractFloat}
-    basematrix!(σ, allocate_basematrix(σ, X, Y), f, X, Y)
-end
-
-function basematrix(
-        f::BaseFunction,
-        X::AbstractMatrix,
-        Y::AbstractMatrix
-    )
-    basematrix(Val(:row), f, X, Y)
-end
-
-
 # ScalarProduct using BLAS/Built-In methods ================================================
 
 @inline function basematrix!(
@@ -314,4 +279,4 @@ function basematrix!(
     yᵀy = dotvectors(σ, Y)
     squared_distance!(P, xᵀx, yᵀy)
     fix_negatives!(σ, P, X, Y)
-end
+end
\ No newline at end of file
diff --git a/src/deprecated.jl b/src/deprecated.jl
index b80ab35..2ad2962 100644
--- a/src/deprecated.jl
+++ b/src/deprecated.jl
@@ -1,67 +1,26 @@
-import Base: @deprecate, depwarn
-
-abstract type MemoryLayout end
-
-struct ColumnMajor <: MemoryLayout
-    function ColumnMajor()
-        depwarn("Use `Val(:col)` instead of `ColumnMajor()`", :ColumnMajor)
-        new()
+# Removal of LinearKernel
+Base.@deprecate LinearKernel(a::Real=1, c::Real=1) PolynomialKernel(a, c, 1)
+
+# Renaming of GammaRational to GammaRationalQuadratic
+Base.@deprecate GammaRationalKernel GammaRationalQuadraticKernel
+
+# Removal of PeriodicKernel
+struct SineSquared <: PreMetric end
+@inline base_aggregate(::SineSquared, s::T, x::T, y::T) where {T} = s + sin(x-y)^2
+@inline isstationary(::SineSquared) = true
+
+struct PeriodicKernel{T<:AbstractFloat} <: MercerKernel{T}
+    α::T
+    function PeriodicKernel{T}(α::Real) where {T<:AbstractFloat}
+        Base.depwarn("PeriodicKernel will be removed in the next major release", :PeriodicKernel)
+        @check_args(PeriodicKernel, α, α > zero(α), "α > 0")
+        new{T}(α)
     end
 end
+PeriodicKernel(α::T₁ = 1.0) where {T₁<:Real} = PeriodicKernel{promote_float(T₁)}(α)
 
-struct RowMajor <: MemoryLayout
-    function RowMajor()
-        depwarn("Use `Val(:row)` instead of `RowMajor()`", :RowMajor)
-        new()
-    end
-end
-
-@deprecate RowMajor Val{:row}
-@deprecate ColumnMajor Val{:col}
-
-layout_map(orient) = typeof(orient) <: RowMajor ? Val(:row) : Val(:col)
-
-col_warn = "Use `Val(:col)` instead of `ColumnMajor()`"
-row_warn = "Use `Val(:row)` instead of `RowMajor()`"
-
-function kernelmatrix!(
-        σ::MemoryLayout,
-        P::Matrix,
-        κ::Kernel,
-        X::AbstractMatrix,
-        symmetrize::Bool = true
-    )
-    orientation = layout_map(σ)
-    kernelmatrix!(orientation, P, κ, X, symmetrize)
-end
-
-function kernelmatrix!(
-        σ::MemoryLayout,
-        P::Matrix,
-        κ::Kernel,
-        X::AbstractMatrix,
-        Y::AbstractMatrix
-    )
-    orientation = layout_map(σ)
-    kernelmatrix!(orientation, P, κ, X, Y)
-end
-
-function kernelmatrix(
-        σ::MemoryLayout,
-        κ::Kernel,
-        X::AbstractMatrix,
-        symmetrize::Bool = true
-    )
-    orientation = layout_map(σ)
-    kernelmatrix(orientation, κ, X, symmetrize)
-end
+@inline basefunction(::PeriodicKernel) = SineSquared()
 
-function kernelmatrix(
-        σ::MemoryLayout,
-        κ::Kernel,
-        X::AbstractMatrix,
-        Y::AbstractMatrix
-    )
-    orientation = layout_map(σ)
-    kernelmatrix(orientation, κ, X, Y)
+@inline function kappa(κ::PeriodicKernel{T}, z::T) where {T}
+    return exp(-κ.α*z)
 end
\ No newline at end of file
diff --git a/src/kernelfunctions.jl b/src/kernelfunctions.jl
index 4dfdc40..9cb29fe 100644
--- a/src/kernelfunctions.jl
+++ b/src/kernelfunctions.jl
@@ -2,20 +2,16 @@
 
 abstract type Kernel{T<:AbstractFloat} end
 
-function string(κ::Kernel)
-    args = [string(getvalue(getfield(κ,θ))) for θ in fieldnames(typeof(κ))]
+function string(κ::Kernel{T}) where {T}
+    args = [string(getfield(κ,θ)) for θ in fieldnames(typeof(κ))]
     kernelname = typeof(κ).name.name
-    string(kernelname, "(", join(args, ","), ")")
+    string(kernelname, "{", string(T), "}(", join(args, ","), ")")
 end
 
 function show(io::IO, κ::Kernel)
     print(io, string(κ))
 end
 
-function basefunction(::Kernel)
-    error("No base function specified for kernel")
-end
-
 @inline eltype(::Type{<:Kernel{E}}) where {E} = E
 @inline eltype(κ::Kernel) = eltype(typeof(κ))
 
@@ -47,57 +43,20 @@ Returns `true` if the kernel `κ` is an isotropic kernel; `false` otherwise.
 """
 isisotropic(κ::Kernel)  = isisotropic(basefunction(κ))
 
-thetafieldnames(κ::Kernel) = fieldnames(typeof(κ))
-
-gettheta(κ::Kernel{T}) where {T} = T[gettheta(getfield(κ,θ)) for θ in thetafieldnames(κ)]
-
-function settheta!(κ::Kernel{T},v::Vector{T}) where {T}
-    fields = thetafieldnames(κ)
-    if length(fields) != length(v)
-        throw(DimensionMismatch("Update vector has invalid length"))
-    end
-    for i in eachindex(fields)
-        settheta!(getfield(κ, fields[i]), v[i])
-    end
-    return κ
-end
-
-function checktheta(κ::Kernel{T},v::Vector{T}) where {T}
-    fields = thetafieldnames(κ)
-    if length(fields) != length(v)
-        throw(DimensionMismatch("Update vector has invalid length"))
-    end
-    for i in eachindex(fields)
-        if !checktheta(getfield(κ, fields[i]), v[i])
-            return false
-        end
-    end
-    return true
-end
-
-function floattype(T_i::DataType...)
-    T_max = promote_type(T_i...)
-    T_max <: AbstractFloat ? T_max : Float64
-end
-
 
 # Mercer Kernels ===========================================================================
 
 abstract type MercerKernel{T<:AbstractFloat} <: Kernel{T} end
+
 @inline ismercer(::MercerKernel) = true
 
 const mercer_kernels = [
     "exponential",
-    "squaredexponential",
-    "gammaexponential",
+    "exponentiated",
     "rationalquadratic",
-    "gammarational",
     "matern",
-    "linear",
-    "polynomial",
-    "exponentiated",
-    "periodic"
-    ]
+    "polynomial"
+]
 
 for kname in mercer_kernels
     include(joinpath("kernelfunctions", "mercer", "$(kname).jl"))
@@ -107,11 +66,12 @@ end
 # Negative Definite Kernels ================================================================
 
 abstract type NegativeDefiniteKernel{T<:AbstractFloat} <: Kernel{T} end
+
 @inline isnegdef(::NegativeDefiniteKernel) = true
 
 const negdef_kernels = [
-    "power",
-    "log"
+    "log",
+    "power"
 ]
 
 for kname in negdef_kernels
@@ -127,45 +87,4 @@ const other_kernels = [
 
 for kname in other_kernels
     include(joinpath("kernelfunctions", "$(kname).jl"))
-end
-
-for κ in [
-        ExponentialKernel,
-        SquaredExponentialKernel,
-        GammaExponentialKernel,
-        RationalQuadraticKernel,
-        GammaRationalKernel,
-        MaternKernel,
-        LinearKernel,
-        PolynomialKernel,
-        ExponentiatedKernel,
-        PeriodicKernel,
-        PowerKernel,
-        LogKernel,
-        SigmoidKernel
-    ]
-    κ_sym = nameof(κ)
-    κ_args = [:(getvalue(κ.$(θ))) for θ in fieldnames(κ)]
-
-    @eval begin
-        function ==(κ1::$(κ_sym), κ2::$(κ_sym))
-            mapreduce(θ -> getfield(κ1,θ) == getfield(κ2,θ), &, fieldnames(typeof(κ1)), init = true)
-        end
-    end
-
-    @eval begin
-        function convert(::Type{$(κ_sym){T}}, κ::$(κ_sym)) where {T}
-            $(Expr(:call, :($(κ_sym){T}), κ_args...))
-        end
-    end
-
-    κs = supertype(κ)
-    while κs != Any
-        @eval begin
-            function convert(::Type{$(nameof(κs)){T}}, κ::$(κ_sym)) where {T}
-                convert($(κ_sym){T}, κ)
-            end
-        end
-        κs = supertype(κs)
-    end
-end
+end
\ No newline at end of file
diff --git a/src/kernelfunctions/mercer/exponential.jl b/src/kernelfunctions/mercer/exponential.jl
index 37e5d73..0e3f6f9 100644
--- a/src/kernelfunctions/mercer/exponential.jl
+++ b/src/kernelfunctions/mercer/exponential.jl
@@ -1,29 +1,158 @@
+# Abstract Exponential Kernel ==============================================================
+
+abstract type AbstractExponentialKernel{T<:AbstractFloat} <: MercerKernel{T} end
+
+@inline basefunction(::AbstractExponentialKernel) = SquaredEuclidean()
+
+
+# Exponential Kernel =======================================================================
 @doc raw"""
     ExponentialKernel([α=1])
 
-The exponential kernel is given by the formula:
+The exponential kernel is an isotropic Mercer kernel given by the formula:
+
+```
+    κ(x,y) = exp(α‖x-y‖)   α > 0
+```
+
+where `α` is a positive scaling parameter. See also [`SquaredExponentialKernel`](@ref) for
+a related form of the kernel or [`GammaExponentialKernel`](@ref) for a generalization.
+
+# Examples
+
+```jldoctest; setup = :(using MLKernels)
+julia> ExponentialKernel()
+ExponentialKernel{Float64}(1.0)
+
+julia> ExponentialKernel(2.0f0)
+ExponentialKernel{Float32}(2.0)
+```
+"""
+struct ExponentialKernel{T<:AbstractFloat} <: AbstractExponentialKernel{T}
+    α::T
+    function ExponentialKernel{T}(α::Real=T(1)) where {T<:AbstractFloat}
+        @check_args(ExponentialKernel, α, α > zero(T), "α > 0")
+        return new{T}(α)
+    end
+end
+ExponentialKernel(α::T=1.0) where {T<:Real} = ExponentialKernel{promote_float(T)}(α)
+
+@inline kappa(κ::ExponentialKernel{T}, d²::T) where {T} = exp(-κ.α*√(d²))
+
+function convert(::Type{K}, κ::ExponentialKernel) where {K>:ExponentialKernel{T}} where T
+    return ExponentialKernel{T}(κ.α)
+end
+
+"""
+    LaplacianKernel([α=1])
+
+Alias for [`ExponentialKernel`](@ref).
+"""
+const LaplacianKernel = ExponentialKernel
+
+
+# Squared Exponential Kernel ===============================================================
+@doc raw"""
+    SquaredExponentialKernel([α=1])
+
+The squared exponential kernel is an isotropic Mercer kernel given by the formula:
+
+```
+    κ(x,y) = exp(α‖x-y‖²)   α > 0
+```
+
+where `α` is a positive scaling parameter. See also [`ExponentialKernel`](@ref) for a
+related form of the kernel or [`GammaExponentialKernel`](@ref) for a generalization.
 
-```math
-\kappa(\mathbf{x},\mathbf{y}) = \exp\left(-\alpha ||\mathbf{x} - \mathbf{y}||\right) 
-\qquad \alpha > 0
+# Examples
+
+```jldoctest; setup = :(using MLKernels)
+julia> SquaredExponentialKernel()
+SquaredExponentialKernel{Float64}(1.0)
+
+julia> SquaredExponentialKernel(2.0f0)
+SquaredExponentialKernel{Float32}(2.0)
 ```
-  
-where ``\alpha`` is a scaling parameter of the Euclidean distance. The exponential kernel, 
-also known as the Laplacian kernel, is an isotropic Mercer kernel. The constructor is 
-aliased by `LaplacianKernel`, so both names may be used:
 """
-struct ExponentialKernel{T<:AbstractFloat} <: MercerKernel{T}
-    alpha::HyperParameter{T}
-    ExponentialKernel{T}(α::Real) where {T<:AbstractFloat} = new{T}(
-        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing))
-    )
+struct SquaredExponentialKernel{T<:AbstractFloat} <: AbstractExponentialKernel{T}
+    α::T
+    function SquaredExponentialKernel{T}(α::Real=T(1)) where {T<:AbstractFloat}
+        @check_args(SquaredExponentialKernel, α, α > zero(T), "α > 0")
+        return new{T}(α)
+    end
+end
+function SquaredExponentialKernel(α::T=1.0) where {T<:Real}
+    return SquaredExponentialKernel{promote_float(T)}(α)
+end
+
+@inline kappa(κ::SquaredExponentialKernel{T}, d²::T) where {T} = exp(-κ.α*d²)
+
+function convert(
+        ::Type{K},
+        κ::SquaredExponentialKernel
+    ) where {K>:SquaredExponentialKernel{T}} where T
+    return SquaredExponentialKernel{T}(κ.α)
+end
+
+"""
+    GaussianKernel([α=1])
+
+Alias of [`SquaredExponentialKernel`](@ref).
+"""
+const GaussianKernel = SquaredExponentialKernel
+
+"""
+    RadialBasisKernel([α=1])
+
+Alias of [`SquaredExponentialKernel`](@ref).
+"""
+const RadialBasisKernel = SquaredExponentialKernel
+
+
+# Gamma Exponential Kernel =================================================================
+@doc raw"""
+    GammaExponentialKernel([α=1 [,γ=1]])
+
+The ``\gamma``-exponential kernel is an isotropic Mercer kernel given by the formula:
+
+```
+    κ(x,y) = exp(α‖x-y‖²ᵞ)   α > 0, γ ∈ (0,1]
+```
+where `α` is a scaling parameter and `γ` is a shape parameter of the Euclidean distance.
+When `γ = 1` use [`SquaredExponentialKernel`](@ref) and [`SquaredExponentialKernel`](@ref)
+when `γ = 0.5` since these are more efficient implementations.
+
+# Examples
+
+```jldoctest; setup = :(using MLKernels)
+julia> GammaExponentialKernel()
+GammaExponentialKernel{Float64}(1.0,1.0)
+
+julia> GammaExponentialKernel(2.0f0)
+GammaExponentialKernel{Float32}(2.0,1.0)
+
+julia> GammaExponentialKernel(2.0, 0.5)
+GammaExponentialKernel{Float64}(2.0,0.5)
+```
+"""
+struct GammaExponentialKernel{T<:AbstractFloat} <: AbstractExponentialKernel{T}
+    α::T
+    γ::T
+    function GammaExponentialKernel{T}(α::Real=T(1), γ::Real=T(1)) where {T<:AbstractFloat}
+        @check_args(GammaExponentialKernel, α, α > zero(T), "α > 0")
+        @check_args(GammaExponentialKernel, γ, one(T) >= γ > zero(T), "γ ∈ (0,1]")
+        return new{T}(α, γ)
+    end
+end
+function GammaExponentialKernel(α::T₁=1.0, γ::T₂=T₁(1)) where {T₁<:Real, T₂<:Real}
+    return GammaExponentialKernel{promote_float(T₁,T₂)}(α, γ)
 end
-ExponentialKernel(α::T=1.0) where {T<:Real} = ExponentialKernel{floattype(T)}(α)
-LaplacianKernel = ExponentialKernel
 
-@inline exponentialkernel(z::T, α::T) where {T<:AbstractFloat} = exp(-α*sqrt(z))
+@inline kappa(κ::GammaExponentialKernel{T}, d²::T) where {T} = exp(-κ.α*d²^κ.γ)
 
-@inline basefunction(::ExponentialKernel) = SquaredEuclidean()
-@inline function kappa(κ::ExponentialKernel{T}, z::T) where {T<:AbstractFloat}
-    exponentialkernel(z, getvalue(κ.alpha))
+function convert(
+        ::Type{K},
+        κ::GammaExponentialKernel
+    ) where {K>:GammaExponentialKernel{T}} where T
+    return GammaExponentialKernel{T}(κ.α, κ.γ)
 end
\ No newline at end of file
diff --git a/src/kernelfunctions/mercer/exponentiated.jl b/src/kernelfunctions/mercer/exponentiated.jl
index e1aac9d..5b45497 100644
--- a/src/kernelfunctions/mercer/exponentiated.jl
+++ b/src/kernelfunctions/mercer/exponentiated.jl
@@ -1,22 +1,42 @@
 @doc raw"""
-    ExponentiatedKernel([a=1])
+    ExponentiatedKernel([α=1])
 
 The exponentiated kernel is a Mercer kernel given by:
 
-```math
-\kappa(\mathbf{x},\mathbf{y}) = \exp\left(a \mathbf{x}^\intercal \mathbf{y} \right) 
-\qquad a > 0
+```
+    κ(x,y) = exp(α⋅xᵀy)   α > 0
+```
+where `α` is a positive scaling parameter.
+
+# Examples
+
+```jldoctest; setup = :(using MLKernels)
+julia> ExponentiatedKernel()
+ExponentiatedKernel{Float64}(1.0)
+
+julia> ExponentiatedKernel(2)
+ExponentiatedKernel{Float64}(2.0)
+
+julia> ExponentiatedKernel(2.0f0)
+ExponentiatedKernel{Float32}(2.0)
 ```
 """
 struct ExponentiatedKernel{T<:AbstractFloat} <: MercerKernel{T}
-    alpha::HyperParameter{T}
-    ExponentiatedKernel{T}(α::Real) where {T<:AbstractFloat} = new{T}(
-        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing))
-    )
+    α::T
+    function ExponentiatedKernel{T}(α::Real=T(1)) where {T<:AbstractFloat}
+        @check_args(ExponentiatedKernel, α, α > zero(T), "α > 0")
+        return new{T}(α)
+    end
 end
-ExponentiatedKernel(α::T1 = 1.0) where {T1<:Real} = ExponentiatedKernel{floattype(T1)}(α)
-
-@inline exponentiatedkernel(z::T, α::T) where {T<:AbstractFloat} = exp(α*z)
+ExponentiatedKernel(α::T=1.0) where {T<:Real} = ExponentiatedKernel{promote_float(T)}(α)
 
 @inline basefunction(::ExponentiatedKernel) = ScalarProduct()
-@inline kappa(κ::ExponentiatedKernel{T}, z::T) where {T} = exponentiatedkernel(z, getvalue(κ.alpha))
\ No newline at end of file
+
+@inline kappa(κ::ExponentiatedKernel{T}, xᵀy::T) where {T} = exp(κ.α*xᵀy)
+
+function convert(
+        ::Type{K},
+        κ::ExponentiatedKernel
+    ) where {K>:ExponentiatedKernel{T}} where T
+    return ExponentiatedKernel{T}(κ.α)
+end
\ No newline at end of file
diff --git a/src/kernelfunctions/mercer/gammaexponential.jl b/src/kernelfunctions/mercer/gammaexponential.jl
deleted file mode 100644
index 6526220..0000000
--- a/src/kernelfunctions/mercer/gammaexponential.jl
+++ /dev/null
@@ -1,30 +0,0 @@
-@doc raw"""
-    GammaExponentialKernel([α=1 [,γ=1]])
-
-The gamma exponential kernel is a generalization of the exponential and squared exponential 
-kernels:
-
-```math
-\kappa(\mathbf{x},\mathbf{y}) = \exp\left(-\alpha ||\mathbf{x} - \mathbf{y}||^{\gamma} 
-\right) \qquad \alpha > 0, \; 0 < \gamma \leq 1
-```
-where ``\alpha`` is a scaling parameter and ``\gamma`` is a shape parameter.
-"""
-struct GammaExponentialKernel{T<:AbstractFloat} <: MercerKernel{T}
-    alpha::HyperParameter{T}
-    gamma::HyperParameter{T}
-    GammaExponentialKernel{T}(α::Real, γ::Real) where {T<:AbstractFloat} = new{T}(
-        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing)),
-        HyperParameter(convert(T,γ), interval(OpenBound(zero(T)), ClosedBound(one(T))))
-    )
-end
-function GammaExponentialKernel(α::T1=1.0, γ::T2=one(T1)) where {T1<:Real,T2<:Real}
-    GammaExponentialKernel{floattype(T1,T2)}(α,γ)
-end
-
-@inline gammaexponentialkernel(z::T, α::T, γ::T) where {T<:AbstractFloat} = exp(-α*z^γ)
-
-@inline basefunction(::GammaExponentialKernel) = SquaredEuclidean()
-@inline function kappa(κ::GammaExponentialKernel{T}, z::T) where {T}
-    gammaexponentialkernel(z, getvalue(κ.alpha), getvalue(κ.gamma))
-end
\ No newline at end of file
diff --git a/src/kernelfunctions/mercer/gammarational.jl b/src/kernelfunctions/mercer/gammarational.jl
deleted file mode 100644
index 1615076..0000000
--- a/src/kernelfunctions/mercer/gammarational.jl
+++ /dev/null
@@ -1,38 +0,0 @@
-@doc raw"""
-    GammaRationalKernel([α [,β [,γ]]])
-  
-The gamma-rational kernel is a generalization of the rational-quadratic kernel with an 
-additional shape parameter:
-
-```math
-\kappa(\mathbf{x},\mathbf{y})
-= \left(1 +\alpha ||\mathbf{x},\mathbf{y}||^{\gamma}\right)^{-\beta} 
-\qquad \alpha > 0, \; \beta > 0, \; 0 < \gamma \leq 1
-```
-
-where ``\alpha`` is a scaling parameter and ``\beta`` and ``\gamma`` are shape parameters.
-"""
-struct GammaRationalKernel{T<:AbstractFloat} <: MercerKernel{T}
-    alpha::HyperParameter{T}
-    beta::HyperParameter{T}
-    gamma::HyperParameter{T}
-    GammaRationalKernel{T}(α::Real, β::Real, γ::Real) where {T<:AbstractFloat} = new{T}(
-        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing)),
-        HyperParameter(convert(T,β), interval(OpenBound(zero(T)), nothing)),
-        HyperParameter(convert(T,γ), interval(OpenBound(zero(T)), ClosedBound(one(T))))
-    )
-end
-function GammaRationalKernel(
-        α::T1 = 1.0,
-        β::T2 = one(T1),
-        γ::T3 = one(floattype(T1,T2))
-    ) where {T1<:Real,T2<:Real,T3<:Real}
-    GammaRationalKernel{floattype(T1,T2,T3)}(α,β,γ)
-end
-
-@inline gammarationalkernel(z::T, α::T, β::T, γ::T) where {T<:AbstractFloat} = (1 + α*(z^γ))^(-β)
-
-@inline basefunction(::GammaRationalKernel) = SquaredEuclidean()
-@inline function kappa(κ::GammaRationalKernel{T}, z::T) where {T}
-    gammarationalkernel(z, getvalue(κ.alpha), getvalue(κ.beta), getvalue(κ.gamma))
-end
\ No newline at end of file
diff --git a/src/kernelfunctions/mercer/linear.jl b/src/kernelfunctions/mercer/linear.jl
deleted file mode 100644
index 716a61f..0000000
--- a/src/kernelfunctions/mercer/linear.jl
+++ /dev/null
@@ -1,24 +0,0 @@
-@doc raw"""
-    LinearKernel([a=1 [,c=1]])
-
-The linear kernel is a Mercer kernel given by:
-
-```math
-\kappa(\mathbf{x},\mathbf{y}) = 
-a \mathbf{x}^\intercal \mathbf{y} + c \qquad \alpha > 0, \; c \geq 0
-```
-"""
-struct LinearKernel{T<:AbstractFloat} <: MercerKernel{T}
-    a::HyperParameter{T}
-    c::HyperParameter{T}
-    LinearKernel{T}(a::Real, c::Real) where {T<:AbstractFloat} = new{T}(
-        HyperParameter(convert(T,a), interval(OpenBound(zero(T)), nothing)),
-        HyperParameter(convert(T,c), interval(ClosedBound(zero(T)), nothing))
-    )
-end
-LinearKernel(a::T1=1.0, c::T2=one(T1)) where {T1<:Real,T2<:Real} = LinearKernel{floattype(T1,T2)}(a,c)
-
-@inline linearkernel(z::T, a::T, c::T) where {T<:AbstractFloat} = a*z + c
-
-@inline basefunction(::LinearKernel) = ScalarProduct()
-@inline kappa(κ::LinearKernel{T}, z::T) where {T} = linearkernel(z, getvalue(κ.a), getvalue(κ.c))
diff --git a/src/kernelfunctions/mercer/matern.jl b/src/kernelfunctions/mercer/matern.jl
index 199b7e1..a81694b 100644
--- a/src/kernelfunctions/mercer/matern.jl
+++ b/src/kernelfunctions/mercer/matern.jl
@@ -1,38 +1,44 @@
 @doc raw"""
-    MaternKernel([ν=1 [,θ=1]])
+    MaternKernel([ν=1 [, θ=1]])
 
-The Matern kernel is a Mercer kernel given by:
+The Matern kernel is a Mercer kernel with parameters `ν > 0` and `ρ > 0`. See the published
+documentation for the full definition of the function.
 
-```math
-\kappa(\mathbf{x},\mathbf{y}) =
-\frac{1}{2^{\nu-1}\Gamma(\nu)}
-\left(\frac{\sqrt{2\nu}||\mathbf{x}-\mathbf{y}||}{\theta}\right)^{\nu}
-K_{\nu}\left(\frac{\sqrt{2\nu}||\mathbf{x}-\mathbf{y}||}{\theta}\right)
-```
+# Examples
+
+```jldoctest; setup = :(using MLKernels)
+julia> MaternKernel()
+MaternKernel{Float64}(1.0,1.0)
 
-where ``\Gamma`` is the gamma function, ``K_{\nu}`` is the modified Bessel function of the
-second kind, ``\nu > 0`` and ``\theta > 0``.
+julia> MaternKernel(2.0f0)
+MaternKernel{Float32}(2.0,1.0)
+
+julia> MaternKernel(2.0f0, 2.0)
+MaternKernel{Float64}(2.0,2.0)
+```
 """
 struct MaternKernel{T<:AbstractFloat} <: MercerKernel{T}
-    nu::HyperParameter{T}
-    rho::HyperParameter{T}
-    MaternKernel{T}(ν::Real, ρ::Real) where {T<:AbstractFloat}  = new{T}(
-        HyperParameter(convert(T,ν), interval(OpenBound(zero(T)), nothing)),
-        HyperParameter(convert(T,ρ), interval(OpenBound(zero(T)), nothing))
-    )
+    ν::T
+    ρ::T
+    function MaternKernel{T}(ν::Real=T(1), ρ::Real=T(1)) where {T<:AbstractFloat}
+        @check_args(MaternKernel, ν, ν > zero(T), "ν > 0")
+        @check_args(MaternKernel, ρ, ρ > zero(T), "ρ > 0")
+        return new{T}(ν, ρ)
+    end
 end
-function MaternKernel(ν::T1=1.0, ρ::T2=one(T1)) where {T1<:Real,T2<:Real}
-    MaternKernel{floattype(T1,T2)}(ν,ρ)
+function MaternKernel(ν::T₁=1.0, ρ::T₂=T₁(1)) where {T₁<:Real,T₂<:Real}
+    MaternKernel{promote_float(T₁,T₂)}(ν,ρ)
 end
 
-@inline function maternkernel(d²::T, ν::T, ρ::T) where {T<:AbstractFloat}
+@inline basefunction(::MaternKernel) = SquaredEuclidean()
+
+@inline function kappa(κ::MaternKernel{T}, d²::T) where {T}
     d = √(d²)
     d = d < eps(T) ? eps(T) : d  # If d is zero, besselk will return NaN
-    tmp = √(2ν)d/ρ
-    (2^(1 - ν))*(tmp^ν)*besselk(ν, tmp)/gamma(ν)
+    tmp = √(2κ.ν)*d/κ.ρ
+    return (convert(T, 2)^(one(T) - κ.ν))*(tmp^κ.ν)*besselk(κ.ν, tmp)/gamma(κ.ν)
 end
 
-@inline basefunction(::MaternKernel) = SquaredEuclidean()
-@inline function kappa(κ::MaternKernel{T}, z::T) where {T}
-    maternkernel(z, getvalue(κ.nu), getvalue(κ.rho))
+function convert(::Type{K}, κ::MaternKernel) where {K>:MaternKernel{T}} where T
+    return MaternKernel{T}(κ.ν, κ.ρ)
 end
\ No newline at end of file
diff --git a/src/kernelfunctions/mercer/periodic.jl b/src/kernelfunctions/mercer/periodic.jl
deleted file mode 100644
index 6aabb1b..0000000
--- a/src/kernelfunctions/mercer/periodic.jl
+++ /dev/null
@@ -1,25 +0,0 @@
-@doc raw"""
-    PeriodicKernel([α=1 [,p=π]])
-
-The periodic kernel is given by:
-
-```math
-\kappa(\mathbf{x},\mathbf{y}) =
-\exp\left(-\alpha \sum_{i=1}^n \sin(p(x_i - y_i))^2\right)
-\qquad p >0, \; \alpha > 0
-```
-
-where ``\mathbf{x}`` and ``\mathbf{y}`` are ``n`` dimensional vectors. The parameters ``p`` 
-and ``\alpha`` are scaling parameters for the periodicity and the magnitude, respectively. 
-This kernel is useful when data has periodicity to it.
-"""
-struct PeriodicKernel{T<:AbstractFloat} <: MercerKernel{T}
-    alpha::HyperParameter{T}
-    PeriodicKernel{T}(α::Real) where {T<:AbstractFloat} = new{T}(
-        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing))
-    )
-end
-PeriodicKernel(α::T1 = 1.0) where {T1<:Real} = PeriodicKernel{floattype(T1)}(α)
-
-@inline basefunction(::PeriodicKernel) = SineSquared()
-@inline kappa(κ::PeriodicKernel{T}, z::T) where {T} = squaredexponentialkernel(z, getvalue(κ.alpha))
\ No newline at end of file
diff --git a/src/kernelfunctions/mercer/polynomial.jl b/src/kernelfunctions/mercer/polynomial.jl
index 737ca76..530154a 100644
--- a/src/kernelfunctions/mercer/polynomial.jl
+++ b/src/kernelfunctions/mercer/polynomial.jl
@@ -3,32 +3,54 @@
 
 The polynomial kernel is a Mercer kernel given by:
 
-```math
-\kappa(\mathbf{x},\mathbf{y}) = 
-(a \mathbf{x}^\intercal \mathbf{y} + c)^d
-\qquad \alpha > 0, \; c \geq 0, \; d \in \mathbb{Z}_{+}
+```
+    κ(x,y) = (a⋅xᵀy + c)ᵈ   α > 0, c ≧ 0, d ∈ ℤ⁺
+```
+
+# Examples
+
+```jldoctest; setup = :(using MLKernels)
+julia> PolynomialKernel(2.0f0)
+PolynomialKernel{Float32}(2.0,1.0,3)
+
+julia> PolynomialKernel(2.0f0, 2.0)
+PolynomialKernel{Float64}(2.0,2.0,3)
+
+julia> PolynomialKernel(2.0f0, 2.0, 2)
+PolynomialKernel{Float64}(2.0,2.0,2)
 ```
 """
-struct PolynomialKernel{T<:AbstractFloat,U<:Integer} <: MercerKernel{T}
-    a::HyperParameter{T}
-    c::HyperParameter{T}
-    d::HyperParameter{U}
-    function PolynomialKernel{T}(a::Real, c::Real, d::U) where {T<:AbstractFloat,U<:Integer}
-        new{T,U}(HyperParameter(convert(T,a), interval(OpenBound(zero(T)), nothing)),
-                 HyperParameter(convert(T,c), interval(ClosedBound(zero(T)), nothing)),
-                 HyperParameter(d, interval(ClosedBound(one(U)), nothing)))
+struct PolynomialKernel{T<:AbstractFloat} <: MercerKernel{T}
+    a::T
+    c::T
+    d::T
+    function PolynomialKernel{T}(
+            a::Real=T(1),
+            c::Real=T(1),
+            d::Real=T(3)
+        ) where {T<:AbstractFloat}
+        @check_args(PolynomialKernel, a, a >  zero(a), "a > 0")
+        @check_args(PolynomialKernel, c, c >= zero(c), "c ≧ 0")
+        @check_args(PolynomialKernel, d, d >= one(d) && d == trunc(d), "d ∈ ℤ₊")
+        return new{T}(a, c, d)
     end
 end
-function PolynomialKernel(a::T1=1.0, c::T2=one(T1), d::Integer=3) where {T1<:Real,T2<:Real}
-    PolynomialKernel{floattype(T1,T2)}(a, c, d)
+
+function PolynomialKernel(
+        a::T₁=1.0,
+        c::T₂=T₁(1),
+        d::T₃=convert(promote_float(T₁,T₂), 3)
+    ) where {T₁<:Real,T₂<:Real,T₃<:Real}
+    T = promote_float(T₁,T₂,T₃)
+    return PolynomialKernel{T}(a, c, d)
 end
 
-@inline eltypes(::Type{<:PolynomialKernel{T,U}}) where {T,U} = (T,U)
-@inline thetafieldnames(κ::PolynomialKernel) = Symbol[:a, :c]
+@inline basefunction(::PolynomialKernel) = ScalarProduct()
 
-@inline polynomialkernel(z::T, a::T, c::T, d::U) where {T<:AbstractFloat,U<:Integer} = (a*z + c)^d
+@inline function kappa(κ::PolynomialKernel{T}, xᵀy::T) where {T}
+    return (κ.a*xᵀy + κ.c)^(κ.d)
+end
 
-@inline basefunction(::PolynomialKernel) = ScalarProduct()
-@inline function kappa(κ::PolynomialKernel{T}, z::T) where {T}
-    polynomialkernel(z, getvalue(κ.a), getvalue(κ.c), getvalue(κ.d))
+function convert(::Type{K}, κ::PolynomialKernel) where {K>:PolynomialKernel{T}} where T
+    return PolynomialKernel{T}(κ.a, κ.c, κ.d)
 end
\ No newline at end of file
diff --git a/src/kernelfunctions/mercer/rationalquadratic.jl b/src/kernelfunctions/mercer/rationalquadratic.jl
index 47690d6..9e1ecc4 100644
--- a/src/kernelfunctions/mercer/rationalquadratic.jl
+++ b/src/kernelfunctions/mercer/rationalquadratic.jl
@@ -1,34 +1,125 @@
+# Abstract Rational-Quadratic Kernel =======================================================
+abstract type AbstractRationalQuadraticKernel{T<:AbstractFloat} <: MercerKernel{T} end
+
+@inline basefunction(::AbstractRationalQuadraticKernel) = SquaredEuclidean()
+
+
+# Rational-Quadratic Kernel ================================================================
 @doc raw"""
-    RationalQuadraticKernel([α=1 [,β=1]])
+    RationalQuadraticKernel([α [,β]])
 
-The rational-quadratic kernel is given by:
+The rational quadratic kernel is an isotropic Mercer kernel given by the formula:
 
-```math
-\kappa(\mathbf{x},\mathbf{y}) 
-= \left(1 +\alpha ||\mathbf{x},\mathbf{y}||^2\right)^{-\beta} 
-\qquad \alpha > 0, \; \beta > 0
 ```
+    κ(x,y) = (1 + α‖x-y‖²)⁻ᵝ   α > 0, β > 0
+```
+where `α` is a scaling parameter and `β` is a shape parameter. The rational quadratic
+kernel is a special of the more general gamma-rational-quadratic kernel (see
+[`GammaRationalQuadraticKernel`](@ref)) with `γ = 1`.
 
-where ``\alpha`` is a scaling parameter and ``\beta`` is a shape parameter. This kernel can 
-be seen as an infinite sum of Gaussian kernels. If one sets ``\alpha = \alpha_0 / \beta``, 
-then taking the limit ``\beta \rightarrow \infty`` results in the Gaussian kernel with 
-scaling parameter ``\alpha_0``.
+# Examples
+
+```jldoctest; setup = :(using MLKernels)
+julia> RationalQuadraticKernel()
+RationalQuadraticKernel{Float64}(1.0,1.0)
+
+julia> RationalQuadraticKernel(2.0f0)
+RationalQuadraticKernel{Float32}(2.0,1.0)
+
+julia> RationalQuadraticKernel(2.0f0, 2.0)
+RationalQuadraticKernel{Float64}(2.0,2.0)
+```
 """
-struct RationalQuadraticKernel{T<:AbstractFloat} <: MercerKernel{T}
-    alpha::HyperParameter{T}
-    beta::HyperParameter{T}
-    RationalQuadraticKernel{T}(α::Real, β::Real) where {T<:AbstractFloat} = new{T}(
-        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing)),
-        HyperParameter(convert(T,β), interval(OpenBound(zero(T)), nothing))
-    )
+struct RationalQuadraticKernel{T<:AbstractFloat} <: AbstractRationalQuadraticKernel{T}
+    α::T
+    β::T
+    function RationalQuadraticKernel{T}(
+            α::Real=T(1),
+            β::Real=T(1)
+        ) where {T<:AbstractFloat}
+        @check_args(RationalQuadraticKernel, α, α > zero(T), "α > 0")
+        @check_args(RationalQuadraticKernel, β, β > zero(T), "β > 0")
+        return new{T}(α, β)
+    end
 end
-function RationalQuadraticKernel(α::T1 = 1.0, β::T2 = one(T1)) where {T1<:Real,T2<:Real}
-    RationalQuadraticKernel{floattype(T1,T2)}(α, β)
+function RationalQuadraticKernel(
+        α::T₁=1.0,
+        β::T₂=T₁(1)
+    ) where {T₁<:Real, T₂<:Real}
+    RationalQuadraticKernel{promote_float(T₁, T₂)}(α, β)
 end
 
-@inline rationalquadratickernel(z::T, α::T, β::T) where {T<:AbstractFloat} = (1 + α*z)^(-β)
+@inline function kappa(κ::RationalQuadraticKernel{T}, d²::T) where {T}
+    return (one(T) + κ.α*d²)^(-κ.β)
+end
+
+function convert(
+        ::Type{K},
+        κ::RationalQuadraticKernel
+    ) where {K>:RationalQuadraticKernel{T}} where T
+    return RationalQuadraticKernel{T}(κ.α, κ.β)
+end
+
+
+# Gamma Rational-Quadratic Kernel ==========================================================
+@doc raw"""
+    GammaRationalKernel([α [,β [,γ]]])
+
+The gamma-rational-quadratic kernel is a generalization of the rational-quadratic kernel
+with an additional shape parameter `γ`:
+
+```
+    κ(x,y) = (1 + α‖x-y‖²ᵞ)⁻ᵝ   α > 0, β > 0, γ ∈ (0,1]
+```
+where ``\alpha`` is a scaling parameter, ``\beta`` is a shape parameter and ``\gamma`` is a
+shape parameter of the Euclidean distance.
+
+# Examples
+
+```jldoctest; setup = :(using MLKernels)
+julia> GammaRationalQuadraticKernel()
+GammaRationalQuadraticKernel{Float64}(1.0,1.0,1.0)
+
+julia> GammaRationalQuadraticKernel(2.0f0)
+GammaRationalQuadraticKernel{Float32}(2.0,1.0,1.0)
+
+julia> GammaRationalQuadraticKernel(2.0f0, 2.0f0)
+GammaRationalQuadraticKernel{Float32}(2.0,2.0,1.0)
+
+julia> GammaRationalQuadraticKernel(2.0f0, 2.0f0, 0.5f0)
+GammaRationalQuadraticKernel{Float32}(2.0,2.0,0.5)
+```
+"""
+struct GammaRationalQuadraticKernel{T<:AbstractFloat} <: AbstractRationalQuadraticKernel{T}
+    α::T
+    β::T
+    γ::T
+    function GammaRationalQuadraticKernel{T}(
+            α::Real=T(1),
+            β::Real=T(1),
+            γ::Real=T(1)
+        ) where {T<:AbstractFloat}
+        @check_args(GammaRationalQuadraticKernel, α, α > zero(T), "α > 0")
+        @check_args(GammaRationalQuadraticKernel, β, β > zero(T), "β > 0")
+        @check_args(GammaRationalQuadraticKernel, γ, one(T) >= γ > zero(T), "γ ∈ (0,1]")
+        return new{T}(α, β, γ)
+    end
+end
+function GammaRationalQuadraticKernel(
+        α::T₁ = 1.0,
+        β::T₂ = T₁(1),
+        γ::T₃ = one(promote_float(T₁, T₂))
+    ) where {T₁<:Real, T₂<:Real, T₃<:Real}
+    GammaRationalQuadraticKernel{promote_float(T₁,T₂,T₃)}(α, β, γ)
+end
+
+@inline function kappa(κ::GammaRationalQuadraticKernel{T}, d²::T) where {T}
+    return (one(T) + κ.α*(d²^κ.γ))^(-κ.β)
+end
 
-@inline basefunction(::RationalQuadraticKernel) = SquaredEuclidean()
-@inline function kappa(κ::RationalQuadraticKernel{T}, z::T) where {T}
-    rationalquadratickernel(z, getvalue(κ.alpha), getvalue(κ.beta))
+function convert(
+        ::Type{K},
+        κ::GammaRationalQuadraticKernel
+    ) where {K>:GammaRationalQuadraticKernel{T}} where T
+    return GammaRationalQuadraticKernel{T}(κ.α, κ.β, κ.γ)
 end
\ No newline at end of file
diff --git a/src/kernelfunctions/mercer/squaredexponential.jl b/src/kernelfunctions/mercer/squaredexponential.jl
deleted file mode 100644
index 8685968..0000000
--- a/src/kernelfunctions/mercer/squaredexponential.jl
+++ /dev/null
@@ -1,32 +0,0 @@
-@doc raw"""
-    SquaredExponentialKernel([α=1])
-
-The squared exponential kernel, or alternatively the Gaussian kernel, is identical to the 
-exponential kernel except that the Euclidean distance is squared:
-
-```math
-\kappa(\mathbf{x},\mathbf{y}) = \exp\left(-\alpha ||\mathbf{x} - \mathbf{y}||^2\right) 
-\qquad \alpha > 0
-```
-
-where ``\alpha`` is a scaling parameter of the squared Euclidean distance.
-Just like the exponential kernel, the squared exponential kernel is an
-isotropic Mercer kernel. The squared exponential kernel is more commonly known
-as the radial basis kernel within machine learning communities.
-"""
-struct SquaredExponentialKernel{T<:AbstractFloat} <: MercerKernel{T}
-    alpha::HyperParameter{T}
-    SquaredExponentialKernel{T}(α::Real) where {T<:AbstractFloat} = new{T}(
-        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing))
-    )
-end
-SquaredExponentialKernel(α::T=1.0) where {T<:Real} = SquaredExponentialKernel{floattype(T)}(α)
-GaussianKernel = SquaredExponentialKernel
-RadialBasisKernel = SquaredExponentialKernel
-
-@inline squaredexponentialkernel(z::T, α::T) where {T<:AbstractFloat} = exp(-α*z)
-
-@inline basefunction(::SquaredExponentialKernel) = SquaredEuclidean()
-@inline function kappa(κ::SquaredExponentialKernel{T}, z::T) where {T}
-    squaredexponentialkernel(z, getvalue(κ.alpha))
-end
\ No newline at end of file
diff --git a/src/kernelfunctions/negativedefinite/log.jl b/src/kernelfunctions/negativedefinite/log.jl
index 699c896..d075984 100644
--- a/src/kernelfunctions/negativedefinite/log.jl
+++ b/src/kernelfunctions/negativedefinite/log.jl
@@ -1,29 +1,43 @@
 @doc raw"""
     LogKernel([α [,γ]])
 
-The Log Kernel is a negative definite kernel given by:
+The Log Kernel is a negative definite kernel given by the formula:
 
-```math
-\kappa(\mathbf{x},\mathbf{y}) = 
-\log \left(1 + \alpha\|\mathbf{x} - \mathbf{y} \|^{2\gamma}\right)
-\qquad \alpha > 0, \; \gamma \in (0,1]
+```
+    κ(x,y) = log(1 + α‖x-y‖²ᵞ)   α > 0, γ ∈ (0,1]
+```
+where `α` is a scaling parameter and `γ` is a shape parameter of the Euclidean distance.
+
+# Examples
+
+```jldoctest; setup = :(using MLKernels)
+julia> LogKernel()
+LogKernel{Float64}(1.0,1.0)
+
+julia> LogKernel(0.5f0)
+LogKernel{Float32}(0.5,1.0)
+
+julia> LogKernel(0.5, 0.5)
+LogKernel{Float64}(0.5,0.5)
 ```
 """
 struct LogKernel{T<:AbstractFloat} <: NegativeDefiniteKernel{T}
-    alpha::HyperParameter{T}
-    gamma::HyperParameter{T}
-    LogKernel{T}(α::Real, γ::Real) where {T<:AbstractFloat} = new{T}(
-        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing)),
-        HyperParameter(convert(T,γ), interval(OpenBound(zero(T)), ClosedBound(one(T))))
-    )
+    α::T
+    γ::T
+    function LogKernel{T}(α::Real=T(1), γ::Real=T(1)) where {T<:AbstractFloat}
+        @check_args(LogKernel, α, α > zero(T), "α > 0")
+        @check_args(LogKernel, γ, one(T) >= γ > zero(T), "γ ∈ (0,1]")
+        return new{T}(α, γ)
+    end
 end
-function LogKernel(α::T1 = 1.0, γ::T2 = one(T1)) where {T1<:Real,T2<:Real}
-    LogKernel{floattype(T1,T2)}(α, γ)
+function LogKernel(α::T₁=1.0, γ::T₂=T₁(1)) where {T₁<:Real,T₂<:Real}
+    LogKernel{promote_float(T₁,T₂)}(α, γ)
 end
 
-@inline logkernel(z::T, α::T, γ::T) where {T<:AbstractFloat} = log(α*z^γ+1)
-
 @inline basefunction(::LogKernel) = SquaredEuclidean()
-@inline function kappa(κ::LogKernel{T}, z::T) where {T}
-    logkernel(z, getvalue(κ.alpha), getvalue(κ.gamma))
+
+@inline kappa(κ::LogKernel{T}, d²::T) where {T} = log(one(T) + κ.α*d²^(κ.γ))
+
+function convert(::Type{K}, κ::LogKernel) where {K>:LogKernel{T}} where T
+    return LogKernel{T}(κ.α, κ.γ)
 end
\ No newline at end of file
diff --git a/src/kernelfunctions/negativedefinite/power.jl b/src/kernelfunctions/negativedefinite/power.jl
index 0cdb86d..f37a9b3 100644
--- a/src/kernelfunctions/negativedefinite/power.jl
+++ b/src/kernelfunctions/negativedefinite/power.jl
@@ -2,22 +2,34 @@
     PowerKernel([γ=1])
 
 The Power Kernel is a negative definite kernel given by:
+```
+    κ(x,y) = ‖x-y‖²ᵞ   γ ∈ (0,1]
+```
+where `γ` is a shape parameter of the Euclidean distance.
+
+# Examples
+
+```jldoctest; setup = :(using MLKernels)
+julia> PowerKernel()
+PowerKernel{Float64}(1.0)
 
-```math
-\kappa(\mathbf{x},\mathbf{y}) = 
-\|\mathbf{x} - \mathbf{y} \|^{2\gamma}
-\qquad \gamma \in (0,1]
+julia> PowerKernel(0.5f0)
+PowerKernel{Float32}(0.5)
 ```
 """
 struct PowerKernel{T<:AbstractFloat} <: NegativeDefiniteKernel{T}
-    gamma::HyperParameter{T}
-    PowerKernel{T}(γ::Real) where {T<:AbstractFloat} = new{T}(
-        HyperParameter(convert(T,γ), interval(OpenBound(zero(T)), ClosedBound(one(T))))
-    )
+    γ::T
+    function PowerKernel{T}(γ::Real=T(1)) where {T<:AbstractFloat}
+        @check_args(PowerKernel, γ, one(T) >= γ > zero(T), "γ ∈ (0,1]")
+        new{T}(γ)
+    end
 end
-PowerKernel(γ::T1 = 1.0) where {T1<:Real} = PowerKernel{floattype(T1)}(γ)
-
-@inline powerkernel(z::T, γ::T) where {T<:AbstractFloat} = z^γ
+PowerKernel(γ::T=1.0) where {T<:Real} = PowerKernel{promote_float(T)}(γ)
 
 @inline basefunction(::PowerKernel) = SquaredEuclidean()
-@inline kappa(κ::PowerKernel{T}, z::T) where {T} = powerkernel(z, getvalue(κ.gamma))
\ No newline at end of file
+
+@inline kappa(κ::PowerKernel{T}, d²::T) where {T} = d²^κ.γ
+
+function convert(::Type{K}, κ::PowerKernel) where {K>:PowerKernel{T}} where T
+    return PowerKernel{T}(κ.γ)
+end
\ No newline at end of file
diff --git a/src/kernelfunctions/sigmoid.jl b/src/kernelfunctions/sigmoid.jl
index 44c1055..3a83f76 100644
--- a/src/kernelfunctions/sigmoid.jl
+++ b/src/kernelfunctions/sigmoid.jl
@@ -1,30 +1,41 @@
 @doc raw"""
     SigmoidKernel([a=1 [,c=1]])
-    
-The Sigmoid Kernel is given by
 
-```math
-\kappa(\mathbf{x},\mathbf{y}) = 
-\tanh(a \mathbf{x}^\intercal \mathbf{y} + c) 
-\qquad \alpha > 0, \; c \geq 0
+The Sigmoid Kernel is given by:
+```
+    κ(x,y) = tanh(a⋅xᵀy + c)
+```
+
+# Examples
+
+```jldoctest; setup = :(using MLKernels)
+julia> SigmoidKernel()
+SigmoidKernel{Float64}(1.0,1.0)
+
+julia> SigmoidKernel(0.5f0)
+SigmoidKernel{Float32}(0.5,1.0)
+
+julia> SigmoidKernel(0.5f0, 0.5)
+SigmoidKernel{Float64}(0.5,0.5)
 ```
-The sigmoid kernel is a not a true kernel, although it has been used in application. 
 """
 struct SigmoidKernel{T<:AbstractFloat} <: Kernel{T}
-    a::HyperParameter{T}
-    c::HyperParameter{T}
-    SigmoidKernel{T}(a::Real, c::Real) where {T<:AbstractFloat} = new{T}(
-        HyperParameter(convert(T,a), interval(OpenBound(zero(T)),   nothing)),
-        HyperParameter(convert(T,c), interval(ClosedBound(zero(T)), nothing))
-    )
+    a::T
+    c::T
+    function SigmoidKernel{T}(a::Real=T(1), c::Real=T(1)) where {T<:AbstractFloat}
+        @check_args(SigmoidKernel, a, a >  zero(T), "a > 0")
+        @check_args(SigmoidKernel, c, c >= zero(T), "c ≧ 0")
+        return new{T}(a, c)
+    end
 end
-function SigmoidKernel(a::T1 = 1.0, c::T2 = one(T1)) where {T1<:Real,T2<:Real}
-    SigmoidKernel{floattype(T1,T2)}(a,c)
+function SigmoidKernel(a::T₁=1.0, c::T₂=T₁(1)) where {T₁<:Real,T₂<:Real}
+    SigmoidKernel{promote_float(T₁,T₂)}(a,c)
 end
 
-@inline sigmoidkernel(z::T, a::T, c::T) where {T<:AbstractFloat} = tanh(a*z + c)
-
 @inline basefunction(::SigmoidKernel) = ScalarProduct()
-@inline function kappa(κ::SigmoidKernel{T}, z::T) where {T}
-    sigmoidkernel(z, getvalue(κ.a), getvalue(κ.c))
+
+@inline kappa(κ::SigmoidKernel{T}, xᵀy::T) where {T} = tanh(κ.a*xᵀy + κ.c)
+
+function convert(::Type{K}, κ::SigmoidKernel) where {K>:SigmoidKernel{T}} where T
+    return SigmoidKernel{T}(κ.a, κ.c)
 end
\ No newline at end of file
diff --git a/src/kernelmatrix.jl b/src/kernelmatrix.jl
index 43546e9..9205b94 100644
--- a/src/kernelmatrix.jl
+++ b/src/kernelmatrix.jl
@@ -37,37 +37,37 @@ function symmetric_kappamatrix!(
 end
 
 """
-    kernelmatrix!(P::Matrix, σ::Orientation, κ::Kernel, X::Matrix, symmetrize::Bool)
+    kernelmatrix!(σ::Orientation, K::Matrix, κ::Kernel, X::Matrix, symmetrize::Bool)
 
 In-place version of `kernelmatrix` where pre-allocated matrix `K` will be overwritten
 with the kernel matrix.
 """
 function kernelmatrix!(
         σ::Orientation,
-        P::Matrix{T},
+        K::Matrix{T},
         κ::Kernel{T},
         X::AbstractMatrix{T},
         symmetrize::Bool
     ) where {T<:AbstractFloat}
-    basematrix!(σ, P, basefunction(κ), X, false)
-    symmetric_kappamatrix!(κ, P, symmetrize)
+    basematrix!(σ, K, basefunction(κ), X, false)
+    symmetric_kappamatrix!(κ, K, symmetrize)
 end
 
 """
-    kernelmatrix!(K::Matrix, σ::Orientation, κ::Kernel, X::Matrix, Y::Matrix)
+    kernelmatrix!(σ::Orientation, K::Matrix, κ::Kernel, X::Matrix, Y::Matrix)
 
 In-place version of `kernelmatrix` where pre-allocated matrix `K` will be overwritten with
 the kernel matrix.
 """
 function kernelmatrix!(
         σ::Orientation,
-        P::Matrix{T},
+        K::Matrix{T},
         κ::Kernel{T},
         X::AbstractMatrix{T},
         Y::AbstractMatrix{T}
     ) where {T<:AbstractFloat}
-    basematrix!(σ, P, basefunction(κ), X, Y)
-    kappamatrix!(κ, P)
+    basematrix!(σ, K, basefunction(κ), X, Y)
+    kappamatrix!(κ, K)
 end
 
 function kernelmatrix(
@@ -76,7 +76,8 @@ function kernelmatrix(
         X::AbstractMatrix{T},
         symmetrize::Bool = true
     ) where {T<:AbstractFloat}
-    symmetric_kappamatrix!(κ, basematrix(σ, basefunction(κ), X, false), symmetrize)
+    K = allocate_basematrix(σ, X)
+    symmetric_kappamatrix!(κ, basematrix!(σ, K, basefunction(κ), X, false), symmetrize)
 end
 
 function kernelmatrix(
@@ -85,7 +86,8 @@ function kernelmatrix(
         X::AbstractMatrix{T},
         Y::AbstractMatrix{T}
     ) where {T<:AbstractFloat}
-    kappamatrix!(κ, basematrix(σ, basefunction(κ), X, Y))
+    K = allocate_basematrix(σ, X, Y)
+    kappamatrix!(κ, basematrix!(σ, K, basefunction(κ), X, Y))
 end
 
 
diff --git a/src/utils.jl b/src/utils.jl
index 903058c..7df545c 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,3 +1,40 @@
+# Check Arguments Macro ====================================================================
+
+macro check_args(K, param, cond, desc=string(cond))
+    quote
+        if !($(esc(cond)))
+            throw(ArgumentError(string(
+                $(string(K)), ": ", $(string(param)), " = ", $(esc(param)), " does not ",
+                "satisfy the constraint ", $(string(desc)), ".")))
+        end
+    end
+end
+
+#macro default_kernel_convert(K)
+#    θ = fieldnames(eval(K))
+#    K_name = esc(K)
+#    conversion_call = Expr(:call, :($K_name{T}), [:(κ.$θₖ) for θₖ in θ]...)
+#    return quote
+#        function test(::Type{K}, κ::$K_name) where {K>:$K_name{T}} where T
+#            return $conversion_call
+#        end
+#    end
+#end
+
+
+# Type Rules ===============================================================================
+
+function promote_float(Tₖ::DataType...)
+    if length(Tₖ) == 0
+        return Float64
+    end
+    T = promote_type(Tₖ...)
+    return T <: AbstractFloat ? T : Float64
+end
+
+
+# Common Functions =========================================================================
+
 for orientation in (:row, :col)
 
     row_oriented = orientation == :row
diff --git a/test/HyperParameters/hyperparameter.jl b/test/HyperParameters/hyperparameter.jl
deleted file mode 100644
index 586cc45..0000000
--- a/test/HyperParameters/hyperparameter.jl
+++ /dev/null
@@ -1,457 +0,0 @@
-@testset "Testing MLK.OpenBound" begin
-    for T in (FloatingPointTypes..., IntegerTypes...)
-        if T <: Integer
-            @test_throws ErrorException MLK.OpenBound(one(T))
-        else
-            if T <: AbstractFloat
-                @test_throws ErrorException MLK.OpenBound(convert(T,  NaN))
-                @test_throws ErrorException MLK.OpenBound(convert(T,  Inf))
-                @test_throws ErrorException MLK.OpenBound(convert(T, -Inf))
-            end
-
-            for x in (1,2)
-                B = MLK.OpenBound(convert(T,x))
-
-                @test B.value  == x
-                @test eltype(B) == T
-            end
-
-            for U in FloatingPointTypes
-                B = MLK.OpenBound(one(T))
-
-                B_u = convert(MLK.OpenBound{U}, B)
-                @test B_u.value == one(U)
-                @test eltype(B_u) == U
-            end
-
-            B = MLK.OpenBound(zero(T))
-
-            @test checkvalue(-one(T), B) == true
-            @test checkvalue(zero(T), B) == false
-            @test checkvalue(one(T),  B) == false
-
-            @test checkvalue(B, -one(T)) == false
-            @test checkvalue(B, zero(T)) == false
-            @test checkvalue(B, one(T))  == true
-
-            @test string(B) == string("OpenBound(", zero(T), ")")
-            @test show(devnull, B) == nothing
-        end
-    end
-end
-
-@testset "Testing MLK.ClosedBound" begin
-    for T in (FloatingPointTypes..., IntegerTypes...)
-        if T <: AbstractFloat
-            @test_throws ErrorException MLK.ClosedBound(convert(T,  NaN))
-            @test_throws ErrorException MLK.ClosedBound(convert(T,  Inf))
-            @test_throws ErrorException MLK.ClosedBound(convert(T, -Inf))
-        end
-
-        for x in (1,2)
-            B = MLK.ClosedBound(convert(T,x))
-
-            @test B.value  == x
-            @test eltype(B) == T
-        end
-
-        for U in (FloatingPointTypes..., IntegerTypes...)
-            B = MLK.ClosedBound(one(T))
-
-            B_u = convert(MLK.ClosedBound{U}, B)
-            @test B_u.value == one(U)
-            @test eltype(B_u) == U
-        end
-
-        B = MLK.ClosedBound(zero(T))
-
-        @test checkvalue(-one(T), B) == true
-        @test checkvalue(zero(T), B) == true
-        @test checkvalue(one(T),  B) == false
-
-        @test checkvalue(B, -one(T)) == false
-        @test checkvalue(B, zero(T)) == true
-        @test checkvalue(B, one(T))  == true
-
-        @test string(B) == string("ClosedBound(", zero(T), ")")
-        @test show(devnull, B) == nothing
-    end
-end
-
-@testset "Testing MLK.NullBound" begin
-    for T in (FloatingPointTypes..., IntegerTypes...)
-        @test eltype(NullBound(T)) == T
-
-        for U in (FloatingPointTypes..., IntegerTypes...)
-            B = MLK.NullBound(T)
-
-            B_u = convert(MLK.NullBound{U}, B)
-            @test eltype(B_u) == U
-        end
-
-        B = MLK.NullBound(T)
-
-        @test checkvalue(-one(T), B) == true
-        @test checkvalue(zero(T), B) == true
-        @test checkvalue(one(T),  B) == true
-
-        @test checkvalue(B, -one(T)) == true
-        @test checkvalue(B, zero(T)) == true
-        @test checkvalue(B, one(T))  == true
-
-        @test string(B) == string("NullBound(", T, ")")
-        @test show(devnull, B) == nothing
-    end
-end
-
-@testset "Testing MLK.Interval" begin
-    for T in FloatingPointTypes
-        @test_throws ErrorException MLK.Interval(MLK.ClosedBound(one(T)), MLK.OpenBound(one(T)))
-        @test_throws ErrorException MLK.Interval(MLK.OpenBound(one(T)),   MLK.ClosedBound(one(T)))
-        @test_throws ErrorException MLK.Interval(MLK.OpenBound(one(T)),   MLK.OpenBound(one(T)))
-
-        a = MLK.ClosedBound(one(T))
-        b = MLK.ClosedBound(one(T))
-        I = MLK.Interval(a,b)
-        @test I.a == a
-        @test I.b == b
-        @test eltype(I) == T
-
-        I = MLK.interval(nothing, nothing)
-        @test I.a == NullBound(Float64)
-        @test I.b == NullBound(Float64)
-        @test eltype(I) == Float64
-
-        for a in (NullBound(T), ClosedBound(-one(T)), OpenBound(-one(T)))
-            for b in (NullBound(T), ClosedBound(one(T)), OpenBound(one(T)))
-                I = MLK.Interval(a, b)
-                @test I.a == a
-                @test I.b == b
-                @test eltype(I) == T
-
-                if typeof(a) <: NullBound
-                    if typeof(b) <: NullBound
-                        @test I == MLK.interval(T)
-                        @test string(I) == string("interval(", T, ")")
-                    else
-                        @test I == MLK.interval(nothing, b)
-                        @test string(I) == string("interval(nothing,", string(b), ")")
-                    end
-                else
-                    if typeof(b) <: NullBound
-                        @test I == MLK.interval(a,nothing)
-                        @test string(I) == string("interval(", string(a), ",nothing)")
-                    else
-                        @test I == MLK.interval(a, b)
-                        @test string(I) == string("interval(", string(a), ",", string(b), ")")
-                    end
-                end
-
-                @test MLK.checkvalue(I, convert(T,-2)) == (typeof(a) <: NullBound ? true  : false)
-                @test MLK.checkvalue(I, -one(T))       == (typeof(a) <: OpenBound ? false : true)
-                @test MLK.checkvalue(I, zero(T))       == true
-                @test MLK.checkvalue(I,  one(T))       == (typeof(b) <: OpenBound ? false : true)
-                @test MLK.checkvalue(I, convert(T,2))  == (typeof(b) <: NullBound ? true  : false)
-            end
-        end
-
-        a = convert(T,7.6)
-        b = convert(T,23)
-        c = convert(T,13)
-
-        @test show(devnull, MLK.interval(MLK.ClosedBound(one(T)), MLK.ClosedBound(one(T)))) == nothing
-    end
-end
-
-@testset "Testing MLK.interval" begin
-    @test typeof(MLK.interval(nothing, nothing)) == MLK.Interval{Float64,NullBound{Float64},NullBound{Float64}}
-    for T in FloatingPointTypes
-        for a in (NullBound(T), ClosedBound(-one(T)), OpenBound(-one(T)))
-            for b in (NullBound(T), ClosedBound(one(T)), OpenBound(one(T)))
-                null_a = typeof(a) <: NullBound
-                null_b = typeof(b) <: NullBound
-
-                I = interval(null_a ? nothing : a, null_b ? nothing : b)
-
-                if null_a && null_b
-                    @test typeof(I) == MLK.Interval{Float64,NullBound{Float64},NullBound{Float64}}
-                else
-                    @test I.a == a
-                    @test I.b == b
-                    @test eltype(T) == T
-                end
-            end
-        end
-    end
-end
-
-
-@testset "Testing MLK.checkvalue" begin
-    for T in FloatingPointTypes
-        l = convert(T,-9)
-        a = convert(T,-3)
-        b = convert(T,3)
-        u = convert(T,9)
-        for A in (NullBound(T), ClosedBound(a), OpenBound(a))
-            T_a = typeof(A)
-            for B in (NullBound(T), ClosedBound(b), OpenBound(b))
-                T_b = typeof(B)
-                I = HP.interval(A,B)
-
-                for x in range(l, stop = a, length = 30)
-                    if T_a <: NullBound || T_a <: ClosedBound && x == a
-                        @test MLK.checkvalue(I,x) == true
-                    else
-                        @test MLK.checkvalue(I,x) == false
-                    end
-                end
-
-                for x in range(a, stop = b, length = 30)
-                    if x == a
-                        @test MLK.checkvalue(I,x) == T_a <: OpenBound ? false : true
-                    elseif x == b
-                        @test MLK.checkvalue(I,x) == T_b <: OpenBound ? false : true
-                    else
-                        @test MLK.checkvalue(I,x) == true
-                    end
-                end
-
-                for x in range(b, stop = u, length = 30)
-                    if T_b <: NullBound || T_b <: ClosedBound && x == b
-                        @test MLK.checkvalue(I,x) == true
-                    else
-                        @test MLK.checkvalue(I,x) == false
-                    end
-                end
-            end
-        end
-
-    end
-end
-
-
-@testset "Testing MLK.lowerboundtheta" begin
-    for T in FloatingPointTypes
-        for A in (NullBound(T), ClosedBound(-one(T)), OpenBound(-one(T)))
-            T_a = typeof(A)
-            for B in (NullBound(T), ClosedBound(one(T)), OpenBound(one(T)))
-                T_b = typeof(B)
-                I = MLK.Interval(A, B)
-
-                if T_a <: NullBound
-                    if T_b <: NullBound
-                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
-                    elseif T_b <: ClosedBound
-                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
-                    else
-                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
-                    end
-                elseif T_a <: ClosedBound
-                    if T_b <: NullBound
-                        @test MLK.lowerboundtheta(I) == I.a.value
-                    elseif T_b <: ClosedBound
-                        @test MLK.lowerboundtheta(I) == I.a.value
-                    else
-                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
-                    end
-                else
-                    if T_b <: NullBound
-                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
-                    elseif T_b <: ClosedBound
-                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
-                    else
-                        @test MLK.lowerboundtheta(I) == convert(T,-Inf)
-                    end
-                end
-            end
-        end
-    end
-end
-
-@testset "Testing MLK.upperboundtheta" begin
-    for T in FloatingPointTypes
-        a = convert(T,-5)
-        b = convert(T,5)
-        for A in (NullBound(T), ClosedBound(a), OpenBound(a))
-            T_a = typeof(A)
-            for B in (NullBound(T), ClosedBound(b), OpenBound(b))
-                T_b = typeof(B)
-                I = MLK.Interval(A, B)
-
-                if T_a <: NullBound
-                    if T_b <: NullBound
-                        @test MLK.upperboundtheta(I) == convert(T,Inf)
-                    elseif T_b <: ClosedBound
-                        @test MLK.upperboundtheta(I) == b
-                    else
-                        @test MLK.upperboundtheta(I) == convert(T,Inf)
-                    end
-                elseif T_a <: ClosedBound
-                    if T_b <: NullBound
-                        @test MLK.upperboundtheta(I) == convert(T,Inf)
-                    elseif T_b <: ClosedBound
-                        @test MLK.upperboundtheta(I) == b
-                    else
-                        @test MLK.upperboundtheta(I) == log(b - a)
-                    end
-                else
-                    if T_b <: NullBound
-                        @test MLK.upperboundtheta(I) == convert(T,Inf)
-                    elseif T_b <: ClosedBound
-                        @test MLK.upperboundtheta(I) == log(b - a)
-                    else
-                        @test MLK.upperboundtheta(I) == convert(T,Inf)
-                    end
-                end
-            end
-        end
-    end
-end
-
-
-@testset "Testing HP.theta" begin
-    for T in FloatingPointTypes
-        l = convert(T,-9)
-        a = convert(T,-3)
-        b = convert(T,3)
-        u = convert(T,9)
-        for A in (NullBound(T), ClosedBound(a), OpenBound(a))
-            T_a = typeof(A)
-            for B in (NullBound(T), ClosedBound(b), OpenBound(b))
-                T_b = typeof(B)
-                I = MLK.Interval(A, B)
-
-                for x in range(l, stop = u, length = 60)
-                    if T_a <: NullBound
-                        if T_b <: NullBound
-                            @test HP.theta(I,x) == x
-                        elseif T_b <: ClosedBound
-                            if x <= b
-                                @test HP.theta(I,x) == x
-                            else
-                                @test_throws DomainError HP.theta(I,x)
-                            end
-                        else
-                            if x < b
-                                @test HP.theta(I,x) == log(b-x)
-                            else
-                                @test_throws DomainError HP.theta(I,x)
-                            end
-                        end
-                    elseif T_a <: ClosedBound
-                        if T_b <: NullBound
-                            if a <= x
-                                @test HP.theta(I,x) == x
-                            else
-                                @test_throws DomainError HP.theta(I,x)
-                            end
-                        elseif T_b <: ClosedBound
-                            if a <= x <= b
-                                @test HP.theta(I,x) == x
-                            else
-                                @test_throws DomainError HP.theta(I,x)
-                            end
-                        else
-                            if a <= x < b
-                                @test HP.theta(I,x) == log(b-x)
-                            else
-                                @test_throws DomainError HP.theta(I,x)
-                            end
-                        end
-                    else
-                        if T_b <: NullBound
-                            if a < x
-                                @test HP.theta(I,x) == log(x-a)
-                            else
-                                @test_throws DomainError HP.theta(I,x)
-                            end
-                        elseif T_b <: ClosedBound
-                            if a < x <= b
-                                @test HP.theta(I,x) == log(x-a)
-                            else
-                                @test_throws DomainError HP.theta(I,x)
-                            end
-                        else
-                            if a < x < b
-                                @test HP.theta(I,x) == log(x-a) - log(b-x)
-                            else
-                                @test_throws DomainError HP.theta(I,x)
-                            end
-                        end
-                    end
-                end
-            end
-        end
-    end
-end
-
-
-@testset "Testing HP.eta" begin
-    for T in FloatingPointTypes
-        a = -one(T)
-        b = one(T)
-        for A in (NullBound(T), ClosedBound(a), OpenBound(a))
-            for B in (NullBound(T), ClosedBound(b), OpenBound(b))
-                I = HP.interval(A,B)
-                l_theta = max(HP.lowerboundtheta(I),convert(T,-10))
-                u_theta = min(HP.upperboundtheta(I),convert(T,10))
-
-                for x in range(convert(T, -10), stop = convert(T, 10), length = 201)
-                    if l_theta <= x <= u_theta
-                        v = HP.eta(I,x)
-                        @test isapprox(HP.theta(I,v), x)
-                    else
-                        @test_throws DomainError HP.eta(I,x)
-                    end
-                end
-            end
-        end
-    end
-end
-
-
-@testset "Testing MLK.HyperParameter" begin
-    for T in (FloatingPointTypes..., IntegerTypes...)
-        I = interval(nothing, ClosedBound(zero(T)))
-
-        @test_throws ErrorException HyperParameter(one(T), I)
-
-        P = HyperParameter(one(T), interval(nothing, ClosedBound(one(T))))
-        @test getindex(P.value) == one(T)
-        @test eltype(P) == T
-
-        @test getvalue(P) == one(T)
-
-        @test checkvalue(P, convert(T,2)) == false
-        @test checkvalue(P, zero(T)) == true
-
-        setvalue!(P, zero(T))
-        @test getindex(P.value) == zero(T)
-        @test getvalue(P) == zero(T)
-
-        if T in FloatingPointTypes
-            for a in (NullBound(T), ClosedBound(-one(T)), OpenBound(-one(T)))
-                for b in (NullBound(T), ClosedBound(one(T)), OpenBound(one(T)))
-                    I = HP.Interval(a, b)
-                    P = HyperParameter(zero(T),I)
-
-                    @test isapprox(HP.gettheta(P), HP.theta(I, zero(T)))
-
-                    HP.settheta!(P, HP.theta(I, convert(T,0.5)))
-                    @test isapprox(HP.gettheta(P), HP.theta(I, convert(T,0.5)))
-                    @test isapprox(HP.getvalue(P), convert(T,0.5))
-                end
-            end
-        end
-
-        P1 = HyperParameter(zero(T), interval(T))
-        P2 = HyperParameter(one(T),  interval(T))
-        for op in (isless, ==, +, -, *, /)
-            @test op(P1, one(T)) == op(getvalue(P1), one(T))
-            @test op(one(T), P1) == op(one(T), getvalue(P1))
-            @test op(P1, P2)     == op(getvalue(P1), getvalue(P2))
-        end
-
-        @test show(devnull, P) == nothing
-    end
-end
diff --git a/test/basefunctions.jl b/test/basefunctions.jl
index eb0edfa..a553d9b 100644
--- a/test/basefunctions.jl
+++ b/test/basefunctions.jl
@@ -1,34 +1,54 @@
-for f in base_functions
-    @testset "Testing $f" begin
-        F = (f)()
+function test_base_function(f)
+    F = (f)()
+
+    @testset "Testing properties" begin
+        test_properties = get(base_functions_properties, f, (false,false))
+        @test MLK.isstationary(F) == test_properties[1]
+        @test MLK.isisotropic(F)  == test_properties[2]
+    end
+
+    # Test that initial value is set correctly compared to reference
+    @testset "Testing $(MLK.base_initiate)" begin
         for T in FloatingPointTypes
             s = MLK.base_initiate(F, T)
             s_tmp = convert(T, get(base_functions_initiate, f, 0))
 
             @test typeof(s) == T
             @test s == s_tmp
+        end
+    end
 
+    # Test that aggregation works correctly
+    @testset "Testing $(MLK.base_aggregate)" begin
+        for T in FloatingPointTypes
             f_tmp = get(base_functions_aggregate, f, (s,x,y) -> convert(T, NaN))
-            for x in (zero(T), one(T), convert(T,2)), y in (zero(T), one(T), convert(T,2))
-                s = MLK.base_aggregate(F, one(T), x, y)
-                s_tmp = f_tmp(one(T), x, y)
+            for x in (T(0), T(1), T(2)), y in (T(0), T(1), T(2))
+                s = MLK.base_aggregate(F, T(1), x, y)
+                s_tmp = f_tmp(T(1), x, y)
 
                 @test typeof(s) == T
                 @test s == s_tmp
             end
+        end
+    end
 
+    # Test that return works correctly
+    @testset "Testing $(MLK.base_return)" begin
+        for T in FloatingPointTypes
             f_tmp = get(base_functions_return, f, s -> s)
-            for x in (zero(T), one(T), convert(T,2))
+            for x in (T(0), T(1), T(2))
                 s = MLK.base_return(F, x)
                 s_tmp = f_tmp(x)
 
                 @test typeof(s) == T
                 @test s == s_tmp
             end
-
-            test_properties = get(base_functions_properties, f, (false,false))
-            @test MLK.isstationary(F) == test_properties[1]
-            @test MLK.isisotropic(F)  == test_properties[2]
         end
     end
 end
+
+for f in base_functions
+    @testset "Testing $f" begin
+        test_base_function(f)
+    end
+end
\ No newline at end of file
diff --git a/test/basematrix.jl b/test/basematrix.jl
index af1e619..289e5ff 100644
--- a/test/basematrix.jl
+++ b/test/basematrix.jl
@@ -117,12 +117,21 @@ end
         P_tst_nn = Array{T}(undef, n, n)
         P_tst_nm = Array{T}(undef, n, m)
 
-        for layout in (Val(:row), Val(:col))
-            X = layout == Val(:row) ? permutedims(hcat(X_set...)) : hcat(X_set...)
-            Y = layout == Val(:row) ? permutedims(hcat(Y_set...)) : hcat(Y_set...)
-
-            for f in base_functions
-                F = (f)()
+        for f in base_functions
+            F = (f)()
+
+            #X = permutedims(hcat(X_set...))
+            #Y = permutedims(hcat(Y_set...))
+            #
+            #P = [MLK.base_evaluate(F,x,y) for x in X_set, y in X_set]
+            #@test isapprox(P, MLK.basematrix(F, X, true))
+            #
+            #P = [MLK.base_evaluate(F,x,y) for x in X_set, y in Y_set]
+            #@test isapprox(P, MLK.basematrix(F, X, Y))
+
+            for layout in (Val(:row), Val(:col))
+                X = layout == Val(:row) ? permutedims(hcat(X_set...)) : hcat(X_set...)
+                Y = layout == Val(:row) ? permutedims(hcat(Y_set...)) : hcat(Y_set...)
 
                 P = [MLK.base_evaluate(F,x,y) for x in X_set, y in X_set]
                 @test isapprox(P, MLK.basematrix!(layout, P_tst_nn, F, X, true))
diff --git a/test/kernelfunctions.jl b/test/kernelfunctions.jl
index daca1ae..8a95b2a 100644
--- a/test/kernelfunctions.jl
+++ b/test/kernelfunctions.jl
@@ -1,92 +1,67 @@
-#= Test Constructors =#
+function test_kernel_function(k)
+    default_args, alt_args = get(kernel_functions_arguments, k, ((), ()))
 
-for k in kernel_functions
-    @testset "Testing $k" begin
-        def_args, alt_args = get(kernel_functions_arguments, k, ((), ()))
+    local n = length(default_args)
+    K = (k)()
 
-        # Test Constructors
-        for T in FloatingPointTypes
-            for i = 1:length(alt_args)
-                K = (k)([typeof(θ) <: AbstractFloat ? convert(T,θ) : θ for θ in alt_args[1:i]]...)
-                for j = 1:length(alt_args)
-                    @test MLK.getvalue(getfield(K,j)) == (j <= i ? alt_args[j] : def_args[j])
-                end
-            end
+    @testset "Testing constructors" begin
+        for j = 1:n
+            @test getfield(K,j) == default_args[j]
         end
+        @test eltype(K) == Float64
 
-        # Test Conversions
-        K = (k)()
-        for psi in (Kernel, MercerKernel, NegativeDefiniteKernel)
-            if k <: psi
-                for T1 in FloatingPointTypes
-                    T2 = T1 == Float64 ? Float32 : T1
-                    @test eltype(convert(psi{T2},K)) == T2
-                end
+        for T in FloatingPointTypes, i = 1:n
+            K = (k)([T(a) for a in alt_args[1:i]]...)
+            for j = 1:n
+                @test getfield(K,j) == (j <= i ? alt_args[j] : default_args[j])
             end
+            @test eltype(K) == T
         end
+    end
 
-        # Test Properties
-        @test MLK.ismercer(K) == (typeof(K) <: MLK.MercerKernel ? true : false)
-        @test MLK.isnegdef(K) == (typeof(K) <: MLK.NegativeDefiniteKernel ? true : false)
+    @testset "Testing properties" begin
+        @test MLK.ismercer(K) == isa(K, MLK.MercerKernel)
+        @test MLK.isnegdef(K) == isa(K, MLK.NegativeDefiniteKernel)
         @test MLK.isstationary(K) == MLK.isstationary(MLK.basefunction(K))
         @test MLK.isisotropic(K) == MLK.isisotropic(MLK.basefunction(K))
-
-        # Test Display
-        @test eval(Meta.parse(string(K))) == K
-        @test show(devnull, K) == nothing
     end
-end
 
-@testset "Testing MLK.kappa" begin
-    for k in kernel_functions
-        k_tmp = get(kernel_functions_kappa, k, x->error(""))
-        for T in FloatingPointTypes
-            K = convert(Kernel{T}, (k)())
-            args = T[MLK.getvalue(getfield(K,theta)) for theta in fieldnames(typeof(K))]
-
-            for z in (zero(T), one(T), convert(T,2))
-                v = MLK.kappa(K, z)
-                v_tmp = (k_tmp)(z, args...)
-                @test isapprox(v, v_tmp)
+    @testset "Testing conversions" begin
+        psi = k
+        while psi != Any
+            for T1 in FloatingPointTypes, T2 in FloatingPointTypes
+                K1 = k{T1}()
+                K2 = convert(psi{T2}, K1)
+                @test eltype(K2) == T2
+                @test isa(K2, k)
             end
+            psi = supertype(psi)
         end
     end
-end
 
-@testset "Testing MLK.gettheta" begin
-    for k in kernel_functions
-        K = (k)()
-        K_theta = MLK.gettheta(K)
-        @test K_theta == [MLK.gettheta(getfield(K,field)) for field in MLK.thetafieldnames(K)]
-        @test MLK.checktheta(K, K_theta) == true
+    @testset "Testing kappa function" begin
+        f = get(kernel_functions_kappa, k, x->error(""))
+        args1, args2 = get(kernel_functions_arguments, k, ((), ()))
+        for i = 0:length(args1)
+            args = [j <= i ? args1[j] : args2[j] for j in eachindex(args1)]
+            for T in FloatingPointTypes
+                K = k{T}(args...)
+                for z in [T(0), T(1), T(2)]
+                    v1 = MLK.kappa(K, z)
+                    v2 = (f)(z, args...)
+                    @test isapprox(v1, v2)
+                end
+            end
+        end
     end
-end
 
-@testset "Testing MLK.settheta!" begin
-    for T in FloatingPointTypes
-        alpha1 = one(T)
-        alpha2 = convert(T,0.6)
-        gamma1 = one(T)
-        gamma2 = convert(T,0.4)
-        K = MLK.GammaExponentialKernel(alpha1,gamma1)
-        MLK.settheta!(K, [log(alpha2); log(gamma2)])
-        @test getvalue(K.alpha) == alpha2
-        @test getvalue(K.gamma) == gamma2
-
-        @test_throws DimensionMismatch MLK.settheta!(K, [one(T)])
-        @test_throws DimensionMismatch MLK.settheta!(K, [one(T); one(T); one(T)])
+    @testset "Testing show function" begin
+        @test show(devnull, K) == nothing
     end
 end
 
-
-@testset "Testing MLK.checktheta" begin
-    for T in FloatingPointTypes
-        K = MLK.GammaExponentialKernel(one(T),one(T))
-
-        @test MLK.checktheta(K, [log(one(T)); log(one(T))]) == true
-        @test MLK.checktheta(K, [log(one(T)); log(convert(T,2))]) == false
-
-        @test_throws DimensionMismatch MLK.checktheta(K, [one(T)])
-        @test_throws DimensionMismatch MLK.checktheta(K, [one(T); one(T); one(T)])
+for k in kernel_functions
+    @testset "Testing $k" begin
+        test_kernel_function(k)
     end
-end
+end
\ No newline at end of file
diff --git a/test/reference.jl b/test/reference.jl
index 1731242..67524b6 100644
--- a/test/reference.jl
+++ b/test/reference.jl
@@ -1,100 +1,98 @@
-# Pairwise Function References
+# Base Function Test
 
-base_functions = (
+struct Euclidean <: Metric end
+
+MLK.base_aggregate(::Euclidean, s::T, x::T, y::T) where {T} = s + (x-y)^2
+MLK.base_return(::Euclidean, s::T) where {T} = sqrt(s)
+
+
+# Base Function References
+
+const base_functions = (
     SquaredEuclidean,
-    SineSquared,
-    ChiSquared,
-    ScalarProduct
+    ScalarProduct,
+    Euclidean
 )
-
-base_functions_initiate = Dict(
+    
+const base_functions_initiate = Dict(
     SquaredEuclidean => 0,
-    SineSquared      => 0,
-    ChiSquared       => 0,
-    ScalarProduct    => 0
+    ScalarProduct    => 0,
+    Euclidean        => 0
 )
 
-base_functions_aggregate = Dict(
+const base_functions_aggregate = Dict(
     SquaredEuclidean => (s,x,y) -> s + (x-y)^2,
-    SineSquared      => (s,x,y) -> s + sin((x-y))^2,
-    ChiSquared       => (s,x,y) -> x == y == 0 ? s : s + ((x-y)^2/(x+y)),
-    ScalarProduct    => (s,x,y) -> s + x*y
+    ScalarProduct    => (s,x,y) -> s + x*y,
+    Euclidean        => (s,x,y) -> s + (x-y)^2
 )
 
-base_functions_return = Dict{DataType,Any}()
+const base_functions_return = Dict(
+    SquaredEuclidean => s -> s,
+    ScalarProduct    => s -> s,
+    Euclidean        => s -> sqrt(s)
+)
+
+const base_functions_properties = Dict(
+                       #|stnry |isotrop
+    ScalarProduct    => (false, false),
+    SquaredEuclidean => (true,  true),
+    Euclidean        => (false, false) # just default values test
+)
 
 
 # Kernel Function References
 
-kernel_functions = (
+const kernel_functions = (
     ExponentialKernel,
     SquaredExponentialKernel,
     GammaExponentialKernel,
     RationalQuadraticKernel,
-    GammaRationalKernel,
+    GammaRationalQuadraticKernel,
     MaternKernel,
-    LinearKernel,
     PolynomialKernel,
     ExponentiatedKernel,
-    PeriodicKernel,
     PowerKernel,
     LogKernel,
     SigmoidKernel
 )
 
-kernel_functions_arguments = Dict(
-    ExponentialKernel        => ((1.0,),        (2.0,)),
-    SquaredExponentialKernel => ((1.0,),        (2.0,)),
-    GammaExponentialKernel   => ((1.0,1.0),     (2.0,0.5)),
-    RationalQuadraticKernel  => ((1.0,1.0),     (2.0,2.0)),
-    GammaRationalKernel      => ((1.0,1.0,1.0), (2.0,2.0,0.5)),
-    MaternKernel             => ((1.0,1.0),     (2.0,2.0)),
-    LinearKernel             => ((1.0,1.0),     (2.0,2.0)),
-    PolynomialKernel         => ((1.0,1.0,3),   (2.0,2.0,2)),
-    ExponentiatedKernel      => ((1.0,),        (2.0,)),
-    PeriodicKernel           => ((1.0,),        (2.0,)),
-    PowerKernel              => ((1.0,),        (0.5,)),
-    LogKernel                => ((1.0,1.0),     (2.0,0.5)),
-    SigmoidKernel            => ((1.0,1.0),     (2.0,2.0))
+const kernel_functions_arguments = Dict(
+    ExponentialKernel            => ((1.0,),        (2.0,)),
+    SquaredExponentialKernel     => ((1.0,),        (2.0,)),
+    GammaExponentialKernel       => ((1.0,1.0),     (2.0,0.5)),
+    RationalQuadraticKernel      => ((1.0,1.0),     (2.0,2.0)),
+    GammaRationalQuadraticKernel => ((1.0,1.0,1.0), (2.0,2.0,0.5)),
+    MaternKernel                 => ((1.0,1.0),     (2.0,2.0)),
+    PolynomialKernel             => ((1.0,1.0,3),   (2.0,2.0,2)),
+    ExponentiatedKernel          => ((1.0,),        (2.0,)),
+    PowerKernel                  => ((1.0,),        (0.5,)),
+    LogKernel                    => ((1.0,1.0),     (2.0,0.5)),
+    SigmoidKernel                => ((1.0,1.0),     (2.0,2.0))
 )
 
-kernel_functions_kappa = Dict(
-    ExponentialKernel        => (z,α)     -> exp(-α*sqrt(z)),
-    SquaredExponentialKernel => (z,α)     -> exp(-α*z),
-    GammaExponentialKernel   => (z,α,γ)   -> exp(-α*z^γ),
-    RationalQuadraticKernel  => (z,α,β)   -> (1 + α*z)^(-β),
-    GammaRationalKernel      => (z,α,β,γ) -> (1 + α*z^γ)^(-β),
-    MaternKernel             => (z,ν,ρ)   -> begin
-                                                d = √(z)
-                                                T = typeof(z)
-                                                d = d < eps(T) ? eps(T) : d
-                                                tmp1 = √(2*ν)*d/ρ
-                                                tmp2 = 2^(1 - ν)
-                                                tmp2*(tmp1^ν)*besselk(ν, tmp1)/gamma(ν)
-                                             end,
-    LinearKernel             => (z,a,c)   -> (a*z+c),
-    PolynomialKernel         => (z,a,c,d) -> (a*z+c)^d,
-    ExponentiatedKernel      => (z,a)     -> exp(a*z),
-    PeriodicKernel           => (z,α)     -> exp(-α*z),
-    PowerKernel              => (z,γ)     -> z^γ,
-    LogKernel                => (z,α,γ)   -> log(α*z^γ+1),
-    SigmoidKernel            => (z,a,c)   -> tanh(a*z+c)
+const kernel_functions_kappa = Dict(
+    ExponentialKernel            => (z,α)     -> exp(-α*sqrt(z)),
+    SquaredExponentialKernel     => (z,α)     -> exp(-α*z),
+    GammaExponentialKernel       => (z,α,γ)   -> exp(-α*z^γ),
+    RationalQuadraticKernel      => (z,α,β)   -> (1 + α*z)^(-β),
+    GammaRationalQuadraticKernel => (z,α,β,γ) -> (1 + α*z^γ)^(-β),
+    MaternKernel                 => (z,ν,ρ)   -> begin
+                                                   d = √(z)
+                                                   T = typeof(z)
+                                                   d = d < eps(T) ? eps(T) : d
+                                                   tmp1 = √(2*ν)*d/ρ
+                                                   tmp2 = 2^(1 - ν)
+                                                   tmp2*(tmp1^ν)*besselk(ν, tmp1)/gamma(ν)
+                                                end,
+    PolynomialKernel             => (z,a,c,d) -> (a*z+c)^d,
+    ExponentiatedKernel          => (z,a)     -> exp(a*z),
+    PowerKernel                  => (z,γ)     -> z^γ,
+    LogKernel                    => (z,α,γ)   -> log(α*z^γ+1),
+    SigmoidKernel                => (z,a,c)   -> tanh(a*z+c)
 )
 
-
-kernel_functions_base = Dict(
-    LinearKernel             => ScalarProduct,
+const kernel_functions_base = Dict(
     PolynomialKernel         => ScalarProduct,
     ExponentiatedKernel      => ScalarProduct,
-    PeriodicKernel           => SineSquared,
     SigmoidKernel            => ScalarProduct
-)
-
-
-base_functions_properties = Dict(
-                       #|stnry |isotrop
-    SineSquared      => (true,  false),
-    ScalarProduct    => (false, false),
-    SquaredEuclidean => (true,  true),
-    ChiSquared       => (false, false)
- )
+)
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 6c830ce..7062cc4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,8 +7,6 @@ using MLKernels:
         InnerProduct,
             ScalarProduct,
         PreMetric,
-            ChiSquared,
-            SineSquared,
             Metric,
                 SquaredEuclidean
 
@@ -18,19 +16,32 @@ using SpecialFunctions:
 import LinearAlgebra
 import Statistics
 
-FloatingPointTypes = (Float32, Float64)
-IntegerTypes = (Int32, Int64)
-MLK = MLKernels
-HP = MLKernels.HyperParameters
+const FloatingPointTypes = (Float32, Float64)
+const IntegerTypes = (Int32, Int64)
+const MLK = MLKernels
 
 include("reference.jl")
 
-include("HyperParameters/hyperparameter.jl")
+@testset "Testing utility functions" begin
+    include("utils.jl")
+end
 
-include("utils.jl")
-include("basefunctions.jl")
-include("basematrix.jl")
+@testset "Testing BaseFunction types" begin
+    include("basefunctions.jl")
+end
 
-include("kernelfunctions.jl")
-include("kernelmatrix.jl")
-include("nystrom.jl")
\ No newline at end of file
+@testset "Testing base function evaluation" begin
+    include("basematrix.jl")
+end
+
+@testset "Testing Kernel types" begin
+    include("kernelfunctions.jl")
+end
+
+@testset "Testing kernel function evaluation" begin
+    include("kernelmatrix.jl")
+end
+
+@testset "Testing nystrom" begin
+    include("nystrom.jl")
+end
\ No newline at end of file
diff --git a/test/utils.jl b/test/utils.jl
index 3c05348..f093442 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -2,6 +2,16 @@ n = 3
 m = 2
 p = 5
 
+@testset "Testing $(MLK.promote_float)" begin
+    @test MLK.promote_float() == Float64
+    @test MLK.promote_float(Float16) == Float16
+    @test MLK.promote_float(Float16, Float32) == Float32
+    @test MLK.promote_float(Float16, Float32, Float64) == Float64
+    @test MLK.promote_float(BigFloat, Float64, Float32, Float16) == BigFloat
+    @test MLK.promote_float(Int64) == Float64
+    @test MLK.promote_float(Int64, Float32) == Float32
+end
+
 @testset "Testing $(MLK.dotvectors!)" begin
     for T in FloatingPointTypes
         Set_X = [rand(T, p) for i = 1:n]