The cosine bump basis functions used in FilterNet seem to be different from the formula in Pillow et.al. J Neurosci 2005

Open
  • Is something wrong?
#334 0 comments 0 reactions 0 assignees View on GitHub

Maintainers usually reply within 1 day

Nobody has claimed this yet.

Assessment

Difficulty
4/5
Estimated time
3-5 days
Newbie friendliness
38/100
Issue type
Documentation
Clarity
Mostly clear
Activity status
Stale
Tech stack
matplotlib, numpy, python
Domain
backend, documentation, performance

Research direction

Start with makeBasis_StimKernel and ff() in bmtk/simulator/filternet/lgnmodel/fitfuns.py, then compare them with the cited Pillow et al. formula and the examples in bmtk/docs/tutorial/07_filter_models.ipynb and filternet_images_helpers.py. Done means documenting b, ylim, db, kpeaks, normalization, and the formula difference, while assessing whether ff() should be vectorized.

Written by the indexing model from the issue text.

Description

The consine bump basis function of LGN cells in Pillow et.al. J Neurosci 2005 is:

B(t) = (cos(log[t+tau] - phi) + 1) / 2,  if phi-pi < log[t+tau] < phi+pi

However, the actual function makeBasis_StimKernel implemented in bmtk/simulator/filternet/lgnmodel/fitfuns.py (from line 28-81 and line 97-122) seems to be more complex than that and lacks documentation (the tutorial in bmtk/docs/tutorial/07_filter_models.ipynb is not enough). For example, the values of variable b (line 33), variable ylim (line 37), and variable db (line 40) are all hard-coded, do not appear in the original equation, and I can't find an explanation on what they mean. Besides, I am also not sure which variable in the original equation does kpeaks represent, and why there is an extra normalization step.

There is also a demo in bmtk/docs/tutorial/helpers/filternet_images_helpers.py that draws a figure of basis function with different parameters. I am able to reproduce that figure with a much simpler code that summarize the actual computation bmtk performs. The computation in get_temporal_kernel() in the following code is indeed different from the original equation in Pillow et.al. J Neurosci 2005. So, what is the reason that bmtk uses a different formula?

import numpy as np
import matplotlib.pyplot as plt

#%% Helper functions
def get_temporal_kernel(t, kpeaks, delay, weights):
    '''
    Inputs:
        t: [length,1] array
        kpeaks: [1,2] array
        delay: [1,2] array
        weights: [2,] array
    Output:
        kernel: [length,] array
    '''
    log_t = np.log(t+1.3-delay) - np.log(kpeaks)  # Use 1.3 not 0.3 here to compensate for the extra 1-point offset in line 65 of `fitfuns.py`.
    log_t_pi = np.clip(np.pi*log_t/get_temporal_kernel.db2, -np.pi, np.pi)
    log_t_pi[np.isnan(log_t_pi)] = -np.pi
    basis = (np.cos(log_t_pi) + 1) / 2.0
    basis_norm = basis / np.linalg.norm(basis, ord=2, axis=0)
    kernel = basis_norm @ weights
    return kernel
get_temporal_kernel.db2 = 2.0 * np.diff(np.log([100.3, 200.3]))  # 2*db in `fitfuns.py`

#%% Plot kernels with different parameters
weights = np.array([[30.0, -20.0], [30.0, -1.0], [15.0, -20.0]])
kpeaks = np.array([[3.0, 5.0], [3.0, 30.0], [20.0, 40.0]])
delays = np.array([[.0, 0.0], [0.0, 60.0], [20.0, 60.0]])

t = np.expand_dims(np.arange(0,150,1), axis=1)  # milliseconds

fig, axes = plt.subplots(3, 3, figsize=(10, 7))
ri = ci = 0

for ci in range(weights.shape[0]):
    kernel = get_temporal_kernel(t, np.array([[9.67, 20.03]]), np.array([[0.0, 1.0]]), weights[ci,:])
    idx = np.abs(kernel) > 0.0
    axes[ri, ci].plot(-t[idx]/1000, kernel[idx])
    axes[ri, ci].set_ylim([-3.5, 10.0])
    axes[ri, ci].text(0.05, 0.90, 'weights={}'.format(weights[ci,:]), horizontalalignment='left', verticalalignment='top', transform=axes[ri, ci].transAxes)
axes[0, 0].set_ylabel('effect of weights')
ri += 1

# kpeaks parameters controll the spread of both peaks, the second peak must have a bigger spread
for ci in range(kpeaks.shape[0]):
    kernel = get_temporal_kernel(t, kpeaks[[ci],:], np.array([[0.0, 1.0]]), np.array([30.0, -20.0]))
    idx = np.abs(kernel) > 0.0
    axes[ri, ci].plot(-t[idx]/1000, kernel[idx])
    axes[ri, ci].set_xlim([-0.15, 0.005])
    axes[ri, ci].text(0.05, 0.90, 'kpeaks={}'.format(kpeaks[ci,:]), horizontalalignment='left', verticalalignment='top', transform=axes[ri, ci].transAxes)
axes[1, 0].set_ylabel('effects of kpeaks')
ri += 1

for ci in range(delays.shape[0]):
    kernel = get_temporal_kernel(t, np.array([[9.67, 20.03]]), delays[[ci],:], np.array([30.0, -20.0]))
    idx = np.abs(kernel) > 0.0
    axes[ri, ci].plot(-t[idx]/1000, kernel[idx])
    axes[ri, ci].set_xlim([-0.125, 0.001])
    axes[ri, ci].text(0.05, 0.90, 'delays={}'.format(delays[ci,:]), horizontalalignment='left', verticalalignment='top', transform=axes[ri, ci].transAxes)
axes[2, 0].set_ylabel('effects of delays')

fig.show()

Besides, the function ff() in line 108-117 of bmtk/simulator/filternet/lgnmodel/fitfuns.py loops over an numpy array one-by-one, but this can be easily vectorized. Why does it choose a slower computation?

Dominant language
Python
Stars
339
Forks
103
Avg merge
2h 7m
Merged PRs (30d)
14

Getting set up

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

More from AllenInstitute/bmtk

All issues in AllenInstitute/bmtk

Similar issues

More Python issues

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.