`get_unit` not supported for callable parameter

[TorkelE]

Also means these cannot be used in Catalyst, as we automatically calls get_unit on all parameters when we perform unit checks.

Example:

using ModelingToolkit

ts = collect(0.0:0.1:10.0)
spline = LinearInterpolation(ts .^ 2, ts)
Tspline = typeof(spline)
@parameters (interp::Tspline)(..)
@parameters p

ModelingToolkit.get_unit(p) # 1.0
ModelingToolkit.get_unit(interp) # ERROR: ArgumentError: Unsupported value interp⋆.

[ChrisRackauckas]

It's not clear this is well defined? A function doesn't have parameters. It could have a parameter map, like tells you the parameters in and parameters out, but doesnt necessarily have parameters. I think Catalyst would need to specialize on this

[TorkelE]

Do you mean a function has units? Or a callable parameter have units?

I am following the tutorial here: https://docs.sciml.ai/ModelingToolkit/stable/tutorials/callable_params/ on how to create a parameter which has some defined values across time values. If you do

using OrdinaryDiffEq
using DataInterpolations
using ModelingToolkit: t_nounits as t, D_nounits as D

# Defines the input.
ts = collect(0.0:0.1:10.0)
spline = LinearInterpolation(ts .^ 2, ts)
Tspline = typeof(spline)
@parameters (interp::Tspline)(..)

# Defines the model.
@variables V(t)
@parameters p d
eqs = [D(V) ~ p - interp(t)*V]
@mtkbuild sys = ODESystem(eqs, t)

Then interp is stored as a parameter in parameters(sys). However, ModelingToolkit.get_unit(interp) fails. In this case, Catalyst performs unit checking for all parameters stored in ReactionSystems. But since this type of parameter is not supported by ModelingToolkit.get_unit (but stored as a parameter) we get an error.

Example where we would want to use this

Imagine we want to make a model where a protein is activated by sunlight (L):

rn = @reaction_network begin
    L, X --> Xa
end

Here, L would be a function of t. I could make this a sinus, but sometimes you want to use more realistic conditions. Then it is common to take data of sunlight from somewhere in the world. In Catalyst, we would want to encode this in a LinearInterpolation and feed it as a callable parameter:

rn = @reaction_network begin
    L(t), X --> Xa # Never mind how the values of `L(t)` is defined right now
end

I would like to make this work, but it is currently not possible.

Ofcourse I could have a check and then disable unit checking if one of these parameters is present. However, it would be rather messy. It also does seem like an omission that ModelingToolkit.get_unit doesn't work here (especially as we actually can assign units to these parameters

[ChrisRackauckas]

Do you mean a function has units? Or a callable parameter have units?

A callable parameter is a function. Functions don't have units.

[TorkelE]

I don't necessarily agree. In my example, the function describes the amount of light somewhere, and light totally has a unit. I'd imagine most places where we put these within MTK they would have some unit.

And if that is the case that they cannot have units, shouldn't

@parameters (interp::Tspline)(..) [unit = u"s"]

throw an error?

[ChrisRackauckas]

It would probably be good to make that throw an error, yes.

In my example, the function describes the amount of light somewhere

The output of the function can have a unit. That unit can be dependent on the units of the input. Take the function f(x)=x^2. If the input is units m then the units of f(x) is m^2. If x has s units then f(x) is s^2. Saying that the units of f are m^2 is simply not true: a function if it has units, that unit would be a function or a map of units.

[hersle]

Would it work to keep the function/spline itself dimensionless (i.e. convert to some fixed unit inside it and strip that unit away), and in the equations multiply up by that unit or a variable with that unit?

[TorkelE]

So I don't disagree with that argument per se, but it kind of obviously sidestep what I am talking about. Algebraic expressions cannot have units as as the unit of X^2+2 depends on the unit of X. Does that mean that we cannot check the units of MTK equations? Of course not because context matters.

Practically here these objects are not used as general Julia functions, and any cases where I can imagine them units makes total sense. But sure, maybe they have more extensive, undocumented, uses that I am aware of which changes things.

It does change quite a lot of code though, it a model can have parameters that cannot have units assigned.

[AayushSabharwal]

Does that mean that we cannot check the units of MTK equations? Of course not because context matters.

Context is not the difference here. We can check the units of MTK equations because we can inspect equations. i.e. if X has the unit m then we can infer that X^2 has a unit of m^2. Callable parameters are black boxes - we can't look inside and need a way for the user to specify the output unit(s) in terms of the input unit(s). If that same algebraic equation was replaced by f(x) where f is a callable parameter suddenly there's nothing we can infer.

Black box functions are hard. Even our current handling of unit verification through registered functions is problematic simply because we can't look inside.

Practically here these objects are not used as general Julia functions, and any cases where I can imagine them units makes total sense.

In my opinion the argument "any cases where I can imagine" is fallible - we can't be aware of all the ways someone uses tools. MTK especially is very generic and practically a compiler. Callable parameters don't have to be interpolations, and e.g. may be used simply as a replaceable aggregator function to swap out max for mean without rebuilding the model. We can't document every use case, only how things can be used.

As I see it, there are two ways for specifying the units for a callable parameter:

• A generic mapping of input units to output units. E.g. given m it returns m^2, given s it returns s^2 and so on.
• a rigid specification of input and output units. E.g. The input has to be s and the output will always be s^2.

Simply saying that "X function has Y unit" is ambiguous at best and incorrect at worst.

[hersle]

I want to add that it could be wise to not only consider units of callable parameters, but more generally of any registered function.

For example, with a spline/interpolation parameter, one might want both the value and derivative of the spline. Those would have different units if the time variable of the spline is dimensionful. Calling the parameter directly will only get you the value.

[AayushSabharwal]

I honestly think that for better unit support we need to lift this to Symbolics rather than implement it at the MTK level. @register_symbolic needs a begin..end block to add metadata (such as unit mapping) and Symbolics should validate units in arithmetic functions.

[ChrisRackauckas]

We do need that, where the unit map can be a function overload just like how derivatives are done. But I don't think that solves the problem here.

I think for this we need the ability to declare unit maps on function variables. That would be a tuple where wildcards could be used, like (x, x^2), or (x, x*m). Then the interpolation case would be either (s, m) like constant units, (0, m) where the input is unit less, or (x, m) where wildcards in are ignored and m units always come out. The last is weird but hey it's expressible. The constant case is then what interpolations and table lookups want.