Adding unit conversion for code generation

[skeating]

Code generation needs to be able to convert units. Currently Maurice is doing

factor = model.units.get_conversion_factor(1 * current_units, desired_units)
eq = sp.relational.Equality(eq.lhs, factor * eq.rhs)

which works so long as lhs isnt used elsewhere.

This issue will overlap with parts of other issues e.g. #8, #21 but adding it here to be clear what I'm working on.

[jonc125]

Principles:

• cellmlmanip should provide a general approach to adding conversions at model/world interfaces (whether the 'world' is Chaste or a WL protocol)
• It will do so by adding a new variable in the units desired by the world (if these are different from the model's units), transferring RDF annotations to the new variable (by just moving the cmeta:id in XML terms), and adding/changing equations as necessary so the new variable behaves appropriately.
• Inputs and outputs are with reference to the model. So an 'input' is when the world wants to change a model variable; 'output' is when the world wants to read a model variable.
• Only state variables, constants, or the free variable can be inputs. A variable calculated by an expression cannot (unless the expression is just a number, as this is really a constant).
• Any variable can be an output.

Some useful bits of pycml code to look at:

• Specific to Web Lab (all code in protocol.py falls in this category):
  • protocol.py:_split_all_odes adjusts all ODEs to account for time being units-converted. It assigns each RHS to a new variable in the original model units, and creates a new derivative in the new units to which this is assigned.
  • protocol.py:specify_as_output deals with creating new variables when needed for outputs
  • protocol.py:process_input_declarations deals with inputs
  • Actually adding the units conversions is done at a later processing stage in the finalize method by _add_units_conversions. With the simpler maths representation I'm hopeful we can do this more easily now!
• Specific to Chaste code generation (the InterfaceGenerator class is just used for Chaste):
  • processors.py:_split_ode does a similar transformation to _split_all_odes for handling a state variable being declared as an input (i.e. Chaste wishes to be able to adjust its value)
  • processors.py:_transform_derivatives_on_rhs will update all references to things like dV/dt on the RHS of any equation in the model to refer to the new variables created above, to ensure they're seeing the units they expect
  • Hopefully we can design cellmlmanip to avoid this duplication of logic...
• Code in the processors.py:ModelModifier class is common to both. Mainly this manages the fact that pycml keeps the notion of components, so adding/removing variables is complicated! It tries to update connections as required and ensure pycml's internal state is consistent.
• processors.py:UnitsConverter is a high-level interface for adding units conversions. Again most of the complexity is due to the maths representation. But there's also the logic for 'special' conversions that can convert between dimensions (notably used for currents).

[jonc125]

@MichaelClerx @MauriceHendrix @skeating let me know if you have any questions on the above :)

[MichaelClerx]

Only state variables, constants, or the free variable can be inputs. A variable calculated by an expression cannot (unless the expression is just a number, as this is really a constant).

So x = 1, y = x means y can't be an input, right? So only literal constants?

[jonc125]

Only state variables, constants, or the free variable can be inputs. A variable calculated by an expression cannot (unless the expression is just a number, as this is really a constant).

So x = 1, y = x means y can't be an input, right? So only literal constants?

That's correct, yes.

[MichaelClerx]

Let's spell out the "adding/changing equations as necessary" bit for inputs:

Literal constant

Example: Set current as an input, in nA

x = 1 [pA]
  in [pA]
  oxmeta: current
y = 1 / x
  in [1/pA]

would become

input_current = 1e-3 [nA]
  in [nA]
  oxmeta: current
x = input_current * 1e3 [pA/nA]
  in [pA]
y = 1 / x
  in [1/pA]

State variable

Example 1: Set voltage as an input in V

dV/dt = -1 / i_ion
  V in [mV]
  oxmeta: voltage
i_na = g * (V - E)
  in [pA]
x = 3 * dV/dt
...

would become

d(input_V)/dt =  original_dV_dt * 1e-3 [V/mV]
  input_V in [V]
  oxmeta: voltage
original_dV_dt =  -1 / i_ion
  in [mV/ms]
V = input_V * 1e3 [mV/V]
  in [mV]
i_na = g * (V - E)
  in [pA]
x = 3 * original_dV_dt

In this case we'd have to update the equation for some_variable, changing it from dV/dt (which is no longer available since V has been demoted to an ordinary variable) to original_dV_dt

Free variable

Set time as input in s

t = 0 [ms] (initial value)
  in [ms]
  oxmeta: time
stimulus = if(t < 5, 1, 0)
dV/dt = -1 / i_ion
  V in [mV]
  oxmeta: voltage
x = 3 * dV/dt

becomes

input_t = 0 [s] (initial value)
  in [s]
  oxmeta: time
t = input_t * 1000 [ms/s]
  in [ms]
stimulus = if(t < 5, 1, 0)
original_dV_dt = -1 / i_ion
  in [mV/ms]
dV/d_input_t = original_dV_dt * 1000 [ms/s]
  V in [mV]
  oxmeta: voltage
x = 3 * original_dV_dt

So in this case we also need to update some equations, changing any d/dt to a d/d_input_t

[jonc125]

Literal would become:

input_current = 1e-3 [nA]

but otherwise as you have it.

The state variable example would be as you've written it, although I'd call input_derivative_V original_dV_dt instead!

For the free variable, you alter every single ODE in a similar way to the state variable example. So something like

input_t = 0 [s] (initial value)
  in [s]
  oxmeta:time
t = input_t * 1000 [ms/s]
  in [ms]
orig_dV_dt_rhs =  -1 / i_ion
  in [mV/ms]
new_dV_dt_rhs = orig_dV_dt_rhs * 1000 [ms/s]
  in [mV/s]
dot(V) = new_dV_dt_rhs
  in [mV]
  oxmeta:voltage

and all references to dV/dt get changed to use orig_dV_dt_rhs instead.

[MichaelClerx]

Thanks! Have updated accordingly.
Why would the input_current RHS be evaluated/simplified? What if it's more complex, e.g. x = 2 * (1 [pA] + 2 [pA]), would that lead to input_current = 6e-3 [nA] ?

[jonc125]

But that x couldn't be an input, as it's a complex computed variable. Because the protocol (or Chaste) needs to be able to set the value of an input, which means it needs to be defined as a simple number. In the model code, it needs to be either in the parameters list, or the state vars list, or time, in order to be settable.

[MichaelClerx]

Oh ok so an input needs to be a number!
...
Turns out I've been thinking of expressions like 1+1 as "literals", since they don't reference any variable. But that's not the definition of a literal at all 😅

[MauriceHendrix]

All this is exactly how I understood it as well.

Moving the cmeta:id is important as otherwise I can't recognize the relevant variables.
I haven't looked at the code in great detail but the new variables like original_dV_dt might need some metadata and maybe a made up cmeta:id as I think various other parts of the code expect variables to have this.

What I'm wondering about is the units: will you put the ones asked for? And what if those are not in the model? It would be very nice if we didn't have to add those units to the model first (like is the case now).
Also what if I ask for say ms and the model already has millisecond. will we just have some bits in ms and some in millisecond? That's fine by me but we should make sure we warn users NOT to do `model. model.units.summarise_units(term1) == model. model.units.summarise_units(term2)

Finally the big question: who is going to implement all this :)

[jonc125]

The new variables shouldn't necessarily need RDF. Adding a cmeta:id should I think only be done at the point of need. So if code decides it wants to add an annotation for a variable, the id should be created at that point (so cellmlmanip will need a Model method to add an annotation, which can handle this, for instance).

I think units should be added transparently as needed too. To some degree I believe pint can do this for us, as it 'creates' unit definition objects as needed from expressions? The possibly complex bit is as you say if you want a 'nice' name for the units, which only matters either when referring to them by name (which I don't think we should do after the parsing - much easier to refer to actual objects) or when reporting to users. In my opinion, given a units definition, summarise_units should check if it matches one defined in the model (or protocol, depending on context) and use that name if so. Otherwise (for derived units) it could make one up?