Type error in SSA after substitution of rates

[tbose1]

The error message is TypeError: cannot determine truth value of Relational

Attached is a minimal example producing the error:

SSAerror.ipynb.zip

In this example ssa1 works whereas ssa2 produces the error message

[jarmarshall]

This seems to have the same root cause as issue #321 and originates in substitute(). See attached notebook
substitution_error.ipynb.zip

[jarmarshall]

Updated notebook showing that it is the equations themselves that are truncated by substitute() in particular circumstances (to be determined)

As a temporary workaround where this bug crops up, constructing models from scratch rather than using substitute() should work.

substitution_error.ipynb.zip

[jarmarshall]

This seems very likely to be an error in simpy's subs() function (at least under v1.1.1) - need to check if pull request #344 also closes this, otherwise raise an issue on the sympy GitHub...

[jarmarshall]

So this is not a Sympy error, but rather an unhandled side effect of substituting in rates that lead to simplification of the equation system (see sympy/sympy#17627) - need to figure out precisely where the logic is that leads to that happening and deal with it...

[jarmarshall]

Alternatively could use xreplace() instead of subs() to leave in all sliders, but leaving sliders in that have no effect on the equations would be confusing for users. A possible short-term fix though...

[jarmarshall]

Also using subs() with evaluate = False would have the same effect...

[jarmarshall]

Having now delved some way down into this I have concluded that this is an error around how SSA() is invoked and am asking @joefresna to take over...
The substitute() logic seems sound... in the example given where the two mass action terms cancel in the ODEs once the substitution r_A=r_B=r is given, it is appropriate to still include the reactions having rate r - these reactions are happening, they just cancel each other out in the mean field limit. Thus the rate r should still be associated with the model. Under the current substitute() logic it is. Furthermore setting the system size for model2 and calling stream() for example will work correctly - the redundant slider for r will be removed and the plot will function as normal.

The issue for SSA() is that the slider for r is not present - this causes an issue for SSA since then no value for it is available (whereas having no value for r is no problem in stream() since the rate has cancelled out in the ODE system). The fix is to ensure that for stochastic views only, rates for sliders are taken from model._rates. This makes sense also because we cannot know a priori which rates may affect the dynamics in finite populations, even if they do not affect the infinite population dynamics, and also the interactions are still meaningful. @joefresna if you'd look into making that change that would be helpful.

[jarmarshall]

N.B. branch missing-rate already created to work on this...

[joefresna]

It seems that all widgets are generated from model._rates already.

MuMoT/mumot/models.py

Line 1713 in 8dcd2c3

for freeParam in self._rates.union(self._constantReactants):

In fact, if I print the model2._rates in the notebook linked above it returns only {g}.

@jarmarshall , from which variable can I get the rate r?

[jarmarshall]

Can you dig around a bit for the best place to take it from? It's definitely in ._stochiometry but that needs some processing...

[joefresna]

@jarmarshall , I think that the only place where to find the missing rate is the stoichiometry...
091cb93 may fix the issue...

but my stoichiometry parser may fail for special cases which I did not test... hopefully not

[jarmarshall]

Thanks - do you want to make a pull request for review?

[joefresna]

yes, but before doing so I wanted to locally run some test notebook to see if there were errors, but I notice that the user-manual hangs at cell 19
stream1 = model4.stream('A','B', showFixedPoints = True, initWidgets={'mu':[3, 1, 5, 0.5]})
and it does not move forward.... it's now a couple of hours it's running at cell 19... I know I've an old laptop ;-) but maybe it's not the expected behaviour. Do you know anything about this, or it has been caused by my changes?

[jarmarshall]

If you’ve pulled all changes in from master then that is almost certainly down to your changes I think... however I don’t believe you should have changed any code that’s called there, so maybe restart the tests?

[joefresna]

My attempted patch is not correct.

I don't know where we can get the rates from...in case they must be taken from the stoichiometry it will require quite an articulated parser because sometimes they can get quite messy.

For instance, in my tests the run was hung up at model4. If I ask model4._stoichiometry I get

{Reaction 1: {'rate': Delta/2 + mu, U: [1, 0, {U: -A - B + N}], A: [0, 1]},
 Reaction 2: {'rate': -Delta/2 + mu, U: [1, 0, {U: -A - B + N}], B: [0, 1]},
 Reaction 3: {'rate': 1/(Delta/2 + mu), A: [1, 0], U: [0, 1, {U: -A - B + N}]},
 Reaction 4: {'rate': 1/(-Delta/2 + mu),
  B: [1, 0],
  U: [0, 1, {U: -A - B + N}]},
 Reaction 5: {'rate': Delta/2 + mu, A: [1, 2], U: [1, 0, {U: -A - B + N}]},
 Reaction 6: {'rate': -Delta/2 + mu, B: [1, 2], U: [1, 0, {U: -A - B + N}]},
 Reaction 7: {'rate': s, A: [1, 1], B: [1, 0], U: [0, 1, {U: -A - B + N}]},
 Reaction 8: {'rate': s, A: [1, 0], B: [1, 1], U: [0, 1, {U: -A - B + N}]}}

So, I don't know if we need to parse this stuff, or there could be an easier way to get the rates.

[jarmarshall]

I think you can use free_symbols() on the rate objects, then take a union over those sets?

[jarmarshall]

This is a kludge for now and we may think of a better way of doing it in the future, at a lower level in the base model class, e.g. keeping separate lists of rates from equations, and from interaction rules.

[joefresna]

it wasn't so easy but it should be now sorted... waiting for tests to finish.
I agree that we should create a new object named _allRates that keeps track of all rates of the model. It may require a bit of code to update such a list for substitute()