cellml / libcellml

Repository for libCellML development.
https://libcellml.org
Apache License 2.0
17 stars 21 forks source link

Analyser: undesired behaviour with external variables #958

Closed hsorby closed 3 years ago

hsorby commented 3 years ago

If I have the following model:

<?xml version='1.0' encoding='UTF-8'?>
<model name="test_model" xmlns="http://www.cellml.org/cellml/2.0#" xmlns:cellml="http://www.cellml.org/cellml/2.0#">
    <units name="centimeter">
        <unit prefix="centi" units="metre"/>
    </units>
    <units name="centimeter_squared">
        <unit exponent="2" units="centimeter"/>
    </units>
     <units name="microlitre">
        <unit prefix="micro" units="litre"/>
    </units>
    <units name="microlitre_per_centimeter_cubed">
        <unit units="microlitre"/>
        <unit exponent="-3" units="centimeter"/>
    </units>
   <component name="cell_geometry">
        <variable initial_value="0.01" name="L" units="centimeter"/>
        <variable initial_value="0.0011" name="rad" units="centimeter"/>
        <variable name="vcell" units="microlitre" interface="public"/>
        <variable name="vss" units="microlitre" interface="public"/>
        <math xmlns="http://www.w3.org/1998/Math/MathML">
            <apply>
                <eq/>
                <ci>vcell</ci>
                <apply>
                    <times/>
                    <cn cellml:units="microlitre_per_centimeter_cubed">1000</cn>
                    <cn cellml:units="dimensionless">3.14</cn>
                    <ci>rad</ci>
                    <ci>rad</ci>
                    <ci>L</ci>
                </apply>
            </apply>
            <apply>
                <eq/>
                <ci>vss</ci>
                <apply>
                    <times/>
                    <cn cellml:units="dimensionless">0.02</cn>
                    <ci>vcell</ci>
                </apply>
            </apply>
        </math>
    </component>
</model>

and I run it through the following Python script


import libcellml as lc

def main():
    with open('test_model.cellml') as f:
        content = f.read()

    p = lc.Parser()
    m = p.parseModel(content)

    v = lc.Validator()
    v.validateModel(m)

    print(v.errorCount())

    a = lc.Analyser()

    a.addExternalVariable(lc.AnalyserExternalVariable(m.component("cell_geometry").variable("L")))
    a.addExternalVariable(lc.AnalyserExternalVariable(m.component("cell_geometry").variable("rad")))

    a.analyseModel(m)

    am = a.model()

    g = lc.Generator()
    g.setModel(am)

    pr = lc.GeneratorProfile(lc.GeneratorProfile.Profile.PYTHON)
    g.setProfile(pr)

    print(g.implementationCode())

if __name__ == "__main__":
    main()

I get the following Python implementation code

# The content of this file was generated using the Python profile of libCellML 0.2.0.

from enum import Enum
from math import *

__version__ = "0.3.0"
LIBCELLML_VERSION = "0.2.0"

VARIABLE_COUNT = 4

class VariableType(Enum):
    CONSTANT = 1
    COMPUTED_CONSTANT = 2
    ALGEBRAIC = 3
    EXTERNAL = 4

VARIABLE_INFO = [
    {"name": "L", "units": "centimeter", "component": "cell_geometry", "type": VariableType.EXTERNAL},
    {"name": "rad", "units": "centimeter", "component": "cell_geometry", "type": VariableType.EXTERNAL},
    {"name": "vcell", "units": "microlitre", "component": "cell_geometry", "type": VariableType.ALGEBRAIC},
    {"name": "vss", "units": "microlitre", "component": "cell_geometry", "type": VariableType.COMPUTED_CONSTANT}
]

def create_variables_array():
    return [nan]*VARIABLE_COUNT

def initialise_constants(variables):
    pass

def compute_computed_constants(variables):
    variables[3] = 0.02*variables[2]

def compute_variables(variables, external_variable):
    variables[0] = external_variable(variables, 0)
    variables[1] = external_variable(variables, 1)
    variables[2] = 1000.0*3.14*variables[1]*variables[1]*variables[0]

Which I think has the problem in the compute_computed_constants function where the equation for setting variable[2] does not appear. Also, I cannot initialise the external variables.

agarny commented 3 years ago

Agreed that something funny is happening here. I believe compute_computed_constants and compute_variables should actually be defined as follows since we cannot guarantee that variables[3] is a computed constant (since variables[0] and variables[1] are external variables):

def compute_computed_constants(variables):
    pass

def compute_variables(variables, external_variable):
    variables[0] = external_variable(variables, 0)
    variables[1] = external_variable(variables, 1)
    variables[2] = 1000.0*3.14*variables[1]*variables[1]*variables[0]
    variables[3] = 0.02*variables[2]

Otherwise, good point about initialising the external variables. I guess we should have a new method for the case where we have external variables. Something like:

def initialise_external_variables(variables, external_variable):
    variables[0] = external_variable(variables, 0)
    variables[1] = external_variable(variables, 1)

Something to discuss at our next libCellML meeting @hsorby and @nickerso.

agarny commented 3 years ago

After our weekly meeting, we agreed that we should have the following (@hsorby and @nickerso to comment if I got it wrong):

Algebraic model + external variables:

def initialise_variables(variables, external_variable): ... variables[index] = external_variable(variables, index)

Differential model:

def initialise_variables(states, variables): ...

Differential model + external variables:

def initialise_variables(voi, states, variables, external_variable): variables[index] = external_variable(voi, states, variables, index)

- `compute_rates()` is to be updated as follows:
```python
# Differential model:
def compute_rates(voi, states, rates, variables):
    ...

# Differential model + external variables:
def compute_rates(voi, states, rates, variables, external_variable):
    variables[index] = external_variable(voi, states, variables, index)

Algebraic model + external variables:

def compute_variables(variables, external_variable): ... variables[index] = external_variable(variables, index)

Differential model:

def compute_variables(voi, states, rates, variables): ...

Differential model + external variables:

def compute_variables(voi, states, rates, variables, external_variable): variables[index] = external_variable(voi, states, variables, index)

hsorby commented 3 years ago

Can you confirm that you need voi for an algebraic equation with external variables. I don't see why it should be there.

agarny commented 3 years ago

That's, I believe, what @nickerso was talking about this morning, i.e. the user may have some recorded data that s/he wants to use as the value of some external variable. To retrieve the correct value from that recorded data, to have access to voi will (likely?) be needed.

nickerso commented 3 years ago

Algebraic models do not have a voi, right? so "time" would just be a variable...

agarny commented 3 years ago

Oops, sorry, yes of course! Was thinking about initialise_variables() for a differential model + external variables. My bad. I have corrected my comment.