Hodgkin-Huxley 5: Interact with generated code

By the time you have worked through this tutorial you will be able to:

• Investigate and understand the contents of files created by the Generator; and

• Integrate generated code into a simple solver to run a simulation.

This tutorial assumes that you are comfortable with:

• Interacting with a model and its entities using the API;

• Using the Generator functionality to output files in C or Python; and

• The basic idea of numerical integration using Euler’s method.

C++ resources

• CMakeLists.txt The CMake file for building this tutorial;

• solveGeneratedCode.cpp Either the skeleton code, or ..

• solveGeneratedCode_completed.cpp the completed tutorial code;

• HodgkinHuxleyModel.cpp Generated implementation code from Tutorial 4; and

• HodgkinHuxleyModel.h Generated interface code from Tutorial 4.

Python resources

• solveGeneratedCode.py Either the skeleton code, or ..

• solveGeneratedCode_completed.py the completed tutorial code;

• HodgkinHuxleyModel.py Generated implementation code from Tutorial 4.

Step 1: Link to the generated code

The first step in interacting with the output from the code generation is including the generated files in the project. There are necessarily big differences between C++ and Python for this tutorial!

C++ instructionsPython instructions

1.a Enter the path to the generated header/interface *.h file in the #include block above.

1.b Add the name and path of the implementation *.cpp file in the CMakeLists.txt file, (or whatever your local toolchain requires).

1.c Open the implementation file *.cpp file and verify that the #include statement in line 3 has the filename of your interface *.h file. Amend if needed and close the file.

1.d Call cmake to create the Makefile. Call make -j to build the executable. Run the code so far to check that the libCellML versions match.

Show C++ snippet

    //  1.a 
    //      Enter the path to the generated header/interface *.h file in the #include block above.

    //  1.b 
    //      Add the name and path of the implementation *.cpp file in the CMakeLists.txt file, 
    //      (or whatever your local toolchain requires).

    //  1.c 
    //      Open the implementation file *.cpp file and verify that the #include statement in
    //      line 3 has the filename of your interface *.h file.  Amend if needed and close the file.

    //  1.d 
    //      Call cmake to create the Makefile.  Call make -j to build the executable.
    //      Run the code so far to check that the libCellML versions match.

Step 2: Access the variables in the generated files

Probably the best way to understand the contents of the generated files is o open them and look! The implementation file *.cpp has two types of items:

• information structures (in all-caps); and

• access functions.

It’s important to remember that in the generated code we don’t have the notion of separate components: they are listed here with the variables only in order to give the correct context to the variable names.

“Variables” are anything which does not require integration as part of the solution, and could have types COMPUTED_CONSTANT (needs to be calculated but doesn’t need integration), CONSTANT (no calculation needed), or ALGEBRAIC (TODO) as defined in the VariableType enum. They are stored in an array of VariableInfoWithType structs called VARIABLE_INFO which is VARIABLE_COUNT long.

The VariableInfoWithType contains:

• name,

• units,

• component, and

• VariableType.

2.a Get the number of variables and iterate through the VARIABLE_INFO structure to retrieve and print each variable’s information to the terminal.

Show C++ snippet

    //  2.a
    //      Get the number of variables and iterate through the VARIABLE_INFO structure to
    //      retrieve and print each variable's information to the terminal.

    std::cout << "VARIABLE_COUNT = " << VARIABLE_COUNT << std::endl;
    for (size_t v = 0; v < VARIABLE_COUNT; ++v) {
        std::cout << "Variable " << v << ": " << std::endl;
        std::cout << "  name = " << VARIABLE_INFO[v].name << std::endl;
        std::cout << "  units = " << VARIABLE_INFO[v].units << std::endl;
        std::cout << "  component = " << VARIABLE_INFO[v].component << std::endl;
        std::cout << "  type = " << VARIABLE_INFO[v].type << std::endl;
    }

Show Python snippet

    #  2.a
    #       Get the number of variables and iterate through the VARIABLE_INFO dict to
    #       retrieve and print each variable's information to the terminal.
    print('VARIABLE_COUNT = {}'.format(model.VARIABLE_COUNT))

    for v in range (0, model.VARIABLE_COUNT):
        print('Variable {}:'.format(v))
        print('  name = {}'.format(model.VARIABLE_INFO[v]['name']))
        print('  units = {}'.format(model.VARIABLE_INFO[v]['units']))
        print('  component = {}'.format(model.VARIABLE_INFO[v]['component']))
        print('  type = {}'.format(model.VARIABLE_INFO[v]['type']))
    print()

“State variables” are those which need integration. They are stored in an array of VariableInfo structs called STATE_INFO which is STATE_COUNT long.

The VariableInfo struct contains:

• name,

• units, and

• component.

2.b Get the number of state variables and iterate through the STATE_INFO structure to retrieve and print each state variable’s information to the terminal.

Show C++ snippet

    //  2.b
    //      Get the number of state variables and iterate through the STATE_INFO structure to
    //      retrieve and print each state variable's information to the terminal.
    std::cout << std::endl;
    std::cout << "STATE_COUNT = " << STATE_COUNT << std::endl;
    for (size_t s = 0; s < STATE_COUNT; ++s) {
        std::cout << "State variable " << s << ": " << std::endl;
        std::cout << "  name = " << STATE_INFO[s].name << std::endl;
        std::cout << "  units = " << STATE_INFO[s].units << std::endl;
        std::cout << "  component = " << STATE_INFO[s].component << std::endl;
    }

Show Python snippet

    #  2.b
    #       Get the number of state variables and iterate through the STATE_INFO dict to
    #       retrieve and print each state variable's information to the terminal.
    print('STATE_COUNT = {}'.format(model.STATE_COUNT))
    for s in range(0, model.STATE_COUNT):
        print('State variable {}:'.format(s))
        print('  name = {}'.format(model.STATE_INFO[s]['name']))
        print('  units = {}'.format(model.STATE_INFO[s]['units']))
        print('  component = {}'.format(model.STATE_INFO[s]['component']))
    print()

2.c Get the integration variable and print its information to the terminal. This is stored in a VariableInfo struct called VOI_INFO.

Show C++ snippet

    //  2.c 
    //      Get the integration variable and print its information to the terminal. This
    //      is stored in a VariableInfo struct called VOI_INFO.
    std::cout << std::endl;
    std::cout << "VOI_INFO" << std::endl;
    std::cout << "  name = " << VOI_INFO.name << std::endl;
    std::cout << "  units = " << VOI_INFO.units << std::endl;
    std::cout << "  component = " << VOI_INFO.component << std::endl;
    std::cout << std::endl;

Show Python snippet

    #  2.c 
    #       Get the integration variable and print its information to the terminal. This
    #       is stored in a dict called VOI_INFO.
    print('VOI_INFO')
    print('  name = {}'.format(model.VOI_INFO['name']))
    print('  units = {}'.format(model.VOI_INFO['units']))
    print('  component = {}'.format(model.VOI_INFO['component']))
    print()

Step 3: Access the functions in the generated files

C++Python

The generated code contains seven functions:

• createStatesArray() to allocate an array of length STATE_COUNT. This can be used to allocate the “rates” or gradient function array too as they’re the same length;

• createVariablesArray() to allocate an array of length VARIABLE_COUNT;

• deleteArray() to free memory used by the given array;

• initialiseStatesAndConstants(states, variables) will do what it says on the tin, and populate the given pre-allocated arrays with the initial values for all of the model’s state variables and constants.

• computeComputedConstants(variables) will fill in values for any variables that do not change in value throughout the solution, but still need to be calculated;

• computeRates(VOI, states, rates, variables) updates the rates array with the gradients of the state variables, given the values of the other variables and the variable of integration (VOI);

• computeVariables(VOI, states, rates, variables) updates any non-integrated variables whose values do not affect the integration. Since this doesn’t affect the solution process it only needs to be called whenever the values need to be output; not necessarily each integration timestep.

3.a Create three arrays representing:

• the variables (which here includes constants);

• the states (the integrated variables); and

• the rates.

Create and initialise a variable of integration, time.

Show C++ snippet

    //  3.a 
    //      Create three arrays and use the functions to allocate them, representing
    //      - variables,
    //      - rates
    //      - states.
    //      Create a variable of integration and set it to 0.
    auto myVariables = createVariablesArray();
    auto myStateVariables = createStatesArray();
    auto myRates = createStatesArray();
    double time = 0;

Show Python snippet

    #  3.a 
    #       Create three arrays representing:
    #       - the variables (which here includes constants)
    #       - the states (the integrated variables)
    #       - the rates 
    #       Create and initialise a variable of integration, time.
    my_variables = model.create_variables_array()
    my_state_variables = model.create_states_array()
    my_rates = model.create_states_array()
    time = 0.0

3.b Use the functions provided to initialise the states array you created, then print them to the screen for checking.

Show C++ snippet

    //  3.b 
    //      Use the functions provided to initialise the arrays you created, then print them 
    //      to the screen for checking.
    initialiseStatesAndConstants(myStateVariables, myVariables);

    std::cout << "The initial conditions for state variables are:" << std::endl;
    for (size_t v = 0; v < STATE_COUNT; ++v) {
        std::cout << "  " << STATE_INFO[v].component << " " << STATE_INFO[v].name << "  = " << myStateVariables[v] << " (" << STATE_INFO[v].units << ")"<< std::endl;
    }
    std::cout << std::endl;

Show Python snippet

    #  3.b 
    #       Use the functions provided to initialise the arrays you created, then print them 
    #       to the screen for checking.
    model.initialise_states_and_constants(my_state_variables, my_variables)
    print('The initial conditions for state variables are:')
    for v in range(0, model.STATE_COUNT):
        print('{} {} = {} ({})'.format(
            model.STATE_INFO[v]['component'],
            model.STATE_INFO[v]['name'],
            my_state_variables[v],
            model.STATE_INFO[v]['units']
            ))
    print()

3.c Compute the constants, compute the variables, and print them to the screen for checking.

Show C++ snippet

    //  3.c 
    //      Compute the constants, compute the variables, and print them to the screen for checking.
    computeComputedConstants(myVariables);
    computeVariables(time, myStateVariables, myRates, myVariables);
    std::cout << "The initial values including all computed constants are:" << std::endl;
    for (size_t v = 0; v < VARIABLE_COUNT; ++v) {
        std::cout << "  " << VARIABLE_INFO[v].name << " = " << myVariables[v] << " (" << VARIABLE_INFO[v].units << ")"<<std::endl;
    }
    std::cout << std::endl;

Show Python snippet

    #  3.c 
    #       Compute the constants, compute the variables,  and print them to the screen for checking.
    print('The initial values including all computed constants are:')
    model.compute_computed_constants(my_variables)
    model.compute_variables(time, my_state_variables, my_rates, my_variables)
    for v in range(0, model.VARIABLE_COUNT):
        print('  {} = {} ({})'.format(
            model.VARIABLE_INFO[v]['name'],
            my_variables[v],
            model.VARIABLE_INFO[v]['units']
            ))
    print()

Step 4: Iterate through the solution

This part will make use of a simple routine to step through the solution iterations using the Euler method to update the state variables.

4.a Create variables which control how the solution will run, representing:

• step size; and

• the number of steps to take.

Show C++ snippet

    //  4.a 
    //      Create variables which control how the solution will run, representing:
    //      - variable of integration (time);
    //      - step size; and
    //      - the number of steps to take.
    double stepSize = 0.01;
    int stepCount = 2000;
    int incr = (int)(stepCount/60) + 1;

Show Python snippet

    #  4.a 
    #       Create variables which control how the solution will run, representing:
    #       - step size and
    #       - the number of steps to take.
    step_size = 0.01
    step_count = 2000
    incr = int(step_count/60) + 1

4.b Create a file for output and open it. You can use the information to name columns with the variables, component, and units so you can keep track later.

The Euler update method is: \(x[n+1] = x[n] + x'[n].dx\)

At each step you will need to:

• Compute the rates;

• Compute the state variables using the update method above;

• Compute the variables; **

• Print to a file.

** We only need to compute these each timestep here because we’re also writing the values to the file at each timestep.

Show C++ snippet

    //  4.b 
    //      Create a file for output and open it. You can use the information to name columns
    //      with the variables, component, and units so you can keep track later.
    std::ofstream outFile("HodgkinHuxleyModelSolution.txt");
    outFile << "iteration";
    outFile << "\t" << VOI_INFO.name << " (" << VOI_INFO.units << ")";
    for (size_t v = 0; v < VARIABLE_COUNT; ++v) {
        outFile << "\t" << VARIABLE_INFO[v].component<<":"<<VARIABLE_INFO[v].name << " (" << VARIABLE_INFO[v].units << ")";
    }
    for (size_t s = 0; s < STATE_COUNT; ++s) {
        outFile << "\t" << STATE_INFO[s].component<<":"<<STATE_INFO[s].name << " (" << STATE_INFO[s].units << ")";
    }
    outFile << std::endl;
   

Show Python snippet

    #  4.b Create a file for output and open it. You can use the information to name columns
    #      with the variables, component, and units so you can keep track later.
    write_file = open('HodgkinHuxleyModelSolution.txt', 'w')
    row = 'iteration\t{}({})'.format(
        model.VOI_INFO['name'], model.VOI_INFO['units'])
    for s in range(0, model.STATE_COUNT):
        row += '\t{}({})'.format(model.STATE_INFO[s]
                                 ['name'], model.STATE_INFO[s]['units'])
    for s in range(0, model.VARIABLE_COUNT):
        row += '\t{}({})'.format(model.VARIABLE_INFO[s]
                                 ['name'], model.VARIABLE_INFO[s]['units'])
    row += '\n'
    write_file.write(row)

4.c Iterate through the time domain, calculate and write the solution at each step.

Show C++ snippet

    //  4.c
    //      Iterate through the time domain and write the solution at each step. 
    for (size_t step = 0; step < stepCount; ++step) {
        time = step * stepSize;

        // Compute the rates at this step using the given function.
        computeRates(time, myStateVariables, myRates, myVariables);
        
        // Compute the solution at the next step using whatever numerical integration 
        // method you choose. In this example we've used Euler, as given above.
        for (size_t s = 0; s < STATE_COUNT; ++s) {
            myStateVariables[s] = myStateVariables[s] + myRates[s] * stepSize;
        }

        // Compute the variables at this step using the given function.
        computeVariables(time, myStateVariables, myRates, myVariables);

        // Write everything to the output file.  Keep the order of columns consistent with
        // whatever you've used in step 4.c.
        outFile << step << "\t " << time;
        for(size_t v = 0; v < VARIABLE_COUNT; ++v) {
            outFile << "\t " << myVariables[v];
        }
        for (size_t s = 0; s < STATE_COUNT; ++s) {
            outFile << "\t" << myStateVariables[s];
        }
        outFile << "\n";
        if(step % incr == 0) {
            std::cout << "." << std::flush;
        }
    }
    std::cout << std::endl << "Finished!" << std::endl;
    outFile.close();

Show Python snippet

    #  4.c 
    #      Iterate through the time domain and write the solution at each step. 
    for step in range(0, step_count):
        time = step * step_size

        # Compute the rates at this step using the given function.
        model.compute_rates(time, my_state_variables, my_rates, my_variables)

        row = '{}\t{}'.format(step, time)

        # Compute the states.
        for s in range(0, model.STATE_COUNT):
            my_state_variables[s] = my_state_variables[s] + my_rates[s] * step_size
            row += '\t{}'.format(my_state_variables[s])

        # Compute the variables.
        model.compute_variables(time, my_state_variables, my_rates, my_variables)
        for s in range(0, model.VARIABLE_COUNT):
            row += '\t{}'.format(my_variables[s])

        row += '\n'
        # Write to file.
        write_file.write(row)

        # Print progress bar.
        if step % incr == 0:
            print('.', end='', flush=True)

    write_file.close()
   

Step 5: Sanity check

The file that results from running the completed version of this tutorial can be downloaded from HodgkinHuxleyModelSolution.txt, and is a tab-delimited file that can be easily read into your favourite plotting program.