Invalid steady-states in Scanner.

[exe0cdc]

In the Run method of Scanner there is code that deals with invalid steady-states using the __StateOK__ flag of PysMod. Unfortunately it only works when the Scanner _MODE_ flag is set to 'state'.

As soon as an invalid state is reached with any other of the modes ('mca', 'elasticity'), RUN exits due to the self.__StateOK__ == True assertion in doMca and doElas.

This can be fixed at the doMca or doElas method level by adding an optional parameter (or setting) that allows the method to ignore the assertion and just skip methods coming after the assertion, thus having the same effect without causing an assertion error. Run will then handle these invalid steady-states in the same way as those produced by doState (as far as I can see).

OR

It can be fixed by running doState in the Analyse method of Scanner regardless of _MODE_ and then just run additional methods as needed and only when __StateOK__ is True.

I'm happy to implement either of these changes depending on which one is decided to be better and if anyone wants this change. Let me know if I have missed anything.

[bgoli]

I would go with something like the second option, except that it incurrs a performance penalty as you run a steady-state solve twice right? I would try and deal with the assertion error at the Analyse level, consider the following:

"""this is the doX method where the solver is failing and raising an error""" def doError(b=1): assert b == False

"""deal with it in Analyse""" try: doError() except AssertionError: print('I am dealing with a steady-state solve error')

instead of wrapping the whole of analyse in a try/except block I would set a boolean control variable "doneState" and use that to run the analyses or not. You will also have to deal with the output, insert nan's or something. Let me know what you think I can also look into it next week or so,

[exe0cdc]

It is possible to do it without solving twice. It could look something like this:

def Analyze(self):
    """
    The analysis method, the mode is automatically set by the self.addUserOutput()
method but can be reset by the user.
    """
    # State is run regardless of mode (except 'null')
    valid_state = False
    if self.__MODE__ !='null':
        self.mod.State()
        valid_state = self.mod.__StateOK__
    # Now to avoid redundancy the methods that each of doElas, doMca and 
    # doEigenMca call (except for State) are called in each conditional.  
    # the steady state must now also be valid before they are called. 
    elif self._MODE_ == 'elasticity' and valid_state:
        self.mod.EvalEvar()
        self.mod.EvalEpar()
    elif self._MODE_ == 'mca'  and valid_state:
        self.mod.EvalEvar()
        self.mod.EvalEpar()
        self.mod.EvalCC()
    elif self._MODE_ == 'stability'  and valid_state:
        self.mod.EvalEvar()
        self.mod.EvalCC()
        self.mod.EvalEigen()
    elif self._MODE_ == 'null' :
        self.HAS_STATE_OUTPUT = False  
    self.mod.User_Function() #What does this do? The docstring says 
    # deprecated and internally it does exec(self.__CODE_User) 

and then StoreData can be modified to store NaNs if a steady state is not reached

def StoreData(self):
    """
    Internal function which concatenates and stores the data generated by Analyze.
    """
    valid_state = self.mod.__StateOK__
    if self.HAS_STATE_OUTPUT and self._MODE_ != 'null':
        if valid_state:
            self.SteadyStateResults.append(scipy.hstack((self.mod.state_species,self.mod.state_flux)))
        else:
            # this may be redundant due to __settings__["mode_state_nan_on_fail"]
            self.SteadyStateResults.append(scipy.hstack(([scipy.NaN for _ in self.mod.state_species],[scipy.NaN for _ in self.mod.state_flux])))    
    if self.HAS_USER_OUTPUT and self._MODE_ != 'null':
        if valid_state:
            self.UserOutputResults.append(scipy.array([getattr(self.mod,res) for res in self.UserOutputList]))
        else:
            self.UserOutputResults.append(scipy.array([scipy.NaN for _ in self.UserOutputList]))

I think by making these changes it should work. (I've just realised that I might as well have implemented them to see for myself, but I'm off to bed - I'll test it in the morning)

[bgoli]

Try it but you will be duplicating the doX functions what about something like this:

    def Analyze(self):
        """
        The analysis method, the mode is automatically set by the self.addUserOutput()
        method but can be reset by the user.
        """
        try:
            if self._MODE_ == 'state':
                self.mod.doState()
            elif self._MODE_ == 'elasticity':
                self.mod.doElas()
            elif self._MODE_ == 'mca':
                self.mod.doMca()
            elif self._MODE_ == 'stability':
                self.mod.doEigenMca()
            elif self._MODE_ == 'null':
                self.HAS_STATE_OUTPUT = False
        except AssertionError:
            self.HAS_STATE_OUTPUT = False
        self.mod.User_Function()

[exe0cdc]

I'm coming back to this issue after more than a year, because previously I just found a workaround for my use case and never really needed to edit this code. However, while refactoring my own code I have run into it again, and I want to fix it for good. But first I want to work through a few things:

When PySCeS solves the steady-state via the state method of State it looks to the setting 'mode_state_nan_on_fail' to determine if an invalid state variable should be replaced with NaN or not. This behaviour is not found in the EvalEvar, EvalEpar, EvalCC, etc. methods, and rather an AssertionError is raised when __StateOK__ != True in the doMca method (and others) .

So when Scanner runs, the invalid steady state results for the 'state' mode is handled at the State method level while the other modes just fail on invalid steady states due to the assertion.

Now, as previously suggested, the problem can be fixed at the Analyse and StoreData level (similar to above) as follows:

def Analyze(self):
        """
        The analysis method, the mode is automatically set by the self.addUserOutput()
        method but can be reset by the user.
        """
        try:
            if self._MODE_ == 'state':
                self.mod.doState()
            elif self._MODE_ == 'elasticity':
                self.mod.doElas()
            elif self._MODE_ == 'mca':
                self.mod.doMca()
            elif self._MODE_ == 'stability':
                self.mod.doEigenMca()
            elif self._MODE_ == 'null':
                self.HAS_STATE_OUTPUT = False
        except AssertionError:
            # Only when _MODE_ = 'elasticity', 'mca' or 'stability' and
            # an invalid state is reached will this AssertaionError be 
            # raised. It should be passed silently as the issue will be
            # dealt with in StoreData
            pass
        self.mod.User_Function()

    def StoreData(self):
        """
        Internal function which concatenates and stores the data generated by Analyze.
        """
        # it seems that checking for both HAS_STATE_OUTPUT and _MODE_!='null' is 
        # redundant as the only place where HAS_STATE_OUTPUT can be set to False
        # is when _MODE_='null' - correct me if I am wrong
        if self.HAS_STATE_OUTPUT and self._MODE_ != 'null': 
            # invalid states are handled on the State level
            # either valid states, invalid states (depending on mode_nan_on_fail)
            # or NaNs will be added to SteadyStateResults
            self.SteadyStateResults.append(scipy.hstack((self.mod.state_species,self.mod.state_flux)))
        if self.HAS_USER_OUTPUT and self._MODE_ != 'null':
            # this should give the same result as above 
            if self.mod.__StateOK__ or not self.mod.__settings__["mode_state_nan_on_fail"]:
                self.UserOutputResults.append(scipy.array([getattr(self.mod,res) for res in self.UserOutputList]))
            elif self.mod.__settings__["mode_state_nan_on_fail"]:
                self.UserOutputResults.append(scipy.array([scipy.NaN for res in self.UserOutputList]))

OR

It can be fixed on the EvalEvar, EvalEpar etc. level by implementing behaviour similar to State. This can be done by removing the assertion in the higher level functions and adding code to check for the 'mode_state_nan_on_fail' setting.

The first option has the least side effects, but the second one is more consistent with the behaviour of State and will mean that no changes to Scanner will be needed. However, I do not know if it will even work, i.e., is it possible to calculate elasticities, control coefficients etc., for invalid steady states (even though these values will obviously also be invalid)? Secondly should access to invalid steady states even be an option? Why does State have the option to return invalid steady states, such as negative concentrations, in the first place (this is actually the default behaviour)? I am sorry for my ignorance regarding this final point, but I am truly curious why PySCeS is designed this way.

EDIT: I have tested the first solution and it works as expected (at least for invalid elasticity coefficients).

[bgoli]

Hi. PySCeS behaves like for three reasons. When it was coded not all SciPy algorithms could handle NaN values, watching where negative concentrations occur can give you useful information on structural errors in your model and it was sometimes useful to let Pitcon initialise on mathematically correct solutions.

The point is that somehow you have to let the user know that the solution is (biologically) invalid, either by raising an error and stopping the calculation or by using NaN's which should just result in all output being NaN.

So I tend to lean towards a solution which raises an error, but returns a results for diagnostic purposes or returns NaN's. I think this is in line with your first solution, why don't you just branch and implement it and we can test how it works.

[bgoli]

Closing as I think this suggestion was implemented in newer releases