In [1]:
%pylab nbagg
from tvb.simulator.lab import *
Populating the interactive namespace from numpy and matplotlib
   INFO  log level set to INFO

Comparing integration methods

Often it is difficult to know which is the best integration method for your simulation. TVB offers a handful of methods, with recently added variable order methods which allow one to safely choose higher time steps, as long as they remain smaller than the smallest time delay. In the following tutorial, we'll perform identical simulations, changing only the integration method in order to check what differences result.

In [2]:
base_dt = 0.1
var_order_dt = 5.0

methods = [
    (integrators.EulerDeterministic, base_dt),
    (integrators.HeunDeterministic, 2*base_dt),
    (integrators.Dop853, var_order_dt),
    (integrators.Dopri5, var_order_dt),
    (integrators.RungeKutta4thOrderDeterministic, 4*base_dt),
    (integrators.VODE, var_order_dt),
]
In [3]:
conn  = connectivity.Connectivity(load_default=True)

figure(figsize=(15, 10))
for i, (method, dt) in enumerate(methods):
    numpy.random.seed(42)
    sim = simulator.Simulator(
        connectivity=conn,
        model=models.Generic2dOscillator(a=0.1),
        coupling=coupling.Linear(a=0.0),
        integrator=method(dt=dt),
        monitors=monitors.TemporalAverage(period=5.0),
        simulation_length=1000.0,
    ).configure()
    (t, raw), = sim.run()
    
    if i == 0:
        euler_raw = raw
    else:
        if raw.shape[0] != euler_raw.shape[0]:
            continue
        raw = abs(raw - euler_raw) / euler_raw.ptp() * 100.0
    
    subplot(3, 2, i + 1)
    plot(t, raw[:, 0, :, 0], 'k', alpha=0.1)
    if i > 0:
        ylabel('% diff')
        plot(t, raw[:, 0, :, 0].mean(axis=1), 'k', linewidth=3)
    title(method._ui_name)
    grid(True)
WARNING  File 'hemispheres' not found in ZIP.

At the moment of writing the RK4 method is not producing the correct number of output samples, and we're investigating the issue. In any case, the different integrators produce qualitatively similar results for this transient.

What is the order of magnitude of these diffences? This time, we compare pair-wise:

In [4]:
conn  = connectivity.Connectivity(load_default=True)

raws = []
names = []
for i, (method, dt) in enumerate(methods):
    numpy.random.seed(42)
    sim = simulator.Simulator(
        connectivity=conn,
        model=models.Generic2dOscillator(a=0.1),
        coupling=coupling.Linear(a=0.0),
        integrator=method(dt=dt),
        monitors=monitors.TemporalAverage(period=5.0),
        simulation_length=200.0,
    ).configure()
    (t, raw), = sim.run()
    raws.append(raw)
    names.append(method._ui_name)

n_raw = len(raws)
figure(figsize=(15, 15))
for i, (raw_i, name_i) in enumerate(zip(raws, names)):
    for j, (raw_j, name_j) in enumerate(zip(raws, names)):
        subplot(n_raw, n_raw, i*n_raw + j + 1)
        if raw_i.shape[0] != t.shape[0] or raw_i.shape[0] != raw_j.shape[0]:
            continue
        if i == j:
            plot(t, raw_i[:, 0, :, 0], 'k', alpha=0.1)
        else:
            semilogy(t, (abs(raw_i - raw_j) / raw_i.ptp())[:, 0, :, 0], 'k', alpha=0.1)
            ylim(exp(r_[-30, 0]))
        
        grid(True)
        if i==0:
            title(name_j)
        if j==0:
            ylabel(name_i)
    
    if i == 0:
        euler_raw = raw
    else:
        raw = abs(raw - euler_raw) / euler_raw.ptp() * 100.0

tight_layout()
WARNING  File 'hemispheres' not found in ZIP.

Along the diagonal, the time series of the simulation are plotted. Off the diagonal, the differences between methods for the same simulation are given. We can see that the Euler & Heun methods have low agreement with the variable order, which is to be expected, while the variable order methods are in good agreement.

So far we've allowed the variable order integrators to choose their dt, bounded by a maximum of 5.0. We might ask, however, about the changes as dt grows. For example,

In [5]:
dts = [float(10**e) for e in r_[-2:0:10j]]
print dts
conn  = connectivity.Connectivity(load_default=True)

raws = []
for i, dt in enumerate(dts):
    numpy.random.seed(42)
    sim = simulator.Simulator(
        connectivity=conn,
        model=models.Generic2dOscillator(a=0.1),
        coupling=coupling.Linear(a=0.0),
        integrator=integrators.VODE(dt=dt),
        monitors=monitors.TemporalAverage(period=1.0),
        simulation_length=1200.0,
    ).configure()
    (t, raw), = sim.run()
    t = t[:1000]
    raw = raw[:1000]
    raws.append(raw)
    
figure(figsize=(10, 10))
for i, dt in enumerate(dts):
    subplot(len(dts)/3, 3+1, i + 1)
    if i == 0:
        dat = raws[i]
    else:
        dat = log10((abs(raws[i] - raws[0]) / raws[0].ptp()))
    plot(t, dat[:, 0, :, 0], 'k', alpha=0.1)
[0.01, 0.016681005372000592, 0.027825594022071243, 0.046415888336127774, 0.0774263682681127, 0.1291549665014884, 0.21544346900318834, 0.3593813663804626, 0.5994842503189409, 1.0]
WARNING  File 'hemispheres' not found in ZIP.

In conclusion, we can see that difference do exist between techniques. Variable order, especially the VODE, provide some speed up over fixed step, but because the connectivity depends on the step size, the step size should remain much lower that the smallest delays, to maintain a releastic simulation of the conduction speed.

When in doubt, compare the variable order solution to the fixed order, keeping in mind that a much lower step size will be required for the fixed step size, because such methods do not internally adapt their step size up or down according to an error criteria.