diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index 4536f845b..636a39ef8 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -131,6 +131,7 @@ include("global_identifiability.jl") include("identifiable_functions.jl") include("observation_field.jl") include("parametrizations.jl") +include("weights.jl") include("lincomp.jl") include("pb_representation.jl") include("submodels.jl") diff --git a/src/weights.jl b/src/weights.jl new file mode 100644 index 000000000..146499a3e --- /dev/null +++ b/src/weights.jl @@ -0,0 +1,62 @@ +""" + cmp_weighted(ode, W) + +Comparator with arbitrary weights. +""" +function cmp_weighted(ode, W::Matrix{Int}) + @assert all(>=(0), W) "Weights must be non-negative." + states, params = ode.x_vars, ode.parameters + df(f) = lie_derivative(f, ode) + features(f) = [ + !iszero(f), + total_degree(denominator(f)), + total_degree(numerator(f)), + length(numerator(f)) + length(denominator(f)), + !issubset(vars(f), params), + total_degree_wrt(numerator(f), states), + total_degree_wrt(denominator(f), states), + number_of_terms(f, states), + + # Lie derivative + !iszero(df(f)), + total_degree(numerator(df(f))), + total_degree(denominator(df(f))), + length(numerator(df(f))) + length(denominator(df(f))), + total_degree_wrt(df(f), states), + total_degree_wrt(denominator(df(f)), states), + number_of_terms(df(f), states), + ] + + fallback = rational_function_cmp + cmp = function (f, g) + f = parent_ring_change(f, parent(ode)) + g = parent_ring_change(g, parent(ode)) + v_f = features(f) + v_g = features(g) + score_f = W * v_f + score_g = W * v_g + score_f != score_g && return score_f < score_g + return fallback(f, g) + end + return cmp +end + +# Matrices of weights for existing named comparators. +const matrix_comparators = let + row(indices...) = [Int(j in indices) for _ in 1:1, j in 1:15] + ordinary = vcat(row(1), row(2, 3), row(4), row(2)) + state = vcat(row(6, 7), row(7), row(8)) + lie = vcat(row(9), row(10, 11), row(12), row(11)) + lie_state = vcat(row(13), row(14), row(15)) + parameter_first = vcat(row(5), ordinary) + + matrices = Dict( + cmp_default => parameter_first, + cmp_lie => lie, + cmp_states => vcat(state, ordinary), + cmp_lie_states => vcat(lie_state, lie), + cmp_pushka => vcat(state, lie_state, lie), + ) + nrows = maximum(size(W, 1) for W in values(matrices)) + Dict(cmp => vcat(W, zeros(Int, nrows - size(W, 1), 15)) for (cmp, W) in matrices) +end