search

cellml / libcellml

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

Importing resources with cellml 1.1 #1169

Closed FinbarArgus closed 1 year ago

FinbarArgus commented 1 year ago

I am trying to import a model with separate module, parameters, and units cellml files. I end up getting an invalid model. When doing resolve_imports I get a warning saying the units can't be found. Not sure what I am doing wrong?


import os
from libcellml import Analyser, AnalyserModel, AnalyserExternalVariable, Generator, GeneratorProfile, \
                      Parser, Importer

def resolve_imports(model, base_dir, strict_mode):
    importer = Importer(strict_mode)
    importer.resolveImports(model, base_dir)
    _dump_issues("resolve_imports", importer)
    if model.hasUnresolvedImports():
        print("unresolved imports?")
    else:
        print("no unresolved imports.")
    return importer

def parse_model(filename, strict_mode):
    cellml_file = open(filename)
    parser = Parser(strict_mode)
    model = parser.parseModel(cellml_file.read())
    _dump_issues("parse_model", parser)
    return model

def flatten_model(model, importer):
    flat_model = importer.flattenModel(model)
    return flat_model

def _dump_issues(source_method_name, logger):
    if logger.issueCount() > 0:
        print('The method "{}" found {} issues:'.format(source_method_name, logger.issueCount()))
        for i in range(0, logger.issueCount()):
            print('    - {}'.format(logger.issue(i).description()))

model_dir = os.path.dirname(__file__) 
model_name = 'cpp_coupling'
model_file_path = os.path.join(model_dir, model_name + '.cellml')

# parse the model in non-strict mode to allow non CellML 2.0 models
model = parse_model(model_file_path, False)
model_dir = os.path.dirname(model_file_path)

# code for importing the modules, parameters, and units, files.

# resolve imports, in non-strict mode
importer = resolve_imports(model, model_dir, False)
print("model.hasUnresolvedImports() = {}".format(model.hasUnresolvedImports()))
# need a flattened model for analysing
flat_model = flatten_model(model, importer)

# create model analyser
a = Analyser()

# add some external variables etc here, for this example case do nothing

# analyse the model
a.analyseModel(flat_model)

analysed_model = a.model()
print(analysed_model.type())
if not analysed_model.isValid():
    print("model is not valid, aborting...")
    exit()

I have attached the directory with the relevant cellml files (with .cellml changed to .txt so I could upload them).

cpp_coupling.txt cpp_coupling_modules.txt cpp_coupling_parameters.txt cpp_coupling_units.txt

I am using the up to date main branch of libcellml

hsorby commented 1 year ago

We have recently discovered an issue with importing of Units, a fix is currently being worked on here #1167.

FinbarArgus commented 1 year ago

Good to know a fix for the units is on the way.

I've tried changing all of the units to dimensionless and not importing a units file as a temporary hack... And I still get an invalid model. However now the analysed_model.type() gives 7, whereas before it gave 3. From [https://libcellml.org/documentation/v0.4.0/api/classlibcellml_1_1AnalyserModel] it looks like 3 means invalid? But 7 shouldn't be possible (assuming numbering starts at 0). Any ideas what I'm doing wrong?

Here are the updated dimensionless files cpp_coupling.txt cpp_coupling_modules.txt cpp_coupling_parameters.txt

agarny commented 1 year ago

The documentation is not up to date. Quite a few changes have been added to the analyser and although all of our tests pass, who knows, maybe something got broken somewhere. In the meantime, regarding the type code of 7 that you are getting, it means that your model is under constrained (see here). Please, have another look at your model. If you reckon that it is fine then I will have a look at why the analyser categorised it as under constrainted.

FinbarArgus commented 1 year ago

Thanks Alan! It is running in OpenCOR so it should be correctly constrained. It seems that either the modules file or the parameters file isn't being read correctly. I might be doing the resolve_imports above incorrectly? Not too sure.

agarny commented 1 year ago

Well, OpenCOR uses the CellML API which is, unfortunately, known to validate malformed models. So, who knows, it might be one of those cases... or libCellML's analyser could indeed be buggy. Ok, I am going to look into it sometime today, using the models you posted here.

agarny commented 1 year ago

Ok, the fact that you mentioned OpenCOR should have made me realised that your model might not have been encoded in CellML 2.0 and... indeed, it is encoded in CellML 1.1! Anyway, at this stage, I can confirm the under constrained type.

agarny commented 1 year ago

Ok, I feel like there is an issue with resolving imports. @hsorby to confirm. Here is the flattened versio of the model:

<?xml version="1.0" encoding="UTF-8"?>
<model xmlns="http://www.cellml.org/cellml/2.0#" name="CardiovascularSystem">
  <component name="parameters">
    <variable name="R_main_vessel" units="dimensionless" initial_value="1e6" interface="public"/>
    <variable name="C_main_vessel" units="dimensionless" initial_value="1e-8" interface="public"/>
    <variable name="u_ext_main_vessel" units="dimensionless" initial_value="0.0" interface="public"/>
    <variable name="I_main_vessel" units="dimensionless" initial_value="1e-6" interface="public"/>
    <variable name="u_out_main_vessel" units="dimensionless" initial_value="10" interface="public"/>
    <variable name="v_in_main_vessel" units="dimensionless" initial_value="1e-5" interface="public"/>
  </component>
  <component name="environment">
    <variable name="time" units="dimensionless" interface="public"/>
  </component>
  <component name="main_vessel_module">
    <variable name="t" units="dimensionless" interface="public"/>
    <variable name="I" units="dimensionless" interface="public"/>
    <variable name="C" units="dimensionless" interface="public"/>
    <variable name="R" units="dimensionless" interface="public"/>
    <variable name="R_v" units="dimensionless" interface="public"/>
    <variable name="u_ext" units="dimensionless" interface="public"/>
    <variable name="v_in" units="dimensionless" interface="public"/>
    <variable name="u" units="dimensionless" interface="public"/>
    <variable name="q_C" units="dimensionless" initial_value="0.0" interface="public"/>
    <variable name="u_C" units="dimensionless" interface="public"/>
    <variable name="v" units="dimensionless" initial_value="0.0" interface="public"/>
    <variable name="u_out" units="dimensionless" interface="public"/>
    <math xmlns="http://www.w3.org/1998/Math/MathML" xmlns:cellml="http://www.cellml.org/cellml/2.0#">
      <apply>
        <eq/>
        <ci>R_v</ci>
        <apply>
          <divide/>
          <cn cellml:units="dimensionless">0.01</cn>
          <ci>C</ci>
        </apply>
      </apply>
      <!-- Conservation Laws-->
      <apply>
        <eq/>
        <apply>
          <diff/>
          <bvar>
            <ci>t</ci>
          </bvar>
          <ci>v</ci>
        </apply>
        <apply>
          <divide/>
          <apply>
            <minus/>
            <apply>
              <minus/>
              <ci>u</ci>
              <ci>u_out</ci>
            </apply>
            <apply>
              <times/>
              <ci>R</ci>
              <ci>v</ci>
            </apply>
          </apply>
          <ci>I</ci>
        </apply>
      </apply>
      <apply>
        <eq/>
        <apply>
          <diff/>
          <bvar>
            <ci>t</ci>
          </bvar>
          <ci>q_C</ci>
        </apply>
        <apply>
          <minus/>
          <ci>v_in</ci>
          <ci>v</ci>
        </apply>
      </apply>
      <apply>
        <eq/>
        <ci>u_C</ci>
        <apply>
          <plus/>
          <apply>
            <divide/>
            <ci>q_C</ci>
            <ci>C</ci>
          </apply>
          <ci>u_ext</ci>
        </apply>
      </apply>
      <apply>
        <eq/>
        <ci>u</ci>
        <apply>
          <plus/>
          <ci>u_C</ci>
          <apply>
            <times/>
            <ci>R_v</ci>
            <apply>
              <minus/>
              <ci>v_in</ci>
              <ci>v</ci>
            </apply>
          </apply>
        </apply>
      </apply>
    </math>
  </component>
  <component name="main_vessel">
    <variable name="u" units="dimensionless" interface="public"/>
    <variable name="v" units="dimensionless" interface="public"/>
    <variable name="q_C" units="dimensionless" interface="public"/>
    <variable name="R" units="dimensionless" interface="public"/>
    <variable name="C" units="dimensionless" interface="public"/>
    <variable name="u_ext" units="dimensionless" interface="public"/>
    <variable name="I" units="dimensionless" interface="public"/>
    <variable name="u_out" units="dimensionless" interface="public"/>
    <variable name="v_in" units="dimensionless" interface="public"/>
  </component>
  <connection component_1="environment" component_2="main_vessel_module">
    <map_variables variable_1="time" variable_2="t"/>
  </connection>
  <connection component_1="main_vessel_module" component_2="main_vessel">
    <map_variables variable_1="I" variable_2="I"/>
    <map_variables variable_1="C" variable_2="C"/>
    <map_variables variable_1="R" variable_2="R"/>
    <map_variables variable_1="u_ext" variable_2="u_ext"/>
    <map_variables variable_1="v_in" variable_2="v_in"/>
    <map_variables variable_1="u" variable_2="u"/>
    <map_variables variable_1="q_C" variable_2="q_C"/>
    <map_variables variable_1="v" variable_2="v"/>
    <map_variables variable_1="u_out" variable_2="u_out"/>
  </connection>
</model>

and the analyser reports the following issues:

If we look at R_v, we can see that it's computed using R_v = 0.01/C and that it is used somewhere else. So, as far as R_V is concerned, it's all good. However, it depends on C and looking at the flattened version of the model, we can see that it used in a couple of equations, but it is never computed and it doesn't have an initial value. However, looking at the original model, we can see that C is connected to C_main_vessel which has an initial value of 1e-8, but this somehow didn't make it through the import process and, therefore, through the flatttening. The same, for instance, with R, which should be initialised to 1e6 in the flattened model, but it is not.

So, yes, it looks like the resolving of imports somehow got broken. @hsorby, do you agree with my "analysis"?

hsorby commented 1 year ago

I have found an error in the code and made a correction that will work for the dimensionless form of the models given here. It won't work with the Unitized version because those models are not valid enough for the importer. The Unitized models must import all required units from the units definitions model. The modules model doesn't do this and thus cannnot be flattened due to insufficient knowledge of the units used by the imported component.

hsorby commented 1 year ago

Here is the dimensionless model I generated with the bug fixed code:

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

from enum import Enum
from math import *

__version__ = "0.4.0"
LIBCELLML_VERSION = "0.4.0"

STATE_COUNT = 2
VARIABLE_COUNT = 9

class VariableType(Enum):
    VARIABLE_OF_INTEGRATION = 0
    STATE = 1
    CONSTANT = 2
    COMPUTED_CONSTANT = 3
    ALGEBRAIC = 4

VOI_INFO = {"name": "time", "units": "dimensionless", "component": "environment", "type": VariableType.VARIABLE_OF_INTEGRATION}

STATE_INFO = [
    {"name": "v", "units": "dimensionless", "component": "main_vessel_module", "type": VariableType.STATE},
    {"name": "q_C", "units": "dimensionless", "component": "main_vessel_module", "type": VariableType.STATE}
]

VARIABLE_INFO = [
    {"name": "R_main_vessel", "units": "dimensionless", "component": "parameters", "type": VariableType.CONSTANT},
    {"name": "C_main_vessel", "units": "dimensionless", "component": "parameters", "type": VariableType.CONSTANT},
    {"name": "u_ext_main_vessel", "units": "dimensionless", "component": "parameters", "type": VariableType.CONSTANT},
    {"name": "I_main_vessel", "units": "dimensionless", "component": "parameters", "type": VariableType.CONSTANT},
    {"name": "u_out_main_vessel", "units": "dimensionless", "component": "parameters", "type": VariableType.CONSTANT},
    {"name": "v_in_main_vessel", "units": "dimensionless", "component": "parameters", "type": VariableType.CONSTANT},
    {"name": "R_v", "units": "dimensionless", "component": "main_vessel_module", "type": VariableType.COMPUTED_CONSTANT},
    {"name": "u", "units": "dimensionless", "component": "main_vessel_module", "type": VariableType.ALGEBRAIC},
    {"name": "u_C", "units": "dimensionless", "component": "main_vessel_module", "type": VariableType.ALGEBRAIC}
]

def create_states_array():
    return [nan]*STATE_COUNT

def create_variables_array():
    return [nan]*VARIABLE_COUNT

def initialise_variables(states, rates, variables):
    variables[0] = 1.0e6
    variables[1] = 1.0e-8
    variables[2] = 0.0
    variables[3] = 1.0e-6
    variables[4] = 10.0
    variables[5] = 1.0e-5
    states[0] = 0.0
    states[1] = 0.0

def compute_computed_constants(variables):
    variables[6] = 0.01/variables[1]

def compute_rates(voi, states, rates, variables):
    variables[8] = states[1]/variables[1]+variables[2]
    variables[7] = variables[8]+variables[6]*(variables[5]-states[0])
    rates[0] = (variables[7]-variables[4]-variables[0]*states[0])/variables[3]

def compute_variables(voi, states, rates, variables):
    variables[8] = states[1]/variables[1]+variables[2]
    variables[7] = variables[8]+variables[6]*(variables[5]-states[0])
                                                                                                                                                                                                            33,1
FinbarArgus commented 1 year ago

Thanks Hugh! The dimensionless model works for me too after pulling your branch.

Is there an easy way to make all units dimensionless in a model? So that I can do that temporarily for now without changing to dimensionless by hand.

Cheers

agarny commented 1 year ago

Is there an easy way to make all units dimensionless in a model? So that I can do that temporarily for now without changing to dimensionless by hand.

This should do what you want (from a BASH prompt):

while IFS= read -r line
do
    if [[ "$line" == *"<variable"* ]]; then
        echo $line | sed 's/units=\".*\"/units=\"dimensionless\"/g'
    else
        echo $line
    fi
done < model.cellml > new.model.cellml

It goes through every single line of a file called model.cellml and check whether it contains <variable. If it does, it will replace the units with dimensionless. If the line doesn't contain <variable (e.g., a line in a units definition) then it will just output the line.

hsorby commented 1 year ago

You still need to worry about cn elements as well. A similar tactic would work for that as well. One of the problems that faces libCellML is the fact that it treats math as a string as much as it can, so it has limited ability to modify any of the math.

agarny commented 1 year ago

Oops, indeed, completely overlooked cn elements. Here goes for a revised version:

while IFS= read -r line
do
    if [[ "$line" != *"<unit"* ]]; then
        echo $line | sed 's/units=\".*\"/units=\"dimensionless\"/g'
    else
        echo $line
    fi
done < model.cellml > new.model.cellml

Note that it only works if an element, that contains a units attribute that you want to change, is not over several lines, i.e.

<variable ... units="xxx" /> <!-- OK -->
<variable ...
units="xxx" /> <!-- NOK -->
hsorby commented 1 year ago

I think I would be more inclined to redefine all the units to be dimensionless, using something like the following:

import os.path

import libcellml
import sys

def main():
    args = sys.argv[:]
    model_file = args.pop()
    with open(model_file) as f:
        content = f.read()

    parser = libcellml.Parser(False)

    model = parser.parseModel(content)

    for units_index in range(model.unitsCount()):
        units = model.units(units_index)
        units.removeAllUnits()
        units.addUnit("dimensionless")

    printer = libcellml.Printer()
    new_content = printer.printModel(model)
    split = os.path.splitext(model_file)
    dimensionless_model_file = f"{split[0]}_dimensionless{split[1]}"
    with open(dimensionless_model_file, "w") as f:
        f.write(new_content)

if __name__ == "__main__":
    main()
hsorby commented 1 year ago

With the fixes from this PR I managed to get a Units appropriate generated python file that was the same as the dimensionless one. I did have one error to fix though:

The method "analyse_model" found 1 issues:
    - Variable 'q_C' in component 'main_vessel_module' has units of 'J_per_m3' and an equivalent variable 'q_C' in component 'main_vessel' with non-matching units of 'm3'. The mismatch is: kilogram^1, metre^-4, second^-2.

Two of the q_C variables are defined with J_per_m3, after changing their units to m3 I was successfully able to generate python code that I could compute.

FinbarArgus commented 1 year ago

Great catch, thanks Hugh!

Could you please send through the python file you used to do the generation?

I have made the same fix to the q_C unit, and included unit imports in the parameters and modules files. (Is this the expected usage going forwards? Or will we be able to import units only once in the main file in the future?)

I am still getting an invalid analysed model type when I try to analyse the model.

Cheers, Finbar

hsorby commented 1 year ago

Sorry, can't do the continual attaching files method. Here are the files available in this repo: https://github.com/hsorby/cpp_coupling_example

When we are importing units and components everything they rely on from the imported model needs to be defined, the units or component are only defined from what is in that model.
This means we need to also import common units/components from a units definition file/components definition file if a units or component makes use of them.
The imported model itself can be invalid or ill-defined in other parts that are not imported by the importing model without effecting your ability to import a units or component from the imported model.

FinbarArgus commented 1 year ago

I had a couple of units defined twice, which broke things. All is working now! Thanks for the fixes!