.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/00_encodingdecoding/fit_prf.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_00_encodingdecoding_fit_prf.py: Different flavors of visual population receptive field models ============================================================== In this example script we will try out increasingly complex models for visual population receptive fields (PRFs). We will start with a simple Gaussian PRF model, and then add more complexity step by step. .. GENERATED FROM PYTHON SOURCE LINES 12-15 Load data --------- First we load in the data. We will use the Szinte (2024)-dataset. .. GENERATED FROM PYTHON SOURCE LINES 15-30 .. code-block:: Python from braincoder.utils.data import load_szinte2024 import matplotlib.pyplot as plt data = load_szinte2024() # This is the visual stimulus ("design matrix") paradigm = data['stimulus'] grid_coordinates = data['grid_coordinates'] # This is the fMRI response data d = data['v1_timeseries'] d.index.name = 'frame' tr = data['tr'] .. GENERATED FROM PYTHON SOURCE LINES 31-34 Simple 2D Gaussian Recetive Field model ------------------------------------- Now we set up a simple Gaussian PRF model .. GENERATED FROM PYTHON SOURCE LINES 34-39 .. code-block:: Python from braincoder.models import GaussianPRF2DWithHRF from braincoder.hrf import SPMHRFModel hrf_model = SPMHRFModel(tr=tr) model_gauss = GaussianPRF2DWithHRF(data=d, paradigm=paradigm, hrf_model=hrf_model, grid_coordinates=grid_coordinates) .. GENERATED FROM PYTHON SOURCE LINES 40-41 And a parameter fitter... .. GENERATED FROM PYTHON SOURCE LINES 41-45 .. code-block:: Python from braincoder.optimize import ParameterFitter par_fitter = ParameterFitter(model=model_gauss, data=d, paradigm=paradigm) .. GENERATED FROM PYTHON SOURCE LINES 46-48 Now we try out a relatively coarse grid search to find the some parameters to start the gradient descent from. .. GENERATED FROM PYTHON SOURCE LINES 48-67 .. code-block:: Python import numpy as np x = np.linspace(-8, 8, 10) y = np.linspace(-4, 4, 10) sd = np.linspace(0.1, 4, 10) # We start the grid search using a correlation cost, so ampltiude # and baseline do not influence those results. # We will optimize them later using OLS. baseline = [0.0] amplitude = [1.0] # Now we can do the grid search pars_gauss_grid = par_fitter.fit_grid(x, y, sd, baseline, amplitude, use_correlation_cost=True) # And refine the baseline and amplitude parameters using OLS pars_gauss_ols = par_fitter.refine_baseline_and_amplitude(pars_gauss_grid) .. rst-class:: sphx-glr-script-out .. code-block:: none 0%| | 0/1 [00:00] .. GENERATED FROM PYTHON SOURCE LINES 95-100 Fit HRFs -------- The standard canonical (SPM) HRF we use is often not a great fit to actual data. To better account for the HRF. We can optimize the HRFs per voxel. We first initialize a GaussianPRF-model with a flexible HRF. .. GENERATED FROM PYTHON SOURCE LINES 100-111 .. code-block:: Python model_hrf = GaussianPRF2DWithHRF(data=d, paradigm=paradigm, hrf_model=hrf_model, grid_coordinates=grid_coordinates, flexible_hrf_parameters=True) par_fitter_hrf = ParameterFitter(model=model_hrf, data=d, paradigm=paradigm) # We set hrf_delay and hrf_dispersion to standard values pars_gauss_gd['hrf_delay'] = 6 pars_gauss_gd['hrf_dispersion'] = 1 pars_gauss_hrf = par_fitter_hrf.fit(init_pars=pars_gauss_gd, max_n_iterations=1000) .. rst-class:: sphx-glr-script-out .. code-block:: none 0%| | 0/1000 [00:00] .. GENERATED FROM PYTHON SOURCE LINES 119-123 Here we plot the predicted time courses of the original model and the model with the optimized HRFs for 9 voxels where the fit improved the most. You can clearly see that, in general, the HRFs have shorter delays than the default setting. .. GENERATED FROM PYTHON SOURCE LINES 123-134 .. code-block:: Python improvement = r2_gauss_hrf - r2_gauss_gd largest_improvements = improvement.sort_values(ascending=False).index[:9] pred_gauss_gd = model_gauss.predict(parameters=pars_gauss_gd) pred_gauss_hrf = model_hrf.predict(parameters=pars_gauss_hrf) pred = pd.concat((d.loc[:, largest_improvements], pred_gauss_gd.loc[:, largest_improvements], pred_gauss_hrf.loc[:, largest_improvements]), axis=1, keys=['data', 'gauss', 'gauss+hrf'], names=['model']) # tmp = pred.stack(['model', 'source']).to_frame('value') sns.relplot(x='frame', y='value', hue='model', col='source', data=tmp.reset_index(), kind='line', col_wrap=3) .. image-sg:: /auto_examples/00_encodingdecoding/images/sphx_glr_fit_prf_005.png :alt: source = 105, source = 275, source = 794, source = 795, source = 1931, source = 2039, source = 2165, source = 2166, source = 2195 :srcset: /auto_examples/00_encodingdecoding/images/sphx_glr_fit_prf_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 137-143 Fit a Difference of Gaussians model ----------------------------------- Now we will try to fit a Difference of Gaussians model. This model has two Gaussian receptive fields, one excitatory and one inhibitory. The inhibitory receptive field is subtracted from the excitatory one. The resulting receptive field is then convolved with the HRF. .. GENERATED FROM PYTHON SOURCE LINES 143-170 .. code-block:: Python from braincoder.models import DifferenceOfGaussiansPRF2DWithHRF model_dog = DifferenceOfGaussiansPRF2DWithHRF(data=d, paradigm=paradigm, hrf_model=hrf_model, grid_coordinates=grid_coordinates, flexible_hrf_parameters=True) pars_dog_init = pars_gauss_hrf.copy() # This is the relative amplitude of the inhibitory receptive field # compared to the excitatory one. pars_dog_init['srf_amplitude'] = 0.1 # This is the relative size of the inhibitory receptive field # compared to the excitatory one. pars_dog_init['srf_size'] = 2. # Let's set up a new parameterfitter par_fitter_dog = ParameterFitter(model=model_dog, data=d, paradigm=paradigm) # Note how, for now, we are not optimizing the HRF parameters. pars_dog = par_fitter_dog.fit(init_pars=pars_dog_init, max_n_iterations=1000, fixed_pars=['hrf_delay', 'hrf_dispersion']) # Now we optimize _with_ the HRF parameters pars_dog_hrf = par_fitter_dog.fit(init_pars=pars_dog, max_n_iterations=1000) r2_dog_hrf = par_fitter_dog.get_rsq(pars_dog_hrf) sns.relplot(x='r2_hrf', y='r2_dog_hrf', data=pd.concat((r2_gauss_hrf, r2_dog_hrf), axis=1, keys=['r2_hrf', 'r2_dog_hrf']).reset_index(), kind='scatter') .. rst-class:: sphx-glr-script-out .. code-block:: pytb Traceback (most recent call last): File "/Users/gdehol/git/braincoder/examples/00_encodingdecoding/fit_prf.py", line 160, in pars_dog = par_fitter_dog.fit(init_pars=pars_dog_init, max_n_iterations=1000, ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/gdehol/git/braincoder/braincoder/optimize/parameter_fitter.py", line 251, in fit loss, gradients = compute_gradients(loss_fn, trainable_variables) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/Users/gdehol/git/braincoder/braincoder/utils/backend.py", line 443, in compute_gradients loss = loss_fn() ^^^^^^^^^ File "/Users/gdehol/git/braincoder/braincoder/optimize/parameter_fitter.py", line 212, in loss_fn pars = build_parameters() ^^^^^^^^^^^^^^^^^^ File "/Users/gdehol/git/braincoder/braincoder/optimize/parameter_fitter.py", line 410, in build return (init_pars_t * (1.0 - voxel_mask) ~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~ File "/Users/gdehol/mambaforge/envs/braincoder/lib/python3.12/site-packages/tensorflow/python/util/traceback_utils.py", line 153, in error_handler raise e.with_traceback(filtered_tb) from None File "/Users/gdehol/mambaforge/envs/braincoder/lib/python3.12/site-packages/tensorflow/python/framework/ops.py", line 6027, in raise_from_not_ok_status raise core._status_to_exception(e) from None # pylint: disable=protected-access ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ tensorflow.python.framework.errors_impl.InvalidArgumentError: {{function_node __wrapped__Mul_device_/job:localhost/replica:0/task:0/device:CPU:0}} Incompatible shapes: [2365,7] vs. [2365,9] [Op:Mul] name: .. GENERATED FROM PYTHON SOURCE LINES 174-177 Here, we plot the predicted time courses of the difference-of-gaussians model versus the original Gaussian model for the 9 voxels where the fit imoproved the most. .. GENERATED FROM PYTHON SOURCE LINES 177-188 .. code-block:: Python improvement = r2_dog_hrf - r2_gauss_hrf largest_improvements = improvement.sort_values(ascending=False).index[:9] pred_dog_hrf = model_dog.predict(parameters=pars_dog_hrf) pred = pd.concat((d.loc[:, largest_improvements], pred_gauss_hrf.loc[:, largest_improvements], pred_dog_hrf.loc[:, largest_improvements]), axis=1, keys=['data', 'gauss+hrf', 'dog+hrf'], names=['model']) tmp = pred.stack(['model', 'source']).to_frame('value') sns.relplot(x='frame', y='value', hue='model', col='source', data=tmp.reset_index(), kind='line', col_wrap=3, palette=['k'] + sns.color_palette(), hue_order=['data', 'gauss+hrf', 'dog+hrf']) .. GENERATED FROM PYTHON SOURCE LINES 189-196 Divisve Normalization PRF model ------------------------------- The most complex model we have is the DN-PRF model (Aqil et al., 2021). This model has a Gaussian excitatory receptive field, and a Gaussian inhibitory receptive field. The excitatory receptive field is divided by the sum of the excitatory and inhibitory receptive fields. The resulting receptive field is then convolved with the HRF. .. GENERATED FROM PYTHON SOURCE LINES 196-217 .. code-block:: Python from braincoder.models import DivisiveNormalizationGaussianPRF2DWithHRF model_dn = DivisiveNormalizationGaussianPRF2DWithHRF(data=d, paradigm=paradigm, hrf_model=hrf_model, grid_coordinates=grid_coordinates, flexible_hrf_parameters=True) pars_dn_init = pars_dog_hrf.copy() pars_dn_init['srf_amplitude'] = 0.01 pars_dn_init['srf_size'] = 2. pars_dn_init['baseline'] = 0.0 pars_dn_init['neural_baseline'] = 1.0 pars_dn_init['surround_baseline'] = 1.0 par_fitter_dn = ParameterFitter(model=model_dn, data=d, paradigm=paradigm) # Without HRF pars_dn = par_fitter_dn.fit(init_pars=pars_dn_init, max_n_iterations=1000, fixed_pars=['hrf_delay', 'hrf_dispersion']) # With HRF pars_dn = par_fitter_dn.fit(init_pars=pars_dn, max_n_iterations=1000) .. GENERATED FROM PYTHON SOURCE LINES 218-219 Again, let's plot the R2 improvements .. GENERATED FROM PYTHON SOURCE LINES 219-225 .. code-block:: Python r2_dn = par_fitter_dn.get_rsq(pars_dn) sns.relplot(x='r2_dog_hrf', y='r2_dn', data=pd.concat((r2_dog_hrf, r2_dn), axis=1, keys=['r2_dog_hrf', 'r2_dn']).reset_index(), kind='scatter') plt.plot([0, 1], [0, 1], 'k--') .. GENERATED FROM PYTHON SOURCE LINES 226-236 .. code-block:: Python improvement = r2_dn - r2_dog_hrf largest_improvements = improvement.sort_values(ascending=False).index[:9] pred_dn = model_dn.predict(parameters=pars_dn) pred = pd.concat((d.loc[:, largest_improvements], pred_dog_hrf.loc[:, largest_improvements], pred_dn.loc[:, largest_improvements]), axis=1, keys=['data', 'dog+hrf', 'dn+hrf'], names=['model']) tmp = pred.stack(['model', 'source']).to_frame('value') sns.relplot(x='frame', y='value', hue='model', col='source', data=tmp.reset_index(), kind='line', col_wrap=3, palette=['k'] + sns.color_palette(), hue_order=['data', 'dog+hrf', 'dn+hrf']) .. GENERATED FROM PYTHON SOURCE LINES 237-267 .. code-block:: Python # Decoding # -------- # We can also use the fitted models to decode the stimulus from the # fMRI response. Let's compare our simplest model versus our most # complex model. # First we fit the noise models from braincoder.optimize import ResidualFitter, StimulusFitter # Let's first get grid coordinates and paradigm at a slightly lower resolution data = load_szinte2024(resize_factor=2.5) grid_coordinates = data['grid_coordinates'] paradigm = data['stimulus'] best_voxels_gauss = r2_gauss_gd[pars_gauss_gd['sd'] > 0.5].sort_values(ascending=False).index[:200] model_gauss = GaussianPRF2DWithHRF(data=d[best_voxels_gauss], hrf_model=hrf_model, grid_coordinates=grid_coordinates.astype(np.float32), parameters=pars_gauss_gd.loc[best_voxels_gauss].astype(np.float32)) resid_fitter_gauss = ResidualFitter(model=model_gauss, data=d.loc[:, best_voxels_gauss], paradigm=paradigm.astype(np.float32), parameters=pars_gauss_gd.loc[best_voxels_gauss].astype(np.float32)) omega_gauss, _ = resid_fitter_gauss.fit() .. GENERATED FROM PYTHON SOURCE LINES 271-283 .. code-block:: Python best_voxels_dn = r2_dn[pars_dn['sd'] > 0.5].sort_values(ascending=False).index[:200] model_dn = DivisiveNormalizationGaussianPRF2DWithHRF(data=d[best_voxels_dn], hrf_model=hrf_model, grid_coordinates=grid_coordinates.astype(np.float32), parameters=pars_dn.loc[best_voxels_dn].astype(np.float32)) resid_fitter_dn = ResidualFitter(model=model_dn, data=d.loc[:, best_voxels_dn], paradigm=paradigm, parameters=pars_dn.loc[best_voxels_dn]) omega_dn, _ = resid_fitter_dn.fit() .. GENERATED FROM PYTHON SOURCE LINES 284-287 Decoded stimulus: Gaussian model =============================== Now we can decode the stimulus from the fMRI responses .. GENERATED FROM PYTHON SOURCE LINES 287-290 .. code-block:: Python stim_fitter_gauss = StimulusFitter(model=model_gauss, data=d.loc[:, best_voxels_gauss], omega=omega_gauss) stim_gauss = stim_fitter_gauss.fit(l2_norm=0.01, learning_rate=0.01, max_n_iterations=1000) .. GENERATED FROM PYTHON SOURCE LINES 291-315 .. code-block:: Python from matplotlib.animation import FuncAnimation from IPython.display import HTML def play_reconstruction(reconstructed_stimulus): # Here we make a movie of the decoded stimulus # Set up a function to draw a single frame vmin, vmax = 0.0, np.quantile(reconstructed_stimulus.values.ravel(), 0.99) def update(frame): plt.clf() # Clear the current figure plt.imshow(reconstructed_stimulus.stack('y').loc[frame].iloc[::-1, :], cmap='viridis', vmin=vmin, vmax=vmax) plt.axis('off') plt.title(f"Frame {frame}") # Create the animation fig = plt.figure() ani = FuncAnimation(fig, update, frames=range(paradigm.shape[0]), interval=100) return HTML(ani.to_html5_video()) play_reconstruction(stim_gauss) .. GENERATED FROM PYTHON SOURCE LINES 316-318 Decoded stimulus: DN model ========================== .. GENERATED FROM PYTHON SOURCE LINES 318-321 .. code-block:: Python stim_fitter_dn = StimulusFitter(model=model_dn, data=d.loc[:, best_voxels_dn], omega=omega_dn) stim_dn = stim_fitter_dn.fit(l2_norm=0.01, learning_rate=0.01, max_n_iterations=1000) .. GENERATED FROM PYTHON SOURCE LINES 322-325 .. code-block:: Python play_reconstruction(stim_dn) # As you can see, the DN model works a lot better than the Gaussian model. ;) .. rst-class:: sphx-glr-timing **Total running time of the script:** (13 minutes 37.291 seconds) .. _sphx_glr_download_auto_examples_00_encodingdecoding_fit_prf.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: fit_prf.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: fit_prf.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: fit_prf.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_