Consider a "auto-normalization-wrapper" in the future

[gregordecristoforo]

[gregordecristoforo]

A possible text case could be a vector with the fields:

a: 1.0 km
b: 1.0 s
c: 2.0 kg

Store these three values without the units in a single array (all values in this array need to be of same type, this means that the values are contiguously in memory which improves performance) and store units in a separate array.

Have a look for examples here: https://github.com/chaobrain/saiunit-experiments/tree/main/fig3-mutiscale-brain-network

[yaccos]

Primary objectives of this project:

• Provide a more efficient way to handle heterogeneous collections of Unitful quantities and other types when the underlying numerical type is identical
• Enable unit-checking for stiff ODE solvers without overhead. The units of measurements should protect the user from making mistakes when using the ODE solver or writing the ODE function, the correctness here is indenpendent of the ODE solver used. Hence, unit-checking should work no matter which ODE solver is chosen.

Top of the blackboard shows the current solution with HeterogeneousVector, the bottom shows the proposed solution using the auto-normalization wrapper:

Preliminary roadmap:

1. Implement structure named CollectionVector{T,S} (named can be changed later) which acts a view over a Vector{T}, which is similar to ValueShapes. T will be a concrete numerical type, whereas S will be a tuple defining the types components of the structure and their lengths. When an CollectionVector x with components a, b, and c is acessed, say x.a = 1.0 km, Julia will automatically return the correct Unitful value to the user. Julia can infer this type at compile-time, making no runtime overhead. However, when giving the structure to an ODE solver or a linear optimizer, the solver will only see the stored numbers and thus not have to care about the underlying units.
2. Implement wrappers around the ODE solver such that the user can provide the system with units of measurement through CollectionVector. Here, the wrapper will strip away the units of measurement (and convert if necessary) before passing the data to the ODE solver. A second wrapper around the ODE solver will re-attach the units of measurement before entering the user-defined ODE solver such that unit-checking can occur within the ODE function. This approach is supposed to work for explicit ODE solvers (eg. RK4).
3. Optimize the wrapper such that unnecessary conversions within the ODE function is avoided.
4. Enable support for implicit solvers (eg. trapezoid) with autodiff = false by providing a unit-aware version of the linsolve argument.
5. Implement automatic differentiation for CollectionVector with Unitful quantities. This may require implementing a two-dimensional structure to store the Jacobian matrix in a unit-aware fashion.
6. Modify the wrappers such that the implicit ODE solvers can use automatic differentiation instead of relying on finite differences to estimate the Jacobian.

[yaccos]

Prototype for new functions of the API:

using Unitful

x = [4.2u"m", 602.2u"kg/m^2/s", 1.42, 201.14u"kg^2*s/m^2", 603.3u"m",  6.45u"ft", 1.5u"lb/(ft)^2/hr"]

function get_dimensions(x::Number)
    return dimension(x)
end

function get_dimensions(x::AbstractArray)
    unique(get_dimensions.(x))
end

function get_units(::Number)
    return unit(x)
end

function get_units(x::Unitful.AbstractQuantity{T,D,U}) where {T,D,U}
    return unit(x)
end

function get_units(x::AbstractArray)
    dims = get_dimensions(x)
    res = (dim => Set{Unitful.FreeUnits}() for dim in dims) |> Dict
    for el in x
        push!(res[get_dimensions(el)], get_units(el))
    end
    return res
end

should print

> get_dimensions(x)
4-element Vector{Unitful.Dimensions}:
 𝐋
 𝐌 𝐋^-2 𝐓^-1
 NoDims
 𝐌^2 𝐓 𝐋^-2
> get_units(x)
Dict{Unitful.Dimensions, Set{Unitful.FreeUnits}} with 4 entries:
  𝐌 𝐋^-2 𝐓^-1 => Set([lb ft^-2 hr^-1, kg m^-2 s^-1])
  𝐋           => Set([m, ft])
  𝐌^2 𝐓 𝐋^-2  => Set([kg^2 s m^-2])
  NoDims      => Set([])

In the same spirit, I think of functions such as:

• get_units(::AbstractArray,::Unitful.Dimension): Returns the unique units of measure for one specific (composite) dimension of the array.
• unit_convert_array(x::AbstractArray, U::Unitful.Unit): Converts all elements of x matching the (composite) dimension of U into U.
• unit_convert_array(x::AbstractArray, D::Dict{Unitful.Dimensions, Unitful.Unit}): Converts all elements of x into the unit specified by its dimension.
• unit_normalize_array(x::AbstractArray, D::Dict{Unitful.Dimensions, Unitful.Unit}): Converts all elements of x into the unit specified by its dimension and strips away the units. Throws an error if there are dimensionful elements in x not covered by the keys in D.
• unit_base_convert_array(x::AbstractArray, U::Unitful.Unit): Converts all elements of x composed of the base dimension of U such that U is used to represent the said base dimension. Requires U to be a base unit.
• unit_base_convert_array(x::AbstractArray, U::Vector{Unitful.Unit}): Converts all elements of x composed of the base dimensions included in U such that the base units U is used to represent the said base dimensions. Requires all the elements of U to be disjoint base units.
• unit_base_normalize_array(x::AbstractArray, U::Vector{Unitful.Unit}): Converts all elements of x composed of the base dimensions included in U such that the base units U is used to represent the said base dimensions and strips away the units. Requires all the elements of U to be disjoint base units. Throws an error if there are base dimensions in x which are not represented by U.

[yaccos]

Also, consider the Buckingham pi theorem: https://en.wikipedia.org/wiki/Buckingham_pi_theorem

[yaccos]

Outline for the inner workings of the front normalization wrapper:

• Find the dimensions encoded in the system variables, simulation time and the model variables.
• Choose one base unit of measure for each of the base dimension encountered
• Determine a mapping between the units of measurement in the user-provided system and the units of measurements used internally based on the units of the base dimensions. We could consider adding some sort of scaling to avoid numeric artifacts if necessary.
• Convert all user-provided input into the internal system of units before stripping away the units when feeding the arguments into the ODE solver
• Consider how to do error handling, as the error message may not be informative to the user when the unit system is modified or units removed
• When the results come back from the ODE solver, convert the respective results for the system variables into their original units of measurement. The same applies to the simulation time.

For the back normalization wrapper, I envision three different modes of function:

• Mode A: The state and time are converted back into the original units of measurements before passed into the user-provided ODE function. Pros: Better compatibility with the user-provided function, internal printouts work as expected, allows for hard-coded parameters, the parameter struct does not need conversion. Cons: Require runtime overhead converting back and forth between units in the wrapper
• Mode B: The internal units of measurement are provided to the ODE solver. Pros: No runtime overhead in units conversion, while still retaining units checking. Cons: The units of the parameters must also be converted, greater care must be taken when writing the ODE function to allow any unit system of measurement to be used, constants cannot be hard-coded, more difficult to debug.
• Mode C: The ODE function only receives plain numbers directly from the ODE solver. Pros: The easiest to implement, requires no special structure to host the system variables (ComponentVector will work fine here), computationally efficient. Cons: Function must be written to support plain numerics, unit checking inside the function (a common mistake) does out the window.

I suggest we start implementing support for Mode A and continue to develop support for Mode B once the solution is mature enough. I don't think Mode C is useful enough for us to prioritize it, but we could keep it as an option in case it requires minimal effort to support it.

[yaccos]

Also, consider the Buckingham pi theorem: https://en.wikipedia.org/wiki/Buckingham_pi_theorem

There is in fact a Julia package centered around this theorem: https://github.com/rmsrosa/UnitfulBuckinghamPi.jl. I would though not recommend adopting the approach of using a global variable to store the parameters.

[yaccos]

Let us assume that the original system has units $U_1,U_2,\ldots,U_n$ where unit $U_i$ has dimensions $D_{1}^{a_{1, i}}D_{2}^{a_{2, i}}\ldots D_{\ell}^{a_{\ell,i}}$ with $D_{1},D_{2},\ldots,D_{\ell}$ being base dimensions. Each of the base dimensions are represented by its corresponding base reference unit $B_1,B_2,\ldots,B_{\ell}$. Let $c_i$ be the conversion factor such that $U_i = c_i B_{1}^{a_{1,i}}B_{2}^{a_{2,i}}\ldots B_{\ell}^{a_{\ell,i}}$. Let $A$ be the ($\ell\times n$) matrix with elements $a_{i,j}$. The matrix $A$ will thus represent the dimensions of the units in terms of the base units. $k=\mathrm{rank}A$ will then be the number of independent dimension of the unit system, this is: How many different units are necessary to represent the systems. For reference, $p=n-k$ is the number of dimensionsless parameters according to the Buckingham $\pi$ theorem.

Let us now consider how the original unit system can be transformed into another one with units $U_{1}^{*},U_{2}^{*},\ldots,U_{k}^{*}$, where everything is defined similarly to the first system expect for the constraint that the columns of the $\ell\times k$ matrix $A^*$ form a basis for the column space of $A$ and hence $A^*$ has full rank $k$. Hence, there exists a $k\times n$ matrix $Q$ with elements $q_{i,j}$ such that $A=A^*Q$. Hence, for each unit $U_i$ in the old system, we can express this in terms of the new system as

$$U_i=c_i\left(\frac{U_{1}^{*}}{c_{1}^*}\right)^{q_{1,i}}\left(\frac{U_{2}^{*}}{c_2^*}\right)^{q_{2,i}}\dots\left(\frac{U_{k}^{*}}{c_k^*}\right)^{q_{k,i}}$$.

In principle, there does not exist any real inverse $Q^{-1}$ for conversion back between the two systems. In practice for the ODE solvers this is not a problem because we are not constructing new kinds of quantities which are returned, but always feed the same system variables back in the same unit which were used to receive them.

[yaccos]

I think this package does something similar: https://github.com/martinkosch/Dimensionless.jl

[yaccos]

https://github.com/JuliaPhysics/NaturallyUnitful.jl also seems to be good inspiration for what we want to do.