Parsing of realsitic models - from Model manipulation

[skeating]

Remaining items from this section of #21

• Check we can parse all the models in Chaste's CellML project
• We possibly don't want to check all the equations manually, but doing a few cases might be wise

@jonc125 Would this be solved by creating CellML models with various MathML in eg #80

[jonc125]

No, this is different from testing just basic MathML, it's testing some actual cardiac models. Again this gets implicitly tested by running full protocols on them, so it's just a question of degree as to how much we need. Perhaps worth checking at least hodgkin-huxley!?

[MichaelClerx]

One way to check is to evaluate the RHS and then compare the calculated derivatives with some reference tool. I think @tamuri did this at some point

[tamuri]

I think @tamuri did this at some point

I did, but looks like it's been removed.

[MichaelClerx]

Still here a few weeks ago:

cellmlmanip/tests/test_hodgkin.py

Lines 139 to 236 in 2419d8b

	def test_evaluation(self, model):
	"""
	From Michael:
	These are my calculations for the derivatives of the state variables
	in the hodgkin-huxley model, evaluated from the initial state (so the
	initial values of the state variables)
	Reading model from hodgkin_huxley_squid_axon_model_1952_modified.mmt...
	Model hodgkin_huxley_squid_axon_model_1952_modified read successfully.
	Evaluating state vector derivatives...
	-------------------------------------------------------------------------------
	Name Initial value Derivative at t=0
	-------------------------------------------------------------------------------
	membrane.V -7.50000000000000000e+01 -6.00768750000000740e-01
	sodium_channel_m_gate.m 5.00000000000000028e-02 1.23855383553985177e-02
	sodium_channel_h_gate.h 5.99999999999999978e-01 -4.55523906540064583e-04
	potassium_channel_n_gate.n 3.25000000000000011e-01 -1.34157228632045961e-03
	-------------------------------------------------------------------------------```
	The free variable time is set to 0
	Numbers are converted to string with `'{:<1.17e}'.format(number)`,
	which I think should give the maximum precision for any double precision float
	(as in, converting to and back from string shouldn't lose
	accuracy, it does _not_ mean all printed digits are accurate)
	"""
	graph = model.graph_with_sympy_numbers
	# initial values for the free and state variables
	initials = {
	'environment$time': 0.0,
	'membrane$V': -7.50000000000000000e+01,
	'sodium_channel_m_gate$m': 5.00000000000000028e-02,
	'sodium_channel_h_gate$h': 5.99999999999999978e-01,
	'potassium_channel_n_gate$n': 3.25000000000000011e-01
	}
	# the calculated derivatives at the initial conditions
	evaluated_derivatives = {
	'membrane$V': -6.00768750000000740e-01,
	'sodium_channel_m_gate$m': 1.23855383553985177e-02,
	'sodium_channel_h_gate$h': -4.55523906540064583e-04,
	'potassium_channel_n_gate$n': -1.34157228632045961e-03
	}
	sorted_symbols = nx.lexicographical_topological_sort(graph, key=str)
	# collects all the evaluated lh-sides of the equations
	evaluated_symbols = dict()
	# saves the evaluated lh-side for the state variables
	state_derivatives = dict()
	# loop over each of the equations in the model (topologically sorted)
	for symbol in sorted_symbols:
	equation = graph.nodes[symbol]['equation']
	# if this symbol is calculated using an equation
	if equation is not None:
	# substitute all symbols in the rhs of the equation
	eq_substituted = equation.rhs.subs(evaluated_symbols)
	# if any symbols remain, we have a problem
	remaining_symbols = eq_substituted.atoms(sympy.Symbol)
	if remaining_symbols:
	for remaining_symbol in remaining_symbols:
	pytest.fail("Unresolved symbol %s in %s" % (remaining_symbol, equation))
	# calculate the result
	eq_evaluated = eq_substituted.evalf()
	# save the calculation for this symbol (to be used in the next equations)
	evaluated_symbols[equation.lhs] = eq_evaluated
	# if the symbol on the lhs is a derivative
	if equation.lhs.is_Derivative:
	# save the calculation for testing
	state_derivatives[equation.lhs.free_symbols.pop()] = eq_evaluated
	# otherwise the symbol doesn't have an equation
	# if the symbol has an initial value
	elif symbol.name in initials:
	# add the symbol's initial value to the substitution dictionary
	evaluated_symbols[symbol] = sympy.Number(initials[symbol.name])
	else:
	# something has gone wrong - no equation or initial value
	pytest.fail("Symbol " + str(symbol) + " does not have equation or initial value.")
	# for each calculated state variable
	for state_symbol, evaluated_deriv in state_derivatives.items():
	# check evaluation against expected
	expected = evaluated_derivatives[state_symbol.name]
	actual = evaluated_deriv
	assert float(actual) == pytest.approx(expected)

[MichaelClerx]

Hmmm, yup disappeared in #116 .

Wonder if we should bring it back in some form, @jonc125 ?

[skeating]

So test_hodgkin disappeared when I removed all test_specific_model files and incorporated the tests in test_model_functions; since we had talked about not having test files named after specific models. I didn't change test coverage when I did this; so clearly the functionality tested here is tested elsewhere but I could see the point of having a test_evaluation file that included the above test for HH and then possibly similar tests for other models if we get round to that.

[jonc125]

I wouldn't worry about adding other models, but bringing that specific test for HH back is probably worthwhile (especially as it's easy to do!). Then we can close this issue.

[MauriceHendrix]

I wouldn't worry about adding other models, but bringing that specific test for HH back is probably worthwhile (especially as it's easy to do!). Then we can close this issue.

I think pretty much all "other models" will be implicitly tested because I'm testing their chaste codegen. If they don't parse with cellmlmanip we'll find out soon enough!

[MichaelClerx]

So test_hodgkin disappeared when I removed all test_specific_model files and incorporated the tests in test_model_functions; since we had talked about not having test files named after specific models. I didn't change test coverage when I did this; so clearly the functionality tested here is tested elsewhere

That makes sense!
Regarding coverage: We've talked about this a lot in PINTS, where we now have an informal rule to try and add value-based tests for every method. I.E. tests that actually check the output is correct.
I've heard about structured ways of doing this but they all sound like a lot of work / extra software.

Does anyone here have experience with those?

[jonc125]

Do you mean tools like hypothesis? Not directly, but I've seen them used.