diff --git a/Project.toml b/Project.toml index 19063a0342..8022fff78a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.54.0" +version = "9.58.0" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" @@ -44,6 +44,7 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" +SCCNonlinearSolve = "9dfe8606-65a1-4bb3-9748-cb89d1561431" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" @@ -121,13 +122,15 @@ NonlinearSolve = "3.14, 4" OffsetArrays = "1" OrderedCollections = "1" OrdinaryDiffEq = "6.82.0" -OrdinaryDiffEqCore = "1.7.0" +OrdinaryDiffEqCore = "1.13.0" +OrdinaryDiffEqNonlinearSolve = "1.3.0" PrecompileTools = "1" REPL = "1" RecursiveArrayTools = "3.26" Reexport = "0.2, 1" RuntimeGeneratedFunctions = "0.5.9" -SciMLBase = "2.57.1" +SCCNonlinearSolve = "1.0.0" +SciMLBase = "2.66" SciMLStructures = "1.0" Serialization = "1" Setfield = "0.7, 0.8, 1" @@ -162,6 +165,7 @@ OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1" OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" +OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -176,4 +180,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEq", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET"] +test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEq", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve"] diff --git a/docs/Project.toml b/docs/Project.toml index 8d0f1667c9..a358455503 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -15,6 +15,7 @@ OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226" +Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" @@ -37,6 +38,7 @@ OptimizationOptimJL = "0.1, 0.4" OrdinaryDiffEq = "6.31" Plots = "1.36" SciMLStructures = "1.1" +Setfield = "1" StochasticDiffEq = "6" SymbolicIndexingInterface = "0.3.1" SymbolicUtils = "3" diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index f425fdce5b..23e1e6d7d1 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -378,3 +378,218 @@ sol.ps[c] # sol[c] will error, since `c` is not a timeseries value ``` It can be seen that the timeseries for `c` is not saved. + +## [(Experimental) Imperative affects](@id imp_affects) + +The `ImperativeAffect` can be used as an alternative to the aforementioned functional affect form. Note +that `ImperativeAffect` is still experimental; to emphasize this, we do not export it and it should be +included as `ModelingToolkit.ImperativeAffect`. `ImperativeAffect` aims to simplify the manipulation of +system state. + +We will use two examples to describe `ImperativeAffect`: a simple heater and a quadrature encoder. +These examples will also demonstrate advanced usage of `ModelingToolkit.SymbolicContinuousCallback`, +the low-level interface of the tuple form converts into that allows control over the SciMLBase-level +event that is generated for a continuous event. + +### [Heater](@id heater_events) + +Bang-bang control of a heater connected to a leaky plant requires hysteresis in order to prevent rapid control oscillation. + +```@example events +@variables temp(t) +params = @parameters furnace_on_threshold=0.5 furnace_off_threshold=0.7 furnace_power=1.0 leakage=0.1 furnace_on(t)::Bool=false +eqs = [ + D(temp) ~ furnace_on * furnace_power - temp^2 * leakage +] +``` + +Our plant is simple. We have a heater that's turned on and off by the time-indexed parameter `furnace_on` +which adds `furnace_power` forcing to the system when enabled. We then leak heat proportional to `leakage` +as a function of the square of the current temperature. + +We need a controller with hysteresis to control the plant. We wish the furnace to turn on when the temperature +is below `furnace_on_threshold` and off when above `furnace_off_threshold`, while maintaining its current state +in between. To do this, we create two continuous callbacks: + +```@example events +using Setfield +furnace_disable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i + @set! x.furnace_on = false + end) +furnace_enable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_on_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i + @set! x.furnace_on = true + end) +``` + +We're using the explicit form of `SymbolicContinuousCallback` here, though +so far we aren't using anything that's not possible with the implicit interface. +You can also write + +```julia +[temp ~ furnace_off_threshold] => ModelingToolkit.ImperativeAffect(modified = (; + furnace_on)) do x, o, i, c + @set! x.furnace_on = false +end +``` + +and it would work the same. + +The `ImperativeAffect` is the larger change in this example. `ImperativeAffect` has the constructor signature + +```julia +ImperativeAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx) +``` + +that accepts the function to call, a named tuple of both the names of and symbolic values representing +values in the system to be modified, a named tuple of the values that are merely observed (that is, used from +the system but not modified), and a context that's passed to the affect function. + +In our example, each event merely changes whether the furnace is on or off. Accordingly, we pass a `modified` tuple +`(; furnace_on)` (creating a `NamedTuple` equivalent to `(furnace_on = furnace_on)`). `ImperativeAffect` will then +evaluate this before calling our function to fill out all of the numerical values, then apply them back to the system +once our affect function returns. Furthermore, it will check that it is possible to do this assignment. + +The function given to `ImperativeAffect` needs to have the signature: + +```julia +f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple +``` + +The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions. +In our example, if `furnace_on` is `false`, then the value of the `x` that's passed in as `modified` will be `(furnace_on = false)`. +The modified values should be passed out in the same format: to set `furnace_on` to `true` we need to return a tuple `(furnace_on = true)`. +The examples does this with Setfield, recreating the result tuple before returning it; the returned tuple may optionally be missing values as +well, in which case those values will not be written back to the problem. + +Accordingly, we can now interpret the `ImperativeAffect` definitions to mean that when `temp = furnace_off_threshold` we +will write `furnace_on = false` back to the system, and when `temp = furnace_on_threshold` we will write `furnace_on = true` back +to the system. + +```@example events +@named sys = ODESystem( + eqs, t, [temp], params; continuous_events = [furnace_disable, furnace_enable]) +ss = structural_simplify(sys) +prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 10.0)) +sol = solve(prob, Tsit5()) +plot(sol) +hline!([sol.ps[furnace_off_threshold], sol.ps[furnace_on_threshold]], + l = (:black, 1), primary = false) +``` + +Here we see exactly the desired hysteresis. The heater starts on until the temperature hits +`furnace_off_threshold`. The temperature then bleeds away until `furnace_on_threshold` at which +point the furnace turns on again until `furnace_off_threshold` and so on and so forth. The controller +is effectively regulating the temperature of the plant. + +### [Quadrature Encoder](@id quadrature) + +For a more complex application we'll look at modeling a quadrature encoder attached to a shaft spinning at a constant speed. +Traditionally, a quadrature encoder is built out of a code wheel that interrupts the sensors at constant intervals and two sensors slightly out of phase with one another. +A state machine can take the pattern of pulses produced by the two sensors and determine the number of steps that the shaft has spun. The state machine takes the new value +from each sensor and the old values and decodes them into the direction that the wheel has spun in this step. + +```@example events +@variables theta(t) omega(t) +params = @parameters qA=0 qB=0 hA=0 hB=0 cnt::Int=0 +eqs = [D(theta) ~ omega + omega ~ 1.0] +``` + +Our continuous-time system is extremely simple. We have two unknown variables `theta` for the angle of the shaft +and `omega` for the rate at which it's spinning. We then have parameters for the state machine `qA, qB, hA, hB` +(corresponding to the current quadrature of the A/B sensors and the historical ones) and a step count `cnt`. + +We'll then implement the decoder as a simple Julia function. + +```@example events +function decoder(oldA, oldB, newA, newB) + state = (oldA, oldB, newA, newB) + if state == (0, 0, 1, 0) || state == (1, 0, 1, 1) || state == (1, 1, 0, 1) || + state == (0, 1, 0, 0) + return 1 + elseif state == (0, 0, 0, 1) || state == (0, 1, 1, 1) || state == (1, 1, 1, 0) || + state == (1, 0, 0, 0) + return -1 + elseif state == (0, 0, 0, 0) || state == (0, 1, 0, 1) || state == (1, 0, 1, 0) || + state == (1, 1, 1, 1) + return 0 + else + return 0 # err is interpreted as no movement + end +end +``` + +Based on the current and old state, this function will return 1 if the wheel spun in the positive direction, +-1 if in the negative, and 0 otherwise. + +The encoder state advances when the occlusion begins or ends. We model the +code wheel as simply detecting when `cos(100*theta)` is 0; if we're at a positive +edge of the 0 crossing, then we interpret that as occlusion (so the discrete `qA` goes to 1). Otherwise, if `cos` is +going negative, we interpret that as lack of occlusion (so the discrete goes to 0). The decoder function is +then invoked to update the count with this new information. + +We can implement this in one of two ways: using edge sign detection or right root finding. For exposition, we +will implement each sensor differently. + +For sensor A, we're using the edge detection method. By providing a different affect to `SymbolicContinuousCallback`'s +`affect_neg` argument, we can specify different behaviour for the negative crossing vs. the positive crossing of the root. +In our encoder, we interpret this as occlusion or nonocclusion of the sensor, update the internal state, and tick the decoder. + +```@example events +qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], + ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, c, i + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 1 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x + end, + affect_neg = ModelingToolkit.ImperativeAffect( + (; qA, hA, hB, cnt), (; qB)) do x, o, c, i + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 0 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x + end) +``` + +The other way we can implement a sensor is by changing the root find. +Normally, we use left root finding; the affect will be invoked instantaneously _before_ +the root is crossed. This makes it trickier to figure out what the new state is. +Instead, we can use right root finding: + +```@example events +qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0], + ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA, theta)) do x, o, c, i + @set! x.hA = o.qA + @set! x.hB = x.qB + @set! x.qB = clamp(sign(cos(100 * o.theta - π / 2)), 0.0, 1.0) + @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) + x + end; rootfind = SciMLBase.RightRootFind) +``` + +Here, sensor B is located `π / 2` behind sensor A in angular space, so we're adjusting our +trigger function accordingly. We here ask for right root finding on the callback, so we know +that the value of said function will have the "new" sign rather than the old one. Thus, we can +determine the new state of the sensor from the sign of the indicator function evaluated at the +affect activation point, with -1 mapped to 0. + +We can now simulate the encoder. + +```@example events +@named sys = ODESystem( + eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) +ss = structural_simplify(sys) +prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) +sol = solve(prob, Tsit5(); dtmax = 0.01) +sol.ps[cnt] +``` + +`cos(100*theta)` will have 200 crossings in the half rotation we've gone through, so the encoder would notionally count 200 steps. +Our encoder counts 198 steps (it loses one step to initialization and one step due to the final state falling squarely on an edge). diff --git a/docs/src/examples/higher_order.md b/docs/src/examples/higher_order.md index 7dafe758dc..fac707525f 100644 --- a/docs/src/examples/higher_order.md +++ b/docs/src/examples/higher_order.md @@ -3,7 +3,7 @@ ModelingToolkit has a system for transformations of mathematical systems. These transformations allow for symbolically changing the representation of the model to problems that are easier to -numerically solve. One simple to demonstrate transformation is the +numerically solve. One simple to demonstrate transformation, is `structural_simplify`, which does a lot of tricks, one being the transformation that turns an Nth order ODE into N coupled 1st order ODEs. @@ -15,16 +15,28 @@ We utilize the derivative operator twice here to define the second order: using ModelingToolkit, OrdinaryDiffEq using ModelingToolkit: t_nounits as t, D_nounits as D -@parameters σ ρ β -@variables x(t) y(t) z(t) - -eqs = [D(D(x)) ~ σ * (y - x), - D(y) ~ x * (ρ - z) - y, - D(z) ~ x * y - β * z] - -@named sys = ODESystem(eqs, t) +@mtkmodel SECOND_ORDER begin + @parameters begin + σ = 28.0 + ρ = 10.0 + β = 8 / 3 + end + @variables begin + x(t) = 1.0 + y(t) = 0.0 + z(t) = 0.0 + end + @equations begin + D(D(x)) ~ σ * (y - x) + D(y) ~ x * (ρ - z) - y + D(z) ~ x * y - β * z + end +end +@mtkbuild sys = SECOND_ORDER() ``` +The second order ODE has been automatically transformed to two first order ODEs. + Note that we could've used an alternative syntax for 2nd order, i.e. `D = Differential(t)^2` and then `D(x)` would be the second derivative, and this syntax extends to `N`-th order. Also, we can use `*` or `∘` to compose @@ -33,28 +45,17 @@ and this syntax extends to `N`-th order. Also, we can use `*` or `∘` to compos Now let's transform this into the `ODESystem` of first order components. We do this by calling `structural_simplify`: -```@example orderlowering -sys = structural_simplify(sys) -``` - Now we can directly numerically solve the lowered system. Note that, following the original problem, the solution requires knowing the -initial condition for `x'`, and thus we include that in our input -specification: +initial condition for both `x` and `D(x)`. +The former already got assigned a default value in the `@mtkmodel`, +but we still have to provide a value for the latter. ```@example orderlowering -u0 = [D(x) => 2.0, - x => 1.0, - y => 0.0, - z => 0.0] - -p = [σ => 28.0, - ρ => 10.0, - β => 8 / 3] - +u0 = [D(sys.x) => 2.0] tspan = (0.0, 100.0) -prob = ODEProblem(sys, u0, tspan, p, jac = true) +prob = ODEProblem(sys, u0, tspan, [], jac = true) sol = solve(prob, Tsit5()) -using Plots; -plot(sol, idxs = (x, y)); +using Plots +plot(sol, idxs = (sys.x, sys.y)) ``` diff --git a/docs/src/examples/modelingtoolkitize_index_reduction.md b/docs/src/examples/modelingtoolkitize_index_reduction.md index 415d5b85ff..b19ea46701 100644 --- a/docs/src/examples/modelingtoolkitize_index_reduction.md +++ b/docs/src/examples/modelingtoolkitize_index_reduction.md @@ -51,6 +51,14 @@ In this tutorial, we will look at the pendulum system: \end{aligned} ``` +These equations can be derived using the [Lagrangian equation of the first kind.](https://en.wikipedia.org/wiki/Lagrangian_mechanics#Lagrangian) +Specifically, for a pendulum with unit mass and length $L$, which thus has +kinetic energy $\frac{1}{2}(v_x^2 + v_y^2)$, +potential energy $gy$, +and holonomic constraint $x^2 + y^2 - L^2 = 0$. +The Lagrange multiplier related to this constraint is equal to half of $T$, +and represents the tension in the rope of the pendulum. + As a good DifferentialEquations.jl user, one would follow [the mass matrix DAE tutorial](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/dae_example/#Mass-Matrix-Differential-Algebraic-Equations-(DAEs)) to arrive at code for simulating the model: diff --git a/ext/MTKHomotopyContinuationExt.jl b/ext/MTKHomotopyContinuationExt.jl index fa5d1bfcd4..c4a090d9a8 100644 --- a/ext/MTKHomotopyContinuationExt.jl +++ b/ext/MTKHomotopyContinuationExt.jl @@ -11,217 +11,6 @@ using ModelingToolkit: iscomplete, parameters, has_index_cache, get_index_cache, const MTK = ModelingToolkit -function contains_variable(x, wrt) - any(y -> occursin(y, x), wrt) -end - -""" -Possible reasons why a term is not polynomial -""" -MTK.EnumX.@enumx NonPolynomialReason begin - NonIntegerExponent - ExponentContainsUnknowns - BaseNotPolynomial - UnrecognizedOperation -end - -function display_reason(reason::NonPolynomialReason.T, sym) - if reason == NonPolynomialReason.NonIntegerExponent - pow = arguments(sym)[2] - "In $sym: Exponent $pow is not an integer" - elseif reason == NonPolynomialReason.ExponentContainsUnknowns - pow = arguments(sym)[2] - "In $sym: Exponent $pow contains unknowns of the system" - elseif reason == NonPolynomialReason.BaseNotPolynomial - base = arguments(sym)[1] - "In $sym: Base $base is not a polynomial in the unknowns" - elseif reason == NonPolynomialReason.UnrecognizedOperation - op = operation(sym) - """ - In $sym: Operation $op is not recognized. Allowed polynomial operations are \ - `*, /, +, -, ^`. - """ - else - error("This should never happen. Please open an issue in ModelingToolkit.jl.") - end -end - -mutable struct PolynomialData - non_polynomial_terms::Vector{BasicSymbolic} - reasons::Vector{NonPolynomialReason.T} - has_parametric_exponent::Bool -end - -PolynomialData() = PolynomialData(BasicSymbolic[], NonPolynomialReason.T[], false) - -abstract type PolynomialTransformationError <: Exception end - -struct MultivarTerm <: PolynomialTransformationError - term::Any - vars::Any -end - -function Base.showerror(io::IO, err::MultivarTerm) - println(io, - "Cannot convert system to polynomial: Found term $(err.term) which is a function of multiple unknowns $(err.vars).") -end - -struct MultipleTermsOfSameVar <: PolynomialTransformationError - terms::Any - var::Any -end - -function Base.showerror(io::IO, err::MultipleTermsOfSameVar) - println(io, - "Cannot convert system to polynomial: Found multiple non-polynomial terms $(err.terms) involving the same unknown $(err.var).") -end - -struct SymbolicSolveFailure <: PolynomialTransformationError - term::Any - var::Any -end - -function Base.showerror(io::IO, err::SymbolicSolveFailure) - println(io, - "Cannot convert system to polynomial: Unable to symbolically solve $(err.term) for $(err.var).") -end - -struct NemoNotLoaded <: PolynomialTransformationError end - -function Base.showerror(io::IO, err::NemoNotLoaded) - println(io, - "ModelingToolkit may be able to solve this system as a polynomial system if `Nemo` is loaded. Run `import Nemo` and try again.") -end - -struct VariablesAsPolyAndNonPoly <: PolynomialTransformationError - vars::Any -end - -function Base.showerror(io::IO, err::VariablesAsPolyAndNonPoly) - println(io, - "Cannot convert convert system to polynomial: Variables $(err.vars) occur in both polynomial and non-polynomial terms in the system.") -end - -struct NotPolynomialError <: Exception - transformation_err::Union{PolynomialTransformationError, Nothing} - eq::Vector{Equation} - data::Vector{PolynomialData} -end - -function Base.showerror(io::IO, err::NotPolynomialError) - if err.transformation_err !== nothing - Base.showerror(io, err.transformation_err) - end - for (eq, data) in zip(err.eq, err.data) - if isempty(data.non_polynomial_terms) - continue - end - println(io, - "Equation $(eq) is not a polynomial in the unknowns for the following reasons:") - for (term, reason) in zip(data.non_polynomial_terms, data.reasons) - println(io, display_reason(reason, term)) - end - end -end - -function is_polynomial!(data, y, wrt) - process_polynomial!(data, y, wrt) - isempty(data.reasons) -end - -""" -$(TYPEDSIGNATURES) - -Return information about the polynmial `x` with respect to variables in `wrt`, -writing said information to `data`. -""" -function process_polynomial!(data::PolynomialData, x, wrt) - x = unwrap(x) - symbolic_type(x) == NotSymbolic() && return true - iscall(x) || return true - contains_variable(x, wrt) || return true - any(isequal(x), wrt) && return true - - if operation(x) in (*, +, -, /) - # `map` because `all` will early exit, but we want to search - # through everything to get all the non-polynomial terms - return all(map(y -> is_polynomial!(data, y, wrt), arguments(x))) - end - if operation(x) == (^) - b, p = arguments(x) - is_pow_integer = symtype(p) <: Integer - if !is_pow_integer - push!(data.non_polynomial_terms, x) - push!(data.reasons, NonPolynomialReason.NonIntegerExponent) - end - if symbolic_type(p) != NotSymbolic() - data.has_parametric_exponent = true - end - - exponent_has_unknowns = contains_variable(p, wrt) - if exponent_has_unknowns - push!(data.non_polynomial_terms, x) - push!(data.reasons, NonPolynomialReason.ExponentContainsUnknowns) - end - base_polynomial = is_polynomial!(data, b, wrt) - return base_polynomial && !exponent_has_unknowns && is_pow_integer - end - push!(data.non_polynomial_terms, x) - push!(data.reasons, NonPolynomialReason.UnrecognizedOperation) - return false -end - -""" -$(TYPEDSIGNATURES) - -Given a `x`, a polynomial in variables in `wrt` which may contain rational functions, -express `x` as a single rational function with polynomial `num` and denominator `den`. -Return `(num, den)`. -""" -function handle_rational_polynomials(x, wrt) - x = unwrap(x) - symbolic_type(x) == NotSymbolic() && return x, 1 - iscall(x) || return x, 1 - contains_variable(x, wrt) || return x, 1 - any(isequal(x), wrt) && return x, 1 - - # simplify_fractions cancels out some common factors - # and expands (a / b)^c to a^c / b^c, so we only need - # to handle these cases - x = simplify_fractions(x) - op = operation(x) - args = arguments(x) - - if op == / - # numerator and denominator are trivial - num, den = args - # but also search for rational functions in numerator - n, d = handle_rational_polynomials(num, wrt) - num, den = n, den * d - elseif op == + - num = 0 - den = 1 - - # we don't need to do common denominator - # because we don't care about cases where denominator - # is zero. The expression is zero when all the numerators - # are zero. - for arg in args - n, d = handle_rational_polynomials(arg, wrt) - num += n - den *= d - end - else - return x, 1 - end - # if the denominator isn't a polynomial in `wrt`, better to not include it - # to reduce the size of the gcd polynomial - if !contains_variable(den, wrt) - return num / den, 1 - end - return num, den -end - """ $(TYPEDSIGNATURES) @@ -289,12 +78,6 @@ end SymbolicIndexingInterface.parameter_values(s::MTKHomotopySystem) = s.p -struct PolynomialTransformationData - new_var::BasicSymbolic - term::BasicSymbolic - inv_term::Vector -end - """ $(TYPEDSIGNATURES) @@ -312,128 +95,37 @@ Keyword arguments: All other keyword arguments are forwarded to `HomotopyContinuation.solver_startsystems`. """ function MTK.HomotopyContinuationProblem( - sys::NonlinearSystem, u0map, parammap = nothing; eval_expression = false, - eval_module = ModelingToolkit, warn_parametric_exponent = true, kwargs...) + sys::NonlinearSystem, u0map, parammap = nothing; kwargs...) + prob = MTK._safe_HomotopyContinuationProblem(sys, u0map, parammap; kwargs...) + prob isa MTK.HomotopyContinuationProblem || throw(prob) + return prob +end + +function MTK._safe_HomotopyContinuationProblem(sys, u0map, parammap = nothing; kwargs...) if !iscomplete(sys) error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `HomotopyContinuationProblem`") end - - dvs = unknowns(sys) - # we need to consider `full_equations` because observed also should be - # polynomials (if used in equations) and we don't know if observed is used - # in denominator. - # This is not the most efficient, and would be improved significantly with - # CSE/hashconsing. - eqs = full_equations(sys) - - polydata = map(eqs) do eq - data = PolynomialData() - process_polynomial!(data, eq.lhs, dvs) - process_polynomial!(data, eq.rhs, dvs) - data + transformation = MTK.PolynomialTransformation(sys) + if transformation isa MTK.NotPolynomialError + return transformation end - - has_parametric_exponents = any(d -> d.has_parametric_exponent, polydata) - - all_non_poly_terms = mapreduce(d -> d.non_polynomial_terms, vcat, polydata) - unique!(all_non_poly_terms) - - var_to_nonpoly = Dict{BasicSymbolic, PolynomialTransformationData}() - - is_poly = true - transformation_err = nothing - for t in all_non_poly_terms - # if the term involves multiple unknowns, we can't invert it - dvs_in_term = map(x -> occursin(x, t), dvs) - if count(dvs_in_term) > 1 - transformation_err = MultivarTerm(t, dvs[dvs_in_term]) - is_poly = false - break - end - # we already have a substitution solving for `var` - var = dvs[findfirst(dvs_in_term)] - if haskey(var_to_nonpoly, var) && !isequal(var_to_nonpoly[var].term, t) - transformation_err = MultipleTermsOfSameVar([t, var_to_nonpoly[var].term], var) - is_poly = false - break - end - # we want to solve `term - new_var` for `var` - new_var = gensym(Symbol(var)) - new_var = unwrap(only(@variables $new_var)) - invterm = Symbolics.ia_solve( - t - new_var, var; complex_roots = false, periodic_roots = false, warns = false) - # if we can't invert it, quit - if invterm === nothing || isempty(invterm) - transformation_err = SymbolicSolveFailure(t, var) - is_poly = false - break - end - # `ia_solve` returns lazy terms i.e. `asin(1.0)` instead of `pi/2` - # this just evaluates the constant expressions - invterm = Symbolics.substitute.(invterm, (Dict(),)) - # RootsOf implies Symbolics couldn't solve the inner polynomial because - # `Nemo` wasn't loaded. - if any(x -> MTK.iscall(x) && MTK.operation(x) == Symbolics.RootsOf, invterm) - transformation_err = NemoNotLoaded() - is_poly = false - break - end - var_to_nonpoly[var] = PolynomialTransformationData(new_var, t, invterm) - end - - if !is_poly - throw(NotPolynomialError(transformation_err, eqs, polydata)) - end - - subrules = Dict() - combinations = Vector[] - new_dvs = [] - for x in dvs - if haskey(var_to_nonpoly, x) - _data = var_to_nonpoly[x] - subrules[_data.term] = _data.new_var - push!(combinations, _data.inv_term) - push!(new_dvs, _data.new_var) - else - push!(combinations, [x]) - push!(new_dvs, x) - end - end - all_solutions = collect.(collect(Iterators.product(combinations...))) - - denoms = [] - eqs2 = map(eqs) do eq - t = eq.rhs - eq.lhs - t = Symbolics.fixpoint_sub(t, subrules; maxiters = length(dvs)) - # the substituted variable occurs outside the substituted term - poly_and_nonpoly = map(dvs) do x - haskey(var_to_nonpoly, x) && occursin(x, t) - end - if any(poly_and_nonpoly) - throw(NotPolynomialError( - VariablesAsPolyAndNonPoly(dvs[poly_and_nonpoly]), eqs, polydata)) - end - - num, den = handle_rational_polynomials(t, new_dvs) - # make factors different elements, otherwise the nonzero factors artificially - # inflate the error of the zero factor. - if iscall(den) && operation(den) == * - for arg in arguments(den) - # ignore constant factors - symbolic_type(arg) == NotSymbolic() && continue - push!(denoms, abs(arg)) - end - elseif symbolic_type(den) != NotSymbolic() - push!(denoms, abs(den)) - end - return 0 ~ num + result = MTK.transform_system(sys, transformation) + if result isa MTK.NotPolynomialError + return result end + MTK.HomotopyContinuationProblem(sys, transformation, result, u0map, parammap; kwargs...) +end - sys2 = MTK.@set sys.eqs = eqs2 - MTK.@set! sys2.unknowns = new_dvs - # remove observed equations to avoid adding them in codegen - MTK.@set! sys2.observed = Equation[] - MTK.@set! sys2.substitutions = nothing +function MTK.HomotopyContinuationProblem( + sys::MTK.NonlinearSystem, transformation::MTK.PolynomialTransformation, + result::MTK.PolynomialTransformationResult, u0map, + parammap = nothing; eval_expression = false, + eval_module = ModelingToolkit, warn_parametric_exponent = true, kwargs...) + sys2 = result.sys + denoms = result.denominators + polydata = transformation.polydata + new_dvs = transformation.new_dvs + all_solutions = transformation.all_solutions _, u0, p = MTK.process_SciMLProblem( MTK.EmptySciMLFunction, sys, u0map, parammap; eval_expression, eval_module) @@ -443,10 +135,11 @@ function MTK.HomotopyContinuationProblem( unpack_solution = MTK.build_explicit_observed_function(sys2, all_solutions) hvars = symbolics_to_hc.(new_dvs) - mtkhsys = MTKHomotopySystem(nlfn.f, p, nlfn.jac, hvars, length(eqs)) + mtkhsys = MTKHomotopySystem(nlfn.f, p, nlfn.jac, hvars, length(new_dvs)) obsfn = MTK.ObservedFunctionCache(sys; eval_expression, eval_module) + has_parametric_exponents = any(d -> d.has_parametric_exponent, polydata) if has_parametric_exponents if warn_parametric_exponent @warn """ diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index de8e69c41f..20b2ada8fa 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -50,6 +50,7 @@ using Distributed import JuliaFormatter using MLStyle using NonlinearSolve +import SCCNonlinearSolve using Reexport using RecursiveArrayTools import Graphs: SimpleDiGraph, add_edge!, incidence_matrix @@ -145,10 +146,12 @@ include("systems/parameter_buffer.jl") include("systems/abstractsystem.jl") include("systems/model_parsing.jl") include("systems/connectors.jl") +include("systems/imperative_affect.jl") include("systems/callbacks.jl") include("systems/problem_utils.jl") include("systems/nonlinear/nonlinearsystem.jl") +include("systems/nonlinear/homotopy_continuation.jl") include("systems/diffeqs/odesystem.jl") include("systems/diffeqs/sdesystem.jl") include("systems/diffeqs/abstractodesystem.jl") diff --git a/src/structural_transformation/bipartite_tearing/modia_tearing.jl b/src/structural_transformation/bipartite_tearing/modia_tearing.jl index cef2f5f6d7..5da873afdf 100644 --- a/src/structural_transformation/bipartite_tearing/modia_tearing.jl +++ b/src/structural_transformation/bipartite_tearing/modia_tearing.jl @@ -62,6 +62,15 @@ function tear_graph_block_modia!(var_eq_matching, ict, solvable_graph, eqs, vars return nothing end +function build_var_eq_matching(structure::SystemStructure, ::Type{U} = Unassigned; + varfilter::F2 = v -> true, eqfilter::F3 = eq -> true) where {U, F2, F3} + @unpack graph, solvable_graph = structure + var_eq_matching = maximal_matching(graph, eqfilter, varfilter, U) + matching_len = max(length(var_eq_matching), + maximum(x -> x isa Int ? x : 0, var_eq_matching, init = 0)) + return complete(var_eq_matching, matching_len), matching_len +end + function tear_graph_modia(structure::SystemStructure, isder::F = nothing, ::Type{U} = Unassigned; varfilter::F2 = v -> true, @@ -78,10 +87,7 @@ function tear_graph_modia(structure::SystemStructure, isder::F = nothing, # find them here [TODO: It would be good to have an explicit example of this.] @unpack graph, solvable_graph = structure - var_eq_matching = maximal_matching(graph, eqfilter, varfilter, U) - matching_len = max(length(var_eq_matching), - maximum(x -> x isa Int ? x : 0, var_eq_matching, init = 0)) - var_eq_matching = complete(var_eq_matching, matching_len) + var_eq_matching, matching_len = build_var_eq_matching(structure, U; varfilter, eqfilter) full_var_eq_matching = copy(var_eq_matching) var_sccs = find_var_sccs(graph, var_eq_matching) vargraph = DiCMOBiGraph{true}(graph, 0, Matching(matching_len)) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 3a9bbd33e1..e44f250a7f 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -162,11 +162,12 @@ object. """ function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys), ps = parameters(sys); wrap_code = nothing, postprocess_fbody = nothing, states = nothing, - expression = Val{true}, eval_expression = false, eval_module = @__MODULE__, kwargs...) + expression = Val{true}, eval_expression = false, eval_module = @__MODULE__, + cachesyms::Tuple = (), kwargs...) if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system.") end - p = reorder_parameters(sys, unwrap.(ps)) + p = (reorder_parameters(sys, unwrap.(ps))..., cachesyms...) isscalar = !(exprs isa AbstractArray) if wrap_code === nothing wrap_code = isscalar ? identity : (identity, identity) @@ -187,7 +188,7 @@ function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys postprocess_fbody, states, wrap_code = wrap_code .∘ wrap_mtkparameters(sys, isscalar) .∘ - wrap_array_vars(sys, exprs; dvs) .∘ + wrap_array_vars(sys, exprs; dvs, cachesyms) .∘ wrap_parameter_dependencies(sys, isscalar), expression = Val{true} ) @@ -199,7 +200,7 @@ function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys postprocess_fbody, states, wrap_code = wrap_code .∘ wrap_mtkparameters(sys, isscalar) .∘ - wrap_array_vars(sys, exprs; dvs) .∘ + wrap_array_vars(sys, exprs; dvs, cachesyms) .∘ wrap_parameter_dependencies(sys, isscalar), expression = Val{true} ) @@ -231,117 +232,47 @@ end function wrap_array_vars( sys::AbstractSystem, exprs; dvs = unknowns(sys), ps = parameters(sys), - inputs = nothing, history = false) + inputs = nothing, history = false, cachesyms::Tuple = ()) isscalar = !(exprs isa AbstractArray) - array_vars = Dict{Any, AbstractArray{Int}}() - if dvs !== nothing - for (j, x) in enumerate(dvs) - if iscall(x) && operation(x) == getindex - arg = arguments(x)[1] - inds = get!(() -> Int[], array_vars, arg) - push!(inds, j) - end - end - for (k, inds) in array_vars - if inds == (inds′ = inds[1]:inds[end]) - array_vars[k] = inds′ - end - end + var_to_arridxs = Dict() - uind = 1 - else + if dvs === nothing uind = 0 - end - # values are (indexes, index of buffer, size of parameter) - array_parameters = Dict{Any, Tuple{AbstractArray{Int}, Int, Tuple{Vararg{Int}}}}() - # If for some reason different elements of an array parameter are in different buffers - other_array_parameters = Dict{Any, Any}() - - hasinputs = inputs !== nothing - input_vars = Dict{Any, AbstractArray{Int}}() - if hasinputs - for (j, x) in enumerate(inputs) - if iscall(x) && operation(x) == getindex - arg = arguments(x)[1] - inds = get!(() -> Int[], input_vars, arg) - push!(inds, j) - end - end - for (k, inds) in input_vars - if inds == (inds′ = inds[1]:inds[end]) - input_vars[k] = inds′ - end - end - end - if has_index_cache(sys) - ic = get_index_cache(sys) else - ic = nothing - end - if ps isa Tuple && eltype(ps) <: AbstractArray - ps = Iterators.flatten(ps) - end - for p in ps - p = unwrap(p) - if iscall(p) && operation(p) == getindex - p = arguments(p)[1] - end - symtype(p) <: AbstractArray && Symbolics.shape(p) != Symbolics.Unknown() || continue - scal = collect(p) - # all scalarized variables are in `ps` - any(isequal(p), ps) || all(x -> any(isequal(x), ps), scal) || continue - (haskey(array_parameters, p) || haskey(other_array_parameters, p)) && continue - - idx = parameter_index(sys, p) - idx isa Int && continue - if idx isa ParameterIndex - if idx.portion != SciMLStructures.Tunable() - continue - end - array_parameters[p] = (vec(idx.idx), 1, size(idx.idx)) + uind = 1 + for (i, x) in enumerate(dvs) + iscall(x) && operation(x) == getindex || continue + arg = arguments(x)[1] + inds = get!(() -> [], var_to_arridxs, arg) + push!(inds, (uind, i)) + end + end + p_start = uind + 1 + history + rps = (reorder_parameters(sys, ps)..., cachesyms...) + if inputs !== nothing + rps = (inputs, rps...) + end + for sym in reduce(vcat, rps; init = []) + iscall(sym) && operation(sym) == getindex || continue + arg = arguments(sym)[1] + + bufferidx = findfirst(buf -> any(isequal(sym), buf), rps) + idxinbuffer = findfirst(isequal(sym), rps[bufferidx]) + inds = get!(() -> [], var_to_arridxs, arg) + push!(inds, (p_start + bufferidx - 1, idxinbuffer)) + end + + viewsyms = Dict() + splitsyms = Dict() + for (arrsym, idxs) in var_to_arridxs + length(idxs) == length(arrsym) || continue + # allequal(first, idxs) is a 1.11 feature + if allequal(Iterators.map(first, idxs)) + viewsyms[arrsym] = (first(first(idxs)), reshape(last.(idxs), size(arrsym))) else - # idx === nothing - idxs = map(Base.Fix1(parameter_index, sys), scal) - if first(idxs) isa ParameterIndex - buffer_idxs = map(Base.Fix1(iterated_buffer_index, ic), idxs) - if allequal(buffer_idxs) - buffer_idx = first(buffer_idxs) - if first(idxs).portion == SciMLStructures.Tunable() - idxs = map(x -> x.idx, idxs) - else - idxs = map(x -> x.idx[end], idxs) - end - else - other_array_parameters[p] = scal - continue - end - else - buffer_idx = 1 - end - - sz = size(idxs) - if vec(idxs) == idxs[begin]:idxs[end] - idxs = idxs[begin]:idxs[end] - elseif vec(idxs) == idxs[begin]:-1:idxs[end] - idxs = idxs[begin]:-1:idxs[end] - end - idxs = vec(idxs) - array_parameters[p] = (idxs, buffer_idx, sz) + splitsyms[arrsym] = reshape(idxs, size(arrsym)) end end - - inputind = if history - uind + 2 - else - uind + 1 - end - params_offset = if history && hasinputs - uind + 2 - elseif history || hasinputs - uind + 1 - else - uind - end if isscalar function (expr) Func( @@ -349,15 +280,11 @@ function wrap_array_vars( [], Let( vcat( - [k ← :(view($(expr.args[uind].name), $v)) for (k, v) in array_vars], - [k ← :(view($(expr.args[inputind].name), $v)) - for (k, v) in input_vars], - [k ← :(reshape( - view($(expr.args[params_offset + buffer_idx].name), $idxs), - $sz)) - for (k, (idxs, buffer_idx, sz)) in array_parameters], - [k ← Code.MakeArray(v, symtype(k)) - for (k, v) in other_array_parameters] + [sym ← :(view($(expr.args[i].name), $idxs)) + for (sym, (i, idxs)) in viewsyms], + [sym ← + MakeArray([expr.args[bufi].elems[vali] for (bufi, vali) in idxs], + expr.args[idxs[1][1]]) for (sym, idxs) in splitsyms] ), expr.body, false @@ -371,15 +298,11 @@ function wrap_array_vars( [], Let( vcat( - [k ← :(view($(expr.args[uind].name), $v)) for (k, v) in array_vars], - [k ← :(view($(expr.args[inputind].name), $v)) - for (k, v) in input_vars], - [k ← :(reshape( - view($(expr.args[params_offset + buffer_idx].name), $idxs), - $sz)) - for (k, (idxs, buffer_idx, sz)) in array_parameters], - [k ← Code.MakeArray(v, symtype(k)) - for (k, v) in other_array_parameters] + [sym ← :(view($(expr.args[i].name), $idxs)) + for (sym, (i, idxs)) in viewsyms], + [sym ← + MakeArray([expr.args[bufi].elems[vali] for (bufi, vali) in idxs], + expr.args[idxs[1][1]]) for (sym, idxs) in splitsyms] ), expr.body, false @@ -392,17 +315,11 @@ function wrap_array_vars( [], Let( vcat( - [k ← :(view($(expr.args[uind + 1].name), $v)) - for (k, v) in array_vars], - [k ← :(view($(expr.args[inputind + 1].name), $v)) - for (k, v) in input_vars], - [k ← :(reshape( - view($(expr.args[params_offset + buffer_idx + 1].name), - $idxs), - $sz)) - for (k, (idxs, buffer_idx, sz)) in array_parameters], - [k ← Code.MakeArray(v, symtype(k)) - for (k, v) in other_array_parameters] + [sym ← :(view($(expr.args[i + 1].name), $idxs)) + for (sym, (i, idxs)) in viewsyms], + [sym ← MakeArray( + [expr.args[bufi + 1].elems[vali] for (bufi, vali) in idxs], + expr.args[idxs[1][1] + 1]) for (sym, idxs) in splitsyms] ), expr.body, false diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index e5db9a71b6..eaf31a9c5f 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -62,8 +62,6 @@ function Base.hash(a::FunctionalAffect, s::UInt) hash(a.ctx, s) end -has_functional_affect(cb) = affects(cb) isa FunctionalAffect - namespace_affect(affect, s) = namespace_equation(affect, s) function namespace_affect(affect::FunctionalAffect, s) FunctionalAffect(func(affect), @@ -75,6 +73,10 @@ function namespace_affect(affect::FunctionalAffect, s) context(affect)) end +function has_functional_affect(cb) + (affects(cb) isa FunctionalAffect || affects(cb) isa ImperativeAffect) +end + #################################### continuous events ##################################### const NULL_AFFECT = Equation[] @@ -88,17 +90,26 @@ By default `affect_neg = affect`; to only get rising edges specify `affect_neg = Assume without loss of generality that the equation is of the form `c(u,p,t) ~ 0`; we denote the integrator state as `i.u`. For compactness, we define `prev_sign = sign(c(u[t-1], p[t-1], t-1))` and `cur_sign = sign(c(u[t], p[t], t))`. A condition edge will be detected and the callback will be invoked iff `prev_sign * cur_sign <= 0`. +The positive edge `affect` will be triggered iff an edge is detected and if `prev_sign < 0`; similarly, `affect_neg` will be +triggered iff an edge is detected and `prev_sign > 0`. + Inter-sample condition activation is not guaranteed; for example if we use the dirac delta function as `c` to insert a sharp discontinuity between integrator steps (which in this example would not normally be identified by adaptivity) then the condition is not guaranteed to be triggered. Once detected the integrator will "wind back" through a root-finding process to identify the point when the condition became active; the method used -is specified by `rootfind` from [`SciMLBase.RootfindOpt`](@ref). Multiple callbacks in the same system with different `rootfind` operations will be resolved -into separate VectorContinuousCallbacks in the enumeration order of `SciMLBase.RootfindOpt`, which may cause some callbacks to not fire if several become -active at the same instant. See the `SciMLBase` documentation for more information on the semantic rules. - -The positive edge `affect` will be triggered iff an edge is detected and if `prev_sign < 0`; similarly, `affect_neg` will be -triggered iff an edge is detected `prev_sign > 0`. +is specified by `rootfind` from [`SciMLBase.RootfindOpt`](@ref). If we denote the time when the condition becomes active as `tc`, +the value in the integrator after windback will be: +* `u[tc-epsilon], p[tc-epsilon], tc` if `LeftRootFind` is used, +* `u[tc+epsilon], p[tc+epsilon], tc` if `RightRootFind` is used, +* or `u[t], p[t], t` if `NoRootFind` is used. +For example, if we want to detect when an unknown variable `x` satisfies `x > 0` using the condition `x ~ 0` on a positive edge (that is, `D(x) > 0`), +then left root finding will get us `x=-epsilon`, right root finding `x=epsilon` and no root finding will produce whatever the next step of the integrator was after +it passed through 0. + +Multiple callbacks in the same system with different `rootfind` operations will be grouped +by their `rootfind` value into separate VectorContinuousCallbacks in the enumeration order of `SciMLBase.RootfindOpt`. This may cause some callbacks to not fire if several become +active at the same instant. See the `SciMLBase` documentation for more information on the semantic rules. Affects (i.e. `affect` and `affect_neg`) can be specified as either: * A list of equations that should be applied when the callback is triggered (e.g. `x ~ 3, y ~ 7`) which must be of the form `unknown ~ observed value` where each `unknown` appears only once. Equations will be applied in the order that they appear in the vector; parameters and state updates will become immediately visible to following equations. @@ -108,6 +119,7 @@ Affects (i.e. `affect` and `affect_neg`) can be specified as either: + `read_parameters` is a vector of the parameters that are *used* by `f!`. Their indices are passed to `f` in `p` similarly to the indices of `unknowns` passed in `u`. + `modified_parameters` is a vector of the parameters that are *modified* by `f!`. Note that a parameter will not appear in `p` if it only appears in `modified_parameters`; it must appear in both `parameters` and `modified_parameters` if it is used in the affect definition. + `ctx` is a user-defined context object passed to `f!` when invoked. This value is aliased for each problem. +* A [`ImperativeAffect`](@ref); refer to its documentation for details. DAEs will be reinitialized using `reinitializealg` (which defaults to `SciMLBase.CheckInit`) after callbacks are applied. This reinitialization algorithm ensures that the DAE is satisfied after the callback runs. The default value of `CheckInit` will simply validate @@ -119,10 +131,10 @@ will run as soon as the solver starts, while finalization affects will be execut """ struct SymbolicContinuousCallback eqs::Vector{Equation} - initialize::Union{Vector{Equation}, FunctionalAffect} - finalize::Union{Vector{Equation}, FunctionalAffect} - affect::Union{Vector{Equation}, FunctionalAffect} - affect_neg::Union{Vector{Equation}, FunctionalAffect, Nothing} + initialize::Union{Vector{Equation}, FunctionalAffect, ImperativeAffect} + finalize::Union{Vector{Equation}, FunctionalAffect, ImperativeAffect} + affect::Union{Vector{Equation}, FunctionalAffect, ImperativeAffect} + affect_neg::Union{Vector{Equation}, FunctionalAffect, ImperativeAffect, Nothing} rootfind::SciMLBase.RootfindOpt reinitializealg::SciMLBase.DAEInitializationAlgorithm function SymbolicContinuousCallback(; @@ -369,7 +381,7 @@ SymbolicDiscreteCallback(cb::SymbolicDiscreteCallback) = cb # passthrough function Base.show(io::IO, db::SymbolicDiscreteCallback) println(io, "condition: ", db.condition) println(io, "affects:") - if db.affects isa FunctionalAffect + if db.affects isa FunctionalAffect || db.affects isa ImperativeAffect # TODO println(io, " ", db.affects) else @@ -722,6 +734,7 @@ function generate_single_rootfinding_callback( else initfn = user_initfun end + return ContinuousCallback( cond, affect_function.affect, affect_function.affect_neg, rootfind = cb.rootfind, initialize = initfn, @@ -755,7 +768,6 @@ function generate_vector_rootfinding_callback( finalize::Union{Function, Nothing}}[ compile_affect_fn(cb, sys, dvs, ps, kwargs) for cb in cbs] - cond = function (out, u, t, integ) rf_ip(out, u, parameter_values(integ), t) end @@ -801,31 +813,21 @@ function generate_vector_rootfinding_callback( initialize = nothing if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing initialize = handle_optional_setup_fn( - map( - (cb, fn) -> begin - if (save_idxs = get(ic.callback_to_clocks, cb, nothing)) !== nothing - let save_idxs = save_idxs - if !isnothing(fn.initialize) - (i) -> begin - fn.initialize(i) - for idx in save_idxs - SciMLBase.save_discretes!(i, idx) - end - end - else - (i) -> begin - for idx in save_idxs - SciMLBase.save_discretes!(i, idx) - end - end + map(cbs, affect_functions) do cb, fn + if (save_idxs = get(ic.callback_to_clocks, cb, nothing)) !== nothing + let save_idxs = save_idxs + custom_init = fn.initialize + (i) -> begin + !isnothing(custom_init) && custom_init(i) + for idx in save_idxs + SciMLBase.save_discretes!(i, idx) end end - else - fn.initialize end - end, - cbs, - affect_functions), + else + fn.initialize + end + end, SciMLBase.INITIALIZE_DEFAULT) else @@ -847,6 +849,13 @@ function compile_affect_fn(cb, sys::AbstractTimeDependentSystem, dvs, ps, kwargs eq_aff = affects(cb) eq_neg_aff = affect_negs(cb) affect = compile_affect(eq_aff, cb, sys, dvs, ps; expression = Val{false}, kwargs...) + function compile_optional_affect(aff, default = nothing) + if isnothing(aff) || aff == default + return nothing + else + return compile_affect(aff, cb, sys, dvs, ps; expression = Val{false}, kwargs...) + end + end if eq_neg_aff === eq_aff affect_neg = affect else @@ -944,10 +953,48 @@ function compile_user_affect(affect::FunctionalAffect, cb, sys, dvs, ps; kwargs. end end +function invalid_variables(sys, expr) + filter(x -> !any(isequal(x), all_symbols(sys)), reduce(vcat, vars(expr); init = [])) +end +function unassignable_variables(sys, expr) + assignable_syms = reduce( + vcat, Symbolics.scalarize.(vcat(unknowns(sys), parameters(sys))); init = []) + written = reduce(vcat, Symbolics.scalarize.(vars(expr)); init = []) + return filter( + x -> !any(isequal(x), assignable_syms), written) +end + +@generated function _generated_writeback(integ, setters::NamedTuple{NS1, <:Tuple}, + values::NamedTuple{NS2, <:Tuple}) where {NS1, NS2} + setter_exprs = [] + for name in NS2 + if !(name in NS1) + missing_name = "Tried to write back to $name from affect; only declared states ($NS1) may be written to." + error(missing_name) + end + push!(setter_exprs, :(setters.$name(integ, values.$name))) + end + return :(begin + $(setter_exprs...) + end) +end + +function check_assignable(sys, sym) + if symbolic_type(sym) == ScalarSymbolic() + is_variable(sys, sym) || is_parameter(sys, sym) + elseif symbolic_type(sym) == ArraySymbolic() + is_variable(sys, sym) || is_parameter(sys, sym) || + all(x -> check_assignable(sys, x), collect(sym)) + elseif sym isa Union{AbstractArray, Tuple} + all(x -> check_assignable(sys, x), sym) + else + false + end +end + function compile_affect(affect::FunctionalAffect, cb, sys, dvs, ps; kwargs...) compile_user_affect(affect, cb, sys, dvs, ps; kwargs...) end - function _compile_optional_affect(default, aff, cb, sys, dvs, ps; kwargs...) if isnothing(aff) || aff == default return nothing diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 3a502d3183..87ee6e823d 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1384,6 +1384,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, initialization_eqs = [], fully_determined = nothing, check_units = true, + use_scc = true, kwargs...) where {iip, specialize} if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`") @@ -1393,16 +1394,18 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, elseif isempty(u0map) && get_initializesystem(sys) === nothing isys = structural_simplify( generate_initializesystem( - sys; initialization_eqs, check_units, pmap = parammap); fully_determined) + sys; initialization_eqs, check_units, pmap = parammap, + guesses, extra_metadata = (; use_scc)); fully_determined) else isys = structural_simplify( generate_initializesystem( - sys; u0map, initialization_eqs, check_units, pmap = parammap); fully_determined) + sys; u0map, initialization_eqs, check_units, + pmap = parammap, guesses, extra_metadata = (; use_scc)); fully_determined) end ts = get_tearing_state(isys) - if warn_initialize_determined && - (unassigned_vars = StructuralTransformations.singular_check(ts); !isempty(unassigned_vars)) + unassigned_vars = StructuralTransformations.singular_check(ts) + if warn_initialize_determined && !isempty(unassigned_vars) errmsg = """ The initialization system is structurally singular. Guess values may \ significantly affect the initial values of the ODE. The problematic variables \ @@ -1424,11 +1427,17 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, neqs = length(equations(isys)) nunknown = length(unknowns(isys)) + if use_scc + scc_message = "`SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here. " + else + scc_message = "" + end + if warn_initialize_determined && neqs > nunknown - @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true" + @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. $(scc_message)To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true" end if warn_initialize_determined && neqs < nunknown - @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true" + @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. $(scc_message)To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true" end parammap = parammap isa DiffEqBase.NullParameters || isempty(parammap) ? @@ -1464,9 +1473,20 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, end for (k, v) in u0map) end - if neqs == nunknown - NonlinearProblem(isys, u0map, parammap; kwargs...) + + TProb = if neqs == nunknown && isempty(unassigned_vars) + if use_scc && neqs > 0 + if is_split(isys) + SCCNonlinearProblem + else + @warn "`SCCNonlinearProblem` can only be used with `split = true` systems. Simplify your `ODESystem` with `split = true` or pass `use_scc = false` to disable this warning" + NonlinearProblem + end + else + NonlinearProblem + end else - NonlinearLeastSquaresProblem(isys, u0map, parammap; kwargs...) + NonlinearLeastSquaresProblem end + TProb(isys, u0map, parammap; kwargs...) end diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl index 011ba6e216..ac47f4c45c 100644 --- a/src/systems/diffeqs/sdesystem.jl +++ b/src/systems/diffeqs/sdesystem.jl @@ -263,6 +263,47 @@ function Base.:(==)(sys1::SDESystem, sys2::SDESystem) all(s1 == s2 for (s1, s2) in zip(get_systems(sys1), get_systems(sys2))) end +""" + function ODESystem(sys::SDESystem) + +Convert an `SDESystem` to the equivalent `ODESystem` using `@brownian` variables instead +of noise equations. The returned system will not be `iscomplete` and will not have an +index cache, regardless of `iscomplete(sys)`. +""" +function ODESystem(sys::SDESystem) + neqs = get_noiseeqs(sys) + eqs = equations(sys) + is_scalar_noise = get_is_scalar_noise(sys) + nbrownian = if is_scalar_noise + length(neqs) + else + size(neqs, 2) + end + brownvars = map(1:nbrownian) do i + name = gensym(Symbol(:brown_, i)) + only(@brownian $name) + end + if is_scalar_noise + brownterms = reduce(+, neqs .* brownvars; init = 0) + neweqs = map(eqs) do eq + eq.lhs ~ eq.rhs + brownterms + end + else + if neqs isa AbstractVector + neqs = reshape(neqs, (length(neqs), 1)) + end + brownterms = neqs * brownvars + neweqs = map(eqs, brownterms) do eq, brown + eq.lhs ~ eq.rhs + brown + end + end + newsys = ODESystem(neweqs, get_iv(sys), unknowns(sys), parameters(sys); + parameter_dependencies = parameter_dependencies(sys), defaults = defaults(sys), + continuous_events = continuous_events(sys), discrete_events = discrete_events(sys), + name = nameof(sys), description = description(sys), metadata = get_metadata(sys)) + @set newsys.parent = sys +end + function __num_isdiag_noise(mat) for i in axes(mat, 1) nnz = 0 diff --git a/src/systems/imperative_affect.jl b/src/systems/imperative_affect.jl new file mode 100644 index 0000000000..2f489913ad --- /dev/null +++ b/src/systems/imperative_affect.jl @@ -0,0 +1,218 @@ + +""" + ImperativeAffect(f::Function; modified::NamedTuple, observed::NamedTuple, ctx) + +`ImperativeAffect` is a helper for writing affect functions that will compute observed values and +ensure that modified values are correctly written back into the system. The affect function `f` needs to have +the signature + +``` + f(modified::NamedTuple, observed::NamedTuple, ctx, integrator)::NamedTuple +``` + +The function `f` will be called with `observed` and `modified` `NamedTuple`s that are derived from their respective `NamedTuple` definitions. +Each declaration`NamedTuple` should map an expression to a symbol; for example if we pass `observed=(; x = a + b)` this will alias the result of executing `a+b` in the system as `x` +so the value of `a + b` will be accessible as `observed.x` in `f`. `modified` currently restricts symbolic expressions to only bare variables, so only tuples of the form +`(; x = y)` or `(; x)` (which aliases `x` as itself) are allowed. + +The argument NamedTuples (for instance `(;x=y)`) will be populated with the declared values on function entry; if we require `(;x=y)` in `observed` and `y=2`, for example, +then the NamedTuple `(;x=2)` will be passed as `observed` to the affect function `f`. + +The NamedTuple returned from `f` includes the values to be written back to the system after `f` returns. For example, if we want to update the value of `x` to be the result of `x + y` we could write + + ImperativeAffect(observed=(; x_plus_y = x + y), modified=(; x)) do m, o + @set! m.x = o.x_plus_y + end + +Where we use Setfield to copy the tuple `m` with a new value for `x`, then return the modified value of `m`. All values updated by the tuple must have names originally declared in +`modified`; a runtime error will be produced if a value is written that does not appear in `modified`. The user can dynamically decide not to write a value back by not including it +in the returned tuple, in which case the associated field will not be updated. +""" +@kwdef struct ImperativeAffect + f::Any + obs::Vector + obs_syms::Vector{Symbol} + modified::Vector + mod_syms::Vector{Symbol} + ctx::Any + skip_checks::Bool +end + +function ImperativeAffect(f::Function; + observed::NamedTuple = NamedTuple{()}(()), + modified::NamedTuple = NamedTuple{()}(()), + ctx = nothing, + skip_checks = false) + ImperativeAffect(f, + collect(values(observed)), collect(keys(observed)), + collect(values(modified)), collect(keys(modified)), + ctx, skip_checks) +end +function ImperativeAffect(f::Function, modified::NamedTuple; + observed::NamedTuple = NamedTuple{()}(()), ctx = nothing, skip_checks = false) + ImperativeAffect( + f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) +end +function ImperativeAffect( + f::Function, modified::NamedTuple, observed::NamedTuple; ctx = nothing, skip_checks = false) + ImperativeAffect( + f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) +end +function ImperativeAffect( + f::Function, modified::NamedTuple, observed::NamedTuple, ctx; skip_checks = false) + ImperativeAffect( + f, observed = observed, modified = modified, ctx = ctx, skip_checks = skip_checks) +end + +function Base.show(io::IO, mfa::ImperativeAffect) + obs_vals = join(map((ob, nm) -> "$ob => $nm", mfa.obs, mfa.obs_syms), ", ") + mod_vals = join(map((md, nm) -> "$md => $nm", mfa.modified, mfa.mod_syms), ", ") + affect = mfa.f + print(io, + "ImperativeAffect(observed: [$obs_vals], modified: [$mod_vals], affect:$affect)") +end +func(f::ImperativeAffect) = f.f +context(a::ImperativeAffect) = a.ctx +observed(a::ImperativeAffect) = a.obs +observed_syms(a::ImperativeAffect) = a.obs_syms +discretes(a::ImperativeAffect) = filter(ModelingToolkit.isparameter, a.modified) +modified(a::ImperativeAffect) = a.modified +modified_syms(a::ImperativeAffect) = a.mod_syms + +function Base.:(==)(a1::ImperativeAffect, a2::ImperativeAffect) + isequal(a1.f, a2.f) && isequal(a1.obs, a2.obs) && isequal(a1.modified, a2.modified) && + isequal(a1.obs_syms, a2.obs_syms) && isequal(a1.mod_syms, a2.mod_syms) && + isequal(a1.ctx, a2.ctx) +end + +function Base.hash(a::ImperativeAffect, s::UInt) + s = hash(a.f, s) + s = hash(a.obs, s) + s = hash(a.obs_syms, s) + s = hash(a.modified, s) + s = hash(a.mod_syms, s) + hash(a.ctx, s) +end + +namespace_affects(af::ImperativeAffect, s) = namespace_affect(af, s) +function namespace_affect(affect::ImperativeAffect, s) + ImperativeAffect(func(affect), + namespace_expr.(observed(affect), (s,)), + observed_syms(affect), + renamespace.((s,), modified(affect)), + modified_syms(affect), + context(affect), + affect.skip_checks) +end + +function compile_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs...) + compile_user_affect(affect, cb, sys, dvs, ps; kwargs...) +end + +function compile_user_affect(affect::ImperativeAffect, cb, sys, dvs, ps; kwargs...) + #= + Implementation sketch: + generate observed function (oop), should save to a component array under obs_syms + do the same stuff as the normal FA for pars_syms + call the affect method + unpack and apply the resulting values + =# + function check_dups(syms, exprs) # = (syms_dedup, exprs_dedup) + seen = Set{Symbol}() + syms_dedup = [] + exprs_dedup = [] + for (sym, exp) in Iterators.zip(syms, exprs) + if !in(sym, seen) + push!(syms_dedup, sym) + push!(exprs_dedup, exp) + push!(seen, sym) + elseif !affect.skip_checks + @warn "Expression $(expr) is aliased as $sym, which has already been used. The first definition will be used." + end + end + return (syms_dedup, exprs_dedup) + end + + obs_exprs = observed(affect) + if !affect.skip_checks + for oexpr in obs_exprs + invalid_vars = invalid_variables(sys, oexpr) + if length(invalid_vars) > 0 + error("Observed equation $(oexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing).") + end + end + end + obs_syms = observed_syms(affect) + obs_syms, obs_exprs = check_dups(obs_syms, obs_exprs) + + mod_exprs = modified(affect) + if !affect.skip_checks + for mexpr in mod_exprs + if !check_assignable(sys, mexpr) + @warn ("Expression $mexpr cannot be assigned to; currently only unknowns and parameters may be updated by an affect.") + end + invalid_vars = unassignable_variables(sys, mexpr) + if length(invalid_vars) > 0 + error("Modified equation $(mexpr) in affect refers to missing variable(s) $(invalid_vars); the variables may not have been added (e.g. if a component is missing) or they may have been reduced away.") + end + end + end + mod_syms = modified_syms(affect) + mod_syms, mod_exprs = check_dups(mod_syms, mod_exprs) + + overlapping_syms = intersect(mod_syms, obs_syms) + if length(overlapping_syms) > 0 && !affect.skip_checks + @warn "The symbols $overlapping_syms are declared as both observed and modified; this is a code smell because it becomes easy to confuse them and assign/not assign a value." + end + + # sanity checks done! now build the data and update function for observed values + mkzero(sz) = + if sz === () + 0.0 + else + zeros(sz) + end + obs_fun = build_explicit_observed_function( + sys, Symbolics.scalarize.(obs_exprs); + mkarray = (es, _) -> MakeTuple(es)) + obs_sym_tuple = (obs_syms...,) + + # okay so now to generate the stuff to assign it back into the system + mod_pairs = mod_exprs .=> mod_syms + mod_names = (mod_syms...,) + mod_og_val_fun = build_explicit_observed_function( + sys, Symbolics.scalarize.(first.(mod_pairs)); + mkarray = (es, _) -> MakeTuple(es)) + + upd_funs = NamedTuple{mod_names}((setu.((sys,), first.(mod_pairs))...,)) + + if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing + save_idxs = get(ic.callback_to_clocks, cb, Int[]) + else + save_idxs = Int[] + end + + let user_affect = func(affect), ctx = context(affect) + function (integ) + # update the to-be-mutated values; this ensures that if you do a no-op then nothing happens + modvals = mod_og_val_fun(integ.u, integ.p, integ.t) + upd_component_array = NamedTuple{mod_names}(modvals) + + # update the observed values + obs_component_array = NamedTuple{obs_sym_tuple}(obs_fun( + integ.u, integ.p, integ.t)) + + # let the user do their thing + upd_vals = user_affect(upd_component_array, obs_component_array, ctx, integ) + + # write the new values back to the integrator + _generated_writeback(integ, upd_funs, upd_vals) + + for idx in save_idxs + SciMLBase.save_discretes!(integ, idx) + end + end + end +end + +scalarize_affects(affects::ImperativeAffect) = affects diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl index ab0dd08764..4fc2f18bfd 100644 --- a/src/systems/index_cache.jl +++ b/src/systems/index_cache.jl @@ -116,7 +116,7 @@ function IndexCache(sys::AbstractSystem) for affect in affs if affect isa Equation is_parameter(sys, affect.lhs) && push!(discs, affect.lhs) - elseif affect isa FunctionalAffect + elseif affect isa FunctionalAffect || affect isa ImperativeAffect union!(discs, unwrap.(discretes(affect))) else error("Unhandled affect type $(typeof(affect))") @@ -594,3 +594,28 @@ function reorder_dimension_by_tunables( reorder_dimension_by_tunables!(buffer, sys, arr, syms; dim) return buffer end + +function subset_unknowns_observed( + ic::IndexCache, sys::AbstractSystem, newunknowns, newobsvars) + unknown_idx = copy(ic.unknown_idx) + empty!(unknown_idx) + for (i, sym) in enumerate(newunknowns) + ttsym = default_toterm(sym) + rsym = renamespace(sys, sym) + rttsym = renamespace(sys, ttsym) + unknown_idx[sym] = unknown_idx[ttsym] = unknown_idx[rsym] = unknown_idx[rttsym] = i + end + observed_syms_to_timeseries = copy(ic.observed_syms_to_timeseries) + empty!(observed_syms_to_timeseries) + for sym in newobsvars + ttsym = default_toterm(sym) + rsym = renamespace(sys, sym) + rttsym = renamespace(sys, ttsym) + for s in (sym, ttsym, rsym, rttsym) + observed_syms_to_timeseries[s] = ic.observed_syms_to_timeseries[sym] + end + end + ic = @set ic.unknown_idx = unknown_idx + @set! ic.observed_syms_to_timeseries = observed_syms_to_timeseries + return ic +end diff --git a/src/systems/nonlinear/homotopy_continuation.jl b/src/systems/nonlinear/homotopy_continuation.jl new file mode 100644 index 0000000000..03aeed1edf --- /dev/null +++ b/src/systems/nonlinear/homotopy_continuation.jl @@ -0,0 +1,538 @@ +""" +$(TYPEDEF) + +A type of Nonlinear problem which specializes on polynomial systems and uses +HomotopyContinuation.jl to solve the system. Requires importing HomotopyContinuation.jl to +create and solve. +""" +struct HomotopyContinuationProblem{uType, H, D, O, SS, U} <: + SciMLBase.AbstractNonlinearProblem{uType, true} + """ + The initial values of states in the system. If there are multiple real roots of + the system, the one closest to this point is returned. + """ + u0::uType + """ + A subtype of `HomotopyContinuation.AbstractSystem` to solve. Also contains the + parameter object. + """ + homotopy_continuation_system::H + """ + A function with signature `(u, p) -> resid`. In case of rational functions, this + is used to rule out roots of the system which would cause the denominator to be + zero. + """ + denominator::D + """ + The `NonlinearSystem` used to create this problem. Used for symbolic indexing. + """ + sys::NonlinearSystem + """ + A function which generates and returns observed expressions for the given system. + """ + obsfn::O + """ + The HomotopyContinuation.jl solver and start system, obtained through + `HomotopyContinuation.solver_startsystems`. + """ + solver_and_starts::SS + """ + A function which takes a solution of the transformed system, and returns a vector + of solutions for the original system. This is utilized when converting systems + to polynomials. + """ + unpack_solution::U +end + +function HomotopyContinuationProblem(::AbstractSystem, _u0, _p; kwargs...) + error("HomotopyContinuation.jl is required to create and solve `HomotopyContinuationProblem`s. Please run `Pkg.add(\"HomotopyContinuation\")` to continue.") +end + +""" + $(TYPEDSIGNATURES) + +Utility function for `safe_HomotopyContinuationProblem`, implemented in the extension. +""" +function _safe_HomotopyContinuationProblem end + +""" + $(TYPEDSIGNATURES) + +Return a `HomotopyContinuationProblem` if the extension is loaded and the system is +polynomial. If the extension is not loaded, return `nothing`. If the system is not +polynomial, return the appropriate `NotPolynomialError`. +""" +function safe_HomotopyContinuationProblem(sys::NonlinearSystem, args...; kwargs...) + if Base.get_extension(ModelingToolkit, :MTKHomotopyContinuationExt) === nothing + return nothing + end + return _safe_HomotopyContinuationProblem(sys, args...; kwargs...) +end + +SymbolicIndexingInterface.symbolic_container(p::HomotopyContinuationProblem) = p.sys +SymbolicIndexingInterface.state_values(p::HomotopyContinuationProblem) = p.u0 +function SymbolicIndexingInterface.set_state!(p::HomotopyContinuationProblem, args...) + set_state!(p.u0, args...) +end +function SymbolicIndexingInterface.parameter_values(p::HomotopyContinuationProblem) + parameter_values(p.homotopy_continuation_system) +end +function SymbolicIndexingInterface.set_parameter!(p::HomotopyContinuationProblem, args...) + set_parameter!(parameter_values(p), args...) +end +function SymbolicIndexingInterface.observed(p::HomotopyContinuationProblem, sym) + if p.obsfn !== nothing + return p.obsfn(sym) + else + return SymbolicIndexingInterface.observed(p.sys, sym) + end +end + +function contains_variable(x, wrt) + any(y -> occursin(y, x), wrt) +end + +""" +Possible reasons why a term is not polynomial +""" +EnumX.@enumx NonPolynomialReason begin + """ + Exponent of an expression involving unknowns is not an integer. + """ + NonIntegerExponent + """ + Exponent is an expression containing unknowns. + """ + ExponentContainsUnknowns + """ + The base of an exponent is not a polynomial in the unknowns. + """ + BaseNotPolynomial + """ + An expression involves a non-polynomial operation involving unknowns. + """ + UnrecognizedOperation +end + +function display_reason(reason::NonPolynomialReason.T, sym) + if reason == NonPolynomialReason.NonIntegerExponent + pow = arguments(sym)[2] + "In $sym: Exponent $pow is not an integer" + elseif reason == NonPolynomialReason.ExponentContainsUnknowns + pow = arguments(sym)[2] + "In $sym: Exponent $pow contains unknowns of the system" + elseif reason == NonPolynomialReason.BaseNotPolynomial + base = arguments(sym)[1] + "In $sym: Base $base is not a polynomial in the unknowns" + elseif reason == NonPolynomialReason.UnrecognizedOperation + op = operation(sym) + """ + In $sym: Operation $op is not recognized. Allowed polynomial operations are \ + `*, /, +, -, ^`. + """ + else + error("This should never happen. Please open an issue in ModelingToolkit.jl.") + end +end + +""" + $(TYPEDEF) + +Information about an expression about its polynomial nature. +""" +mutable struct PolynomialData + """ + A list of all non-polynomial terms in the expression. + """ + non_polynomial_terms::Vector{BasicSymbolic} + """ + Corresponding to `non_polynomial_terms`, a list of reasons why they are + not polynomial. + """ + reasons::Vector{NonPolynomialReason.T} + """ + Whether the polynomial contains parametric exponents of unknowns. + """ + has_parametric_exponent::Bool +end + +PolynomialData() = PolynomialData(BasicSymbolic[], NonPolynomialReason.T[], false) + +abstract type PolynomialTransformationError <: Exception end + +struct MultivarTerm <: PolynomialTransformationError + term::Any + vars::Any +end + +function Base.showerror(io::IO, err::MultivarTerm) + println(io, + "Cannot convert system to polynomial: Found term $(err.term) which is a function of multiple unknowns $(err.vars).") +end + +struct MultipleTermsOfSameVar <: PolynomialTransformationError + terms::Any + var::Any +end + +function Base.showerror(io::IO, err::MultipleTermsOfSameVar) + println(io, + "Cannot convert system to polynomial: Found multiple non-polynomial terms $(err.terms) involving the same unknown $(err.var).") +end + +struct SymbolicSolveFailure <: PolynomialTransformationError + term::Any + var::Any +end + +function Base.showerror(io::IO, err::SymbolicSolveFailure) + println(io, + "Cannot convert system to polynomial: Unable to symbolically solve $(err.term) for $(err.var).") +end + +struct NemoNotLoaded <: PolynomialTransformationError end + +function Base.showerror(io::IO, err::NemoNotLoaded) + println(io, + "ModelingToolkit may be able to solve this system as a polynomial system if `Nemo` is loaded. Run `import Nemo` and try again.") +end + +struct VariablesAsPolyAndNonPoly <: PolynomialTransformationError + vars::Any +end + +function Base.showerror(io::IO, err::VariablesAsPolyAndNonPoly) + println(io, + "Cannot convert convert system to polynomial: Variables $(err.vars) occur in both polynomial and non-polynomial terms in the system.") +end + +struct NotPolynomialError <: Exception + transformation_err::Union{PolynomialTransformationError, Nothing} + eq::Vector{Equation} + data::Vector{PolynomialData} +end + +function Base.showerror(io::IO, err::NotPolynomialError) + if err.transformation_err !== nothing + Base.showerror(io, err.transformation_err) + end + for (eq, data) in zip(err.eq, err.data) + if isempty(data.non_polynomial_terms) + continue + end + println(io, + "Equation $(eq) is not a polynomial in the unknowns for the following reasons:") + for (term, reason) in zip(data.non_polynomial_terms, data.reasons) + println(io, display_reason(reason, term)) + end + end +end + +function is_polynomial!(data, y, wrt) + process_polynomial!(data, y, wrt) + isempty(data.reasons) +end + +""" +$(TYPEDSIGNATURES) + +Return information about the polynmial `x` with respect to variables in `wrt`, +writing said information to `data`. +""" +function process_polynomial!(data::PolynomialData, x, wrt) + x = unwrap(x) + symbolic_type(x) == NotSymbolic() && return true + iscall(x) || return true + contains_variable(x, wrt) || return true + any(isequal(x), wrt) && return true + + if operation(x) in (*, +, -, /) + # `map` because `all` will early exit, but we want to search + # through everything to get all the non-polynomial terms + return all(map(y -> is_polynomial!(data, y, wrt), arguments(x))) + end + if operation(x) == (^) + b, p = arguments(x) + is_pow_integer = symtype(p) <: Integer + if !is_pow_integer + push!(data.non_polynomial_terms, x) + push!(data.reasons, NonPolynomialReason.NonIntegerExponent) + end + if symbolic_type(p) != NotSymbolic() + data.has_parametric_exponent = true + end + + exponent_has_unknowns = contains_variable(p, wrt) + if exponent_has_unknowns + push!(data.non_polynomial_terms, x) + push!(data.reasons, NonPolynomialReason.ExponentContainsUnknowns) + end + base_polynomial = is_polynomial!(data, b, wrt) + return base_polynomial && !exponent_has_unknowns && is_pow_integer + end + push!(data.non_polynomial_terms, x) + push!(data.reasons, NonPolynomialReason.UnrecognizedOperation) + return false +end + +""" + $(TYPEDEF) + +Information about how an unknown in the system is substituted for a non-polynomial +expression to turn the system into a polynomial. Used in `PolynomialTransformation`. +""" +struct PolynomialTransformationData + """ + The new variable to use as an unknown of the transformed system. + """ + new_var::BasicSymbolic + """ + The non-polynomial expression being substituted. + """ + term::BasicSymbolic + """ + A vector of expressions corresponding to the solutions of + the non-polynomial expression `term` in terms of the new unknown `new_var`, + used to backsolve for the original unknown of the system. + """ + inv_term::Vector{BasicSymbolic} +end + +""" + $(TYPEDEF) + +Information representing how to transform a `NonlinearSystem` into a polynomial +system. +""" +struct PolynomialTransformation + """ + Substitutions mapping non-polynomial terms to temporary unknowns. The system + is a polynomial in the new unknowns. Currently, each non-polynomial term is a + function of a single unknown of the original system. + """ + substitution_rules::Dict{BasicSymbolic, BasicSymbolic} + """ + A vector of expressions involving unknowns of the transformed system, mapping + back to solutions of the original system. + """ + all_solutions::Vector{Vector{BasicSymbolic}} + """ + The new unknowns of the transformed system. + """ + new_dvs::Vector{BasicSymbolic} + """ + The polynomial data for each equation. + """ + polydata::Vector{PolynomialData} +end + +function PolynomialTransformation(sys::NonlinearSystem) + # we need to consider `full_equations` because observed also should be + # polynomials (if used in equations) and we don't know if observed is used + # in denominator. + # This is not the most efficient, and would be improved significantly with + # CSE/hashconsing. + eqs = full_equations(sys) + dvs = unknowns(sys) + + # Collect polynomial information about all equations + polydata = map(eqs) do eq + data = PolynomialData() + process_polynomial!(data, eq.lhs, dvs) + process_polynomial!(data, eq.rhs, dvs) + data + end + + # Get all unique non-polynomial terms + # NOTE: + # Is there a better way to check for uniqueness? `simplify` is relatively slow + # (maybe use the threaded version?) and `expand` can blow up expression size. + # Could metatheory help? + all_non_poly_terms = mapreduce(d -> d.non_polynomial_terms, vcat, polydata) + unique!(all_non_poly_terms) + + # each variable can only be replaced by one non-polynomial expression involving + # that variable. Keep track of this mapping. + var_to_nonpoly = Dict{BasicSymbolic, PolynomialTransformationData}() + + is_poly = true + transformation_err = nothing + for t in all_non_poly_terms + # if the term involves multiple unknowns, we can't invert it + dvs_in_term = map(x -> occursin(x, t), dvs) + if count(dvs_in_term) > 1 + transformation_err = MultivarTerm(t, dvs[dvs_in_term]) + is_poly = false + break + end + # we already have a substitution solving for `var` + var = dvs[findfirst(dvs_in_term)] + if haskey(var_to_nonpoly, var) && !isequal(var_to_nonpoly[var].term, t) + transformation_err = MultipleTermsOfSameVar([t, var_to_nonpoly[var].term], var) + is_poly = false + break + end + # we want to solve `term - new_var` for `var` + new_var = gensym(Symbol(var)) + new_var = unwrap(only(@variables $new_var)) + invterm = Symbolics.ia_solve( + t - new_var, var; complex_roots = false, periodic_roots = false, warns = false) + # if we can't invert it, quit + if invterm === nothing || isempty(invterm) + transformation_err = SymbolicSolveFailure(t, var) + is_poly = false + break + end + # `ia_solve` returns lazy terms i.e. `asin(1.0)` instead of `pi/2` + # this just evaluates the constant expressions + invterm = Symbolics.substitute.(invterm, (Dict(),)) + # RootsOf implies Symbolics couldn't solve the inner polynomial because + # `Nemo` wasn't loaded. + if any(x -> iscall(x) && operation(x) == Symbolics.RootsOf, invterm) + transformation_err = NemoNotLoaded() + is_poly = false + break + end + var_to_nonpoly[var] = PolynomialTransformationData(new_var, t, invterm) + end + + # return the error instead of throwing it, so the user can choose what to do + # without having to catch the exception + if !is_poly + return NotPolynomialError(transformation_err, eqs, polydata) + end + + subrules = Dict{BasicSymbolic, BasicSymbolic}() + # corresponding to each unknown in `dvs`, the list of its possible solutions + # in terms of the new unknown. + combinations = Vector{BasicSymbolic}[] + new_dvs = BasicSymbolic[] + for x in dvs + if haskey(var_to_nonpoly, x) + _data = var_to_nonpoly[x] + # map term to new unknown + subrules[_data.term] = _data.new_var + push!(combinations, _data.inv_term) + push!(new_dvs, _data.new_var) + else + push!(combinations, BasicSymbolic[x]) + push!(new_dvs, x) + end + end + all_solutions = vec(collect.(collect(Iterators.product(combinations...)))) + return PolynomialTransformation(subrules, all_solutions, new_dvs, polydata) +end + +""" + $(TYPEDEF) + +A struct containing the result of transforming a system into a polynomial system +using the appropriate `PolynomialTransformation`. Also contains the denominators +in the equations, to rule out invalid roots. +""" +struct PolynomialTransformationResult + sys::NonlinearSystem + denominators::Vector{BasicSymbolic} +end + +""" + $(TYPEDSIGNATURES) + +Transform the system `sys` with `transformation` and return a +`PolynomialTransformationResult`, or a `NotPolynomialError` if the system cannot +be transformed. +""" +function transform_system(sys::NonlinearSystem, transformation::PolynomialTransformation) + subrules = transformation.substitution_rules + dvs = unknowns(sys) + eqs = full_equations(sys) + polydata = transformation.polydata + new_dvs = transformation.new_dvs + all_solutions = transformation.all_solutions + + eqs2 = Equation[] + denoms = BasicSymbolic[] + for eq in eqs + t = eq.rhs - eq.lhs + t = Symbolics.fixpoint_sub(t, subrules; maxiters = length(dvs)) + # the substituted variable occurs outside the substituted term + poly_and_nonpoly = map(dvs) do x + all(!isequal(x), new_dvs) && occursin(x, t) + end + if any(poly_and_nonpoly) + return NotPolynomialError( + VariablesAsPolyAndNonPoly(dvs[poly_and_nonpoly]), eqs, polydata) + end + num, den = handle_rational_polynomials(t, new_dvs) + # make factors different elements, otherwise the nonzero factors artificially + # inflate the error of the zero factor. + if iscall(den) && operation(den) == * + for arg in arguments(den) + # ignore constant factors + symbolic_type(arg) == NotSymbolic() && continue + push!(denoms, abs(arg)) + end + elseif symbolic_type(den) != NotSymbolic() + push!(denoms, abs(den)) + end + push!(eqs2, 0 ~ num) + end + + sys2 = @set sys.eqs = eqs2 + @set! sys2.unknowns = new_dvs + # remove observed equations to avoid adding them in codegen + @set! sys2.observed = Equation[] + @set! sys2.substitutions = nothing + return PolynomialTransformationResult(sys2, denoms) +end + +""" +$(TYPEDSIGNATURES) + +Given a `x`, a polynomial in variables in `wrt` which may contain rational functions, +express `x` as a single rational function with polynomial `num` and denominator `den`. +Return `(num, den)`. +""" +function handle_rational_polynomials(x, wrt) + x = unwrap(x) + symbolic_type(x) == NotSymbolic() && return x, 1 + iscall(x) || return x, 1 + contains_variable(x, wrt) || return x, 1 + any(isequal(x), wrt) && return x, 1 + + # simplify_fractions cancels out some common factors + # and expands (a / b)^c to a^c / b^c, so we only need + # to handle these cases + x = simplify_fractions(x) + op = operation(x) + args = arguments(x) + + if op == / + # numerator and denominator are trivial + num, den = args + # but also search for rational functions in numerator + n, d = handle_rational_polynomials(num, wrt) + num, den = n, den * d + elseif op == + + num = 0 + den = 1 + + # we don't need to do common denominator + # because we don't care about cases where denominator + # is zero. The expression is zero when all the numerators + # are zero. + for arg in args + n, d = handle_rational_polynomials(arg, wrt) + num += n + den *= d + end + else + return x, 1 + end + # if the denominator isn't a polynomial in `wrt`, better to not include it + # to reduce the size of the gcd polynomial + if !contains_variable(den, wrt) + return num / den, 1 + end + return num, den +end diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index eff19afb07..2344727920 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -11,7 +11,7 @@ function generate_initializesystem(sys::ODESystem; default_dd_guess = 0.0, algebraic_only = false, check_units = true, check_defguess = false, - name = nameof(sys), kwargs...) + name = nameof(sys), extra_metadata = (;), kwargs...) trueobs, eqs = unhack_observed(observed(sys), equations(sys)) vars = unique([unknowns(sys); getfield.(trueobs, :lhs)]) vars_set = Set(vars) # for efficient in-lookup @@ -30,7 +30,8 @@ function generate_initializesystem(sys::ODESystem; # 1) process dummy derivatives and u0map into initialization system eqs_ics = eqs[idxs_alge] # start equation list with algebraic equations defs = copy(defaults(sys)) # copy so we don't modify sys.defaults - guesses = merge(get_guesses(sys), todict(guesses)) + additional_guesses = anydict(guesses) + guesses = merge(get_guesses(sys), additional_guesses) schedule = getfield(sys, :schedule) if !isnothing(schedule) for x in filter(x -> !isnothing(x[1]), schedule.dummy_sub) @@ -178,7 +179,8 @@ function generate_initializesystem(sys::ODESystem; for k in keys(defs) defs[k] = substitute(defs[k], paramsubs) end - meta = InitializationSystemMetadata(Dict{Any, Any}(u0map), Dict{Any, Any}(pmap)) + meta = InitializationSystemMetadata( + anydict(u0map), anydict(pmap), additional_guesses, extra_metadata) return NonlinearSystem(eqs_ics, vars, pars; @@ -193,6 +195,8 @@ end struct InitializationSystemMetadata u0map::Dict{Any, Any} pmap::Dict{Any, Any} + additional_guesses::Dict{Any, Any} + extra_metadata::NamedTuple end function is_parameter_solvable(p, pmap, defs, guesses) @@ -208,17 +212,16 @@ function is_parameter_solvable(p, pmap, defs, guesses) _val1 === nothing && _val2 !== nothing)) && _val3 !== nothing end -function SciMLBase.remake_initializeprob(sys::ODESystem, odefn, u0, t0, p) +function SciMLBase.remake_initialization_data(sys::ODESystem, odefn, u0, t0, p, newu0, newp) if u0 === missing && p === missing - return odefn.initializeprob, odefn.update_initializeprob!, odefn.initializeprobmap, - odefn.initializeprobpmap + return odefn.initialization_data end if !(eltype(u0) <: Pair) && !(eltype(p) <: Pair) oldinitprob = odefn.initializeprob - if oldinitprob === nothing || !SciMLBase.has_sys(oldinitprob.f) || - !(oldinitprob.f.sys isa NonlinearSystem) - return oldinitprob, odefn.update_initializeprob!, odefn.initializeprobmap, - odefn.initializeprobpmap + oldinitprob === nothing && return nothing + if !SciMLBase.has_sys(oldinitprob.f) || !(oldinitprob.f.sys isa NonlinearSystem) + return SciMLBase.OverrideInitData(oldinitprob, odefn.update_initializeprob!, + odefn.initializeprobmap, odefn.initializeprobpmap) end pidxs = ParameterIndex[] pvals = [] @@ -259,79 +262,86 @@ function SciMLBase.remake_initializeprob(sys::ODESystem, odefn, u0, t0, p) newp = remake_buffer( oldinitprob.f.sys, parameter_values(oldinitprob), pidxs, pvals) end - initprob = remake(oldinitprob; u0 = newu0, p = newp) - return initprob, odefn.update_initializeprob!, odefn.initializeprobmap, - odefn.initializeprobpmap - end - if u0 === missing || isempty(u0) - u0 = Dict() - elseif !(eltype(u0) <: Pair) - u0 = Dict(unknowns(sys) .=> u0) - end - if p === missing - p = Dict() - end - if t0 === nothing - t0 = 0.0 + if oldinitprob.f.resid_prototype === nothing + newf = oldinitprob.f + else + newf = NonlinearFunction{ + SciMLBase.isinplace(oldinitprob.f), SciMLBase.specialization(oldinitprob.f)}( + oldinitprob.f; + resid_prototype = calculate_resid_prototype( + length(oldinitprob.f.resid_prototype), newu0, newp)) + end + initprob = remake(oldinitprob; f = newf, u0 = newu0, p = newp) + return SciMLBase.OverrideInitData(initprob, odefn.update_initializeprob!, + odefn.initializeprobmap, odefn.initializeprobpmap) end - u0 = todict(u0) + dvs = unknowns(sys) + ps = parameters(sys) + u0map = to_varmap(u0, dvs) + symbols_to_symbolics!(sys, u0map) + pmap = to_varmap(p, ps) + symbols_to_symbolics!(sys, pmap) + guesses = Dict() defs = defaults(sys) - varmap = merge(defs, u0) - for k in collect(keys(varmap)) - if varmap[k] === nothing - delete!(varmap, k) + cmap, cs = get_cmap(sys) + use_scc = true + + if SciMLBase.has_initializeprob(odefn) + oldsys = odefn.initializeprob.f.sys + meta = get_metadata(oldsys) + if meta isa InitializationSystemMetadata + u0map = merge(meta.u0map, u0map) + pmap = merge(meta.pmap, pmap) + merge!(guesses, meta.additional_guesses) + use_scc = get(meta.extra_metadata, :use_scc, true) end - end - varmap = canonicalize_varmap(varmap) - missingvars = setdiff(unknowns(sys), collect(keys(varmap))) - setobserved = filter(keys(varmap)) do var - has_observed_with_lhs(sys, var) || has_observed_with_lhs(sys, default_toterm(var)) - end - p = todict(p) - guesses = ModelingToolkit.guesses(sys) - solvablepars = [par - for par in parameters(sys) - if is_parameter_solvable(par, p, defs, guesses)] - pvarmap = merge(defs, p) - setparobserved = filter(keys(pvarmap)) do var - has_parameter_dependency_with_lhs(sys, var) - end - if (((!isempty(missingvars) || !isempty(solvablepars) || - !isempty(setobserved) || !isempty(setparobserved)) && - ModelingToolkit.get_tearing_state(sys) !== nothing) || - !isempty(initialization_equations(sys))) - if SciMLBase.has_initializeprob(odefn) - oldsys = odefn.initializeprob.f.sys - meta = get_metadata(oldsys) - if meta isa InitializationSystemMetadata - u0 = merge(meta.u0map, u0) - p = merge(meta.pmap, p) + else + # there is no initializeprob, so the original problem construction + # had no solvable parameters and had the differential variables + # specified in `u0map`. + if u0 === missing + # the user didn't pass `u0` to `remake`, so they want to retain + # existing values. Fill the differential variables in `u0map`, + # initialization will either be elided or solve for the algebraic + # variables + diff_idxs = isdiffeq.(equations(sys)) + for i in eachindex(dvs) + diff_idxs[i] || continue + u0map[dvs[i]] = newu0[i] end end - for k in collect(keys(u0)) - if u0[k] === nothing - delete!(u0, k) + if p === missing + # the user didn't pass `p` to `remake`, so they want to retain + # existing values. Fill all parameters in `pmap` so that none of + # them are solvable. + for p in ps + pmap[p] = getp(sys, p)(newp) end end - for k in collect(keys(p)) - if p[k] === nothing - delete!(p, k) - end + # all non-solvable parameters need values regardless + for p in ps + haskey(pmap, p) && continue + is_parameter_solvable(p, pmap, defs, guesses) && continue + pmap[p] = getp(sys, p)(newp) end + end + if t0 === nothing + t0 = 0.0 + end + filter_missing_values!(u0map) + filter_missing_values!(pmap) - initprob = InitializationProblem(sys, t0, u0, p) - initprobmap = getu(initprob, unknowns(sys)) - punknowns = [p for p in all_variable_symbols(initprob) if is_parameter(sys, p)] - getpunknowns = getu(initprob, punknowns) - setpunknowns = setp(sys, punknowns) - initprobpmap = GetUpdatedMTKParameters(getpunknowns, setpunknowns) - reqd_syms = parameter_symbols(initprob) - update_initializeprob! = UpdateInitializeprob( - getu(sys, reqd_syms), setu(initprob, reqd_syms)) - return initprob, update_initializeprob!, initprobmap, initprobpmap - else - return nothing, nothing, nothing, nothing + op, missing_unknowns, missing_pars = build_operating_point( + u0map, pmap, defs, cmap, dvs, ps) + kws = maybe_build_initialization_problem( + sys, op, u0map, pmap, t0, defs, guesses, missing_unknowns; use_scc) + initprob = get(kws, :initializeprob, nothing) + if initprob === nothing + return nothing end + return SciMLBase.OverrideInitData(initprob, get(kws, :update_initializeprob!, nothing), + get(kws, :initializeprobmap, nothing), + get(kws, :initializeprobpmap, nothing)) end """ diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index 94b44e0a6f..3cb68853aa 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -283,6 +283,16 @@ function hessian_sparsity(sys::NonlinearSystem) unknowns(sys)) for eq in equations(sys)] end +function calculate_resid_prototype(N, u0, p) + u0ElType = u0 === nothing ? Float64 : eltype(u0) + if SciMLStructures.isscimlstructure(p) + u0ElType = promote_type( + eltype(SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]), + u0ElType) + end + return zeros(u0ElType, N) +end + """ ```julia SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(sys), @@ -337,13 +347,7 @@ function SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(s if length(dvs) == length(equations(sys)) resid_prototype = nothing else - u0ElType = u0 === nothing ? Float64 : eltype(u0) - if SciMLStructures.isscimlstructure(p) - u0ElType = promote_type( - eltype(SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]), - u0ElType) - end - resid_prototype = zeros(u0ElType, length(equations(sys))) + resid_prototype = calculate_resid_prototype(length(equations(sys)), u0, p) end NonlinearFunction{iip}(f, @@ -492,10 +496,14 @@ end function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem, u0map, parammap = DiffEqBase.NullParameters(); - check_length = true, kwargs...) where {iip} + check_length = true, use_homotopy_continuation = true, kwargs...) where {iip} if !iscomplete(sys) error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `NonlinearProblem`") end + prob = safe_HomotopyContinuationProblem(sys, u0map, parammap; check_length, kwargs...) + if prob isa HomotopyContinuationProblem + return prob + end f, u0, p = process_SciMLProblem(NonlinearFunction{iip}, sys, u0map, parammap; check_length, kwargs...) pt = something(get_metadata(sys), StandardNonlinearProblem()) @@ -531,6 +539,186 @@ function DiffEqBase.NonlinearLeastSquaresProblem{iip}(sys::NonlinearSystem, u0ma NonlinearLeastSquaresProblem{iip}(f, u0, p; filter_kwargs(kwargs)...) end +struct CacheWriter{F} + fn::F +end + +function (cw::CacheWriter)(p, sols) + cw.fn(p.caches[1], sols, p...) +end + +function CacheWriter(sys::AbstractSystem, exprs, solsyms, obseqs::Vector{Equation}; + eval_expression = false, eval_module = @__MODULE__) + ps = parameters(sys) + rps = reorder_parameters(sys, ps) + obs_assigns = [eq.lhs ← eq.rhs for eq in obseqs] + cmap, cs = get_cmap(sys) + cmap_assigns = [eq.lhs ← eq.rhs for eq in cmap] + fn = Func( + [:out, DestructuredArgs(DestructuredArgs.(solsyms)), + DestructuredArgs.(rps)...], + [], + SetArray(true, :out, exprs) + ) |> wrap_assignments(false, obs_assigns)[2] |> + wrap_parameter_dependencies(sys, false)[2] |> + wrap_array_vars(sys, exprs; dvs = nothing, inputs = [])[2] |> + wrap_assignments(false, cmap_assigns)[2] |> toexpr + return CacheWriter(eval_or_rgf(fn; eval_expression, eval_module)) +end + +struct SCCNonlinearFunction{iip} end + +function SCCNonlinearFunction{iip}( + sys::NonlinearSystem, _eqs, _dvs, _obs, cachesyms; eval_expression = false, + eval_module = @__MODULE__, kwargs...) where {iip} + ps = parameters(sys) + rps = reorder_parameters(sys, ps) + + obs_assignments = [eq.lhs ← eq.rhs for eq in _obs] + + cmap, cs = get_cmap(sys) + cmap_assignments = [eq.lhs ← eq.rhs for eq in cmap] + rhss = [eq.rhs - eq.lhs for eq in _eqs] + wrap_code = wrap_assignments(false, cmap_assignments) .∘ + (wrap_array_vars(sys, rhss; dvs = _dvs, cachesyms)) .∘ + wrap_parameter_dependencies(sys, false) .∘ + wrap_assignments(false, obs_assignments) + f_gen = build_function( + rhss, _dvs, rps..., cachesyms...; wrap_code, expression = Val{true}) + f_oop, f_iip = eval_or_rgf.(f_gen; eval_expression, eval_module) + + f(u, p) = f_oop(u, p) + f(u, p::MTKParameters) = f_oop(u, p...) + f(resid, u, p) = f_iip(resid, u, p) + f(resid, u, p::MTKParameters) = f_iip(resid, u, p...) + + subsys = NonlinearSystem(_eqs, _dvs, ps; observed = _obs, + parameter_dependencies = parameter_dependencies(sys), name = nameof(sys)) + if get_index_cache(sys) !== nothing + @set! subsys.index_cache = subset_unknowns_observed( + get_index_cache(sys), sys, _dvs, getproperty.(_obs, (:lhs,))) + @set! subsys.complete = true + end + + return NonlinearFunction{iip}(f; sys = subsys) +end + +function SciMLBase.SCCNonlinearProblem(sys::NonlinearSystem, args...; kwargs...) + SCCNonlinearProblem{true}(sys, args...; kwargs...) +end + +function SciMLBase.SCCNonlinearProblem{iip}(sys::NonlinearSystem, u0map, + parammap = SciMLBase.NullParameters(); eval_expression = false, eval_module = @__MODULE__, kwargs...) where {iip} + if !iscomplete(sys) || get_tearing_state(sys) === nothing + error("A simplified `NonlinearSystem` is required. Call `structural_simplify` on the system before creating an `SCCNonlinearProblem`.") + end + + if !is_split(sys) + error("The system has been simplified with `split = false`. `SCCNonlinearProblem` is not compatible with this system. Pass `split = true` to `structural_simplify` to use `SCCNonlinearProblem`.") + end + + ts = get_tearing_state(sys) + var_eq_matching, var_sccs = StructuralTransformations.algebraic_variables_scc(ts) + + if length(var_sccs) == 1 + return NonlinearProblem{iip}( + sys, u0map, parammap; eval_expression, eval_module, kwargs...) + end + + condensed_graph = MatchedCondensationGraph( + DiCMOBiGraph{true}(complete(ts.structure.graph), + complete(var_eq_matching)), + var_sccs) + toporder = topological_sort_by_dfs(condensed_graph) + var_sccs = var_sccs[toporder] + eq_sccs = map(Base.Fix1(getindex, var_eq_matching), var_sccs) + + dvs = unknowns(sys) + ps = parameters(sys) + eqs = equations(sys) + obs = observed(sys) + + _, u0, p = process_SciMLProblem( + EmptySciMLFunction, sys, u0map, parammap; eval_expression, eval_module, kwargs...) + + explicitfuns = [] + nlfuns = [] + prevobsidxs = Int[] + cachesize = 0 + for (i, (escc, vscc)) in enumerate(zip(eq_sccs, var_sccs)) + # subset unknowns and equations + _dvs = dvs[vscc] + _eqs = eqs[escc] + # get observed equations required by this SCC + obsidxs = observed_equations_used_by(sys, _eqs) + # the ones used by previous SCCs can be precomputed into the cache + setdiff!(obsidxs, prevobsidxs) + _obs = obs[obsidxs] + + # get all subexpressions in the RHS which we can precompute in the cache + banned_vars = Set{Any}(vcat(_dvs, getproperty.(_obs, (:lhs,)))) + for var in banned_vars + iscall(var) || continue + operation(var) === getindex || continue + push!(banned_vars, arguments(var)[1]) + end + state = Dict() + for i in eachindex(_obs) + _obs[i] = _obs[i].lhs ~ subexpressions_not_involving_vars!( + _obs[i].rhs, banned_vars, state) + end + for i in eachindex(_eqs) + _eqs[i] = _eqs[i].lhs ~ subexpressions_not_involving_vars!( + _eqs[i].rhs, banned_vars, state) + end + + # cached variables and their corresponding expressions + cachevars = Any[obs[i].lhs for i in prevobsidxs] + cacheexprs = Any[obs[i].lhs for i in prevobsidxs] + for (k, v) in state + push!(cachevars, unwrap(v)) + push!(cacheexprs, unwrap(k)) + end + cachesize = max(cachesize, length(cachevars)) + + if isempty(cachevars) + push!(explicitfuns, Returns(nothing)) + else + solsyms = getindex.((dvs,), view(var_sccs, 1:(i - 1))) + push!(explicitfuns, + CacheWriter(sys, cacheexprs, solsyms, obs[prevobsidxs]; + eval_expression, eval_module)) + end + f = SCCNonlinearFunction{iip}( + sys, _eqs, _dvs, _obs, (cachevars,); eval_expression, eval_module, kwargs...) + push!(nlfuns, f) + append!(cachevars, _dvs) + append!(cacheexprs, _dvs) + for i in obsidxs + push!(cachevars, obs[i].lhs) + push!(cacheexprs, obs[i].rhs) + end + append!(prevobsidxs, obsidxs) + end + + if cachesize != 0 + p = rebuild_with_caches(p, BufferTemplate(eltype(u0), cachesize)) + end + + subprobs = [] + for (f, vscc) in zip(nlfuns, var_sccs) + prob = NonlinearProblem(f, u0[vscc], p) + push!(subprobs, prob) + end + + new_dvs = dvs[reduce(vcat, var_sccs)] + new_eqs = eqs[reduce(vcat, eq_sccs)] + @set! sys.unknowns = new_dvs + @set! sys.eqs = new_eqs + sys = complete(sys) + return SCCNonlinearProblem(subprobs, explicitfuns, p, true; sys) +end + """ $(TYPEDSIGNATURES) @@ -682,72 +870,3 @@ function Base.:(==)(sys1::NonlinearSystem, sys2::NonlinearSystem) _eq_unordered(get_ps(sys1), get_ps(sys2)) && all(s1 == s2 for (s1, s2) in zip(get_systems(sys1), get_systems(sys2))) end - -""" -$(TYPEDEF) - -A type of Nonlinear problem which specializes on polynomial systems and uses -HomotopyContinuation.jl to solve the system. Requires importing HomotopyContinuation.jl to -create and solve. -""" -struct HomotopyContinuationProblem{uType, H, D, O, SS, U} <: - SciMLBase.AbstractNonlinearProblem{uType, true} - """ - The initial values of states in the system. If there are multiple real roots of - the system, the one closest to this point is returned. - """ - u0::uType - """ - A subtype of `HomotopyContinuation.AbstractSystem` to solve. Also contains the - parameter object. - """ - homotopy_continuation_system::H - """ - A function with signature `(u, p) -> resid`. In case of rational functions, this - is used to rule out roots of the system which would cause the denominator to be - zero. - """ - denominator::D - """ - The `NonlinearSystem` used to create this problem. Used for symbolic indexing. - """ - sys::NonlinearSystem - """ - A function which generates and returns observed expressions for the given system. - """ - obsfn::O - """ - The HomotopyContinuation.jl solver and start system, obtained through - `HomotopyContinuation.solver_startsystems`. - """ - solver_and_starts::SS - """ - A function which takes a solution of the transformed system, and returns a vector - of solutions for the original system. This is utilized when converting systems - to polynomials. - """ - unpack_solution::U -end - -function HomotopyContinuationProblem(::AbstractSystem, _u0, _p; kwargs...) - error("HomotopyContinuation.jl is required to create and solve `HomotopyContinuationProblem`s. Please run `Pkg.add(\"HomotopyContinuation\")` to continue.") -end - -SymbolicIndexingInterface.symbolic_container(p::HomotopyContinuationProblem) = p.sys -SymbolicIndexingInterface.state_values(p::HomotopyContinuationProblem) = p.u0 -function SymbolicIndexingInterface.set_state!(p::HomotopyContinuationProblem, args...) - set_state!(p.u0, args...) -end -function SymbolicIndexingInterface.parameter_values(p::HomotopyContinuationProblem) - parameter_values(p.homotopy_continuation_system) -end -function SymbolicIndexingInterface.set_parameter!(p::HomotopyContinuationProblem, args...) - set_parameter!(parameter_values(p), args...) -end -function SymbolicIndexingInterface.observed(p::HomotopyContinuationProblem, sym) - if p.obsfn !== nothing - return p.obsfn(sym) - else - return SymbolicIndexingInterface.observed(p.sys, sym) - end -end diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl index 3425216b5f..43e9294dd3 100644 --- a/src/systems/optimization/optimizationsystem.jl +++ b/src/systems/optimization/optimizationsystem.jl @@ -168,7 +168,8 @@ function OptimizationSystem(objective; constraints = [], kwargs...) push!(new_ps, p) end end - return OptimizationSystem(objective, collect(allunknowns), collect(new_ps); constraints, kwargs...) + return OptimizationSystem( + objective, collect(allunknowns), collect(new_ps); constraints, kwargs...) end function flatten(sys::OptimizationSystem) diff --git a/src/systems/parameter_buffer.jl b/src/systems/parameter_buffer.jl index 72af9f594d..7d64054acc 100644 --- a/src/systems/parameter_buffer.jl +++ b/src/systems/parameter_buffer.jl @@ -3,11 +3,12 @@ symconvert(::Type{T}, x) where {T} = convert(T, x) symconvert(::Type{Real}, x::Integer) = convert(Float64, x) symconvert(::Type{V}, x) where {V <: AbstractArray} = convert(V, symconvert.(eltype(V), x)) -struct MTKParameters{T, D, C, N} +struct MTKParameters{T, D, C, N, H} tunable::T discrete::D constant::C nonnumeric::N + caches::H end """ @@ -181,11 +182,18 @@ function MTKParameters( mtkps = MTKParameters{ typeof(tunable_buffer), typeof(disc_buffer), typeof(const_buffer), - typeof(nonnumeric_buffer)}(tunable_buffer, disc_buffer, const_buffer, - nonnumeric_buffer) + typeof(nonnumeric_buffer), typeof(())}(tunable_buffer, + disc_buffer, const_buffer, nonnumeric_buffer, ()) return mtkps end +function rebuild_with_caches(p::MTKParameters, cache_templates::BufferTemplate...) + buffers = map(cache_templates) do template + Vector{template.type}(undef, template.length) + end + @set p.caches = buffers +end + function narrow_buffer_type(buffer::AbstractArray) type = Union{} for x in buffer @@ -297,7 +305,8 @@ end for (Portion, field, recurse) in [(SciMLStructures.Discrete, :discrete, 1) (SciMLStructures.Constants, :constant, 1) - (Nonnumeric, :nonnumeric, 1)] + (Nonnumeric, :nonnumeric, 1) + (SciMLStructures.Caches, :caches, 1)] @eval function SciMLStructures.canonicalize(::$Portion, p::MTKParameters) as_vector = buffer_to_arraypartition(p.$field) repack = let p = p @@ -324,11 +333,13 @@ function Base.copy(p::MTKParameters) discrete = Tuple(eltype(buf) <: Real ? copy(buf) : copy.(buf) for buf in p.discrete) constant = Tuple(eltype(buf) <: Real ? copy(buf) : copy.(buf) for buf in p.constant) nonnumeric = copy.(p.nonnumeric) + caches = copy.(p.caches) return MTKParameters( tunable, discrete, constant, - nonnumeric + nonnumeric, + caches ) end @@ -640,7 +651,7 @@ end # getindex indexes the vectors, setindex! linearly indexes values # it's inconsistent, but we need it to be this way @generated function Base.getindex( - ps::MTKParameters{T, D, C, N}, idx::Int) where {T, D, C, N} + ps::MTKParameters{T, D, C, N, H}, idx::Int) where {T, D, C, N, H} paths = [] if !(T <: SizedVector{0, Float64}) push!(paths, :(ps.tunable)) @@ -654,6 +665,9 @@ end for i in 1:fieldcount(N) push!(paths, :(ps.nonnumeric[$i])) end + for i in 1:fieldcount(H) + push!(paths, :(ps.caches[$i])) + end expr = Expr(:if, :(idx == 1), :(return $(paths[1]))) curexpr = expr for i in 2:length(paths) @@ -663,12 +677,12 @@ end return Expr(:block, expr, :(throw(BoundsError(ps, idx)))) end -@generated function Base.length(ps::MTKParameters{T, D, C, N}) where {T, D, C, N} +@generated function Base.length(ps::MTKParameters{T, D, C, N, H}) where {T, D, C, N, H} len = 0 if !(T <: SizedVector{0, Float64}) len += 1 end - len += fieldcount(D) + fieldcount(C) + fieldcount(N) + len += fieldcount(D) + fieldcount(C) + fieldcount(N) + fieldcount(H) return len end @@ -691,7 +705,10 @@ end function Base.:(==)(a::MTKParameters, b::MTKParameters) return a.tunable == b.tunable && a.discrete == b.discrete && - a.constant == b.constant && a.nonnumeric == b.nonnumeric + a.constant == b.constant && a.nonnumeric == b.nonnumeric && + all(Iterators.map(a.caches, b.caches) do acache, bcache + eltype(acache) == eltype(bcache) && length(acache) == length(bcache) + end) end # to support linearize/linearization_function diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 1cf110df13..6b75fb2c6d 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -4,12 +4,13 @@ const AnyDict = Dict{Any, Any} $(TYPEDSIGNATURES) If called without arguments, return `Dict{Any, Any}`. Otherwise, interpret the input -as a symbolic map and turn it into a `Dict{Any, Any}`. Handles `SciMLBase.NullParameters` -and `nothing`. +as a symbolic map and turn it into a `Dict{Any, Any}`. Handles `SciMLBase.NullParameters`, +`missing` and `nothing`. """ anydict() = AnyDict() anydict(::SciMLBase.NullParameters) = AnyDict() anydict(::Nothing) = AnyDict() +anydict(::Missing) = AnyDict() anydict(x::AnyDict) = x anydict(x) = AnyDict(x) @@ -51,6 +52,42 @@ function add_toterms(varmap::AbstractDict; toterm = default_toterm) return cp end +""" + $(TYPEDSIGNATURES) + +Turn any `Symbol` keys in `varmap` to the appropriate symbolic variables in `sys`. Any +symbols that cannot be converted are ignored. +""" +function symbols_to_symbolics!(sys::AbstractSystem, varmap::AbstractDict) + if is_split(sys) + ic = get_index_cache(sys) + for k in collect(keys(varmap)) + k isa Symbol || continue + newk = get(ic.symbol_to_variable, k, nothing) + newk === nothing && continue + varmap[newk] = varmap[k] + delete!(varmap, k) + end + else + syms = all_symbols(sys) + for k in collect(keys(varmap)) + k isa Symbol || continue + idx = findfirst(syms) do sym + hasname(sym) || return false + name = getname(sym) + return name == k + end + idx === nothing && continue + newk = syms[idx] + if iscall(newk) && operation(newk) === getindex + newk = arguments(newk)[1] + end + varmap[newk] = varmap[k] + delete!(varmap, k) + end + end +end + """ $(TYPEDSIGNATURES) @@ -388,6 +425,15 @@ function evaluate_varmap!(varmap::AbstractDict, vars; limit = 100) end end +""" + $(TYPEDSIGNATURES) + +Remove keys in `varmap` whose values are `nothing`. +""" +function filter_missing_values!(varmap::AbstractDict) + filter!(kvp -> kvp[2] !== nothing, varmap) +end + struct GetUpdatedMTKParameters{G, S} # `getu` functor which gets parameters that are unknowns during initialization getpunknowns::G @@ -431,12 +477,99 @@ end $(TYPEDEF) A simple utility meant to be used as the `constructor` passed to `process_SciMLProblem` in -case constructing a SciMLFunction is not required. +case constructing a SciMLFunction is not required. The arguments passed to it are available +in the `args` field, and the keyword arguments in the `kwargs` field. """ -struct EmptySciMLFunction end +struct EmptySciMLFunction{A, K} + args::A + kwargs::K +end function EmptySciMLFunction(args...; kwargs...) - return nothing + return EmptySciMLFunction{typeof(args), typeof(kwargs)}(args, kwargs) +end + +""" + $(TYPEDSIGNATURES) + +Construct the operating point of the system from the user-provided `u0map` and `pmap`, system +defaults `defs`, constant equations `cmap` (from `get_cmap(sys)`), unknowns `dvs` and +parameters `ps`. Return the operating point as a dictionary, the list of unknowns for which +no values can be determined, and the list of parameters for which no values can be determined. +""" +function build_operating_point( + u0map::AbstractDict, pmap::AbstractDict, defs::AbstractDict, cmap, dvs, ps) + op = add_toterms(u0map) + missing_unknowns = add_fallbacks!(op, dvs, defs) + for (k, v) in defs + haskey(op, k) && continue + op[k] = v + end + merge!(op, pmap) + missing_pars = add_fallbacks!(op, ps, defs) + for eq in cmap + op[eq.lhs] = eq.rhs + end + return op, missing_unknowns, missing_pars +end + +""" + $(TYPEDSIGNATURES) + +Build and return the initialization problem and associated data as a `NamedTuple` to be passed +to the `SciMLFunction` constructor. Requires the system `sys`, operating point `op`, +user-provided `u0map` and `pmap`, initial time `t`, system defaults `defs`, user-provided +`guesses`, and list of unknowns which don't have a value in `op`. The keyword `implicit_dae` +denotes whether the `SciMLProblem` being constructed is in implicit DAE form (`DAEProblem`). +All other keyword arguments are forwarded to `InitializationProblem`. +""" +function maybe_build_initialization_problem( + sys::AbstractSystem, op::AbstractDict, u0map, pmap, t, defs, + guesses, missing_unknowns; implicit_dae = false, kwargs...) + guesses = merge(ModelingToolkit.guesses(sys), todict(guesses)) + has_observed_u0s = any( + k -> has_observed_with_lhs(sys, k) || has_parameter_dependency_with_lhs(sys, k), + keys(op)) + solvablepars = [p + for p in parameters(sys) + if is_parameter_solvable(p, pmap, defs, guesses)] + has_dependent_unknowns = any(unknowns(sys)) do sym + val = get(op, sym, nothing) + val === nothing && return false + return symbolic_type(val) != NotSymbolic() || is_array_of_symbolics(val) + end + if (((implicit_dae || has_observed_u0s || !isempty(missing_unknowns) || + !isempty(solvablepars) || has_dependent_unknowns) && + get_tearing_state(sys) !== nothing) || + !isempty(initialization_equations(sys))) && t !== nothing + initializeprob = ModelingToolkit.InitializationProblem( + sys, t, u0map, pmap; guesses, kwargs...) + initializeprobmap = getu(initializeprob, unknowns(sys)) + + punknowns = [p + for p in all_variable_symbols(initializeprob) + if is_parameter(sys, p)] + getpunknowns = getu(initializeprob, punknowns) + setpunknowns = setp(sys, punknowns) + initializeprobpmap = GetUpdatedMTKParameters(getpunknowns, setpunknowns) + + reqd_syms = parameter_symbols(initializeprob) + update_initializeprob! = UpdateInitializeprob( + getu(sys, reqd_syms), setu(initializeprob, reqd_syms)) + for p in punknowns + p = unwrap(p) + stype = symtype(p) + op[p] = get_temporary_value(p) + end + + for v in missing_unknowns + op[v] = zero_var(v) + end + empty!(missing_unknowns) + return (; + initializeprob, initializeprobmap, initializeprobpmap, update_initializeprob!) + end + return (;) end """ @@ -491,6 +624,8 @@ Keyword arguments: Only applicable if `warn_cyclic_dependency == true`. - `substitution_limit`: The number times to substitute initial conditions into each other to attempt to arrive at a numeric value. +- `use_scc`: Whether to use `SCCNonlinearProblem` for initialization if the system is fully + determined. All other keyword arguments are passed as-is to `constructor`. """ @@ -504,7 +639,7 @@ function process_SciMLProblem( symbolic_u0 = false, warn_cyclic_dependency = false, circular_dependency_max_cycle_length = length(all_symbols(sys)), circular_dependency_max_cycles = 10, - substitution_limit = 100, kwargs...) + substitution_limit = 100, use_scc = true, kwargs...) dvs = unknowns(sys) ps = parameters(sys) iv = has_iv(sys) ? get_iv(sys) : nothing @@ -516,73 +651,26 @@ function process_SciMLProblem( pType = typeof(pmap) _u0map = u0map u0map = to_varmap(u0map, dvs) + symbols_to_symbolics!(sys, u0map) _pmap = pmap pmap = to_varmap(pmap, ps) + symbols_to_symbolics!(sys, pmap) defs = add_toterms(recursive_unwrap(defaults(sys))) cmap, cs = get_cmap(sys) kwargs = NamedTuple(kwargs) - op = add_toterms(u0map) - missing_unknowns = add_fallbacks!(op, dvs, defs) - for (k, v) in defs - haskey(op, k) && continue - op[k] = v - end - merge!(op, pmap) - missing_pars = add_fallbacks!(op, ps, defs) - for eq in cmap - op[eq.lhs] = eq.rhs - end - if sys isa ODESystem - guesses = merge(ModelingToolkit.guesses(sys), todict(guesses)) - has_observed_u0s = any( - k -> has_observed_with_lhs(sys, k) || has_parameter_dependency_with_lhs(sys, k), - keys(op)) - solvablepars = [p - for p in parameters(sys) - if is_parameter_solvable(p, pmap, defs, guesses)] - has_dependent_unknowns = any(unknowns(sys)) do sym - val = get(op, sym, nothing) - val === nothing && return false - return symbolic_type(val) != NotSymbolic() || is_array_of_symbolics(val) - end - if build_initializeprob && - (((implicit_dae || has_observed_u0s || !isempty(missing_unknowns) || - !isempty(solvablepars) || has_dependent_unknowns) && - get_tearing_state(sys) !== nothing) || - !isempty(initialization_equations(sys))) && t !== nothing - initializeprob = ModelingToolkit.InitializationProblem( - sys, t, u0map, pmap; guesses, warn_initialize_determined, - initialization_eqs, eval_expression, eval_module, fully_determined, - warn_cyclic_dependency, check_units = check_initialization_units, - circular_dependency_max_cycle_length, circular_dependency_max_cycles) - initializeprobmap = getu(initializeprob, unknowns(sys)) - - punknowns = [p - for p in all_variable_symbols(initializeprob) - if is_parameter(sys, p)] - getpunknowns = getu(initializeprob, punknowns) - setpunknowns = setp(sys, punknowns) - initializeprobpmap = GetUpdatedMTKParameters(getpunknowns, setpunknowns) - - reqd_syms = parameter_symbols(initializeprob) - update_initializeprob! = UpdateInitializeprob( - getu(sys, reqd_syms), setu(initializeprob, reqd_syms)) - for p in punknowns - p = unwrap(p) - stype = symtype(p) - op[p] = get_temporary_value(p) - delete!(missing_pars, p) - end + op, missing_unknowns, missing_pars = build_operating_point( + u0map, pmap, defs, cmap, dvs, ps) - for v in missing_unknowns - op[v] = zero_var(v) - end - empty!(missing_unknowns) - kwargs = merge(kwargs, - (; initializeprob, initializeprobmap, - initializeprobpmap, update_initializeprob!)) - end + if sys isa ODESystem && build_initializeprob + kws = maybe_build_initialization_problem( + sys, op, u0map, pmap, t, defs, guesses, missing_unknowns; + implicit_dae, warn_initialize_determined, initialization_eqs, + eval_expression, eval_module, fully_determined, + warn_cyclic_dependency, check_units = check_initialization_units, + circular_dependency_max_cycle_length, circular_dependency_max_cycles, use_scc) + + kwargs = merge(kwargs, kws) end if t !== nothing && !(constructor <: Union{DDEFunction, SDDEFunction}) diff --git a/src/systems/systems.jl b/src/systems/systems.jl index a54206d1dd..47acd81a82 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -26,7 +26,7 @@ topological sort of the observed equations in `sys`. + `fully_determined=true` controls whether or not an error will be thrown if the number of equations don't match the number of inputs, outputs, and equations. """ function structural_simplify( - sys::AbstractSystem, io = nothing; simplify = false, split = true, + sys::AbstractSystem, io = nothing; additional_passes = [], simplify = false, split = true, allow_symbolic = false, allow_parameter = true, conservative = false, fully_determined = true, kwargs...) isscheduled(sys) && throw(RepeatedStructuralSimplificationError()) @@ -46,6 +46,9 @@ function structural_simplify( not yet supported. """) end + for pass in additional_passes + newsys = pass(newsys) + end if newsys isa ODESystem || has_parent(newsys) @set! newsys.parent = complete(sys; split, flatten = false) end @@ -72,6 +75,10 @@ function __structural_simplify(sys::JumpSystem, args...; kwargs...) return sys end +function __structural_simplify(sys::SDESystem, args...; kwargs...) + return __structural_simplify(ODESystem(sys), args...; kwargs...) +end + function __structural_simplify(sys::AbstractSystem, io = nothing; simplify = false, kwargs...) sys = expand_connections(sys) @@ -154,8 +161,9 @@ function __structural_simplify(sys::AbstractSystem, io = nothing; simplify = fal end noise_eqs = StructuralTransformations.tearing_substitute_expr(ode_sys, noise_eqs) - return SDESystem(full_equations(ode_sys), noise_eqs, + return SDESystem(Vector{Equation}(full_equations(ode_sys)), noise_eqs, get_iv(ode_sys), unknowns(ode_sys), parameters(ode_sys); - name = nameof(ode_sys), is_scalar_noise, observed = observed(ode_sys)) + name = nameof(ode_sys), is_scalar_noise, observed = observed(ode_sys), defaults = defaults(sys), + parameter_dependencies = parameter_dependencies(sys)) end end diff --git a/src/utils.jl b/src/utils.jl index 1555cd624e..dd7113cbe6 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -569,7 +569,6 @@ function collect_vars!(unknowns, parameters, p::Pair, iv; depth = 0, op = Differ return nothing end - function collect_var!(unknowns, parameters, var, iv; depth = 0) isequal(var, iv) && return nothing check_scope_depth(getmetadata(var, SymScope, LocalScope()), depth) || return nothing @@ -1002,9 +1001,143 @@ end diff2term_with_unit(x, t) = _with_unit(diff2term, x, t) lower_varname_with_unit(var, iv, order) = _with_unit(lower_varname, var, iv, iv, order) +""" + $(TYPEDSIGNATURES) + +Check if `sym` represents a symbolic floating point number or array of such numbers. +""" function is_variable_floatingpoint(sym) sym = unwrap(sym) T = symtype(sym) return T == Real || T <: AbstractFloat || T <: AbstractArray{Real} || T <: AbstractArray{<:AbstractFloat} end + +""" + $(TYPEDSIGNATURES) + +Return the `DiCMOBiGraph` denoting the dependencies between observed equations `eqs`. +""" +function observed_dependency_graph(eqs::Vector{Equation}) + for eq in eqs + if symbolic_type(eq.lhs) == NotSymbolic() + error("All equations must be observed equations of the form `var ~ expr`. Got $eq") + end + end + graph, assigns = observed2graph(eqs, getproperty.(eqs, (:lhs,))) + matching = complete(Matching(Vector{Union{Unassigned, Int}}(assigns))) + return DiCMOBiGraph{false}(graph, matching) +end + +""" + $(TYPEDSIGNATURES) + +Return the indexes of observed equations of `sys` used by expression `exprs`. +""" +function observed_equations_used_by(sys::AbstractSystem, exprs) + obs = observed(sys) + + obsvars = getproperty.(obs, :lhs) + graph = observed_dependency_graph(obs) + + syms = vars(exprs) + + obsidxs = BitSet() + for sym in syms + idx = findfirst(isequal(sym), obsvars) + idx === nothing && continue + parents = dfs_parents(graph, idx) + for i in eachindex(parents) + parents[i] == 0 && continue + push!(obsidxs, i) + end + end + + obsidxs = collect(obsidxs) + sort!(obsidxs) + return obsidxs +end + +""" + $(TYPEDSIGNATURES) + +Given an expression `expr`, return a dictionary mapping subexpressions of `expr` that do +not involve variables in `vars` to anonymous symbolic variables. Also return the modified +`expr` with the substitutions indicated by the dictionary. If `expr` is a function +of only `vars`, then all of the returned subexpressions can be precomputed. + +Note that this will only process subexpressions floating point value. Additionally, +array variables must be passed in both scalarized and non-scalarized forms in `vars`. +""" +function subexpressions_not_involving_vars(expr, vars) + expr = unwrap(expr) + vars = map(unwrap, vars) + state = Dict() + newexpr = subexpressions_not_involving_vars!(expr, vars, state) + return state, newexpr +end + +""" + $(TYPEDSIGNATURES) + +Mutating version of `subexpressions_not_involving_vars` which writes to `state`. Only +returns the modified `expr`. +""" +function subexpressions_not_involving_vars!(expr, vars, state::Dict{Any, Any}) + expr = unwrap(expr) + symbolic_type(expr) == NotSymbolic() && return expr + iscall(expr) || return expr + is_variable_floatingpoint(expr) || return expr + symtype(expr) <: Union{Real, AbstractArray{<:Real}} || return expr + Symbolics.shape(expr) == Symbolics.Unknown() && return expr + haskey(state, expr) && return state[expr] + vs = ModelingToolkit.vars(expr) + intersect!(vs, vars) + if isempty(vs) + sym = gensym(:subexpr) + stype = symtype(expr) + var = similar_variable(expr, sym) + state[expr] = var + return var + end + op = operation(expr) + args = arguments(expr) + if (op == (+) || op == (*)) && symbolic_type(expr) !== ArraySymbolic() + indep_args = [] + dep_args = [] + for arg in args + _vs = ModelingToolkit.vars(arg) + intersect!(_vs, vars) + if !isempty(_vs) + push!(dep_args, subexpressions_not_involving_vars!(arg, vars, state)) + else + push!(indep_args, arg) + end + end + indep_term = reduce(op, indep_args; init = Int(op == (*))) + indep_term = subexpressions_not_involving_vars!(indep_term, vars, state) + dep_term = reduce(op, dep_args; init = Int(op == (*))) + return op(indep_term, dep_term) + end + newargs = map(args) do arg + symbolic_type(arg) != NotSymbolic() || is_array_of_symbolics(arg) || return arg + subexpressions_not_involving_vars!(arg, vars, state) + end + return maketerm(typeof(expr), op, newargs, metadata(expr)) +end + +""" + $(TYPEDSIGNATURES) + +Create an anonymous symbolic variable of the same shape, size and symtype as `var`, with +name `gensym(name)`. Does not support unsized array symbolics. +""" +function similar_variable(var::BasicSymbolic, name = :anon) + name = gensym(name) + stype = symtype(var) + sym = Symbolics.variable(name; T = stype) + if size(var) !== () + sym = setmetadata(sym, Symbolics.ArrayShapeCtx, map(Base.OneTo, size(var))) + end + return sym +end diff --git a/src/variables.jl b/src/variables.jl index d11c6c1834..3b0d64e7ea 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -489,7 +489,9 @@ $(SIGNATURES) Define one or more Brownian variables. """ macro brownian(xs...) - all(x -> x isa Symbol || Meta.isexpr(x, :call) && x.args[1] == :$, xs) || + all( + x -> x isa Symbol || Meta.isexpr(x, :call) && x.args[1] == :$ || Meta.isexpr(x, :$), + xs) || error("@brownian only takes scalar expressions!") Symbolics._parse_vars(:brownian, Real, diff --git a/test/dde.jl b/test/dde.jl index 2030a90d06..c7561e6c24 100644 --- a/test/dde.jl +++ b/test/dde.jl @@ -76,12 +76,13 @@ prob = SDDEProblem(hayes_modelf, hayes_modelg, [1.0], h, tspan, pmul; constant_lags = (pmul[1],)); sol = solve(prob, RKMil(), seed = 100) -@variables x(..) +@variables x(..) delx(t) @parameters a=-4.0 b=-2.0 c=10.0 α=-1.3 β=-1.2 γ=1.1 @brownian η τ = 1.0 -eqs = [D(x(t)) ~ a * x(t) + b * x(t - τ) + c + (α * x(t) + γ) * η] +eqs = [D(x(t)) ~ a * x(t) + b * x(t - τ) + c + (α * x(t) + γ) * η, delx ~ x(t - τ)] @mtkbuild sys = System(eqs, t) +@test ModelingToolkit.has_observed_with_lhs(sys, delx) @test ModelingToolkit.is_dde(sys) @test !is_markovian(sys) @test equations(sys) == [D(x(t)) ~ a * x(t) + b * x(t - τ) + c] diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index 81c252c84a..554f9e1e1d 100644 --- a/test/extensions/homotopy_continuation.jl +++ b/test/extensions/homotopy_continuation.jl @@ -1,6 +1,18 @@ using ModelingToolkit, NonlinearSolve, SymbolicIndexingInterface +import ModelingToolkit as MTK using LinearAlgebra using Test + +@testset "Safe HCProblem" begin + @variables x y z + eqs = [0 ~ x^2 + y^2 + 2x * y + 0 ~ x^2 + 4x + 4 + 0 ~ y * z + 4x^2] + @mtkbuild sys = NonlinearSystem(eqs) + prob = MTK.safe_HomotopyContinuationProblem(sys, [x => 1.0, y => 1.0, z => 1.0], []) + @test prob === nothing +end + import HomotopyContinuation @testset "No parameters" begin @@ -9,12 +21,19 @@ import HomotopyContinuation 0 ~ x^2 + 4x + 4 0 ~ y * z + 4x^2] @mtkbuild sys = NonlinearSystem(eqs) - prob = HomotopyContinuationProblem(sys, [x => 1.0, y => 1.0, z => 1.0], []) + u0 = [x => 1.0, y => 1.0, z => 1.0] + prob = HomotopyContinuationProblem(sys, u0) @test prob[x] == prob[y] == prob[z] == 1.0 @test prob[x + y] == 2.0 sol = solve(prob; threading = false) @test SciMLBase.successful_retcode(sol) @test norm(sol.resid)≈0.0 atol=1e-10 + + prob2 = NonlinearProblem(sys, u0) + @test prob2 isa HomotopyContinuationProblem + sol = solve(prob2; threading = false) + @test SciMLBase.successful_retcode(sol) + @test norm(sol.resid)≈0.0 atol=1e-10 end struct Wrapper @@ -78,30 +97,45 @@ end @test_throws ["Cannot convert", "Unable", "symbolically solve", "Exponent", "not an integer", "not a polynomial"] HomotopyContinuationProblem( sys, []) + @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError + @test NonlinearProblem(sys, []) isa NonlinearProblem + @mtkbuild sys = NonlinearSystem([x^x - x ~ 0]) @test_throws ["Cannot convert", "Unable", "symbolically solve", "Exponent", "unknowns", "not a polynomial"] HomotopyContinuationProblem( sys, []) + @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError + @test NonlinearProblem(sys, []) isa NonlinearProblem @mtkbuild sys = NonlinearSystem([((x^2) / sin(x))^2 + x ~ 0]) @test_throws ["Cannot convert", "both polynomial", "non-polynomial", "recognized", "sin", "not a polynomial"] HomotopyContinuationProblem( sys, []) + @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError + @test NonlinearProblem(sys, []) isa NonlinearProblem @variables y = 2.0 @mtkbuild sys = NonlinearSystem([x^2 + y^2 + 2 ~ 0, y ~ sin(x)]) @test_throws ["Cannot convert", "recognized", "sin", "not a polynomial"] HomotopyContinuationProblem( sys, []) + @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError + @test NonlinearProblem(sys, []) isa NonlinearProblem @mtkbuild sys = NonlinearSystem([x^2 + y^2 - 2 ~ 0, sin(x + y) ~ 0]) @test_throws ["Cannot convert", "function of multiple unknowns"] HomotopyContinuationProblem( sys, []) + @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError + @test NonlinearProblem(sys, []) isa NonlinearProblem @mtkbuild sys = NonlinearSystem([sin(x)^2 + 1 ~ 0, cos(y) - cos(x) - 1 ~ 0]) @test_throws ["Cannot convert", "multiple non-polynomial terms", "same unknown"] HomotopyContinuationProblem( sys, []) + @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError + @test NonlinearProblem(sys, []) isa NonlinearProblem @mtkbuild sys = NonlinearSystem([sin(x^2)^2 + sin(x^2) - 1 ~ 0]) @test_throws ["import Nemo"] HomotopyContinuationProblem(sys, []) + @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError + @test NonlinearProblem(sys, []) isa NonlinearProblem end import Nemo @@ -126,6 +160,9 @@ end @test prob[x] ≈ 0.25 @test prob[y] ≈ 0.125 sol = solve(prob; threading = false) + # can't replicate the solve failure locally, so CI logs might help + @show sol.u sol.original.path_results + @test SciMLBase.successful_retcode(sol) @test sol[a]≈0.5 atol=1e-6 @test sol[b]≈0.25 atol=1e-6 end diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index 770f9f12bd..9c06dd2030 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -28,7 +28,7 @@ sol = solve(initprob) initprob = ModelingToolkit.InitializationProblem(pend, 0.0, [x => 1, y => 0], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) -@test initprob isa NonlinearProblem +@test initprob isa NonlinearLeastSquaresProblem sol = solve(initprob) @test SciMLBase.successful_retcode(sol) @test sol.u == [0.0, 0.0, 0.0, 0.0] @@ -975,3 +975,77 @@ end @test integ.ps[p] ≈ 1.0 @test integ.ps[q]≈cbrt(2) rtol=1e-6 end + +@testset "Guesses provided to `ODEProblem` are used in `remake`" begin + @variables x(t) y(t)=2x + @parameters p q=3x + @mtkbuild sys = ODESystem([D(x) ~ x * p + q, x^3 + y^3 ~ 3], t) + prob = ODEProblem( + sys, [], (0.0, 1.0), [p => 1.0]; guesses = [x => 1.0, y => 1.0, q => 1.0]) + @test prob[x] == 0.0 + @test prob[y] == 0.0 + @test prob.ps[p] == 1.0 + @test prob.ps[q] == 0.0 + integ = init(prob) + @test integ[x] ≈ 1 / cbrt(3) + @test integ[y] ≈ 2 / cbrt(3) + @test integ.ps[p] == 1.0 + @test integ.ps[q] ≈ 3 / cbrt(3) + prob2 = remake(prob; u0 = [y => 3x], p = [q => 2x]) + integ2 = init(prob2) + @test integ2[x] ≈ cbrt(3 / 28) + @test integ2[y] ≈ 3cbrt(3 / 28) + @test integ2.ps[p] == 1.0 + @test integ2.ps[q] ≈ 2cbrt(3 / 28) +end + +@testset "Remake problem with no initializeprob" begin + @variables x(t) [guess = 1.0] y(t) [guess = 1.0] + @parameters p [guess = 1.0] q [guess = 1.0] + @mtkbuild sys = ODESystem( + [D(x) ~ p * x + q * y, y ~ 2x], t; parameter_dependencies = [q ~ 2p]) + prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [p => 1.0]) + @test prob.f.initialization_data === nothing + prob2 = remake(prob; u0 = [x => 2.0]) + @test prob2[x] == 2.0 + @test prob2.f.initialization_data === nothing + prob3 = remake(prob; u0 = [y => 2.0]) + @test prob3.f.initialization_data !== nothing + @test init(prob3)[x] ≈ 1.0 + prob4 = remake(prob; p = [p => 1.0]) + @test prob4.f.initialization_data === nothing + prob5 = remake(prob; p = [p => missing, q => 2.0]) + @test prob5.f.initialization_data !== nothing + @test init(prob5).ps[p] ≈ 1.0 +end + +@testset "Variables provided as symbols" begin + @variables x(t) [guess = 1.0] y(t) [guess = 1.0] + @parameters p [guess = 1.0] q [guess = 1.0] + @mtkbuild sys = ODESystem( + [D(x) ~ p * x + q * y, y ~ 2x], t; parameter_dependencies = [q ~ 2p]) + prob = ODEProblem(sys, [:x => 1.0], (0.0, 1.0), [p => 1.0]) + @test prob.f.initialization_data === nothing + prob2 = remake(prob; u0 = [:x => 2.0]) + @test prob2.f.initialization_data === nothing + prob3 = remake(prob; u0 = [:y => 1.0]) + @test prob3.f.initialization_data !== nothing + @test init(prob3)[x] ≈ 0.5 +end + +@testset "Issue#3246: type promotion with parameter dependent initialization_eqs" begin + @variables x(t)=1 y(t)=1 + @parameters a = 1 + @named sys = ODESystem([D(x) ~ 0, D(y) ~ x + a], t; initialization_eqs = [y ~ a]) + + ssys = structural_simplify(sys) + prob = ODEProblem(ssys, [], (0, 1), []) + + @test SciMLBase.successful_retcode(solve(prob)) + + seta = setsym_oop(prob, [a]) + (newu0, newp) = seta(prob, ForwardDiff.Dual{ForwardDiff.Tag{:tag, Float64}}.([1.0], 1)) + newprob = remake(prob, u0 = newu0, p = newp) + + @test SciMLBase.successful_retcode(solve(newprob)) +end diff --git a/test/mtkparameters.jl b/test/mtkparameters.jl index 603426aaf7..ce524fdb76 100644 --- a/test/mtkparameters.jl +++ b/test/mtkparameters.jl @@ -299,7 +299,7 @@ end # Parameter timeseries ps = MTKParameters(([1.0, 1.0],), (BlockedArray(zeros(4), [2, 2]),), - (), ()) + (), (), ()) ps2 = SciMLStructures.replace(Discrete(), ps, ones(4)) @test typeof(ps2.discrete) == typeof(ps.discrete) with_updated_parameter_timeseries_values( @@ -316,7 +316,7 @@ with_updated_parameter_timeseries_values( ps = MTKParameters( (), (BlockedArray([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], [3, 3]), BlockedArray(falses(1), [1, 0])), - (), ()) + (), (), ()) @test SciMLBase.get_saveable_values(sys, ps, 1).x isa Tuple{Vector{Float64}, Vector{Bool}} tsidx1 = 1 tsidx2 = 2 diff --git a/test/runtests.jl b/test/runtests.jl index eaa87e6407..d240aaafa9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -78,6 +78,7 @@ end @safetestset "SDESystem Test" include("sdesystem.jl") @safetestset "DDESystem Test" include("dde.jl") @safetestset "NonlinearSystem Test" include("nonlinearsystem.jl") + @safetestset "SCCNonlinearProblem Test" include("scc_nonlinear_problem.jl") @safetestset "PDE Construction Test" include("pde.jl") @safetestset "JumpSystem Test" include("jumpsystem.jl") @safetestset "BVProblem Test" include("bvproblem.jl") diff --git a/test/scc_nonlinear_problem.jl b/test/scc_nonlinear_problem.jl new file mode 100644 index 0000000000..fdf1646343 --- /dev/null +++ b/test/scc_nonlinear_problem.jl @@ -0,0 +1,163 @@ +using ModelingToolkit +using NonlinearSolve, SCCNonlinearSolve +using OrdinaryDiffEq +using SciMLBase, Symbolics +using LinearAlgebra, Test +using ModelingToolkit: t_nounits as t, D_nounits as D + +@testset "Trivial case" begin + function f!(du, u, p) + du[1] = cos(u[2]) - u[1] + du[2] = sin(u[1] + u[2]) + u[2] + du[3] = 2u[4] + u[3] + 1.0 + du[4] = u[5]^2 + u[4] + du[5] = u[3]^2 + u[5] + du[6] = u[1] + u[2] + u[3] + u[4] + u[5] + 2.0u[6] + 2.5u[7] + 1.5u[8] + du[7] = u[1] + u[2] + u[3] + 2.0u[4] + u[5] + 4.0u[6] - 1.5u[7] + 1.5u[8] + du[8] = u[1] + 2.0u[2] + 3.0u[3] + 5.0u[4] + 6.0u[5] + u[6] - u[7] - u[8] + end + @variables u[1:8] [irreducible = true] + eqs = Any[0 for _ in 1:8] + f!(eqs, u, nothing) + eqs = 0 .~ eqs + @named model = NonlinearSystem(eqs) + @test_throws ["simplified", "required"] SCCNonlinearProblem(model, []) + _model = structural_simplify(model; split = false) + @test_throws ["not compatible"] SCCNonlinearProblem(_model, []) + model = structural_simplify(model) + prob = NonlinearProblem(model, [u => zeros(8)]) + sccprob = SCCNonlinearProblem(model, [u => zeros(8)]) + sol1 = solve(prob, NewtonRaphson()) + sol2 = solve(sccprob, NewtonRaphson()) + @test SciMLBase.successful_retcode(sol1) + @test SciMLBase.successful_retcode(sol2) + @test sol1[u] ≈ sol2[u] +end + +@testset "With parameters" begin + function f!(du, u, (p1, p2), t) + x = (*)(p1[4], u[1]) + y = (*)(p1[4], (+)(0.1016, (*)(-1, u[1]))) + z1 = ifelse((<)(p2[1], 0), + (*)((*)(457896.07999999996, p1[2]), sqrt((*)(1.1686468413521012e-5, p1[3]))), + 0) + z2 = ifelse((>)(p2[1], 0), + (*)((*)((*)(0.58, p1[2]), sqrt((*)(1 // 86100, p1[3]))), u[4]), + 0) + z3 = ifelse((>)(p2[1], 0), + (*)((*)(457896.07999999996, p1[2]), sqrt((*)(1.1686468413521012e-5, p1[3]))), + 0) + z4 = ifelse((<)(p2[1], 0), + (*)((*)((*)(0.58, p1[2]), sqrt((*)(1 // 86100, p1[3]))), u[5]), + 0) + du[1] = p2[1] + du[2] = (+)(z1, (*)(-1, z2)) + du[3] = (+)(z3, (*)(-1, z4)) + du[4] = (+)((*)(-1, u[2]), (*)((*)(1 // 86100, y), u[4])) + du[5] = (+)((*)(-1, u[3]), (*)((*)(1 // 86100, x), u[5])) + end + p = ( + [0.04864391799335977, 7.853981633974484e-5, 1.4034843205574914, + 0.018241469247509915, 300237.05, 9.226186337232914], + [0.0508]) + u0 = [0.0, 0.0, 0.0, 789476.0, 101325.0] + tspan = (0.0, 1.0) + mass_matrix = [1.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0] + dt = 1e-3 + function nlf(u1, (u0, p)) + resid = Any[0 for _ in u0] + f!(resid, u1, p, 0.0) + return mass_matrix * (u1 - u0) - dt * resid + end + + prob = NonlinearProblem(nlf, u0, (u0, p)) + @test_throws Exception solve(prob, SimpleNewtonRaphson(), abstol = 1e-9) + sol = solve(prob, TrustRegion(); abstol = 1e-9) + + @variables u[1:5] [irreducible = true] + @parameters p1[1:6] p2 + eqs = 0 .~ collect(nlf(u, (u0, (p1, p2)))) + @mtkbuild sys = NonlinearSystem(eqs, [u], [p1, p2]) + sccprob = SCCNonlinearProblem(sys, [u => u0], [p1 => p[1], p2 => p[2][]]) + sccsol = solve(sccprob, SimpleNewtonRaphson(); abstol = 1e-9) + @test SciMLBase.successful_retcode(sccsol) + @test norm(sccsol.resid) < norm(sol.resid) +end + +@testset "Transistor amplifier" begin + C = [k * 1e-6 for k in 1:5] + Ub = 6 + UF = 0.026 + α = 0.99 + β = 1e-6 + R0 = 1000 + R = 9000 + Ue(t) = 0.1 * sin(200 * π * t) + + function transamp(out, du, u, p, t) + g(x) = 1e-6 * (exp(x / 0.026) - 1) + y1, y2, y3, y4, y5, y6, y7, y8 = u + out[1] = -Ue(t) / R0 + y1 / R0 + C[1] * du[1] - C[1] * du[2] + out[2] = -Ub / R + y2 * 2 / R - (α - 1) * g(y2 - y3) - C[1] * du[1] + C[1] * du[2] + out[3] = -g(y2 - y3) + y3 / R + C[2] * du[3] + out[4] = -Ub / R + y4 / R + α * g(y2 - y3) + C[3] * du[4] - C[3] * du[5] + out[5] = -Ub / R + y5 * 2 / R - (α - 1) * g(y5 - y6) - C[3] * du[4] + C[3] * du[5] + out[6] = -g(y5 - y6) + y6 / R + C[4] * du[6] + out[7] = -Ub / R + y7 / R + α * g(y5 - y6) + C[5] * du[7] - C[5] * du[8] + out[8] = y8 / R - C[5] * du[7] + C[5] * du[8] + end + + u0 = [0, Ub / 2, Ub / 2, Ub, Ub / 2, Ub / 2, Ub, 0] + du0 = [ + 51.338775, + 51.338775, + -Ub / (2 * (C[2] * R)), + -24.9757667, + -24.9757667, + -Ub / (2 * (C[4] * R)), + -10.00564453, + -10.00564453 + ] + daeprob = DAEProblem(transamp, du0, u0, (0.0, 0.1)) + daesol = solve(daeprob, DImplicitEuler()) + + t0 = daesol.t[5] + t1 = daesol.t[6] + u0 = daesol.u[5] + u1 = daesol.u[6] + dt = t1 - t0 + + @variables y(t)[1:8] + eqs = Any[0 for _ in 1:8] + transamp(eqs, collect(D(y)), y, nothing, t) + eqs = 0 .~ eqs + subrules = Dict(Symbolics.unwrap(D(y[i])) => ((y[i] - u0[i]) / dt) for i in 1:8) + eqs = substitute.(eqs, (subrules,)) + @mtkbuild sys = NonlinearSystem(eqs) + prob = NonlinearProblem(sys, [y => u0], [t => t0]) + sol = solve(prob, NewtonRaphson(); abstol = 1e-12) + + sccprob = SCCNonlinearProblem(sys, [y => u0], [t => t0]) + sccsol = solve(sccprob, NewtonRaphson(); abstol = 1e-12) + + @test sol.u≈sccsol.u atol=1e-10 +end + +@testset "Expression caching" begin + @variables x[1:4] = rand(4) + val = Ref(0) + function func(x, y) + val[] += 1 + x + y + end + @register_symbolic func(x, y) + @mtkbuild sys = NonlinearSystem([0 ~ x[1]^3 + x[2]^3 - 5 + 0 ~ sin(x[1] - x[2]) - 0.5 + 0 ~ func(x[1], x[2]) * exp(x[3]) - x[4]^3 - 5 + 0 ~ func(x[1], x[2]) * exp(x[4]) - x[3]^3 - 4]) + sccprob = SCCNonlinearProblem(sys, []) + sccsol = solve(sccprob, NewtonRaphson()) + @test SciMLBase.successful_retcode(sccsol) + @test val[] == 1 +end diff --git a/test/sdesystem.jl b/test/sdesystem.jl index 749aca86a7..8069581dcc 100644 --- a/test/sdesystem.jl +++ b/test/sdesystem.jl @@ -809,3 +809,62 @@ end prob = SDEProblem(sys, [x => 1.0, y => 1.0], (0.0, 1.0)) @test prob[z] ≈ 2.0 end + +@testset "SDESystem to ODESystem" begin + @variables x(t) y(t) z(t) + @testset "Scalar noise" begin + @named sys = SDESystem([D(x) ~ x, D(y) ~ y, z ~ x + y], [x, y, 3], + t, [x, y, z], [], is_scalar_noise = true) + odesys = ODESystem(sys) + @test odesys isa ODESystem + vs = ModelingToolkit.vars(equations(odesys)) + nbrownian = count( + v -> ModelingToolkit.getvariabletype(v) == ModelingToolkit.BROWNIAN, vs) + @test nbrownian == 3 + for eq in equations(odesys) + ModelingToolkit.isdiffeq(eq) || continue + @test length(arguments(eq.rhs)) == 4 + end + end + + @testset "Non-scalar vector noise" begin + @named sys = SDESystem([D(x) ~ x, D(y) ~ y, z ~ x + y], [x, y, 0], + t, [x, y, z], [], is_scalar_noise = false) + odesys = ODESystem(sys) + @test odesys isa ODESystem + vs = ModelingToolkit.vars(equations(odesys)) + nbrownian = count( + v -> ModelingToolkit.getvariabletype(v) == ModelingToolkit.BROWNIAN, vs) + @test nbrownian == 1 + for eq in equations(odesys) + ModelingToolkit.isdiffeq(eq) || continue + @test length(arguments(eq.rhs)) == 2 + end + end + + @testset "Matrix noise" begin + noiseeqs = [x+y y+z z+x + 2y 2z 2x + z+1 x+1 y+1] + @named sys = SDESystem([D(x) ~ x, D(y) ~ y, D(z) ~ z], noiseeqs, t, [x, y, z], []) + odesys = ODESystem(sys) + @test odesys isa ODESystem + vs = ModelingToolkit.vars(equations(odesys)) + nbrownian = count( + v -> ModelingToolkit.getvariabletype(v) == ModelingToolkit.BROWNIAN, vs) + @test nbrownian == 3 + for eq in equations(odesys) + @test length(arguments(eq.rhs)) == 4 + end + end +end + +@testset "`structural_simplify(::SDESystem)`" begin + @variables x(t) y(t) + @mtkbuild sys = SDESystem( + [D(x) ~ x, y ~ 2x], [x, 0], t, [x, y], []; is_scalar_noise = true) + @test sys isa SDESystem + @test length(equations(sys)) == 1 + @test length(ModelingToolkit.get_noiseeqs(sys)) == 1 + @test length(observed(sys)) == 1 +end diff --git a/test/split_parameters.jl b/test/split_parameters.jl index 22c90edf7a..2f8667faa8 100644 --- a/test/split_parameters.jl +++ b/test/split_parameters.jl @@ -206,7 +206,7 @@ S = get_sensitivity(closed_loop, :u) BlockedArray([[1 2; 3 4], [2 4; 6 8]], [1, 1])), # (BlockedArray([[true, false], [false, true]]), BlockedArray([[[1 2; 3 4]], [[2 4; 6 8]]])), ([5, 6],), - (["hi", "bye"], [:lie, :die])) + (["hi", "bye"], [:lie, :die]), ()) @test ps[ParameterIndex(Tunable(), 1)] == 1.0 @test ps[ParameterIndex(Tunable(), 2:4)] == collect(2.0:4.0) @test ps[ParameterIndex(Tunable(), reshape(4:7, 2, 2))] == reshape(4.0:7.0, 2, 2) diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl index 2704559f72..863e091aad 100644 --- a/test/structural_transformation/utils.jl +++ b/test/structural_transformation/utils.jl @@ -152,3 +152,12 @@ end end end end + +@testset "additional passes" begin + @variables x(t) y(t) + @named sys = ODESystem([D(x) ~ x, y ~ x + t], t) + value = Ref(0) + pass(sys; kwargs...) = (value[] += 1; return sys) + structural_simplify(sys; additional_passes = [pass]) + @test value[] == 1 +end diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 61690acdce..ccf0a17a40 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -8,6 +8,7 @@ using ModelingToolkit: SymbolicContinuousCallback, using StableRNGs import SciMLBase using SymbolicIndexingInterface +using Setfield rng = StableRNG(12345) @variables x(t) = 0 @@ -227,6 +228,119 @@ affect_neg = [x ~ 1] @test e[].affect == affect end +@testset "ImperativeAffect constructors" begin + fmfa(o, x, i, c) = nothing + m = ModelingToolkit.ImperativeAffect(fmfa) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test m.obs == [] + @test m.obs_syms == [] + @test m.modified == [] + @test m.mod_syms == [] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa, (;)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test m.obs == [] + @test m.obs_syms == [] + @test m.modified == [] + @test m.mod_syms == [] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa, (; x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, []) + @test m.obs_syms == [] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:x] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa, (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, []) + @test m.obs_syms == [] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa; observed = (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:y] + @test m.modified == [] + @test m.mod_syms == [] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa; modified = (; x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, []) + @test m.obs_syms == [] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:x] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa; modified = (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, []) + @test m.obs_syms == [] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa, (; x), (; x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:x] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:x] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect(fmfa, (; y = x), (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:y] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect( + fmfa; modified = (; y = x), observed = (; y = x)) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:y] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === nothing + + m = ModelingToolkit.ImperativeAffect( + fmfa; modified = (; y = x), observed = (; y = x), ctx = 3) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:y] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:y] + @test m.ctx === 3 + + m = ModelingToolkit.ImperativeAffect(fmfa, (; x), (; x), 3) + @test m isa ModelingToolkit.ImperativeAffect + @test m.f == fmfa + @test isequal(m.obs, [x]) + @test m.obs_syms == [:x] + @test isequal(m.modified, [x]) + @test m.mod_syms == [:x] + @test m.ctx === 3 +end + ## @named sys = ODESystem(eqs, t, continuous_events = [x ~ 1]) @@ -822,7 +936,7 @@ end @named trigsys = ODESystem(eqs, t; continuous_events = [evt1, evt2]) trigsys_ss = structural_simplify(trigsys) prob = ODEProblem(trigsys_ss, [], (0.0, 2π)) - sol = solve(prob, Tsit5()) + sol = solve(prob, Tsit5(); dtmax = 0.01) required_crossings_c1 = [π / 2, 3 * π / 2] required_crossings_c2 = [π / 6, π / 2, 5 * π / 6, 7 * π / 6, 3 * π / 2, 11 * π / 6] @test maximum(abs.(first.(cr1) .- required_crossings_c1)) < 1e-4 @@ -969,6 +1083,173 @@ end @test sol[b] == [2.0, 5.0, 5.0] @test sol[c] == [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0] end +@testset "Heater" begin + @variables temp(t) + params = @parameters furnace_on_threshold=0.5 furnace_off_threshold=0.7 furnace_power=1.0 leakage=0.1 furnace_on::Bool=false + eqs = [ + D(temp) ~ furnace_on * furnace_power - temp^2 * leakage + ] + + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c + @set! x.furnace_on = false + end) + furnace_enable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_on_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, i, c + @set! x.furnace_on = true + end) + @named sys = ODESystem( + eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) + ss = structural_simplify(sys) + prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + sol = solve(prob, Tsit5(); dtmax = 0.01) + @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) + + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i + @set! x.furnace_on = false + end; initialize = ModelingToolkit.ImperativeAffect(modified = (; + temp)) do x, o, c, i + @set! x.temp = 0.2 + end) + furnace_enable = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_on_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on)) do x, o, c, i + @set! x.furnace_on = true + end) + @named sys = ODESystem( + eqs, t, [temp], params; continuous_events = [furnace_off, furnace_enable]) + ss = structural_simplify(sys) + prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + sol = solve(prob, Tsit5(); dtmax = 0.01) + @test all(sol[temp][sol.t .> 1.0] .<= 0.79) && all(sol[temp][sol.t .> 1.0] .>= 0.49) + @test all(sol[temp][sol.t .!= 0.0] .<= 0.79) && all(sol[temp][sol.t .!= 0.0] .>= 0.2) +end + +@testset "ImperativeAffect errors and warnings" begin + @variables temp(t) + params = @parameters furnace_on_threshold=0.5 furnace_off_threshold=0.7 furnace_power=1.0 leakage=0.1 furnace_on::Bool=false + eqs = [ + D(temp) ~ furnace_on * furnace_power - temp^2 * leakage + ] + + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect( + modified = (; furnace_on), observed = (; furnace_on)) do x, o, c, i + @set! x.furnace_on = false + end) + @named sys = ODESystem(eqs, t, [temp], params; continuous_events = [furnace_off]) + ss = structural_simplify(sys) + @test_logs (:warn, + "The symbols Any[:furnace_on] are declared as both observed and modified; this is a code smell because it becomes easy to confuse them and assign/not assign a value.") prob=ODEProblem( + ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + + @variables tempsq(t) # trivially eliminated + eqs = [tempsq ~ temp^2 + D(temp) ~ furnace_on * furnace_power - temp^2 * leakage] + + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect( + modified = (; furnace_on, tempsq), observed = (; furnace_on)) do x, o, c, i + @set! x.furnace_on = false + end) + @named sys = ODESystem( + eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) + ss = structural_simplify(sys) + @test_throws "refers to missing variable(s)" prob=ODEProblem( + ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + + @parameters not_actually_here + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on), + observed = (; furnace_on, not_actually_here)) do x, o, c, i + @set! x.furnace_on = false + end) + @named sys = ODESystem( + eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) + ss = structural_simplify(sys) + @test_throws "refers to missing variable(s)" prob=ODEProblem( + ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + + furnace_off = ModelingToolkit.SymbolicContinuousCallback( + [temp ~ furnace_off_threshold], + ModelingToolkit.ImperativeAffect(modified = (; furnace_on), + observed = (; furnace_on)) do x, o, c, i + return (; fictional2 = false) + end) + @named sys = ODESystem( + eqs, t, [temp, tempsq], params; continuous_events = [furnace_off]) + ss = structural_simplify(sys) + prob = ODEProblem( + ss, [temp => 0.0, furnace_on => true], (0.0, 100.0)) + @test_throws "Tried to write back to" solve(prob, Tsit5()) +end + +@testset "Quadrature" begin + @variables theta(t) omega(t) + params = @parameters qA=0 qB=0 hA=0 hB=0 cnt::Int=0 + eqs = [D(theta) ~ omega + omega ~ 1.0] + function decoder(oldA, oldB, newA, newB) + state = (oldA, oldB, newA, newB) + if state == (0, 0, 1, 0) || state == (1, 0, 1, 1) || state == (1, 1, 0, 1) || + state == (0, 1, 0, 0) + return 1 + elseif state == (0, 0, 0, 1) || state == (0, 1, 1, 1) || state == (1, 1, 1, 0) || + state == (1, 0, 0, 0) + return -1 + elseif state == (0, 0, 0, 0) || state == (0, 1, 0, 1) || state == (1, 0, 1, 0) || + state == (1, 1, 1, 1) + return 0 + else + return 0 # err is interpreted as no movement + end + end + qAevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta) ~ 0], + ModelingToolkit.ImperativeAffect((; qA, hA, hB, cnt), (; qB)) do x, o, c, i + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 1 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x + end, + affect_neg = ModelingToolkit.ImperativeAffect( + (; qA, hA, hB, cnt), (; qB)) do x, o, c, i + @set! x.hA = x.qA + @set! x.hB = o.qB + @set! x.qA = 0 + @set! x.cnt += decoder(x.hA, x.hB, x.qA, o.qB) + x + end; rootfind = SciMLBase.RightRootFind) + qBevt = ModelingToolkit.SymbolicContinuousCallback([cos(100 * theta - π / 2) ~ 0], + ModelingToolkit.ImperativeAffect((; qB, hA, hB, cnt), (; qA)) do x, o, c, i + @set! x.hA = o.qA + @set! x.hB = x.qB + @set! x.qB = 1 + @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) + x + end, + affect_neg = ModelingToolkit.ImperativeAffect( + (; qB, hA, hB, cnt), (; qA)) do x, o, c, i + @set! x.hA = o.qA + @set! x.hB = x.qB + @set! x.qB = 0 + @set! x.cnt += decoder(x.hA, x.hB, o.qA, x.qB) + x + end; rootfind = SciMLBase.RightRootFind) + @named sys = ODESystem( + eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) + ss = structural_simplify(sys) + prob = ODEProblem(ss, [theta => 1e-5], (0.0, pi)) + sol = solve(prob, Tsit5(); dtmax = 0.01) + @test getp(sol, cnt)(sol) == 198 # we get 2 pulses per phase cycle (cos 0 crossing) and we go to 100 cycles; we miss a few due to the initial state +end @testset "Initialization" begin @variables x(t) @@ -1063,7 +1344,7 @@ end cb = [x ~ 0.0] => [x ~ 0, y ~ 1] @mtkbuild pend = ODESystem(eqs, t; continuous_events = [cb]) prob = ODEProblem(pend, [x => 1], (0.0, 3.0), guesses = [y => x]) - @test_throws "CheckInit specified but initialization" solve(prob, Rodas5()) + @test_throws "DAE initialization failed" solve(prob, Rodas5()) cb = [x ~ 0.0] => [y ~ 1] @mtkbuild pend = ODESystem(eqs, t; continuous_events = [cb])