Analyser unit check raises the multiplier of a unit to its exponent, in contrast to the specification and Units::scalingFactor

[matthiaskoenig]

The CellML 2.0 specification (section 3.3, rule 1.4) interprets a unit element as

multiplier * (10^prefix * units)^exponent

i.e., the exponent applies to the prefix, not to the multiplier. Units::scalingFactor follows this for a reference to standard units, but the unit check of the Analyser also raises the multiplier to the exponent, so it reports consistent equations as not equivalent. A second, related difference: for a reference to non-standard units Units::scalingFactor does not raise the prefix to the exponent.

Found with libcellml 0.7.1 (python, Linux) while writing CellML units in sbml2cellml; the code on main (33f65ba) is the same.

Example

import libcellml

CELLML = """<?xml version="1.0" encoding="UTF-8"?>
<model xmlns="http://www.cellml.org/cellml/2.0#" name="per_minute">
  <units name="minute">
    <unit units="second" multiplier="60"/>
  </units>
  <!-- 1/60 second^-1, the multiplier is outside of the exponent -->
  <units name="per_minute">
    <unit units="second" exponent="-1" multiplier="0.0166666666666667"/>
  </units>
  <!-- the same units, by reference -->
  <units name="per_minute_by_reference">
    <unit units="minute" exponent="-1"/>
  </units>
  <component name="main">
    <variable name="t" units="minute"/>
    <variable name="x" units="mole" initial_value="1"/>
    <variable name="k" units="per_minute" initial_value="0.1"/>
    <math xmlns="http://www.w3.org/1998/Math/MathML">
      <apply><eq/>
        <apply><diff/><bvar><ci>t</ci></bvar><ci>x</ci></apply>
        <apply><times/><ci>k</ci><ci>x</ci></apply>
      </apply>
    </math>
  </component>
</model>
"""

print("libcellml", libcellml.versionString())
model = libcellml.Parser().parseModel(CELLML)

print(
    "scalingFactor(per_minute, per_minute_by_reference):",
    libcellml.Units.scalingFactor(
        model.units("per_minute"), model.units("per_minute_by_reference")
    ),
)

analyser = libcellml.Analyser()
analyser.analyseModel(model)
for i in range(analyser.issueCount()):
    print(analyser.issue(i).description())

# second case: the prefix of a reference to non-standard units
units_model = libcellml.Model("prefix")
b = libcellml.Units("b")
b.addUnit("metre")
a = libcellml.Units("a")  # (kilo b)^2
a.addUnit("b", "kilo", 2.0, 1.0)
c = libcellml.Units("c")  # (kilo metre)^2
c.addUnit("metre", "kilo", 2.0, 1.0)
for units in (b, a, c):
    units_model.addUnits(units)
units_model.linkUnits()
print("scalingFactor((kilo b)^2, (kilo metre)^2):", libcellml.Units.scalingFactor(a, c))

Output:

libcellml 0.7.1
scalingFactor(per_minute, per_minute_by_reference): 0.999999999999998
The units in 'dx/dt = k*x' in component 'main' are not equivalent. 'dx/dt' is in 'minute^-1 x mole' (i.e. '10^-1.77815 x mole x second^-1') while 'k*x' is in 'mole x per_minute' (i.e. '10^1.77815 x mole x second^-1').
scalingFactor((kilo b)^2, (kilo metre)^2): 1000.0

1. Analyser: the multiplier is raised to the exponent

per_minute is 1/60 second^-1, which scalingFactor confirms (it equals minute^-1), so dx/dt = k*x is consistent. The analyser takes per_minute as 10^1.77815 second^-1, i.e. (1/60 second)^-1 = 60 second^-1, and warns. With k in per_minute_by_reference there is no warning.

The cause is in Analyser::AnalyserImpl::updateUnitsMultiplier, where log10(multiplier) is inside the factor exponent:

libcellml/src/analyser.cpp

Line 1401 in 33f65ba

newUnitsMultiplier += unitsMultiplier + (standardMultiplierList.at(reference) + std::log10(multiplier) + convertPrefixToInt(prefix)) * exponent * unitsExponent;

newUnitsMultiplier += unitsMultiplier + (standardMultiplierList.at(reference) + std::log10(multiplier) + convertPrefixToInt(prefix)) * exponent * unitsExponent;

while updateUnitMultiplier in units.cpp keeps it outside, as the specification does:

libcellml/src/units.cpp

Line 256 in 33f65ba

localMultiplier += mult + (standardMult + prefixMult) * exp;

localMultiplier += mult + (standardMult + prefixMult) * exp;

I would expect (std::log10(multiplier) + (standardMultiplierList.at(reference) + convertPrefixToInt(prefix)) * exponent) * unitsExponent in the analyser.

2. Units::scalingFactor: the prefix of a reference to non-standard units is not raised to the exponent

(kilo b)^2 with b = metre should equal (kilo metre)^2, the scaling factor is 1000 instead of 1. In the branch for non-standard units the prefix is outside of the exponent:

libcellml/src/units.cpp

Line 269 in 33f65ba

localMultiplier += mult + branchMult * exp + prefixMult;

localMultiplier += mult + branchMult * exp + prefixMult;

I would expect mult + (branchMult + prefixMult) * exp, like for standard units. (The analyser raises the prefix to the exponent in this case.)

Both are about the same formula in two places, so I report them together; happy to split the issue if you prefer.

[agarny]

Good catch @matthiaskoenig! I have just fixed it and created a PR for it here.

Note that I used the result of 1.0/60.0 rather than 0.0166666666666667 to avoid tolerance issues (see areNearlyEqual() and Units::equivalent(), currently used by the Analyser and Units classes respectively).