Skip to content
Snippets Groups Projects
Commit 411cacdf authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'feature/cell-model-translate' into 'master'

Cell model translator for cardiac solver.

This adds a python utility for generating cell model code from CellML models.

See merge request !723
parents 2db0ae69 c5d533fc
No related branches found
No related tags found
No related merge requests found
Showing
with 22087 additions and 3 deletions
......@@ -58,6 +58,9 @@ v4.4.0
- Added a convective like outflow boundary condition from Dong (!713)
- Added the ability to specifiy Womersley boundary conditions for pulsatile flow (!472)
**CardiacEPSolver:**
- Added a Python translator utility to generate cell models from CellML (!723)
**FieldConvert:**
- Allow equi-spaced output for 1D and 2DH1D fields (!613)
- Update quality metric to include scaled Jacobian output (!695)
......
......@@ -484,6 +484,7 @@ set(CPACK_SOURCE_IGNORE_FILES
"/library/Demos/LocalRegions/XmlFiles/"
"/solvers/ImageWarpingSolver/"
"/solvers/VortexWaveInteraction/"
"/solvers/CardiacEPSolver/Utilities/CellMLToNektar/"
"/solvers/IncNavierStokesSolver/Utilities/HybridPerformanceModel/"
"/utilities/Extras/"
)
......
......@@ -12,7 +12,7 @@
#include <iostream>
#include <string>
//#include <LibUtilities/BasicUtils/Vmath.hpp>
//#inc lude <LibUtilities/BasicUtils/Vmath.hpp>
#include <CardiacEPSolver/CellModels/LuoRudy91.h>
namespace Nektar
......
......@@ -60,7 +60,7 @@ std::string TenTusscher06::lookupIds[4] = {
std::string TenTusscher06::def =
LibUtilities::SessionReader::RegisterDefaultSolverInfo(
"CellModelVariant", "eEpicardium");
"CellModelVariant", "Epicardium");
/**
*
......
#!/usr/bin/env python
"""Copyright (c) 2005-2016, University of Oxford.
All rights reserved.
University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.
This file is part of Chaste.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
* Neither the name of the University of Oxford nor the names of its
contributors may be used to endorse or promote products derived from this
software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
"""
"""
This part of PyCml deals with converting CellML models into programming language code.
It is a thin executable wrapper around translators.py.
"""
import os #importing functions related to operating system
import sys #importing functions related to system (e.g. path directory etc)
# Make sure PyCml is on sys.path
pycml_path = os.path.dirname(os.path.realpath(__file__))
sys.path[0:0] = [pycml_path]
import translators #importing the main translator class
from translators import CellMLTranslator #again, just to make sure
import CellMLToNektarTranslator #importing the nektar sub-class
#This part actually runs the code
if __name__ == '__main__': #this part of the code is only run when CellMLToNektar.py is run directly, not imported
CellMLTranslator.register(CellMLTranslator, 'C++') #this registers the main translator and it's name when called as an option in modified_translators
CellMLTranslator.register(CellMLToNektarTranslator.CellMLToNektarTranslator, 'Nektar') #this registers the nektar subclass and it's name when called as an option in modified_translators
if '--profile' in sys.argv: #sys.argv is the list of command line items passed to the script, this checks if a specific translator is requested
import time, cProfile
profile_name = '/tmp/pycml-profile-%f-%d' % (time.time(), os.getpid())
cProfile.run('modified_translators.run()', profile_name) #the requested translator is run
else:
translators.run() #if no specific translator is requested, the default run function in modified_translators is run
#Below here seem to be variables that do some form of testing
# For use in testing
def euler(doc, t, nsteps=1000, dt=0.01):
global tvar, state_vars, exprs #assigns tvar, state_vars and exprs to be global variables, meaning they can be used by any function
tvar = t.free_vars[0] #takes the t-input and does "free_vars[0]" to it?
state_vars = t.state_vars #defines state_vars as the state_vars of the t-input?
for var in state_vars: #cycles through all the entries in "state_vars"
var.set_value(float(var.initial_value)) #for all the entries in "state_vars" it sets the value to a float of the initial value
tvar.set_value(0.0) #this sets the value of tvars to (0.0)?
exprs = [e for e in doc.model.get_assignments() #I don't understand this part
if isinstance(e, modified_translators.mathml_apply)] #what is the mathml_apply function? Can't find it in modified_translators
for _ in range(nsteps): #this is a for-loop through all the values up to the input of nsteps
for expr in exprs: #this cycles through all the entries in exprs
expr.evaluate()
tvar.set_value(tvar.get_value() + dt)
for var in state_vars:
var.set_value(var.get_value() +
dt * var.get_value(ode=tvar))
return
def writefile(doc, outfn='test.cml'):
# Write out CellML file
st = modified_translators.open_output_stream(outfn)
doc.xml(indent=1, stream=st)
st.close()
return
def show_usage(doc):
for comp in doc.model.component:
for var in comp.variable:
print var.fullname(), var._cml_usage_count
def fix_divide_by_zero(doc):
"""
Several models have equations of a form that may give rise to
a divide by zero error on simulation, especially when lookup
tables are used. The general form is:
(a * (V - v0)) / (exp(b * (V - v0)) - 1)
When V = v0 this is undefined, however the limit of the
function as V approaches v0 from either side is well-defined,
and each limit is the same. We approximate the limit by
linear interpolation between values of the expression for
(V-v0) = +/- 1e-10.
"""
divides = [d.xml_parent
for d in doc.xml_xpath(u'//m:apply/m:divide')]
for divide in divides:
pass
return
This diff is collapsed.
This diff is collapsed.
"""Copyright (c) 2005-2016, University of Oxford.
All rights reserved.
University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.
This file is part of Chaste.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice,
this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
* Neither the name of the University of Oxford nor the names of its
contributors may be used to endorse or promote products derived from this
software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
"""
"""
This module abstracts the interface to RDF metadata about CellML models.
The RdfProcessor class below pretends to be the module itself, so all its properties
are available at module-level, and these should typically be called by users.
It also provides sets METADATA_NAMES and STIMULUS_NAMES, which contain the local names
of terms in the ontology that can annotate variables, and the subset of those names
which define properties of the stimulus current (but not the current itself), respectively.
"""
import logging
import os
import sys
import types
from cStringIO import StringIO
# We now only support rdflib for RDF processing
import rdflib
def __init__(module):
# Import pycml here, to avoid circular import surprises
import pycml
module.pycml = pycml
class RdfProcessor(object):
"""Implements CellML metadata functionality using the RDFLib library."""
def __init__(self, name):
"""Create the wrapper."""
# Magic for pretending to be a module
self._module = sys.modules[name]
sys.modules[name] = self
self._initializing = True
# Map from cellml_model instances to RDF stores
self._models = {}
# Oxford metadata will be loaded lazily
self._metadata_names = self._stimulus_names = None
# Cope with differences in API between library versions
rdflib_major_version = int(rdflib.__version__[0])
if rdflib_major_version >= 3:
self.Graph = rdflib.Graph
else:
self.Graph = rdflib.ConjunctiveGraph
def __getattribute__(self, name):
"""Provide access to real module-level variables as though they're class properties."""
# call module.__init__ after import introspection is done
baseget = super(RdfProcessor, self).__getattribute__
module = baseget('_module')
if baseget('_initializing') and not name[:2] == '__' == name[-2:]:
setattr(self, '_initializing', False)
__init__(module)
try:
return baseget(name)
except AttributeError:
return getattr(module, name)
def _debug(*args):
pycml.DEBUG('cellml-metadata', *args)
def _load_ontology(self):
"""Load the Oxford metadata ontology the first time it's needed."""
pycml_path = os.path.dirname(os.path.realpath(__file__))
oxmeta_ttl = os.path.join(pycml_path, 'oxford-metadata.ttl')
oxmeta_rdf = os.path.join(pycml_path, 'oxford-metadata.rdf')
g = self.Graph()
# We allow a difference in modification time of 10s, so we don't get confused when checking out!
if os.stat(oxmeta_ttl).st_mtime > os.stat(oxmeta_rdf).st_mtime + 10.0:
# Try to regenerate RDF/XML version of ontology
try:
g.parse(oxmeta_ttl, format='turtle')
except Exception, e:
print >> sys.stderr, 'Unable to convert metadata from Turtle format to RDF/XML.'
print >> sys.stderr, 'Probably you need to upgrade rdflib to version 4.\nDetails of error:'
raise
g.serialize(oxmeta_rdf, format='xml')
else:
# Just parse the RDF/XML version
g.parse(oxmeta_rdf, format='xml')
annotation_terms = list(g.subjects(rdflib.RDF.type, rdflib.URIRef(pycml.NSS['oxmeta']+u'Annotation')))
self._metadata_names = frozenset(map(lambda node: self.namespace_member(node, pycml.NSS['oxmeta']), annotation_terms))
# Parameters for the stimulus current
self._stimulus_names = frozenset(filter(lambda name: name.startswith('membrane_stimulus_current_'), self._metadata_names))
@property
def METADATA_NAMES(self):
"""Fake a module-level constant as a property for lazy loading."""
if self._metadata_names is None:
self._load_ontology()
return self._metadata_names
@property
def STIMULUS_NAMES(self):
"""Fake a module-level constant as a property for lazy loading."""
if self._stimulus_names is None:
self._load_ontology()
return self._stimulus_names
def _create_new_store(self, cellml_model):
"""Create a new RDF store for the given CellML model.
The new store will be available as self._models[cellml_model].
"""
self._models[cellml_model] = self.Graph()
def _add_rdf_element(self, cellml_model, rdf_text):
"""Add statements to the model's graph from the given serialized RDF."""
g = self.Graph()
g.parse(StringIO(rdf_text))
rdf_model = self._models[cellml_model]
for stmt in g:
rdf_model.add(stmt)
def _serialize(self, cellml_model):
"""Serialize the RDF model for this CellML model to XML."""
return self._models[cellml_model].serialize()
def get_rdf_from_model(self, cellml_model):
"""Get the RDF graph of the given CellML model.
If this model is already in our map, return the existing RDF store.
Otherwise, extract metadata from all RDF elements in the cellml_model,
create a new RDF graph from these, and delete the original elements.
"""
if not cellml_model in self._models:
rdf_blocks = cellml_model.xml_xpath(u'//rdf:RDF')
self._create_new_store(cellml_model)
for rdf_block in rdf_blocks:
rdf_text = rdf_block.xml()
self._add_rdf_element(cellml_model, rdf_text)
rdf_block.xml_parent.xml_remove_child(rdf_block)
return self._models[cellml_model]
def remove_model(self, cellml_model):
"""The given model is being deleted / no longer needed."""
if cellml_model in self._models:
del self._models[cellml_model]
self._debug('Clearing RDF state for model', cellml_model.name)
def update_serialized_rdf(self, cellml_model):
"""Ensure the RDF serialized into the given CellML model is up-to-date.
If we have done any metadata processing on the given model, will serialize
our RDF store into the rdf:RDF element child of the model.
"""
if cellml_model in self._models:
# Paranoia: ensure it doesn't already contain serialized RDF
rdf_blocks = cellml_model.xml_xpath(u'//rdf:RDF')
if rdf_blocks:
pycml.LOG('cellml-metadata', logging.WARNING, 'Removing existing RDF in model.')
for rdf_block in rdf_blocks:
rdf_block.xml_parent.xml_remove_child(rdf_block)
# Serialize the RDF model into cellml_model.RDF
rdf_text = self._serialize(cellml_model)
rdf_doc = pycml.amara.parse(rdf_text)
cellml_model.xml_append(rdf_doc.RDF)
# Remove the RDF model
self.remove_model(cellml_model)
def create_rdf_node(self, node_content=None, fragment_id=None):
"""Create an RDF node.
node_content, if given, must either be a tuple (qname, namespace_uri),
or a string, in which case it is interpreted as a literal RDF node.
Alternatively, fragment_id may be given to refer to a cmeta:id within the
current model.
If neither are given, a blank node is created.
"""
if fragment_id:
node = rdflib.URIRef(str('#'+fragment_id))
elif node_content:
if type(node_content) == types.TupleType:
qname, nsuri = node_content
if nsuri[-1] not in ['#', '/']:
nsuri = nsuri + '#'
ns = rdflib.Namespace(nsuri)
prefix, local_name = pycml.SplitQName(qname)
node = ns[local_name]
elif type(node_content) in types.StringTypes:
node = rdflib.Literal(node_content)
else:
raise ValueError("Don't know how to make a node from " + str(node_content)
+ " of type " + type(node_content))
else:
node = rdflib.BNode()
return node
def create_unique_id(self, cellml_model, base_id):
"""Create a fragment identifier that hasn't already been used.
If base_id hasn't been used, it will be returned. Otherwise, underscores will
be added until a unique id is obtained.
"""
while True:
node = self.create_rdf_node(fragment_id=base_id)
if not self.get_targets(cellml_model, node, None):
break
base_id += u'_'
return base_id
def add_statement(self, cellml_model, source, property, target):
"""Add a statement to the model."""
self._debug("add_statement(", source, ",", property, ",", target, ")")
rdf_model = self.get_rdf_from_model(cellml_model)
rdf_model.add((source, property, target))
def replace_statement(self, cellml_model, source, property, target):
"""Add a statement to the model, avoiding duplicates.
Any existing statements with the same source and property will first be removed.
"""
self._debug("replace_statement(", source, ",", property, ",", target, ")")
rdf_model = self.get_rdf_from_model(cellml_model)
rdf_model.set((source, property, target))
def remove_statements(self, cellml_model, source, property, target):
"""Remove all statements matching (source,property,target).
Any of these may be None to match anything.
"""
self._debug("remove_statements(", source, ",", property, ",", target, ")")
rdf_model = self.get_rdf_from_model(cellml_model)
rdf_model.remove((source, property, target))
def get_target(self, cellml_model, source, property):
"""Get the target of property from source.
Returns None if no such target exists. Throws if there is more than one match.
If the target is a literal node, returns its string value. Otherwise returns an RDF node.
"""
rdf_model = self.get_rdf_from_model(cellml_model)
try:
target = rdf_model.value(subject=source, predicate=property, any=False)
except rdflib.exceptions.UniquenessError:
raise ValueError("Too many targets for source " + str(source) + " and property " + str(property))
if isinstance(target, rdflib.Literal):
target = str(target)
self._debug("get_target(", source, ",", property, ") -> ", "'" + str(target) + "'")
return target
def get_targets(self, cellml_model, source, property):
"""Get a list of all targets of property from source.
If no such targets exist, returns an empty list.
If property is None, targets of any property will be returned.
Alternatively if source is None, targets of the given property from any source will be found.
For each target, if it is a literal node then its string value is given.
Otherwise the list will contain an RDF node.
"""
rdf_model = self.get_rdf_from_model(cellml_model)
targets = list(rdf_model.objects(subject=source, predicate=property))
for i, target in enumerate(targets):
if isinstance(target, rdflib.Literal):
targets[i] = str(target)
return targets
def find_variables(self, cellml_model, property, value=None):
"""Find variables in the cellml_model with the given property, and optionally value.
property (and value if given) should be a suitable input for create_rdf_node.
Will return a list of cellml_variable instances.
"""
self._debug("find_variables(", property, ",", value, ")")
rdf_model = self.get_rdf_from_model(cellml_model)
property = self.create_rdf_node(property)
if value:
value = self.create_rdf_node(value)
vars = []
for result in rdf_model.subjects(property, value):
assert isinstance(result, rdflib.URIRef), "Non-resource annotated."
uri = str(result)
assert uri[0] == '#', "Annotation found on non-local URI"
var_id = uri[1:] # Strip '#'
var_objs = cellml_model.xml_xpath(u'*/cml:variable[@cmeta:id="%s"]' % var_id)
assert len(var_objs) == 1, "Didn't find a unique variable with ID " + var_id
vars.append(var_objs[0])
return vars
def get_all_rdf(self, cellml_model):
"""Return an iterator over all RDF triples in the model."""
rdf_model = self.get_rdf_from_model(cellml_model)
for triple in rdf_model:
yield triple
def namespace_member(self, node, nsuri, not_uri_ok=False, wrong_ns_ok=False):
"""Given a URI reference RDF node and namespace URI, return the local part.
Will raise an exception if node is not a URI reference unless not_uri_ok is True.
Will raise an exception if the node doesn't live in the given namespace, unless
wrong_ns_ok is True. In both cases, if the error is suppressed the empty string
will be returned instead.
"""
local_part = ""
if not isinstance(node, rdflib.URIRef):
if not not_uri_ok:
raise ValueError("Cannot extract namespace member for a non-URI RDF node.")
if node.startswith(nsuri):
local_part = node[len(nsuri):]
elif not wrong_ns_ok:
raise ValueError("Node is not in correct namespace.")
self._debug("namespace_member(", node, ",", nsuri, ") = ", local_part)
return local_part
####################################################################################
# Instantiate a processor instance that pretends to be this module
####################################################################################
p = RdfProcessor(__name__)
if __name__ == '__main__':
# Just load the ontology to trigger TTL -> RDF/XML conversion if needed
p._load_ontology()
# Proper enums in Python
# Taken from http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/413486
# Submitter: Zoran Isailovski
# Version: 1.6
# Updated: 2005/05/09
# License: http://www.python.org/psf/license/
def Enum(*names):
##assert names, "Empty enums are not supported" # <- Don't like empty enums? Uncomment!
class EnumClass(object):
__slots__ = names
def __iter__(self): return iter(constants)
def __len__(self): return len(constants)
def __getitem__(self, i): return constants[i]
def __repr__(self): return 'Enum' + str(names)
def __str__(self): return 'enum ' + str(constants)
class EnumValue(object):
__slots__ = ('__value')
def __init__(self, value): self.__value = value
Value = property(lambda self: self.__value)
EnumType = property(lambda self: EnumType)
def __hash__(self): return hash(self.__value)
def __cmp__(self, other):
# C fans might want to remove the following assertion
# to make all enums comparable by ordinal value {;))
assert self.EnumType is other.EnumType, "Only values from the same enum are comparable"
return cmp(self.__value, other.__value)
def __invert__(self): return constants[maximum - self.__value]
def __nonzero__(self): return bool(self.__value)
def __repr__(self): return str(names[self.__value])
maximum = len(names) - 1
constants = [None] * len(names)
for i, each in enumerate(names):
val = EnumValue(i)
setattr(EnumClass, each, val)
constants[i] = val
constants = tuple(constants)
EnumType = EnumClass()
return EnumType
if __name__ == '__main__':
print '\n*** Enum Demo ***'
print '--- Days of week ---'
Days = Enum('Mo', 'Tu', 'We', 'Th', 'Fr', 'Sa', 'Su')
print Days
print Days.Mo
print Days.Fr
print Days.Mo < Days.Fr
print list(Days)
for each in Days:
print 'Day:', each
print '--- Yes/No ---'
Confirmation = Enum('No', 'Yes')
answer = Confirmation.No
print 'Your answer is not', ~answer
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -4,7 +4,7 @@ SET(PP_SOURCES ./Prepacing.cpp
../../CellModels/CellModel.cpp
../../CellModels/CourtemancheRamirezNattel98.cpp
../../CellModels/FentonKarma.cpp
../../CellModels/TenTusscher06.cpp
../../CellModels/TenTusscher06.cpp
../../Stimuli/Stimulus.cpp
../../Stimuli/StimulusPoint.cpp
../../Stimuli/Protocol.cpp
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment