ARD Implementation and Automatic Differentation compatibility

Manifest.toml

@@ -9,15 +9,15 @@ version = "0.8.10"
 
 [[BinaryProvider]]
 deps = ["Libdl", "Pkg", "SHA", "Test"]
-git-tree-sha1 = "9930c1a6cd49d9fcd7218df6be417e6ae4f1468a"
+git-tree-sha1 = "055eb2690182ebc31087859c3dd8598371d3ef9e"
 uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
-version = "0.5.2"
+version = "0.5.3"
 
 [[Compat]]
 deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
-git-tree-sha1 = "2d9e14d19bad3f9ad5cc5e4cffabc3cfa59de825"
+git-tree-sha1 = "ec61a16eed883ad0cfa002d7489b3ce6d039bb9a"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "1.3.0"
+version = "1.4.0"
 
 [[Dates]]
 deps = ["Printf"]

src/basefunctions/scalarproduct.jl

@@ -1,11 +1,11 @@
 @doc raw"""
     ScalarProduct()
-    
+
 The scalar product is an inner product of the form:
 
 ```math
 f(\mathbf{x}, \mathbf{y}) = \mathbf{x}^{\intercal}\mathbf{y}
 ```
 """
 struct ScalarProduct <: InnerProduct end
-@inline base_aggregate(::ScalarProduct, s::T, x::T, y::T) where {T} = s + x*y
+@inline base_aggregate(::ScalarProduct, s::T, scale::T, x::T, y::T) where {T} = s + scale*x*y

src/basefunctions/squaredeuclidean.jl

@@ -9,7 +9,7 @@ f(\mathbf{x}, \mathbf{y}) = (\mathbf{x} - \mathbf{y})^{\intercal}(\mathbf{x} - \
 """
 struct SquaredEuclidean <: Metric end
 
-@inline base_aggregate(::SquaredEuclidean, s::T, x::T, y::T) where {T} = s + (x-y)^2
+@inline base_aggregate(::SquaredEuclidean, s::T, scale::T, x::T, y::T) where {T} = s + scale*(x-y)^2
 
 @inline isstationary(::SquaredEuclidean) = true
 @inline isisotropic(::SquaredEuclidean)  = true

src/basematrix.jl

@@ -3,36 +3,75 @@
 @inline base_initiate(::BaseFunction, ::Type{T}) where {T} = zero(T)
 @inline base_return(::BaseFunction, s::T) where {T} = s
 
-function base_evaluate(f::BaseFunction, x::T, y::T) where {T<:AbstractFloat}
-    base_return(f, base_aggregate(f, base_initiate(f,T), x, y))
+function base_evaluate(f::BaseFunction, scale::T, x::T, y::T) where {T<:Real}
+    base_return(f, base_aggregate(f, base_initiate(f,T), scale, x, y))
 end
 
+function base_evaluate(f::BaseFunction, scale::AbstractVector{T}, x::T, y::T) where {T<:Real}
+    base_return(f, base_aggregate(f, base_initiate(f,T), scale[1], x, y))
+end
+
+
 # Note: no checks, assumes length(x) == length(y) >= 1
 function unsafe_base_evaluate(
         f::BaseFunction,
+        scale::AbstractArray{T},
         x::AbstractArray{T},
-        y::AbstractArray{T}
-    ) where {T<:AbstractFloat}
+        y::AbstractArray{T},
+    ) where {T<:Real}
+    s = base_initiate(f, T)
+    @simd for I in eachindex(x, y)
+        @inbounds xi = x[I]
+        @inbounds yi = y[I]
+        @inbounds si = scale[I]
+        s = base_aggregate(f, s, si, xi, yi)
+    end
+    base_return(f, s)
+end
+
+function unsafe_base_evaluate(
+        f::BaseFunction,
+        scale::T,
+        x::AbstractArray{T},
+        y::AbstractArray{T},
+    ) where {T<:Real}
     s = base_initiate(f, T)
     @simd for I in eachindex(x, y)
         @inbounds xi = x[I]
         @inbounds yi = y[I]
-        s = base_aggregate(f, s, xi, yi)
+        s = base_aggregate(f, s, scale, xi, yi)
     end
     base_return(f, s)
 end
 
 function base_evaluate(
         f::BaseFunction,
+        scale::AbstractArray{T},
         x::AbstractArray{T},
         y::AbstractArray{T}
-    ) where {T<:AbstractFloat}
+    ) where {T<:Real}
     if (n = length(x)) != length(y)
         throw(DimensionMismatch("Arrays x and y must have the same length."))
+    elseif n != length(scale)
+        throw(DimensionMismatch("Arrays x and y must have the same length as the scaling vector of the kernel"))
     elseif n == 0
         throw(DimensionMismatch("Arrays x and y must be at least of length 1."))
     end
-    unsafe_base_evaluate(f, x, y)
+    unsafe_base_evaluate(f, scale, x, y)
+end
+
+function base_evaluate(
+        f::BaseFunction,
+        scale::T,
+        x::AbstractArray{T},
+        y::AbstractArray{T}
+    ) where {T<:Real}
+    if (n = length(x)) != length(y)
+        throw(DimensionMismatch("Arrays x and y must have the same length."))
+    elseif n == 0
+        throw(DimensionMismatch("Arrays x and y must be at least of length 1."))
+    end
+    unsafe_base_evaluate(f, scale, x, y)
 end
 
 
@@ -57,15 +96,15 @@ for orientation in (:row, :col)
         @inline function allocate_basematrix(
                 ::Val{$(Meta.quot(orientation))},
                 X::AbstractMatrix{T}
-            ) where {T<:AbstractFloat}
+            ) where {T<:Real}
             Array{T}(undef, size(X,$dim_obs), size(X,$dim_obs))
         end
 
         @inline function allocate_basematrix(
                 ::Val{$(Meta.quot(orientation))},
                 X::AbstractMatrix{T},
                 Y::AbstractMatrix{T}
-            ) where {T<:AbstractFloat}
+            ) where {T<:Real}
             Array{T}(undef, size(X,$dim_obs), size(Y,$dim_obs))
         end
 
@@ -112,15 +151,35 @@ function basematrix!(
         σ::Orientation,
         P::Matrix{T},
         f::BaseFunction,
+        scale::AbstractArray{T},
+        X::AbstractMatrix{T},
+        symmetrize::Bool
+    ) where {T<:Real}
+    n = checkdimensions(σ, P, X)
+    for j = 1:n
+        xj = subvector(σ, X, j)
+        for i = 1:j
+            xi = subvector(σ, X, i)
+            @inbounds P[i,j] = unsafe_base_evaluate(f, scale, xi, xj)
+        end
+    end
+    symmetrize ? LinearAlgebra.copytri!(P, 'U', false) : P
+end
+
+function basematrix!(
+        σ::Orientation,
+        P::Matrix{T},
+        f::BaseFunction,
+        scale::T,
         X::AbstractMatrix{T},
         symmetrize::Bool
-    ) where {T<:AbstractFloat}
+    ) where {T<:Real}
     n = checkdimensions(σ, P, X)
     for j = 1:n
         xj = subvector(σ, X, j)
         for i = 1:j
             xi = subvector(σ, X, i)
-            @inbounds P[i,j] = unsafe_base_evaluate(f, xi, xj)
+            @inbounds P[i,j] = unsafe_base_evaluate(f, scale, xi, xj)
         end
     end
     symmetrize ? LinearAlgebra.copytri!(P, 'U', false) : P
@@ -130,20 +189,39 @@ function basematrix!(
         σ::Orientation,
         P::Matrix{T},
         f::BaseFunction,
+        scale::AbstractArray{T},
         X::AbstractMatrix{T},
         Y::AbstractMatrix{T},
-    ) where {T<:AbstractFloat}
+    ) where {T<:Real}
     n, m = checkdimensions(σ, P, X, Y)
     for j = 1:m
         yj = subvector(σ, Y, j)
         for i = 1:n
             xi = subvector(σ, X, i)
-            @inbounds P[i,j] = unsafe_base_evaluate(f, xi, yj)
+            @inbounds P[i,j] = unsafe_base_evaluate(f, scale, xi, yj)
         end
     end
     P
 end
 
+function basematrix!(
+        σ::Orientation,
+        P::Matrix{T},
+        f::BaseFunction,
+        scale::T,
+        X::AbstractMatrix{T},
+        Y::AbstractMatrix{T},
+    ) where {T<:Real}
+    n, m = checkdimensions(σ, P, X, Y)
+    for j = 1:m
+        yj = subvector(σ, Y, j)
+        for i = 1:n
+            xi = subvector(σ, X, i)
+            @inbounds P[i,j] = unsafe_base_evaluate(f, scale, xi, yj)
+        end
+    end
+    P
+end
 
 # ScalarProduct using BLAS/Built-In methods ================================================
 
@@ -174,7 +252,7 @@ function squared_distance!(
         G::Matrix{T},
         xᵀx::Vector{T},
         symmetrize::Bool
-    ) where {T<:AbstractFloat}
+    ) where {T<:Real}
     if !((n = size(G,1)) == size(G,2))
         throw(DimensionMismatch("Gramian matrix must be square."))
     end
@@ -195,7 +273,7 @@ function squared_distance!(
         G::Matrix{T},
         xᵀx::Vector{T},
         yᵀy::Vector{T}
-    ) where {T<:AbstractFloat}
+    ) where {T<:Real}
     n, m = size(G)
     if length(xᵀx) != n
         throw(DimensionMismatch("Length of xᵀx must match rows of G"))
@@ -217,7 +295,7 @@ function fix_negatives!(
         X::Matrix{T},
         symmetrize::Bool,
         ϵ::T=zero(T)
-    ) where {T<:AbstractFloat}
+    ) where {T<:Real}
     if !((n = size(D,1)) == size(D,2))
         throw(DimensionMismatch("Distance matrix must be square."))
     end
@@ -239,7 +317,7 @@ function fix_negatives!(
         X::Matrix{T},
         Y::Matrix{T},
         ϵ::T=zero(T)
-    ) where {T<:AbstractFloat}
+    ) where {T<:Real}
     n, m = size(D)
     for j = 1:m
         yj = subvector(σ, Y, j)
@@ -279,4 +357,4 @@ function basematrix!(
     yᵀy = dotvectors(σ, Y)
     squared_distance!(P, xᵀx, yᵀy)
     fix_negatives!(σ, P, X, Y)
-end
+end

src/deprecated.jl

@@ -9,9 +9,9 @@ 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}
+struct PeriodicKernel{T<:Real} <: MercerKernel{T}
     α::T
-    function PeriodicKernel{T}(α::Real) where {T<:AbstractFloat}
+    function PeriodicKernel{T}(α::Real) where {T<:Real}
         Base.depwarn("PeriodicKernel will be removed in the next major release", :PeriodicKernel)
         @check_args(PeriodicKernel, α, α > zero(α), "α > 0")
         new{T}(α)
@@ -23,4 +23,4 @@ PeriodicKernel(α::T₁ = 1.0) where {T₁<:Real} = PeriodicKernel{promote_float
 
 @inline function kappa(κ::PeriodicKernel{T}, z::T) where {T}
     return exp(-κ.α*z)
-end
+end

src/kernelfunctions.jl

@@ -1,6 +1,6 @@
 # Kernel Functions =========================================================================
 
-abstract type Kernel{T<:AbstractFloat} end
+abstract type Kernel{T<:Real} end
 
 function string(κ::Kernel{T}) where {T}
     args = [string(getfield(κ,θ)) for θ in fieldnames(typeof(κ))]
@@ -46,7 +46,7 @@ isisotropic(κ::Kernel)  = isisotropic(basefunction(κ))
 
 # Mercer Kernels ===========================================================================
 
-abstract type MercerKernel{T<:AbstractFloat} <: Kernel{T} end
+abstract type MercerKernel{T<:Real} <: Kernel{T} end
 
 @inline ismercer(::MercerKernel) = true
 
@@ -65,7 +65,7 @@ end
 
 # Negative Definite Kernels ================================================================
 
-abstract type NegativeDefiniteKernel{T<:AbstractFloat} <: Kernel{T} end
+abstract type NegativeDefiniteKernel{T<:Real} <: Kernel{T} end
 
 @inline isnegdef(::NegativeDefiniteKernel) = true
 
@@ -87,4 +87,4 @@ const other_kernels = [
 
 for kname in other_kernels
     include(joinpath("kernelfunctions", "$(kname).jl"))
-end
+end