Merge pull request #98 from DanielVandH/deps

Up the package compats. Use SymbolicIndexingInterface

Project.toml

@@ -1,7 +1,7 @@
 name = "ProfileLikelihood"
 uuid = "3fca794e-44fa-49a6-b6a0-8cd45572ba0a"
 authors = ["Daniel VandenHeuvel <[email protected]>"]
-version = "0.3.4"
+version = "0.4.0"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
@@ -17,6 +17,7 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 
 [weakdeps]
 DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
@@ -30,14 +31,15 @@ ProfileLikelihoodMakieExt = "Makie"
 ChunkSplitters = "1.0, 2.0"
 Contour = "0.6"
 FunctionWrappers = "1.1"
-Interpolations = "0.14"
+Interpolations = "0.14, 0.15"
 InvertedIndices = "1.2, 1.3"
 OffsetArrays = "1.12"
 PolygonInbounds = "0.2"
 PreallocationTools = "0.4"
 Requires = "1.3"
-SciMLBase = "=2.9.1"
-SimpleNonlinearSolve = "0.1"
+SciMLBase = "2.0"
+SimpleNonlinearSolve = "1"
+SymbolicIndexingInterface = "0.3"
 StatsFuns = "1.1, 1.3"
 julia = "1"
 

ext/ProfileLikelihoodDelaunayTriangulationExt.jl

@@ -13,8 +13,8 @@ function ProfileLikelihood._get_confidence_regions_delaunay!(confidence_regions,
     grid_xy = vec([(x, y) for x in range_1, y in range_2])
     tri = DelaunayTriangulation.triangulate(grid_xy)
     conf_contour = NTuple{2,T}[]
-    DelaunayTriangulation.delete_boundary_vertices_from_graph!(tri)
-    for (u, v) in DelaunayTriangulation.get_edges(tri)
+    for e in DelaunayTriangulation.each_solid_edge(tri)
+        u, v = DelaunayTriangulation.edge_vertices(e)
         u₁, u₂ = profile_values[u], profile_values[v]
         if threshold_intersection_exists(threshold, u₁, u₂)
             p₁ = grid_xy[u]

ext/ProfileLikelihoodMakieExt.jl

@@ -1,6 +1,6 @@
 module ProfileLikelihoodMakieExt
 
- using ProfileLikelihood
+using ProfileLikelihood
 @static if isdefined(Base, :get_extension)
     using Makie
 else
@@ -21,7 +21,7 @@ function ProfileLikelihood.plot_profiles(prof::ProfileLikelihood.ProfileLikeliho
     axis_kwargs=nothing,
     show_points=false,
     markersize=9,
-    latex_names=Dict(vars .=> [L"\theta_{%$i}" for i in ProfileLikelihood.SciMLBase.sym_to_index.(vars, Ref(prof))]))
+    latex_names=Dict(vars .=> [L"\theta_{%$i}" for i in ProfileLikelihood.variable_index.((prof,), vars)]))
     num_plots = vars isa Symbol ? 1 : length(vars)
     _, _, plot_positions = ProfileLikelihood.choose_grid_layout(num_plots, ncol, nrow)
     if fig_kwargs !== nothing
@@ -51,7 +51,7 @@ function ProfileLikelihood.plot_profiles(prof::ProfileLikelihood.BivariateProfil
     interpolation=false,
     smooth_confidence_boundary=false,
     close_contour=true,
-    latex_names=Dict(1:ProfileLikelihood.number_of_parameters(ProfileLikelihood.get_likelihood_problem(prof)) .=> ProfileLikelihood.get_syms(prof)),
+    latex_names=Dict(1:ProfileLikelihood.number_of_parameters(ProfileLikelihood.get_likelihood_problem(prof)) .=> ProfileLikelihood.variable_symbols(prof)),
     xlim_tuples=nothing,
     ylim_tuples=nothing)
     vars = ProfileLikelihood.convert_symbol_tuples(vars, prof)
@@ -120,26 +120,26 @@ function ProfileLikelihood.plot_profiles!(prof::ProfileLikelihood.ProfileLikelih
     Makie.ylims!(ax, threshold - 1, 0.1)
     if !spline
         Makie.lines!(ax, θ_vals, ℓ_vals)
-        Makie.hlines!(ax, [threshold], color=:red, linetype=:dashed)
+        Makie.hlines!(ax, [threshold], color=:red, linestyle=:dash)
         CI_range = lower_ci .< θ_vals .< upper_ci
         shade_ci && Makie.band!(ax, θ_vals[CI_range], ℓ_vals[CI_range], repeat([threshold], count(CI_range)), color=(:blue, 0.35))
     else
         val_range = extrema(θ_vals)
         Δθ₁ = (val_range[2] - val_range[1]) / max(length(θ_vals), 1000)
         data_vals = val_range[1]:Δθ₁:val_range[2]
         Makie.lines!(ax, data_vals, prof(data_vals))
-        Makie.hlines!(ax, [threshold], color=:red, linetype=:dashed)
+        Makie.hlines!(ax, [threshold], color=:red, linestyle=:dash)
         Δθ₂ = (upper_ci - lower_ci) / max(length(θ_vals), 1000)
         if Δθ₂ ≠ 0.0
             ci_vals = lower_ci:Δθ₂:upper_ci
             shade_ci && Makie.band!(ax, ci_vals, prof(ci_vals), repeat([threshold], length(ci_vals)), color=(:blue, 0.35))
         end
     end
     if !isnothing(true_vals)
-        Makie.vlines!(ax, [true_vals], color=:black, linetype=:dashed)
+        Makie.vlines!(ax, [true_vals], color=:black, linestyle=:dash)
     end
     if !isnothing(mle_val)
-        Makie.vlines!(ax, [mle_val], color=:red, linetype=:dashed)
+        Makie.vlines!(ax, [mle_val], color=:red, linestyle=:dash)
     end
     if show_points
         Makie.scatter!(ax, θ_vals, ℓ_vals, color=:black, markersize=markersize)

src/ProfileLikelihood.jl

@@ -12,6 +12,7 @@ using OffsetArrays
 using PolygonInbounds
 using Contour
 using ChunkSplitters
+using SymbolicIndexingInterface
 
 include("utils.jl")
 include("problem_updates.jl")

src/abstract_type_definitions.jl

@@ -31,24 +31,10 @@ finite_bounds(prob) = finite_lower_bounds(prob) && finite_upper_bounds(prob)
 number_of_parameters(prob::OptimizationProblem) = length(prob.u0)
 number_of_parameters(::AbstractLikelihoodProblem{N, L}) where {N, L} = N
 
-Base.getindex(prob::AbstractLikelihoodProblem, i::Integer) = get_θ₀(prob, i)
-SciMLBase.sym_to_index(sym, prob::AbstractLikelihoodProblem) = SciMLBase.sym_to_index(sym, get_syms(prob))
-SciMLBase.sym_to_index(i::Integer, ::AbstractLikelihoodProblem) = i
+SymbolicIndexingInterface.symbolic_container(prob::AbstractLikelihoodProblem) = get_problem(prob)
+SymbolicIndexingInterface.state_values(prob::AbstractLikelihoodProblem) = get_θ₀(prob)
 function Base.getindex(prob::AbstractLikelihoodProblem, sym)
-    if SciMLBase.issymbollike(sym)
-        i = SciMLBase.sym_to_index(sym, prob)
-    else
-        i = sym
-    end
-    if i === nothing
-        throw(BoundsError(get_θ₀(prob), sym))
-    else
-        get_θ₀(prob, i)
-    end
-end
-function Base.getindex(prob::AbstractLikelihoodProblem, sym::AbstractVector{Symbol})
-    idx = SciMLBase.sym_to_index.(sym, Ref(prob))
-    return get_θ₀(prob, idx)
+    return getu(prob, sym)(prob)
 end
 
 function parameter_is_inbounds(prob, θ)
@@ -80,23 +66,10 @@ get_syms(sol::AbstractLikelihoodSolution) = get_syms(get_problem(sol))
 
 number_of_parameters(::AbstractLikelihoodSolution{N,P}) where {N,P}=N
 
-Base.getindex(sol::AbstractLikelihoodSolution, i::Integer) = get_mle(sol, i)
-SciMLBase.sym_to_index(sym, sol::AbstractLikelihoodSolution) = SciMLBase.sym_to_index(sym, get_syms(sol))
+SymbolicIndexingInterface.symbolic_container(sol::AbstractLikelihoodSolution) = get_problem(sol)
+SymbolicIndexingInterface.state_values(sol::AbstractLikelihoodSolution) = get_mle(sol)
 function Base.getindex(sol::AbstractLikelihoodSolution, sym)
-    if SciMLBase.issymbollike(sym)
-        i = SciMLBase.sym_to_index(sym, sol)
-    else
-        i = sym
-    end
-    if i === nothing
-        throw(BoundsError(get_mle(sol), sym))
-    else
-        get_mle(sol, i)
-    end
-end
-function Base.getindex(sol::AbstractLikelihoodSolution, sym::AbstractVector{Symbol})
-    idx = SciMLBase.sym_to_index.(sym, Ref(sol))
-    return get_mle(sol, idx)
+    return getu(sol, sym)(sol)
 end
 
 update_initial_estimate(prob::AbstractLikelihoodProblem, sol::AbstractLikelihoodSolution) = update_initial_estimate(prob, get_mle(sol))

src/display.jl

@@ -6,7 +6,7 @@ function Base.show(io::IO, ::MIME"text/plain", prob::T) where {T<:AbstractLikeli
         type_color, isinplace(get_problem(prob)))
     println(io,
         no_color, "θ₀: ", summary(prob.θ₀))
-    sym_param_names = get_syms(prob)
+    sym_param_names = variable_symbols(prob)
     for i in 1:number_of_parameters(prob)
         i < number_of_parameters(prob) ? println(io,
             no_color, "     $(sym_param_names[i]): $(prob[i])") : print(io,
@@ -35,7 +35,7 @@ function Base.show(io::IO, mime::MIME"text/plain", sol::T) where {T<:AbstractLik
     println(io)
     println(io,
         no_color, "Maximum likelihood estimates: ", summary(get_mle(sol)))
-    sym_param_names = get_syms(sol)
+    sym_param_names = variable_symbols(sol)
     for i in 1:number_of_parameters(sol)
         i < number_of_parameters(sol) ? println(io,
             no_color, "     $(sym_param_names[i]): $(sol[i])") : print(io,
@@ -61,7 +61,7 @@ function Base.show(io::IO, ::MIME"text/plain", prof::T) where {T<:ProfileLikelih
     println(io,
         no_color, "Confidence intervals: ")
     CIs = get_confidence_intervals(prof)
-    sym_param_names = get_syms(prof)
+    sym_param_names = variable_symbols(prof)
     for i in profiled_parameters(prof)
         i ≠ last(profiled_parameters(prof)) ? println(io,
             no_color, "     $(100get_level(CIs[i]))% CI for $(sym_param_names[i]): ($(get_lower(CIs[i])), $(get_upper(CIs[i])))") : print(io,
@@ -71,7 +71,7 @@ end
 
 function Base.show(io::IO, ::MIME"text/plain", prof::ProfileLikelihoodSolutionView{N,PLS}) where {N,PLS}
     type_color, no_color = SciMLBase.get_colorizers(io)
-    param_name = get_syms(prof)
+    param_name = variable_symbols(prof)[N]
     CIs = get_confidence_intervals(prof)
     println(io,
         type_color, "Profile likelihood",
@@ -94,7 +94,7 @@ function Base.show(io::IO, ::MIME"text/plain", prof::T) where {T<:BivariateProfi
     println(io,
         no_color, "Profile info: ")
     pairs = profiled_parameters(prof)
-    sym_param_names = get_syms(prof)
+    sym_param_names = variable_symbols(prof)
     for (i, j) in pairs
         a, b, c, d = get_bounding_box(prof, i, j)
         num_layers = number_of_layers(prof, i, j)

src/likelihood_problem.jl

@@ -26,7 +26,7 @@ The log-likelihood function, taking the form `ℓ(θ, p)`.
 Initial estimates for the MLE `θ`.
 - `syms::S`
 
-Variable names for the parameters.
+Variable names for the parameters. This is wrapped into a `SymbolCache` from SymbolicIndexingInterface.jl.
     
 The extra parameter `N` is the number of parameters.
 
@@ -106,16 +106,20 @@ Returns the [`LikelihoodProblem`](@ref) problem object.
 """
 function LikelihoodProblem end
 
+_to_symbolcache(syms, θ₀) = SymbolCache(syms, defaults = Dict(syms .=> θ₀))
+_to_symbolcache(syms::SymbolCache, θ₀) = syms
+
 function LikelihoodProblem(loglik::Function, θ₀;
     syms=eachindex(θ₀), data=SciMLBase.NullParameters(),
     f_kwargs=nothing, prob_kwargs=nothing)
     Base.require_one_based_indexing(θ₀)
+    _syms = _to_symbolcache(syms, θ₀)
     negloglik = negate_loglik(loglik)
-    opt_f = isnothing(f_kwargs) ? construct_optimisation_function(negloglik, syms) : construct_optimisation_function(negloglik, syms; f_kwargs...)
+    opt_f = isnothing(f_kwargs) ? construct_optimisation_function(negloglik, _syms) : construct_optimisation_function(negloglik, _syms; f_kwargs...)
     opt_prob = isnothing(prob_kwargs) ? construct_optimisation_problem(opt_f, θ₀, data) : construct_optimisation_problem(opt_f, θ₀, data; prob_kwargs...)
     return LikelihoodProblem{length(θ₀),typeof(opt_prob),
         typeof(data),typeof(loglik),
-        typeof(θ₀),typeof(syms)}(opt_prob, data, loglik, θ₀, syms)
+        typeof(θ₀),typeof(_syms)}(opt_prob, data, loglik, θ₀, _syms)
 end
 
 function LikelihoodProblem(loglik::Function, θ₀, integrator;
@@ -140,9 +144,9 @@ end
 
 function construct_optimisation_function(negloglik, syms; f_kwargs...)
     if :adtype ∈ keys(f_kwargs)
-        return OptimizationFunction(negloglik, f_kwargs[:adtype]; syms=syms, f_kwargs[Not(:adtype)]...)
+        return OptimizationFunction(negloglik, f_kwargs[:adtype]; sys=syms, f_kwargs[Not(:adtype)]...)
     else
-        return OptimizationFunction(negloglik, SciMLBase.NoAD(); syms=syms, f_kwargs...)
+        return OptimizationFunction(negloglik, SciMLBase.NoAD(); sys=syms, f_kwargs...)
     end
 end
 function construct_optimisation_problem(opt_f, θ₀, data; prob_kwargs...)

src/problem_updates.jl

@@ -1,35 +1,19 @@
 update_initial_estimate(prob::OptimizationProblem, θ) = remake(prob; u0=θ)
 update_initial_estimate(prob::OptimizationProblem, sol::SciMLBase.OptimizationSolution) = update_initial_estimate(prob, sol.u)
 
+@generated function _to_namedtuple(obj)
+    A = (Expr(:(=), n, :(obj.$n)) for n in setdiff(fieldnames(obj), (:f, :adtype, :grad, :hess)))
+    Expr(:tuple, A...)
+end
+
 # replace obj with a new obj
-@inline function replace_objective_function(f::OF, new_f::FF, new_grad!, new_hess!) where {OF, FF}
-    # scimlbase needs to add a constructorof method for OptimizationFunction before we can just do @set with accessors.jl. 
-    # g = @set f.f = new_f
-    # return g
+@inline function replace_objective_function(f::OF, new_f::FF, new_grad!::GG, new_hess!::HH) where {OF, FF, GG, HH}
     return OptimizationFunction(
         new_f,
         f.adtype;
         grad=new_grad!,
         hess=new_hess!,
-        hv=f.hv,
-        cons=f.cons,
-        cons_j=f.cons_j,
-        cons_h=f.cons_h,
-        hess_prototype=f.hess_prototype,
-        cons_jac_prototype=f.cons_jac_prototype,
-        f.cons_hess_prototype,
-        f.syms,
-        f.paramsyms,
-        f.observed,
-        f.expr,
-        f.cons_expr,
-        f.sys,
-        f.lag_h,
-        f.lag_hess_prototype,
-        f.hess_colorvec,
-        f.cons_jac_colorvec,
-        f.cons_hess_colorvec,
-        f.lag_hess_colorvec,
+        _to_namedtuple(f)...
     )
 end
 

src/profile_likelihood.jl

@@ -719,8 +719,8 @@ function _get_confidence_intervals_spline!(confidence_intervals, n, param_vals,
     right_bracket = (mles[n], param_vals[end])
     left_prob = IntervalNonlinearProblem(itp_f, left_bracket)
     right_prob = IntervalNonlinearProblem(itp_f, right_bracket)
-    ℓ = solve(left_prob, Falsi()).u
-    u = solve(right_prob, Falsi()).u
+    ℓ = solve(left_prob, ITP()).u # Used to be Falsi(), but Falsi() just sometimes doesn't want to give good results anymore...
+    u = solve(right_prob, ITP()).u
     confidence_intervals[n] = ConfidenceInterval(ℓ, u, conf_level)
     return nothing
 end