feat: implement BVProblem for System

Author: AayushSabharwal

src/problems/bvproblem.jl

@@ -0,0 +1,34 @@
+@fallback_iip_specialize function SciMLBase.BVProblem{iip, spec}(
+        sys::System, u0map, tspan, parammap = SciMLBase.NullParameters();
+        check_compatibility = true, cse = true, checkbounds = false, eval_expression = false,
+        eval_module = @__MODULE__, guesses = Dict(), kwargs...) where {iip, spec}
+    check_complete(sys, BVProblem)
+    check_compatibility && check_compatible_system(BVProblem, sys)
+
+    # ODESystems without algebraic equations should use both fixed values + guesses
+    # for initialization.
+    _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses))
+    fode, u0, p = process_SciMLProblem(
+        ODEFunction{iip, spec}, sys, _u0map, parammap; guesses,
+        t = tspan !== nothing ? tspan[1] : tspan, check_compatibility = false, cse, checkbounds,
+        time_dependent_init = false, kwargs...)
+
+    dvs = unknowns(sys)
+    stidxmap = Dict([v => i for (i, v) in enumerate(dvs)])
+    u0_idxs = has_alg_eqs(sys) ? collect(1:length(dvs)) : [stidxmap[k] for (k, v) in u0map]
+    fbc = generate_boundary_conditions(
+        sys, u0, u0_idxs, tspan; expression = Val{false}, cse, checkbounds)
+    kwargs = process_kwargs(sys; kwargs...)
+    # Call `remake` so it runs initialization if it is trivial
+    return remake(BVProblem{iip}(fode, fbc, u0, tspan[1], p; kwargs...))
+end
+
+function check_compatible_system(T::Union{Type{BVPFunction}, Type{BVProblem}}, sys::System)
+    check_time_dependent(sys, T)
+    check_not_dde(sys)
+    check_no_cost(sys, T)
+    check_has_constraints(sys, T)
+    check_no_jumps(sys, T)
+    check_no_noise(sys, T)
+    check_is_continuous(sys, T)
+end

src/problems/compatibility.jl

@@ -69,6 +69,15 @@ function check_no_constraints(sys::System, T)
     end
 end
 
+function check_has_constraints(sys::System, T)
+    if isempty(constraints(sys))
+        throw(SystemCompatibilityError("""
+        A system without constraints cannot be used to construct a `$T`. Consider an \
+        `ODEProblem` instead.
+        """))
+    end
+end
+
 function check_no_jumps(sys::System, T)
     if !isempty(jumps(sys))
         throw(SystemCompatibilityError("""

src/systems/codegen.jl

@@ -317,3 +317,55 @@ function isautonomous(sys::System)
     tgrad = calculate_tgrad(sys; simplify = true)
     all(iszero, tgrad)
 end
+
+function get_bv_solution_symbol(ns)
+    only(@variables BV_SOLUTION(..)[1:ns])
+end
+
+function get_constraint_unknown_subs!(subs::Dict, cons::Vector, stidxmap::Dict, iv, sol)
+    vs = vars(cons)
+    for v in vs
+        iscall(v) || continue
+        op = operation(v)
+        args = arguments(v)
+        issym(op) && length(args) == 1 || continue
+        newv = op(iv)
+        haskey(stidxmap, newv) || continue
+        subs[v] = sol(args[1])[stidxmap[newv]]
+    end
+end
+
+function generate_boundary_conditions(sys::System, u0, u0_idxs, t0; expression = Val{true},
+        eval_expression = false, eval_module = @__MODULE__, kwargs...)
+    iv = get_iv(sys)
+    sts = unknowns(sys)
+    ps = parameters(sys)
+    np = length(ps)
+    ns = length(sts)
+    stidxmap = Dict([v => i for (i, v) in enumerate(sts)])
+    pidxmap = Dict([v => i for (i, v) in enumerate(ps)])
+
+    sol = get_bv_solution_symbol(ns)
+
+    cons = [con.lhs - con.rhs for con in constraints(sys)]
+    conssubs = Dict()
+    get_constraint_unknown_subs!(conssubs, cons, stidxmap, iv, sol)
+    cons = map(x -> fast_substitute(x, conssubs), cons)
+
+    init_conds = Any[]
+    for i in u0_idxs
+        expr = sol(t0)[i] - u0[i]
+        push!(init_conds, expr)
+    end
+
+    exprs = vcat(init_conds, cons)
+    _p = reorder_parameters(sys, ps)
+
+    res = build_function_wrapper(sys, exprs, sol, _p..., iv; output_type = Array, kwargs...)
+    if expression == Val{true}
+        return res
+    end
+
+    f_oop, f_iip = eval_or_rgf.(res; eval_expression, eval_module)
+    return GeneratedFunctionWrapper{(2, 3, is_split(sys))}(f_oop, f_iip)
+end