Skip to content

Conversation

@bradcarman
Copy link
Contributor

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

Additional context

Closes issue #3823

@github-actions
Copy link
Contributor

github-actions bot commented Oct 28, 2025

Benchmark Results (Julia vlts)

Time benchmarks
master 9f04501... master / 9f04501...
ODEProblem 0.0919 ± 0.0068 s 0.0928 ± 0.0064 s 0.99 ± 0.1
init 0.0584 ± 0.0015 ms 0.0579 ± 0.0017 ms 1.01 ± 0.039
large_parameter_init/ODEProblem 0.982 ± 0.0063 s 0.974 ± 0.01 s 1.01 ± 0.013
large_parameter_init/init 0.108 ± 0.0022 ms 0.109 ± 0.0021 ms 0.99 ± 0.027
mtkcompile 0.0326 ± 0.00094 s 0.0338 ± 0.0011 s 0.965 ± 0.042
time_to_load 4.97 ± 0.036 s 5.07 ± 0.076 s 0.98 ± 0.016
Memory benchmarks
master 9f04501... master / 9f04501...
ODEProblem 0.667 M allocs: 22.4 MB 0.667 M allocs: 22.4 MB 1
init 0.892 k allocs: 0.0461 MB 0.892 k allocs: 0.0461 MB 1
large_parameter_init/ODEProblem 4.15 M allocs: 0.161 GB 4.15 M allocs: 0.161 GB 1
large_parameter_init/init 1.74 k allocs: 0.0857 MB 1.74 k allocs: 0.0857 MB 1
mtkcompile 0.225 M allocs: 8.34 MB 0.225 M allocs: 8.34 MB 1
time_to_load 0.153 k allocs: 14.5 kB 0.153 k allocs: 14.5 kB 1

@github-actions
Copy link
Contributor

github-actions bot commented Oct 28, 2025

Benchmark Results (Julia v1)

Time benchmarks
master 9f04501... master / 9f04501...
ODEProblem 0.0944 ± 0.0058 s 0.0981 ± 0.006 s 0.962 ± 0.083
init 0.0438 ± 0.011 ms 0.0443 ± 0.011 ms 0.988 ± 0.35
large_parameter_init/ODEProblem 1.15 ± 0.026 s 1.12 ± 0.035 s 1.02 ± 0.039
large_parameter_init/init 0.077 ± 0.018 ms 0.0772 ± 0.016 ms 0.997 ± 0.3
mtkcompile 0.033 ± 0.00075 s 0.0342 ± 0.00075 s 0.966 ± 0.031
time_to_load 4.59 ± 0.026 s 4.59 ± 0.015 s 0.999 ± 0.0065
Memory benchmarks
master 9f04501... master / 9f04501...
ODEProblem 0.681 M allocs: 22.5 MB 0.681 M allocs: 22.5 MB 1
init 0.797 k allocs: 30.6 kB 0.797 k allocs: 30.6 kB 1
large_parameter_init/ODEProblem 4.39 M allocs: 0.162 GB 4.39 M allocs: 0.162 GB 1
large_parameter_init/init 1.64 k allocs: 0.0596 MB 1.64 k allocs: 0.0596 MB 1
mtkcompile 0.23 M allocs: 7.88 MB 0.23 M allocs: 7.88 MB 1
time_to_load 0.149 k allocs: 11.1 kB 0.149 k allocs: 11.1 kB 1

Comment on lines 1 to 4
using SymbolicIndexingInterface
using Setfield
using StaticArrays
using OrdinaryDiffEqCore
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not top level

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fixed

Comment on lines +38 to +39
events::Tuple{SymbolicDiscreteCallback}
vars::Tuple{SymbolicUtils.BasicSymbolic{Real}}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not concrete

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fixed

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not fixed. Tuples are only concrete if you also choose a size for them

return sys
end

function DiffEqBase.solve(prob::SciMLBase.AbstractDEProblem, inputs::Vector{Input}, args...; kwargs...)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this violates the contract, and its not DiffEqBase

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fixed

Copy link
Member

@AayushSabharwal AayushSabharwal left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The added field isn't handled in namespacing. It looks like the intended behavior is for input_functions to only be populated on a simplified system? Consequently, I would recommend that build_input_functions should either be a model transformation in basic_transformations.jl or an optional simplification pass that can be included in additional_passes keyword of mtkcompile (similar to IfLifting).

Comment on lines +28 to +29
data::SVector
time::SVector
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why specifically SVector? Not only is this not concrete, I don't think it doesn't offer any performance gain given how it is used.

Comment on lines +33 to +34
n = length(data)
return Input(var, SVector{n}(data), SVector{n}(time))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is type-unstable. n isn't inferrable.

return Input(var, SVector{n}(data), SVector{n}(time))
end

struct InputFunctions
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do the fields of this struct need to be wrapped in Tuple{...}?


struct InputFunctions
events::Tuple{SymbolicDiscreteCallback}
vars::Tuple{SymbolicUtils.BasicSymbolic{Real}}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fairly minor, but vars is a bit of a misnomer for a field only containing one variable.

struct InputFunctions
events::Tuple{SymbolicDiscreteCallback}
vars::Tuple{SymbolicUtils.BasicSymbolic{Real}}
setters::Tuple{SymbolicIndexingInterface.ParameterHookWrapper}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

SII.ParameterHookWrapper is not public API.

return nothing
end
function set_input!(integrator, var, value::Real)
set_input!(get_input_functions(integrator.f.sys), integrator, var, value)
Copy link
Member

@AayushSabharwal AayushSabharwal Oct 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Accessing sys during a solve is discouraged. The solve should not involve symbolic quantities. This is a lot of unnecessary, type-unstable overhead. Solves should also work without the system - anything possible symbolically should be possible without the symbolic infrastructure.


# Here we ensure the inputs have metadata marking the discrete variables as parameters. In some
# cases the inputs can be fed to this function before they are converted to parameters by mtkcompile.
vars = SymbolicUtils.BasicSymbolic[isparameter(x) ? x : toparam(x)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think build_input_functions should be called as a model transformation after mtkcompile, and verify that the list of inputs it gets are inputs (and part of the list of parameters). This is more robust than trusting that the user passed values around correctly and trusting that there won't be bugs in the future.

end
end

@set! sys.discrete_events = events
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This would overwrite any pre-existing discrete events in the system.


# ensure that the ODEProblem does not complain about missing parameter map
if !haskey(defaults, x)
push!(defaults, x => 0.0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is mutating the defaults stored in sys. Model transformations should ideally be out of place.

end

function DiffEqBase.solve(prob::SciMLBase.AbstractDEProblem, inputs::Vector{Input}, args...; kwargs...)
tstops = Float64[]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's not guaranteed that the time is Float64.

@AayushSabharwal
Copy link
Member

I thought about this a bit. In my opinion, the best way to phrase this is a model transformation in basic_transformations.jl. It would operate only on simplified models with inputs. For every input, a parameter representing the input data is added. A callback is also generated for each input which sets the input at the appropriate time and sets the tstop to the next time the input value has to be set. Then simply generating and solving an ODEProblem like normal would replicate this. The generated parameters can be specially named using operators. So InputValues(x) would be a parameter for the data, and InputTimes(x) would be a parameter for the time points.


ModelingToolkit supports setting input values during simulation for variables marked with the `[input=true]` metadata. This is useful for real-time simulations, hardware-in-the-loop testing, interactive simulations, or any scenario where input values need to be determined during integration rather than specified beforehand.

To use this functionality, variables must be marked as inputs using the `[input=true]` metadata and specified in the `inputs` keyword argument of `@mtkcompile`.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

variables must be marked as inputs using the [input=true] metadata

is this actually required? The inputs are specified

in the inputs keyword argument of @mtkcompile

anyway. No other input-using functionality, like linearization or generation of input-dependent functions, rely on the input metadata

@bradcarman
Copy link
Contributor Author

I thought about this a bit. In my opinion, the best way to phrase this is a model transformation in basic_transformations.jl. It would operate only on simplified models with inputs. For every input, a parameter representing the input data is added. A callback is also generated for each input which sets the input at the appropriate time and sets the tstop to the next time the input value has to be set. Then simply generating and solving an ODEProblem like normal would replicate this. The generated parameters can be specially named using operators. So InputValues(x) would be a parameter for the data, and InputTimes(x) would be a parameter for the time points.

I really don't want to bake the input data into the System itself. The system should be able to be used in determinate or indeterminate forms. Putting the data into the System would work just for the determinate form, meaning we'd need a different System for the indeterminate form. We should be able to serialize/precompile a single System and ODEProblem, and then quickly run that system with different data sets. If the data gets put into the system as parameters, then each new data set needs to make (or remake) an ODEProblem.

@ChrisRackauckas
Copy link
Member

That's not true. If you do it with FunctionWrappers or DataInterpolations you can just swap the interpolation via remake. It would be type-stable and faster while also matching the normal semantics.

return sys
end

function CommonSolve.solve(prob::SciMLBase.AbstractDEProblem, inputs::Vector{Input}, args...; kwargs...)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It still violates the dispatch contract that SciML does. The second argument should be the algorithm, if we do this. But I don't think it makes sense to do this as a solve, it violates a few other principles of solve

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I want to keep the data fully separate from the System and the ODEProblem. We need some way to be able to solve a problem where data is fed in at the top level. I have benchmarks that show that this is much more efficient, and follows the norm for what other modeling platforms offer. If we move inputs as the 3rd argument, will that be acceptable?

Comment on lines +171 to +186
for input::Input in inputs
tstops = union(tstops, input.time)
condition = (u, t, integrator) -> any(t .== input.time)
affect! = function (integrator)
@inbounds begin
i = findfirst(integrator.t .== input.time)
set_input!(integrator, input.var, input.data[i])
end
end
push!(callbacks, DiscreteCallback(condition, affect!))

# DiscreteCallback doesn't hit on t==0, workaround...
if input.time[1] == 0
prob.ps[input.var] = input.data[1]
end
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This shouldn't be done here. This should instead be done by registering discontinuities for the interpolation function if it is a discontinuous one.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is no interpolation here, this is explicitly a discrete input. The point is to not have interpolation.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is ConstantInterpolation

@bradcarman
Copy link
Contributor Author

That's not true. If you do it with FunctionWrappers or DataInterpolations you can just swap the interpolation via remake. It would be type-stable and faster while also matching the normal semantics.

It's not faster, I have benchmarks that show this.

@ChrisRackauckas
Copy link
Member

No, you're not understanding... it's not faster you're just doing search reduction so you're not benchmarking the right thing.

Ultimately this just isn't the way to implement it. We just need a specialization on ConstantInterpolation that lowers it to a callback form to remove search. That needs none of the other extra interface breaking machinery here. This form is just impossible to maintain and breaks a lot of contracts, and most of it isn't necessary (and it's half broken in many scenarios 😅). That would also fix the type instabilities that are rampant in here (which end up breaking a lot of precompilation improvements)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants