From ec5aa1ffa62cb9ff356dea93d1ffe0734fc63813 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan <35105271+sathvikbhagavan@users.noreply.github.com> Date: Sat, 2 Sep 2023 17:23:31 +0530 Subject: [PATCH 1/6] fix: kernel functions `abs(t) > 0` will give all points except 0. This means K(t) will result in 0 for all points except 0. That is incorrect. --- src/collocation.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/collocation.jl b/src/collocation.jl index 105d62b02..6f1bd0831 100644 --- a/src/collocation.jl +++ b/src/collocation.jl @@ -36,7 +36,7 @@ function calckernel(::TriangularKernel,t) end function calckernel(::QuarticKernel,t) - if abs(t)>0 + if abs(t) > 1 return 0 else return (15*(1-t^2)^2)/16 @@ -44,7 +44,7 @@ function calckernel(::QuarticKernel,t) end function calckernel(::TriweightKernel,t) - if abs(t)>0 + if abs(t) > 1 return 0 else return (35*(1-t^2)^3)/32 @@ -52,7 +52,7 @@ function calckernel(::TriweightKernel,t) end function calckernel(::TricubeKernel,t) - if abs(t)>0 + if abs(t) > 1 return 0 else return (70*(1-abs(t)^3)^3)/80 @@ -64,7 +64,7 @@ function calckernel(::GaussianKernel,t) end function calckernel(::CosineKernel,t) - if abs(t)>0 + if abs(t) > 1 return 0 else return (π*cos(π*t/2))/4 From 36c95c501de3d6db1fafb37edbed71ec0f16bf01 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Mon, 4 Sep 2023 06:08:55 +0000 Subject: [PATCH 2/6] test: add unit tests for kernel functions --- test/collocation.jl | 45 +++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 9 +++++++-- 2 files changed, 52 insertions(+), 2 deletions(-) create mode 100644 test/collocation.jl diff --git a/test/collocation.jl b/test/collocation.jl new file mode 100644 index 000000000..9f5a770a0 --- /dev/null +++ b/test/collocation.jl @@ -0,0 +1,45 @@ +using DiffEqFlux +using Test + +ts = collect(-5.0:0.1:5.0) + +@testset "Kernels with support from -1 to 1" begin + minus_one_index = findfirst(x -> ==(x, -1.0), ts) + plus_one_index = findfirst(x -> ==(x, 1.0), ts) + @testset "$kernel" for (kernel, x0) in zip( + [ + EpanechnikovKernel(), + UniformKernel(), + TriangularKernel(), + QuarticKernel(), + TriweightKernel(), + TricubeKernel(), + CosineKernel(), + ], + [0.75, 0.50, 1.0, 15.0 / 16.0, 35.0 / 32.0, 70.0 / 81.0, pi / 4.0], + ) + ws = DiffEqFlux.calckernel.((kernel,), ts) + + # t < -1 + @test all(ws[1:minus_one_index-1] .== 0.0) + + # t > 1 + @test all(ws[plus_one_index+1:end] .== 0.0) + + # -1 < t <1 + @test all(ws[minus_one_index+1:plus_one_index-1] .> 0.0) + + # t = 0 + @test DiffEqFlux.calckernel(kernel, 0.0) == x0 + end +end + +@testset "Kernels with unconstrained support" begin + @testset "$kernel" for (kernel, x0) in zip( + [GaussianKernel(), LogisticKernel(), SigmoidKernel(), SilvermanKernel()], + [1 / (sqrt(2 * pi)), 0.25, 1 / pi, 1 / (2 * sqrt(2))], + ) + # t = 0 + @test DiffEqFlux.calckernel(kernel, 0.0) == x0 + end +end diff --git a/test/runtests.jl b/test/runtests.jl index e76a92d59..61e594945 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,8 +6,13 @@ const is_CI = haskey(ENV, "CI") @time begin if GROUP == "All" || GROUP == "DiffEqFlux" || GROUP == "Layers" - @safetestset "Collocation Regression" begin - include("collocation_regression.jl") + @safetestset "Collocation" begin + @safetestset "Unit tests" begin + include("collocation.jl") + end + @safetestset "Regression" begin + include("collocation_regression.jl") + end end @safetestset "Stiff Nested AD Tests" begin include("stiff_nested_ad.jl") From 8dfd4aa23374688c2a3579b6be14d29e30f447b8 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Mon, 4 Sep 2023 06:09:20 +0000 Subject: [PATCH 3/6] test: modify and expand kernel collocation tests --- test/collocation_regression.jl | 37 ++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 17 deletions(-) diff --git a/test/collocation_regression.jl b/test/collocation_regression.jl index d30721520..4c7102fa4 100644 --- a/test/collocation_regression.jl +++ b/test/collocation_regression.jl @@ -3,25 +3,28 @@ using DiffEqFlux, OrdinaryDiffEq, Test function f(u, p, t) p .* u end -rc=62 + +rc = 62 ps = repeat([-0.001], rc) -tspan = (7.0, 84.0) +tspan = (0.0, 50.0) u0 = 3.4 .+ ones(rc) -t = collect(range(minimum(tspan), stop=maximum(tspan), length=157)) +t = collect(range(minimum(tspan), stop = maximum(tspan), length = 1000)) prob = ODEProblem(f, u0, tspan, ps) -data = Array(solve(prob, Tsit5(), saveat = t, abstol=1e-12, reltol=1e-12)) -ptest = ones(rc) - -u′,u = collocate_data(data,t,SigmoidKernel()) -@test sum(abs2,u - data) < 1e-8 +data = Array(solve(prob, Tsit5(), saveat = t, abstol = 1e-12, reltol = 1e-12)) -function loss(p) - cost = zero(first(p)) - for i in 1:size(u′,2) - _du = f(@view(u[:,i]),p,t[i]) - u′i = @view u′[:,i] - cost += sum(abs2,u′i .- _du) - end - sqrt(cost) +@testset "$kernel" for kernel in [ + EpanechnikovKernel(), + UniformKernel(), + TriangularKernel(), + QuarticKernel(), + TriweightKernel(), + TricubeKernel(), + CosineKernel(), + GaussianKernel(), + LogisticKernel(), + SigmoidKernel(), + SilvermanKernel(), +] + u′, u = collocate_data(data, t, kernel) + @test sum(abs2, u - data) < 1e-8 end -@test loss(ptest) ≈ 418.3400017500223 From 61cc59113792f70e8da799968ca383460f88ad1e Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Mon, 4 Sep 2023 06:09:34 +0000 Subject: [PATCH 4/6] fix: `TricubeKernel` kernel --- src/collocation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/collocation.jl b/src/collocation.jl index 6f1bd0831..6fb94986c 100644 --- a/src/collocation.jl +++ b/src/collocation.jl @@ -55,7 +55,7 @@ function calckernel(::TricubeKernel,t) if abs(t) > 1 return 0 else - return (70*(1-abs(t)^3)^3)/80 + return (70*(1-abs(t)^3)^3)/81 end end From 724962d8f8b0e13d7192189f153fb4943b29ae55 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Tue, 19 Sep 2023 05:11:02 +0000 Subject: [PATCH 5/6] test: update tests for collocation 1. All tests in one file 2. Check collocation of data --- test/collocation.jl | 96 +++++++++++++++++++++------------- test/collocation_regression.jl | 30 ----------- test/runtests.jl | 7 +-- 3 files changed, 60 insertions(+), 73 deletions(-) delete mode 100644 test/collocation_regression.jl diff --git a/test/collocation.jl b/test/collocation.jl index 9f5a770a0..14cac4017 100644 --- a/test/collocation.jl +++ b/test/collocation.jl @@ -1,45 +1,67 @@ -using DiffEqFlux -using Test +using DiffEqFlux, OrdinaryDiffEq, Test -ts = collect(-5.0:0.1:5.0) +bounded_support_kernels = [ + EpanechnikovKernel(), + UniformKernel(), + TriangularKernel(), + QuarticKernel(), + TriweightKernel(), + TricubeKernel(), + CosineKernel(), +] -@testset "Kernels with support from -1 to 1" begin - minus_one_index = findfirst(x -> ==(x, -1.0), ts) - plus_one_index = findfirst(x -> ==(x, 1.0), ts) - @testset "$kernel" for (kernel, x0) in zip( - [ - EpanechnikovKernel(), - UniformKernel(), - TriangularKernel(), - QuarticKernel(), - TriweightKernel(), - TricubeKernel(), - CosineKernel(), - ], - [0.75, 0.50, 1.0, 15.0 / 16.0, 35.0 / 32.0, 70.0 / 81.0, pi / 4.0], - ) - ws = DiffEqFlux.calckernel.((kernel,), ts) +unbounded_support_kernels = + [GaussianKernel(), LogisticKernel(), SigmoidKernel(), SilvermanKernel()] - # t < -1 - @test all(ws[1:minus_one_index-1] .== 0.0) - - # t > 1 - @test all(ws[plus_one_index+1:end] .== 0.0) - - # -1 < t <1 - @test all(ws[minus_one_index+1:plus_one_index-1] .> 0.0) - - # t = 0 - @test DiffEqFlux.calckernel(kernel, 0.0) == x0 +@testset "Kernel Functions" begin + ts = collect(-5.0:0.1:5.0) + @testset "Kernels with support from -1 to 1" begin + minus_one_index = findfirst(x -> ==(x, -1.0), ts) + plus_one_index = findfirst(x -> ==(x, 1.0), ts) + @testset "$kernel" for (kernel, x0) in zip( + bounded_support_kernels, + [0.75, 0.50, 1.0, 15.0 / 16.0, 35.0 / 32.0, 70.0 / 81.0, pi / 4.0], + ) + ws = DiffEqFlux.calckernel.((kernel,), ts) + # t < -1 + @test all(ws[1:minus_one_index-1] .== 0.0) + # t > 1 + @test all(ws[plus_one_index+1:end] .== 0.0) + # -1 < t <1 + @test all(ws[minus_one_index+1:plus_one_index-1] .> 0.0) + # t = 0 + @test DiffEqFlux.calckernel(kernel, 0.0) == x0 + end + end + @testset "Kernels with unbounded support" begin + @testset "$kernel" for (kernel, x0) in zip( + unbounded_support_kernels, + [1 / (sqrt(2 * pi)), 0.25, 1 / pi, 1 / (2 * sqrt(2))], + ) + # t = 0 + @test DiffEqFlux.calckernel(kernel, 0.0) == x0 + end end end -@testset "Kernels with unconstrained support" begin - @testset "$kernel" for (kernel, x0) in zip( - [GaussianKernel(), LogisticKernel(), SigmoidKernel(), SilvermanKernel()], - [1 / (sqrt(2 * pi)), 0.25, 1 / pi, 1 / (2 * sqrt(2))], - ) - # t = 0 - @test DiffEqFlux.calckernel(kernel, 0.0) == x0 +@testset "Collocation of data" begin + function f(u, p, t) + p .* u + end + rc = 2 + ps = repeat([-0.001], rc) + tspan = (0.0, 50.0) + u0 = 3.4 .+ ones(rc) + t = collect(range(minimum(tspan), stop = maximum(tspan), length = 1000)) + prob = ODEProblem(f, u0, tspan, ps) + data = Array(solve(prob, Tsit5(), saveat = t, abstol = 1e-12, reltol = 1e-12)) + @testset "$kernel" for kernel in + [bounded_support_kernels..., unbounded_support_kernels...] + u′, u = collocate_data(data, t, kernel, 0.003) + @test sum(abs2, u - data) < 1e-8 + end + @testset "$kernel" for kernel in [bounded_support_kernels...] + # Errors out as the bandwidth is too low + @test_throws ErrorException collocate_data(data, t, kernel, 0.001) end end diff --git a/test/collocation_regression.jl b/test/collocation_regression.jl deleted file mode 100644 index 4c7102fa4..000000000 --- a/test/collocation_regression.jl +++ /dev/null @@ -1,30 +0,0 @@ -using DiffEqFlux, OrdinaryDiffEq, Test - -function f(u, p, t) - p .* u -end - -rc = 62 -ps = repeat([-0.001], rc) -tspan = (0.0, 50.0) -u0 = 3.4 .+ ones(rc) -t = collect(range(minimum(tspan), stop = maximum(tspan), length = 1000)) -prob = ODEProblem(f, u0, tspan, ps) -data = Array(solve(prob, Tsit5(), saveat = t, abstol = 1e-12, reltol = 1e-12)) - -@testset "$kernel" for kernel in [ - EpanechnikovKernel(), - UniformKernel(), - TriangularKernel(), - QuarticKernel(), - TriweightKernel(), - TricubeKernel(), - CosineKernel(), - GaussianKernel(), - LogisticKernel(), - SigmoidKernel(), - SilvermanKernel(), -] - u′, u = collocate_data(data, t, kernel) - @test sum(abs2, u - data) < 1e-8 -end diff --git a/test/runtests.jl b/test/runtests.jl index 61e594945..0bf6c449e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,12 +7,7 @@ const is_CI = haskey(ENV, "CI") @time begin if GROUP == "All" || GROUP == "DiffEqFlux" || GROUP == "Layers" @safetestset "Collocation" begin - @safetestset "Unit tests" begin - include("collocation.jl") - end - @safetestset "Regression" begin - include("collocation_regression.jl") - end + include("collocation.jl") end @safetestset "Stiff Nested AD Tests" begin include("stiff_nested_ad.jl") From 41d10318b5719475898af820d6510f4e0ac32832 Mon Sep 17 00:00:00 2001 From: Sathvik Bhagavan Date: Tue, 19 Sep 2023 05:12:04 +0000 Subject: [PATCH 6/6] refactor: make bandwidth a user passable argument and normalize time into kernel functions --- src/collocation.jl | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/collocation.jl b/src/collocation.jl index 6fb94986c..93f7a16b1 100644 --- a/src/collocation.jl +++ b/src/collocation.jl @@ -92,14 +92,14 @@ function construct_t2(t,tpoints) end function construct_w(t,tpoints,h,kernel) - W = @. calckernel((kernel,),(tpoints-t)/h)/h + W = @. calckernel((kernel,),((tpoints-t)/(tpoints[end]-tpoints[begin]))/h)/h Diagonal(W) end """ ```julia -u′,u = collocate_data(data,tpoints,kernel=SigmoidKernel()) +u′,u = collocate_data(data,tpoints,kernel=TriangularKernel(),bandwidth=nothing) u′,u = collocate_data(data,tpoints,tpoints_sample,interp,args...) ``` @@ -128,24 +128,29 @@ Additionally, we can use interpolation methods from data from intermediate timesteps. In this case, pass any of the methods like `QuadraticInterpolation` as `interp`, and the timestamps to sample from as `tpoints_sample`. """ -function collocate_data(data,tpoints,kernel=TriangularKernel()) +function collocate_data(data, tpoints, kernel=TriangularKernel(), bandwidth=nothing) _one = oneunit(first(data)) _zero = zero(first(data)) e1 = [_one;_zero] e2 = [_zero;_one;_zero] n = length(tpoints) - h = (n^(-1/5))*(n^(-3/35))*((log(n))^(-1/16)) + bandwidth = isnothing(bandwidth) ? (n^(-1/5))*(n^(-3/35))*((log(n))^(-1/16)) : bandwidth Wd = similar(data, n, size(data,1)) WT1 = similar(data, n, 2) WT2 = similar(data, n, 3) + T2WT2 = similar(data, 3, 3) + T1WT1 = similar(data, 2, 2) x = map(tpoints) do _t T1 = construct_t1(_t,tpoints) T2 = construct_t2(_t,tpoints) - W = construct_w(_t,tpoints,h,kernel) + W = construct_w(_t,tpoints,bandwidth,kernel) mul!(Wd,W,data') mul!(WT1,W,T1) mul!(WT2,W,T2) + mul!(T2WT2,T2',WT2) + mul!(T1WT1,T1',WT1) + (det(T2WT2) ≈ 0.0 || det(T1WT1) ≈ 0.0) && error("Collocation failed with bandwidth $bandwidth. Please choose a higher bandwidth") (e2'*((T2'*WT2)\T2'))*Wd,(e1'*((T1'*WT1)\T1'))*Wd end estimated_derivative = reduce(hcat,transpose.(first.(x))) @@ -163,7 +168,7 @@ function collocate_data(data::AbstractMatrix{T},tpoints::AbstractVector{T}, tpoints_sample::AbstractVector{T},interp,args...) where T u = zeros(T,size(data, 1),length(tpoints_sample)) du = zeros(T,size(data, 1),length(tpoints_sample)) - for d1 in 1:size(data,1) + for d1 in axes(data, 1) interpolation = interp(data[d1,:],tpoints,args...) u[d1,:] .= interpolation.(tpoints_sample) du[d1,:] .= DataInterpolations.derivative.((interpolation,), tpoints_sample)