Chapter 13 Step-by-step tutorial 𝑰: Re-creating the TBT module

13.1 Create the module

New modules are created in the maven module EvoLudoCore. In order to get started, create a new java file org.evoludo.simulator.modules.TxT in the EvoLudoCore maven module. In VS Code this is done by right-clicking on the modules folder in EvoLudoCore/src/main/java/org/evoludo/simulator and selecting New Java File ¿ Class…. Enter the name of the new file as TxT. The new file should look like this:

1 package org.evoludo.simulator.modules;

2

3 public class TxT {

4

5 }

Note, that the following code snippets in this tutorial contain little to no comments because the TBT class serves as a reference file. Besides, the motivation and purpose of each code snippet is explained in detail in the accompanying text.

The new TxT class must extend the Discrete class because we are dealing with a discrete set of traits, say $A$ and $B$ . This is done by changing the class definition from TxT to TxT extends Discrete. The Discrete class provides a number of useful methods and fields to handle discrete traits (and their interactions) as opposed to continuous traits. The Discrete class is located in the same package as the TxT class, so no import statement is needed.

However, extending Discrete triggers two errors in VS Code because (i) the Discrete class does not provide a default constructor and (ii) TxT does not implement the abstract method getTitle() of the Discrete class.

.

Note: The Discrete class has no default constructor and therefore subclasses must provide their own constructor using one of the constructors of the super class. All single species modules use the constructor Discrete(EvoLudo engine).

.

Both errors are easily addressed in VS Code by hovering over the error location and selecting Quick fix…, followed by Add constructor TxT(EvoLudo) as well as Add unimplemented methods. Now the content of the file looks as follows:

1 package org.evoludo.simulator.modules;

2

3 import org.evoludo.simulator.EvoLudo;

4

5 public class TxT extends Discrete {

6

7 public TxT(EvoLudo engine) {

8 super(engine);

9 //TODO Auto-generated constructor stub

10 }

11

12 @Override

13 public String getTitle() {

14 // TODO Auto-generated method stub

15 throw new UnsupportedOperationException("Unimplemented method ’getTitle’");

16 }

17 }

In the constructor there is nothing else we need to do, so simply delete the TODO comment.

.

Note: If you follow the steps in this tutorial the line numbers in the code snippets should always match the line numbers in your file.

.

The method getTitle() returns a short title of the module:

11 @Override

12 public String getTitle() {

13 return "2x2 Games Tutorial";

14 }

This information will be displayed in the list of available modules on the EvoLudo help screen. Before we can verify this we need to make EvoLudo aware of the new module.

13.1.1 Adding the module to EvoLudo

In order to run the new module in EvoLudo, we first need to add it to the list of available modules. This is done in the EvoLudo class: open org.evoludo.simulator.EvoLudo in the EvoLudoCore maven module and located in the parent directory of our TxT class. The EvoLudo(boolean loadModules) constructor at the end of the file includes several addModule(...) calls to add all available modules. Add the following line to the list of modules:

addModule(new TxT(this));

together with the necessary import statement at the top:

import org.evoludo.simulator.modules.TxT;

13.2 Running (and debugging) the module

Open the terminal in VS Code. Make sure the current working directory is the root directory of the EvoLudo project. Compile the project with the command ’mvn clean install’. This may take a few moments. Next, run ’mvn gwt:devmode’. This starts the GWT development web server.

Now switch to the Run and Debug view in VS Code and launch the GWT EvoLudo (localhost) target. If everything is properly configured, this opens a window in the Google Chrome browser and displays the message Compiling evoludoweb. This shouldn’t take long and then display the familiar EvoLudo GUI.

In the GUI click on Settings and then Help to see the list of available modules. The newly minted TxT module is listed as: TxT: 2x2 Games Tutorial. The first part, TxT, is the key to load the module and the second part the descriptive string returned by getTitle().

.

Note: The key to load a module can be customized by overriding the method getKey(). The key defaults to the class name of the module, here TxT.

.

The module is loaded by entering --module TxT in the settings field and then clicking Apply. Loading the module merely adds a line to the console: Module loaded: 2x2 Games Tutorial and reports an error that no model was found:

EvoLudo: v2.0.2 © Christoph Hauert

GWT Version: 2.11.0

GUI features: WebGL XML keyboard mouse

Module loaded: 2x2 Games Tutorial

ERROR: No model found!

.

Note: Similarly, the method getAuthors() can be overridden to claim the authorship of the module. By default, no authors are returned.

.

More importantly, note that running ’mvn gwt:devmode’ actually set up a powerful development environment: after saving changes to the source code all that is needed to recompile and run the revised code is to reload the browser window! However, reloading the browser window also resets EvoLudo and, more specifically, the option string in the settings field.

Therefore, to make our lives easier by not having to adjust the settings all the time, it’s a good idea to edit the file EvoLudoGWT/src/main/webapp/TestEvoLudo.html and look for the string <div class="evoludo-simulation" $\ldots$ > and change its data-clo attribute to --module TxT. This way the TxT module is loaded by default when the browser window is reloaded.

.

Important: The development server that was started with ’mvn gwt:devmode’ serves pages out of the EvoLudoGWT/target/EvoLudoGWT-1.0-SNAPSHOT folder. Hence, changes to EvoLudoGWT/src/main/webapp/TestEvoLudo.html are ignored until the project is recompiled with ’mvn clean install’. Alternatively, the file EvoLudoGWT/target/EvoLudoGWT-1.0-SNAPSHOT/TestEvoLudo.html can be edited directly for immediate effect, but those changes are eventually overridden by ’mvn clean install’.

.

Moreover, the development environment allows setting breakpoints in the Java code of VS Code and debug the JavaScript code running in the browser. How cool is that?!

.

Important: Step-by-step debugging is possible by setting break points in VS Code.

.

While debugging becomes more important later on and when adding your own customizations, it doesn’t hurt to give it a try now. For example, set a break point in VS Code at the return of the method getTitle() (line $13$ in the TxT class) by clicking the line number on the left side of the editor. The break point is marked by a $\bullet$ . Clicking again on the line number removes the break point.

The method getTitle() is called when displaying the list of available options. Click on Help again. This will trigger a call to the method getTitle() and the execution of the code stops at the requested line. The Debug view in VS Code shows the current stack trace and all local variables.

Now you can explore the values of all variables and step through the code. Once you have convinced yourself of how seamlessly Google Chrome and VS Code work together, you can remove the breakpoint by clicking on the $\bullet$ . In order to resume execution of the EvoLudo code, simply click on the Continue, F5 button in the debug view.

13.3 Implementing the module

Many features and capabilities of EvoLudo modules are advertised by implementing certain interfaces. More specifically, with our TxT module we want to model pairwise interactions in $2\times 2$ games, which means individuals not only carry a trait but also a score (or fitness) that is determined through payoffs from interactions with other members of the population. In order to advertise that our module uses payoffs we need to implement the interface Payoffs:

6 public class TxT extends Discrete implements Payoffs {

Now VS Code alerts us that the interface Payoffs requires the implementation of two methods: getMinGameScore() and getMaxGameScore(). As before, just hover over the location of the error (marked by a red squiggly line) and select Quick fix… followed by Add unimplemented methods.

The two methods, getMinGameScore() and getMaxGameScore(), must return the minimum and maximum scores, respectively, that an individual can achieve in a single interaction. Those bounds on the scores are important to convert payoffs to probabilities and properly scale the range of views in the GUI. We will get back to those methods later. That is, once we have implemented the payoff matrix.

In order to model $2\times 2$ games we need to define the payoff matrix as well as the indices of the different traits. For readability, we define two constants for the trait indices. In TBT the two traits are called COOPERATE and DEFECT for historical reasons. Let’s call them A with index 0 and B with index 1. Next, we define a two-dimensional array payoffs, where payoffs[x][y] holds the payoff to an individual with trait x against another with trait y, where x and y refer to either A or B.

After adding the definition of constants and the payoff array to the beginning of our module, we get:

6 public class TxT extends Discrete implements Payoffs {

7

8 public static final int A = 0;

9

10 public static final int B = 1;

11

12 double[][] payoffs;

Note, to conserve memory, do not allocate the payoffs array. The gain in this case is minimal, of course, but it is still a good habit. Also, you are encouraged to add comments to explain the purpose of all variables and methods. Now, before implementing the two remaining, auto-generated method stubs, let’s make the newly created module available in EvoLudo.

13.3.1 Loading and unloading of the module

The EvoLudo code has a highly modular structure with modules for different evolutionary setups, models for different kinds of analysis and modes for analyzing the evolutionary process (see section 11.2.1: Modules, Models and Modes). The new module TxT is loaded into the EvoLudo engine with the setting --module TxT, where TxT denotes the key to identify a module. The key defaults to the class name but can be customized by overriding the optional method getKey(). Whenever a module gets loaded, the method load() is called. If any initialization or memory allocations are required, this is the place to do so. In our case, this means allocating (and freeing) the memory for the payoff matrix. In addition, the names and colours for the different traits can be set, but this is optional. Before a new module is loaded the current one (if any) gets unloaded by calling the method unload() to free up resources.

.

Important: The number of traits nTraits must be either set during loading by overriding the method load() (and include a call super.load()), or in a static block of the class.

.

23 @Override

24 public void load() {

25 super.load();

26 nTraits = 2;

27 // trait names (optional)

28 String[] names = new String[nTraits];

29 names[A] = "A";

30 names[B] = "B";

31 setTraitNames(names);

32 // trait colors (optional)

33 Color[] colors = new Color[2 * nTraits];

34 colors[A] = Color.BLUE;

35 colors[B] = Color.RED;

36 colors[A + nTraits] = Color.GREEN;

37 colors[B + nTraits] = Color.YELLOW;

38 setTraitColors(colors);

39 // allocate memory for payoffs

40 payoffs = new double[nTraits][nTraits];

41 }

The Color class requires another import statement at the top:

import java.awt.Color;

.

Note: The colors-array can hold nTraits or 2 * nTraits colours. The first nTraits colours are used to mark the traits of individuals in the GUI views and tooltips (or the output when exporting graphics). The second set of nTraits colours mark the traits of individuals that have changed since the previous update. If no (or insufficient numbers of) colours are set, the default colours of Color.BLUE, Color.RED, Color.GREEN, Color.YELLOW, Color.MAGENTA, Color.ORANGE, Color.PINK, Color.CYAN are used for up to eight traits. Further trait colours are generated at random. Newly changed traits are marked by a paler shade of the trait colour.

.

Note: The names-array holds the names for all nTraits traits. If no names are set, the default names are "Trait A", "Trait B" etc. This affects the axis labels and tooltips in the GUI as well as the output when exporting data.

.

Similarly, when unloading the module (before potentially loading a different one), the method unload() is called. This is primarily an opportunity to release all resources that were allocated in load().

45 @Override

46 public void unload() {

47 super.unload();

48 payoffs = null;

49 }

.

Important: Always include a call to super.unload() to release all resources allocated by parent classes. In most cases the garbage collector is unable to identify unused resources because the shell of every module is initialized at launch time and remains in memory.

.

13.4 Module specific fields

The TxT module introduces only a single custom field payoffs to store the $2\times 2$ game matrix. Restricting access to fields of an object is not only a good habit but failure to do so would defeat the basic tenets of object-oriented programming. Therefore, we provide setter and getter methods to access the payoff matrix from outside the module:

63 public void setPayoffs(double[][] payoffs) {

64 this.payoffs = payoffs;

65 }

66

67 public double[][] getPayoffs() {

68 return payoffs;

69 }

.

Note: Another good habit is to include Javadoc descriptions of all fields, constants and especially methods.

.

Since we can always refer to the corresponding documentation in the original TBT module, we will not include the Javadoc descriptions in this tutorial.

For quick and dirty tests it is possible to hard code parameter settings. For example, we could fix the payoff matrix to Axelrod’s famous payoffs by replacing the line

42 payoffs = new double[nTraits][nTraits];

with

42 payoffs = new double[][]{{3, 0}, {5, 1}};

at the end of the method load(). This is perfectly fine during code development and testing. However, in the long run a much better and more flexible approach is to add with very little extra effort a new command line option to set the payoff matrix.

13.4.1 Adding command-line options

All command line options are instances of the CLOption class. The class provides a number of constructors for different kinds of command line options.

Using the constructor CLOption(…)

The most common constructor is: CLOption(String name, String defaultArg, Category category, String description, CLODelegate delegate) and takes five arguments:

name:: the name of the command line option (without the ’--’ prefix): "paymatrix".
defaultArg:: the default value of the command line option (Axelrod’s famous payoffs in the prisoner’s dilemma): "3,0;5,1".
category:: the command line option falls into the module category (for organizing options in the help screen): Category.Module.
description:: a short description of the command line option (for the help screen): "--paymatrix <a,b;c,d> 2x2 payoff matrix".
delegate:: the delegate to parse the command line option and configure the module accordingly.

As the last argument the CLOption(...) constructor expects an implementation of the CLODelegate interface.

Creating the delegate

The delegate must implement the CLODelegate interface and provide the method parse(String arg) to process the command line option and its arguments. The delegate defines three methods, but the implementation of only one is required:

boolean parse(String arg):: (required) parse the arg string and set the parameters of the module accordingly. If parsing is successful, the payoff matrix is set and the method returns true. If parsing fails, for example if the matrix has wrong dimensions or non-numerical entries, the method must return false. This triggers a warning message with the name and argument of the option that caused trouble, which is logged (and shown in the GUI). Finally, if the option is missing, or parsing failed, the method parse(String) is called with the defaultArg provided in the CLOption constructor (see above).
String getDescription():: (optional) return the description of the option. In most cases the String description argument to the CLOption(...) constructor is sufficient. Implementing the method getDescription() is only required if the description depends on other parameter settings. Most commonly this is the number of traits in a module, which may affect the formatting of the argument arg.
int getKeyPos():: (optional) some options accept a selection of key values. For example, the options --module or --model. The getKeyPos() method returns the position of the key in the argument string for options that accept keys and additional arguments. By default, the key is expected in position 0 and may be followed by additional parameters. The option --init is an example. In contrast, --mutation, for example, expects the key in position 1.

In Java this CLODelegate can be easily implemented through anonymous inner classes. Instead of writing a custom class that implements the CLODelegate interface, we can simply pass the following code as the last argument to the above CLOption(...) constructor:

new CLODelegate() {

@Override

public boolean parse(String arg) {

// parse arg and set parameters

return true;

}

This creates an anonymous inner class (one without a name) that implements the CLODelegate interface and overrides the method parse(String arg). Assembling all pieces yields the following code snippet to define the command line option --paymatrix:

74 public final CLOption cloPayoffs = new CLOption("paymatrix", "3,0;5,1", Category.Module,

75 "--paymatrix <a,b;c,d> 2x2 payoff matrix", new CLODelegate() {

76

77 @Override

78 public boolean parse(String arg) {

79 double[][] payMatrix = CLOParser.parseMatrix(arg);

80 if (payMatrix == null || payMatrix.length != 2

81 || payMatrix[0].length != 2)

82 return false;

83 setPayoffs(payMatrix);

84 return true;

85 }

86 });

The code snippet above requires importing the following classes:

7 import org.evoludo.util.CLOParser;

8 import org.evoludo.util.CLOption;

9 import org.evoludo.util.CLOption.CLODelegate;

and has already been taken into account with the line numbers.

.

Important: Make sure the default arguments are valid. Invalid default arguments may result in unexpected or undefined behaviour.

.

Finally, register the new command line option with the option parser:

88 @Override

89 public void collectCLO(CLOParser parser) {

90 super.collectCLO(parser);

91 parser.addCLO(cloPayoffs);

92 }

13.4.2 Tying up loose ends

After dealing with the payoff matrix, we can finally implement the methods getMinGameScore() and getMaxGameScore(). The minimum score is the smallest entry in the payoff matrix and the maximum score is the largest entry. We can use the ArrayMath utility class to find the minimum and maximum values in the payoff matrix. Moreover, we can optionally override the method getMonoGameScore(int type) to return the payoffs in a monomorphic population of type (e.g., A or B), which corresponds to the diagonal entries in the matrix payoffs:

54 @Override

55 public double getMinPayoff() {

56 return ArrayMath.min(payoffs);

57 }

58

59 @Override

60 public double getMaxPayoff() {

61 return ArrayMath.max(payoffs);

62 }

63

64 @Override

65 public double getMonoPayoff(int type) {

66 return payoffs[type][type];

67 }

together with the necessary import statement

import org.evoludo.math.ArrayMath;

The minimum and maximum scores are essential for converting fitness to probabilities in individual based simulations, and to rates in differential equation models. In addition, the range of scores is used to set the scale of fitness in several data views as well as to scale the colour gradient for individual fitness in the population.

.

Note: The implementation of getMonoGameScore(int) is optional. Monomorphic scores are used as visual references in the GUI.

.

For example, in the fitness view, the monomorphic scores are shown as horizontal lines, or used to mark bins in the fitness histogram, and to colour individuals that obtain the monomorphic scores in the structure view.

Optionally, we may implement the methods getPayoff(int me, int you) to return the payoff of an individual with trait me against an individual with trait you and similarly getPayoff(double payoff, int me, int you) to set the payoff to an individual with trait me against one with trait you to payoff. Those methods are available in the TBT module even though they are presently unused.

13.5 Adding IBS model: Individual based simulations

In order to run individual based simulations (IBS) we need to implement the interface HasIBS.DPairs, where DPairs indicates has discrete traits and deals with interactions in pairs. Thus, import the interface and implement it in the class definition:

8 import org.evoludo.simulator.models.IBS.HasIBS;

13 public class TxT extends Discrete implements Payoffs,

14 HasIBS.DPairs {

This requires implementing two additional methods:

pairScores(int me, int[] traitCount, double[] traitScore): to calculate the scores of a focal individual with trait me against all other traits given the counts of each trait in the array traitCount. The method returns the accumulated score of the focal individual. The payoffs that each trait scores against trait me are returned in the array traitScore.
mixedScores(int[] traitCount, double[] traitScore): to calculate the scores of each individual against all other individuals in the population. This only applies to unstructured (or well-mixed) populations where every individual may interact with every other. The state of the population is given by the counts of each trait in the array traitCount and the average scores of each trait are returned in the second array traitScore.

First, let VS Code create the stubs by hovering over the red squiggly line and selecting Quick fix… followed by Add unimplemented methods. After moving the methods to the middle of the class (after the getMonoPayoff(int) method), we have:

72 @Override

73 public double pairScores(int me, int[] traitCount, double[] traitScore) {

74 // TODO Auto-generated method stub

75 throw new UnsupportedOperationException("Unimplemented method ’pairScores’");

76 }

77

78 @Override

79 public void mixedScores(int[] traitCount, double[] traitScore) {

80 // TODO Auto-generated method stub

81 throw new UnsupportedOperationException("Unimplemented method ’mixedScores’");

82 }

The field me in the method pairScores(int me, int[] traitCount, double[] traitScore) denotes the index of the focal individual’s trait. Thus, to determine the payoffs of other traits against me, we simply need to retrieve the corresponding entries in the payoff matrix payoffs. Similarly, the payoff to me is the accumulated score of all interactions against each trait with counts given in traitCount.

72 @Override

73 public double pairScores(int me, int[] traitCount, double[] traitScore) {

74 traitScore[A] = payoffs[A][me];

75 traitScore[B] = payoffs[B][me];

76 return traitCount[A] * payoffs[me][A] + traitCount[B] * payoffs[me][B];

77 }

Implementing the mixedScores(...) method is similarly straight forward to calculate the average payoff in a well-mixed population, where every individual interacts with everyone else, with payoffs returned in traitScore.

.

Important: Individuals do not interact with themselves. Hence, individuals with trait

A

see one less individual with trait

A

among their peers than individuals with trait

B

, and vice versa.

.

79 @Override

80 public void mixedScores(int[] traitCount, double[] traitScore) {

81 int x = traitCount[A];

82 int y = traitCount[B];

83 double in = 1.0 / (x + y - 1);

84 traitScore[A] = ((x - 1) * payoffs[A][A] + y * payoffs[A][B]) * in;

85 traitScore[B] = (x * payoffs[B][A] + (y - 1) * payoffs[B][B]) * in;

86 }

This completes the first stage of our new module, and we are ready to run our first simulation!

Reload the browser window to recompile the new code and load the TxT module. The error message ERROR: No model found! is gone and replaced by the more encouraging message:

Model ’IBS: individual based simulations’ loaded

Model reset

Next is to visualize the state of the simulation.

13.6 Adding GUI views

Currently, the only information that our TxT module provides on the state of the individual based simulation is a console displaying the state of the module and model as well as to log errors, warnings and miscellaneous information about the state of the module. Adding more sophisticated views to visualize the state of the population is often as simple as implementing a particular interface.

13.6.1 Traits and fitness in 2D and 3D

In order to add a 2D visualization of the structure of the population as well as the traits of all individuals, TxT merely needs to implement the interface HasPop2D.Traits:

13 public class TxT extends Discrete implements Payoffs,

14 HasIBS.DPairs,

15 HasPop2D.Traits {

with

import org.evoludo.simulator.views.HasPop2D;

No further actions are required! Now try out the new view by reloading the browser window. The colours of the nodes indicates the trait (or strategy) of the individual. Alternatively, the colours could indicate the fitness of each individual. This is easily achieved by implementing another interface: HasPop2D.Fitness. After reloading the browser window there are now three views available: the Console, the view of the traits, Traits – 2D Structure and now the view of the fitness, Fitness – 2D Structure.

With little more additional effort views of the population structure in 3D can be added as well. Simply implement the interfaces HasPop3D.Traits and HasPop3D.Fitness in the TxT class. Exploring the network structures in 3D is visually interesting and a particularly pleasing extension for the nerdy folks out there that still have a pair of cyan-magenta 3D glasses lying around or, even better, a Google cardboard VR viewer for your phone. Check out the options in the context menu. Don’t be shy and give it a try!

This brings the number of visualizations available in the TxT module to a total of four views(plus the console), each triggered by implementing the corresponding interface:

14 public class TxT extends Discrete implements Payoffs,

15 HasIBS.DPairs,

16 HasPop2D.Traits, HasPop2D.Fitness, HasPop3D.Traits, HasPop3D.Fitness {

with

import org.evoludo.simulator.views.HasPop3D;

The fitness and trait views, each in 2D and 3D, are the most important ones to visualize the structure of the population as well as its current state. In addition, the trait views also allow to interactively manipulate the current state (change traits through double-clicks or double taps). A detailed description of all the ways you can interact with the different kinds of views is given in section 11.1.2: Views, including the 2D Structure and 3D Structure views.

Exercise 13.1:

Consider a population of $20\%$ $A$ and $80\%$ $B$ types on a random graph with $N=200$ nodes and $k=4$ neighbours, on average. The two types engage in a donation game (c.f. section 4.7.4: Donation game) with cost-to-benefit ratio $r=0.1$ , which results in the rescaled payoff matrix

\displaystyle\bordermatrix{&A&B\cr A&1&0\cr B&1+r&r},

(13.1)

where $A$ denotes cooperation and $B$ defection. Plot the population structure in 2D. This is the first in a series of exercises to explore the $b>c\,k$ -rule in structured populations (Ohtsuki et al., 2006).

.

Note: Click Settings to open a text field where the command line options can be entered. Activate the new settings by clicking Apply. A brief description of all available options is available by clicking Help or consulting section 11.2: Parameters.

.

Hint Solution

13.6.2 Mean traits and fitness

While information on the current state and structure of the population is important, interesting and illuminating, it’s the changes over time that rank among the most relevant and informative aspects of the simulated evolutionary process. Luckily, all that is needed to visualize the mean frequency and the mean payoff of each trait over time is to implement the interfaces HasMean.Traits and HasMean.Fitness in the TxT class. This is done by adding the following lines to the class definition:

15 public class TxT extends Discrete implements Payoffs,

16 HasIBS.DPairs,

17 HasPop2D.Traits, HasPop2D.Fitness, HasPop3D.Traits, HasPop3D.Fitness,

18 HasMean.Traits, HasMean.Fitness {

with

import org.evoludo.simulator.views.HasMean;

The two additional views, Traits – Mean and Fitness – Mean, are available in the GUI to visualize the mean frequency and the mean fitness of each trait as well as the mean fitness of the population over time.

Exercise 13.2:

Continuing with Exercise 13.1, draw the frequency of $A$ and $B$ types for Moran death-birth updating over time until one type goes extinct. Use a convex payoff-to-fitness mapping $f_{i}=(1-w)b+w\pi_{i}$ for the fitness of an individual with trait $i$ and payoff $\pi_{i}$ , where $b=1$ denotes the baseline fitness, and $w=0.01$ the selection strength.
Hint Note, your trajectory looks almost certainly different due to the stochastic nature of the evolutionary process. In order to generate identical results, set the seed of the random number generator by adding --seed 0 to the command line options. Traits – Mean
Solution Set command line options to: --module TxT --popsize 200 --geometry R4 --popupdate dB --init frequency 2,8 --paymatrix 1,0;1.1,0.1 --fitnessmap convex 1,0.01 --view 3 and click Apply followed by Start. . . Note: Alternatively, add the option --run to the command line to immediately launch the simulations. Instead of setting --view 3 you can also switch manually to the Traits – Mean view. .

Note that the protection against exploitation provided by the spatial structure enables cooperators ( $A$ types) to reach fixation most of the time. However, because selection is fairly weak the process is fairly noisy with only a small bias in favour of cooperation based on the limited local interactions on the graph.

13.6.3 Fitness histogram

The mean fitness values are sufficient to describe the dynamics in unstructured (or well-mixed) populations. However, in structured populations the distribution of fitness provides additional insights into the evolutionary dynamics. To visualize the distribution of fitness values for each trait, a histogram of the respective fitness distributions is available when implementing the interface HasHistogram.Fitness:

17 public class TxT extends Discrete implements Payoffs,

18 HasIBS.DPairs,

19 HasPop2D.Traits, HasPop2D.Fitness, HasPop3D.Traits, HasPop3D.Fitness,

20 HasMean.Traits, HasMean.Fitness, HasHistogram.Fitness {

with

import org.evoludo.simulator.views.HasHistogram;

In each histogram the marked bin corresponds to the fitness of a population that is monomorphic in the respective trait.

Exercise 13.3:

Continuing with Exercise 13.2, reset the simulation and draw the fitness distribution after ten generations.
Hint Note, your distribution looks likely different due to the stochastic nature of the evolutionary process. In order to generate identical results, set the seed of the random number generator by adding --seed 0 to the command line options. Histogram - Fitness
Solution Set command line options to: --module TxT --popsize 200 --geometry R4 --popupdate dB --init frequency 2,8 --paymatrix 1,0;1.1,0.1 --fitnessmap convex 1,0.01 --view 7 and click Apply followed by Start. Quickly click Stop to halt the simulations again. . . Note: Add the option --run on the command line to immediately launch the simulation. Alternatively, you can simply click Step ten times. Moreover, the simulations can be stopped automatically after ten generations by adding --generations 10 to the command line options. Finally, instead of setting --view 7 you can switch manually to the Histogram - Fitness view. .

13.6.4 Degree histogram

An interesting and relevant measure to characterize heterogeneous population structures is to consider the degree distribution of nodes in the population. The degree of a node is the number of its neighbours. In directed graphs the in-degree is the number of incoming edges and the out-degree is the number of outgoing edges. In undirected graphs the in- and out-degree are the same. Moreover, the structure of the interaction network may be different from the reproduction network.

The degree distribution depicts the proportion of nodes with a particular degree and visualizes the heterogeneity of the population structure. For example, in scale-free networks the degree distribution follows a power-law, which means that the probability of finding a node with a degree $k$ is proportional to $k^{-\gamma}$ , where $\gamma$ is the power-law exponent. For random graphs the degree distribution is a binomial distribution, which is peaked around the average degree $\bar{k}$ . Finally, for lattices (with periodic boundaries) or random regular graphs the degree distribution has a single peak at $k$ because all nodes have the same degree. In general, the degree distribution provides information about the heterogeneity of the population structure.

In order to visualize the degree distribution of the population structure, all that is needed is to implement the interface HasHistogram.Degree:

17 public class TxT extends Discrete implements Payoffs,

18 HasIBS.DPairs,

19 HasPop2D.Traits, HasPop2D.Fitness, HasPop3D.Traits, HasPop3D.Fitness,

20 HasMean.Traits, HasMean.Fitness, HasHistogram.Fitness, HasHistogram.Degree {

Exercise 13.4:

Extending on Exercise 13.1, plot the degree distribution of the population structure.
Hint Note, your degree distribution looks likely slightly different because the generation of the network structure involves drawing links between randomly drawn nodes. In order to generate an identical structure, set the seed of the random number generator by adding --seed 0 to the command line options. Histogram - Degree
Solution Set command line options to: --module TxT --popsize 200 --geometry R4 --view 8 [--popupdate dB --init frequency 2,8 --paymatrix 1,0;1.1,0.1 --fitnessmap convex 1,0.01] and click Apply. . . Note: The parameters in brackets are part of Exercise 13.1 and Exercise 13.2 but do not affect the degree distribution. Instead of setting --view 8 you can switch manually to the Histogram - Degree view. .

13.7 Adding statistics modes

The EvoLudo GUI supports two kinds of running modes: dynamics and statistics. The default is the dynamics mode, which simply reports the progress of the evolutionary process over time and is visualized in intervals set by the --timestep command line option. Statistics mode works slightly different in that samples are collected for statistical evaluations. At present two complementary sampling types are supported: fixation statistics and stationary distributions. For fixation statistics, each sample represents a run starting from a well-defined initial configuration that resulted in the fixation of one trait in the population. Naturally this requires that fixation is possible, that is, absorbing states must exist. Fixation is not possible, for example, in the presence of mutations. The initial configuration is typically a monomorphic resident population with a single mutant placed in a random location set by the --init command line option.

In contrast, stationary distributions are obtained by sampling the population at regular intervals, again set by the --timestep command line option, and recording the frequencies of each trait. The long-term average over the samples converges to the stationary distribution. This requires that no absorbing states exist, which can be easily achieved by adding small mutation rates (see --mutation option).

The statistical results are visualized as histograms that are continuously updated as more samples are collected. Neither statistics mode requires any particular consideration in the code of the module. Instead, all that is required is to implement the corresponding interfaces:

HasHistogram.StatisticsProbability:: record the probability of fixation for each trait.
HasHistogram.StatisticsTime:: record the expected time until fixation for each trait plus the time to absorption, that is, the time until the population turns monomorphic, regardless of which trait fixated. For a detailed discussion about the connection between fixation probabilities as well as fixation and absorption times see sections 2.5 & 5.4.
HasHistogram.StatisticsStationary:: record the probability to find the population in any particular state, specified by the number of mutants. This probability approximates the relative time spent in this particular state.

.

Important: Whether and what kind of statistics modes are available in the GUI is not only determined by the interfaces implemented by a module but also by the command line options. For example, with mutations no absorbing states exist and hence fixation can never occur. Accordingly, only stationary distributions are available. Conversely, without mutations absorbing states exist and fixation is inevitable, eventually. Thus, fixation statistics are available provided that a compatible initial configuration is set.

.

After this long introduction, the actual implementation of the statistics modes is as easy as adding another view to the GUI:

17 public class TxT extends Discrete implements Payoffs,

18 HasIBS.DPairs,

19 HasPop2D.Traits, HasPop2D.Fitness, HasPop3D.Traits, HasPop3D.Fitness,

20 HasMean.Traits, HasMean.Fitness, HasHistogram.Fitness,

21 HasHistogram.StatisticsProbability, HasHistogram.StatisticsTime,

22 HasHistogram.StatisticsStationary {

Exercise 13.5:

Building on Exercise 13.2, calculate the distribution of fixation probabilities for mutants arising in all nodes based on $10^{\prime}000$ samples. Note that this requires changing the settings for the initial configuration to a single $B$ mutant in a resident $A$ population. Also, be aware that generating the samples may require a few minutes.
Hint Note, your distribution looks likely different due to the stochastic nature of the evolutionary process. In order to generate identical results, set the seed of the random number generator by adding --seed 0 to the command line options. Statistics - Fixation probabilities
Solution Set command line options to: --module TxT --popsize 200 --geometry R4 --popupdate dB --init mutant 0,1 --paymatrix 1,0;1.1,0.1 --fitnessmap convex 1,0.01 --view 9 --samples 10000 and click Apply followed by Start. . . Note: Compiling $10^{\prime}000$ samples may take a few minutes. Instead of setting --view 9 you can switch manually to the Statistics - Fixation probabilities view. .

While, overall, the results may not look in favour of cooperation, things need to be put into proper perspective. The average fixation probability of a single cooperator, $A$ , in a homogeneous defector population, $B$ , is approximately $0.59\%$ . In contrast, the fixation probability of a single trait under neutral selection is $1/N$ and hence for a population of size $N=200$ amounts to $0.5\%$ . Indeed, in comparison, the spatial structure promotes cooperation. Moreover, keep in mind that we are considering weak selection with small fitness differences and hence the results are expected to be close to the neutral case. Of course, to gain more confidence in the significance of this difference we would need to increase the sample size to $10^{6}$ or more. Nevertheless, the finding is in line with the $b>ck$ -rule reported in Ohtsuki et al. (2006). More specifically, in our case, the cost-to-benefit ratio is $r=c/b=0.1$ and the (average) number of neighbours is $\bar{k}=4$ and hence the rule applies.

.

Note: The statistical data of the simulation runs can be exported in CSV format by choosing Export… and Statistics (csv) from the context menu of any statistics view.

.

The histogram for the fixation probabilities shows significant variation, which has two independent sources: first, the location of the initial mutant on the regular graph and second the structure of the random graph itself. In heterogeneous graphs the fixation probability sensitively depends on the location of the initial mutant as well as the overall structure of the graph. In Exercise 13.5 the random graph was regenerated for each sample before placing the mutant uniformly at random. As a consequence the distribution of fixation probabilities across nodes should eventually converge to a constant for a sufficiently large number of samples.

Exercise 13.6:

Repeat the simulations in Exercise 13.5, but this time use a single graph for all samples. Now any systematic variability in the fixation probabilities is due to the structural features of the regular graph.
Hint Note, your distribution looks likely different due to the stochastic nature of the evolutionary process. In order to generate identical results, set the seed of the random number generator by adding --seed 0 to the command line options. Statistics - Fixation probabilities
Solution Set command line options to: --module TxT --popsize 200 --geometry R4 --popupdate dB --init mutant 0,1 --paymatrix 1,0;1.1,0.1 --fitnessmap convex 1,0.01 --samples 10000 --statistics reset 0 --view 9 and click Apply followed by Start.
Exercise 13.7:

Run the simulations in Exercise 13.6 but make sure to include the --seed 0 option to guarantee reproducible results. Moreover, in order to reduce the noise in the evolutionary process, increase the selection strength to the maximum of $w=1$ . This is the maximum because the minimum payoff is zero. Note that the fixation probability sensitively depends on the location of the initial mutant. In particular, with the above settings it turns out that mutants arising in nodes with indices $8,10,22,57,68,105,135,198$ have all remarkable fixation probabilities of $20\%$ or higher, with a maximum of $30.43\%$ for node $135$ .

. . Note: Hover with the mouse over the bars (or tap on a bar) to view a tooltip with more information .

Keep in mind that the statistical power remains limited because even for $10^{\prime}000$ samples this leaves only an average of $50$ samples per node. Nevertheless, this should be enough to identify relevant patterns and mechanisms.

Now turn to the Traits – 2D Structure view and try to locate these nodes. Even more importantly, figure out why they are special and what their characteristics are. What are your observations? How could those features help to promote cooperation?
Hint The random regular graph used for the fixation probability histogram in Exercise 13.6. Because of the --seed 0 option your graph should look identical. Can you identify which locations might be conducive to the fixation of cooperators and why? Traits – 2D Structure
Solution The key determining feature is that all nodes resulting in high fixation probabilities have (at least) one neighbour which is a leaf node (i.e. has no further neighbours). The leaf node provides a safeguard against accidental extinction of the mutant. At the same time, leaf nodes themselves are too isolated and the mutant unlikely to spread. The neighbours of leaf nodes strike a balance between protection from extinction through their leaf neighbours and opportunities for proliferating.
Exercise 13.8:

An alternate and complementing view on cooperation in structured populations is to consider the stationary distribution of the trait frequencies. Stationary distributions require an ergodic system, which means no absorbing states. In the presence of absorbing states the stationary distribution is either trivial or depends on the initial configuration. The Moran process is easily turned into an ergodic process by simply adding a small mutation rate, say $\mu=0.001$ . This translates into, on average, one mutation every five generations in a population of size $N=200$ . In addition, change the initial condition to start with a monomorphic defector population, $B$ types. Derive the stationary distribution after $10^{\prime}000$ generations (samples).
Hint Note, your distribution may look slightly different due to the stochastic nature of the evolutionary process, but the differences should be essentially imperceptible. In order to generate identical results, set the seed of the random number generator by adding --seed 0 to the command line options. Statistics - Stationary distribution
Solution Set the command line options to: --module TxT --popsize 200 --geometry R4 --popupdate dB --init mono 1 --mutation 0.001 --paymatrix 1,0;1.1,0.1 --fitnessmap convex 1,0.01 --timestop 10000 and click Apply followed by Start . . Note: In order to collect the samples at maximum speed, move the slider to the far right or set --delay 0. .

13.8 Advanced topics

This section covers several, more advanced topics, which go beyond the basic functionalities of an EvoLudo module but serve to illustrate the flexibility and extensibility of the EvoLudo framework.

13.8.1 Implementing evolutionary kaleidoscopes

As illustrated at the very beginning of this book (see section 1.2: Gallery) some types of interactions spatially structured populations give rise to intriguing spatio-temporal patterns (Nowak and May, 1992). In particular, given the appropriate settings, $2\times 2$ games exhibit fascinating evolutionary kaleidoscopes, or patterns that are reminiscent of fractal structures or Persian carpets and even change dynamically over time (see e.g. Lab. 1.1).

In order to create evolutionary kaleidoscopes very specific settings for the simulations are required. Most importantly, the population must be structured, the updating of the population must be deterministic and requires very particular initial configurations that are symmetric. Structured populations ensure that different individuals experience distinct environments, that is, have sets of neighbours with different traits. The simplest case are square lattices with von Neumann, $k=4$ , or Moore, $k=8$ , neighbourhoods. Deterministic updates are achieved if:

1.

the population is updated synchronously, that is, all individuals update their trait at the same time;
2.

every individual interacts with all its neighbours;
3.

every individual adopts the trait that performs best among all its neighbours, including itself.

For any symmetrical initial configuration the deterministic updating guarantees that the trait distribution in the population remains symmetrical at all times. The simplest case is given by a monomorphic $A$ population with a single $B$ placed at the centre.

Finally, kaleidoscopes only emerge for special types of interactions, that is particular payoff matrices. The most prominent example is the so-called weak prisoner’s dilemma (see section 6.1.1: Evolutionary Kaleidoscopes for more details) with a payoff matrix given by:

\displaystyle\bordermatrix{&A&B\cr A&1&0\cr B&b&0},

(13.2)

such that a player is indifferent between traits $A$ and $B$ against a co-player opting for $B$ . For $b>1$ the trait (or strategy) $B$ is weakly dominant and $b$ denotes the temptation to defect.

On square lattices with Moore neighbourhood evolutionary kaleidoscopes arise for $5/3>b>8/5$ , while those with von Neumann neighbourhood they emerge for $3/2>b>4/3$ (Nowak and May, 1992).

Given this long list of requirements it bears the question of why bother with evolutionary kaleidoscopes at all? The answer is simple: they may not carry any biological relevance, but they are fascinating and prime examples of the beauty in mathematics.

Assuming that the reader agrees with this assessment, let’s add the necessary code to allow for customized initial configurations. This requires a deeper access to the inner workings of individual based simulation models.

Customizing IBS populations

In order to get started, let us create our own implementation of an individual based simulation. Subclassing IBSDPopulation as an inner class of TxT will do the trick.

.

Note: Using an inner class has the advantage that it maintains access to all fields and methods of the outer class TxT without the need to pass them as parameters.

.

Add the following code to the end of the TxT class:

124 public class IBSPop extends IBSDPopulation {

125 }

with

import org.evoludo.simulator.models.IBSDPopulation;

This triggers an error because the IBSDPopulation does not provide a default constructor, i.e. one without arguments. Thus, we need to explicitly implement one, which refers to the constructor IBSDPopulation(EvoLudo engine, Discrete module). Let VS Code create the constructor stub by hovering over the error, marked by a red squiggly line, and select Quick fix… followed by Add constructor ’IBSPop(EvoLudo engine, Discrete module)’ to obtain:

125 public class IBSPop extends IBSDPopulation {

126

127 public IBSPop(EvoLudo engine, Discrete module) {

128 super(engine, module);

129 //TODO Auto-generated constructor stub

130 }

131 }

.

Note: Because we are using an inner class, the constructor of IBSPop can be simplified to only take the TxT module as a parameter. The field engine is known by the outer class TxT and hence can be accessed directly using TxT.this.engine. this by itself refers to the current instance of IBSPop. Thus, in order to access the instance of the outer class we need to use TxT.this.

.

Our simplified constructor now reads:

127 public IBSPop(TxT module) {

128 super(TxT.this.engine, module);

129 }

As the name of the super class suggests (reading from the back), it provides the framework to model a Population of individuals with Discrete traits through individual based simulations, IBS.

Now we are ready to link our newly minted subclass into the TxT module. The next step is to replace the default EvoLudo implementation for individual based simulation with our own. This is done by overriding the method createPopulation() in the TxT class and return our own implementation:

125 @Override

126 public TxT.IBSPop createIBSPop() {

127 return new TxT.IBSPop(this);

128 }

It is easy to check that our subclass gets indeed loaded by setting a break point on the constructor. After re-loading your browser window execution will stop at this line. Remove the break point again and click continue to proceed.

Customizing initial configurations

Before linking our newly minted subclass into the TxT module, let’s add a method to set the initial configuration of the population. For evolutionary kaleidoscopes some groundwork is already implemented, and we only need to implement the method initKaleidoscope(). This method should first check whether the population has an appropriate lattice structure and report a warning if not. On any lattice it then places a single $B$ mutant in the centre of a monomorphic $A$ population:

131 @Override

132 protected void initKaleidoscope() {

133 if (!interaction.isLattice()) {

134 logger.warning("kaleidoscopes require lattices - using default initialization.");

135 initUniform();

136 return;

137 }

138 int mid = -1;

139 initMono(TxT.A);

140 switch (interaction.getType()) {

141 case CUBE:

142 int side = (int) (Math.pow(nPopulation, 1.0 / 3.0) + 0.5);

143 int side2 = side * side;

144 mid = (side / 2) * (side2 + side + 1);

145 break;

146

147 case SQUARE_NEUMANN:

148 case SQUARE_NEUMANN_2ND:

149 case SQUARE_MOORE:

150 case SQUARE:

151 case HONEYCOMB:

152 case TRIANGULAR:

153 side = (int) Math.sqrt(nPopulation);

154 mid = (side / 2) * (side + 1);

155 break;

156

157 case LINEAR:

158 mid = nPopulation / 2;

159 break;

160

161 default:

162 throw new Error("geometry incompatible with kaleidoscopes!");

163 }

164 traitsCount[TxT.A]--;

165 setTraitAt(mid, TxT.B);

166 traitsCount[TxT.B]++;

167 }

.

Note: In principle, the method initKaleidoscope() can be used to provide any kind of special initial configuration – regardless of whether they are capable of producing evolutionary kaleidoscopes.

.

The last step is to modify the --init command line option to include a key that triggers the initial configuration for kaleidoscopes. This is done by adjusting the command line options by overriding the adjustCLO(CLOParser) method in the TxT class:

125 @Override

126 public void adjustCLO(CLOParser parser) {

127 super.adjustCLO(parser);

128 if (model instanceof IBSD) {

129 CLOption clo = ((IBSDPopulation) getIBSPopulation()).getInit().clo;

130 clo.addKey(Init.Type.KALEIDOSCOPE);

131 }

132 }

with

import org.evoludo.simulator.models.IBSD;

import org.evoludo.simulator.models.IBSD.Init;

This adds the key kaleidoscope to the --init command line option. Reload your browser window to recompile the module and click Help to verify the addition of the new key.

Exercise 13.9:

Create your first evolutionary kaleidoscope with EvoLudo on a $101\times 101$ square lattice with Moore neighbourhood.
Hint Note, the snapshot of the kaleidoscope is taken after $49$ generations. Traits – 2D Structure
Solution Set the command line options to: --module TxT --popsize 101x101 --geometry m --popupdate sync --playerupdate best --interactions all --references all --init kaleidoscope --paymatrix 1,0;1.65,0 and click Apply.

13.8.2 Adding square lattice geometry with second neighbours

On a square lattice every node has four nearest neighbours (to the north, east, south and west), the von Neumann neighbourhood, and also four second-nearest neighbours (on either side along the two diagonals). Combining first- and second-nearest neighbours results in the Moore neighbourhood, which includes all eight neighbours reachable by a chess-kings move, see Fig. 13.2.

Refer to caption — Figure 13.2: Square lattice. a The nearest neighbours of the focal node, $\bullet$ (in the centre), are shown as $\bullet$ and the second neighbours as $\bullet$ . b The nearest neighbours form the standard 2D square lattice with degree $k=4$ and includes all members of the population. c The second neighbours also form a square lattice with the same degree. However, in contrast to the first neighbours, they form two disjoint sub-lattices (black and red dashed lines, respectively) of equal sizes, provided that $N$ is even on an $N\times N$ lattice.

Interestingly, however, the second neighbours form two disjoint sub-lattices, which can give rise to intriguing dynamical features. More specifically, in the spatial donation game where individuals interact with their four nearest neighbours but compete with, or imitate, their four second-nearest neighbours, spontaneous symmetry breaking may occur below a critical cost-to-benefit ratio $r$ . As a result, cooperation is more abundant in one sub-lattice than in the other. Rather surprisingly, this asymmetry increases with decreasing $r$ – even though conditions for cooperation become less challenging (smaller costs or increased benefits) – and eventually leads to complete asymmetry with the fixation of one trait in each sub-lattice. This is a neat example of a critical phase transition in evolutionary game theory (Szabó and Hauert, 2002, Hauert and Szabó, 2024, 2025).

The easy part is implementing the geometry with second neighbours. This is already accomplished in the Geometry class (see SQUARE_NEUMANN_2ND geometry). Similarly, a customized graphical visualization is available, which reveals the two sub-lattices. We will explore this further once the trickier part is dealt with for dealing with the two sub-lattices when reporting the state of the population. In principle, simulations using second neighbour lattices for interaction or competition structures work out of the box. However, without proper reporting it would be hard if not impossible to quantify the dynamics in each sub-lattice, in general, and the asymmetry between sub-lattices, in particular. For this reason the SQUARE_NEUMANN_2ND geometry is disabled by default.

The goal of this section is to enable the SQUARE_NEUMANN_2ND geometry and to report the state of each sub-lattice separately in terms of traits as well as fitness. First off, to enable the SQUARE_NEUMANN_2ND geometry, we simply need to add the corresponding key to the parser for the command line option --geometry for IBS models. If you skipped the previous section on evolutionary kaleidoscopes, you will need to add the following code:

127 @Override

128 public void adjustCLO(CLOParser parser) {

129 super.adjustCLO(parser);

130 if (model instanceof IBSD) {

131 cloGeometry.addKey(Geometry.Type.SQUARE_NEUMANN_2ND);

132 }

133 }

with

import org.evoludo.simulator.models.Geometry;

Otherwise, simply insert the line

134 cloGeometry.addKey(Geometry.Type.SQUARE_NEUMANN_2ND);

after clo.addKey(Init.Type.KALEIDOSCOPE) in the adjustCLO(CLOParser). This will make the SQUARE_NEUMANN_2ND geometry available in the GUI, which can be verified by reloading the browser window, clicking on Help and checking the list of available geometries.

Customizing trait names for graphics and tooltips

Now, let us define different names for the traits $A$ and $B$ in each sub-lattice, say $A_{1},B_{1}$ and $A_{2},B_{2}$ , respectively. This requires two changes: First we need to extend the method getTraitName(int) in our TxT module for the SQUARE_NEUMANN_2ND geometry to return modified trait names for one sub-lattice and add two additional trait names (with indices $2$ and $3$ ) for the corresponding traits in the other sub-lattice:

100 @Override

101 public String getTraitName(int idx) {

102 String idxname = super.getTraitName(idx % nTraits);

103 if (competition.getType() != Geometry.Type.SQUARE_NEUMANN_2ND)

104 return idxname;

105 if (idx >= nTraits)

106 return idxname + "<sub>2</sub>";

107 return idxname + "<sub>1</sub>";

108 }

.

Note: HTML formatting is acceptable for most strings, including trait names. Exceptions are the keys and descriptions of the command line options, which must be plain text.

.

Second, when reporting the trait names, the tooltips of the 2D and 3D lattice views need to take the sub-lattice into account on which the individual resides in the SQUARE_NEUMANN_2ND geometry. If you skipped the previous section on evolutionary kaleidoscopes, you will need to first create a custom IBS population as described in the section section 13.8.1: Customizing IBS populations above. The inner class TxT.IBSPop then needs to override the method getTraitNameAt(int):

159 @Override

160 public String getTraitNameAt(int idx) {

161 if (competition.getType() != Geometry.Type.SQUARE_NEUMANN_2ND)

162 return super.getTraitNameAt(idx);

163 int side = (int) Math.sqrt(nPopulation);

164 int trait = getTraitAt(idx);

165 if ((idx / side) % 2 == (idx % side) % 2)

166 return module.getTraitName(trait);

167 return module.getTraitName(nTraits + trait);

168 }

Now is a good time to check our progress. Reload your browser window and click on Help to verify that the SQUARE_NEUMANN_2ND geometry is available, enter the settings --module TxT --geominter n --geomcomp n2 and click Apply to check the visualization of the geometry and the tooltips. In the 2D the nodes in one sub-lattice are represented by squares and those in the other by disks. In 3D the cubes representing the nodes of the sub-lattices are vertically displaced such that when looking at the structure at a low angle only one or the other sub-lattice is visible.

Customizing the status bar

Next, let’s update the status bar to report the number of individuals in each sub-lattice. In preparation, we need to override the method getMeanTraits(double[]) and calculate the mean frequency of traits in each sub-lattice separately. Since this procedure requires additional CPU resources and the mean frequencies are also used by the Traits – Mean let us store the results in a local field, together with the time stamp of the last update. This allows us to only update the values if the time stamp is outdated. This is done by adding the following code to the TxT.IBSPop class:

171 @Override

172 public void getMeanTraits(double[] mean) {

173 if (competition.getType() != Geometry.Type.SQUARE_NEUMANN_2ND) {

174 super.getMeanTraits(mean);

175 return;

176 }

177

178 double newtime = model.getTime();

179 if (Math.abs(tsMean - newtime) < model.getTimeStep()) {

180 System.arraycopy(mean2nd, 0, mean, 0, mean.length);

181 return;

182 }

183 int n = 0;

184 Arrays.fill(mean, 0);

185 int side = (int) Math.sqrt(nPopulation);

186 while (n < nPopulation) {

187 if ((n / side) % 2 == 0) {

188 mean[getTraitAt(n++)]++;

189 mean[nTraits + getTraitAt(n++)]++;

190 continue;

191 }

192 mean[nTraits + getTraitAt(n++)]++;

193 mean[getTraitAt(n++)]++;

194 }

195 ArrayMath.multiply(mean, 2.0 / nPopulation);

196 System.arraycopy(mean, 0, mean2nd, 0, mean.length);

197 tsMean = newtime;

198 }

Of course this also requires defining the local double array mean2nd to store the mean frequencies and the double tsMean to save the time stamp of the last calculation:

156 double tsMean = -1.0;

157 double[] mean2nd;

Their initialization is done by overriding the check() method, which is called after the initial loading of the module or after clicking Init, Reset, or Apply:

156 @Override

157 public boolean check() {

158 boolean doReset = super.check();

159 if (competition.getType() != Geometry.Type.SQUARE_NEUMANN_2ND) {

160 mean2nd = null;

161 } else {

162 mean2nd = new double[2 * nTraits];

163 }

164 tsMean = -1.0;

165 return doReset;

166 }

Now we are ready to show the information in the status bar. Override the method getStatus():

215 @Override

216 public String getStatus() {

217 if (competition.getType() != Geometry.Type.SQUARE_NEUMANN_2ND)

218 return super.getStatus();

219

220 getMeanTraits(mean2nd);

221 String status = "";

222 for (int i = 0; i < 2 * nTraits; i++)

223 status += (status.length() > 0 ? ", " : "") + module.getTraitName(i) + ": "

224 + Formatter.formatPercent(mean2nd[i], 1);

225 return status;

226 }

with

import org.evoludo.util.Formatter;

Verify that the status bar works by reloading your browser window. Make sure that the settings include --module TxT --geominter n --geomcomp n2. Click Step to update the status information.

Customizing the views for mean trait frequencies and fitness

With these changes in place we are almost ready to also visualize the frequencies of the traits in each sub-lattice over time in the Traits – Mean view. All that is needed is to advertise that our simulation, the TxT.IBSPop class, provides four mean values instead of the default two in case of the SQUARE_NEUMANN_2ND geometry:

187 @Override

188 public int getNMean() {

189 if (competition.getType() != Geometry.Type.SQUARE_NEUMANN_2ND)

190 return super.getNMean();

191 return 2 * nTraits;

192 }

As you can easily verify, the Traits – Mean view now indeed reports four mean values for the frequencies of $A$ and $B$ traits in each sub-lattice. However, you will also notice that now the colours of the lines are off. It would be nicer to use paler versions of the $A$ and $B$ colours to distinguish their frequencies in the two sub-lattices. In order to do this, override the method getMeanColor(int) in the TxT module and simply add some transparency to the existing trait colours:

64 @Override

65 public Color[] getMeanColors() {

66 Color[] colors = super.getMeanColors();

67 if (competition.getType() != Geometry.Type.SQUARE_NEUMANN_2ND)

68 return colors;

69 Color[] color2nd = new Color[2 * nTraits];

70 for (int n = 0; n < nTraits; n++) {

71 Color cn = colors[n];

72 color2nd[n] = cn;

73 color2nd[nTraits + n] = ColorMap.addAlpha(cn, 150);

74 }

75 return color2nd;

76 }

with

import org.evoludo.util.ColorMap;

Note, at this point our view Fitness – Mean is slightly messed up because, based on the return value of getNMean(), it now expects four mean values instead of two (plus the overall mean fitness). This is addressed by overriding the method for the fitness values, getMeanFitness(double[]). Similar to getMeanTraits(double[]) we calculate the mean fitness for $A$ and $B$ types in the two sub-lattices separately.

238 @Override

239 public void getMeanFitness(double[] mean) {

240 if (competition.getType() != Geometry.Type.SQUARE_NEUMANN_2ND) {

241 super.getMeanFitness(mean);

242 return;

243 }

244

245 double newtime = model.getTime();

246 if (Math.abs(tsFit - newtime) < model.getTimeStep()) {

247 System.arraycopy(fit2nd, 0, mean, 0, mean.length);

248 return;

249 }

250 int n = 0;

251 Arrays.fill(mean, 0);

252 int side = (int) Math.sqrt(nPopulation);

253 while (n < nPopulation) {

254 if ((n / side) % 2 == 0) {

255 mean[getTraitAt(n)] += getFitnessAt(n++);

256 mean[nTraits + getTraitAt(n)] += getFitnessAt(n++);

257 continue;

258 }

259 mean[nTraits + getTraitAt(n)] += getFitnessAt(n++);

260 mean[getTraitAt(n)] += getFitnessAt(n++);

261 }

262 mean[2 * nTraits] = sumFitness * 0.25;

263 ArrayMath.multiply(mean, 2.0 / nPopulation);

264 System.arraycopy(mean, 0, fit2nd, 0, mean.length);

265 tsFit = newtime;

266 }

As before, this requires introducing two new fields, the time stamp for fitness calculations, tsFit, and the double array fit2nd to store the fitness,

174 double tsFit = -1.0;

175 double[] fit2nd;

as well as initializing them in the check() method:

169 if (competition.getType() != Geometry.Type.SQUARE_NEUMANN_2ND) {

170 mean2nd = null;

171 fit2nd = null;

172 } else {

173 mean2nd = new double[2 * nTraits];

174 fit2nd = new double[2 * nTraits + 1];

175 }

176 tsMean = -1.0;

177 tsFit = -1.0;

with the additions highlighted in red.

Exercise 13.10:

Explore the spatial dynamics of cooperation in the two sub-lattices in a population of size $200\times 200$ . The population is updated asynchronously and individuals compete with a single, randomly selected neighbour among the four along the diagonals. The focal individual updates its trait following a probabilistic comparison of their payoff to that of the competitor according to the Fermi update with a noise term of $K=0.5$ . The payoffs are accumulated over interactions with all four nearest neighbours from donation games. All payoffs are ephemeral and solely used for the probabilistic comparison and discarded afterwards. Fairly long relaxation times are required to reach the equilibrium configuration. Produce a snapshot of a typical equilibrium configuration for a cost-to-benefit ratio of $r=0.035$ after $10^{\prime}000$ generations with symmetrical levels of cooperation in the two sub-lattices.
Hint For the symmetrical setup, the equilibrium level of cooperation is around $22\%$ in each sub-lattice. Traits – 2D Structure
Solution Set the command line options to: --module TxT --timestep 1 --geominter n --geomcomp n2 --popsize 200x --popupdate async --interactions all --references random 1 --playerupdate thermal 0.5 --init frequency 1,1 --paymatrix 0.965,-0.035;1,0 --accuscores --resetscores ephemeral --timestop 10000 and click Apply. The same distribution as shown above can be replicated by adding the option --seed 0.
Exercise 13.11:

Explore the dynamics of cooperation in the two sub-lattices for the donation game with cost-to-benefit ratio $r=0.005$ . Use otherwise the same settings as in Exercise 13.10 to produce a snapshot of a typical equilibrium configuration with highly asymmetrical levels of cooperation in the two sub-lattices after $10^{\prime}000$ generations.
Hint For the asymmetrical setup, the equilibrium level of cooperation is around $12\%$ in one sub-lattice and around $88\%$ in the other. Thus, the trait frequencies are largely complementary between lattices. Traits – 2D Structure
Solution Set the command line options to: --module TxT --timestep 1 --geominter n --geomcomp n2 --popsize 200x --popupdate async --interactions all --references random 1 --playerupdate thermal 0.5 --init frequency 1,1 --paymatrix 0.995,-0.005;1,0 --accuscores --resetscores ephemeral --timestop 10000 and click Apply. The same distribution as shown above can be replicated by adding the option --seed 0.

13.9 Adding ODE model:
ordinary differential equations

The traditional way to investigate $2\times 2$ games is based on the replicator dynamics (Hofbauer and Sigmund, 1998). In order to do this in EvoLudo, we first need to implement the interface HasDE.DPairs, to advertise that we intend to offer a model based on differential equations (the HasDE interface) for discrete traits and interactions in pairs (the DPairs interface).

25 public class TxT extends Discrete implements Payoffs,

26 HasIBS.DPairs, HasDE.DPairs,

27 HasPop2D.Traits, HasPop2D.Fitness, HasPop3D.Traits, HasPop3D.Fitness,

28

\vdots

with

15 import org.evoludo.simulator.models.Model.HasDE;

This advertises that our module can handle ordinary differential equations, which requires the implementation of two methods: getDependent() and avgScores(double[] density, double[] avgscores) in our TxT class. Again, you may let VS Code create the method stubs by hovering over the error, marked by a red squiggly line, and select Quick fix… followed by Add unimplemented methods to obtain (after moving the stubs to more appropriate locations):

74 @Override

75 public int getDependent() {

76 // TODO Auto-generated method stub

77 throw new UnsupportedOperationException("Unimplemented method ’getDependent’");

78 }

and

110 @Override

111 public void avgScores(double[] density, double[] avgscores) {

112 // TODO Auto-generated method stub

113 throw new UnsupportedOperationException("Unimplemented method ’avgScores’");

114 }

The second method calculates the average payoffs from pairwise interactions and returns them in the array avgscores for each trait given their frequencies in the state array:

104 @Override

105 public void avgScores(double[] state, double[] scores) {

106 scores[A] = state[A] * payoffs[A][A] + state[B] * payoffs[A][B];

107 scores[B] = state[A] * payoffs[B][A] + state[B] * payoffs[B][B];

108 }

The state of the population is fully determined by the frequencies of the two traits, $A$ and $B$ . More specifically this means that none of the entries in the state array are negative, and they all must sum up to one. Thus, the state space corresponds to a simplex, $S_{2}$ in the present case, and the number of dynamical variables is reduced by one. This is where the first method getDependent() comes in and is used to identify the dependent trait. Let’s declare $B$ as the dependent trait, and $A$ the independent one:

73 @Override

74 public int getDependent() {

75 return B;

76 }

.

Note: The method getDependent() helps reduce numerical rounding errors by ensuring that frequencies always lie in

[0,1]

and sum up to one, i.e. remain confined to the simplex. This results in computational redundancy that can be exploited to reduce numerical rounding errors and improve accuracy.

.

In contrast, models based on densities have no dependent trait and thus numerical precautions are limited to enforcing non-negative values. In those cases the method getDependent() returns $-1$ to indicate that there is no dependent trait.

Now we have finished the groundwork for implementing differential equations and are ready to specify the kinds of ordinary differential equation models that our model should support. EvoLudo supports two kinds of ODE models, again in the form of interfaces:

HasDE.EM:: Simple integrator based on the Euler method.
HasDE.RK5:: More advanced fifth-order Runge-Kutta method with variable step size.

There is no reason not to support both methods, especially because all we need to do is to implement the corresponding interfaces:

25 public class TxT extends Discrete implements Payoffs,

26 HasIBS.DPairs, HasDE.DPairs, HasDE.EM, HasDE.RK5,

27 HasPop2D.Traits, HasPop2D.Fitness, HasPop3D.Traits, HasPop3D.Fitness,

28

\vdots

.

Note: The models are selected by adding the option --model EM or --model RK5 to the command line, respectively. In addition, the generic option --model ODE is also available, which defaults to RK5, if supported, and EM otherwise.

.

Now we are almost done with the implementation of the ODE model. If you skipped section 13.8: Advanced topics you are good to go and may continue with the exercise. However, in case you did complete the advanced section, there is one important detail to consider:

.

Attention: Two of the methods that were added to the TxT class in the advanced section 13.8.2: Adding square lattice geometry with second neighbours, require further considerations for ODE models. More specifically, getTraitName(int idx) and getMeanColors(), implicitly assume that the competition geometry exists, i.e. is not null. While this is the case for IBS models, it does not apply to ODE models (as well as stochastic, SDE, and partial differential equations, PDE, models, see below).

.

In order to prevent access to an undefined field, we need to take suitable precautions by either checking that the field competition is not null, or even better, check the type of the current model:

118 public String getTraitName(int idx) {

119 String idxname = super.getTraitName(idx % nTraits);

120 if (model.isDE() || competition.getType() != Geometry.Type.SQUARE_NEUMANN_2ND)

121

\vdots

and

128 public Color[] getMeanColors() {

129 Color[] colors = super.getMeanColors();

130 if (model.isDE() || competition.getType() != Geometry.Type.SQUARE_NEUMANN_2ND)

131

\vdots

These last considerations complete the implementation of the ODE model for $2\times 2$ games.

Exercise 13.12:

Plot the demise of cooperators ( $A$ types) in the donation game with a cost-to-benefit ratio of $r=0.1$ in a population with initially $99.9\%$ $A$ ’s using the ODE model and a time increment of $2$ generations.
Hint Note, after $156$ generations all cooperators have disappeared and their frequencies converged to zero. Of course, in principle, extinction would take forever with infinite precision. However, given the deterministic nature of this model you should observe the exact same outcome. Traits – Mean
Solution Set the command line options to: --module TxT --model ODE --timestep 2 --init frequency 0.999,0.001 --paymatrix 1,0;1.1,0.1 and click Apply.

.

Note: Not all the visualizations that we implemented before are available for the ODE model. In particular, when loading the ODE model with the options --model ODE, --model RK5, or --model EM only Traits – Mean, Fitness – Mean and Histogram – Fitness and, of course, the Console log are available.

.

13.10 Adding SDE model:
stochastic differential equations

Stochastic differential equations, SDE’s, establish an insightful link between the stochastic processes in individual based simulations, IBS, and the deterministic dynamics of ordinary differential equations, ODE’s. For increasing population sizes the effects of demographic fluctuations gradually decrease and eventually disappear in the deterministic limit of infinite populations (see section 2.8: From finite to infinite populations & section 5.5: From finite to infinite populations). At the same time, running IBS models becomes increasingly expensive for increasing population sizes in terms of CPU time. For large but finite populations SDE’s provide an effective and numerically efficient alternative to IBS models. In particular, the performance of SDE’s remains unaffected by the population size and hence the computational cost is constant.

The implementation of SDE’s in EvoLudo is based on the Euler-Maruyama method, which is a simple and efficient numerical scheme to solve SDE’s. The Euler-Maruyama method is a generalization of the Euler method for ordinary differential equations to stochastic differential equations and simply adds a stochastic term arising from demographic fluctuations in finite populations.

In particular, more sophisticated ODE solvers, like Runge-Kutta methods (XXX), may, for example, adjust the increments of the independent variable to allow for bigger steps if changes are small. Such optimizations are often based on evaluating the rates of change at different points in time. While this is not a problem for a deterministic system, it is not possible for SDE’s because the stochastic term is not replicable and not known in advance.

In order to add the SDE model to our EvoLudo module, all we need to do is implement another interface, HasSDE.DPairs:

26 public class TxT extends Discrete implements Payoffs,

27 HasIBS.DPairs, HasDE.DPairs, HasDE.EM, HasDE.RK5, HasDE.SDE,

28 HasPop2D.Traits, HasPop2D.Fitness, HasPop3D.Traits, HasPop3D.Fitness,

29

\vdots

.

Important: Implementing the interface HasDE.SDE also requires the implementation of either HasDE.DPairs or HasDE.DGroups for pairwise or group interactions, respectively. This requirement can only be enforced at run time without overcomplicating the interface structure of EvoLudo.

.

Since we have already implemented HasDE.DPairs in the previous section (see section 13.9: Adding ODE model: ordinary differential equations), nothing else is required. However, in case you skipped the ODE section, you need to go back now before continuing and taking the SDE model for a test drive with the exercise.

.

Note: In contrast to the ODE model, declaring a dependent trait is now essential. If no dependent trait is declared the demographic noise may result in invalid states outside the simplex!

.

Problem:

Investigate the dynamics of cooperation in the snowdrift game with a cost-to-benefit ratio of $r=0.5$ using the payoff matrix

$\displaystyle\bordermatrix{&A&B\cr A&1-r/2&1-r\cr B&1&0}.$ (13.3)
Exercise 13.13:

Calculate the equilibrium frequency of $A$ types analytically.
Solution Solve $x(1-r/2)+(1-x)(1-r)=x$ for $x$ to obtain $x=(1-r)/(1-r/2)$ or $x=2/3$ for $r=0.5$ .
Exercise 13.14:

Determine the stationary distribution after $10^{\prime}000$ generations in a population of size $N=10^{\prime}000$ with initially $50\%$ $A$ types using an SDE model and the default time increment of $1$ generation.
Hint The stationary distribution is roughly a normal distribution with mean around $67\%$ $A$ types and $33\%$ $B$ types. Statistics – Stationary distribution
Solution Set the command line options to: --module TxT --view 4 --model SDE --timestep 1 --init frequency 0.5,0.5 --paymatrix 0.75,0.5;1,0 --mutation 0.001 --timestop 10000 and click Apply and Start.
Exercise 13.15:

For comparison, repeat the simulation with an IBS model and note the difference in speed.
Solution Set the command line options to: --module TxT --view Statistics_-_Stationary_distribution --model IBS --timestep 1 --init frequency 0.5,0.5 --paymatrix 0.75,0.5;1,0 --mutation 0.001 --timestop 10000 and click Apply and Start.

13.11 Adding PDE model:
partial differential equations

Partial differential equations, PDE, provide an alternative to individual based simulations, IBS, for exploring the effect of spatial extensions on the evolutionary process. However, in contrast to the explicit modelling of population structures in IBS models, PDE models track the frequency distribution of traits in continuous space. The implementation of PDE models in EvoLudo is based on discretizing the spatial domain into a grid of discrete elements. Each element is updated following a reaction-diffusion process, or, in an evolutionary context maybe better termed a selection-diffusion process. Trait frequencies in each element change following a locally well-mixed dynamic in the selection step. During the diffusion step traits then diffuse between neighbouring elements depending on the differences in trait frequencies between adjacent elements as well as the diffusivity of each trait. Effectively this transforms the PDE into a system of coupled ordinary differential equations, where the coupling arises from the diffusion between the discrete elements.

EvoLudo supports two kinds of PDE models, again in the form of interfaces:

HasDE.PDERD:: reaction-diffusion model.
HasDE.PDEADV:: reaction-diffusion-advection model.

Just as for ODE and SDE models, both interfaces additionally require implementing either the HasDE.DPairs or HasDE.DGroups interface for pairwise or group interactions, respectively, which we completed already (see section 13.9: Adding ODE model: ordinary differential equations). Thus, the PDE models are simply added by implementing the two interfaces:

27 public class TxT extends Discrete implements Payoffs,

28 HasIBS.DPairs, HasDE.DPairs, HasDE.EM, HasDE.RK5, HasDE.SDE, HasDE.PDERD, HasDE.PDEADV,

29 HasPop2D.Traits, HasPop2D.Fitness, HasPop3D.Traits, HasPop3D.Fitness,

30

\vdots

.

Note: The models are selected by adding the option --model PDERD or --model PDEADV to the command line, respectively. In addition, the generic option --model PDE is also available, which defaults to PDERD.

.

Unfortunately, for $2\times 2$ games the results based on PDE models are a bit underwhelming. At first glance this may come as a surprise given that spatial structures are instrumental to enable cooperators to persist in the donation game (see e.g. Exercise 13.1 or Exercise 13.5) and can develop rich and complex dynamics as exemplified by evolutionary kaleidoscopes. Diffusion results in non-vanishing frequencies of both traits everywhere. As a consequence the spatial dimensions have only a transient impact on the dynamics. In the long term, the outcome is the same as in the well-mixed case. More specifically, in the donation game the spatial dimensions fail to protect cooperators from exploitation, because in continuous space diffusion undermines clustering of cooperators to offset losses against defectors. In analytical terms, the PDE for the donation game is

\displaystyle\partial_{t}u=D\nabla^{2}u+u(f_{C}-\bar{f})=D\nabla^{2}u-cu(1-u),

(13.4)

where $u$ is the frequency of cooperators at time $t$ and location $(x,y)$ . The corresponding fitness of cooperators is $f_{C}$ and $\bar{f}$ denotes the average fitness in the population. The diffusion constant $D$ specifies the migration rate of cooperators and $\nabla^{2}$ the diffusion operator. Finally, $c$ denotes the cost of cooperation. For $c>0$ the only steady state solution is $u=0$ , which is globally stable.

Exercise 13.16:

Returning to Exercise 13.1, draw the frequency distribution of $A$ and $B$ types using a PDE model for a homogeneous initial population of $A$ ’s with a single perturbation of only $B$ ’s in the centre element of the default $71\times 71$ elements. At what time does the model stop for the donation game with a cost-to-benefit ratio of $r=0.1$ and Fermi updating with noise $K=0.1$ ?
Hint After $81.1$ generations the PDE integration stops because cooperators have effectively disappeared and their frequencies converged to zero. Of course, in principle, extinction takes forever with infinite precision. However, given the deterministic nature of this model you should observe the exact same outcome. Traits – 2D Structure
. . Note: Even though the frequency distribution as shown does not appear to be homogeneous, it is actually easy to verify that the frequencies of $A$ ’s are zero everywhere by hovering over the elements in EvoLudo. It only appears as if the density of cooperators remains high in the corners because the colour gradient is automatically adjusted such that the lowest frequency of cooperators is red and the highest is blue. This helps to visually emphasize and amplify spatial patterns. However, in the present case this ‘feature’ backfires because the actual differences are negligibly small and controlled by the numerical precision (see option --accuracy). . Solution Set command line options to: --module TxT --model PDE --timestep 1 --init perturbation 1.0,0.0 --paymatrix 1,0;1.1,0.1 --playerupdate thermal 0.1 and click Apply followed by Start. . . Note: Alternatively, add the option --run to the command line to immediately launch the simulations. .

One possibility to avert this outcome, but keep the continuous spatial dimensions, is to explicitly consider the population dynamics and independently model the cooperator and defector densities rather than just their frequencies. Even though this may seem like an innocent change, the fact that birth and death events need to be explicitly modelled to allow for variable population densities requires a fundamentally different approach and a new module. Maybe an idea to write your own EvoLudo module?

13.12 Simulations using the command line

Without any further modifications the newly completed TxT module is now also available to be run from the command line as a classical Java application. For an introduction on running simulations from the command line see section 12.2.2: Command line. This avoids any security restrictions imposed on JavaScript running in a browser. Moreover, the command line allows running multiple simulations in parallel and automating the process of data collection and analysis.

For a test drive, first quit the Google Chrome browser window. This also ends the GWT EvoLudo (localhost) debug target in VS Code. Now switch to the terminal and execute ’mvn clean install’ in the top EvoLudo directory. This compiles the Java code and creates a self-contained, executable jar-archive in the ./EvoLudoJRE/target directory. It can be copied to any location and executed with ’java -jar EvoLudo.<version>.jar --module TxT --help’, which loads the TxT module and displays all the available options.

The most important, but rather technical, difference between running JavaScript simulations in the browser and a java application in the terminal, is that JavaScript is single threaded, whereas java provides a multithreaded environment. For example, in order to maintain a responsive user interface, JavaScript requires elaborate scheduling of short running tasks to prevent the browser from freezing. In contrast, in java the user interface and the simulations run in separate threads.

.

Important: The semicolon, ’;’, used to separate rows of the payoff matrix, for example, can cause troubles on the command line because it is a special character and prematurely ends the option string, which results parsing failures. Either the ’;’ must be escaped with a backslash, ’\;’, or the entire option string needs to be protected by enclosing it in quotes, ".

.

13.12.1 The --data option

The option --data is only available when running simulations from the command line. It requests and collects specific data from simulations. For a list of all available data types see section 11.2.5: Simulation options. If present, the jar-archive runs in headless mode (without a GUI) and returns the requested data in a structured format. This is particularly useful for running extended simulations or batch processing for further data analysis.

Exercise 13.17:

Repeat Exercise 13.5 but now using the jar-archive and the --data option to obtain the fixation probabilities.

. . Note: Don’t forget to remove the --view option. This triggers launching the java GUI instead of running the simulations on the command line. .

Compare the results to the histogram obtained in Exercise 13.5. When using the --seed 0 option in both cases the JavaScript and java results are identical!
Solution Run the command: ’java -jar EvoLudo.<version>.jar --module TxT --popsize 200 --geometry R4 --popupdate dB --init mutant 0,1 --paymatrix 1,0;1.1,0.1 --fitnessmap convex 1,0.01 --samples 10000 --data statprob’.

13.12.2 The --export option

Another option that is only available when running simulations from the command line is --export <filename>. Once the simulation finishes, it writes the final configuration to a file named <filename>. The exported configuration is stored as a PLIST file, which is a simple, human-readable, and machine-parsable format that can be easily imported into the EvoLudo GUI for visualization or used to restore the configuration using the --restore <filename> option and resume execution, see section 11.1.6: Export & Restore.

Exercise 13.18:

The extinction of cooperators on lattices exhibits a critical phase transition, see section LABEL:sect:XXX: LABEL:sect:XXX. This requires long relaxation times for the population to reach its dynamical equilibrium, where cooperators form isolated clusters to reduce exploitation by defectors. Consider the following string of options:
--module TxT --accuscores --geometry n --init frequency 8,2
--interactions all --paymatrix 1,0;1.037,0.037
--playerupdate thermal 0.4 --popsize 240x --popupdate async
--references random 1 --resetscores ephemeral --timestop 10000
--seed 0 --export dgr037K4.plist
What scenario does this describe?
Solution The spatial donation game with a cost-to-benefit ratio of $r=0.037$ on a $240\times 240$ square lattice with von Neumann neighbourhood and asynchronous updating. Each individual accumulates its payoff from interactions with all four nearest neighbours and probabilistically compares its payoff to a randomly chosen neighbour according to the Fermi update with a noise term of $K=0.4$ . The payoffs are ephemeral and solely used for the probabilistic comparison and discarded afterwards. The simulation stops after $10^{\prime}000$ generations and the final configuration is written to a file named dgr037K4.plist.

Exercise 13.19:

Use the EvoLudo.<version>.jar to run the simulation in Exercise 13.18 in java and generate a PLIST-file describing the relaxed configuration. Now visualize the configuration in the browser by simply dragging the generated PLIST-file onto the EvoLudo GUI in a browser window.

.

Note: To resume the simulation in the browser, simply click Start. For comparison, it is also possible to click Reset and redo the simulations in the browser. If the option string includes --seed 0 the results are identical.

.

Solution

13.12.3 The panic file

Relaxing the spatial configuration in Exercise 13.18 requires some time to finish, but that was only the preparation for the actual measurements, such as determining the size of fluctuations in equilibrium. Typically, this requires even longer execution times. All together this provides opportunities for something to go wrong – whether a power failure, a system crash, or simply killing the running process by mistake, for example installing a software update that requires rebooting the system. In an attempt to minimize the risk of data loss, all running simulations try to save their current state to a file named panic-t<time>-<date>.plist, where <time> is replaced by the current simulation time of the running model and <date> indicates the date and time when the crash happened.

.

Important: The panic file is just a last resort and there is no guarantee that saving the configuration succeeds or restoring the configuration is possible.

.

Exercise 13.20:

In order to test the panic recovery, run the simulation in Exercise 13.18 but change the name of the export file to dgr037K4a.plist. Kill the process after a few second, e.g. with Ctrl-C, and check if a file named panic-t<time>-<date>.plist has been created in the current working directory. Now resume the simulation by executing ’java -jar EvoLudo.<version>.jar --restore panic-t<time>-<date>.plist’. Let the simulation finish and compare the two final configurations dgr037K4.plist and dgr037K4a.plist. Apart from time stamps they are identical.

13.12.4 The java GUI

If neither the --data nor the --export options are present, the jar-archive launches the java GUI.

.

Note: The java GUI is not as sophisticated as the JavaScript GUI. The java GUI is still maintained but no longer under active development. Consequently, it provides only a subset of the visualizations and interactive capabilities.

.

Exercise 13.21:

For PDE models the ability to run multiple processes in parallel provides a significant speed boost. In the java implementation of PDE models both the reaction and diffusion steps are parallelized. This performance difference can be nicely illustrated by returning to Exercise 13.16. However, instead of the default $71\times 71$ elements, use a grid of $241\times 241$ elements. Open EvoLudo in the browser (as in Exercise 13.16). At the same time run the java application in the terminal with the exact same options. This launches another copy of EvoLudo but with a (more rudimentary) java GUI. Now click Step or Start in the browser as well as the java GUI and compare the impressive difference in speed of the two simulation environments.

13.13 Custom simulations

This final section covers the steps for running custom simulations with the TxT module. In order to understand an evolutionary process in greater detail it is often necessary to run simulations for a range of parameters. In principle, this is possible by manually changing the parameter of interest in the EvoLudo GUI and recording the results. However, this is not only tedious and error-prone, but also suffers from the slower execution speed of JavaScript on top of it.

A more efficient way is to automate the process through a batch file, which executes EvoLudo.<version>.jar multiple times with appropriate parameters, see section 13.12: Simulations using the command line. This allows to run simulations in the background and even to run multiple simulations in parallel. The results can be written to one or several output files. This output can then be processed for visualization or further analysis. However, for more complex investigations, custom simulations may be required.

13.13.1 Verifying the $\bm{b>c\,k\,}$ - rule

Our goal here is to create a custom simulation that calculates the fixation probability of a single cooperator in the donation game for a range of cost-to-benefit ratios $r$ , see Eq. (13.1). This requires running a series of simulations for different values of $r$ and collecting the fixation probabilities. In principle, this can be done by either method outlined above, but also serves as an illustrative example on how to build custom simulations. At the end we will have an executable jar-archive that can be run from the command line and, in addition to the parameters of TxT, provides further parameters to control the simulations. In the end this allows to validate the $b>c\,k$ -rule in structured populations (Ohtsuki et al., 2006).

All custom simulations reside in the maven module EvoLudoSims folder. Several sample files can be found in the EvoLudoSims/src/main/java/org/evoludo/simulator/exec folder, including one for $2\times 2$ games, simTBT.java but let us start from scratch.

13.13.2 Create the simulation module

To get started, create a new file simDG.java in the EvoLudoSims/src/main/java/org/evoludo/simulator/exec folder with the following code:

1 package org.evoludo.simulator.exec;

2

3 import org.evoludo.simulator.modules.TxT;

4 import org.evoludo.simulator.EvoLudo;

5

6 /**

7 * Custom EvoLudo simulations to reproduce the {@code b>c*k} rule for the

8 * donation game.

9 */

10 public class simDG extends TxT {

11

12 /**

13 * Create a new module for custom simulations for the {@code TxT} model.

14 *

15 * @param engine the pacemaker for managing the module and running the model

16 */

17 public simDG(EvoLudo evoLudo) {

18 super(evoLudo);

19 }

20 }

.

Note: Since the simDG class extends the TxT class, it simply represents another module that can be loaded by the EvoLudo engine and has access to all features of TxT.

.

An executable jar-archive requires the implementation of an entry method: public static void main(String[] args), where args refers to the array of options provided on the command line. This method is called when the jar-archive is executed. The key tasks of the main(...) method are to instantiate the simDG class and start the EvoLudo engine. Custom simulations are handled by the method custom(Module module, String[] args) in the EvoLudo engine:

21 /**

22 * Main entry point for running the custom simulation.

23 *

24 * @param args the command line arguments

25 */

26 public static void main(String[] args) {

27 EvoLudoJRE engine = new EvoLudoJRE(false);

28 engine.custom(new simDG(engine), args);

29 }

with

5 import org.evoludo.simulator.EvoLudoJRE;

In contrast to re-creating the TBT-module there is no reference file available for simDG and hence a few comments are included that help explain the logic of the code. The custom(...) method completes a number of important tasks:

1.

registers itself as a supplier of custom command line options (by overriding the collectCLO(CLOParser parser) method);
2.

loads the module and processes all command line options, which includes loading a model;
3.

resets the module, which also initializes the model and sets the seed of the random number generator (if requested);

4.

registers simDG as a MilestoneListener to monitor the progress of the simulation through milestone events, see section 12.10.1: Milestone listeners;

.

Important: Some simulations may require notifications of every update through ChangeListener events. However, this easily gets computationally expensive and risks slowing down the simulations. For this reason, the ChangeListener is not registered by default. It only gets registered if the module implements the ChangeListener interface, see section 12.10.2: Change listeners.

.

5.

returns control to simDG by calling the run() method of the parental Module class, which is inherited via the TxT and Discrete super classes.

By overriding the run() method we have full control over the simulation process. However, let us first finish the administrative part and build the custom simulation, albeit it’s just an empty shell that does nothing.

13.13.3 Building an executable jar-archive

A number of executable jar-archives are built when running ’mvn clean install’ in the top EvoLudo directory and are located in the EvoLudoSims/target directory. The maven instructions to build the jar-archives are in EvoLudoSims/pom.xml. Open the file and look for the configuration of the maven-assembly-plugin. This includes a number of <executions>. Each one builds one executable jar-archive. Simply copy and paste one <execution> and adjust the <id>, <finalName>, and <mainClass> entries to build the simDG module. The relevant snippet of the pom.xml file is:

⋮

<id>simDG</id>

<finalName>simDG.${evoludo.commit}</finalName>

<mainClass>org.evoludo.simulator.exec.simDG</mainClass>

</manifest>

</archive>

<descriptorRef>jar-with-dependencies</descriptorRef>

</descriptorRefs>

</configuration>

<phase>package</phase>

<goals>

<goal>single</goal>

</goals>

</execution>

⋮

</executions>

⋮

with the modifications highlighted. During development with VS Code, it is recommended to create a new launch configuration. The launch.json configuration file can be opened by selecting the Run & Debug tab and clicking on the settings button. The following lines add a new target for testing the DG custom simulations.

{

"name": "DG custom simulations",

"type": "java",

"request": "launch",

"mainClass": "${workspaceFolder}/EvoLudoSims/src/main/java/org/evoludo/simulator/exec/simDG.java",

"args": "--help"

},

In order to set or change the command line options, simply edit the string "args" field in the target configuration. Now breakpoints can be set and when running the target, execution is interrupted at the first specified point with all the joys and sorrows of step-by-step debugging.

Alternatively, or in order to run the custom simulations in production, run ’mvn clean install’ in the top EvoLudo directory for a full build of the EvoLudo project, which includes EvoLudoSims. In order to restrict the build to the required maven modules, run

’mvn install -pl EvoLudoCore,EvoLudoJRE,EvoLudoSims’

to build the jar-archives. The executable simDG.<version>.jar-archive is located in the EvoLudoSims/target directory. Try it out by running ’java -jar EvoLudoSims/target/simDG.<version>.jar’. This loads the custom simulations and exits. Add --help to display all available options.

.

Note: The --module option is not available for custom simulations and is ignored if provided.

.

13.13.4 Overriding run()

Now we are ready to start implementing the method run() and gaining full control over the simulation process. The skeletal setup is to request statistics mode, print the parameter settings, loop over the range of cost-to-benefit ratios $r$ and to collect statistical samples for each value of $r$ . Eventually we want to add a command line option to set the range of the cost-to-benefit ratio $r$ as well as the number of steps. However, for now let’s just add a double array cbRatio with [start, end, steps]:

9 /**

10 * The array for the range and number of steps of the cost-to-benefit ratio.

11 */

12 double[] cbRatio = new double[] { 0.0, 1.0, 5.0 };

and the structure of the run() method:

16 @Override

17 public void run() {

18 // set statistics mode

19 model.requestMode(Mode.STATISTICS_SAMPLE);

20 // print simulation header with parameters

21 engine.dumpParameters();

22 double r = cbRatio[0];

23 double end = cbRatio[1];

24 int steps = (int) cbRatio[2];

25 double incr = (end - r) / steps;

26 double payoffs[][] = new double[2][2];

27 // the payoffs to A types is independent of r

28 payoffs[A][A] = 1.0;

29 payoffs[A][B] = 0.0;

30 // print data header

31 // TODO

32 for (int s = 0; s <= steps; s++) {

33 // update the payoffs to B types

34 payoffs[B][A] = 1.0 + r;

35 payoffs[B][B] = r;

36 setPayoffs(payoffs);

37 // collect statistics samples

38 // TODO

39 // report statistics for current r

40 // TODO

41 // update r

42 r += incr;

43 }

44 engine.dumpEnd();

45 }

The custom simulation compiles and runs fine but, apart from printing the parameters, nothing much happens. Next, we set everything up for the statistics and print the results for each value of $r$ . To simplify printing output to the console (or any output file that is specified by --output or --append) we introduce the convenience variables out. For dealing with statistics, another field jrengine will be useful later. We declare the two fields at the beginning of the class:

12 /**

13 * The output stream. Defaults to {@code System.out}.

14 */

15 PrintStream out;

16

17 /**

18 * The EvoLudoJRE engine for running the simulation. This is a convenience field

19 * that saves us casting engine to EvoLudoJRE every time we need to access its

20 * methods.

21 */

22 EvoLudoJRE jrengine;

and initialize them in the constructor:

33 public simDG(EvoLudoJRE engine) {

34 super(engine);

35 out = engine.getOutput();

36 jrengine = engine;

37 }

For the statistics, all we need is an array to store the number of times that A or B types have fixed. With these additions we can provide feedback to the user:

52 double[] nFix = new double[2];

53 // print data header

54 out.printf("# r\trhoA\trhoB\n");

62 // report statistics for current r

63 ArrayMath.normalize(nFix);

64 out.printf("%.3f\t%.6f\t%.6f\n", r, nFix[A], nFix[B]);

with

7 import org.evoludo.math.ArrayMath;

The normalize(double[]) conveniently converts the number of fixations into probabilities. Compiling and running the custom module now returns the following output:

EvoLudo % java -jar EvoLudoSims/target/simDG.<version>.jar

\vdots

# data:

# r rhoA rhoB

0.000 NaN NaN

0.200 NaN NaN

0.400 NaN NaN

0.600 NaN NaN

0.800 NaN NaN

1.000 NaN NaN

# runningtime: 0:00:00.53

The NaN values are of no importance and simply a consequence of all elements of nFix being zero. However, at least we confirmed that our loop is working as expected. The key remaining piece is to implement the actual collection of statistical samples. The number of samples is set by the --samples option from the Model class. The method generateSample() in the class EvoLudoJRE generates and returns one statistics sample with the results stored in a FixationData object:

61 // collect statistics samples

62 for (int i = 0; i < model.getNSamples(); i++) {

63 FixationData fixData = jrengine.generateSample();

64 nFix[fixData.typeFixed]++;

65 }

Now let’s see what we get with $1^{\prime}000$ samples for each $r$ :

EvoLudo % java -jar EvoLudoSims/target/simDG.<version>.jar --samples 1000

\vdots

# r rhoA rhoB

Exception in thread "main" java.lang.NullPointerException: Cannot read field "mutantNode" because "fix" is null

at org.evoludo.simulator.EvoLudoJRE.generateSample(EvoLudoJRE.java:746)

at org.evoludo.simulator.exec.simDG.run(simDG.java:63)

at org.evoludo.simulator.EvoLudoJRE.custom(EvoLudoJRE.java:349)

at org.evoludo.simulator.exec.simDG.main(simDG.java:82)

Oops! Recall that fixation statistics are only available for appropriate initial configurations. Let’s try again and add the --init mutant 0,1 to place a single A type in a homogeneous B population. Also, don’t forget to set the Moran death-birth update, --popupdate dB, together with the convex payoff to fitness mapping, --fitnessmap convex 1,0.01 with baseline fitness, $1$ , and weak selection, $0.01$ . The output now reads:

EvoLudo % java -jar EvoLudoSims/target/simDG.<version>.jar --popupdate dB --init mutant 0,1 --fitnessmap convex 1,0.01 --samples 1000

\vdots

# r rhoA rhoB

0.000 0.012000 0.988000

0.200 0.009003 0.990997

0.400 0.005004 0.994996

0.600 0.008996 0.991004

0.800 0.003006 0.996994

1.000 0.007995 0.992005

# runningtime: 0:00:00.247

Looks rather promising – possibly with some minor variation due to the random initialization.

.

Note: In order to obtain reproducible simulations, it is essential to set the seed of the random number generator using the --seed option. This ensures that the same sequence of random numbers is generated for each simulation. Every time the model is reset, the random number generator is reset as well, which not only results in identical initial configurations, but also re-generates the same network (if any) for the population structure.

.

Important: For statistical evaluations it is imperative to clear the seed after the first reset, which has already happened before returning control to simDG through calling the method run(). Otherwise, the same data point is generated over and over and over…

.

If the option --seed has been set, clear it before starting the statistics:

46 // print simulation header with parameters

47 engine.dumpParameters();

48 if (engine.cloSeed.isSet()) {

49 // RNG seed is set. now clear seed to obtain reproducible statistics

50 // rather repeatedly the same data point

51 engine.getRNG().clearRNGSeed();

52 }

.

Note: Clear the seed only after printing the parameters to the console to ensure that the seed settings are correctly reported.

.

Earlier, when using the command model.requestMode(Mode.STATISTICS_SAMPLE); we tacitly assumed that the model actually supports statistics mode. While this is true for IBS and SDE models it does not extend to ODE or PDE models. There are different approaches to handle this situation:

1.

check the model type (for example with model.isIBS() or model instanceof IBS etc. at the beginning of the run() method and exit if not);
2.

check the return value of model.requestMode(Mode.STATISTICS_SAMPLE); and exit if false is returned (indicates that the model does not support statistics mode);
3.

modify the --model option to only allow compatible models;
4.

restrict the available model types for simDG.

The last option is clearly the most elegant solution – and the easiest to implement:

47 @Override

48 public Type[] getModelTypes() {

49 return new Type[] { Type.IBS, Type.SDE };

50 }

13.13.5 Command line options

The range of the cost-to-benefit ratio $r$ is at present hard-coded to $[0,1]$ with $5$ steps. Instead of setting these default values when initializing the field cbRatio it would be more convenient to adjust the settings with a command line option, for example, --ratio <s[,e[,n]]>, where the array of arguments follows the same convention as the cbRatio, only that the endpoint, <e>, and the number of steps, <n>, are optional:

87 /**

88 * The command line option for specifying the range and increment of the

89 * cost-to-benefit ratio.

90 */

91 public final CLOption cloCBRange = new CLOption("ratio", "0,1,5", Category.Simulation,

92 "--ratio <s[,e[,n]]> range [s,e], n steps", new CLODelegate() {

93

94 @Override

95 public boolean parse(String arg) {

96 double[] rVec = CLOParser.parseVector(arg);

97 switch (rVec == null ? -1 : rVec.length) {

98 default:

99 return false;

100 case 1:

101 // run simulation for r=rVec[0]

102 rVec = new double[] { rVec[0], rVec[0], 0.0 };

103 break;

104 case 2:

105 // run simulations for r=rVec[0] and r=rVec[1]

106 rVec = new double[] { rVec[0], rVec[1], 1.0 };

107 break;

108 case 3:

109 // run simulations for r=rVec[0] to r=rVec[1] in rVec[2] steps

110 break;

111 }

112 cbRatio = rVec;

113 return true;

114 }

115 });

116

117 @Override

118 public void collectCLO(CLOParser parser) {

119 super.collectCLO(parser);

120 parser.addCLO(cloCBRange);

121 }

with

4 import org.evoludo.util.CLOParser;

5 import org.evoludo.util.CLOption;

6 import org.evoludo.util.CLOption.CLODelegate;

The initial value of the field cbRatio is no longer used and can be safely removed. If --ratio is not provided on the command line, its default argument, "0,1,5", is still processed.

When running our custom simulations we do not want to allow the user to set the payoff matrix! We can easily prevent this by removing the cloPayoffs option. After all command line options have been collected, the method adjustCLO(CLOParser) is called to make any adjustments. This is the place to remove the --paymatrix option:

130 @Override

131 public void adjustCLO(CLOParser parser) {

132 super.adjustCLO(parser);

133 // remove payoffs option

134 parser.removeCLO(cloPayoffs);

135 }

Now we are ready to verify the $b>c\,k$ -rule.

.

Note: When running statistical evaluations using structures that involve random elements, such as random (regular) graphs, the structure is regenerated for every new sample. Computationally this is a fairly expensive and unnecessary at the same time. In a structured population of size

N

, there are

N

possible configurations for the placement of the initial mutant. Thus, it is safe to use the same structure for at least

N

samples. Even if the same initial configuration is used multiple times by chance, the bias is negligible as long as the evolutionary trajectory is stochastic, i.e. the initial condition does not uniquely determine the end state. The interval for resetting the structure can be controlled with the option --statistics.

.

Exercise 13.22:

Calculate the fixation probabilities on random graphs with $\bar{k}=4$ and $N=200$ , scanning the interval $r\in[0,0.5]$ using $10$ steps and $10^{\prime}000$ samples.
Solution The resulting output is:

EvoLudo % java -jar EvoLudoSims/target/simDG.<version>.jar
--popsize 200 --popupdate dB --init mutant 0,1
--fitnessmap convex 1,0.01 --geometry R4 --ratio 0,0.5,10 --seed 0
--statistics reset 200 --samples 100000
$\vdots$
# data:
# r rhoA rhoB
0.000 0.006020 0.993980
0.050 0.005990 0.994010
0.100 0.005730 0.994270
0.150 0.005390 0.994610
0.200 0.004950 0.995050
0.250 0.005150 0.994850
0.300 0.005000 0.995000
0.350 0.004380 0.995620
0.400 0.004520 0.995480
0.450 0.004240 0.995760
0.500 0.004090 0.995910
# runningtime: 0:01:58.90

With the --seed option the numbers should be identical. Without the --statistics option the results would show insignificant differences, but with a longer execution time. In this case $\approx 25$ seconds longer, which translates to a performance improvement of roughly $20\%$ when using --statistics. More importantly, overall the trend is clear, but the results remain noisy even at $10^{\prime}000$ samples. According to the $b>c\,k$ -rule the threshold for $\varphi_{A}>1/N=0.005$ should be near $1/\bar{k}=0.25$ .

Exercise 13.23:

Equip yourself with some patience and repeat the above calculations for $10^{6}$ samples, while also limiting the range to $r\in[0.15,0.35]$ .

.

Note: With

10^{6}

samples, the statistical error is of the order

1/\sqrt{10^{6}}=10^{-3}

such that we end up with merely two significant digits. In order to add another significant digit,

10^{8}

samples are required, and the running time increase

100

-fold.

.

Solution

Long-running simulations that get accidentally interrupted attempt to dump their final state to a file named panic-t<t>-<date>.plist, where <t> is replaced by the evolutionary time at which the model stopped, and <date> by the actual time when the crash happened, see section 13.12: Simulations using the command line. Unfortunately, this is only of limited use for statistical evaluations. A better approach to reduce the risk of data loss, is to report the current fixation probability in regular intervals together with the number of samples collected thus far. Possibly add a leading # character to simplify post-processing of the output.

Chapter 13 Step-by-step tutorial 𝑰: Re-creating the TBT module

13.1 Create the module

13.1.1 Adding the module to EvoLudo

13.2 Running (and debugging) the module

13.3 Implementing the module

13.3.1 Loading and unloading of the module

13.4 Module specific fields

13.4.1 Adding command-line options

Using the constructor CLOption(…)

Creating the delegate

13.4.2 Tying up loose ends

13.5 Adding IBS model: Individual based simulations

13.6 Adding GUI views

13.6.1 Traits and fitness in 2D and 3D

13.6.2 Mean traits and fitness

13.6.3 Fitness histogram

13.6.4 Degree histogram

13.7 Adding statistics modes

13.8 Advanced topics

13.8.1 Implementing evolutionary kaleidoscopes

Customizing IBS populations

Customizing initial configurations

13.8.2 Adding square lattice geometry with second neighbours

Customizing trait names for graphics and tooltips

Customizing the status bar

Customizing the views for mean trait frequencies and fitness

13.9 Adding ODE model: ordinary differential equations

13.10 Adding SDE model: stochastic differential equations

13.11 Adding PDE model: partial differential equations

13.12 Simulations using the command line

13.12.1 The --data option

13.12.2 The --export option

13.12.3 The panic file

13.12.4 The java GUI

13.13 Custom simulations

13.13.1 Verifying the 𝒃>𝒄⁢𝒌 - rule

13.13.2 Create the simulation module

13.13.3 Building an executable jar-archive

13.13.4 Overriding run()

13.13.5 Command line options

References

Chapter 13 Step-by-step tutorial $\bm{I}$ :
Re-creating the TBT module

13.9 Adding ODE model:
ordinary differential equations

13.10 Adding SDE model:
stochastic differential equations

13.11 Adding PDE model:
partial differential equations

13.13.1 Verifying the $\bm{b>c\,k\,}$ - rule