Frequency problem in template & recording

[ZoeChen96]

Hi there, I am using MEArec to generate neuropixel1.0 (128 channel) simulated dataset, but I can't understand the template length in the recording. These are my steps: 1) I generate a template with 32kHz, and I find the template size is (3900, 128, 224). It makes sense because the default 'cut-out' is [2,5], and 32kHz(2+5)ms = 224. 2) I tried to generate a recording with 32kHz templates. The template size in the recording is (60, 10, 128, 416), which also makes sense because the 'pad_len' is [3,3], and 32kHz(2+5+3+3) = 416.

Then weird things happens: 1) I generate a recording with 32kHz, but I specify the fs to be 30kHz. This time the template size in the recording is (60, 10, 128, 378), but the 378 is weird. I assume the length should be 416*30kHz/32kHz = 390. 2) Based on 1), I also generate a recording with fs to be 20kHz and 40kHz. The length of them, to sum up, are: 20kHz -- 215 30kHz --378 40kHz -- 580

3) I generate a template with 40kHz ('dt = 0.0025'), and the template size is still (3900, 128, 224). I used the template to generate the recording in 20kHz/30kHz/40kHz, and the length of them are: 20kHz -- 172 30kHz --303 40kHz --464.

Besides, I ran sorters in spikeinterface of the above cases( 20k/30k/40k recordings with 40kHz template+ 20k/30k/40k recordings with 32kHz template), and I found the former ones (20k/30k/40k recordings with 40kHz template) have very poor accuracy; while the latter ones (20k/30k/40k recordings with 32kHz template) seems reasonable. Could you have a look at it? I will paste my scripts for generating the template and recording below.

[ZoeChen96]

The template generating script is:

import MEArec as mr
#settings
probe = 'Neuropixels2-128'              #available: CLI: mearec available-probes
freq  = 40                              #40kHz
#folders
dir_full_templates = '/xxx/mearec/temp_datasets/templates300_'+str(freq)+'kHz_minampnull_'+ probe +'.h5'
templates_params = mr.get_default_templates_params()
cell_models_folder = mr.get_default_cell_models_folder()
#parameters
templates_params['probe'] = probe
templates_params['n'] = 300                 #for every cell model (in total 13 cell model), generate 300 templates
templates_params['seed']=0
templates_params['min_amp']=0
templates_params['dt']=1/freq
tempgen = mr.gen_templates(cell_models_folder=cell_models_folder, params=templates_params, verbose=True)
mr.save_template_generator(tempgen=tempgen, filename=dir_full_templates)

[ZoeChen96]

The recording generation script is:

import MEArec as mr
import scipy.io
import spikeinterface.extractors as se
import warnings
warnings.simplefilter("ignore")
############ set global parameters ####################
probe       = 'Neuropixels2-128'
n_exc       = 48
n_inh       = 12
duration    = 6 #seconds
filtering   = False #filtering True/false
noise_type  = 'far-neurons' #uncorrelated/distance-correlated/far-neurons
recording_params = mr.get_default_recordings_params()
recording_params['spiketrains']['n_exc'] = n_exc
recording_params['spiketrains']['n_inh'] = n_inh
recording_params['spiketrains']['duration'] = duration
recording_params['spiketrains']['seed'] = 0
recording_params['templates']['min_amp'] = 40
recording_params['templates']['max_amp'] = 300
recording_params['templates']['seed'] = 0
recording_params['recordings']['modulation'] = 'electrode'
recording_params['recordings']['noise_mode'] = noise_type
recording_params['recordings']['noise_level'] = 10
recording_params['recordings']['chunk_conv_duration'] = 20
recording_params['recordings']['chunk_noise_duration'] = 20
recording_params['recordings']['chunk_filter_duration'] = 20
recording_params['recordings']['seed'] = 0
#filter settings
recording_params['recordings']['filter'] = filtering
recording_params['recordings']['filter_cutoff'] = [300, 6000] # highpass from 0
#seeds settings
recording_params['seeds']['spiketrains'] = 0
recording_params['seeds']['templates'] = 0
recording_params['seeds']['convolution'] = 0
recording_params['seeds']['noise'] = 0
########## load templates ###################
dir_full_templates = '/xxx/mearec/temp_datasets/templates300_32kHz_minampnull_'+ probe +'.h5'
tempgen = mr.load_templates(dir_full_templates)
############## set looping parameters and generate datasets ############
fs          = [20000 30000 40000]
for i in fs:
    print('generating recording for fs', str(i))
    recording_params['recordings']['fs'] = i
    dir_full_recordings = '/xxx/mearec/temp_datasets/recording_'+str(n_exc+n_inh) +'cells_' + '_noise10' +str(filtering)+'filter_'+str(i)+'fs_'+probe +'.h5'
    dir_full_gttimes = '/xxx/mearec/temp_datasets/gttimes_'+str(n_exc+n_inh) +'cells_' + '_noise10' +str(filtering)+'filter_'+str(i)+'fs_'+probe +'.mat'
    recgen = mr.gen_recordings(templates=dir_full_templates, params=recording_params, verbose=True)
    # save recording
    mr.save_recording_generator(recgen, dir_full_recordings)
    #save gt_times
    if (n_exc+n_inh >0):
        sorting_GT = se.MEArecSortingExtractor(dir_full_recordings)
        scipy.io.savemat(dir_full_gttimes, {'gt_times': sorting_GT.to_spike_vector()})

[ZoeChen96]

Besides, as I did not specify the 'cut-out' and 'pad-len', I check that in the .info. I can paste the info for recording (30k), for example:

{'cell_types': {'excitatory': ['PC', 'SS', 'SP'],
  'inhibitory': ['AC', 'BP', 'BC', 'BTC', 'ChC', 'DBC', 'MC', 'NGC']},
 'electrodes': {'description': 'Neuropixels 2,0 probe. 384 square contacts in 2 columns.',
  'dim': [64, 2],
  'electrode_name': 'Neuropixels2-128',
  'pitch': [15.0, 32.0],
  'plane': 'yz',
  'shape': 'square',
  'size': 6.0,
  'sortlist': None,
  'type': 'mea'},
 'recordings': {'adc_bit_depth': None,
  'angle_tol': 15,
  'bursting': False,
  'bursting_units': None,
  'chunk_conv_duration': 20,
  'chunk_duration': 10,
  'chunk_filter_duration': 20,
  'chunk_noise_duration': 20,
  'color_noise_floor': 1,
  'color_peak': 300,
  'color_q': 2,
  'drift_fs': 100,
  'drift_mode_probe': 'rigid',
  'drift_mode_speed': 'slow',
  'drifting': False,
  'dtype': 'float32',
  'duration': 60.0,
  'exp_decay': 0.2,
  'extract_waveforms': False,
  'far_neurons_exc_inh_ratio': 0.8,
  'far_neurons_max_amp': 10,
  'far_neurons_n': 300,
  'far_neurons_noise_floor': 0.5,
  'fast_drift_max_jump': 20,
  'fast_drift_min_jump': 5,
  'fast_drift_period': 10,
  'filter': False,
  'filter_cutoff': [300, 6000],
  'filter_order': 3,
  'fs': 30000.0,
  'gain': None,
  'lsb': None,
  'max_burst_duration': 100,
  'modulation': 'electrode',
  'mrand': 1,
  'n_burst_spikes': 10,
  'n_bursting': None,
  'n_drifting': None,
  'n_neurons': 60,
  'noise_color': False,
  'noise_half_distance': 30,
  'noise_level': 10,
  'noise_mode': 'far-neurons',
  'non_rigid_gradient_mode': 'linear',
  'overlap': False,
  'preferred_dir': [0, 0, 1],
  'sdrand': 0.05,
  'seed': 0,
  'shape_mod': False,
  'shape_stretch': 30.0,
  'slow_drift_amplitude': None,
  'slow_drift_velocity': 5,
  'slow_drift_waveform': 'triangluar',
  'sync_jitt': 1,
  'sync_rate': None,
  't_end_drift': None,
  't_start_drift': 0},
 'seeds': {'convolution': 0, 'noise': 0, 'spiketrains': 0, 'templates': 0},
 'spiketrains': {'duration': 60,
  'f_exc': 5,
  'f_inh': 15,
  'gamma_shape': 2,
  'min_rate': 0.5,
  'n_exc': 240,
  'n_inh': 60,
  'process': 'poisson',
  'ref_per': 2,
  'seed': 0,
  'st_exc': 1,
  'st_inh': 3,
  't_start': 0},
 'templates': {'cut_out': [2, 5],
  'max_amp': 300,
  'min_amp': 40,
  'min_dist': 25,
  'n_jitters': 10,
  'n_overlap_pairs': None,
  'overlap_threshold': 0.9,
  'overlapping': [array([ 0, 34]),
   array([ 0, 40]),
   array([ 7, 15]),
   array([8, 9]),
   array([11, 27]),
   array([12, 23]),
   array([13, 53]),
   array([19, 55]),
   array([27, 50]),
   array([29, 36]),
   array([24, 31]),
   array([31, 45]),
   array([33, 46]),
   array([34, 40]),
   array([50, 52]),
   array([27, 52]),
   array([ 9, 54]),
   array([47, 57])],
  'pad_len': [3, 3],
  'seed': 0,
  'smooth_percent': 0.5,
  'smooth_strength': 1,
  'upsample': 8,
  'xlim': None,
  'ylim': None,
  'zlim': None}}

[ZoeChen96]

And here is the .info for templates (40k):

{'electrodes': {'description': 'Neuropixels 2,0 probe. 384 square contacts in '
                               '2 columns.',
                'dim': [64, 2],
                'electrode_name': 'Neuropixels2-128',
                'pitch': [15.0, 32.0],
                'plane': 'yz',
                'shape': 'square',
                'size': 6.0,
                'sortlist': None,
                'type': 'mea'},
 'params': {'beta_distr_params': [1.5, 5.0],
            'cell_models_folder': '/imec/users/chen14/.config/mearec/1.8.0/cell_models/bbp',
            'check_for_drift_amp': False,
            'cut_out': [2, 5],
            'delay': 10,
            'det_thresh': 30,
            'drift_steps': 31,
            'drift_within_bounds': False,
            'drift_xlim': [-10, 10],
            'drift_ylim': [-10, 10],
            'drift_zlim': [30, 80],
            'drifting': False,
            'dt': 0.025,
            'max_drift': 100,
            'max_iterations': 1000,
            'min_amp': 0,
            'min_drift': 30,
            'n': 300,
            'ncontacts': 10,
            'offset': 0,
            'overhang': 30,
            'probe': 'Neuropixels2-128',
            'rot': 'physrot',
            'seed': 0,
            'sim_time': 1,
            'target_spikes': [3, 50],
            'templates_folder': '/imec/users/chen14/.config/mearec/1.8.0/templates',
            'timeout': None,
            'weights': [0.25, 1.75],
            'x_distr': 'uniform',
            'xlim': [10, 80],
            'ylim': None,
            'zlim': None}}

[alejoe91]

HI @ZoeChen96

Thanks for looking into this. I'm a bit busy thes days but I'll try to take a look next week.

[ZoeChen96]

Ok, take your time. Thanks, Alessio!

[alejoe91]

@ZoeChen96 I have to double check, but I think that:

1. the problem with generating new templates with different dt is that if an intracellular simulation is already present, the simulator will skip the NEURON sim and directly load the currents and compute extracellular spikes.

2. This behavior is wrong, because the simulator thinks it has 40kHz templates, and maybe downstream something breaks (affecting spike sorting).

I'll run a small test, and push a fix in case this is the issue ;)

[ZoeChen96]

Hi Alessio, thanks for your response! I will test it out. BTW, I encountered a problem before when I generate templates (maybe because of np version:) Traceback (most recent call last):

 File "/imec/other/macaw/projectdata/mearec_datasets/new_templates_datasets/template_gen.py", line 16, in <module>

   mr.save_template_generator(tempgen=tempgen, filename=dir_full_templates)

 File "/imec/other/macaw/chen14/gitrepository/MEArec/MEArec/tools.py", line 536, in save_template_generator

   save_dict_to_hdf5(tempgen.info, f, 'info/')

 File "/imec/other/macaw/chen14/gitrepository/MEArec/MEArec/tools.py", line 639, in save_dict_to_hdf5

  recursively_save_dict_contents_to_group(h5file, path, dic)

 File "/imec/other/macaw/chen14/gitrepository/MEArec/MEArec/tools.py", line 656, in recursively_save_dict_contents_to_group

   if isinstance(item, (int, float, np.integer, np.float, str, bytes, np.bool_)):

                                               ^^^^^^^^

 File "/imec/other/macaw/chen14/anaconda43/envs/mearec_git_411/lib/python3.11/site-packages/numpy/__init__.py", line 305, in __getattr__

   raise AttributeError(__former_attrs__[attr])

AttributeError: module 'numpy' has no attribute 'float'.

`np.float` was a deprecated alias for the builtin `float`. To avoid this error in existing code, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.

The aliases was originally deprecated in NumPy 1.20; for more details and guidance see the original release note at:

   https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations. Did you mean: 'cfloat'?

So I delete np.float in Line 656 of file tools.py. https://github.com/alejoe91/MEArec/blob/708be8eed4656fad660aac6d3e0760646a3abf29/MEArec/tools.py#L656

FYI.

[alejoe91]

Uh let me fix that!

[ZoeChen96]

Sorry I put the comment in the wrong place... It's actually for #140

Hi Alessio, thanks for your response! I will test it out. BTW, I encountered a problem before when I generate templates (maybe because of np version:) Traceback (most recent call last):

 File "/imec/other/macaw/projectdata/mearec_datasets/new_templates_datasets/template_gen.py", line 16, in <module>

   mr.save_template_generator(tempgen=tempgen, filename=dir_full_templates)

 File "/imec/other/macaw/chen14/gitrepository/MEArec/MEArec/tools.py", line 536, in save_template_generator

   save_dict_to_hdf5(tempgen.info, f, 'info/')

 File "/imec/other/macaw/chen14/gitrepository/MEArec/MEArec/tools.py", line 639, in save_dict_to_hdf5

  recursively_save_dict_contents_to_group(h5file, path, dic)

 File "/imec/other/macaw/chen14/gitrepository/MEArec/MEArec/tools.py", line 656, in recursively_save_dict_contents_to_group

   if isinstance(item, (int, float, np.integer, np.float, str, bytes, np.bool_)):

                                               ^^^^^^^^

 File "/imec/other/macaw/chen14/anaconda43/envs/mearec_git_411/lib/python3.11/site-packages/numpy/__init__.py", line 305, in __getattr__

   raise AttributeError(__former_attrs__[attr])

AttributeError: module 'numpy' has no attribute 'float'.

`np.float` was a deprecated alias for the builtin `float`. To avoid this error in existing code, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.

The aliases was originally deprecated in NumPy 1.20; for more details and guidance see the original release note at:

   https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations. Did you mean: 'cfloat'?

So I delete np.float in Line 656 of file tools.py.

https://github.com/alejoe91/MEArec/blob/708be8eed4656fad660aac6d3e0760646a3abf29/MEArec/tools.py#L656

FYI.

[alejoe91]

Fixed in the same PR #141