Tutorial 3: Create a model and generate code using the API

By the time you have worked through Tutorial 3 you will be able to:

• Create a new model and its child entities from scratch using the API;

• Write your own MathML syntax to construct governing equations;

• Define custom combinations of built-in units;

• Define your own custom units independent from the built-in units; and

• Use the Generator functionality to transform the model into other languages.

This tutorial assumes that you are comfortable with:

• Accessing and adjusting names of items inside a model hierarchy (see Tutorial 2);

• Creating a validator and using it to check a model for errors (see Tutorial 2);

• Accessing the errors produced by a validator and using them to correct the model (see Tutorial 2); and

• Serialising and printing a model to a CellML2 file (see Tutorial 1).

Requirements

C++ resources

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

• tutorial3.cpp The skeleton code; or

• tutorial3_complete.cpp the completed tutorial code; and

• utilities.h and utilities.cpp Utility functions for use in the tutorials.

Python resources

• tutorial3.py Either the skeleton code; or

• tutorial3_complete.py the completed tutorial code; and

• utilities.py Utility functions for use in the tutorials.

Overview

During this tutorial you will create a simple model representing the population dynamics of two species - one a predator (sharks), and the other their prey (fish). The population of fish can only grow when they are not being constantly eaten by sharks, and the rate at which is grows will depend on how many fish are available for breeding. At the same time, the population of sharks will depend on how much food is available in the fish population. In maths this relationship can be written:

\[ \begin{align}\begin{aligned}\frac{dy_s}{dt} =f(sharks, fishes, time) = a y_s + b y_s y_f\\\frac{dy_f}{dt} =f(sharks, fishes, time) = c y_f + d y_s y_f\end{aligned}\end{align} \]

where the constants \((a, b, c, d)=(-0.8, 0.3, 1.2, -0.6)\) and we’ll use the initial conditions of \(y_s(t=0)=1.0\) and \(y_f(t=0)=2.0\).

In order to model these unusual populations you’ll need to create your own custom units, to enter and check these governing equations in MathML syntax, and to use the Generator functionality to create files able to be solved using a numerical integrator in C++ or Python.

Step 1: Set up the governing equations

Just as you did in Tutorial 2, we need to start by setting up a Model instance, and creating a Component inside it.

1.a Create a new Model and give it a name. This can be done in a single step using the overloaded constructor with a name string as its argument.

1.b Create a new Component with a name, and add it to the model you created in 1.a.

Show C++ snippet

    //  1.a
    //      Create a Model instance, set its name and id.
    auto model = libcellml::Model::create("tutorial_3_model");
    model->setId("tutorial_3_model_id");

    //  Check that it worked.
    std::cout << "Model has name: '" << model->name() << "'" << std::endl;
    std::cout << "Model has id: '" << model->id() << "'" << std::endl;

    //  1.b   
    //      Create a Component instance to use as an integrator, set its attributes and
    //      add it to the model.
    auto component = libcellml::Component::create("predator_prey_component");
    model->addComponent(component);

    //  Check that it worked.
    std::cout << "Model has " << model->componentCount()
              << " components:" << std::endl;
    for (size_t c = 0; c < model->componentCount(); ++c) {
        std::cout << "  Component [" << c << "] has name: '"
                  << model->component(c)->name() << "'" << std::endl;
        std::cout << "  Component [" << c << "] has id: '"
                  << model->component(c)->id() << "'" << std::endl;
    }

Show Python snippet

    #  1.a   
    #      Create a Model and name it.
    model = Model()
    model.setName("tutorial3_model")

    #  1.b   
    #      Create a component to use as an integrator, set its attributes and
    #      add it to the model.
    component = Component()
    component.setName("component")
    model.addComponent(component)

    #  Checking that it worked
    print_model(model)

Now for the real bit. In order to actually model anything, we need to include the mathematical equations which represent the physical situation of interest. As you saw in Tutorial 2, the maths and the Variable items which it references live inside a parent Component item.

At this point it should be noted that the order in which you add your components, or models, or variables (or anything) is not important to the final product, but it can affect how quickly you’re able to find and fix bugs along the way. In these tutorials, we have suggested that you add the mathematics first and use a Validator to notify you of the outstanding items, but you can really do this in any order you choose.

The system of equations which describe the populations are given by:

\[ \begin{align}\begin{aligned}c = a + 2.0\\\frac{dy_s}{dt} =f(sharks, fish, time) = a y_s + b y_s y_f\\\frac{dy_f}{dt} =f(sharks, fish, time) = c y_f + d y_s y_f\end{aligned}\end{align} \]

where \(y_s\) and \(y_f\) are the number of sharks and thousands of fish respectively, and the constants \((a, b, d)=(-0.8, 0.3, -0.6)\) govern their behaviour. It’s clear that the value of constant \(c\) is easily calculable from the first equation, but we will leave it in this form to better illustrate the operation of the Analyser later on.

In order to use this in our model we need to write it as a MathML2 string. The basic structure for these is described in the W3 resource pages regarding MathML2.

It’s highly unlikely that you - the user - will actually be required to write MathML2 code directly, so this part of the tutorial is more about understanding what’s going on under the hood than practising fundamentally necessary skills.

Note that libCellML will only accept MathML2 markup, even though later versions (3 and 4) are now available.

Looking at the top equation first, the MathML2 representation of \(c = a - 2.0\) is:

<apply><eq/>
   <ci>c</ci>
   <apply><plus/>
       <ci>a</ci>
       <cn>2.0</cn>
   </apply>
</apply>

Four things can be seen here:

• The <apply> opening and </apply> closing tags which surround the operations;

• The operations tags like <eq/> and <plus/> (or <minus/>, <times/>, <divide/>) which stand alone rather than in an open/close pair;

• The <ci> opening and </ci> closing tags which surround the variable names; and

• The <cn> opening and </cn> closing tags which surround the constant \(2.0\) value.

1.c Create a string containing the MathML which represents equation 1 above.

Show C++ snippet

    //  1.c
    //      Create the MathML2 strings representing the governing equations.
    std::string equation1 =
        "  <apply><eq/>\n"
        "    <ci>c</ci>\n"
        "    <apply><plus/>\n"
        "      <ci>a</ci>\n"
        "      <cn>2.0</cn>\n"
        "    </apply>\n"
        "  </apply>\n";

Show Python snippet

    #  1.c
    #      Create the MathML2 string representing the governing equations.  
    equation1 = \
        "  <apply><eq/>"\
        "    <ci>c</ci>"\
        "    <apply><plus/>"\
        "      <ci>a</ci>"\
        "      <cn>2.0</cn>"\
        "    </apply>"\
        "  </apply>"

Differential terms, such as those on the left-hand side of equations 2 and 3 \(\frac{dx}{dt}\) in MathML become:

<apply><diff/>
    <bvar>
        <ci>t</ci>
    </bvar>
    <ci>x</ci>
</apply>

Two further items to note:

• The base variable for the integration is identified by the <bvar> ... </bvar> tags. These variables are referred to as variables of integration, VOI or base variables.

• The <diff/> operation signifies differentiation with respect to the base variable.

The right-hand side becomes a collection of nested operations, all bracketed by <apply>...</apply> tags for each operation:

<apply><plus/>
  <apply><times/>
    <ci>a</ci>
    <ci>y_s</ci>
  </apply>
  <apply><times/>
    <ci>b</ci>
    <ci>y_s</ci>
    <ci>y_f</ci>
  </apply>
</apply>

When both sides are defined we need to equate them by <apply> -ing the <eq/> equals operator, and turn it into a string.

1.d Create (or copy from the snippet below) the string representing equation 2 into your code.

Show C++ snippet

    //  1.d
    std::string equation2 =
        "  <apply><eq/>\n"
        "    <apply><diff/>\n"
        "      <bvar><ci>time</ci></bvar>\n"
        "      <ci>y_s</ci>\n"
        "    </apply>\n"
        "    <apply><plus/>\n"
        "      <apply><times/>\n"
        "        <ci>a</ci>\n"
        "        <ci>y_s</ci>\n"
        "      </apply>\n"
        "      <apply><times/>\n"
        "        <ci>b</ci>\n"
        "        <ci>y_s</ci>\n"
        "        <ci>y_f</ci>\n"
        "      </apply>\n"
        "    </apply>\n"
        "  </apply>\n";

Show Python snippet

    #  1.d
    equation2 = \
        "  <apply><eq/>"\
        "    <apply><diff/>"\
        "      <bvar><ci>time</ci></bvar>"\
        "      <ci>y_s</ci>"\
        "    </apply>"\
        "    <apply><plus/>"\
        "      <apply><times/>"\
        "        <ci>a</ci>"\
        "        <ci>y_s</ci>"\
        "      </apply>"\
        "      <apply><times/>"\
        "        <ci>b</ci>"\
        "        <ci>y_s</ci>"\
        "        <ci>y_f</ci>"\
        "      </apply>"\
        "    </apply>"\
        "  </apply>"

1.e Create (or copy from the snippet below) a third string representing equation 3.

Show C++ snippet

    //  1.e
    std::string equation3 =
        "  <apply><eq/>\n"
        "    <apply><diff/>\n"
        "      <bvar><ci>time</ci></bvar>\n"
        "      <ci>y_f</ci>\n"
        "    </apply>\n"
        "    <apply><plus/>\n"
        "      <apply><times/>\n"
        "        <ci>c</ci>\n"
        "        <ci>y_f</ci>\n"
        "      </apply>\n"
        "      <apply><times/>\n"
        "        <ci>d</ci>\n"
        "        <ci>y_s</ci>\n"
        "        <ci>y_f</ci>\n"
        "      </apply>\n"
        "    </apply>\n"
        "  </apply>\n";

Show Python snippet

    #  1.e
    equation3 = \
        "  <apply><eq/>"\
        "    <apply><diff/>"\
        "      <bvar><ci>time</ci></bvar>"\
        "      <ci>y_f</ci>"\
        "    </apply>"\
        "    <apply><plus/>"\
        "      <apply><times/>"\
        "        <ci>c</ci>"\
        "        <ci>y_f</ci>"\
        "      </apply>"\
        "      <apply><times/>"\
        "        <ci>d</ci>"\
        "        <ci>y_s</ci>"\
        "        <ci>y_f</ci>"\
        "      </apply>"\
        "    </apply>"\
        "  </apply>"

Next you need to define the namespace in which the maths will be applied by enclosing it in the <math> ... </math> tags with the two namespaces:

• the MathML2 namespace xmlns at http://www.w3.org/1998/Math/MathML

• the CellML2 namespace xmlns:cellml at http://www.cellml.org/cellml/2.0#

It’s simple to do this once in your code using a string to represent the opening attributes and namespaces; this string can be reused easily throughout your code as needed later.

1.f Create the strings representing the opening and closing tags of the maths block.

Show C++ snippet

    //  1.f
    //      Create the header and footer strings.
    std::string mathHeader = "<math xmlns=\"http://www.w3.org/1998/Math/MathML\" xmlns:cellml=\"http://www.cellml.org/cellml/2.0#\">\n";
    std::string mathFooter = "</math>";

Show Python snippet

    #  1.f
    #    Add the header and footer strings.
    math_header = '<math xmlns="http://www.w3.org/1998/Math/MathML" xmlns:cellml="http://www.cellml.org/cellml/2.0#">'
    math_footer = '</math>'

Component class

• setMath

• appendMath

Utility functions (C++)

• printIssues will write information on any issues from a given Logger to the terminal. Logger classes include the Validator, Analyser, and Parser.

Utility functions (Python)

• print_issues will write information on any issues from a given Logger to the terminal. Logger classes include the Validator, Analyser, and Parser.

Our last step in defining the mathematics is to link it into the component. The functions available to manipulate maths are:

• A setMath(yourMathsHere) function, which overwrites any existing MathML strings stored in the Component item; and

• An appendMath(yourMathsHere) function, which performs a straightforward string concatenation with the current contents of the maths string in the Component.

• There is no specific function to remove maths from a component, but this can be accomplished by using the setMath function with an empty string.

1.g Use the functions above to include the MathML strings you’ve created into your component.

Show C++ snippet

    //  1.g 
    //      Add the maths strings in to the component.
    component->setMath(mathHeader);
    component->appendMath(equation1);
    component->appendMath(equation2);
    component->appendMath(equation3);
    component->appendMath(mathFooter);

Show Python snippet

    #  1.g
    #    Include the MathML strings in the component.
    component.setMath(math_header)
    component.appendMath(equation1)
    component.appendMath(equation2)
    component.appendMath(equation3)
    component.appendMath(math_footer)

Note that the order in which strings are added might be important, as they are stored as a single concatenated string. However, the order in which complete MathML equations occur in the overall MathML string is not important.

1.h Create a Validator and use it to check for errors in the model at this point. Use the utility function printIssues to output the messages to the terminal.

Show C++ snippet

    //  1.h 
    //      Create a Validator instance and use it to check for issues so far.
    //      We expect there to be 18 errors found, related to missing variables
    //      in the component.  You can use the utility printIssues function 
    //      to print them to the terminal.
    auto validator = libcellml::Validator::create();
    validator->validateModel(model);
    printIssues(validator);

Show Python snippet

    #  1.h   
    #    Create a validator and use it to check the model so far.
    validator = Validator()
    validator.validateModel(model)
    print_issues(validator)

You should see an output similar to that shown below:

Recorded 17 issues:

Issue 0 is an ERROR:
    description: MathML ci element has the child text 'c' which does not correspond with any variable names present in component 'predator_prey_component'.
    see section 2.12.3 in the CellML specification.
    more information at: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specB12.html?issue=2.12.3
    stored item type: MATH

Issue 1 is an ERROR:
    description: MathML ci element has the child text 'a' which does not correspond with any variable names present in component 'predator_prey_component'.
    see section 2.12.3 in the CellML specification.
    more information at: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specB12.html?issue=2.12.3
    stored item type: MATH

Issue 2 is an ERROR:
    description: Math cn element with the value '2.0' does not have a valid cellml:units attribute. CellML identifiers must contain one or more basic Latin alphabetic characters.
    see section 2.13.4 in the CellML specification.
    more information at: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specB13.html?issue=2.13.4
    stored item type: MATH

... etc ...

Running the validator will alert you variables in the MathML that don’t (yet) exist in your component. This was explained in Tutorial 2, and we’ll add them below.

Step 2: Create the variables

Variable class

• create

• setName

Component class

• addVariable

2.a Create Variable items for each of the missing variables discovered above. Remember that:

• Each must have a name, either using the naming constructor or by manually calling the setName function; and

• Each name must match what’s used inside your MathML string.

Show C++ snippet

    //  2.a 
    //      Create the variables listed by the validator: d, a, b, c, time, y_s, y_f.
    auto sharks = libcellml::Variable::create("y_s");
    auto fish = libcellml::Variable::create("y_f");
    auto time = libcellml::Variable::create("time");
    auto a = libcellml::Variable::create("a");
    auto b = libcellml::Variable::create("b");
    auto c = libcellml::Variable::create("c");
    auto d = libcellml::Variable::create("d");

Show Python snippet

    #  2.a
    #      Create the variables listed by the validator: d, a, b, c, time, y_s, y_f.
    sharks = Variable("y_s")
    fish = Variable("y_f")
    time = Variable("time")
    a = Variable("a")
    b = Variable("b")
    c = Variable("c")
    d = Variable("d")

2.b Add each of your new variables to the component using the addVariable function.

2.c Call the validator again to check for issues. At this stage you can expect errors like those below relating to missing units for the variables.

Show C++ snippet

    //  2.b 
    //      Add the variables into the component.
    component->addVariable(a);
    component->addVariable(b);
    component->addVariable(c);
    component->addVariable(d);
    component->addVariable(sharks);
    component->addVariable(fish);
    component->addVariable(time);

    //  2.c 
    //      Call the validator again to check.
    validator->validateModel(model);
    printIssues(validator);

Show Python snippet

    #  2.b
    #      Add the variables into the component.
    component.addVariable(sharks)
    component.addVariable(fish)
    component.addVariable(time)
    component.addVariable(a)
    component.addVariable(b)
    component.addVariable(c)
    component.addVariable(d)

    #  2.c  
    #      Call the validator again to check the model.
    validator.validateModel(model)
    print_issues(validator)

Recorded 8 issues:

Issue 0 is an ERROR:
    description: Variable 'a' in component 'predator_prey_component' does not have any units specified.
    see section 2.8.1.2 in the CellML specification.
    more information at: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specB08.html?issue=2.8.1.2
    stored item type: VARIABLE

... etc ...

Step 3: Built-in and customised units

Linking variables to the name of their units is straightforward, but in order to be able to use them we need to also define what the name actually means by creating the units themselves. Some common units have been defined and built into libCellML, others you can define by combining the built-in ones using scaling factors and exponents, or you can define your own from scratch if need be. You can read more about units on the Understanding units page.

Units class

• create

• setName

• addUnit (note singular)

Model class

• addUnits (note plural)

• linkUnits

Variable class

• setUnits (note plural)

To create a Units item you need will follow the same basic steps as other entities: declare it, name it, define it, and then add it in. For example:

C++Python

// Declare, name, and define a "millisecond" unit pointer.
auto ms = libcellml::Units::create("millisecond");

// The manner of specification here is agnostic: all three definitions are identical.
ms->addUnit("second", "milli");  reference unit and built-in prefix
// OR
ms->addUnit("second", 1.0, -3);  reference unit, prefix, exponent
// OR
ms->addUnit("second", 1.0, 0, 0.001);  reference unit, prefix, exponent, multiplier

3.a Use the example above to create, name and define the units of “month” which will represent your time variable. This should be defined as a multiple of the built-in unit second.

Show C++ snippet

    //  3.a 
    //      Create units representing a month, or 2592000 seconds.
    auto month = libcellml::Units::create("month");
    month->addUnit("second", 0, 1, 2592000);  // base unit, prefix, exponent, multiplier
    model->addUnits(month);

Show Python snippet

    #  3.a  
    #      Define the relationship between our custom units and the built-in
    #      units. There is a list of built-in units and their definitions
    #      available in section 19.2 of the CellML2 specification.
    #      First we create the "month" units, which will be equivalent to
    #      60*60*24*30 = 2,592,000 seconds.
    month = Units("month")
    month.addUnit("second", 1, 2592000)  # Setting a month to be 2592000 seconds.
    model.addUnits(month)

    #  "second" is a built-in unit, used with a multiplier of 2592000.
    #  Note that this could have been written:
    #    month.addUnit("second", "mega", 1, 2.592)
    #    month.addUnit("second", 5, 25.92)

Units can be defined based on one another as well. For example, after defining our millisecond units, we could then use this definition to define the per_millisecond units by simply including it with an exponent of -1:

Define a per_millisecond unit based on millisecond^-1:
per_ms->addUnit(ms, -1.0);

3.b Create a Units item called “per_month” based on the one you just created, as shown above.

Show C++ snippet

    //  3.b 
    //      Create the per_month unit based on the month defined in 3.a.
    auto per_month = libcellml::Units::create("per_month");
    per_month->addUnit("month", -1);  // base unit, exponent
    model->addUnits(per_month);

Show Python snippet

    #  3.b  
    #      Create units which represent "per_month", which
    #      is simply the inverse of the "month" unit above.
    per_month = Units()
    per_month.setName("per_month")
    per_month.addUnit("month", -1)
    model.addUnits(per_month)

3.c Create the irreducible units needed by the shark and fish populations. Call these “number_of_sharks” and “thousands_of_fish” respectively.

Show C++ snippet

    //  3.c 
    //      Create the sharks and fishes base units, "number_of_sharks" and "thousands_of_fish".
    auto number_of_sharks = libcellml::Units::create("number_of_sharks");
    auto thousands_of_fish = libcellml::Units::create("thousands_of_fish");
    model->addUnits(number_of_sharks);
    model->addUnits(thousands_of_fish);

Show Python snippet

    #  3.c      
    #      Create the sharks and fishes base units.
    number_of_sharks = Units()
    number_of_sharks.setName("number_of_sharks")
    model.addUnits(number_of_sharks)
    thousands_of_fish = Units()
    thousands_of_fish.setName("thousands_of_fish")
    model.addUnits(thousands_of_fish)

Finally we need to create the units for the constants b and d. These will be combinations of those which we’ve already created, as defined by the need for dimensional consistency in our governing equations.

3.d Create two units representing “per shark month” (for the b variable) and “per fish month” (for the d variable).

Show C++ snippet

    //  3.d 
    //      Create the combined units for the constants, "per_shark_month" and "per_fish_month".
    auto b_units = libcellml::Units::create("per_shark_month");
    b_units->addUnit("per_month");
    b_units->addUnit("number_of_sharks", -1);
    model->addUnits(b_units);

    auto d_units = libcellml::Units::create("per_1000fish_month");
    d_units->addUnit("per_month");
    d_units->addUnit("thousands_of_fish", -1);
    model->addUnits(d_units);

Show Python snippet

    #  3.d  
    #      Create the combined units for the constants.  Note that each item included
    #      with the addUnit command is multiplied to create the final Units definition.
    b_units = Units()
    b_units.setName("per_shark_month")
    b_units.addUnit("per_month")
    b_units.addUnit("number_of_sharks", -1)
    model.addUnits(b_units)

    d_units = Units()
    d_units.setName("per_fish_month")
    d_units.addUnit("per_month")
    d_units.addUnit("thousands_of_fish", -1)
    model.addUnits(d_units)

The final two steps are to associate each variable with its appropriate units, and to include the units in the model.

• When you add different sub-unit parts into a Units item, the function is addUnit (singular), and it takes as argument the name of the sub-unit as a string (eg: "second" used above).

• When you add the final created combination into the Model item, the function is addUnits (plural), and it takes as argument the reference of the combined units (eg: ms).

3.e Add the units to their variables using myVariable->setUnits(myUnits). Add the units to the model using myModel->addUnits(myUnits).

Show C++ snippet

    //  3.e 
    //      Add the units to their variables using the setUnits function.
    time->setUnits(month);
    a->setUnits(per_month);
    b->setUnits(b_units);
    c->setUnits(per_month);
    d->setUnits(d_units);
    sharks->setUnits(number_of_sharks);
    fish->setUnits(thousands_of_fish);

Show Python snippet

    #  3.e  
    #      Set the units to their respective variables.
    time.setUnits(month)
    sharks.setUnits(number_of_sharks)
    fish.setUnits(thousands_of_fish)
    a.setUnits(per_month)
    b.setUnits(b_units)
    c.setUnits(per_month)
    d.setUnits(d_units)

Gotcha! When you specify the Units for a Variable using its name then you may need to call the model’s linkUnits function before validating the model. If you see errors related to missing units which do in fact exist, this indicates that a call to the linkUnits function is needed.

3.f Call the validator to check your model for errors. You should see an output similar to that shown below.

Show C++ snippet

    //  3.f 
    //      Call the validator to check the model.  We expect one error regarding the missing units in the MathML.
    validator->validateModel(model);
    printIssues(validator);

Show Python snippet

    #  3.f  
    #      Call the validator again to check the model.
    #      Expect one error regarding a missing unit in the MathML.
    validator.validateModel(model)
    print_issues(validator)

Recorded 1 issues:

Issue 0 is an ERROR:
    description: Math cn element with the value '2.0' does not have a valid cellml:units attribute. CellML identifiers must contain one or more basic Latin alphabetic characters.
    see section 2.13.4 in the CellML specification.
    more information at: https://cellml-specification.readthedocs.io/en/latest/reference/formal_and_informative/specB13.html?issue=2.13.4
    stored item type: MATH

In the first MathML equation we used a real number <cn>2.0</cn> without specifying any units for it. Because the dimensionality of the equation needs to be valid, all real numbers must be associated with units, just the same way that variables are. These are defined within the tags of the MathML, and must also refer to the cellml namespace. For example:

<cn cellml:units="bunch_of_bananas">1</cn>

… which gives us one bunch of bananas, without needing to create a corresponding Variable item. Of course, you may need to create the corresponding Units item and add it to the model, if it is not already present.

3.g Create a copy of the MathML statement from step 1.c and add the namespace and units definition as in the example above into the string. Recall that using the setMath() function will overwrite the existing maths, and repeat the process you did in step 1.e to include the new MathML instead. Remember that you will need to re-include the opening and closing <math> tags as well as the other equations.

3.h Check that the model is now free of validation errors.

Show C++ snippet

    //  3.g 
    //      Units for constants inside the MathML must be specified at the time.  This means we need to adjust
    //      equation1 to include the per_month units.  We have to wipe all the existing MathML and replace it.
    component->removeMath();
    component->setMath(mathHeader);
    equation1 =
        "  <apply><eq/>\n"
        "    <ci>c</ci>\n"
        "    <apply><plus/>\n"
        "      <ci>a</ci>\n"
        "      <cn cellml:units=\"per_month\">2.0</cn>\n"
        "    </apply>\n"
        "  </apply>\n";

    component->appendMath(equation1);
    component->appendMath(equation2);
    component->appendMath(equation3);
    component->appendMath(mathFooter);

    //  3.h 
    //      Revalidate your model and expect there to be no errors.
    validator->validateModel(model);
    printIssues(validator);
    assert(validator->errorCount() == 0);

Show Python snippet

    #  3.g  
    #      Units for constants inside the MathML must be specified at the time.
    #      This means we need to adjust equation1 to include the per_month units.
    #      We have to wipe all the existing MathML and replace it.
    component.removeMath()
    component.setMath(math_header)
    equation1 = \
        "  <apply><eq/>"\
        "    <ci>c</ci>"\
        "    <apply><plus/>"\
        "      <ci>a</ci>"\
        "      <cn cellml:units=\"per_month\">2.0</cn>"\
        "    </apply>"\
        "  </apply>"
    component.appendMath(equation1)
    component.appendMath(equation2)
    component.appendMath(equation3)
    component.appendMath(math_footer)

    #  3.h  
    #      Validate once more, and expect there to be no errors this time.
    validator.validateModel(model)
    print_issues(validator)

Step 4: Analyse the mathematical model

The Analyser class checks that the underlying mathematical model represented by the entire combination of variables, components, and mathematics, makes sense. The Validator checks your “spelling” and syntax, and the Analyser checks for the ability to find a solution.

Analyser class

• create

• analyseModel

• model

4.a Create an Analyser instance and pass it the model using the analyseModel function.

4.b Check for issues found in the analyser. You should expect 6 errors, related to variables whose values are not computed or initialised. Note that you can use the same utility function printIssues to output issues from the analyser as from the validator.

Show C++ snippet

    //  4.a 
    //      Create an Analyser instance and pass it the model using the
    //      analyseModel function.  
    auto analyser = libcellml::Analyser::create();
    analyser->analyseModel(model);

    //  4.b 
    //      Check for errors found in the analyser. You should expect 6 errors,
    //      related to variables whose values are not computed or initialised.
    printIssues(analyser);

Show Python snippet

    #  4.a 
    #     Create an Analyser instance and pass it the model using the
    #     analyseModel function.  
    analyser = Analyser()
    analyser.analyseModel(model)

    #  4.b 
    #     Check for errors found in the analyser. You should expect 6 errors,
    #     related to variables whose values are not computed or initialised.
    print_issues(analyser)

The messages above refer to the fact that though our model has passed validation tests, it’s not yet sufficiently constrained to allow it to be solved, which is what the Generator checks for. We need to set initial conditions for the variables we’re solving for, the populations of sharks and fish, using the setInitialValue function. The values of the constants a, b, c, d are just that - constant - and their values are set using the same setInitialValue function.

4.c Set the values of the constants \((a, b, d)=(-0.8, 0.3, -0.6)\) and the initial conditions such that \(y_f(t=0)=2.0\) and \(y_s(t=0)=1.0\). Note that:

• The constant \(c\) will be calculated by our equation 1, so does not need to be specified; and

• The base variable (or “variable of integration”, or “voi”) \(t\) must not have an initial condition set.

4.d Reprocess the model and check that the analyser is now free of issues.

Show C++ snippet

    //  4.c 
    //      Add initial conditions to all variables except the base variable, time
    //      and the constant c which will be computed. Reprocess the model.
    a->setInitialValue(-0.8);
    b->setInitialValue(0.3);
    d->setInitialValue(-0.6);
    sharks->setInitialValue(1.0);
    fish->setInitialValue(2.0);

    //  4.d 
    //      Reprocess the model and check that the generator is now free of errors.
    analyser->analyseModel(model);
    printIssues(analyser);

Show Python snippet

    #  4.c 
    #     Add initial conditions to all variables except the base variable, time
    #     and the constant c which will be computed. Reprocess the model.
    a.setInitialValue(-0.8)
    b.setInitialValue(0.3)
    d.setInitialValue(-0.6)
    sharks.setInitialValue(1.0)
    fish.setInitialValue(2.0)

    #  4.d 
    #     Reprocess the model and check that the analyser is now free of errors.
    analyser.analyseModel(model)
    print_issues(analyser)

Step 5: Generate code and output

Some exciting new functionality in libCellML is the ability to generate a runnable file from a model description. This means that if you already have a solver in either C or Python, you can simply translate your model from here into that language.

The Generator has to re-interpret all of the maths, including the variables, their interaction with each other in different equations, values, initial conditions and units before it can output your model in your choice of language. For the maths to make sense, the definitions in your model’s variables, maths blocks and units need to be solvable too. There are several requirements that need to be satisfied in order for the code generation functionality to be able to work, beyond the CellML syntax requirements. These are:

• The mathematical model definition must be appropriately constrained (not over- or under-constrained);

• Initial conditions must be specified for variables which are integrated;

• Initial conditions must not be specified for variables which are the base of integration;

• The values of constants must be specified or calculable; and

• TODO get full list of stuff here …

Generator class

• create

• profile

• setModel

• interfaceCode

• implementationCode

GeneratorProfile class

• create

• setInterfaceFileNameString

The GeneratorProfile class contains an enum indicating the language of profile to set. In C++ this is GeneratorProfile::Profile. In Python this is GeneratorProfile.Profile.

At the time of writing two profiles are available:

• C (default)

• PYTHON

5.a Create a Generator instance. Instead of giving it the model to process, the generator needs an analysed model as its argument. Retrieve the analysed model using the analyser’s model function and pass it to the generator using the setModel function.

Show C++ snippet

    //  5.a  
    //      Create a Generator instance.  Instead of giving it the Model item to process, 
    //      the generator takes the output from the analyser.  
    //      Retrieve the analysed model using the Analyser::model() function and pass it
    //      to the generator using the Generator::setModel function.
    auto generator = libcellml::Generator::create();
    generator->setModel(analyser->model());

    //  The generator takes the CellML model and turns it into procedural code in another 
    //  language.  The default is C, but Python is available too.  This language choice is
    //  called the "profile", and is stored in a GeneratorProfile item.

    //  If you're using the C profile then you have the option at this stage 
    //  to specify the file name of the interface file you'll create in the
    //  next step.  This means that the two files will be prepared to link to
    //  one another without manual editing later.  

Show Python snippet

    #  5.a  
    #     Create a Generator instance.  Instead of giving it the Model item to process, 
    #     the generator takes the output from the analyser.  
    #     Retrieve the analysed model using the Analyser.model() function and pass it
    #     to the generator using the Generator.setModel function.
    generator = Generator()
    generator.setModel(analyser.model())

    #  The generator takes the CellML model and turns it into procedural code in another 
    #  language.  The default is C, but Python is available too.  This language choice is
    #  called the "profile", and is stored in a GeneratorProfile item.

    #  The default profile has a C flavour, and it already exists inside the Generator you've just created.
    #  We need to edit that profile a little, but only to tell it the file name where they interface
    #  (header file) code will be written.  This is so that the implementation code (source file) knows
    #  where to look when it tries to #include it.  

The generator takes the CellML model and turns it into procedural code in another language. The default is C, but Python is available too. This language choice is called the “profile”, and is stored in a GeneratorProfile item.

The default profile already exists inside the Generator you’ve just created. We need to edit that profile a little, but only to tell it the file name where they interface (header file) code will be written. This is so that the implementation code (source file) knows where to look when it tries to include it.

5.b Retrieve the C profile from the generator, and use its setInterfaceFileNameString function to pass in the same filename that you’ll use in 5.c below for the interface code.

Show C++ snippet

    //  5.b
    //      You can do this by specifying the header file name in the GeneratorProfile item
    //      using the setInterfaceFileNameString("yourHeaderFileNameHere.h") function.
    //      This will need to be the same as the file which you write to in step 5.c below.
    auto profileC = generator->profile();
    profileC->setInterfaceFileNameString("PredatorPrey.h");

Show Python snippet

    #  5.b
    #     Retrieve the C profile from the generator, and use its setInterfaceFileNameString function
    #     to pass in the same filename that you'll use in 5.c for the interface code.
    generator.profile().setInterfaceFileNameString('PredatorPrey.h')

5.c Since we’re using the default profile (C), we need to output both the interface code (the header file) and the implementation code (the source file) from the generator and write them to their respective files.

Show C++ snippet

    //  5.c 
    //      First we'll use the default profile (C), so we need to output both the
    //      interfaceCode (the header file) and the implementationCode (source file)
    //      from the generator and write them to their respective files.
    std::ofstream outFile("PredatorPrey.h");
    outFile << generator->interfaceCode();
    outFile.close();

    outFile.open("PredatorPrey.c");
    outFile << generator->implementationCode();
    outFile.close();

Show Python snippet

    #  5.c
    #     First we'll use the default profile (C), so we need to output both the
    #     interfaceCode (the header file) and the implementationCode (source file)
    #     from the generator and write them to their respective files.
    implementation_code_C = generator.implementationCode()
    write_file = open("PredatorPrey.c", "w")
    write_file.write(implementation_code_C)
    write_file.close()

    interface_code = generator.interfaceCode()
    write_file = open("PredatorPrey.h", "w")
    write_file.write(interface_code)
    write_file.close()

The generator takes the CellML model and turns it into procedural code in another language. The default is C, but Python is available too. This language choice is called the “profile”, and is stored in a GeneratorProfile item.

5.d Create a GeneratorProfile item using the Profile::PYTHON enum value in the constructor. Pass this profile to the setProfile function in the generator.

Show C++ snippet

    //  5.d
    //      Create a GeneratorProfile item using the libcellml::GeneratorProfile::Profile::PYTHON
    //      enum value in the constructor.  Pass this profile to the setProfile function in the
    //      generator.
    auto profilePython =
        libcellml::GeneratorProfile::create(libcellml::GeneratorProfile::Profile::PYTHON);
    generator->setProfile(profilePython);

Show Python snippet

    #  5.d
    #     Create a GeneratorProfile item using the libcellml.GeneratorProfile.Profile.PYTHON
    #     enum value in the constructor.  Pass this profile to the setProfile function in the
    #     generator.
    profile = GeneratorProfile(GeneratorProfile.Profile.PYTHON)
    generator.setProfile(profile)

5.e Retrieve the Python implementation code (there is no header file) and write to a *.py file.

Show C++ snippet

    //  5.e
    //      Retrieve the Python implementation code (there is no header file) and write to a *.py file.
    outFile.open("PredatorPrey.py");
    outFile << generator->implementationCode();
    outFile.close();

Show Python snippet

    #  5.d
    #     Create a GeneratorProfile item using the libcellml.GeneratorProfile.Profile.PYTHON
    #     enum value in the constructor.  Pass this profile to the setProfile function in the
    #     generator.
    profile = GeneratorProfile(GeneratorProfile.Profile.PYTHON)
    generator.setProfile(profile)

    #  5.e
    #     Retrieve the Python implementation code (there is no header file) and write to a *.py file.
    implementation_code_python = generator.implementationCode()
    write_file = open("PredatorPrey.py", "w")
    write_file.write(implementation_code_python)

Go and have a cuppa, you’re done!