updated LFP plots (#229)

* Accept iterable for bin sizes for pairwise spike-train correlations
Fixes #216

* corrcoef plot fixes

* --mpi=pmi2 not working anymore

* introduce t_sim_lfp<=t_sim parameter to independently control LFP calculations

* correlation analysis

* Backup files before updating master

* prep merge

* roll back some edits

* some minor plot updates

* add lfp script

* improved figure legibility

* Update base_analysis_params.py

* remove spaces

* removed unused files

* roll back change

* add back deleted ms_figures.correlation() function

* roll back edit

* line change

mesocircuit/analysis/spike_analysis.py

@@ -951,7 +951,7 @@ def _compute_ccs_distances(X, circuit, sptrains_X, binsize_time, ccs_time_interv
     binsize_time
         Temporal resolution of sptrains_X (in ms).
     ccs_time_interval
-        Temporal bin size(s) for correlation coefficients calculation (in ms).
+        Temporal bin size for correlation coefficients calculation (in ms)
     positions_X
         Positions of population X.
     """
@@ -1018,7 +1018,7 @@ def _compute_ccs_distances(X, circuit, sptrains_X, binsize_time, ccs_time_interv
                            extent=[circuit.net_dict['extent']] * 2,
                            edge_wrap=True)
 
-    ccs_dic = {f'ccs_{ccs_time_interval}ms': ccs,
+    ccs_dic = {f'ccs_{ccs_time_interval}': ccs,
                'distances_mm': distances}
     return ccs_dic
 

mesocircuit/lfp/lfp_parameters.py

@@ -212,7 +212,7 @@ def get_parameters(path_lfp_data=None, sim_dict=None, net_dict=None):
         savefolder=path_lfp_data,
         # derived parameters for CachedTopoNetwork instance
         network_params=dict(
-            simtime=sim_dict['t_presim'] + sim_dict['t_sim'],
+            simtime=sim_dict['t_presim'] + sim_dict['t_sim_lfp'],
             dt=2**-2,
             # sim_dict['sim_resolution'],  # run at 4kHz; we're downsampling
             # outputs anyway.
@@ -288,7 +288,7 @@ def get_parameters(path_lfp_data=None, sim_dict=None, net_dict=None):
             'lambda_f': 100,
             'dt': 2**-2,  # sim_dict['sim_resolution'],
             'tstart': 0,
-            'tstop': sim_dict['t_presim'] + sim_dict['t_sim'],
+            'tstop': sim_dict['t_presim'] + sim_dict['t_sim_lfp'],
             'verbose': False,
         }
     ))

mesocircuit/lfp/plotting.py

@@ -462,7 +462,8 @@ def layout_illustration(
 def plot_single_channel_lfp_data(
     ax, PS, fname, title='LFP', ylabel=r'$\Phi$ (mV)', T=[
         500, 550], CONTACTPOS=(
-            (-200, 200), (200, -200)), subtract_mean=True):
+            (-200, 200), (200, -200)), subtract_mean=True, 
+            report_corrcoefs=False):
     '''
     Arguments
     ---------
@@ -476,6 +477,8 @@ def plot_single_channel_lfp_data(
         x and y coordinate of electrode contact points in PS.electrodeParams
     subtract_mean: bool
         if True, subtract mean value trace
+    report_corrcoefs: bool
+        if True, print pairwise correlations in axes title
     '''
     with h5py.File(fname, 'r') as f:
         srate = f['srate'][()]
@@ -491,6 +494,20 @@ def plot_single_channel_lfp_data(
                 ax.plot(tvec, data[tinds] - data[tinds].mean(), f'C{i}')
             else:
                 ax.plot(tvec, data[tinds], f'C{i}')
+
+    if report_corrcoefs:
+        with h5py.File(fname, 'r') as f:
+            tind = int(T[0] * srate / 1000)
+            CONTACTS = []
+            for i, CPOS in enumerate(CONTACTPOS):
+                CONTACT = (PS.electrodeParams['x'] == CPOS[0]
+                            ) & (PS.electrodeParams['y'] == CPOS[1])
+                CONTACTS += np.where(CONTACT)
+            CONTACTS = np.array(CONTACTS).ravel()
+            data = f['data'][()][CONTACTS, tind:]
+        cc = np.corrcoef(data)
+        print(f'{title} correlations', CONTACTS + 1, CONTACTPOS, cc)
+    #else:
     ax.set_title(title)
     ax.set_ylabel(ylabel)
     remove_axis_junk(ax)
@@ -501,7 +518,8 @@ def plot_single_channel_csd_data(
         title='CSD',
         ylabel=r'CSD ($\frac{\mathrm{nA}}{\mathrm{µm}^3}$)',
         T=[500, 550],
-        CONTACTPOS=((-200, 200), (200, -200))):
+        CONTACTPOS=((-200, 200), (200, -200)),
+        report_corrcoefs=False):
     # CSD bin edges
     X, Y, Z = np.meshgrid(PS.CSDParams['x'],
                           PS.CSDParams['y'],
@@ -523,6 +541,16 @@ def plot_single_channel_csd_data(
             CONTACT = (Xmid == CPOS[0]) & (Ymid == CPOS[1])
             data = f['data'][()][CONTACT, ].flatten()
             ax.plot(tvec, data[tinds] - data[tinds].mean(), f'C{i}')
+    if report_corrcoefs:
+        datas = []
+        with h5py.File(fname, 'r') as f:
+            tind = int(T[0] * srate / 1000)
+            for i, CPOS in enumerate(CONTACTPOS):
+                CONTACT = (Xmid == CPOS[0]) & (Ymid == CPOS[1])
+                datas += [f['data'][()][CONTACT, ].flatten()[tind:]]
+        cc = np.corrcoef(np.array(datas))
+        print(f'{title} correlations', cc)
+
     ax.set_title(title)
     ax.set_ylabel(ylabel)
     remove_axis_junk(ax)
@@ -1244,7 +1272,8 @@ def plot_coherence_vs_distance_vs_frequency(
         method='mlab',
         phase_coherence=False,
         title='LFP',
-        show_cbar=True
+        show_cbar=True,
+        clim=[0, 0.25]
 ):
 
     with h5py.File(fname, 'r') as f:
@@ -1279,8 +1308,8 @@ def plot_coherence_vs_distance_vs_frequency(
         unique * 1E-3,  # [µm] -> [mm] unit conversion
         chfreqs,
         means,
-        vmin=0,
-        vmax=0.25,
+        vmin=clim[0],
+        vmax=clim[1],
         cmap='viridis',
         shading='auto')
 

mesocircuit/parameterization/base_simulation_params.py

@@ -30,6 +30,10 @@
     # resolution of the simulation (in ms)
     'sim_resolution': 0.1,
 
+    # simulation time for LFP predictions (in ms),
+    # should be equal to or shorter than 't_sim' above.
+    't_sim_lfp': 1000.0,
+
     # list of recording devices, default is 'spike_recorder'. A 'voltmeter' can
     # be added to record membrane voltages of the neurons. Nothing will be
     # recorded if an empty list is given.

mesocircuit/plotting/figures.py

@@ -4,14 +4,12 @@
 Definition of default figures plotted with plotting.py.
 """
 
-import mesocircuit.plotting.plotting as plot
+from mesocircuit.plotting import plotting as plot
 from mesocircuit.parameterization import helpers_analysis as helpana
 import matplotlib.gridspec as gridspec
 import matplotlib.pyplot as plt
 import numpy as np
 import h5py
-import matplotlib
-matplotlib.use('Agg')
 
 
 def parameters(circuit):
@@ -313,14 +311,11 @@ def statistics_overview(circuit, all_FRs, all_LVs, all_CCs_distances, all_PSDs):
     for X in all_CCs_distances:
         if isinstance(all_CCs_distances[X], h5py._hl.group.Group):
             try:
-                iter(circuit.ana_dict['ccs_time_interval'])
-                ccs_time_interval = circuit.ana_dict['ccs_time_interval'][0]
-                all_CCs[X] = all_CCs_distances[X][f'ccs_{ccs_time_interval}ms']
-
-                print('CCs in statistics_overview use only first value in list of time intervals: ' +
-                      f'{ccs_time_interval}')
-            except TypeError:
-                all_CCs[X] = all_CCs_distances[X][f'ccs_{circuit.ana_dict["ccs_time_interval"]}ms']
+                for i, ccs_time_interval in enumerate(iter(circuit.ana_dict['ccs_time_interval'])):
+                    all_CCs[X] = all_CCs_distances[X][f'ccs_{ccs_time_interval}']
+                    break  # plot statistics overview using only first value in list for now
+            except TypeError as _:
+                all_CCs[X] = all_CCs_distances[X][f'ccs_{circuit.ana_dict["ccs_time_interval"]}']
         else:
             all_CCs[X] = np.array([])
 
@@ -365,20 +360,20 @@ def _corrcoef_distance(key_ccs):
             xlabel='distance (mm)',
             ylabel=r'$CC$')
         return fig
-
+    
     try:
-        # one figure per time interval if a list is provided
+        # stats plot for different time bins 
         for i, ccs_time_interval in enumerate(iter(circuit.ana_dict['ccs_time_interval'])):
             key_ccs = f'ccs_{ccs_time_interval}'
             fig = _corrcoef_distance(key_ccs)
             plot.savefig(circuit.data_dir_circuit, circuit.plot_dict['extension'],
-                         f'corrcoef_distance_{key_ccs}ms')
+                    f'corrcoef_distance_{key_ccs}ms')
     except TypeError as _:
         ccs_time_interval = circuit.ana_dict['ccs_time_interval']
         key_ccs = f'ccs_{ccs_time_interval}'
         fig = _corrcoef_distance(key_ccs)
         plot.savefig(circuit.data_dir_circuit, circuit.plot_dict['extension'],
-                     f'corrcoef_distance_{key_ccs}ms')
+                    f'corrcoef_distance_{key_ccs}ms')
     return
 
 

mesocircuit/plotting/plotting.py

@@ -23,7 +23,7 @@
 import matplotlib
 from scipy.optimize import curve_fit
 from mesocircuit.helpers.io import load_h5_to_sparse_X
-matplotlib.use('Agg')
+
 matplotlib.rcParams.update(rcParams)
 
 # initialize MPI
@@ -1368,7 +1368,7 @@ def plotfunc_CCs_distance(
         return
 
     distances = data[X]['distances_mm'][:max_num_pairs]
-    ccs = data[X][key_ccs + 'ms'][:max_num_pairs]
+    ccs = data[X][key_ccs][:max_num_pairs]
 
     # loop for reducing zorder-bias
     blocksize = int(len(distances) / nblocks)
@@ -1517,7 +1517,7 @@ def savefig(
     data_dir_circuit
         Path to data directory
     extension
-        File extension.
+        file extension
     filename
         Name of the file.
     eps_conv