Added new HRF kernels (gamma, double exponential, mixture of gammas)

[brucsk]

Adapted from TVB (https://github.com/the-virtual-brain/tvb-root/blob/master/tvb_library/tvb/datatypes/equations.py)

[Copilot]

Pull request overview

Adds additional hemodynamic response function (HRF) kernels to the BOLD monitor module, adapting formulations from TVBSim so HRFBold can be driven by alternative, literature-based HRF shapes.

Changes:

• Added GammaHRFKernel (Boynton-style gamma HRF) with TVBSim-like normalization and amplitude scaling.
• Added DoubleExponentialHRFKernel (difference of two damped sinusoids) with TVBSim-like normalization and amplitude scaling.
• Added MixtureOfGammasHRFKernel (Glover-style mixture of gammas) HRF definition.

Comments suppressed due to low confidence (2)

src/tvboptim/observations/tvb_monitors/bold.py:189

• Normalization divides by jnp.max(kernel) without guarding against a zero/invalid maximum. For parameter choices that produce an all-zero kernel (or NaNs/Infs), this will create NaNs and propagate into the BOLD output. Consider guarding with an epsilon / conditional and/or validating parameters (e.g., tau > 0, finite max).

        # Replicate TVBSim's normalization and amplitude scaling from evaluate()
        kernel = kernel / jnp.max(kernel)
        kernel = kernel * self.a

src/tvboptim/observations/tvb_monitors/bold.py:255

• Normalizing by jnp.max(kernel) is fragile here because this kernel is oscillatory and can have negative peaks; if the maximum is <= 0 (or very close to 0), the normalization can flip signs or produce NaNs/Infs. Consider normalizing by max(abs(kernel)) (with a zero guard) or using a more explicit TVB-compatible normalization rule.

        # Replicate TVBSim's normalization + amplitude scaling from evaluate()
        kernel = kernel / jnp.max(kernel)
        kernel = kernel * self.a

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[mapi1]

Thanks for this contribution, looks great so far!

I added some basic unit test for HRF kernels which you can test against by adding your kernels to the list here, after pulling the latest changes from main into this PR.

Final remark, you can get rid of the __init__ functions, as the kernels are a equinox.Module, a JAX specific data class that populates it automatically from the class field definitions. This saves some LOC and redundant definition of defaults.

[brucsk]

Thank you very much for your feedback!

About the __init__, I had realized that later on (hence why it's not on the MixtureOfGammas) but forgot to take it off the earlier ones. Fixed it on the last commit.

About the unit tests: the new kernels failed some of them (cf screenshot below):

I'll try to work on it soon, but if you have any immediate ideas on what might be causing these failures, please let me know!

[mapi1]

Looking at the code, the first two error probably come from the line: kernel = kernel / jnp.max(kernel) which if we only evaluate at one time point and its actually zero gives a nan. The third error looks like the Exponential has not yet decayed with ~0.009, the test currently accepts 5% of the peak value, so you probably need to increase the default duration here. My assumption for the jit error is, that the kernels actually jit fine, but then some minor numerical differences around critical points cause the test to fail.

I adopted the tests accordingly, Most of your kernels should pass now.

[brucsk]

You were right about everything; I updated the normalization in Gamma and Double Exponential kernels to avoid divsion by zero and set the default duration for the Double Exponential to 40000 ms. Indeed, relaxing the jittable test tolerance solved it as well. All the tests are okay now!

[mapi1]

Very nice! All tests are also passing on CI, so I see no reason not to merge this. I will cut a new release later so you can access your changes from PyPI from version v0.2.10 onwards.

Thanks again for your contribution 🙌