Add RBF estimation

Project.toml

@@ -1,7 +1,7 @@
 name = "KissSmoothing"
 uuid = "23b0397c-cd08-4270-956a-157331f0528f"
 authors = ["Francesco Alemanno <francescoalemanno710@gmail.com>"]
-version = "1.0.0"
+version = "1.0.1"
 
 [deps]
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

README.md

@@ -1,6 +1,6 @@
 # KissSmoothing.jl
-
-This package implements a smoothing procedure
+This package implements a denoising procedure, and a Radial Basis Function estimation procedure
+## Denoising
 
     denoise(V::Array; factor=1.0, rtol=1e-12, dims=ndims(V), verbose = false)
 
@@ -25,7 +25,7 @@ returns a tuple (S,N) where:
 in particular `S + N` reconstructs the original data `V`.
 
 
-# Example
+### Example
 
 ```julia
 using KissSmoothing, Statistics, LinearAlgebra
@@ -54,7 +54,7 @@ savefig("test.png")
 ```
 ![test.png](test.png "Plot of 1D signal smoothing")
 
-## Multidimensional example
+### Multidimensional example
 ```julia
 using KissSmoothing, Statistics, LinearAlgebra
 using PyPlot
@@ -81,3 +81,36 @@ tight_layout()
 savefig("test_multi.png")
 ```
 ![test_multi.png](test_multi.png "Plot of multidim smoothing")
+
+## RBF Estimation
+
+    fit_rbf(xv::Array, yv::Array, cp::Array)
+
+fit thin-plate radial basis function according to:
+
+    `xv` : array NxP, N number of training points, P number of input variables
+
+    `yv` : array NxQ, N number of training points, Q number of output variables
+
+    `cp` : array KxP, K number of control points, P number of input variables
+
+returns a callable RBF object.
+
+### Example
+
+```julia
+using PyPlot
+t = LinRange(0,2pi,1000)
+ty = sin.(t)
+y = ty .+ randn(length(t)) .*0.05
+fn = fit_rbf(t,y,LinRange(0,2pi,20))
+scatter(t, y, color="gray",s=2,label="noisy")
+plot(t, fn(t), color="red",lw=1.5,label="rbf estimate")
+plot(t,ty, color="blue",lw=1.0,label="true")
+xlabel("X")
+ylabel("Y")
+legend()
+tight_layout()
+savefig("rbf.png")
+```
+![rbf.png](rbf.png "Plot of rbf estimation")

src/KissSmoothing.jl

@@ -31,7 +31,7 @@ returns a tuple (S,N) where:
 in particular `S + N` reconstructs the original data `V`.
 """
 function denoise(
-    V::Array{Float64,N};
+    V::AbstractArray{Float64,N};
     factor::Float64 = 1.0,
     rtol::Float64 = 1e-12,
     dims::Int64 = N,
@@ -48,7 +48,7 @@ function denoise(
     end
     iV = dct(V)
     stri = map(i -> ifelse(i == dims, lV, 1), 1:ndims(V))
-    X = map(abs2,reshape(LinRange(0.0,1.0,lV), stri...))
+    X = map(abs2, reshape(LinRange(0.0, 1.0, lV), stri...))
     d = factor * mean(abs, diff(V, dims = dims)) * (K1 / K2)
     σt = 0.5
     σd = 0.25
@@ -66,7 +66,7 @@ function denoise(
         if verbose
             println(buf, iter, "  ", σ, "  ", Δ)
         end
-        if abs(d-c) < rtol*d
+        if abs(d - c) < rtol * d
             break
         end
     end
@@ -76,5 +76,59 @@ function denoise(
     f, V .- f
 end
 
-export denoise
+function tps(r)
+    if iszero(r)
+        zero(r)
+    else
+        r * r * log(r)
+    end
+end
+
+function dist(x, y)
+    mapreduce(+, x, y) do a, b
+        abs2(a - b)
+    end
+end
+
+struct RBF{G<:AbstractArray{Float64},C<:AbstractArray{Float64}}
+    Γ::G
+    C::C
+end
+
+function evalPhi(xs::AbstractArray{Float64}, cp::AbstractArray{Float64})
+    Phi = zeros(size(xs, 1), size(cp, 1) + 1)
+    for i = 1:size(xs, 1), j = 1:size(cp, 1)
+        Phi[i, j] = tps(dist(xs[i, :], cp[j, :]))
+    end
+    Phi[:, end] .= 1
+    Phi
+end
+
+function (net::RBF)(X::AbstractArray{Float64})
+    evalPhi(X, net.C) * net.Γ
+end
+
+
+"""
+    fit_rbf(xv::Array, yv::Array, cp::Array)
+
+fit thin-plate radial basis function according to:
+
+    `xv` : array NxP, N number of training points, P number of input variables
+
+    `yv` : array NxQ, N number of training points, Q number of output variables
+
+    `cp` : array KxP, K number of control points, P number of input variables
+
+returns a callable RBF object.
+"""
+function fit_rbf(
+    xv::AbstractArray{Float64},
+    yv::AbstractArray{Float64},
+    cp::AbstractArray{Float64},
+)
+    RBF(evalPhi(xv, cp) \ yv, collect(cp))
+end
+
+export denoise, fit_rbf, RBF
 end # module

test/runtests.jl

@@ -14,32 +14,41 @@ using Random
             c1 = sum(abs2, yr .- y)
             c2 = sum(abs2, ys .- y)
             @test c2 < c1
-            s1 = sum(abs2,n)
-            s2 = sum(abs2,yn)
-            @test abs(log2(s2/s1)) < 1
+            s1 = sum(abs2, n)
+            s2 = sum(abs2, yn)
+            @test abs(log2(s2 / s1)) < 1
         end
     end
 end
 
 @testset "Too few points" begin
     rng = Random.MersenneTwister(1337)
-    O=[1.0,2.0,10.0]
-    S,N = denoise(O)
-    @test all(O.==S)
-    @test all(N.==0)
+    O = [1.0, 2.0, 10.0]
+    S, N = denoise(O)
+    @test all(O .== S)
+    @test all(N .== 0)
     @test length(S) == 3
     @test length(N) == 3
 end
 
 @testset "Infinite smoothing" begin
-    O=sign.(sin.(1:1000)).+1
-    S,N = denoise(O,factor=Inf)
+    O = sign.(sin.(1:1000)) .+ 1
+    S, N = denoise(O, factor = Inf)
     @test all(S .≈ 1)
     @test all(abs.(N) .≈ 1)
 end
 
 @testset "Verbose" begin
-    O=sign.(sin.(1:1000)).+1
-    S,N = denoise(O, factor=0.0, verbose=true)
-    @test all(abs.(S.-O) .< 1e-10)
+    O = sign.(sin.(1:1000)) .+ 1
+    S, N = denoise(O, factor = 0.0, verbose = true)
+    @test all(abs.(S .- O) .< 1e-10)
+end
+
+@testset "Fit 1D RBF" begin
+    t = LinRange(0,2pi,1000)
+    y = sin.(t)
+    fn = fit_rbf(t,y,LinRange(0,2pi,20))
+    pred_y = fn(t)
+    error = sqrt(sum(abs2, pred_y .- y)/length(t))
+    @test error < 0.0005
 end