GLM Basis Functions

[1]:

import cedalion
import cedalion.models.glm.basis_functions as bf
import cedalion.models.glm.design_matrix as dm

import cedalion.dataclasses as cdc

import matplotlib.pyplot as p
import numpy as np
import xarray as xr
import matplotlib.pyplot as p
import cedalion.datasets

units = cedalion.units

xr.set_options(display_expand_data=False)

[1]:

<xarray.core.options.set_options at 0x7fdf0470ebd0>

[2]:

# dummy time series
fs = 8.0
ts = cdc.build_timeseries(
    np.random.random((100, 1, 2)),
    dims=["time", "channel", "chromo"],
    time=np.arange(100) / fs,
    channel=["S1D1"],
    value_units=units.uM,
    time_units=units.s,
    other_coords={'chromo' : ["HbO", "HbR"]}
)
display(ts)

<xarray.DataArray (time: 100, channel: 1, chromo: 2)> Size: 2kB
[µM] 0.5818 0.4453 0.458 0.5398 0.9459 ... 0.06786 0.2011 0.9574 0.4412 0.6611
Coordinates:
  * time     (time) float64 800B 0.0 0.125 0.25 0.375 ... 12.0 12.12 12.25 12.38
    samples  (time) int64 800B 0 1 2 3 4 5 6 7 8 ... 91 92 93 94 95 96 97 98 99
  * channel  (channel) <U4 16B 'S1D1'
  * chromo   (chromo) <U3 24B 'HbO' 'HbR'
[3]:

basis = bf.GaussianKernels(
    t_pre=5 * units.s,
    t_post=30 * units.s,
    t_delta=3 * units.s,
    t_std=3 * units.s,
)
hrf = basis(ts)

p.figure()
for i_comp, comp in enumerate(hrf.component.values):
    p.plot(hrf.time, hrf[:, i_comp], label=comp)

p.axvline(-5, c="r", ls=":")
p.axvline(30, c="r", ls=":")
p.legend(ncols=3)

[3]:

<matplotlib.legend.Legend at 0x7fdf04757490>
../../_images/examples_modeling_31_glm_basis_functions_3_1.png 
[4]:

basis = bf.GaussianKernelsWithTails(
    t_pre=5 * units.s,
    t_post=30 * units.s,
    t_delta=3 * units.s,
    t_std=3 * units.s,
)
hrf = basis(ts)

p.figure()
for i_comp, comp in enumerate(hrf.component.values):
    p.plot(hrf.time, hrf[:, i_comp], label=comp)
p.axvline(-5, c="r", ls=":")
p.axvline(30, c="r", ls=":")
p.legend(ncols=3)

[4]:

<matplotlib.legend.Legend at 0x7fde9017ba10>
../../_images/examples_modeling_31_glm_basis_functions_4_1.png 
[5]:

basis = bf.Gamma(
    tau={"HbO": 0 * units.s, "HbR": 1 * units.s},
    sigma=3 * units.s,
    T=0 * units.s,
)
hrf = basis(ts)
display(hrf)
p.figure()
for i_comp, comp in enumerate(hrf.component.values):
    for i_chromo, chromo in enumerate(hrf.chromo.values):
        p.plot(hrf.time, hrf[:, i_comp, i_chromo], label=f"{comp} {chromo}")

p.legend()

<xarray.DataArray (time: 107, component: 1, chromo: 2)> Size: 2kB
0.0 0.0 0.004711 0.0 0.01875 ... 2.533e-07 3.574e-06 1.789e-07 2.601e-06
Coordinates:
  * time       (time) float64 856B 0.0 0.125 0.25 0.375 ... 13.0 13.12 13.25
  * chromo     (chromo) <U3 24B 'HbO' 'HbR'
  * component  (component) <U5 20B 'gamma'
[5]:

<matplotlib.legend.Legend at 0x7fde8ef897d0>
../../_images/examples_modeling_31_glm_basis_functions_5_2.png 
[6]:

basis = bf.Gamma(
    tau={"HbO": 0 * units.s, "HbR": 1 * units.s},
    sigma=2 * units.s,
    T=5 * units.s,
)
hrf = basis(ts)
display(hrf)

p.figure()
for i_comp, comp in enumerate(hrf.component.values):
    for i_chromo, chromo in enumerate(hrf.chromo.values):
        p.plot(hrf.time, hrf[:, i_comp, i_chromo], label=f"{comp} {chromo}")

p.legend()

<xarray.DataArray (time: 114, component: 1, chromo: 2)> Size: 2kB
0.0 0.0 0.0005512 0.0 0.00273 ... 6.285e-09 4.08e-07 3.599e-09 2.491e-07
Coordinates:
  * time       (time) float64 912B 0.0 0.125 0.25 0.375 ... 13.88 14.0 14.12
  * chromo     (chromo) <U3 24B 'HbO' 'HbR'
  * component  (component) <U5 20B 'gamma'
[6]:

<matplotlib.legend.Legend at 0x7fde8ef8a1d0>
../../_images/examples_modeling_31_glm_basis_functions_6_2.png 
[7]:

basis = bf.GammaDeriv(
    tau=1 * units.s,
    sigma=1 * units.s,
    T=10 * units.s,
)
hrf = basis(ts)
display(hrf)
p.figure()
for i_comp, comp in enumerate(hrf.component.values):
    for i_chromo, chromo in enumerate(["HbO"]):
        p.plot(hrf.time, hrf[:, i_comp, i_chromo], label=f"{comp} {chromo}")

p.legend()

<xarray.DataArray (time: 121, component: 2, chromo: 2)> Size: 4kB
0.0 0.0 0.0 0.0 0.0 0.0 ... -7.624e-06 3.086e-07 3.086e-07 -2.978e-06 -2.978e-06
Coordinates:
  * time       (time) float64 968B 0.0 0.125 0.25 0.375 ... 14.75 14.88 15.0
  * chromo     (chromo) <U3 24B 'HbO' 'HbR'
  * component  (component) <U11 88B 'gamma' 'gamma_deriv'
[7]:

<matplotlib.legend.Legend at 0x7fde8efff550>
../../_images/examples_modeling_31_glm_basis_functions_7_2.png 
[8]:

basis = bf.AFNIGamma(
    p=1,
    q=0.7 * units.s,
    T=0 * units.s,
)
hrf = basis(ts)
display(hrf)
p.figure()
for i_comp, comp in enumerate(hrf.component.values):
    for i_chromo, chromo in enumerate(["HbO"]):
        p.plot(hrf.time, hrf[:, i_comp, i_chromo], label=f"{comp} {chromo}")

p.legend()

<xarray.DataArray (time: 31, component: 1, chromo: 2)> Size: 496B
0.0 0.0 0.407 0.407 0.6809 0.6809 ... 0.0918 0.07953 0.07953 0.06882 0.06882
Coordinates:
  * time       (time) float64 248B 0.0 0.125 0.25 0.375 ... 3.375 3.5 3.625 3.75
  * chromo     (chromo) <U3 24B 'HbO' 'HbR'
  * component  (component) <U10 40B 'afni_gamma'
[8]:

<matplotlib.legend.Legend at 0x7fde8cb2fa10>
../../_images/examples_modeling_31_glm_basis_functions_8_2.png 
[9]:

basis = bf.DiracDelta()
hrf = basis(ts)
display(hrf)
p.figure()
for i_comp, comp in enumerate(hrf.component.values):
    for i_chromo, chromo in enumerate(["HbO"]):
        p.stem(hrf.time, hrf[:, i_comp, i_chromo], label=f"{comp} {chromo}")

p.legend()

<xarray.DataArray (time: 2, component: 1, chromo: 2)> Size: 32B
1.0 1.0 0.0 0.0
Coordinates:
  * time       (time) float64 16B 0.0 0.125
  * chromo     (chromo) <U3 24B 'HbO' 'HbR'
  * component  (component) <U6 24B 'square'
[9]:

<matplotlib.legend.Legend at 0x7fde8dde4290>
../../_images/examples_modeling_31_glm_basis_functions_9_2.png