ModellingWebLab / cellmlmanip

CellML loading and model equation manipulation
Other
3 stars 1 forks source link

Allow more complex unit transformation rules #282

Open MichaelClerx opened 4 years ago

MichaelClerx commented 4 years ago

This issue has two goals:

#!/usr/bin/env python
import pint

r = pint.UnitRegistry()

r.define('uA=ampere * 1e-6')
r.define('pA=ampere * 1e-12')
r.define('uF=farad * 1e-6')
r.define('pF=farad * 1e-12')
r.define('cm2=centimetre ** 2')
r.define('uA_per_cm2=uA / cm2')
r.define('uF_per_cm2=uF / cm2')
r.define('A_per_F=ampere / farad')

Cm = r.Quantity(12, r.pF)
Cs = r.Quantity(1.1, r.uF_per_cm2)

# Rules
c = pint.Context()
c.add_transformation(r.A_per_F, r.uA_per_cm2, lambda ureg, rhs: rhs * Cs)
c.add_transformation(r.uA_per_cm2, r.A_per_F, lambda ureg, rhs: rhs / Cs)
c.add_transformation(r.uA, r.uA_per_cm2, lambda ureg, rhs: rhs * Cs / Cm)
c.add_transformation(r.uA_per_cm2, r.uA, lambda ureg, rhs: rhs * Cm / Cs)
r.enable_contexts(c)

# Test
x = r.Quantity(1, r.pA)
print(x.to(r.A_per_F))

The above is an example of pint transformations.

There's some subtle gotchas:

# Gotchas
r.disable_contexts()
c = pint.Context()
c.add_transformation(r.A_per_F, r.uA_per_cm2, lambda ureg, rhs: rhs * Cs)
c.add_transformation(r.uA_per_cm2, r.A_per_F, lambda ureg, rhs: rhs / Cs)
f = Cs / Cm
c.add_transformation(r.uA, r.uA_per_cm2, lambda ureg, rhs: rhs * f)
f = Cm / Cs
c.add_transformation(r.uA_per_cm2, r.uA, lambda ureg, rhs: rhs * f)
r.enable_contexts(c)

# Test
x = r.Quantity(1, r.pA)
print(x.to(r.A_per_F))

this will break, because the lambda doesn't store a copy of f, it just grabs an f from the local scope whenever it's used, so both functions actually use the same f. So some tricks needed there (e.g. define a callable class).

But the main thing to figure out is whether we can come up with some cellmlmanip API that can be re-used by FC (Web Lab) and chaste codegen, for specifying rules with lambdas (python lambdas? sympy expressions?) so that conversion of currents is easier

MauriceHendrix commented 4 years ago

What I currently use for codegen (outsode pints / cellmlmanip) would need to be able to generate a function call that ends up being in the output. E.g. lambda eq: sympy.Eq(eq.lhs, eq.rhs * sp.Function('HeartConfig::Instance()->GetCapacitance', real=True)() so I think that may not be possible with this.

MichaelClerx commented 4 years ago

Could still work! @jonc125 wrapped sympy objects in pint.Quantity, and I don't see why calling sp.Function(..) couldn't be part of the lambda you passed in? Also, you could just store the output of the function in a local var?

MauriceHendrix commented 4 years ago

The function has no output, it needs to be there as a cal. In other woulds the C++ output is of the for: const double xyz = .... * HeartConfig::Instance()->GetCapacitance()

If you can wrap that in a quantity maybe it's possible, where could I find this wrapping in action to look at?

jonc125 commented 4 years ago

There's an example of the wrapping in the description.

I also wondered about using Sympy's lambdify to convert from a Sympy expression to an equivalent Python that could operate on the Pint units. If you had Pint Quantity instances with Sympy expressions as magnitudes (including UnevaluatedFunction for your capacitance example, and our Variable/Quantity classes for variables & numbers) then applying such a Python-ised function to them should I think result in a Quantity with the desired units, and the Sympy expression you need as the magnitude.

MauriceHendrix commented 4 years ago

There's an example of the wrapping in the description.

I might be blind but I can't see it?

jonc125 commented 4 years ago

Ah, good point, sorry, didn't realise Michael had done something different. My example was at https://github.com/ModellingWebLab/weblab-fc/issues/156#issuecomment-614675629

jonc125 commented 4 years ago

The second point is handled by #286.