Compute BOLD online in Numba backend

[maedoc]

Describe the new feature or enhancement

When simulating BOLD with the Numba backend, one currently has to first do the simulation then process the time series with the BOLD analyzer. For long time series, it is more memory-efficient to compute the BOLD in an online fashion.

Describe your proposed implementation

import numpy as np
import numba as nb

def make_bold():
    tau_s = nb.float32(0.65)
    tau_f = nb.float32(0.41)
    tau_o = nb.float32(0.98)
    alpha = nb.float32(0.32)
    te = nb.float32(0.04)
    v0 = nb.float32(4.0)
    e0 = nb.float32(0.4)
    epsilon = nb.float32(0.5)
    nu_0 = nb.float32(40.3)
    r_0 = nb.float32(25.0)
    recip_tau_s = nb.float32(1.0 / tau_s)
    recip_tau_f = nb.float32(1.0 / tau_f)
    recip_tau_o = nb.float32(1.0 / tau_o)
    recip_alpha = nb.float32(1.0 / alpha)
    recip_e0 = nb.float32(1.0 / e0)
    k1 = nb.float32(4.3 * nu_0 * e0 * te)
    k2 = nb.float32(epsilon * r_0 * e0 * te)
    k3 = nb.float32(1.0 - epsilon)
    
    @numba.njit
    def bold_step(state, x, dt):
        s = state[0]
        f = state[1]
        v = state[2]
        q = state[3]
        ds = x - recip_tau_s * s - recip_tau_f * (f - 1)
        df = s
        dv = recip_tau_o * (f - v**recip_alpha)
        dq = recip_tau_o * (f * (1 - (1 - e0)**(1 / f)) * recip_e0 - v**recip_alpha * (q / v))
        s += dt * ds
        f += dt * df
        v += dt * dv
        q += dt * dq
        state[0] = s
        state[1] = f
        state[2] = v
        state[3] = q
        return v0 * (k1*(1 - q) + k2*(1 - q / v) + k3*(1 - v))

    return bold_step

This should be integrated with the Numba backend:

• check if BOLD is requested
• allocate extra areas for BOLD
• include calls to bold_step at each time step in the simulation, over each node & svar
• return BOLD in addition to time series

[maedoc]

Some of the BOLD derivatives are subject to loss of precision: open Herbie and enter f * (1 - pow(1 - e0, 1 / f)) * recip_e0 - pow(v,recip_alpha) * (q / v), and click v for bit error,

Their analysis suggest how to rewrite to drop that bit error from the red line to the blue line.

Fortunately the 1st ODE where x appears is linear wrt. x, so maybe these changes aren't needed.