CellMLToNektarTranslator.py

Go to the documentation of this file.

#Importing various functions used by this class (i have outcommented modules that we don't need)
import os
import optparse
import re
import time
import sys
from cStringIO import StringIO

# Common CellML processing stuff below
import pycml
from pycml import *  # Put contents in the local namespace as well
import optimize
import processors
import validator

#Importing functions from modified_translators.py
import translators
#from modified_translators import CellMLTranslator

class CellMLToNektarTranslator(translators.CellMLTranslator):
    """
    As CellMLTranslator, but targets more recent Chaste style.

    Includes the ability to output a cell that can solve itself using
    backward Euler, if the appropriate analyses have been done on the
    model.  (See the -J and -j options to translate.py.)
    """

    # We want separate .cpp/.hpp files
    USES_SUBSIDIARY_FILE = True

    # Type of (a reference to) the state variable vector
    TYPE_VECTOR = 'std::vector<double> '
    TYPE_VECTOR_REF = 'std::vector<double>& '
    
    NODESET = []

    def writeln_hpp(self, *args, **kwargs):
        """Convenience wrapper for writing to the header file."""
        kwargs['subsidiary'] = True
        self.writeln(*args, **kwargs)

    def translate(self, *args, **kwargs):
        """Generate code for the given model."""
        our_kwargs = {'use_chaste_stimulus': False,
                      'separate_lut_class': True,
                      'convert_interfaces': False,
                      'use_modifiers': False,
                      'use_data_clamp': False,
                      'dynamically_loadable': False,
                      'use_protocol': False
                      }
        for key, default in our_kwargs.iteritems():
            setattr(self, key, kwargs.get(key, default))
            if key in kwargs:
                del kwargs[key]
        # Some other default settings
        self.use_backward_euler = False
        self.include_serialization = False
        # Last method's access specification
        self._last_method_access = 'private'
        return super(CellMLToNektarTranslator, self).translate(*args, **kwargs)

    def final_configuration_hook(self):
        """Set the LT method prefix (requires self.class_name to be set)."""
        if self.separate_lut_class:
            self.lt_class_name = self.class_name + '_LookupTables'
            self.lookup_method_prefix = self.lt_class_name + '::Instance()->'
        return super(CellMLToNektarTranslator, self).final_configuration_hook()
        
    def output_includes(self, base_class=None):
        """Output the start of each output file.
        
        As well as the #include lines, it also outputs the include guard for
        the .hpp file, and doxygen comment.
        
        If base_class is not None (and self.use_backward_euler isn't set)
        then includes that class' header instead of AbstractCardiacCell.
        
        If self.dynamically_loadable is set, includes extra headers needed
        for that case.
        
        Reads self.include_serialization and self.use_backward_euler.
        Sets self.base_class_name and self.class_inheritance.
        """

        
        from translators import version_comment
        for sub in [False, True]:
            self.output_doxygen('@file\n\n',
                                'This source file was generated from CellML.\n\n',
                                'Model: ', self.model.name, '\n\n',
                                version_comment(self.add_timestamp),
                                '\n\n<autogenerated>',
                                subsidiary=sub)
            self.writeln(subsidiary=sub)


        self.writeln('#include <iostream>')
        self.writeln('#include <string>')
            
        # .cpp should include .hpp
        #writes path
        self.writeln('#include <CardiacEPSolver/CellModels/', os.path.basename(self.subsidiary_filename), '>')
        # if self.include_serialization:
        #     self.writeln_hpp('#include "ChasteSerialization.hpp"')
        #     self.writeln_hpp('#include <boost/serialization/base_object.hpp>')
        #self.writeln('#include <cmath>')
        #self.writeln('#include <cassert>')
        #self.writeln('#include <memory>')

        #write the namespace (and the open bracket)
        self.writeln()
        self.writeln('namespace Nektar')
        self.writeln('{')


        if self.use_backward_euler:
            self.writeln_hpp('#include "AbstractBackwardEulerCardiacCell.hpp"')
            self.writeln('#include "CardiacNewtonSolver.hpp"')
            self.base_class_name = 'AbstractBackwardEulerCardiacCell<' + \
                str(self.nonlinear_system_size) + '>'
        elif self.options.rush_larsen:
            self.base_class_name = 'AbstractRushLarsenCardiacCell'
            self.writeln_hpp('#include "' + self.base_class_name + '.hpp"')
            if not self.doc._cml_rush_larsen:
                self.writeln('#include "Warnings.hpp"')
        elif self.options.grl1:
            self.base_class_name = 'AbstractGeneralizedRushLarsenCardiacCell'
            self.writeln_hpp('#include "' + self.base_class_name + '.hpp"')
        elif self.options.grl2: #1992 TODO: merge with above case
            self.base_class_name = 'AbstractGeneralizedRushLarsenCardiacCell'
            self.writeln_hpp('#include "' + self.base_class_name + '.hpp"')
        elif base_class:
            self.base_class_name = base_class
            self.writeln_hpp('#include "' + self.base_class_name + '.hpp"')
        else:
            self.base_class_name = 'AbstractCardiacCell'
            # self.writeln_hpp('#include "' + self.base_class_name + '.hpp"')
        if self.use_modifiers:
            self.writeln_hpp('#include "AbstractCardiacCellWithModifiers.hpp"')
            self.writeln_hpp('#include "AbstractModifier.hpp"')
            # Modify the base class name
            self.base_class_name = 'AbstractCardiacCellWithModifiers<' + self.base_class_name + ' >'
        self.class_inheritance = ' : public CellModel'# + self.base_class_name
        if self.dynamically_loadable:
            self.writeln_hpp('#include "AbstractDynamicallyLoadableEntity.hpp"')
            self.class_inheritance += ', public AbstractDynamicallyLoadableEntity'
        if self.use_protocol:
            self.writeln_hpp('#include "AbstractTemplatedSystemWithOutputs.hpp"')
            self.class_inheritance += ', public AbstractTemplatedSystemWithOutputs<' + self.TYPE_VECTOR + '>'
        #self.writeln('#include "Exception.hpp"')
        #self.writeln('#include "OdeSystemInformation.hpp"')
        #self.writeln('#include "RegularStimulus.hpp"')
        # self.writeln_hpp('#include "AbstractStimulusFunction.hpp"')
        #self.writeln('#include "HeartConfig.hpp"')
        #self.writeln('#include "IsNan.hpp"')
        #self.writeln('#include "MathsCustomFunctions.hpp"')
        #self.writeln()

        self.writeln_hpp('#ifndef NEKTAR_SOLVERS_ADRSOLVER_EQUATIONSYSTEMS_', self.include_guard)
        self.writeln_hpp('#define NEKTAR_SOLVERS_ADRSOLVER_EQUATIONSYSTEMS_', self.include_guard, '\n')

        self.writeln_hpp('#include <CardiacEPSolver/CellModels/CellModel.h>')
        self.writeln_hpp('namespace Nektar')

        self.writeln_hpp()
        
    def set_access(self, access):
        """Set the access specification for subsequent output.
        
        We keep track of the last access set, either via this method or
        output_method_start, and only output a new declaration to the
        header file if it changes.
        """
        if access != self._last_method_access:
            self._last_method_access = access
            self.writeln_hpp()
            # self.writeln_hpp(access, ':', indent_offset=-1)

    def output_method_start(self, method_name, args, ret_type, access=None, defaults=[]):
        """Output the start of a method declaration/definition.
        
        Will write to both the .hpp and .cpp file.
        
        We keep track of the access of the last method, and only output a new
        declaration to the header file if it changes.  The default is to use
        the same access specification as last time.
        """
        DEBUG('translator', 'Generating code for method', method_name)
        if access:
            self.set_access(access)
        if ret_type:
            if ret_type[-1] != ' ':
                ret_type = ret_type + ' '
        else:
            ret_type = ''
        args_string_cpp = ', '.join(filter(None, map(str, args)))
        if defaults:
            assert len(defaults) == len(args)
            args_with_default = []
            for (arg, default) in zip(map(str, args), map(str, defaults)):
                if arg:
                    if default:
                        args_with_default.append(arg + '=' + default)
                    else:
                        args_with_default.append(arg)
            args_string_hpp = ', '.join(args_with_default)
        else:
            args_string_hpp = args_string_cpp
        self.writeln_hpp(ret_type, method_name, '(', args_string_hpp, ')', self.STMT_END)
        self.writeln(ret_type, self.class_name, '::', method_name, '(', args_string_cpp, ')')

    def output_derived_quantities(self):
        """Output a ComputeDerivedQuantities method if any such quantities exist.
        
        Looks for variables annotated with pycml:derived-quantity=yes, and generates
        a method to compute all these variables from a given state.
        """
        dqs = self.derived_quantities
        if dqs:
            self.output_method_start('ComputeDerivedQuantities',
                                     [self.TYPE_DOUBLE + self.code_name(self.free_vars[0]),
                                      'const ' + self.TYPE_VECTOR + '& rY'], # We need it to really be a reference
                                     self.TYPE_VECTOR, access='public')
            self.open_block()
            self.output_comment('Inputs:')
            self.output_comment('Time units: ', self.free_vars[0].units)
            # Work out what equations are needed
            if self.use_chaste_stimulus:
                i_stim = [self.doc._cml_config.i_stim_var]            
            else:
                i_stim = []
            if self.use_data_clamp:
                prune = [self.config.i_data_clamp_data]
            else:
                prune = []
            nodeset = self.calculate_extended_dependencies(dqs, prune_deps=i_stim, prune=prune)


            # State variable inputs
            self.output_state_assignments(assign_rY=False, nodeset=nodeset)
            self.writeln()
            table_index_nodes_used = self.calculate_lookup_table_indices(nodeset, self.code_name(self.free_vars[0]))
            # Output equations
            self.output_comment('Mathematics')
            self.output_equations(nodeset - table_index_nodes_used)
            self.writeln()
            # Assign to results vector
            self.writeln(self.vector_create('dqs', len(dqs)))
            for i, var in enumerate(dqs):
                self.writeln(self.vector_index('dqs', i), self.EQ_ASSIGN, self.code_name(var), self.STMT_END)
            self.writeln('return dqs', self.STMT_END)
            self.close_block(blank_line=True)
    
    def output_serialize_method(self):
        """This method outputs the boost serialize method for the 
        header files that need it."""
        # Serialization
        if self.include_serialization:
            # self.writeln_hpp('friend class boost::serialization::access;')
            self.writeln_hpp('public:')
            # self.writeln_hpp('template<class Archive>')
            self.writeln_hpp('/// Creates an instance of this class', indent_offset=1)
            # self.writeln_hpp('void serialize(Archive & archive, const unsigned int version)')
            self.writeln_hpp('static CellModelSharedPtr create(', indent_offset=1)
            self.writeln_hpp('const LibUtilities::SessionReaderSharedPtr& pSession,', indent_offset=3)
            self.writeln_hpp('const MultiRegions::ExpListSharedPtr& pField)', indent_offset=3)    
            self.open_block(subsidiary=True,indent_offset=1)
            self.writeln_hpp('return MemoryManager<',self.class_name,'>::AllocateSharedPtr(pSession, pField);',indent_offset=1)
            # self.writeln_hpp('archive & boost::serialization::base_object<', self.base_class_name,
                             # ' >(*this);')
            if self.dynamically_loadable:
                self.writeln_hpp('archive & boost::serialization::base_object<AbstractDynamicallyLoadableEntity>(*this);')
            if self.use_modifiers:
                self.output_comment('Despite this class having modifier member variables, they are all added to the', subsidiary=True)
                self.output_comment('abstract class by the constructor, and archived via that, instead of here.', subsidiary=True) 
            self.close_block(subsidiary=True,indent_offset=1)

            self.writeln_hpp('/// Name of class',indent_offset=1)
            self.writeln_hpp('static std::string className;\n',indent_offset=1)


            class_name_length = len(self.class_name)+1
            self.writeln_hpp('/// Constructor',indent_offset=1)  
            self.writeln_hpp(self.class_name,'(const LibUtilities::SessionReaderSharedPtr& pSession,',indent_offset=1)    
            self.writeln_hpp(' '*class_name_length, 'const MultiRegions::ExpListSharedPtr& pField);\n',indent_offset=1) 

            self.writeln_hpp('/// Destructor',indent_offset=1)
            self.writeln_hpp('virtual ~',self.class_name,'() {}\n',indent_offset=1)

            self.writeln_hpp('protected:')
            self.writeln_hpp('/// Computes the reaction terms $f(u,v)$ and $g(u,v)$.',indent_offset=1)
            self.writeln_hpp('virtual void v_Update(',indent_offset=1)
            self.writeln_hpp('const Array<OneD, const  Array<OneD, NekDouble> >&inarray,',indent_offset=3)
            self.writeln_hpp('      Array<OneD,        Array<OneD, NekDouble> >&outarray,',indent_offset=3)    
            self.writeln_hpp('const NekDouble var_chaste_interface__environment__time);\n',indent_offset=3)

            self.writeln_hpp('/// Prints a summary of the model parameters.',indent_offset=1)
            self.writeln_hpp('virtual void v_GenerateSummary(SummaryList& s);\n',indent_offset=1)

            self.writeln_hpp('/// Set initial conditions for the cell model',indent_offset=1)
            self.writeln_hpp('virtual void v_SetInitialConditions();',indent_offset=1)
            self.writeln_hpp('};')
            self.close_block(subsidiary=True)
            self.writeln_hpp('#endif')

    def output_cell_parameters(self):
        """Output declarations, set & get methods for cell parameters.
        
        Sets self.cell_parameters to be those constant variables annotated with
        pycml:modifiable-parameter.  These use the mParameters functionality in
        Chaste.
        
        Also collects any variables annotated with an RDF oxmeta name into
        self.metadata_vars. Only constants and state variables are included.
        """
        # Find annotated parameters
        self.cell_parameters = filter(
            lambda v: v.is_modifiable_parameter,
            cellml_metadata.find_variables(self.model,
                                           ('pycml:modifiable-parameter', NSS['pycml']),
                                           'yes'))

        # Reduce intra-run variation
        self.cell_parameters.sort(key=self.var_display_name)

        for i, var in enumerate(self.cell_parameters):
            # Remember the var's index
            var._cml_param_index = i

        # Create set of all oxmeta-annotated variables
        vars = cellml_metadata.find_variables(self.model, ('bqbiol:is', NSS[u'bqbiol']))
        # Keep only the variables with an oxmeta name
        vars = filter(lambda v: v.oxmeta_name, vars)
        # We're interested in anything that isn't time or the stimulus
        self.metadata_vars = set([v for v in vars if v.get_type() != VarTypes.Free])
        self.metadata_vars.discard(self.doc._cml_config.i_stim_var)
        self.metadata_vars = list(self.metadata_vars)
        self.metadata_vars.sort(key=self.var_display_name)
        
        # #1464 Create a set of metadata variables that will have modifiers
        # We want to avoid writing out metadata for stimulus current as it is used once and then discarded.
        # \todo - use protocol information to put only the required modifiers into this list.
        self.modifier_vars = [v for v in self.metadata_vars if v.oxmeta_name not in cellml_metadata.STIMULUS_NAMES and v.oxmeta_name != 'membrane_capacitance']
        self.modifier_vars.sort(key=self.var_display_name)

        # Generate member variable declarations
        self.set_access('private')
        if self.metadata_vars:
            self.output_comment('\nSettable parameters and readable variables\n', subsidiary=True)
        
        # Write out the modifier member variables. 
        if self.use_modifiers:
            for var in self.modifier_vars:
                self.writeln_hpp('boost::shared_ptr<AbstractModifier> mp_' + var.oxmeta_name + '_modifier', self.STMT_END)    
        
        # Methods associated with oxmeta annotated variables
        # Don't use LT & modifiers for the const methods
        use_modifiers = self.use_modifiers
        self.use_modifiers = False
        use_lt = self.use_lookup_tables
        self.use_lookup_tables = False
        for var in self.metadata_vars:
            if var.is_statically_const(ignore_annotations=True):
#                 self.output_method_start('Get_' + var.oxmeta_name + '_constant', [], self.TYPE_DOUBLE)
#                 self.open_block()
#                 self.output_comment('Constant value given in CellML')
#                 nodeset = self.calculate_extended_dependencies([var])
#                 self.output_equations(nodeset)
#                 self.writeln('return ', self.code_name(var), self.STMT_END)
#                 self.close_block()
#                 self.writeln()
                if var in self.cell_parameters and var in self.modifier_vars:
                    # 'Forget' its index, so normal code generation occurs (#1647)
                    var._cml_has_modifier = True
        self.use_lookup_tables = use_lt
        self.use_modifiers = use_modifiers
        self.output_default_stimulus()
        self.output_intracellular_calcium()
        
        # Find & store derived quantities, for use elsewhere
        self.derived_quantities = cellml_metadata.find_variables(self.model,
                                                                 ('pycml:derived-quantity', NSS['pycml']),
                                                                 'yes')
        # Reduce intra-run variation
        self.derived_quantities.sort(key=self.var_display_name)
                
    def output_default_stimulus(self):
        """
        Output a default cell stimulus from the metadata specification
        as long as the following metadata exists:
         * membrane_stimulus_current_amplitude
         * membrane_stimulus_current_duration
         * membrane_stimulus_current_period
        and optionally:
         * membrane_stimulus_current_offset
         * membrane_stimulus_current_end
        
        Ensures that the amplitude of the generated RegularStimulus is negative.
        """
        vars = dict()
        for n in ['duration', 'amplitude', 'period', 'offset', 'end']:
            vars[n] = self.model.get_variable_by_oxmeta_name('membrane_stimulus_current_'+n, throw=False)
        if not (vars['duration'] and vars['amplitude'] and vars['period']):
            self.has_default_stimulus = False
            return
        self.has_default_stimulus = True
        nodeset = self.calculate_extended_dependencies(filter(None, vars.values()))

        self.output_method_start('UseCellMLDefaultStimulus', [], 'boost::shared_ptr<RegularStimulus>', 'public')
        self.open_block()
        self.output_comment('Use the default stimulus specified by CellML metadata')
        self.output_equations(nodeset)
        self.writeln('boost::shared_ptr<RegularStimulus> p_cellml_stim(new RegularStimulus(')
        self.writeln('        -fabs(', self.code_name(vars['amplitude']), '),')
        self.writeln('        ', self.code_name(vars['duration']), ',')
        self.writeln('        ', self.code_name(vars['period']), ',')
        if vars['offset']:
            self.writeln('        ', self.code_name(vars['offset']))
        else:
            self.writeln('        0.0')
        if vars['end']:
            self.writeln('      , ', self.code_name(vars['end']))
        self.writeln('        ))', self.STMT_END)
        self.writeln('mpIntracellularStimulus = p_cellml_stim', self.STMT_END)
        self.writeln('return p_cellml_stim', self.STMT_END)
        self.close_block(blank_line=True)
    
    def output_intracellular_calcium(self):
        """
        If a (state) variable has been annotated as cytosolic_calcium_concentration,
        generate a GetIntracellularCalciumConcentration method.
        """
        # Find cytosolic_calcium_concentration
        cai = self.doc.model.get_variable_by_oxmeta_name('cytosolic_calcium_concentration', throw=False)
        if cai and cai in self.state_vars:
            i = self.state_vars.index(cai[0])
            self.output_method_start('GetIntracellularCalciumConcentration', [], self.TYPE_DOUBLE, 'public')
            self.open_block()
            self.writeln('return ', self.vector_index('mStateVariables', i), self.STMT_END)
            self.close_block(blank_line=True)
        
    def code_name(self, var, *args, **kwargs):
        """
        Return the full name of var in a form suitable for inclusion in a source file.
        
        Overrides the base class version to access mParameters for parameters.
        """
        #self.writeln(var)
        if hasattr(var, '_cml_param_index') and not (self.use_modifiers and getattr(var, '_cml_has_modifier', False)):
            return self.vector_index('mParameters', var._cml_param_index)
        elif var is getattr(self.model, u'_cml_Chaste_Cm', None):
            return 'HeartConfig::Instance()->GetCapacitance()'
        elif hasattr(var, '_cml_code_name'):
            return var._cml_code_name % {'time': self.code_name(self.free_vars[0])}
        else:
            # self.writeln(var)
            # self.writeln(super(CellMLToNektarTranslator, self).code_name(var, *args, **kwargs))
            return super(CellMLToNektarTranslator, self).code_name(var, *args, **kwargs)


    def output_top_boilerplate(self):
        """Output top boilerplate.
        
        This method outputs the constructor and destructor of the cell
        class, and also lookup table declarations and lookup methods.
        It also calls output_verify_state_variables.
        """
        self.include_serialization = True

        # Check if we're generating a Backward Euler model
        self.use_backward_euler = self.model.get_option('backward_euler')
        self.use_analytic_jacobian = (self.model.get_option('maple_output') and hasattr(self.model.solver_info, u'jacobian'))
        if self.use_backward_euler:
            assert hasattr(self.model, u'solver_info')
            # Find the size of the nonlinear system
            num_linear_odes = len(self.model.solver_info.xml_xpath(u'solver:linear_odes/m:math/m:apply'))
            self.nonlinear_system_size = len(self.state_vars) - 1 - num_linear_odes
            nonlinear_entries = self.model.solver_info.xml_xpath(u'solver:jacobian/solver:entry/@var_j')
            self.nonlinear_system_vars = map(self.varobj, nonlinear_entries[:self.nonlinear_system_size])
        # Start output
        self.output_includes()
        
        if self.use_backward_euler or self.options.rush_larsen or self.options.grl1 or self.options.grl2:
            # Keep the same signature as forward cell models, but note that the solver isn't used
            solver1 = 'boost::shared_ptr<AbstractIvpOdeSolver> /* unused; should be empty */'
            solver2 = ''
            #solver1 = solver2 = ''
        else: #this currently outputs the boilerplate stuff
            solver1 = 'boost::shared_ptrALPHA<AbstractIvpOdeSolver> pSolverBRAVO'
            solver2 = '"' + self.class_name + '"'

        if self.use_lookup_tables and self.separate_lut_class:
            self.output_lut_class()

        # Cell model class
        self.open_block(subsidiary=True)
        self.writeln_hpp('class ', self.class_name, self.class_inheritance)
        self.writeln_hpp('{\n')
        
        # Put the boost serialize() method in if requested.
        self.output_serialize_method()

        # Parameter declarations, and set & get methods (#666)
        self.output_cell_parameters()
        # Constructor
        self.set_access('public')

        self.output_constructor([solver1, 'boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus'],
                                [solver2, self.class_name + '::create','"Description of the model?"'])
        # Destructor
        #self.output_method_start(''+self.class_name, [], '')
        self.writeln(self.class_name + '::' + self.class_name + '(')
        self.writeln('const LibUtilities::SessionReaderSharedPtr& pSession,' , indent_offset=3)
        self.writeln('const MultiRegions::ExpListSharedPtr& pField):', indent_offset=3)
        self.writeln('CellModel(pSession, pField)', indent_offset=2)


        self.open_block()

        #writing out all the state variables as a comment
        self.writeln('/*')
        self.writeln('State Variables:')
        for var in self.state_vars:
            self.writeln(var)
        self.writeln()
        self.writeln('Free Variables:')
        for var in self.free_vars:
            self.writeln(var)
        self.writeln()
        # for var in self.used_vars:
        #     self.writeln(var)
        self.writeln('*/\n')

        #outputting the number of variables in the state_vars list
        self.writeln('m_nq   = pField->GetNpoints();\n')
        self.writeln('m_nvar = ', len(self.state_vars), ';')

   



        #create a list of tau and infs, if present
        #setting up the arrays to contain the taus, infs, alphas and betas
        self.taus = []
        self.infs = []
        self.alphas = []
        self.betas = []

        # creating the nodeset that contains the mathematics
        state_vars = self.doc.model.find_state_vars()
        derivs = set(map(lambda v: (v, self.free_vars[0]), state_vars))
        extra_nodes=set()
        i_stim = []    

        nonv_nodeset = self.calculate_extended_dependencies(derivs|extra_nodes, prune_deps=i_stim)
        prune = nonv_nodeset
        dvdt = (self.v_variable, self.free_vars[0])
        v_nodeset = self.calculate_extended_dependencies([dvdt], prune=prune, prune_deps=i_stim)
        all_nodes = nonv_nodeset|v_nodeset
        extra_table_nodes=set()
        table_index_nodes_used = self.calculate_lookup_table_indices(all_nodes|extra_table_nodes, self.code_name(self.free_vars[0]))
        self.NODESET = nonv_nodeset - table_index_nodes_used

        #searching through the nodeset for all the infs, taus, alphas and betas and adding them to their arrays
        for expr in (e for e in self.model.get_assignments() if e in self.NODESET):
            # self.writeln('test')
            if isinstance(expr, cellml_variable):
                codename = str(self.code_name(expr))
                # self.writeln(self.code_name(expr))
                if 'gate' in codename:
                    if '__tau_'in codename:
                        self.taus.append(codename)
                    if '_inf' in codename:
                        self.infs.append(codename)
                    if '_alpha' in codename:
                        self.alphas.append(codename)
                    if '_beta' in codename:
                        self.betas.append(codename)

        # debugging help to print the arrays
        # self.writeln('Taus: ' + str(taus))
        # self.writeln('Infs: ' + str(infs))
        # self.writeln('Alphas: ' + str(self.alphas))
        # self.writeln('Betas: ' + str(self.betas) + '\n')
        
        #initialising the gate_vars and concentration_vars counters
        gate_vars = 0
        concentration_vars = 0
        self.state_var_type = {}

        #finding the actual variable name for each of the state variables
        for var in self.state_vars:
            # self.writeln()
            # self.writeln(var)
            var_actual = str(var)[33:str(var).rfind('__')]
            # self.writeln(var_actual)


            #writing the gate and concentration variables
            if str(var).find('membrane__V') != -1:
                self.state_var_type[str(var)] = 'voltage'
            elif filter(lambda element: var_actual in element,self.infs) and filter(lambda element: var_actual in element,self.taus) or filter(lambda element: var_actual in element,self.alphas) and filter(lambda element: var_actual in element,self.betas):
                gate_vars += 1
                # self.writeln(var, ' gating variable')
                self.writeln('m_gates.push_back(',str(gate_vars),');')
                self.state_var_type[str(var)] = 'm_gates'
                # self.writeln()
            else:
                gate_vars += 1
                # self.writeln(var, ' concentration variable')
                self.writeln('m_concentrations.push_back(',str(gate_vars),');')
                self.state_var_type[str(var)] = 'm_concentrations'
                # self.writeln()

        # self.writeln(self.state_var_type)

        self.close_block()






        # Other declarations & methods
        self.output_chaste_lut_methods()
        self.output_verify_state_variables()
        return
    
    @property
    def unsigned_v_index(self):
        if self.v_index == -1:
            return 'UNSIGNED_UNSET'
        else:
            return str(self.v_index)
    
    def output_verify_state_variables(self):
        """Output the VerifyStateVariables method.
        
        This will look for state variables annotated with pycml:range-low and/or pycml:range-high,
        which specify allowable ranges for these variables.  The generated method will check that
        they are within the range.  Both limits are included, i.e. they specify a closed interval.
        """
        
        # First work out if there are any constraints on state variables
    
        low_prop = ('pycml:range-low', NSS['pycml'])
        high_prop = ('pycml:range-high', NSS['pycml'])
        low_range_vars = filter(
            lambda v: v.get_type() == VarTypes.State,
            cellml_metadata.find_variables(self.model, low_prop))
        high_range_vars = filter(
            lambda v: v.get_type() == VarTypes.State,
            cellml_metadata.find_variables(self.model, high_prop))
        nodeset = set(low_range_vars + high_range_vars)
        
        # If not, don't bother writing the method, an empty implementation is in the abstract classes.
        if nodeset:
            # It's not appropriate to apply modifiers here - we want to check the actual values of the state
            use_modifiers = self.use_modifiers
            self.use_modifiers = False
        
            self.output_method_start('VerifyStateVariables', [], 'void')
            self.open_block()
            
            using_cvode = (self.TYPE_VECTOR_REF == CellMLToCvodeTranslator.TYPE_VECTOR_REF)
            if using_cvode:
                self.writeln('/* We only expect CVODE to keep state variables to within its tolerances,')
                self.writeln(' * not exactly the bounds prescribed to each variable that are checked here.')
                self.writeln(' *')
                self.writeln(' * For 99.99% of paces this->mAbsTol works,')
                self.writeln(' * For 99.999% of paces 10*this->mAbsTol is fine,')
                self.writeln(' * but unfortunately 100x seems to be required on rare occasions for upstrokes.')
                self.writeln(' * This sounds bad, but is probably typically only 1e-5 or 1e-6.')
                self.writeln(' */')
                self.writeln('const double tol = 100*this->mAbsTol;')
                            
            self.output_state_assignments(nodeset=nodeset)
            error_template = 'EXCEPTION(DumpState("State variable {0} has gone out of range. Check numerical parameters, for example time and space stepsizes, and/or solver tolerances"));'
            additional_tolerance_adjustment = ''
            for var in low_range_vars:
                if using_cvode:
                    additional_tolerance_adjustment = ' - tol'
                self.writeln('if (', self.code_name(var), ' < ', var.get_rdf_annotation(low_prop), additional_tolerance_adjustment, ')')
                self.open_block()
                #self.writeln('std::cout << "Too small: ', self.code_name(var), ' = " << ', self.code_name(var) , ' << std::endl << std::flush;')
                self.writeln(error_template.format(self.var_display_name(var)))
                self.close_block(False)
            for var in high_range_vars:
                if using_cvode:
                    additional_tolerance_adjustment = ' + tol'            
                self.writeln('if (', self.code_name(var), ' > ', var.get_rdf_annotation(high_prop), additional_tolerance_adjustment, ')')
                self.open_block()
                #self.writeln('std::cout << "Too large: ', self.code_name(var), ' = " << ', self.code_name(var) , ' << std::endl << std::flush;')
                self.writeln(error_template.format(self.var_display_name(var)))
                self.close_block(False)
            self.close_block(True)
            
            self.use_modifiers = use_modifiers
            

    def output_constructor(self, params, base_class_params):
        """Output a cell constructor.
        
        params is a list of constructor parameters, entries of which should be strings
        including both type and parameter name, which will be included verbatim in the
        generated code.
        
        base_class_params is a list of parameters to be supplied to the base class
        constructor.  Entries will be converted to strings.
        """
        #self.output_method_start(self.class_name, params, '', access='public')
        self.writeln('std::string ', self.class_name, '::className')
        self.writeln('= GetCellModelFactory().RegisterCreatorFunction(', indent_offset=2)
        #self.writeln('    : ', self.base_class_name, '(')

        # Filter out empty params, to make backward Euler happy
        base_class_params = filter(None, map(str, base_class_params))
        for i, param in enumerate(base_class_params):
            if i == len(base_class_params)-1: comma = ');'
            else: comma = ','
            self.writeln(param, comma, indent_offset=3)
       
        #self.open_block()
        self.writeln()
        self.writeln('/**')
        self.writeln('*')
        self.writeln('*/')
        # self.writeln('this->mpSystemInfo = OdeSystemInformation<',
        #             self.class_name, '>::Instance();')
        # if self.v_index == -1 and self.v_variable:
        #    self.writeln('this->mVoltageIndex = GetAnyVariableIndex("',
        #                 self.var_display_name(self.v_variable), '");')
        # if self.config.options.include_dt_in_tables:
        #    self.writeln(self.lt_class_name, '::Instance()->SetTimestep(mDt);')
        # self.writeln('Init();\n')
        
        #1861 - Rush-Larsen
        if self.options.rush_larsen and not self.doc._cml_rush_larsen:
            self.writeln('WARNING("No elligible gating variables found for this Rush-Larsen cell model; using normal forward Euler.");')
        
        #1463 - default cellML stimulus
        if self.has_default_stimulus:
            self.output_comment('We have a default stimulus specified in the CellML file metadata')
            self.writeln('this->mHasDefaultStimulusFromCellML = true', self.STMT_END)
            
        #1464 - cleverer modifiers...
        if self.use_modifiers and self.modifier_vars:
            self.output_comment('These will get initialised to DummyModifiers in the base class method.')
            for var in self.modifier_vars:
                self.writeln('this->AddModifier("' + var.oxmeta_name + '",')
                self.writeln('                  mp_' + var.oxmeta_name + '_modifier)', self.STMT_END)        
        
        #666 - initialise parameters
        for var in self.cell_parameters:
            if var.get_type() == VarTypes.Constant:
                self.writeln(self.vector_index('this->mParameters', var._cml_param_index),
                             self.EQ_ASSIGN, var.initial_value, self.STMT_END, ' ',
                             self.COMMENT_START, var.fullname(), ' [', var.units, ']')
        #1354 - specify protocol outputs
        
        if self.use_protocol:
            outputs = cellml_metadata.find_variables(self.model,
                                                     ('pycml:output-variable', NSS['pycml']),
                                                     'yes')
            def write_output_info(output):
                if output.get_type() in [VarTypes.Free, VarTypes.Unknown]:
                    self.writeln('UNSIGNED_UNSET, FREE', indent=False, nl=False)
                elif output.get_type() == VarTypes.State:
                    self.writeln(self.state_vars.index(output), ', STATE', indent=False, nl=False)
                elif output.is_derived_quantity:
                    self.writeln(self.derived_quantities.index(output), ', DERIVED', indent=False, nl=False)
                elif output.is_modifiable_parameter:
                    self.writeln(self.cell_parameters.index(output), ', PARAMETER', indent=False, nl=False)
                else:
                    raise ValueError('Unexpected protocol output: ' + str(output))
            if outputs:
                outputs.sort(key=lambda v: self.var_display_name(v))
                self.output_comment('Protocol outputs')
                self.writeln('this->mOutputsInfo.resize(', len(outputs), ');')
                for i, output in enumerate(outputs):
                    self.writeln('this->mOutputsInfo[', i, ']', self.EQ_ASSIGN,
                                 'std::make_pair(', nl=False)
                    write_output_info(output)
                    self.writeln(')', self.STMT_END, indent=False)
                self.writeln()
            outputs = set(outputs)
            #1925 - outputs that are vectors
            prop = ('pycml:output-vector', NSS['pycml'])
            vector_names = set(cellml_metadata.get_targets(self.model, None,
                                                           cellml_metadata.create_rdf_node(prop)))
            self.writeln('this->mVectorOutputsInfo.resize(', len(vector_names), ');')
            self.writeln('this->mVectorOutputNames.resize(', len(vector_names), ');')
            for i, name in enumerate(sorted(vector_names)):
                self.writeln('this->mVectorOutputNames[', i, ']', self.EQ_ASSIGN, '"', name, '"', self.STMT_END)
                vector_outputs = cellml_metadata.find_variables(self.model, prop, name)
                assert len(vector_outputs) > 0
                vector_outputs.sort(key=lambda v: self.var_display_name(v))
                self.writeln('this->mVectorOutputsInfo[', i, '].resize(', len(vector_outputs), ');')
                for j, output in enumerate(vector_outputs):
                    self.writeln('this->mVectorOutputsInfo[', i, '][', j, ']', self.EQ_ASSIGN,
                                 'std::make_pair(', nl=False)
                    write_output_info(output)
                    self.writeln(')', self.STMT_END, indent=False)
                self.writeln()
                outputs.update(vector_outputs)
            #1910 - SED-ML name mappings
            prop = ('pycml:alias', NSS['pycml'])
            aliased_vars = cellml_metadata.find_variables(self.model, prop, None)
            prop = cellml_metadata.create_rdf_node(prop)
            for var in aliased_vars:
                assert var in outputs
                source = cellml_metadata.create_rdf_node(fragment_id=var.cmeta_id)
                for alias in cellml_metadata.get_targets(self.model, source, prop):
                    name = self.var_display_name(var)
                    self.writeln('this->mNameMap["', alias, '"] = "', name, '";')
            #2178 - set up model outputs environment from above info
            self.writeln()
            self.writeln('ProcessOutputsInfo();')
            self.writeln()
            #2428 - also record protocol inputs
            inputs = cellml_metadata.find_variables(self.model, ('pycml:input-variable', NSS['pycml']), 'yes')
            if inputs:
                inputs.sort(key=lambda v: self.var_display_name(v))
                self.writeln('this->mInputNames.reserve(', len(inputs), ');')
                for input in inputs:
                    self.writeln('this->mInputNames.push_back("', self.var_display_name(input), '");')

        # Lookup table generation, if not in a singleton
        if self.use_lookup_tables and not self.separate_lut_class:
            self.output_lut_generation()
        self.output_extra_constructor_content()
        #self.close_block()

        return
    
    def output_extra_constructor_content(self):
        """Hook for subclasses to add further content to the constructor."""
        pass
    
    def output_chaste_lut_methods(self):
        """
        Output lookup table declarations & methods, if not using a separate class,
        or output the method to get a pointer to the lookup table collection.
        """
        if self.use_lookup_tables:
            if self.separate_lut_class:
                self.output_method_start('GetLookupTableCollection', [], 'AbstractLookupTableCollection*')
                self.open_block()
                self.writeln('return ', self.lt_class_name, '::Instance();')
                self.close_block()
            else:
                self.send_main_output_to_subsidiary()
                self.output_lut_declarations()
                self.output_lut_row_lookup_memory()
                self.output_lut_methods()
                self.send_main_output_to_subsidiary(False)
    
    def lut_parameters(self, key):
        """Get the bounds and step size for a particular table.
        
        key should be a key into self.lookup_table_indices.
        Returns (min, max, step) suitable for putting in generated code.
        """
        if self.separate_lut_class:
            idx = self.doc.lookup_table_indexes[key]
            return map(lambda s: 'mTable%s[%s]' % (s, idx), ['Mins', 'Maxs', 'Steps', 'StepInverses'])
        else:
            return super(CellMLToNektarTranslator, self).lut_parameters(key)
    
    def output_lut_indexing_methods(self):
        """Output methods in the LT class for indexing the tables, and checking index bounds.
        
        These will be methods like
            const double * const IndexTable0(double index_var);
        if self.row_lookup_method, or like
            void IndexTable0(double index_var, unsigned& index, double& factor);
        otherwise, with
            bool CheckIndex0(double& index_var);
        for checking the bounds.
        """
        for key, idx in self.doc.lookup_table_indexes.iteritems():
            varname = self.code_name(key[-1])
            method_name = 'IndexTable' + str(idx)
            if self.row_lookup_method:
                method = 'const double * %s(double %s)' % (method_name, varname)
            else:
                factor = self.lut_factor(idx)
                idx_var = '_table_index_' + str(idx)
                if factor:
                    factor = ', double& ' + factor
                method = 'void %s(double %s, unsigned& %s%s)' % (method_name, varname, idx_var, factor)
            self.writeln(method)
            self.open_block()
            self.output_table_index_generation_code(key, idx, call_method=False)
            if self.row_lookup_method:
                self.writeln('return _lt_', idx, '_row;')
            self.close_block()
            # And check the indexes
            if self.config.options.check_lt_bounds:
                self.writeln('#define COVERAGE_IGNORE', indent=False)
                self.writeln('bool CheckIndex', idx, '(double& ', varname, ')')
                self.open_block()
                self.output_table_index_checking(key, idx, call_method=False)
                self.writeln('return _oob_', idx, self.STMT_END)
                self.close_block(blank_line=False)
                self.writeln('#undef COVERAGE_IGNORE\n', indent=False)
    
    def output_table_index_checking(self, key, idx, call_method=True):
        """Override base class method to call the methods on the lookup table class if needed."""
        if self.separate_lut_class and call_method:
            if self.config.options.check_lt_bounds:
                var = key[-1]
                varname = self.code_name(var)
                self.writeln('const bool _oob_', idx, self.EQ_ASSIGN, self.lt_class_name,
                             '::Instance()->CheckIndex', idx, '(', varname, ')', self.STMT_END)
        else:
            super(CellMLToNektarTranslator, self).output_table_index_checking(key, idx)
    
    def output_table_index_generation_code(self, key, idx, call_method=True):
        """Override base class method to call the methods on the lookup table class if needed."""
        if self.separate_lut_class and call_method:
            var = key[-1]
            varname = self.code_name(var)
            method_name = self.lt_class_name + '::Instance()->IndexTable' + str(idx)
            if self.row_lookup_method:
                self.writeln('const double* const _lt_', idx, '_row = ', method_name, '(', varname, ');')
            else:
                factor = self.lut_factor(idx, include_comma=True)
                idx_var = '_table_index_' + str(idx)
                self.writeln(method_name, '(', varname, ', ', idx_var, factor, ');')
        else:
            super(CellMLToNektarTranslator, self).output_table_index_generation_code(key, idx)

    def output_lut_class(self):
        """Output a separate class for lookup tables.
        
        This will live entirely in the .cpp file."""
        # Lookup tables class
        self.writeln('class ', self.lt_class_name, ' : public AbstractLookupTableCollection')
        self.writeln('{')
        self.writeln('public:')
        self.set_indent(1)
        # Method to get the table instance object
        self.writeln('static ', self.lt_class_name, '* Instance()')
        self.open_block()
        self.writeln('if (mpInstance.get() == NULL)')
        self.writeln('{')
        self.writeln('mpInstance.reset(new ', self.lt_class_name, ');', indent_offset=1)
        self.writeln('}')
        self.writeln('return mpInstance.get();')
        self.close_block()
        # Method to free the table memory
        self.writeln('void FreeMemory()')
        self.open_block()
        self.output_lut_deletion()
        self.writeln('mNeedsRegeneration.assign(mNeedsRegeneration.size(), true);')
        self.close_block()
        # Table lookup methods
        self.output_lut_methods()
        self.output_lut_indexing_methods()
        # Destructor
        self.writeln('~', self.lt_class_name, '()')
        self.open_block()
        self.output_lut_deletion()
        self.close_block()
        # Make the class a singleton
        self.writeln('protected:', indent_level=0)
        self.writeln(self.lt_class_name, '(const ', self.lt_class_name, '&);')
        self.writeln(self.lt_class_name, '& operator= (const ', self.lt_class_name, '&);')
        # Constructor
        self.writeln(self.lt_class_name, '()')
        self.open_block()
        self.writeln('assert(mpInstance.get() == NULL);')
        if self.config.options.include_dt_in_tables:
            self.writeln('mDt = HeartConfig::Instance()->GetOdeTimeStep();')
            self.writeln('assert(mDt > 0.0);')
        num_indexes = len(self.doc.lookup_table_indexes)
        self.writeln('mKeyingVariableNames.resize(', num_indexes, ');')
        self.writeln('mNumberOfTables.resize(', num_indexes, ');')
        self.writeln('mTableMins.resize(', num_indexes, ');')
        self.writeln('mTableSteps.resize(', num_indexes, ');')
        self.writeln('mTableStepInverses.resize(', num_indexes, ');')
        self.writeln('mTableMaxs.resize(', num_indexes, ');')
        self.writeln('mNeedsRegeneration.resize(', num_indexes, ');')
        for key, idx in self.doc.lookup_table_indexes.iteritems():
            min, max, step, var = key
            num_tables = unicode(self.doc.lookup_tables_num_per_index[idx])
            self.writeln('mKeyingVariableNames[', idx, '] = "', self.var_display_name(var), '";')
            self.writeln('mNumberOfTables[', idx, '] = ', num_tables, self.STMT_END)
            self.writeln('mTableMins[', idx, '] = ', min, self.STMT_END)
            self.writeln('mTableSteps[', idx, '] = ', step, self.STMT_END)
            self.writeln('mTableStepInverses[', idx, '] = ', str(1/float(step)), self.STMT_END)
            self.writeln('mTableMaxs[', idx, '] = ', max, self.STMT_END)
            self.writeln('mNeedsRegeneration[', idx, '] = true;')
            self.writeln('_lookup_table_', idx, self.EQ_ASSIGN, 'NULL', self.STMT_END)
        self.writeln(self.lt_class_name, '::RegenerateTables();')
        self.close_block()
        # Table (re-)generation
        self.writeln('void RegenerateTables()')
        self.open_block()
        event_handler = 'AbstractLookupTableCollection::EventHandler::'
        self.writeln(event_handler, 'BeginEvent(', event_handler, 'GENERATE_TABLES);')
        if self.config.options.include_dt_in_tables:
            self.writeln(self.TYPE_CONST_DOUBLE, self.code_name(self.config.dt_variable), ' = mDt;')
            # Hack: avoid unused variable warning
            self.writeln('double _unused = ', self.code_name(self.config.dt_variable), ';')
            self.writeln('_unused = _unused;\n')
        for idx in self.doc.lookup_table_indexes.itervalues():
            self.writeln('if (mNeedsRegeneration[', idx, '])')
            self.open_block()
            self.output_lut_deletion(only_index=idx)
            self.output_lut_generation(only_index=idx)
            self.writeln('mNeedsRegeneration[', idx, '] = false;')
            self.close_block(blank_line=True)
        self.writeln(event_handler, 'EndEvent(', event_handler, 'GENERATE_TABLES);')
        self.close_block()
        # Private data
        self.writeln('private:', indent_level=0)
        self.writeln('/** The single instance of the class */')
        self.writeln('static std::auto_ptr<', self.lt_class_name, '> mpInstance;\n')
        if self.row_lookup_method:
            self.output_lut_row_lookup_memory()
        self.output_lut_declarations()
        # Close the class
        self.set_indent(0)
        self.writeln('};\n')
        # Define the instance pointer
        self.writeln('std::auto_ptr<', self.lt_class_name, '> ', self.lt_class_name, '::mpInstance;')
        self.writeln()
        return

    def output_state_assignments(self, exclude_nonlinear=False,
                                 assign_rY=True,
                                 nodeset=None,
                                 pointer=''):
        """Output statements extracting state variables from their vector.

        If exclude_nonlinear is set to true, state variables appearing
        in the nonlinear system will not be included.

        If nodeset is given, only state variables appearing in nodeset
        will be included.
        
        If pointer is given, then the state variables actually appear in the
        variable given by pointer, which is of type const std::vector<double>*.
        """
        used_vars = set()
        for var in self.state_vars:
            if ((not exclude_nonlinear or var not in self.nonlinear_system_vars)
                and (nodeset is None or var in nodeset)):
                used_vars.add(var)
        if assign_rY and used_vars:
            if pointer:
                self.output_comment('For state variable interpolation (SVI) we read in interpolated state variables,')
                self.output_comment('otherwise for ionic current interpolation (ICI) we use the state variables of this model (node).')
                if self.TYPE_VECTOR_REF == CellMLToNektarTranslator.TYPE_VECTOR_REF:
                    self.writeln('if (!%s) %s = &rGetStateVariables();' % (pointer, pointer))
                    self.writeln('const ', self.TYPE_VECTOR_REF, 'rY = *', pointer, self.STMT_END)
                else:
                    self.writeln(self.TYPE_VECTOR_REF, 'rY;')
                    self.writeln('bool made_new_cvode_vector = false;')
                    self.writeln('if (!%s)' % (pointer))
                    self.open_block()
                    self.writeln('rY = rGetStateVariables();')
                    self.close_block(False)
                    self.writeln('else')
                    self.open_block()
                    self.writeln('made_new_cvode_vector = true;')
                    self.writeln('rY = MakeNVector(*%s);' % (pointer))
                    self.close_block()
            else:
                self.writeln(self.TYPE_VECTOR_REF, 'rY = rGetStateVariables();')
        if self.options.protocol:
            low_prop = ('pycml:range-low', NSS['pycml'])
            high_prop = ('pycml:range-high', NSS['pycml'])
            def check_bound(prop, reln, var, value):
                prop_value = var.get_rdf_annotation(prop)
                if prop_value:
                    value = '(%s %s %s ? %s : %s)' % (value, reln, prop_value, prop_value, value)
                return value
        for i, var in enumerate(self.state_vars):
            if var in used_vars:
                if self.use_modifiers and var in self.modifier_vars:
                    value = self.modifier_call(var, self.vector_index('rY', i))
                else:
                    value = self.vector_index('inarray', i) + '[i]'
                if self.options.protocol:
                    value = check_bound(low_prop, '<', var, value)
                    value = check_bound(high_prop, '>', var, value)
                #2116 - use supplied fixed voltage if we're clamping
                # if var is self.v_variable:
                #     value = '(mSetVoltageDerivativeToZero ? this->mFixedVoltage : %s)' % value
                self.writeln(self.TYPE_DOUBLE, self.code_name(var),
                             self.EQ_ASSIGN, value, self.STMT_END)
                self.writeln(self.COMMENT_START, 'Units: ', var.units,
                             '; Initial value: ',
                             getattr(var, u'initial_value', 'Unknown'))
                #621 TODO: maybe convert if state var dimensions include time
        self.writeln()
        return
    
    def modifier_call(self, var, current_value):
        """Return code for a call to a modifier function for an oxmeta-annotated variable.
        
        The modifier function takes 2 parameters: the current value of the variable,
        and the current time.  It returns a modified value for the variable.
        """
        return ('mp_' + var.oxmeta_name + '_modifier->Calc(' +
                current_value + ', ' + self.code_name(self.free_vars[0]) + ')')
    
    def vector_index(self, vector, i):
        """Return code for accessing the i'th index of vector."""
        return vector + '[' + str(i) + ']'
    
    def vector_create(self, vector, size):
        """Return code for creating a new vector with the given size."""
        return ''.join(map(str, [self.TYPE_VECTOR, vector, '(', size, ')', self.STMT_END]))
    
    def vector_initialise(self, vector, size):
        """Return code for creating an already-declared vector with the given size."""
        return ''.join(map(str, [vector, '.resize(', size, ')', self.STMT_END]))
    
    def output_nonlinear_state_assignments(self, nodeset=None):
        """Output assignments for nonlinear state variables."""
        for i, var in enumerate(self.nonlinear_system_vars):
            if not nodeset or var in nodeset:
                self.writeln(self.TYPE_DOUBLE, self.code_name(var), self.EQ_ASSIGN,
                             self.vector_index('rCurrentGuess', i), self.STMT_END)
                #621 TODO: maybe convert if state var dimensions include time
        self.writeln()
        return
    
    def get_stimulus_assignment(self):
        """Return code for getting Chaste's stimulus current."""
        expr = self.doc._cml_config.i_stim_var
        output = self.code_name(expr) + self.EQ_ASSIGN
        get_stim = 'GetIntracellularAreaStimulus(' + self.code_name(self.free_vars[0]) + ')'
        if self.doc._cml_config.i_stim_negated:
            get_stim = '-' + get_stim
        return output + get_stim + self.STMT_END

    def output_equations(self, nodeset, zero_stimulus=False):
        """Output the mathematics described by nodeset.

        nodeset represents a subset of the assignments in the model.
        Output assignments in the order given by a topological sort,
        but only include those in nodeset.
        """
        # Special case for the stimulus current
        

        if self.doc._cml_config.i_stim_var in nodeset:
            if zero_stimulus:
                i_stim = self.doc._cml_config.i_stim_var
                stim_assignment = self.code_name(i_stim) + self.EQ_ASSIGN + '0.0' + self.STMT_END
            else:
                stim_assignment = self.get_stimulus_assignment()            

        for expr in (e for e in self.model.get_assignments() if e in nodeset):
            
            # Special-case the stimulus current
            # self.writeln(expr)
            # if self.use_chaste_stimulus or zero_stimulus:
            #     if isinstance(expr, cellml_variable) and expr is self.doc._cml_config.i_stim_var:
            #         self.writeln(self.TYPE_CONST_DOUBLE, stim_assignment)
            #     elif not (isinstance(expr, mathml_apply) and
            #               isinstance(expr.operator(), mathml_eq) and
            #               isinstance(expr.eq.lhs, mathml_ci) and
            #               expr.eq.lhs.variable is self.doc._cml_config.i_stim_var):
            #         self.output_assignment(expr)
            
            
            self.output_assignment(expr)
        return

    def output_assignment(self, expr):
        """Output an assignment statement.

        Has overrides for various special cases.
        """
        clear_type = False
        writing_data_clamp_current = False
        # Figure out what is being assigned to
        if isinstance(expr, cellml_variable):
            assigned_var = expr
        else:
            if expr.eq.lhs.localName == 'ci':
                assigned_var = expr.eq.lhs.variable
                if assigned_var is self.config.i_data_clamp_current:
                    writing_data_clamp_current = True
                    self.output_comment('Special handling of data clamp current here (see #2708)')
                    self.output_comment('(we want to save expense of calling the interpolation method if possible.)')
                    self.writeln(self.TYPE_DOUBLE, self.code_name(assigned_var), self.EQ_ASSIGN, '0.0' , self.STMT_END)
                    self.writeln('if (mDataClampIsOn)')
                    self.open_block()
                    clear_type = True
            else:
                assigned_var = None # We don't store derivatives as members
                #907: Check if this is the derivative of the transmembrane potential
                if not self.use_backward_euler and expr.eq.lhs.diff.dependent_variable == self.v_variable:
                    clear_type = True
                    
        # Parameters don't need assigning
        has_modifier = self.use_modifiers and getattr(assigned_var, '_cml_has_modifier', False)
        # self.writeln(has_modifier)
        if assigned_var in self.cell_parameters and not has_modifier:
            return
        # Is the variable declared elsewhere?
        if clear_type:
            self.TYPE_DOUBLE = self.TYPE_CONST_DOUBLE = ''
        elif getattr(assigned_var, '_cml_modifiable', False):
            # Override const-ness, e.g. for a lookup table index
            self.TYPE_CONST_DOUBLE = self.TYPE_DOUBLE
        if (assigned_var and self.use_modifiers and assigned_var in self.modifier_vars
            and assigned_var.get_type() != VarTypes.State):
            # "Constant" oxmeta-annotated parameters may be modified at run-time
            if has_modifier:
                # Turn off the modifier to figure out the base value
                assigned_var._cml_has_modifier = False
                rhs = self.code_name(assigned_var)
                assigned_var._cml_has_modifier = True
            else:
                self.capture_output()
                super(CellMLToNektarTranslator, self).output_assignment(expr)
                assignment = self.get_captured_output()
                eq_pos = assignment.find(self.EQ_ASSIGN)
                end_pos = assignment.find(self.STMT_END)
                rhs = assignment[eq_pos+len(self.EQ_ASSIGN):end_pos]
            if rhs:
                # If assigned_var is computed, it'll 'appear' twice - once with expr==assigned_var,
                # and once for the assignment mathml_apply.  The former will result in an empty rhs.
                # self.writeln('Test')
                self.writeln(self.TYPE_CONST_DOUBLE, self.code_name(assigned_var), self.EQ_ASSIGN,
                             self.modifier_call(assigned_var, rhs), self.STMT_END, nl=False)

                self.output_comment(assigned_var.units, indent=False, pad=True)
        else:
            # self.writeln('Test')
            super(CellMLToNektarTranslator, self).output_assignment(expr)
#        if assigned_var:
#            # Debug
#            self.writeln('EXCEPT_IF_NOT(!std::isinf(', self.code_name(assigned_var), '));')
#            self.writeln('EXCEPT_IF_NOT(!std::isnan(', self.code_name(assigned_var), '));')
        if clear_type:
            # Remove the instance attributes, thus reverting to the class members
            del self.TYPE_DOUBLE
            del self.TYPE_CONST_DOUBLE
        elif getattr(assigned_var, '_cml_modifiable', False):
            del self.TYPE_CONST_DOUBLE
            
        if writing_data_clamp_current:
            self.close_block(False)
            
        return

    def output_mathematics(self):
        """Output the mathematics in this model.

        When backward Euler is used, we do so in 5 methods:
         * UpdateTransmembranePotential  does a forward Euler step for V
         * ComputeOneStepExceptVoltage  co-ordinates a backward Euler step
         * ComputeResidual and ComputeJacobian are used in the Newton iteration
         * GetIIonic returns the total ionic current
        
        Rush-Larsen is implemented similarly, with:
         * EvaluateEquations  evaluate the model derivatives and alpha/beta terms
         * ComputeOneStepExceptVoltage  does a Rush-Larsen update for eligible variables,
           and a forward Euler step for other non-V state variables
        Generalised Rush-Larsen methods also have specialised handling; see the
        individual methods for details.

        For other solvers, only 2 methods are needed:
         * EvaluateYDerivatives computes the RHS of the ODE system
         * GetIIonic is as above
        
        Where derived-quantity annotations are present, we also generate a
        ComputeDerivedQuantities method.
        """
        #self.output_get_i_ionic()
        if self.options.rush_larsen:
            self.output_rush_larsen_mathematics()
        elif self.use_backward_euler:
            self.output_backward_euler_mathematics()
        elif self.options.grl1:
            self.output_grl1_mathematics()
        elif self.options.grl2:
            self.output_grl2_mathematics()
        else:
            self.output_evaluate_y_derivatives()
        self.output_derived_quantities()
    
    def calculate_lookup_table_indices(self, nodeset, time_name=None):
        """Output the lookup table index calculations needed for the given equations, if tables are enabled.
        
        If time_name is given, it may be used in exception messages for tables out of bounds.
        Note that it is needed to be passed in, since GetIIonic does not know the time.
        
        Returns the subset of nodeset used in calculating the indices.
        """
        if self.use_lookup_tables:
            nodes_used = self.output_table_index_generation(time_name, nodeset=nodeset)
        else:
            nodes_used = set()
        return nodes_used

    def output_get_i_ionic(self):
        """Output the GetIIonic method."""
        use_modifiers = self.use_modifiers
        self.use_modifiers = False
        self.output_method_start('GetIIonic', ['const std::vector<double>* pStateVariables'],
                                 self.TYPE_DOUBLE, access='public', defaults=['NULL'])
        self.open_block()
        # Output mathematics to calculate ionic current, using solver_info.ionic_current.
        if (hasattr(self.model, u'solver_info') and hasattr(self.model.solver_info, u'ionic_current')):
            if not hasattr(self.model.solver_info.ionic_current, u'var'):
                raise ValueError('No ionic currents found; check your configuration file')
            nodes = map(lambda elt: self.varobj(unicode(elt)),
                        self.model.solver_info.ionic_current.var)
            # GetIIonic must not include the stimulus current
            i_stim = self.doc._cml_config.i_stim_var
            nodeset = self.calculate_extended_dependencies(nodes, prune_deps=[i_stim])
            #print map(lambda v: v.fullname(), nodes)
            #print filter(lambda p: p[2]>0, map(debugexpr, nodeset))
            # Output main part of maths
            self.output_state_assignments(nodeset=nodeset, pointer='pStateVariables')
            table_index_nodes_used = self.calculate_lookup_table_indices(nodeset)
            self.output_equations(nodeset - table_index_nodes_used, zero_stimulus=True)
            self.writeln()
            # Assign the total current to a temporary so we can check for NaN
            self.writeln(self.TYPE_CONST_DOUBLE, 'i_ionic', self.EQ_ASSIGN, nl=False)
            if self.doc._cml_config.i_ionic_negated:
                self.writeln('-(', nl=False, indent=False)
            plus = False
            for varelt in self.model.solver_info.ionic_current.var:
                if plus: self.write('+')
                else: plus = True
                self.output_variable(varelt)
            if self.doc._cml_config.i_ionic_negated:
                self.writeln(')', nl=False, indent=False)
            self.writeln(self.STMT_END, indent=False)
            """if self.TYPE_VECTOR_REF == CellMLToCvodeTranslator.TYPE_VECTOR_REF:
                self.writeln('if (made_new_cvode_vector)')
                self.open_block()
                self.writeln('DeleteVector(rY);')
                self.close_block(False)"""
            self.writeln('EXCEPT_IF_NOT(!std::isnan(i_ionic));')
            self.writeln('return i_ionic', self.STMT_END)
        else:
            self.writeln('return 0.0;')
        self.close_block()
        self.use_modifiers = use_modifiers

    def output_evaluate_y_derivatives(self, method_name='EvaluateYDerivatives'):
        """Output the EvaluateYDerivatives method."""
        # Start code output

        #write the method start
        self.writeln(self.TYPE_VOID, self.class_name, '::v_Update(')
        self.writeln('const Array<OneD, const  Array<OneD, NekDouble> >&inarray,', indent_offset=4)
        self.writeln('  Array<OneD,        Array<OneD, NekDouble> >&outarray,', indent_offset=5)
        self.writeln(' const NekDouble var_chaste_interface__environment__time)', indent_offset=13)
        self.open_block()
        self.writeln('for (unsigned int i = 0; i < m_nq; ++i)')

        self.open_block()
        if not self.state_vars:
            # This isn't an ODE model!
            self.close_block()
            return
        self.output_comment('Inputs:')
        self.output_comment('Time units: ', self.free_vars[0].units)
        self.output_derivative_calculations(self.state_vars)

        # debugging help by writing out the arrays of taus, infs, alphas, betas and dictionary of variable types
        # self.writeln('Taus: ' + str(self.taus))
        # self.writeln('Infs: ' + str(self.infs))
        # self.writeln('Alphas: ' + str(self.alphas))
        # self.writeln('Betas: ' + str(self.betas))
        # self.writeln('Var Types: ' + str(self.state_var_type)+ '\n')
       
        #write the taus
        self.output_comment('Calculating the tau-values:')
        for i,var in enumerate(self.state_vars):
            var_name = self.code_name(var, True)
            var_actual = str(var)[33:str(var).rfind('__')]
            if 'gate' in str(var_name):
                after_underscore = var_name[var_name.rfind('__')+2:]
        

            if 'membrane__V' in var:
                continue
            
            if self.state_var_type[str(var)] == 'm_gates':                  
                if filter(lambda element: var_actual in element,self.taus):
                    self.output_comment('The tau value for ' + str(var_actual) + ' should already be calculated in the mathematics')
                else:
                    if filter(lambda element: var_actual in element,self.alphas) and filter(lambda element: var_actual in element,self.betas):
                        alpha = filter(lambda element: var_actual in element,self.alphas)[0]
                        beta = filter(lambda element: var_actual in element,self.betas)[0]
                        new_tau = 'var_' + str(var_actual) + '__tau_' + str(after_underscore)
                        self.writeln('const NekDouble ' + str(new_tau) + ' = 1.0 / (' + str(alpha) + ' + ' + str(beta) + ');')  
                        self.taus.append(new_tau)
            elif self.state_var_type[str(var)] == 'm_concentrations':
                continue
        self.writeln()  

       #write the infs
        self.output_comment('Calculating the inf-values:')
        for i,var in enumerate(self.state_vars):
            var_name = self.code_name(var, True)
            var_actual = str(var)[33:str(var).rfind('__')]
            if 'gate' in str(var_name):
                after_underscore = var_name[var_name.rfind('__')+2:]
            

            if 'membrane__V' in str(var):
                continue
            
            if self.state_var_type[str(var)] == 'm_gates':                  
                if filter(lambda element: var_actual in element,self.infs):
                    self.output_comment('The inf value for ' + str(var_actual) + ' should already be calculated in the mathematics')
                else:
                    if filter(lambda element: var_actual in element,self.alphas) and filter(lambda element: var_actual in element,self.betas):
                        alpha = filter(lambda element: var_actual in element,self.alphas)[0]
                        beta = filter(lambda element: var_actual in element,self.betas)[0]
                        new_inf = 'var_' + str(var_actual) + '__' + str(after_underscore) + '_inf'
                        self.writeln('const NekDouble ' + str(new_inf) + ' = ' + str(alpha) + ' / (' + str(alpha) + ' + ' + str(beta) + ');')   
                        self.infs.append(new_inf)
            elif self.state_var_type[str(var)] == 'm_concentrations':
                continue
        self.writeln()

        #debugging help by writing out the taus, infs and dictionary of variable types
        # self.writeln('Taus: ' + str(self.taus))
        # self.writeln('Infs: ' + str(self.infs))
        # self.writeln('Var Types: ' + str(self.state_var_type)+ '\n')

 
        #writing the outarrays
        self.output_comment('Writing the outarrays:')
        m_gate_tau_counter = 0

        for i,var in enumerate(self.state_vars):
            var_name = self.code_name(var, True)
            var_actual = str(var)[33:str(var).rfind('__')]
            if 'gate' in str(var_name):
                after_underscore = var_name[var_name.rfind('__')+2:]



            if 'membrane__V' in str(var):
                self.writeln('outarray[',i,'][i] = ' + str(var_name) + ';')

            if self.state_var_type[str(var)] == 'm_gates': 
                inf_value = filter(lambda element: var_actual in element,self.infs)[0]
                tau_value = filter(lambda element: var_actual in element,self.taus)[0]
                self.writeln('outarray[',i,'][i] = ' + inf_value + ';')
                self.writeln('m_gates_tau[',m_gate_tau_counter,'][i] = ' + tau_value + ';')
                m_gate_tau_counter += 1

            if self.state_var_type[str(var)] == 'm_concentrations': 
                self.writeln('outarray[',i,'][i] = ' + str(var_name) + ';')

            # self.writeln()



        self.close_block()
        self.close_block()

        
    def output_derivative_calculations(self, state_vars, assign_rY=False, extra_nodes=set(),
                                       extra_table_nodes=set()):
        """
        This is used by self.output_evaluate_y_derivatives and self.output_rush_larsen_mathematics
        to compute the derivatives (and any extra nodes, if given).  It contains the special logic
        to obey the mSetVoltageDerivativeToZero member variable in the generated code.
        Returns a nodeset containing the equations output.
        """
        # Work out what equations are needed to compute the derivatives
        
        #this is important and might be where the variables come from
        #derivs creates a set of lists of tuples, where each tuple is a state_var with the free_var
        derivs = set(map(lambda v: (v, self.free_vars[0]), state_vars))
        # for var in derivs:
        #     self.writeln(var)


        # self.writeln(derivs) #this line shows that set

        if self.v_variable in state_vars: #what this does is take out the v_variable into its own tuple and calls it dvdt

            dvdt = (self.v_variable, self.free_vars[0])
            derivs.remove(dvdt) #907: Consider dV/dt separately
        else:
            dvdt = None

        #this if function has been outcommented because i think it might not actually be needed in nektar    
        # if self.use_chaste_stimulus:
        #     i_stim = [self.doc._cml_config.i_stim_var]
        # else:
        i_stim = []
    

        nonv_nodeset = self.calculate_extended_dependencies(derivs|extra_nodes, prune_deps=i_stim)

        if dvdt:
            if self.use_data_clamp:
                prune = set([self.config.i_data_clamp_data]) | nonv_nodeset
            else:
                prune = nonv_nodeset
            # self.writeln('uses dvdt')
            v_nodeset = self.calculate_extended_dependencies([dvdt], prune=prune, prune_deps=i_stim)
        else:
            v_nodeset = set()
        # State variable inputs
        all_nodes = nonv_nodeset|v_nodeset
        self.output_state_assignments(assign_rY=assign_rY, nodeset=all_nodes)
        self.writeln()
        table_index_nodes_used = self.calculate_lookup_table_indices(all_nodes|extra_table_nodes, self.code_name(self.free_vars[0]))
        self.output_comment('Mathematics')
        #907: Declare dV/dt
        if dvdt:
            self.writeln(self.TYPE_DOUBLE, self.code_name(self.v_variable, ode=True), self.STMT_END)
        # Output mathematics required for non-dV/dt derivatives (which may include dV/dt)
        # self.writeln('Test Begin')
        self.NODESET = nonv_nodeset - table_index_nodes_used
        self.output_equations(nonv_nodeset - table_index_nodes_used)
        # self.writeln('Test End')



        self.writeln()
        
        #907: Calculation of dV/dt
        # if dvdt:
        # self.writeln('if (mSetVoltageDerivativeToZero)')
        # self.open_block()
        # self.writeln(self.code_name(self.v_variable, ode=True), self.EQ_ASSIGN, '0.0', self.STMT_END)
        # self.close_block(blank_line=False)
        # self.writeln('else')
        # self.open_block()

        # self.writeln('Test')
        # self.writeln(v_nodeset)
        
        # to_remove = set()
        # for var in v_nodeset:
        #     # self.writeln(var)
        #     if ',stim_' in str(var) or 'membrane,time' in str(var) and not 'membrane,i_st' in str(var):
        #         to_remove.add(var)
        
        # # self.writeln(to_remove)
        
        # for var in to_remove:
        #     v_nodeset.remove(var)

        self.writeln(self.TYPE_CONST_DOUBLE,'var_chaste_interface__membrane__I_stim = 0.0;')
        # self.writeln(self.TYPE_CONST_DOUBLE,'var_membrane__I_stim = var_chaste_interface__membrane__I_stim; // microA_per_cm2')
        self.output_equations(v_nodeset - table_index_nodes_used)
        # self.writeln('Test End')

        # self.close_block()
        
        self.writeln()
        
        return all_nodes | table_index_nodes_used

    def output_backward_euler_mathematics(self):
        """Output the mathematics methods used in a backward Euler cell.

        Outputs ComputeResidual, ComputeJacobian,
        UpdateTransmembranePotential and ComputeOneStepExceptVoltage.
        """
        dt_name = 'mDt'
        #model_dt = self.varobj(self.model.solver_info.dt)
        if self.nonlinear_system_size > 0:
            # Residual
            ##########
            argsize = '[' + str(self.nonlinear_system_size) + ']'
            self.output_method_start('ComputeResidual',
                                     [self.TYPE_DOUBLE + self.code_name(self.free_vars[0]),
                                      self.TYPE_CONST_DOUBLE + 'rCurrentGuess' + argsize,
                                      self.TYPE_DOUBLE + 'rResidual' + argsize],
                                     'void', access='public')
            self.open_block()
            # Output mathematics for computing du/dt for each nonlinear state var u
            nodes = map(lambda u: (u, self.free_vars[0]), self.nonlinear_system_vars)
            nodeset = self.calculate_extended_dependencies(nodes, prune_deps=[self.doc._cml_config.i_stim_var])
            self.output_state_assignments(exclude_nonlinear=True, nodeset=nodeset)
            self.output_nonlinear_state_assignments(nodeset=nodeset)
            table_index_nodes_used = self.calculate_lookup_table_indices(nodeset, self.code_name(self.free_vars[0]))
            self.output_equations(nodeset - table_index_nodes_used)
            self.writeln()
            # Fill in residual
            for i, var in enumerate(self.state_vars):
                try:
                    j = self.nonlinear_system_vars.index(var)
                except ValueError:
                    j = -1
                if j != -1:
                    self.writeln('rResidual[', j, '] = rCurrentGuess[', j, '] - rY[', i, '] - ',
                                 dt_name, '*', self.code_name(var, ode=True), self.STMT_END)
            self.close_block()
            
            # Jacobian
            ##########
            self.output_method_start('ComputeJacobian',
                                     [self.TYPE_DOUBLE + self.code_name(self.free_vars[0]),
                                      self.TYPE_CONST_DOUBLE + 'rCurrentGuess' + argsize,
                                      self.TYPE_DOUBLE + 'rJacobian' + argsize + argsize],
                                     'void', access='public')
            self.open_block()
            # Mathematics that the Jacobian depends on
            used_vars = set()
            for entry in self.model.solver_info.jacobian.entry:
                used_vars.update(self._vars_in(entry.math))
            nodeset = self.calculate_extended_dependencies(used_vars, prune_deps=[self.doc._cml_config.i_stim_var])
            self.output_state_assignments(exclude_nonlinear=True, nodeset=nodeset)
            self.output_nonlinear_state_assignments(nodeset=nodeset)
            self.writeln(self.TYPE_CONST_DOUBLE, self.code_name(self.config.dt_variable), self.EQ_ASSIGN, dt_name, self.STMT_END, '\n')
            table_index_nodes_used = self.calculate_lookup_table_indices(nodeset|set(map(lambda e: e.math, self.model.solver_info.jacobian.entry)), self.code_name(self.free_vars[0]))
            self.output_equations(nodeset - table_index_nodes_used)
            self.writeln()
            # Jacobian entries
            for entry in self.model.solver_info.jacobian.entry:
                var_i, var_j = entry.var_i, entry.var_j
                i = self.nonlinear_system_vars.index(self.varobj(var_i))
                j = self.nonlinear_system_vars.index(self.varobj(var_j))
                self.writeln('rJacobian[', i, '][', j, '] = ', nl=False)
                entry_content = list(entry.math.xml_element_children())
                assert len(entry_content) == 1, "Malformed Jacobian matrix entry: " + entry.xml()
                self.output_expr(entry_content[0], False)
                self.writeln(self.STMT_END, indent=False)
#            self.output_comment('Debugging')
#            self.writeln('#ifndef NDEBUG', indent=False)
#            self.writeln('for (unsigned i=0; i<', len(self.nonlinear_system_vars), '; i++)')
#            self.writeln('for (unsigned j=0; j<', len(self.nonlinear_system_vars), '; j++)', indent_offset=1)
#            self.writeln('EXCEPT_IF_NOT(!std::isnan(rJacobian[i][j]));', indent_offset=2)
#            self.writeln('//DumpJacobianToFile(', self.code_name(self.free_vars[0]),
#                         ', rCurrentGuess, rJacobian, rY);')
#            self.writeln('#endif // NDEBUG', indent=False)
            self.close_block()
        # The other methods are protected
        self.writeln_hpp('protected:', indent_offset=-1)
        
        # UpdateTransmembranePotential
        ##############################
        self.output_method_start('UpdateTransmembranePotential',
                                 [self.TYPE_DOUBLE + self.code_name(self.free_vars[0])],
                                 'void', access='public')
        self.open_block()
        self.output_comment('Time units: ', self.free_vars[0].units)
        # Output mathematics to compute dV/dt
        nodes = [(self.state_vars[self.v_index], self.free_vars[0])]
        nodeset = self.calculate_extended_dependencies(nodes, prune_deps=[self.doc._cml_config.i_stim_var])
        self.output_state_assignments(nodeset=nodeset)
        table_index_nodes_used = self.calculate_lookup_table_indices(nodeset, self.code_name(self.free_vars[0]))
        self.output_equations(nodeset - table_index_nodes_used)
        # Update V
        self.writeln()
        self.writeln('rY[', self.v_index, '] += ', dt_name, '*',
                     self.code_name(self.state_vars[self.v_index], ode=True), self.STMT_END)
        self.close_block()

        # ComputeOneStepExceptVoltage
        #############################
        self.output_method_start('ComputeOneStepExceptVoltage',
                                 [self.TYPE_DOUBLE + self.code_name(self.free_vars[0])],
                                 'void', access='public')
        self.open_block()
        self.writeln(self.COMMENT_START, 'Time units: ',
                     self.free_vars[0].units)
        # Output mathematics to update linear state variables, using solver_info.linear_odes.
        # Also need to use output_equations for variables used in the update equations.
        linear_vars, update_eqns = [], {}
        used_vars = set() # NB: Also contains update equation if is a mathml_apply so table index generation works
        for u, t, update_eqn in SolverInfo(self.model).get_linearised_odes():
            assert t == self.free_vars[0]
            assert len(update_eqn) == 1
            update_eqn = update_eqn[0]
            linear_vars.append(u)
            update_eqns[id(u)] = update_eqn
            if not isinstance(update_eqn, mathml_cn): used_vars.add(update_eqn)
            used_vars.update(self._vars_in(update_eqn))
        # Output required equations for used variables
        nodeset = self.calculate_extended_dependencies(used_vars, prune_deps=[self.doc._cml_config.i_stim_var])
        self.output_state_assignments(nodeset=nodeset)
        if self.config.dt_variable in nodeset:
            self.writeln(self.TYPE_CONST_DOUBLE, self.code_name(self.config.dt_variable), self.EQ_ASSIGN,
                         dt_name, self.STMT_END, '\n')
        table_index_nodes_used = self.calculate_lookup_table_indices(nodeset, self.code_name(self.free_vars[0]))
        self.output_equations(nodeset - table_index_nodes_used)
        # Update state variables:
        #   rY[i] = (rY[i] + _g_j*dt) / (1 - _h_j*dt)
        self.writeln()
        linear_vars.sort(key=lambda v: v.fullname())
        for i, u in enumerate(linear_vars):
            j = self.state_vars.index(u)
            self.writeln('rY[', j, ']', self.EQ_ASSIGN, nl=False)
            self.output_expr(update_eqns[id(u)], False)
            self.writeln(self.STMT_END, indent=False)
        # Set up the Newton iteration, if needed
        self.writeln()
        if self.nonlinear_system_size > 0:
            self.writeln('double _guess[', self.nonlinear_system_size, '] = {', nl=False)
            comma = False
            idx_map = [0] * self.nonlinear_system_size
            for i, var in enumerate(self.state_vars):
                try:
                    j = self.nonlinear_system_vars.index(var)
                    idx_map[j] = i
                except ValueError:
                    pass
            for i in idx_map:
                if comma: self.write(',')
                else: comma = True
                self.write('rY[', i, ']')
            self.writeln('};', indent=False)
            # Solve
            CNS = 'CardiacNewtonSolver<%d,%s>' % (self.nonlinear_system_size, self.class_name)
            self.writeln(CNS, '* _p_solver = ', CNS, '::Instance();')
            self.writeln('_p_solver->Solve(*this, ', self.code_name(self.free_vars[0]), ', _guess);')
            # Update state
            for j, i in enumerate(idx_map):
                self.writeln('rY[', i, '] = _guess[', j, '];')
        self.close_block()
    
    def output_rush_larsen_mathematics(self):
        """Output the special methods needed for Rush-Larsen style cell models.
        
        We generate:
         * EvaluateEquations  evaluate the model derivatives and alpha/beta terms
         * ComputeOneStepExceptVoltage  does a Rush-Larsen update for eligible variables,
           and a forward Euler step for other non-V state variables
        """
        rl_vars = self.doc._cml_rush_larsen
        # EvaluateEquations
        ###################
        self.output_method_start('EvaluateEquations',
                                 [self.TYPE_DOUBLE + self.code_name(self.free_vars[0]),
                                  'std::vector<double> &rDY',
                                  'std::vector<double> &rAlphaOrTau',
                                  'std::vector<double> &rBetaOrInf'],
                                 'void', access='public')
        self.open_block()
        normal_vars = [v for v in self.state_vars if not v in rl_vars]
        nodes, table_nodes = set(), set()
        for _, alpha_or_tau, beta_or_inf, _ in rl_vars.itervalues():
            table_nodes.add(alpha_or_tau)
            nodes.update(self._vars_in(alpha_or_tau))
            table_nodes.add(beta_or_inf)
            nodes.update(self._vars_in(beta_or_inf))
        self.output_derivative_calculations(normal_vars, True, nodes, table_nodes)
        # Now assign input vectors
        for i, var in enumerate(self.state_vars):
            if var in rl_vars:
                # Fill in rAlphaOrTau & rBetaOrInf
                self.writeln(self.vector_index('rAlphaOrTau', i), self.EQ_ASSIGN, nl=False)
                self.output_expr(rl_vars[var][1], False)
                self.writeln(self.STMT_END, indent=False)
                self.writeln(self.vector_index('rBetaOrInf', i), self.EQ_ASSIGN, nl=False)
                self.output_expr(rl_vars[var][2], False)
                self.writeln(self.STMT_END, indent=False)
            else:
                # Fill in rDY
                self.writeln(self.vector_index('rDY', i), self.EQ_ASSIGN, self.code_name(var, True), self.STMT_END)
        self.close_block()
        
        # ComputeOneStepExceptVoltage
        #############################
        self.output_method_start('ComputeOneStepExceptVoltage',
                                 ['const std::vector<double> &rDY',
                                  'const std::vector<double> &rAlphaOrTau',
                                  'const std::vector<double> &rBetaOrInf'],
                                 'void', access='public')
        self.open_block()
        self.writeln('std::vector<double>& rY = rGetStateVariables();')
        for i, var in enumerate(self.state_vars):
            if var in rl_vars:
                # Rush-Larsen update
                conv = rl_vars[var][3] or ''
                if conv: conv = '*' + str(conv)
                if rl_vars[var][0] == 'ab':
                    # Alpha & beta formulation
                    self.open_block()
                    self.writeln(self.TYPE_CONST_DOUBLE, 'tau_inv = rAlphaOrTau[', i, '] + rBetaOrInf[', i, '];')
                    self.writeln(self.TYPE_CONST_DOUBLE, 'y_inf = rAlphaOrTau[', i, '] / tau_inv;')
                    self.writeln('rY[', i, '] = y_inf + (rY[', i, '] - y_inf)*exp(-mDt', conv, '*tau_inv);')
                    self.close_block(blank_line=False)
                else:
                    # Tau & inf formulation
                    self.writeln('rY[', i, '] = rBetaOrInf[', i, '] + (rY[', i, '] - rBetaOrInf[', i, '])',
                                 '*exp(-mDt', conv, '/rAlphaOrTau[', i, ']);')
            elif var is not self.v_variable:
                # Forward Euler update
                self.writeln('rY[', i, '] += mDt * rDY[', i, '];')
        self.close_block()
    
    #Megan E. Marsh, Raymond J. Spiteri 
    #Numerical Simulation Laboratory 
    #University of Saskatchewan 
    #December 2011 
    #Partial support provided by research grants from 
    #the National Science and Engineering Research 
    #Council (NSERC) of Canada and the MITACS/Mprime 
    #Canadian Network of Centres of Excellence.
    def output_derivative_calculations_grl(self, var, assign_rY=False, extra_nodes=set(), extra_table_nodes=set()):
        """This is used by self.output_grl?_mathematics to get equations for each variable separately.

        Returns a node set with the equations output.
        """
        # Work out what equations are needed to compute the derivative of var
        if var in self.state_vars:
            dvardt = (var, self.free_vars[0])
            var_nodeset = self.calculate_extended_dependencies([dvardt])
        else:
            var_nodeset = set()
        # State variable inputs
        self.output_state_assignments(nodeset=var_nodeset, assign_rY=assign_rY)
        self.writeln()
        table_index_nodes_used = self.calculate_lookup_table_indices(var_nodeset, self.code_name(self.free_vars[0]))
        self.output_comment('Mathematics')
        self.output_equations(var_nodeset - table_index_nodes_used)
        return var_nodeset | table_index_nodes_used

    def find_grl_partial_derivatives(self):
        """If we have analytic Jacobian information available from Maple, find the terms needed for GRL methods.

        This caches where the diagonal entries are in the matrix, indexed by the state variable objects currently in use,
        since the entries in the matrix may reference non-partially-evaluated variables.
        """
        if not hasattr(self, 'jacobian_diagonal'):
            self.jacobian_diagonal = {}
        if self.use_analytic_jacobian and not self.jacobian_diagonal:
            for entry in self.model.solver_info.jacobian.entry:
                if entry.var_i == entry.var_j:
                    # It's a diagonal entry
                    var = self.varobj(entry.var_i).get_source_variable(recurse=True)
                    assert var in self.state_vars, "Jacobian diagonal entry is not in the state vector: " + entry.xml()
                    entry_content = list(entry.math.xml_element_children())
                    assert len(entry_content) == 1, "Malformed Jacobian entry: " + entry.xml()
                    self.jacobian_diagonal[var] = entry_content[0]

    def output_grl_compute_partial(self, i, var):
        """Compute the partial derivative of f(var) wrt var, the i'th variable in the state vector.

        This uses an analytic Jacobian if available; otherwise it approximates using finite differences.
        """
        self.output_method_start('EvaluatePartialDerivative'+str(i),
                                 [self.TYPE_DOUBLE + self.code_name(self.free_vars[0]),
                                  'std::vector<double>& rY', 'double delta', 'bool forceNumerical'],
                                 'double', access='public', defaults=['', '', '', 'false'])
        self.open_block()
        self.writeln('double partialF;')
        if self.jacobian_diagonal:
            # Work out what equations are needed to compute the analytic derivative
            self.writeln('if (!forceNumerical && this->mUseAnalyticJacobian)')
            self.open_block()
            entry = self.jacobian_diagonal[var]
            nodeset = self.calculate_extended_dependencies(self._vars_in(entry))
            self.output_state_assignments(nodeset=nodeset, assign_rY=False)
            table_index_nodes_used = self.calculate_lookup_table_indices(nodeset|set([entry]), self.code_name(self.free_vars[0]))
            self.output_equations(nodeset)
            # Calculate the derivative
            self.writeln('partialF = ', nl=False)
            self.output_expr(entry, paren=False)
            self.writeln(self.STMT_END, indent=False)
            self.close_block(blank_line=False)
            self.writeln('else')
            self.open_block()
        # Numerical approximation
        self.writeln('const double y_save = rY[', i, '];')
        self.writeln('rY[', i, '] += delta;')
        self.writeln('const double temp = EvaluateYDerivative', i, '(', self.code_name(self.free_vars[0]), ', rY);')
        self.writeln('partialF = (temp-mEvalF[', i, '])/delta;')
        self.writeln('rY[', i, '] = y_save;')
        if self.jacobian_diagonal:
            self.close_block(blank_line=False)
        self.writeln('return partialF;')
        self.close_block()

    #Megan E. Marsh, Raymond J. Spiteri 
    #Numerical Simulation Laboratory 
    #University of Saskatchewan 
    #December 2011 
    #Partial support provided by research grants from 
    #the National Science and Engineering Research 
    #Council (NSERC) of Canada and the MITACS/Mprime 
    #Canadian Network of Centres of Excellence.
    def output_grl1_mathematics(self):
        """Output the special methods needed for GRL1 style cell models.

        We generate:
         * UpdateTransmembranePotential update V_m
         * ComputeOneStepExceptVoltage  does a GRL1 update for variables except voltage
         * EvaluateYDerivativeI for each variable I
        """
        self.find_grl_partial_derivatives()
        ########################################################UpdateTransmembranePotential
        self.output_method_start('UpdateTransmembranePotential',
                                 [self.TYPE_DOUBLE + self.code_name(self.free_vars[0])],
                                 'void', access='public')
        self.open_block()
        self.writeln('std::vector<double>& rY = rGetStateVariables();')
        self.writeln('unsigned v_index = GetVoltageIndex();')
        self.writeln('const double delta = 1e-8;')
        self.writeln()
        # Compute partial derivative of dV wrt V
        self.writeln(self.TYPE_DOUBLE, self.code_name(self.v_variable, ode=True), self.STMT_END)
        self.output_derivative_calculations_grl(self.v_variable)
        self.writeln()
        self.writeln('double evalF = ', self.code_name(self.v_variable, ode=True), self.STMT_END)
        self.writeln('mEvalF[', self.v_index, '] = ', self.code_name(self.v_variable, ode=True), self.STMT_END)
        self.writeln('double partialF = EvaluatePartialDerivative', self.v_index, '(', self.code_name(self.free_vars[0]), ', rY, delta, true);')
        self.writeln('if (fabs(partialF) < delta)')
        self.open_block()
        self.writeln('rY[v_index] += evalF*mDt;')
        self.close_block(False)
        self.writeln('else')
        self.open_block()
        self.writeln('rY[v_index] += (evalF/partialF)*(exp(partialF*mDt)-1.0);')
        self.close_block()
        self.close_block()

        #########################################################ComputeOneStepExceptVoltage
        self.output_method_start('ComputeOneStepExceptVoltage',
                                 [self.TYPE_DOUBLE + self.code_name(self.free_vars[0])],
                                 'void', access='public')
        self.open_block()
        # Set up variables
        self.writeln('std::vector<double>& rY = rGetStateVariables();')
        self.writeln('const double delta = 1e-8;')
        self.writeln()

        # Evaluate RHS of equations (except dV/dt)
        non_v_vars = self.state_vars[:]
        if self.v_variable in non_v_vars:
            non_v_vars.remove(self.v_variable)
        self.output_derivative_calculations(non_v_vars)

        # Compute partial derivatives (for non-V)
        for i, var in enumerate(self.state_vars):
            if var is not self.v_variable:
                self.writeln('mEvalF[', i, '] = ', self.code_name(var, ode=True), self.STMT_END)
                self.writeln('mPartialF[', i, '] = EvaluatePartialDerivative', i, '(', self.code_name(self.free_vars[0]), ', rY, delta);')

        # Do the GRL updates
        for i, var in enumerate(self.state_vars):
            if var is not self.v_variable:
                self.open_block()
                self.writeln('if (fabs(mPartialF[', i, ']) < delta)')
                self.open_block()
                self.writeln('rY[', i, '] += mDt*', self.code_name(var, True), ';')
                self.close_block(False)
                self.writeln('else')
                self.open_block()
                self.writeln('rY[', i, '] += (', self.code_name(var, True), '/mPartialF[', i, '])*(exp(mPartialF[', i, ']*mDt)-1.0);')
                self.close_block()
                self.close_block()
        self.close_block()

        #########################################################Evaluate each equation
        for i, var in enumerate(self.state_vars):
            self.output_method_start('EvaluateYDerivative'+str(i),
                                     [self.TYPE_DOUBLE + self.code_name(self.free_vars[0]),
                                      'std::vector<double>& rY'],
                                     'double', access='public')
            self.open_block()
            if var is self.v_variable:
                self.writeln(self.TYPE_DOUBLE, self.code_name(self.v_variable, ode=True), self.STMT_END)
            self.output_derivative_calculations_grl(var)
            self.writeln()
            self.writeln('return ', self.code_name(var, True), ';')
            self.close_block()

            self.output_grl_compute_partial(i, var)

    #Megan E. Marsh, Raymond J. Spiteri 
    #Numerical Simulation Laboratory 
    #University of Saskatchewan 
    #December 2011 
    #Partial support provided by research grants from 
    #the National Science and Engineering Research 
    #Council (NSERC) of Canada and the MITACS/Mprime 
    #Canadian Network of Centres of Excellence.     
    def output_grl2_mathematics(self):
        """Output the special methods needed for GRL2 style cell models.

        We generate:
         * Update TransmembranePotential update V_m
         * ComputeOneStepExceptVoltage  does a GRL2 update for variables except voltage
         * EvaluateYDerivativeI for each variable I
        """
        self.find_grl_partial_derivatives()
        ########################################################UpdateTransmembranePotential
        self.output_method_start('UpdateTransmembranePotential',
                                 [self.TYPE_DOUBLE + self.code_name(self.free_vars[0])],
                                 'void', access='public')
        self.open_block()
        self.writeln('std::vector<double>& rY = rGetStateVariables();')
        self.writeln('const unsigned v_index = GetVoltageIndex();')
        self.writeln('const double delta = 1e-8;')
        self.writeln('const double yinit = rY[v_index];')
        self.writeln()

        # Do the first half step
        self.writeln(self.TYPE_DOUBLE, self.code_name(self.v_variable, ode=True), self.STMT_END)
        self.output_derivative_calculations_grl(self.v_variable)
        self.writeln()
        self.writeln('double evalF = ', self.code_name(self.v_variable, ode=True), self.STMT_END)
        self.writeln('mEvalF[', self.v_index, '] = ', self.code_name(self.v_variable, ode=True), self.STMT_END)
        self.writeln('double partialF = EvaluatePartialDerivative', self.v_index, '(', self.code_name(self.free_vars[0]), ', rY, delta, true);')
        self.writeln('if (fabs(partialF) < delta)')
        self.open_block()
        self.writeln('rY[v_index] += 0.5*evalF*mDt;')
        self.close_block(False)
        self.writeln('else')
        self.open_block()
        self.writeln('rY[v_index] += (evalF/partialF)*(exp(partialF*0.5*mDt)-1.0);')
        self.close_block()

        # Do the second half step
        self.writeln('rY[v_index] = yinit;')
        self.writeln('evalF = EvaluateYDerivative', self.v_index, '(', self.code_name(self.free_vars[0]), ', rY);')
        self.writeln('mEvalF[', self.v_index, '] = evalF;')
        self.writeln('partialF = EvaluatePartialDerivative', self.v_index, '(', self.code_name(self.free_vars[0]), ', rY, delta, true);')
        self.writeln('if (fabs(partialF) < delta)')
        self.open_block()
        self.writeln('rY[v_index] = yinit + evalF*mDt;')
        self.close_block(False)
        self.writeln('else')
        self.open_block()
        self.writeln('rY[v_index] = yinit + (evalF/partialF)*(exp(partialF*mDt)-1.0);')
        self.close_block()
        self.close_block() # End method

        #########################################################ComputeOneStepExceptVoltage
        self.output_method_start('ComputeOneStepExceptVoltage',
                                 [self.TYPE_DOUBLE + self.code_name(self.free_vars[0])],
                                 'void', access='public')
        self.open_block()
        # Set up variables
        self.writeln('std::vector<double>& rY = rGetStateVariables();')
        self.writeln('const double delta=1e-8;')
        self.writeln('const unsigned size = GetNumberOfStateVariables();')
        self.writeln('mYInit = rY;')
        self.writeln('double y_save;')
        self.writeln()

        # Calculate partial derivatives
        self.output_derivative_calculations(self.state_vars)
        for i, var in enumerate(self.state_vars):
            self.writeln(self.vector_index('mEvalF', i), self.EQ_ASSIGN, self.code_name(var, True), self.STMT_END)
        self.writeln()
        for i, var in enumerate(self.state_vars):
            if var is not self.v_variable:
                self.writeln('mPartialF[', i, '] = EvaluatePartialDerivative', i, '(', self.code_name(self.free_vars[0]), ', rY, delta);')

        # Update all variables
        self.writeln('for (unsigned var=0; var<size; var++)')
        self.open_block()
        self.writeln('if (var == ', self.v_index, ') continue;')
        self.writeln('if (fabs(mPartialF[var]) < delta)')
        self.open_block()
        self.writeln('rY[var] = mYInit[var] + 0.5*mDt*mEvalF[var];')
        self.close_block(False)
        self.writeln('else')
        self.open_block()
        self.writeln('rY[var] = mYInit[var] + (mEvalF[var]/mPartialF[var])*(exp(mPartialF[var]*0.5*mDt)-1.0);')
        self.close_block()
        self.close_block()
        self.writeln()

        # Determine new partial derivatives
        for i, var in enumerate(self.state_vars):
            if var is not self.v_variable:
                self.writeln()
                self.writeln('y_save = rY[', i, '];')
                self.writeln('rY[', i, '] = mYInit[', i, '];')
                self.writeln('mEvalF[', i, '] = EvaluateYDerivative', i, '(', self.code_name(self.free_vars[0]), ', rY);')
                self.writeln('mPartialF[', i, '] = EvaluatePartialDerivative', i, '(', self.code_name(self.free_vars[0]), ', rY, delta);')
                self.writeln('rY[', i, '] = y_save;')

        # Update all variables
        self.writeln('for (unsigned var=0; var<size; var++)')
        self.open_block()
        self.writeln('if (var == ', self.v_index, ') continue;')
        self.writeln('if (fabs(mPartialF[var]) < delta)')
        self.open_block()
        self.writeln('rY[var] = mYInit[var] + mDt*mEvalF[var];')
        self.close_block(False)
        self.writeln('else')
        self.open_block()
        self.writeln('rY[var] = mYInit[var] + (mEvalF[var]/mPartialF[var])*(exp(mPartialF[var]*mDt)-1.0);')
        self.close_block()
        self.close_block()
        self.writeln()
        self.close_block() # End method

        #########################################################Evaluate each equation
        for i, var in enumerate(self.state_vars):
            self.output_method_start('EvaluateYDerivative'+str(i),
                                     [self.TYPE_DOUBLE + self.code_name(self.free_vars[0]),
                                      'std::vector<double>& rY'],
                                     'double', access='public')
            self.open_block()
            if var is self.v_variable:
                self.writeln(self.TYPE_DOUBLE, self.code_name(self.v_variable, ode=True), self.STMT_END)
            self.output_derivative_calculations_grl(var)
            self.writeln()
            self.writeln('return '+self.code_name(var, True)+';')
            self.close_block()

            self.output_grl_compute_partial(i, var)


    def output_model_attributes(self):
        """Output any named model attributes defined in metadata.

        Such attributes are given by compound RDF annotations:
          model --pycml:named-attribute--> bnode
          bnode --pycml:name--> Literal(Attribute name, string)
          bnode --pycml:value--> Literal(Attribute value, double)
        """
        model = self.model
        meta_id = model.cmeta_id
        attrs = []
        if meta_id:
            property = cellml_metadata.create_rdf_node(('pycml:named-attribute', NSS['pycml']))
            name_prop = cellml_metadata.create_rdf_node(('pycml:name', NSS['pycml']))
            value_prop = cellml_metadata.create_rdf_node(('pycml:value', NSS['pycml']))
            source = cellml_metadata.create_rdf_node(fragment_id=meta_id)
            attr_nodes = cellml_metadata.get_targets(model, source, property)
            for node in attr_nodes:
                name = cellml_metadata.get_target(model, node, name_prop)
                value = cellml_metadata.get_target(model, node, value_prop)
                attrs.append((name, value))
        for name, value in attrs:
            self.writeln('this->mAttributes["', name, '"] = ', value, ';')
        if attrs:
            self.writeln()

    def output_bottom_boilerplate(self):
        """Output bottom boilerplate.

        End class definition, output ODE system information (to .cpp) and
        serialization code (to .hpp), and end the file.
        """
        # End main class
        self.set_indent(offset=0)
        # self.writeln_hpp('};\n\n')
        # ODE system information

        self.writeln(self.TYPE_VOID, self.class_name, '::v_GenerateSummary(SummaryList& s)')
        self.open_block()
        self.writeln('SolverUtils::AddSummaryItem(s, "', 'Cell model' , '", "', self.class_name , '");')
        self.close_block()    

        #self.writeln('template<>')
        self.writeln(self.TYPE_VOID, self.class_name,
                     '::v_SetInitialConditions()')
        self.open_block()
        # self.writeln('this->mSystemName', self.EQ_ASSIGN, '"', self.model.name, '"', self.STMT_END)
        # self.writeln('this->mFreeVariableName', self.EQ_ASSIGN,
        #              '"', self.var_display_name(self.free_vars[0]), '"', self.STMT_END)
        # self.writeln('this->mFreeVariableUnits', self.EQ_ASSIGN,
        #              '"', self.free_vars[0].units, '"', self.STMT_END)
        # self.writeln()
        def output_var(vector, var):
            self.writeln('this->m', vector, 'Names.push_back("', self.var_display_name(var), '");')
            self.writeln('this->m', vector, 'Units.push_back("', var.units, '");')
        
        for i,var in enumerate(self.state_vars):
            #output_var('Variable', var)
            init_val = getattr(var, u'initial_value', None)
            
            if init_val is None:
                init_comm = ' // Value not given in model'
                # Don't want compiler error, but shouldn't be a real number
                init_val = self.NOT_A_NUMBER
            else:
                init_comm = ''
                #21
            self.writeln('Vmath::Fill(m_nq, ', float(init_val), ',' ,  ' '*(21-len(init_val)) , 'm_cellSol[', i ,'],  1);',
                       init_comm)
            
        # Model parameters
        for var in self.cell_parameters:
            if var.get_type() == VarTypes.Constant:
                output_var('Parameter', var)
                self.writeln()
        # Derived quantities
        for var in self.derived_quantities:
            output_var('DerivedQuantity', var)
            self.writeln()
        self.output_model_attributes()
        #self.writeln('this->mInitialised = true;')
        self.writeln('}',indent_level=1)
        #self.writeln()
        # Serialization
        # if self.include_serialization:
        #     self.output_comment('Needs to be included last', subsidiary=True)
        #     self.writeln_hpp('#include "SerializationExportWrapper.hpp"')
        #     self.writeln_hpp('CHASTE_CLASS_EXPORT(', self.class_name, ')')
        #     self.output_comment('Serialization for Boost >= 1.36')
        #     self.writeln('#include "SerializationExportWrapperForCpp.hpp"')
        #     self.writeln('CHASTE_CLASS_EXPORT(', self.class_name, ')')
        #     self.writeln_hpp()
        #     self.writeln_hpp('namespace boost')
        #     self.open_block(subsidiary=True)
        #     self.writeln_hpp('namespace serialization')
        #     self.open_block(subsidiary=True)
        #     # Save
        #     self.writeln_hpp('template<class Archive>')
        #     self.writeln_hpp('inline void save_construct_data(')
        #     self.writeln_hpp('Archive & ar, const ', self.class_name,
        #                      ' * t, const unsigned int fileVersion)',
        #                      indent_offset=1)
        #     self.open_block(subsidiary=True)
        #     self.writeln_hpp('const boost::shared_ptr<AbstractIvpOdeSolver> p_solver = t->GetSolver();')
        #     self.writeln_hpp('const boost::shared_ptr<AbstractStimulusFunction> p_stimulus = t->GetStimulusFunction();')
        #     self.writeln_hpp('ar << p_solver;')
        #     self.writeln_hpp('ar << p_stimulus;')
        #     self.close_block(subsidiary=True)
        #     # Load
        #     self.writeln_hpp('template<class Archive>')
        #     self.writeln_hpp('inline void load_construct_data(')
        #     self.writeln_hpp('Archive & ar, ', self.class_name,
        #                      ' * t, const unsigned int fileVersion)',
        #                      indent_offset=1)
        #     self.open_block(subsidiary=True)
        #     self.writeln_hpp('boost::shared_ptr<AbstractIvpOdeSolver> p_solver;')
        #     self.writeln_hpp('boost::shared_ptr<AbstractStimulusFunction> p_stimulus;')
        #     self.writeln_hpp('ar >> p_solver;')
        #     self.writeln_hpp('ar >> p_stimulus;')
        #     self.writeln_hpp('::new(t)', self.class_name, '(p_solver, p_stimulus);')
        #     self.close_block(subsidiary=True)
        #     self.close_block(subsidiary=True)
        #     self.close_block(subsidiary=True)
        # if self.dynamically_loadable:
        #     # Write the C function to create instances of this cell model
        #     self.writeln('extern "C"')
        #     self.open_block()
        #     self.writeln('AbstractCardiacCellInterface* MakeCardiacCell(')
        #     self.writeln('boost::shared_ptr<AbstractIvpOdeSolver> pSolver,', indent_offset=2)
        #     self.writeln('boost::shared_ptr<AbstractStimulusFunction> pStimulus)', indent_offset=2)
        #     self.open_block()
        #     self.writeln('return new ', self.class_name, '(pSolver, pStimulus);')
        #     self.close_block()
        #     self.close_block()
        # End file
        #self.writeln_hpp('#endif // ', self.include_guard)
        self.writeln('}', indent_level=0)
        return

    def output_lhs(self, expr):
        """Output the left hand side of an assignment expression."""
        if expr.localName == 'ci':
            self.output_variable(expr)
        elif expr.operator().localName == 'diff':
            ci_elt = expr.operands().next()
            self.output_variable(ci_elt, ode=True)
        return

    def output_variable(self, ci_elt, ode=False):
        """Output a ci element, i.e. a variable lookup."""
        if hasattr(ci_elt, '_cml_variable') and ci_elt._cml_variable:
            self.write(self.code_name(ci_elt.variable, ode=ode))
        else:
            # This ci element doesn't have all the extra annotations.  It is a fully
            # qualified name though.  This is typically because PE has been done.
            prefix = ['var_', 'd_dt_'][ode]
            varname = unicode(ci_elt)
            try:
                var = self.varobj(varname)
            except KeyError:
                var = None
            if var:
                self.write(self.code_name(var, ode=ode))
            else:
                # Assume it's a suitable name
                self.write(prefix + varname)
        return

    def output_function(self, func_name, args, *posargs, **kwargs):
        """Override base class method for special case of abs with 2 arguments.
        
        This comes from Maple's Jacobians, and should generate signum of the second argument.
        """
        args = list(args)
        if func_name == 'fabs' and len(args) == 2:
            super(CellMLToNektarTranslator, self).output_function('Signum', [args[1]], *posargs, **kwargs)
        else:
            super(CellMLToNektarTranslator, self).output_function(func_name, args, *posargs, **kwargs)
    
    @staticmethod
    def get_current_units_options(model):
        """
        Return a list of units objects that give the possibilities for the dimensions
        of transmembrane ionic currents.
        """
        chaste_units = cellml_units.create_new(
            model, 'uA_per_cm2',
            [{'units': 'ampere', 'prefix': 'micro'},
             {'units': 'metre', 'prefix': 'centi', 'exponent': '-2'}])
        microamps = cellml_units.create_new(model, u'microamps',
                                            [{'units':'ampere', 'prefix':'micro'}])
        A_per_F = cellml_units.create_new(model, 'A_per_F',
                                          [{'units': 'ampere'},
                                           {'units': 'farad', 'exponent': '-1'}])
        return [chaste_units, microamps, A_per_F]

    # Name in CellML for the variable representing Chaste's membrane capacitance
    MEMBRANE_CAPACITANCE_NAME = u'chaste_membrane_capacitance'
    
    # Name of the component added to interface the model to Chaste
    INTERFACE_COMPONENT_NAME = u'chaste_interface'

    @staticmethod
    def add_special_conversions(converter, comp):
        """Add special units conversions for ionic currents.
        
        Adds conversions for the two other common conventions to/from the units expected by Chaste,
        uA/cm^2.  The cases are:
        
        1. Current in amps/farads.
           In this case we convert to uA/uF then multiply by Chaste's value
           for the membrane capacitance (in uF/cm^2).
        2. Current in amps, capacitance in farads.
           We assume the cell model conceptually represents a cell, and hence
           that its membrane capacitance is supposed to represent the same
           thing as Chaste's.  Thus convert current to uA, capacitance to uF,
           and return current/capacitance * Chaste's capacitance.
        
        comp is a component to which we should add any necessary variables, i.e. Chaste's capacitance.
        """
        klass = CellMLToNektarTranslator
        model = converter.model
        # Variables needed by some conversions
        model_Cm = model.get_config('Cm_variable')
        uF_per_cm2 = cellml_units.create_new(model, 'uF_per_cm2',
                                             [{'units': 'farad', 'prefix': 'micro'},
                                              {'units': 'metre', 'prefix': 'centi', 'exponent': '-2'}])
        Chaste_Cm = converter.add_variable(comp, klass.MEMBRANE_CAPACITANCE_NAME, uF_per_cm2)
        model._cml_Chaste_Cm = Chaste_Cm # Record for use in code_name
        # Add the conversions
        chaste_units, microamps, A_per_F = klass.get_current_units_options(model)
        converter.add_special_conversion(A_per_F, chaste_units,
                                         lambda expr: converter.times_rhs_by(expr, Chaste_Cm))
        converter.add_special_conversion(chaste_units, A_per_F,
                                         lambda expr: converter.divide_rhs_by(expr, Chaste_Cm))
        if model_Cm:
            converter.add_special_conversion(microamps, chaste_units,
                    lambda expr: converter.times_rhs_by(converter.divide_rhs_by(expr, model_Cm),
                                                        Chaste_Cm))
            converter.add_special_conversion(chaste_units, microamps,
                    lambda expr: converter.divide_rhs_by(converter.times_rhs_by(expr, model_Cm),
                                                         Chaste_Cm))
        
    @staticmethod
    def generate_interface(doc, solver_info):
        """Generate an interface component connecting the model to Chaste.
        
        On return from this method, Chaste code will only need to interact with variables in
        the new interface component.  It will contain the transmembrane potential, the ionic
        and stimulus currents, the simulation time, and the derivatives.
        
        It may also contain other variables depending on the model, for example the intracellular
        calcium concentration (if annotated), modifiable parameters, and derived quantities.
        
        If the --convert-interfaces option has been supplied, units conversion will then be
        performed on this component, ensuring that all these variables are in the units expected
        by Chaste and linked by suitable conversions to the rest of the model.
        
        Note that if partial evaluation is then performed, the model will be collapsed into a
        single component.  However, the interface will still be preserved in the correct units.
        """
        model = doc.model
        config = doc._cml_config
        klass = CellMLToNektarTranslator
        # Create Chaste units definitions
        ms = cellml_units.create_new(model, 'millisecond',
                                     [{'units': 'second', 'prefix': 'milli'}])
        mV = cellml_units.create_new(model, 'millivolt',
                                     [{'units': 'volt', 'prefix': 'milli'}])
        current_units, microamps = klass.get_current_units_options(model)[0:2]
        # The interface generator
        generator = processors.InterfaceGenerator(model, name=klass.INTERFACE_COMPONENT_NAME)
        iface_comp = generator.get_interface_component()
        # In case we need to convert initial values, we create the units converter here
        if config.options.convert_interfaces:
            warn_only = not config.options.fully_automatic and config.options.warn_on_units_errors
            notifier = NotifyHandler(level=logging.WARNING)
            logging.getLogger('units-converter').addHandler(notifier)
            converter = processors.UnitsConverter(model, warn_only, show_xml_context_only=True)
            klass.add_special_conversions(converter, iface_comp)
            generator.set_units_converter(converter)
        # And now specify the interface
        t = model.find_free_vars()[0]
        if not ms.dimensionally_equivalent(t.get_units()):
            # Oops!
            raise TranslationError('Time does not have dimensions of time')
        generator.add_input(t, ms)
        if doc.model.get_option('backward_euler'):
            # Backward Euler code generation requires access to the time step
            model_dt = solver_info.create_dt(generator, t.component, t.get_units())
            config.dt_variable = generator.add_input(model_dt, ms)
            config.dt_variable.set_pe_keep(True)
        elif doc.model.get_option('maple_output'):
            # CVODE Jacobians need to be able to scale for time too
            fake_dt = generator.add_variable(t.component, 'fake_dt', ms, initial_value='1.0')
            fake_dt._set_type(VarTypes.Constant)
            config.dt_variable = generator.add_input(fake_dt, t.get_units())
            config.dt_variable.set_is_modifiable_parameter(False)
            config.dt_variable.set_pe_keep(True)

        if config.options.use_chaste_stimulus and config.i_stim_var:
            # We need to make it a constant so add_input doesn't complain, then make it computed
            # again so that exposing metadata-annotated variables doesn't make it a parameter!
            generator.make_var_constant(config.i_stim_var, 0)
            config.i_stim_var = generator.add_input(config.i_stim_var, current_units,
                                                    annotate=False, convert_initial_value=False)
            generator.make_var_computed_constant(config.i_stim_var, 0)
            # Also convert variables that make up the default stimulus
            # Note: we vary in/out-put primarily to test units conversion of initial values
            def add_oxmeta_ioput(oxmeta_name, units, inout):
                var = doc.model.get_variable_by_oxmeta_name(oxmeta_name, throw=False)
                if var:
                    meth = getattr(generator, 'add_%sput' % inout)
                    newvar = meth(var, units, annotate=False)
                    newvar.set_pe_keep(True)
            for n in ['duration', 'period', 'offset', 'end']:
                add_oxmeta_ioput('membrane_stimulus_current_'+n, ms, 'in')
            add_oxmeta_ioput('membrane_stimulus_current_amplitude', current_units, 'out')

        if config.V_variable:
            config.V_variable = generator.add_input(config.V_variable, mV)
        ionic_vars = config.i_ionic_vars
        if ionic_vars:
            i_ionic = generator.add_output_function('i_ionic', 'plus', ionic_vars, current_units)
            config.i_ionic_vars = [i_ionic]

        if doc.model.get_option('use_data_clamp'):
            assert config.V_variable and ionic_vars
            # Create g_clamp
            conductance_units = current_units.quotient(mV).simplify()
            i_data_clamp_conductance = generator.add_variable(iface_comp, 'membrane_data_clamp_current_conductance', conductance_units, initial_value='0.0')
            i_data_clamp_conductance._set_type(VarTypes.Constant)
            i_data_clamp_conductance.set_pe_keep(True) # This prevents it becoming 'chaste_interface__membrane_data_clamp_current_conductance'
            config.i_data_clamp_conductance = generator.add_input(i_data_clamp_conductance, conductance_units)
            # Create V_clamp
            data_var = config.i_data_clamp_data = generator.add_variable(iface_comp, 'experimental_data_voltage', mV, initial_value='0.0')
            data_var._set_type(VarTypes.Constant)
            data_var.set_pe_keep(True)
            data_var._cml_code_name = 'GetExperimentalVoltageAtTimeT(%(time)s)'
            # Create the current: I = g_clamp * (V - V_clamp)
            current_var = config.i_data_clamp_current = generator.add_variable(iface_comp, 'membrane_data_clamp_current', current_units)
            current_var._set_type(VarTypes.Computed)
            current_var.set_is_derived_quantity(True)
            sub = mathml_apply.create_new(model, u'minus', [config.V_variable.name, data_var.name])
            times = mathml_apply.create_new(model, u'times', [config.i_data_clamp_conductance.name, sub])
            assign = mathml_apply.create_new(model, u'eq', [current_var.name, times])
            generator.add_expr_to_comp(iface_comp, assign)
            # Make dV/dt use the new current
            def process_ci(elt):
                # Add reference to new current after first existing ionic current
                ref = mathml_ci.create_new(model, local_current_var.name)
                elt.xml_parent.xml_insert_after(elt, ref)
            if hasattr(ionic_vars[0], '_cml_ref_in_dvdt'):
                local_current_var = generator.connect_variables(current_var, (ionic_vars[0]._cml_ref_in_dvdt.component.name, current_var.name))
                process_ci(ionic_vars[0]._cml_ref_in_dvdt)
            else:
                dVdt = config.V_variable.get_all_expr_dependencies()[0]
                local_current_var = generator.connect_variables(current_var, (config.V_variable.component.name, current_var.name))
                def process_ci_elts(elt):
                    """Recursively process any ci elements in the tree rooted at elt."""
                    if isinstance(elt, mathml_ci):
                        if elt.variable is ionic_vars[0]:
                            process_ci(elt)
                    else:
                        for child in getattr(elt, 'xml_children', []):
                            process_ci_elts(child)
                process_ci_elts(dVdt)

        # Finish up
        def errh(errors):
            raise TranslationError("Creation of Chaste interface component failed:\n  " + str(errors))
        generator.finalize(errh, check_units=False)
        # Apply units conversions just to the interface, if desired
        if config.options.convert_interfaces:
            converter.add_conversions_for_component(iface_comp)
            converter.finalize(errh, check_units=False)
            notifier.flush()
            logging.getLogger('units-converter').removeHandler(notifier)
            if notifier.messages:
                msg = 'Problems occurred converting model variables to Chaste units.\n'
                if ionic_vars and ionic_vars[0].get_units().dimensionally_equivalent(microamps):
                    msg += 'To convert the ionic currents for this model, '\
                           'the model membrane capacitance needs to be identified.'
                if config.options.fully_automatic:
                    raise TranslationError(msg)
                else:
                    print >>sys.stderr, msg