Add Convergence information to profiles

src/DataStructures/DataModel.jl

@@ -156,7 +156,7 @@ function (::Type{T})(DM::DataModel; kwargs...) where T<:Number
         @warn "Was unable to convert $(typeof(Data(DM))) to $T due to: $err"
         Data(DM)
     end
-    DataModel(D, Predictor(DM), dPredictor(DM), T.(MLE(DM)); kwargs...)
+    DataModel(D, Predictor(DM), dPredictor(DM), T.(MLE(DM)), true; kwargs...)
 end
 
 

src/InformationGeometry.jl

@@ -18,6 +18,7 @@ using TreeViews, PrettyTables
 
 import DataInterpolations: AbstractInterpolation
 import SciMLBase: AbstractODESolution, AbstractODEFunction, AbstractODEAlgorithm, AbstractADType
+import ModelingToolkit: AbstractODESystem
 
 import SciMLBase: remake
 export remake

src/ModelMaps.jl

@@ -637,7 +637,7 @@ Log10Transform(Sys::ODESystem, idxs::AbstractVector{<:Bool}=trues(length(paramet
     SystemTransform(Sys::ODESystem, F::Function, idxs::AbstractVector{<:Bool}=trues(length(parameters(Sys)))) -> ODESystem
 Transforms the parameters of a `ODESystem` according to `F`.
 """
-function SystemTransform(Sys::ODESystem, F::Function, idxs::AbstractVector{<:Bool}=trues(length(parameters(Sys))))
+function SystemTransform(Sys::AbstractODESystem, F::Function, idxs::AbstractVector{<:Bool}=trues(length(parameters(Sys))))
     SubstDict = Dict(parameters(Sys) .=> [(idxs[i] ? F(x) : x) for (i,x) in enumerate(parameters(Sys))])
     NewEqs = [(equations(Sys)[i].lhs ~ substitute(equations(Sys)[i].rhs, SubstDict)) for i in eachindex(equations(Sys))]
     # renamed "states" to "unknowns": https://github.com/SciML/ModelingToolkit.jl/pull/2432

src/ProfileLikelihood.jl

@@ -2,6 +2,7 @@
 
 # Returns a copy of type `Vector`, i.e. is not typesafe!
 SafeCopy(X::AbstractVector) = copy(X)
+SafeCopy(X::AbstractVector{<:Num}) = X
 SafeCopy(X::AbstractRange) = collect(X)
 SafeCopy(X::Union{SVector,MVector}) = convert(Vector,X)
 
@@ -210,6 +211,15 @@ end
 
 # USE NelderMead for ODEmodels!!!!!
 
+GetMinimizer(Res::LsqFit.LsqFitResult) = Res.param
+GetMinimum(Res::LsqFit.LsqFitResult, L::Function) = L(GetMinimizer(Res))
+HasConverged(Res::LsqFit.LsqFitResult) = Res.converged
+GetMinimizer(Res::Optim.OptimizationResults) = Optim.minimizer(Res)
+GetMinimum(Res::Optim.OptimizationResults, L::Function) = Res.minimum
+HasConverged(Res::Optim.OptimizationResults) = Optim.converged(Res)
+GetMinimizer(Res::SciMLBase.OptimizationSolution) = Res.u
+GetMinimum(Res::SciMLBase.OptimizationSolution, L::Function) = Res.objective
+HasConverged(Res::SciMLBase.OptimizationSolution) = Res.retcode === ReturnCode.Success
 
 """
     GetProfile(DM::AbstractDataModel, Comp::Int, dom::Tuple{<:Real, <:Real}; N::Int=50, dof::Int=DOF(DM), SaveTrajectories::Bool=false, SavePriors::Bool=false)
@@ -227,9 +237,9 @@ function GetProfile(DM::AbstractDataModel, Comp::Int, ps::AbstractVector{<:Real}
 
     FitFunc = if !isnothing(meth) || general || !isnothing(LogPrior(DM)) || Data(DM) isa AbstractUnknownUncertaintyDataSet
         Meth = isnothing(meth) ? NewtonTrustRegion() : meth
-        ((args...; Kwargs...)->InformationGeometry.minimize(args...; Full=false, tol=tol, meth=Meth, Kwargs...))
+        ((args...; Kwargs...)->InformationGeometry.minimize(args...; tol=tol, meth=Meth, Kwargs..., Full=true))
     else 
-        ((args...; Kwargs...)->curve_fit(args...; tol=tol, Kwargs...).param)
+        ((args...; Kwargs...)->curve_fit(args...; tol=tol, Kwargs...))
     end
 
     # Could use variable size array instead to cut off computation once Confnum+0.1 is reached?
@@ -238,6 +248,7 @@ function GetProfile(DM::AbstractDataModel, Comp::Int, ps::AbstractVector{<:Real}
     path = SaveTrajectories ? Vector{eltype(MLE(DM))}[] : nothing
     priors = SavePriors ? eltype(MLE(DM))[] : nothing
     Converged = BitVector()
+
     if pdim(DM) == 1    # Cannot drop dims if pdim already 1
         visitedps = ps
         Res = map(loglikelihood(DM), [[x] for x in ps])
@@ -262,16 +273,19 @@ function GetProfile(DM::AbstractDataModel, Comp::Int, ps::AbstractVector{<:Real}
             PerformStep!!! = if general || Data(DM) isa AbstractUnknownUncertaintyDataSet
                 function PerformStepGeneral!(Res, MLEstash, Converged, p)
                     L = Negloglikelihood(DM)∘ValInserter(Comp, p, MLE(DM))
-                    # R = FitFunc(L, MLEstash; kwargs...)
-                    copyto!(MLEstash, FitFunc(L, MLEstash; kwargs...))
-                    push!(Res, -L(MLEstash))
+                    R = FitFunc(L, MLEstash; kwargs...)
+                    copyto!(MLEstash, GetMinimizer(R))
+                    push!(Res, -GetMinimum(R,L))
+                    push!(Converged, HasConverged(R))
                 end
             else
                 function PerformStepManual!(Res, MLEstash, Converged, p)
                     NewModel = ProfilePredictor(DM, Comp, p)
                     DroppedLogPrior = EmbedLogPrior(DM, ValInserter(Comp, p, MLE(DM)))
-                    copyto!(MLEstash, FitFunc(Data(DM), NewModel, MLEstash, DroppedLogPrior; kwargs...))
-                    push!(Res, loglikelihood(Data(DM), NewModel, MLEstash, DroppedLogPrior))
+                    R = FitFunc(Data(DM), NewModel, MLEstash, DroppedLogPrior; kwargs...)
+                    copyto!(MLEstash, GetMinimizer(R))
+                    push!(Res, -GetMinimum(R,x->loglikelihood(Data(DM), NewModel, x, DroppedLogPrior)))
+                    push!(Converged, HasConverged(R))
                 end
             end
             # Adaptive instead of ps grid here?
@@ -305,7 +319,7 @@ function GetProfile(DM::AbstractDataModel, Comp::Int, ps::AbstractVector{<:Real}
         end
     end
 
-    ResMat = SavePriors ? [visitedps Res priors] : [visitedps Res]
+    ResMat = SavePriors ? [visitedps Res priors Converged] : [visitedps Res Converged]
     SaveTrajectories ? (ResMat, path) : ResMat
 end
 
@@ -318,7 +332,7 @@ end
 function GetLocalProfileDir(DM::AbstractDataModel, Comp::Int, p::AbstractVector{<:Number}=MLE(DM))
     F = FisherMetric(DM, p)
     F[Comp, :] .= [(j == Comp) for j in eachindex(p)]
-    isinf(logdet(F)) && @warn "Using pseudo-inverse to determine profile direction for parameter $Comp due to local non-identifiability."
+    det(F) == 0 && @warn "Using pseudo-inverse to determine profile direction for parameter $Comp due to local non-identifiability."
     pinv(F)[:, Comp]
 end
 
@@ -352,14 +366,21 @@ end
 PlotSingleProfile(DM::AbstractDataModel, Prof::Tuple{<:AbstractMatrix, <:Any}, i::Int; kwargs...) = PlotSingleProfile(DM, Prof[1], i; kwargs...)
 function PlotSingleProfile(DM::AbstractDataModel, Prof::AbstractMatrix, i::Int; kwargs...)
     P = RecipesBase.plot(view(Prof, :,1), view(Prof, :,2); leg=false, label="Profile", kwargs...)
-    size(Prof,2) == 3 && RecipesBase.plot!(P, view(Prof, :,1), view(Prof, :,3); label="Prior", color=:red, line=:dash)
+    HasPriors(Prof) && RecipesBase.plot!(P, view(Prof, :,1), view(Prof, :,3); label="Prior", color=:red, line=:dash)
     P
 end
 
+
+GetConverged(M::AbstractMatrix) = @view M[:, end]
+
 HasTrajectories(M::AbstractVector) = any(HasTrajectories, M)
 HasTrajectories(M::Tuple) = true
 HasTrajectories(M::AbstractMatrix) = false
 
+HasPriors(M::AbstractVector) = any(HasPriors, M)
+HasPriors(M::Tuple) = HasPriors(M[1])
+HasPriors(M::AbstractMatrix) = size(M,2) > 3
+
 function ProfilePlotter(DM::AbstractDataModel, Profiles::AbstractVector;
     Pnames::AbstractVector{<:AbstractString}=(Predictor(DM) isa ModelMap ? pnames(Predictor(DM)) : CreateSymbolNames(pdim(DM), "θ")), idxs=nothing, kwargs...)
     @assert length(Profiles) == length(Pnames)
@@ -463,7 +484,7 @@ mutable struct ParameterProfiles <: AbstractProfiles
         Trajs = SaveTrajectories ? getindex.(FullProfs,2) : fill(nothing, length(Profs))
         if !(inds == 1:pdim(DM))
             for i in 1:pdim(DM) # Profs and Trajs already sorted by sorting inds
-                i ∉ inds && (insert!(Profs, i, fill(NaN, size(Profs[1])));    SaveTrajectories ? insert!(Trajs, i, fill(NaN, pdim(DM))) : insert!(Trajs, i, nothing))
+                i ∉ inds && (insert!(Profs, i, fill(NaN, size(Profs[1])));  fill!(Profs[i][:,end], 0.);    SaveTrajectories ? insert!(Trajs, i, fill(NaN, pdim(DM))) : insert!(Trajs, i, nothing))
             end
         end
         P = ParameterProfiles(DM, Profs, Trajs; IsCost=IsCost)
@@ -560,7 +581,7 @@ end
             # P(i)
         end
         # Draw prior contribution
-        if size(Profiles(P)[i],2) == 3
+        if HasPriors(Profiles(P)[i])
             @series begin
                 label --> "Prior contribution"
                 color --> :red

test/runtests.jl

@@ -77,7 +77,7 @@ end
     @test ApproxInRegion(sol, MLE(ToyDME)) && !ApproxInRegion(sol, sol.u[1] + 1e-5BasisVector(1,2))
 
     #Check that bounding box from ProfileLikelihood coincides roughly with exact box.
-    Mats = ProfileLikelihood(ToyDME, 2; plot=false)
+    Mats = ParameterProfiles(ToyDME, 2; plot=false)
     ProfBox = ProfileBox(ToyDME, InterpolatedProfiles(Mats), 1)
     ExactBox = ConstructCube(sol)
     @test norm(Center(ProfBox) - Center(ExactBox)) < 3e-5 && norm(CubeWidths(ProfBox) - CubeWidths(ExactBox)) < 3e-4
@@ -250,7 +250,7 @@ end
     @test norm(Score(GDM, P) - Score(CDM, P)) < 2e-4
     @test norm(FisherMetric(GDM, P) - FisherMetric(CDM, P)) < 2e-4
 
-    UDS = DataSetUncertain(1:5, (1:5) + [rand(Normal(0,0.4)) for i in 1:5], (x,y,p)->1/abs(p[1]), [0.4])
+    UDS = DataSetUncertain(1:5, (1:5) + [rand(Normal(0,0.4)) for i in 1:5], (x,y,p)->1/exp(p[1]), [log(0.4)])
     UDM = DataModel(UDS, (x,p)->p[1]*x + p[2], [1, 1, 1.])
 
     # Test Type conversions for Datasets
@@ -294,9 +294,9 @@ end
 
     DS = DataSetUncertain(1:5, (1:5) + [rand(Normal(0,0.4)) for i in 1:5], (x,y,p)->1/abs(p[1]), [0.4]; xnames=["Time"], ynames=["Signal"])
     tester(
-        DataModel(DS, (x,p)->p[1]*x + p[2], [1, 1, 1.]; meth=Newton()),
-        DataModel(DataSet(xdata(DS), ydata(DS), 0.4; xnames=["Time"], ynames=["Signal"]), (x,p)->p[1].*x .+ p[2], [1, 1.]; meth=Newton()),
-        DataModel(DataSetExact(xdata(DS), ydata(DS), 0.4), (x,p)->p[1]*x + p[2], [1, 1.]; meth=Newton()),
+        DataModel(DS, (x,p)->p[1]*x + p[2], [1, 1, 1.]; meth=Optim.Newton()),
+        DataModel(DataSet(xdata(DS), ydata(DS), 0.4; xnames=["Time"], ynames=["Signal"]), (x,p)->p[1].*x .+ p[2], [1, 1.]; meth=Optim.Newton()),
+        DataModel(DataSetExact(xdata(DS), ydata(DS), 0.4), (x,p)->p[1]*x + p[2], [1, 1.]; meth=Optim.Newton()),
         [1,1,0.4]
     )
 
@@ -341,15 +341,15 @@ end
     rSTAT5A_rel = [14.72316822,33.76234229,36.79985129,49.71760229,46.9281201,47.83657456,46.92872725,40.59775294,43.78366389,44.45738765,41.32715926,41.06273321,39.23583003,36.61946054,34.8937144,32.21107716]
 
     BöhmDS = DataSet(BöhmT, [pSTAT5A_rel pSTAT5B_rel rSTAT5A_rel], [[4.12, 7.04, 3.37] for i in eachindex(BöhmT)]; xnames=["Time"], ynames=["pSTAT5A_rel", "pSTAT5B_rel", "rSTAT5A_rel"])
-    BöhmModel = GetModel(Exp10Transform(BöhmSys), BöhmInitial, BöhmObservation; meth=AutoTsit5(Rodas5()), tol=1e-6, Domain=HyperCube(-6ones(6), 6ones(6)))
-    BöhmDM = DataModel(BöhmDS, BöhmModel, [-1.5690, 4.1978, 4.99, -1.7859, -2.2, -5.], meth=LBFGS())
+    BöhmModel = GetModel(Exp10Transform(BöhmSys), BöhmInitial, BöhmObservation; tol=1e-4, Domain=HyperCube(-6ones(6), 6ones(6)))
+    BöhmDM = DataModel(BöhmDS, BöhmModel, [-1.5690, 4.1978, 4.99, -1.7859, -2.2, -5.]; tol=1e-6)
 
     Böhmds = DataSetUncertain(BöhmT, [pSTAT5A_rel pSTAT5B_rel rSTAT5A_rel], (x,y,c)->Diagonal(1 ./ abs.(c)), [4.12, 7.04, 3.37]; xnames=["Time"], ynames=["pSTAT5A_rel", "pSTAT5B_rel", "rSTAT5A_rel"])
-    Böhmdm = DataModel(Böhmds, GetModel(Exp10Transform(BöhmSys), BöhmInitial, BöhmObservation; meth=AutoTsit5(Rodas5()), tol=1e-6, Domain=vcat(HyperCube(-6ones(6), 6ones(6)), HyperCube(1e-2ones(3), 20ones(3)))),
-                            vcat([-1.5690, 4.1978, 4.99, -1.7859, -2.2, -5.], [4.12, 7.04, 3.37]), meth=NewtonTrustRegion())
+    Böhmdm = DataModel(Böhmds, GetModel(Exp10Transform(BöhmSys), BöhmInitial, BöhmObservation; tol=1e-4, Domain=vcat(HyperCube(-6ones(6), 6ones(6)), HyperCube(1e-2ones(3), 20ones(3)))),
+                            vcat([-1.5690, 4.1978, 4.99, -1.7859, -2.2, -5.], [4.12, 7.04, 3.37]); tol=1e-6)
     tester(
         Böhmdm, BöhmDM,
-        DataModel(DataSetExact(BöhmDS), Predictor(BöhmDM), MLE(BöhmDM)),
+        DataModel(DataSetExact(BöhmDS), Predictor(BöhmDM), MLE(BöhmDM); tol=1e-6),
         vcat(MLE(BöhmDM), [4.12, 7.04, 3.37])
     )
 end
@@ -481,19 +481,19 @@ end
 
     @test norm(minimize(F, initial, Cube; tol=1e-5, meth=NelderMead())) < 5e-1
     @test norm(minimize(F, initial, Cube; tol=1e-5, meth=LBFGS())) < 5e-2
-    @test norm(minimize(F, initial, Cube; tol=1e-5, meth=Newton())) < 5e-2
+    @test norm(minimize(F, initial, Cube; tol=1e-5, meth=Optim.Newton())) < 5e-2
 
     @test norm(minimize(F, initial; tol=1e-5, meth=NelderMead())) < 5e-1
     @test norm(minimize(F, initial; tol=1e-5, meth=LBFGS())) < 5e-2
-    @test norm(minimize(F, initial; tol=1e-5, meth=Newton())) < 5e-2
+    @test norm(minimize(F, initial; tol=1e-5, meth=Optim.Newton())) < 5e-2
 
     # Check optimization with non-linear constraints and box constraints
 
     # Check in-place and out-of-place optimization
 
     # Check type stability of optimization
     using ComponentArrays
-    @test minimize(X->X.A[1]^2 + 0.5X.B[1]^4, ComponentVector(A=[initial[1]], B=[initial[1]]); tol=1e-5, meth=Newton()) isa ComponentVector
+    @test minimize(X->X.A[1]^2 + 0.5X.B[1]^4, ComponentVector(A=[initial[1]], B=[initial[1]]); tol=1e-5, meth=Optim.Newton()) isa ComponentVector
 end