Catalyst with units

[vchuravy]

I originally asked this on Slack, but I was looking for examples with Catalyst and units.

It wasn't immediately clear how to do this, so improved documentation would be appreciated.
Looking at some tests

I arrived at:

using Catalyst
using DynamicQuantities
using OrdinaryDiffEqDefault

@parameters t [unit=u"s"]
@species X(t) [unit=u"mol/m^3"]
@parameters d [unit=u"1/s"]
rxs = [
	Reaction(d, [X], [])
]
example = ReactionSystem(rxs, t, name=:example) |> complete

But now the question arises: How do I run this model.

My goal here is that u0 and p are also provided in units and that the solution also remembers the units (and the plotting and so forth)

tspan_example = (0.0, 14*24.0) .* u"s"
p_example = [
	:d => -1.0 * u"s^-1"
]
u0_example = [
	:X => 10 * u"mol/m^3"
]
ode_example = ODEProblem(example, u0_example, tspan_example, p_example)

fails with:

Cannot create an additive identity from `Type{<:DynamicQuantities.Quantity}`, as the dimensions are unknown. Please use `zero(::DynamicQuantities.Quantity)` instead.

Without units

tspan_example = (0.0, 14*60)
p_example = [
	:d => -1.0
]
u0_example = [
	:X => 10
]
ode_example = ODEProblem(example, u0_example, tspan_example, p_example)
sol_example = solve(ode_example, saveat=0.1)

does work.

function ustrip(varmap)
	map(varmap) do (name, value)
		name => DynamicQuantities.ustrip( ModelingToolkit.get_unit(name), value)
	end
end

function from_units(rn, u, tspan, p)
	u = symmap_to_varmap(rn, u)
	p = symmap_to_varmap(rn, p)
	tspan = map(tspan) do t
		DynamicQuantities.ustrip(ModelingToolkit.get_unit(rn.t), t)
	end
	
	ustrip(u), tspan, ustrip(p)
end

ode_example = ODEProblem(example, from_units(example, u0_example, tspan_example, p_example)...)

Allows me to do 50% of what I want, but I still need to add the units manually back on to the solution,
and of course ideally Catalyst would handle that for me.

julia> sol_example[:X] |> typeof
Vector{Float64} (alias for Array{Float64, 1})

[TorkelE]

So, I think the expected intended behaviour here (sorry, should have gone through properly in Slack) is that you only define your symbolci variables with units

@parameters t [unit=u"s"]
@species X(t) [unit=u"mol/m^3"]
@parameters d [unit=u"1/s"]

and your values normally (no units)

tspan_example = (0.0, 14*60)
p_example = [
	:d => -1.0
]
u0_example = [
	:X => 10
]

However, the extent with which carries through in different ways I am a little bit unsure of. I.e. I know ModelingToolkit (which does most of teh actual implementation) does unit checking to see that these add up in the equation. And I think it can add units to plots, but I don't actually think it is intended to store the units in the actual solution objects. However, I would need to check this with someone who knows more about units and ModelingToolkit.
(at some point I think we dropped the units from the docs because this was changing around and we weren't fully confident how it was meant to work around the switch from Unitful to DynamicQuantities)

[isaacsas]

To be honest, I'm not sure I'd trust the units stuff. It is basically all handled via ModelingToolkit, and I don't think is being developed much these days beyond keeping tests passing. I know there were issues with registered functions and units, and we had issues with certain types of expressions giving the wrong units. Also, since the conversion to DynamicQuantites I don't think standard chemistry units like Molars have worked right.

But to your specific question, I don't think you can attach units in the parameter and initial condition maps. See also Chris' recent comments about units in ModelingToolkit at:

https://discourse.julialang.org/t/mtk-units/133795/2

[isaacsas]

(Part of the reason we don't say much about units in the docs is due to the issues I mention -- I think we actually removed units as one of the Catalyst bullet point features at some point.)

[isaacsas]

@AayushSabharwal do you know if there is a way to reattach symbolic units to ODE solution objects for plotting and such?

[vchuravy]

and your values normally (no units)

So one of the reasons why units are interesting to me is that it makes for better input validation and output presentation.

I was showing someone Catalyst and one of the first questions was "what do the values correspond to"

[isaacsas]

I agree it is a nice feature, but unit handling is really all lower level than Catalyst (i.e. in ModelingToolkit for symbolic aspects and the base validation we reuse, and SciMLBase for plot recipes and solution indexing and such). I think the general conclusion @ChrisRackauckas came to over time is having units in the type is a bad idea within the actual SciML numerical solvers as it breaks lots of stuff (maybe numerical linear algebra was an example he ran into?). That said, it seems like the plot recipes could perhaps be modified in SciMLBase to handle automatically reattaching units from the symbolic metadata for ModelingToolkit generated models (or just labeling axes appropriately), and maybe convenience functions could be provided to reattach units to solution objects given the symbolic model the solution was generated from.

@ChrisRackauckas perhaps these might make a reasonable SciML small grant opportunity?

[vchuravy]

In particular Makie@0.24.7 (released a few days ago) comes with support for DynamicalQuantities and axis conversion.

https://docs.makie.org/dev/explanations/dim-converts#Experimental-DynamicQuantities.jl-support
MakieOrg/Makie.jl#5280

[vchuravy]

The DQ in Makie is kinda nice. I can set an axis to days and after manually reattaching units to the time axis it automatically adjusts the axis

[TorkelE]

That's actually rather neat