RuleWorld / nfsim

A general-purpose, stochastic, biochemical reaction simulator for large reaction networks
http://michaelsneddon.net/nfsim/
MIT License
14 stars 9 forks source link

NFSim passes over the wrong subpattern for local function evaluation #13

Closed jjtapia closed 8 years ago

jjtapia commented 8 years ago

IN the example linked here

https://gist.github.com/jjtapia/23bf8426b8a3c4ae15fdd7e3da18951b

NFSim reports

Internal error in LocalFunction::evaluateOn()! trying to evaluate a function with unknown scope. It seems the problem is that NFSim is passing the wrong molecule subpattern for local function evaluaiton, so in the rule

Bug1: %z:X(a!1,a!+).Y(b!1,d~0)->X(a!1,a!+).Y(b!1,d~1) factor(z)*kmod

For whatever reason it is evaluating Ytot on the X molecule, which produces the wrong scope error we are getting. Funny enough, if you change the rule to

NowItWorks1: %z:Y(b!1,d~0).X(a!1,a!+)->Y(b!1,d~1).X(a!1,a!+) factor(z)*kmod

NFsim will pass the correct 'Y' subpattern. I think the error is happening during the DORReaction constructor function where it creates the wrong mapping information on the 'argIndexIntoMappingSet' array. This happens both in NFSim v1.11 and v1.12. The solution would then to take better care on which pattern is passed for evaluation during the mapping of a species pattern to a local function.

It is also work looking into how this relates to validation example v17.bngl. I think they are different issues but it's worth looking over it.

jjtapia commented 8 years ago

I looked a little bit more into this issue and unfortunately the fix does not seem to be simple. The issue has to do with how local functions using Species labels are evaluated in DORRxnClass:evaluateLocalFunctions, specifically this line

this->argMappedMolecule[i] = ms->get(this->argIndexIntoMappingSet[i])->getMolecule();

tells nfsim which molecule agent will be sent for evaluation into a local function, like this

double value = this->cf->evaluateOn(argMappedMolecule,argScope, reactantCounts, n_reactants);

Which molecular pattern is selected to be sent for evaluation (that is, the molecules that argMappedMolecule is mapping) is not chosen based on any labeling scheme, but based on which is the 'root molecule type' of sorts for the mappings contained inside this->argIndexIntoMappingSet. This 'root molecule type' is chosen dependant on how a complex is written.

in X(a!1,a!+).Y(b!1,d~0)

X is selected as the root molecule type

in Y(b!1,d~0).X(a!1,a!+)

Y is selected as the root molecule type

The issue is that when this so-called root molecule types in this case is chosen to be 'X' in the case reported as a bug. When it is compared against Ytot nfsim complains that it doesn't know what to do with it. The real solution would be to make it so that nfsim just doesnt send a mapping's root molecule type for evaluation to a local function but actually selects the right one. (I don't know why it is done this way anyway, but this is the bug right here. when the intersection between the root molecule agent and the observable in question is empty nfsim balks). Alternatively keep walking down a complex till you find something you can evaluate. (this would be simpler to implement although potentially very expensive).

The solution I suggest for now is that when you have a pattern that only partially matches against an observable and you can choose an specific molecule type instead of the whole species you are more specific with the labeling, like so

X(a!1,a!+).Y(b!1,d~0)%z->X(a!1,a!+).Y(b!1,d~1)%z factor(z)*kmod

And avoid species labeling if possible.

jjtapia commented 8 years ago

Addressed in commit ae4237d5d7cdedccc35bafc3f99c7aa4fec3eec8

jjtapia commented 8 years ago

So I opted for the second solution I mentioned, where once it finds a species label-local function parameter matching problem it falls back to trying other molecule patterns included in a reaction's mapping object. I included a test case indexed v19.bngl.