From e5f2613308f0b1d09a4bb9e773aaec95654a6fe7 Mon Sep 17 00:00:00 2001 From: ArnoStrouwen Date: Sat, 15 Oct 2022 17:04:55 +0200 Subject: [PATCH] remove duplicate code for collocation --- Project.toml | 2 + docs/src/methods/collocation_loss.md | 2 +- src/DiffEqParamEstim.jl | 5 +- src/kernels.jl | 86 ---------------------------- src/two_stage_method.jl | 66 +++++++-------------- test/out_of_place_odes.jl | 4 +- 6 files changed, 28 insertions(+), 137 deletions(-) delete mode 100644 src/kernels.jl diff --git a/Project.toml b/Project.toml index 035f4296..b4a89313 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "1.26.0" Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" +DiffEqFlux = "aae7a2af-3d4f-5e19-a356-7da93b79d9d0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -21,6 +22,7 @@ SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" Calculus = "0.5" Dierckx = "0.4, 0.5" DiffEqBase = "6" +DiffEqFlux = "1" Distributions = "0.25" ForwardDiff = "0.10" LsqFit = "0.8, 0.9, 0.10, 0.11, 0.12" diff --git a/docs/src/methods/collocation_loss.md b/docs/src/methods/collocation_loss.md index 00dfe33f..ad0eabf7 100644 --- a/docs/src/methods/collocation_loss.md +++ b/docs/src/methods/collocation_loss.md @@ -10,7 +10,7 @@ but is much faster, and is a good method to try first to get in the general "good parameter" region, to then finish using one of the other methods. ```julia -function two_stage_method(prob::DEProblem,tpoints,data;kernel= :Epanechnikov, +function two_stage_method(prob::DEProblem,tpoints,data;kernel= EpanechnikovKernel(), loss_func = L2DistLoss,mpg_autodiff = false, verbose = false,verbose_steps = 100) ``` diff --git a/src/DiffEqParamEstim.jl b/src/DiffEqParamEstim.jl index 822c5b59..db516e94 100644 --- a/src/DiffEqParamEstim.jl +++ b/src/DiffEqParamEstim.jl @@ -1,9 +1,9 @@ module DiffEqParamEstim using DiffEqBase, LsqFit, PenaltyFunctions, RecursiveArrayTools, ForwardDiff, Calculus, Distributions, - LinearAlgebra, SciMLSensitivity, Dierckx, - SciMLBase + LinearAlgebra, SciMLSensitivity, Dierckx, DiffEqFlux, SciMLBase +import DiffEqFlux.CollocationKernel import PreallocationTools STANDARD_PROB_GENERATOR(prob, p) = remake(prob; u0 = eltype(p).(prob.u0), p = p) function STANDARD_PROB_GENERATOR(prob::EnsembleProblem, p) @@ -26,7 +26,6 @@ include("cost_functions.jl") include("lm_fit.jl") include("build_loss_objective.jl") include("build_lsoptim_objective.jl") -include("kernels.jl") include("two_stage_method.jl") include("multiple_shooting_objective.jl") diff --git a/src/kernels.jl b/src/kernels.jl deleted file mode 100644 index a59b2c3d..00000000 --- a/src/kernels.jl +++ /dev/null @@ -1,86 +0,0 @@ -# Kernel definition is taken from here: https://en.wikipedia.org/wiki/Kernel_(statistics)#Kernel_functions_in_common_use - -abstract type CollocationKernel end -struct EpanechnikovKernel <: CollocationKernel end -struct UniformKernel <: CollocationKernel end -struct TriangularKernel <: CollocationKernel end -struct QuarticKernel <: CollocationKernel end -struct TriweightKernel <: CollocationKernel end -struct TricubeKernel <: CollocationKernel end -struct GaussianKernel <: CollocationKernel end -struct CosineKernel <: CollocationKernel end -struct LogisticKernel <: CollocationKernel end -struct SigmoidKernel <: CollocationKernel end -struct SilvermanKernel <: CollocationKernel end - -function calckernel(::EpanechnikovKernel, t) - if abs(t) > 1 - return 0 - else - return 0.75 * (1 - t^2) - end -end - -function calckernel(::UniformKernel, t) - if abs(t) > 1 - return 0 - else - return 0.5 - end -end - -function calckernel(::TriangularKernel, t) - if abs(t) > 1 - return 0 - else - return (1 - abs(t)) - end -end - -function calckernel(::QuarticKernel, t) - if abs(t) > 0 - return 0 - else - return (15 * (1 - t^2)^2) / 16 - end -end - -function calckernel(::TriweightKernel, t) - if abs(t) > 0 - return 0 - else - return (35 * (1 - t^2)^3) / 32 - end -end - -function calckernel(::TricubeKernel, t) - if abs(t) > 0 - return 0 - else - return (70 * (1 - abs(t)^3)^3) / 80 - end -end - -function calckernel(::GaussianKernel, t) - exp(-0.5 * t^2) / (sqrt(2 * π)) -end - -function calckernel(::CosineKernel, t) - if abs(t) > 0 - return 0 - else - return (π * cos(π * t / 2)) / 4 - end -end - -function calckernel(::LogisticKernel, t) - 1 / (exp(t) + 2 + exp(-t)) -end - -function calckernel(::SigmoidKernel, t) - 2 / (π * (exp(t) + exp(-t))) -end - -function calckernel(::SilvermanKernel, t) - sin(abs(t) / 2 + π / 4) * 0.5 * exp(-abs(t) / sqrt(2)) -end diff --git a/src/two_stage_method.jl b/src/two_stage_method.jl index 6c6ddd29..121d73b8 100644 --- a/src/two_stage_method.jl +++ b/src/two_stage_method.jl @@ -1,4 +1,7 @@ export TwoStageCost, two_stage_method +export EpanechnikovKernel, UniformKernel, TriangularKernel, QuarticKernel +export TriweightKernel, TricubeKernel, GaussianKernel, CosineKernel +export LogisticKernel, SigmoidKernel, SilvermanKernel struct TwoStageCost{F, F2, D} <: Function cost_function::F @@ -10,7 +13,6 @@ end (f::TwoStageCost)(p) = f.cost_function(p) (f::TwoStageCost)(p, g) = f.cost_function2(p, g) -decide_kernel(kernel::CollocationKernel) = kernel function decide_kernel(kernel::Symbol) if kernel == :Epanechnikov return EpanechnikovKernel() @@ -32,51 +34,19 @@ function decide_kernel(kernel::Symbol) return LogisticKernel() elseif kernel == :Sigmoid return SigmoidKernel() - else + elseif kernel == :Silverman return SilvermanKernel() + else + return error("Kernel name not recognized") end end -function construct_t1(t, tpoints) - hcat(ones(eltype(tpoints), length(tpoints)), tpoints .- t) -end -function construct_t2(t, tpoints) - hcat(ones(eltype(tpoints), length(tpoints)), tpoints .- t, (tpoints .- t) .^ 2) -end -function construct_w(t, tpoints, h, kernel) - W = @. calckernel((kernel,), (tpoints - t) / h) / h - Diagonal(W) -end -function construct_estimated_solution_and_derivative!(data, kernel, tpoints) - _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)) - - Wd = similar(data, n, size(data, 1)) - WT1 = similar(data, n, 2) - WT2 = similar(data, n, 3) - x = map(tpoints) do _t - T1 = construct_t1(_t, tpoints) - T2 = construct_t2(_t, tpoints) - W = construct_w(_t, tpoints, h, kernel) - mul!(Wd, W, data') - mul!(WT1, W, T1) - mul!(WT2, W, T2) - (e2' * ((T2' * WT2) \ T2')) * Wd, (e1' * ((T1' * WT1) \ T1')) * Wd - end - estimated_derivative = reduce(hcat, transpose.(first.(x))) - estimated_solution = reduce(hcat, transpose.(last.(x))) - estimated_derivative, estimated_solution -end function construct_iip_cost_function(f, du, preview_est_sol, preview_est_deriv, tpoints) function (p) _du = PreallocationTools.get_tmp(du, p) vecdu = vec(_du) cost = zero(first(p)) - for i in 1:length(preview_est_sol) + for i in eachindex(preview_est_sol) est_sol = preview_est_sol[i] f(_du, est_sol, p, tpoints[i]) vecdu .= vec(preview_est_deriv[i]) .- vec(_du) @@ -89,7 +59,7 @@ end function construct_oop_cost_function(f, du, preview_est_sol, preview_est_deriv, tpoints) function (p) cost = zero(first(p)) - for i in 1:length(preview_est_sol) + for i in eachindex(preview_est_sol) est_sol = preview_est_sol[i] _du = f(est_sol, p, tpoints[i]) cost += sum(abs2, vec(preview_est_deriv[i]) .- vec(_du)) @@ -102,23 +72,27 @@ get_chunksize(cs) = cs get_chunksize(cs::Type{Val{CS}}) where {CS} = CS function two_stage_method(prob::DiffEqBase.DEProblem, tpoints, data; - kernel = EpanechnikovKernel(), + kernel::Union{CollocationKernel, Symbol} = EpanechnikovKernel(), loss_func = L2Loss, mpg_autodiff = false, verbose = false, verbose_steps = 100, autodiff_chunk = length(prob.p)) - f = prob.f - kernel_function = decide_kernel(kernel) - estimated_derivative, estimated_solution = construct_estimated_solution_and_derivative!(data, - kernel_function, - tpoints) + if kernel isa Symbol + @warn "Passing kernels as Symbols will be deprecated" + kernel = decide_kernel(kernel) + end + + # Step - 1 + + estimated_derivative, estimated_solution = collocate_data(data, tpoints, kernel) # Step - 2 du = PreallocationTools.dualcache(similar(prob.u0), autodiff_chunk) preview_est_sol = [@view estimated_solution[:, i] - for i in 1:size(estimated_solution, 2)] + for i in axes(estimated_solution, 2)] preview_est_deriv = [@view estimated_derivative[:, i] - for i in 1:size(estimated_solution, 2)] + for i in axes(estimated_solution, 2)] + f = prob.f if DiffEqBase.isinplace(prob) cost_function = construct_iip_cost_function(f, du, preview_est_sol, preview_est_deriv, tpoints) diff --git a/test/out_of_place_odes.jl b/test/out_of_place_odes.jl index d9536209..1e84dc79 100644 --- a/test/out_of_place_odes.jl +++ b/test/out_of_place_odes.jl @@ -46,7 +46,9 @@ prob_oop = ODEProblem{false}(ff, u0, tspan, ps) data = Array(solve(prob, Tsit5(), saveat = t)) ptest = ones(rc) -obj_ts = two_stage_method(prob, t, data; kernel = :Sigmoid) +obj_ts = two_stage_method(prob, t, data; kernel = SigmoidKernel()) +@test obj_ts(ptest) ≈ 418.3400017500223^2 +obj_ts = two_stage_method(prob_oop, t, data; kernel = SigmoidKernel()) @test obj_ts(ptest) ≈ 418.3400017500223^2 obj_ts = two_stage_method(prob_oop, t, data; kernel = :Sigmoid) @test obj_ts(ptest) ≈ 418.3400017500223^2