diff --git a/Project.toml b/Project.toml
index 9f286dfb2..d69ff6532 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "KernelFunctions"
 uuid = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
-version = "0.10.26"
+version = "0.10.27"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -32,7 +32,7 @@ IrrationalConstants = "0.1"
 LogExpFunctions = "0.2.1, 0.3"
 ParameterHandling = "0.4"
 Requires = "1.0.1"
-SpecialFunctions = "0.8, 0.9, 0.10, 1"
+SpecialFunctions = "0.8, 0.9, 0.10, 1, 2"
 StatsBase = "0.32, 0.33"
 TensorCore = "0.1"
 ZygoteRules = "0.2"
diff --git a/src/basekernels/fbm.jl b/src/basekernels/fbm.jl
index 6db795a1a..11cba1088 100644
--- a/src/basekernels/fbm.jl
+++ b/src/basekernels/fbm.jl
@@ -35,16 +35,16 @@ function (κ::FBMKernel)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
     modX = sum(abs2, x)
     modY = sum(abs2, y)
     modXY = sqeuclidean(x, y)
-    h = first(κ.h)
+    h = only(κ.h)
     return (modX^h + modY^h - modXY^h) / 2
 end
 
 function (κ::FBMKernel)(x::Real, y::Real)
-    return (abs2(x)^first(κ.h) + abs2(y)^first(κ.h) - abs2(x - y)^first(κ.h)) / 2
+    return (abs2(x)^only(κ.h) + abs2(y)^only(κ.h) - abs2(x - y)^only(κ.h)) / 2
 end
 
 function Base.show(io::IO, κ::FBMKernel)
-    return print(io, "Fractional Brownian Motion Kernel (h = ", first(κ.h), ")")
+    return print(io, "Fractional Brownian Motion Kernel (h = ", only(κ.h), ")")
 end
 
 _fbm(modX, modY, modXY, h) = (modX^h + modY^h - modXY^h) / 2
diff --git a/src/basekernels/matern.jl b/src/basekernels/matern.jl
index d792b409f..7f255aeac 100644
--- a/src/basekernels/matern.jl
+++ b/src/basekernels/matern.jl
@@ -38,7 +38,7 @@ function ParameterHandling.flatten(::Type{T}, k::MaternKernel{S}) where {T<:Real
 end
 
 @inline function kappa(κ::MaternKernel, d::Real)
-    result = _matern(first(κ.ν), d)
+    result = _matern(only(κ.ν), d)
     return ifelse(iszero(d), one(result), result)
 end
 
@@ -50,7 +50,7 @@ end
 metric(k::MaternKernel) = k.metric
 
 function Base.show(io::IO, κ::MaternKernel)
-    return print(io, "Matern Kernel (ν = ", first(κ.ν), ", metric = ", κ.metric, ")")
+    return print(io, "Matern Kernel (ν = ", only(κ.ν), ", metric = ", κ.metric, ")")
 end
 
 ## Matern12Kernel = ExponentialKernel aliased in exponential.jl
diff --git a/src/distances/sinus.jl b/src/distances/sinus.jl
index 4bcf4bdf0..51d14c47d 100644
--- a/src/distances/sinus.jl
+++ b/src/distances/sinus.jl
@@ -2,10 +2,12 @@ struct Sinus{T} <: Distances.UnionSemiMetric
     r::Vector{T}
 end
 
+Sinus(r::Real) = Sinus([r])
+
 Distances.parameters(d::Sinus) = d.r
 @inline Distances.eval_op(::Sinus, a::Real, b::Real, p::Real) = abs2(sinpi(a - b) / p)
 @inline (dist::Sinus)(a::AbstractArray, b::AbstractArray) = Distances._evaluate(dist, a, b)
-@inline (dist::Sinus)(a::Number, b::Number) = abs2(sinpi(a - b) / first(dist.r))
+@inline (dist::Sinus)(a::Number, b::Number) = abs2(sinpi(a - b) / only(dist.r))
 
 Distances.result_type(::Sinus{T}, Ta::Type, Tb::Type) where {T} = promote_type(T, Ta, Tb)
 
diff --git a/src/kernels/transformedkernel.jl b/src/kernels/transformedkernel.jl
index 57c2a1d0c..0de2b4f2b 100644
--- a/src/kernels/transformedkernel.jl
+++ b/src/kernels/transformedkernel.jl
@@ -43,10 +43,10 @@ function (k::TransformedKernel{<:SimpleKernel,<:ScaleTransform})(
 end
 
 function _scale(t::ScaleTransform, metric::Euclidean, x, y)
-    return first(t.s) * evaluate(metric, x, y)
+    return only(t.s) * evaluate(metric, x, y)
 end
 function _scale(t::ScaleTransform, metric::Union{SqEuclidean,DotProduct}, x, y)
-    return first(t.s)^2 * evaluate(metric, x, y)
+    return only(t.s)^2 * evaluate(metric, x, y)
 end
 _scale(t::ScaleTransform, metric, x, y) = evaluate(metric, t(x), t(y))
 
diff --git a/src/transform/ardtransform.jl b/src/transform/ardtransform.jl
index 72bc8f9f4..7abb5928a 100644
--- a/src/transform/ardtransform.jl
+++ b/src/transform/ardtransform.jl
@@ -35,7 +35,7 @@ end
 
 dim(t::ARDTransform) = length(t.v)
 
-(t::ARDTransform)(x::Real) = first(t.v) * x
+(t::ARDTransform)(x::Real) = only(t.v) * x
 (t::ARDTransform)(x) = t.v .* x
 
 _map(t::ARDTransform, x::AbstractVector{<:Real}) = t.v' .* x
diff --git a/src/utils.jl b/src/utils.jl
index b7fd48a09..cfcd72b24 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -178,7 +178,7 @@ function validate_inputs(x, y)
     if dim(x) != dim(y) # Passes by default if `dim` is not defined
         throw(
             DimensionMismatch(
-                "Dimensionality of x ($(dim(x))) not equality to that of y ($(dim(y)))"
+                "dimensionality of x ($(dim(x))) is not equal to that of y ($(dim(y)))"
             ),
         )
     end
diff --git a/test/Project.toml b/test/Project.toml
index ef3a56dc2..de1dbd0d6 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9"
+Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
@@ -17,6 +18,7 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 AxisArrays = "0.4.3"
+Compat = "3"
 Distances = "= 0.10.0, = 0.10.1, = 0.10.2, = 0.10.3, = 0.10.4"
 Documenter = "0.25, 0.26, 0.27"
 FiniteDifferences = "0.10.8, 0.11, 0.12"
diff --git a/test/basekernels/constant.jl b/test/basekernels/constant.jl
index fb9f9d0d7..1e7811512 100644
--- a/test/basekernels/constant.jl
+++ b/test/basekernels/constant.jl
@@ -38,6 +38,6 @@
         # Standardised tests.
         TestUtils.test_interface(k, Float64)
         test_params(k, ([log(c)],))
-        test_ADs(c -> ConstantKernel(; c=first(c)), [c])
+        test_ADs(c -> ConstantKernel(; c=only(c)), [c])
     end
 end
diff --git a/test/basekernels/exponential.jl b/test/basekernels/exponential.jl
index fcb3339b2..ecf42c114 100644
--- a/test/basekernels/exponential.jl
+++ b/test/basekernels/exponential.jl
@@ -58,7 +58,7 @@
         @test metric(k2) isa WeightedEuclidean
         @test k2(v1, v2) ≈ k(v1, v2)
 
-        test_ADs(γ -> GammaExponentialKernel(; gamma=first(γ)), [1 + 0.5 * rand()])
+        test_ADs(γ -> GammaExponentialKernel(; gamma=only(γ)), [1 + 0.5 * rand()])
         test_params(k, ([logit(γ / 2)],))
         TestUtils.test_interface(GammaExponentialKernel(; γ=1.36))
 
diff --git a/test/distances/sinus.jl b/test/distances/sinus.jl
index 91ba1b028..d903e765d 100644
--- a/test/distances/sinus.jl
+++ b/test/distances/sinus.jl
@@ -5,5 +5,6 @@
     d = KernelFunctions.Sinus(p)
     @test Distances.parameters(d) == p
     @test evaluate(d, A, B) == sum(abs2.(sinpi.(A - B) ./ p))
-    @test d(3.0, 2.0) == abs2(sinpi(3.0 - 2.0) / first(p))
+    d1 = KernelFunctions.Sinus(first(p))
+    @test d1(3.0, 2.0) == abs2(sinpi(3.0 - 2.0) / first(p))
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 63c8e8f89..c9bc932f3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -14,6 +14,7 @@ using Zygote: Zygote
 using ForwardDiff: ForwardDiff
 using ReverseDiff: ReverseDiff
 using FiniteDifferences: FiniteDifferences
+using Compat: only
 
 using KernelFunctions: SimpleKernel, metric, kappa, ColVecs, RowVecs, TestUtils
 
diff --git a/test/test_utils.jl b/test/test_utils.jl
index bd3c40b50..dce16a969 100644
--- a/test/test_utils.jl
+++ b/test/test_utils.jl
@@ -45,7 +45,7 @@ const FDM = FiniteDifferences.central_fdm(5, 1)
 gradient(f, s::Symbol, args) = gradient(f, Val(s), args)
 
 function gradient(f, ::Val{:Zygote}, args)
-    g = first(Zygote.gradient(f, args))
+    g = only(Zygote.gradient(f, args))
     if isnothing(g)
         if args isa AbstractArray{<:Real}
             return zeros(size(args)) # To respect the same output as other ADs
@@ -66,7 +66,11 @@ function gradient(f, ::Val{:ReverseDiff}, args)
 end
 
 function gradient(f, ::Val{:FiniteDiff}, args)
-    return first(FiniteDifferences.grad(FDM, f, args))
+    return only(FiniteDifferences.grad(FDM, f, args))
+end
+
+function compare_gradient(f, ::Val{:FiniteDiff}, args)
+    @test_nowarn gradient(f, :FiniteDiff, args)
 end
 
 function compare_gradient(f, AD::Symbol, args)
@@ -88,7 +92,7 @@ testdiagfunction(k::MOKernel, A, B) = sum(kernelmatrix_diag(k, A, B))
 function test_ADs(
     kernelfunction, args=nothing; ADs=[:Zygote, :ForwardDiff, :ReverseDiff], dims=[3, 3]
 )
-    test_fd = test_FiniteDiff(kernelfunction, args, dims)
+    test_fd = test_AD(:FiniteDiff, kernelfunction, args, dims)
     if !test_fd.anynonpass
         for AD in ADs
             test_AD(AD, kernelfunction, args, dims)
@@ -114,70 +118,6 @@ function test_ADs(
     end
 end
 
-function test_FiniteDiff(kernelfunction, args=nothing, dims=[3, 3])
-    # Init arguments :
-    k = if args === nothing
-        kernelfunction()
-    else
-        kernelfunction(args)
-    end
-    rng = MersenneTwister(42)
-    @testset "FiniteDifferences" begin
-        if k isa SimpleKernel
-            for d in log.([eps(), rand(rng)])
-                @test_nowarn gradient(:FiniteDiff, [d]) do x
-                    kappa(k, exp(first(x)))
-                end
-            end
-        end
-        ## Testing Kernel Functions
-        x = rand(rng, dims[1])
-        y = rand(rng, dims[1])
-        @test_nowarn gradient(:FiniteDiff, x) do x
-            k(x, y)
-        end
-        if !(args === nothing)
-            @test_nowarn gradient(:FiniteDiff, args) do p
-                kernelfunction(p)(x, y)
-            end
-        end
-        ## Testing Kernel Matrices
-        A = rand(rng, dims...)
-        B = rand(rng, dims...)
-        for dim in 1:2
-            @test_nowarn gradient(:FiniteDiff, A) do a
-                testfunction(k, a, dim)
-            end
-            @test_nowarn gradient(:FiniteDiff, A) do a
-                testfunction(k, a, B, dim)
-            end
-            @test_nowarn gradient(:FiniteDiff, B) do b
-                testfunction(k, A, b, dim)
-            end
-            if !(args === nothing)
-                @test_nowarn gradient(:FiniteDiff, args) do p
-                    testfunction(kernelfunction(p), A, B, dim)
-                end
-            end
-
-            @test_nowarn gradient(:FiniteDiff, A) do a
-                testdiagfunction(k, a, dim)
-            end
-            @test_nowarn gradient(:FiniteDiff, A) do a
-                testdiagfunction(k, a, B, dim)
-            end
-            @test_nowarn gradient(:FiniteDiff, B) do b
-                testdiagfunction(k, A, b, dim)
-            end
-            if args !== nothing
-                @test_nowarn gradient(:FiniteDiff, args) do p
-                    testdiagfunction(kernelfunction(p), A, B, dim)
-                end
-            end
-        end
-    end
-end
-
 function test_FiniteDiff(k::MOKernel, dims=(in=3, out=2, obs=3))
     rng = MersenneTwister(42)
     @testset "FiniteDifferences" begin
@@ -224,68 +164,68 @@ end
 
 function test_AD(AD::Symbol, kernelfunction, args=nothing, dims=[3, 3])
     @testset "$(AD)" begin
-        # Test kappa function
         k = if args === nothing
             kernelfunction()
         else
             kernelfunction(args)
         end
         rng = MersenneTwister(42)
+
         if k isa SimpleKernel
-            for d in log.([eps(), rand(rng)])
-                compare_gradient(AD, [d]) do x
-                    kappa(k, exp(x[1]))
+            @testset "kappa function" begin
+                for d in log.([eps(), rand(rng)])
+                    compare_gradient(AD, [d]) do x
+                        kappa(k, exp(x[1]))
+                    end
                 end
             end
         end
-        # Testing kernel evaluations
-        x = rand(rng, dims[1])
-        y = rand(rng, dims[1])
-        compare_gradient(AD, x) do x
-            k(x, y)
-        end
-        compare_gradient(AD, y) do y
-            k(x, y)
-        end
-        if !(args === nothing)
-            compare_gradient(AD, args) do p
-                kernelfunction(p)(x, y)
-            end
-        end
-        # Testing kernel matrices
-        A = rand(rng, dims...)
-        B = rand(rng, dims...)
-        for dim in 1:2
-            compare_gradient(AD, A) do a
-                testfunction(k, a, dim)
-            end
-            compare_gradient(AD, A) do a
-                testfunction(k, a, B, dim)
+
+        @testset "kernel evaluations" begin
+            x = rand(rng, dims[1])
+            y = rand(rng, dims[1])
+            compare_gradient(AD, x) do x
+                k(x, y)
             end
-            compare_gradient(AD, B) do b
-                testfunction(k, A, b, dim)
+            compare_gradient(AD, y) do y
+                k(x, y)
             end
             if !(args === nothing)
-                compare_gradient(AD, args) do p
-                    testfunction(kernelfunction(p), A, dim)
+                @testset "hyperparameters" begin
+                    compare_gradient(AD, args) do p
+                        kernelfunction(p)(x, y)
+                    end
                 end
             end
+        end
 
-            compare_gradient(AD, A) do a
-                testdiagfunction(k, a, dim)
-            end
-            compare_gradient(AD, A) do a
-                testdiagfunction(k, a, B, dim)
-            end
-            compare_gradient(AD, B) do b
-                testdiagfunction(k, A, b, dim)
-            end
-            if args !== nothing
-                compare_gradient(AD, args) do p
-                    testdiagfunction(kernelfunction(p), A, dim)
+        @testset "kernel matrices" begin
+            A = rand(rng, dims...)
+            B = rand(rng, dims...)
+            @testset "$(_testfn)" for _testfn in (testfunction, testdiagfunction)
+                for dim in 1:2
+                    compare_gradient(AD, A) do a
+                        _testfn(k, a, dim)
+                    end
+                    compare_gradient(AD, A) do a
+                        _testfn(k, a, B, dim)
+                    end
+                    compare_gradient(AD, B) do b
+                        _testfn(k, A, b, dim)
+                    end
+                    if !(args === nothing)
+                        @testset "hyperparameters" begin
+                            compare_gradient(AD, args) do p
+                                _testfn(kernelfunction(p), A, dim)
+                            end
+                            compare_gradient(AD, args) do p
+                                _testfn(kernelfunction(p), A, B, dim)
+                            end
+                        end
+                    end
                 end
             end
-        end
+        end # kernel matrices
     end
 end
 
diff --git a/test/transform/chaintransform.jl b/test/transform/chaintransform.jl
index 6e3d8a44f..f8b19cffe 100644
--- a/test/transform/chaintransform.jl
+++ b/test/transform/chaintransform.jl
@@ -24,6 +24,8 @@
     @test repr(tp ∘ tf) == "Chain of 2 transforms:\n\t - $(tf) |> $(tp)"
     test_ADs(
         x -> SEKernel() ∘ (ScaleTransform(exp(x[1])) ∘ ARDTransform(exp.(x[2:4]))),
-        randn(rng, 4),
+        randn(rng, 4);
+        ADs=[:ForwardDiff, :ReverseDiff],  # explicitly pass ADs to exclude :Zygote
     )
+    @test_broken "test_AD of chain transform is currently broken in Zygote, see GitHub issue #263"
 end