# Checks that various model options (observables, events, defaults and metadata, differential equations,
# non-default_iv) works for hybrid models.
let
    # Creates the model (X species is pure jump, Y is pure ODE, and  Z1,Z2 are mixed).
    # Hybrid species have non-constant rates containing the two other species.
    rn = @reaction_network begin
        @ivs τ
        @differentials Δ = Differential(τ)
        @parameters X0 Y0
        @species X(τ) = X0 [description = "A jump only species"] Y(τ) = Y0 [description = "An ODE only species"]
        @variables Ztot(τ) = 0.0  # Default needed for MTK v10+ callback compilation
        @observables Ztot ~ Z1 + Z2
        @discrete_events [1.0] => [Z1 ~ Pre(Z1) + 1.0]
        @continuous_events [Y ~ 1.0] => [Y ~ 5.0]
        @equations Δ(V) ~ Z1 + X^2 - V
        (p,d), 0 <--> X
        d, Y --> 0, [physical_scale = PhysicalScale.ODE]
        (k*X, k*Y), Z1 <--> Z2, ([physical_scale = PhysicalScale.ODE], [physical_scale = PhysicalScale.Jump])
    end

    # Simulates the model.
    u0 = [:Z1 => 1.0, :Z2 => 2.0, :V => 1.5]
    ps = [:p => 2.0, :d => 0.5, :k => 1.0, :X0 => 0.1, :Y0 => 4.0]
    prob = HybridProblem(rn, u0, (0.0, 10.0), ps)
    sol = solve(prob, Tsit5(); tstops = [1.0 - eps(), 1.0 + eps()])

    # Tests that the model contain the correct stuff.
    @test ModelingToolkitBase.getdescription(rn.X) == "A jump only species"
    @test ModelingToolkitBase.getdescription(rn.Y) == "An ODE only species"
    @test sol[:X][1] == sol.ps[:X0] == 0.1
    @test sol[:Y][1] == sol.ps[:Y0] == 4.0
    @test sol[:V][1] == 1.5
    @test sol[:Ztot] ≈ sol[rn.Z1 + rn.Z2]
    @test minimum(sol[:Y]) ≈ 1.0
    @test maximum(sol[:Y]) ≈ 5.0 atol = 1e-1 rtol = 1e-1
    @test all(isequal([rn.τ], Symbolics.arguments(Symbolics.value(u))) for u in unknowns(rn))
    @test sol(1.0 - eps(); idxs = :Z1) + 1 ≈ sol(1.0 + eps(); idxs = :Z1)
end