2713 - Chaste_codegen Add a gating variable annotation to the metadata about chaste

[MauriceHendrix]

Comment by @mirams:
When you do a normal Chaste build it makes ones like this:

    void CellLuoRudy1991FromCellMLCvode::VerifyStateVariables()
    {
        /* We only expect CVODE to keep state variables to within its tolerances,
         * not exactly the bounds prescribed to each variable that are checked here.
         *
         * For 99.99% of paces this->mAbsTol works,
         * For 99.999% of paces 10*this->mAbsTol is fine,
         * but unfortunately 100x seems to be required on rare occasions for upstrokes.
         * This sounds bad, but is probably typically only 1e-5 or 1e-6.
         */
        const double tol = 100*this->mAbsTol;
        N_Vector rY = rGetStateVariables();
        double var_chaste_interface__fast_sodium_current_m_gate__m = NV_Ith_S(rY, 2);
        // Units: dimensionless; Initial value: 0.00187018
        
        if (var_chaste_interface__fast_sodium_current_m_gate__m < 0.0 - tol || var_chaste_interface__fast_sodium_current_m_gate__m > 1.0 + tol)
        {
            EXCEPTION(DumpState("State variable fast_sodium_current_m_gate__m has gone out of range. Check numerical parameters, for example time and space stepsizes, and/or solver tolerances"));
        }
        
    }

for CVODE (which is "allowed" to go a little bit outside physical bounds as long as its within tolerances according to its manual)
and like this

    void CellLuoRudy1991FromCellML::VerifyStateVariables()
    {
        std::vector<double>& rY = rGetStateVariables();
        double var_chaste_interface__fast_sodium_current_m_gate__m = rY[2];
        // Units: dimensionless; Initial value: 0.00187018
        
        if (var_chaste_interface__fast_sodium_current_m_gate__m < 0.0 || var_chaste_interface__fast_sodium_current_m_gate__m > 1.0)
        {
            EXCEPTION(DumpState("State variable fast_sodium_current_m_gate__m has gone out of range. Check numerical parameters, for example time and space stepsizes, and/or solver tolerances"));
        }
        
    }

for non-CVODE cells

I can't remember what it is that was tagged in that model to produce that like that, but it is the one in the main Chaste source that is producing that (via codegen already :)

So definitely an Exception for the non-CVODE cells, and worth reading the CVODE documentation to see if there is anything about the "recoverable error flag" instead that Jon mentioned

[MauriceHendrix]

"The SUNDIALS solvers allow for the RHS/RES functions to return a positive flag indicating a "recoverable error"." https://computing.llnl.gov/projects/sundials/usage-notes
sounds like it's something the solving step should return not something you set.

[MauriceHendrix]

comment by @mirams
Perhaps https://chaste.cs.ox.ac.uk/trac/ticket/2715 gives a nasty ODE system (you will need to login to trac to see my project svn where it lives probably)

[MauriceHendrix]

right so I can build & run the system & put some loops in to run the various conditions but the numbers I get don't look like the table in the ticket.
my numbers:

Scaling Factor	CVODE No Jacobian	CVODE With Jacobian
1.0	2.80875e-05	4.26515e-07
1.5	0.00175806	0.00132057
2.0	0.00175806	0.000294663
2.5	0.0211667	0.00170412
3.0	CV_ERR_FAILURE	CV_ERR_FAILURE
3.5	CV_ERR_FAILURE	CV_ERR_FAILURE

When the RHS returns a recoverable error flag I get:

Scaling Factor	CVODE No Jacobian	CVODE With Jacobian
1.0	3.55766e-05	4.76901e-07
1.5	0.000312776	CV_REPTD_RHSFUNC_ERR
2.0	CV_REPTD_RHSFUNC_ERR	0.000294663
2.5	CV_REPTD_RHSFUNC_ERR	CV_REPTD_RHSFUNC_ERR
3.0	CV_REPTD_RHSFUNC_ERR	CV_REPTD_RHSFUNC_ERR
3.5	CV_REPTD_RHSFUNC_ERR	CV_REPTD_RHSFUNC_ERR

So it appears to not help and we decided to park teh issue for now.

[mirams]

Yeah, parking the recoverable flag is a bit different to the gating variable annotations existing and generating useful things in the VerifyStateVariables method though, still useful to have, particularly for non-CVODE runs. We probably still want to be able to do that. (We may have done this already and I have forgotten).

But the idea would be to have some metadata structure like

HodgkinHuxleyGatingVariable "is a" Probability
Markov Model State "is a" Probability

and then any variable that is a "Probability" automatically gets [0,1] or for CVODE [0-tol,1+tol] checks in the VerifyStateVartiables() method.

[mirams]

see also the comment on concentrations in the ticket description, are they all getting >0 checks at the moment?

[MauriceHendrix]

see also the comment on concentrations in the ticket description, are they all getting >0 checks at the moment?

no that's indeed a seperate thing to the recoverable error as you say (re-opened ticket)

[mirams]

Hmmm, reading through this, not sure it is all finished (reopening to check!)

[mirams]

OK, yeah, the generated C++ files get these in the .hpp

private:
const bool is_concentration[8] = {false, true, false, false, false, false, false, false};
const bool is_probability[8] = {false, false, false, false, false, false, false, false};

and a method like this in the .cpp

    void CellDiFrancescoNoble1985FromCellML::VerifyStateVariables()
    {
        std::vector<double>& rY = rGetStateVariables();
        std::string error_message = "";
        
        for (unsigned i=0; i < 16; i++)
        {
            if(std::isnan(rY[i]))
            {
                error_message += "State variable " + this->rGetStateVariableNames()[i] + " is not a number\n";
            }
            if(std::isinf(rY[i]))
            {
                error_message += "State variable " + this->rGetStateVariableNames()[i] + " has become INFINITE\n";
            }
            if(this->is_concentration[i] && rY[i] < 0)
            {
                error_message += "Concentration " + this->rGetStateVariableNames()[i] + " below 0\n";
            }
            if(this->is_probability[i] && rY[i] < 0)
            {
                error_message += "Probability " + this->rGetStateVariableNames()[i] + " below 0\n";
            }
            if(this->is_probability[i] && rY[i] > 1)
            {
                error_message += "Probability " + this->rGetStateVariableNames()[i] + " above 1\n";
            }
        }
        if (error_message != ""){
            EXCEPTION(DumpState(error_message));
        }
    }

if it is showing a false when it should be true for a certain state variable, it is probably just a matter of tagging the CellML appropriately with ontology terms