Model and experimental timespace data in one figure

Author: m-kowalska

figures/npx/figure_traub_timespace.py

@@ -9,10 +9,12 @@
 import os
 from traub_data_kcsd_column_figure import (prepare_electrodes, prepare_pots,
                                            do_kcsd, set_axis)
+import kCSD2D_reconstruction_from_npx as npx
+from scipy.signal import filtfilt, butter
 
 
 def make_plot_spacetime(ax, val, cut=9, title='True CSD',
-                        cmap=plt.cm.bwr, letter='A'):
+                        cmap=plt.cm.bwr, letter='A', ylabel=True):
     yy = np.linspace(-3500, 500, val.shape[1])
     xx = np.linspace(-50, 250, val[cut, :, :].shape[1])
     max_val = np.max(np.abs(val[cut, :, :]))
@@ -26,23 +28,50 @@ def make_plot_spacetime(ax, val, cut=9, title='True CSD',
                         color='k', ls='--')
             plt.text(110, value+165, name[i], fontsize=10, va='top', ha='center')
     ax.set_xlabel('Time (ms)')
-    ax.set_ylabel('Y ($\mu$m)')
+    if ylabel:
+        ax.set_ylabel('Y ($\mu$m)')
     ax.set_title(title, fontsize=20, pad=30)
+    ax.set_xlim(-50, 100)
     ticks = np.linspace(-max_val, max_val, 3, endpoint=True)
-    plt.colorbar(im, orientation='vertical', format='%.3f', ticks=ticks)
+    plt.colorbar(im, orientation='horizontal', format='%.3f', ticks=ticks)
     set_axis(ax, letter=letter)
-    plt.tight_layout()
+    #plt.tight_layout()
+
+
+def make_plot_1D_pics(ax, k, est_val, tp, Fs, cut=9, title='Experimental data',
+                      cmap=plt.cm.bwr, letter='A', ylabel=True):
+
+    set_axis(ax, letter=letter)
+    npx.make_plot_spacetime(ax, k.estm_x, k.estm_y, est_val[cut,:,:], Fs,
+                            title=title, cmap=cmap, ylabel=ylabel)
+    if letter == 'D':
+        for lvl, name in zip([-500,-850,-2000], ['II/III', 'IV', 'V/VI']):
+            plt.axhline(lvl, ls='--', color='grey')
+            plt.text(340, lvl+20, name)
+    elif letter == 'C':
+        plt.axvline(tp/Fs*1000, ls='--', color ='grey', lw=2)
+    
+    plt.xlim(250, 400)
+    plt.xticks([250, 300, 350, 400], [-50, 0, 50, 100])
+    
+    #plt.tight_layout()
+    plt.savefig('figure_1D_pics', dpi=300)
 
 
-def make_figure_spacetime(val_pots, val_csd, cut=13, titl1='POT', titl2='CSD',
-                          fig_title='Pot and CSD in time'):
-    fig = plt.figure(figsize=(12, 8))
-    ax2 = plt.subplot(121)
-    make_plot_spacetime(ax2, val_pots, cut=cut,
+def make_figure_spacetime(val_pots_m, val_csd_m, kcsd_obj, val_pots_e, val_csd_e, tp, Fs, cut1=13, cut2=15,
+                          titl1='POT', titl2='CSD', fig_title='Pot and CSD in time'):
+    fig = plt.figure(figsize=(16, 8))
+    ax2 = plt.subplot(141)
+    make_plot_spacetime(ax2, val_pots_m, cut=cut1,
               title=titl1, cmap=plt.cm.PRGn, letter='A')
-    ax1 = plt.subplot(122)
-    make_plot_spacetime(ax1, val_csd, cut=cut, 
-              title=titl2, cmap=plt.cm.bwr, letter='B')
+    ax1 = plt.subplot(142)
+    make_plot_spacetime(ax1, val_csd_m, cut=cut1, 
+              title=titl2, cmap=plt.cm.bwr, letter='B', ylabel=False)
+    ax3 = plt.subplot(143)
+    make_plot_1D_pics(ax3, kcsd_obj, val_pots_e, tp, Fs, cut=cut2, title=titl1, cmap=plt.cm.PRGn, letter='C', ylabel=False)
+    ax4 = plt.subplot(144)
+    make_plot_1D_pics(ax4, kcsd_obj, val_csd_e, tp, Fs, cut=cut2, title=titl2, cmap=plt.cm.bwr, letter='D', ylabel=False)
+    plt.tight_layout()
     fig.savefig(os.path.join(fig_title + '.png'), dpi=300)
 
 
@@ -62,14 +91,35 @@ def make_figure_spacetime(val_pots, val_csd, cut=13, titl1='POT', titl2='CSD',
     elec_pos_list, names_list = prepare_electrodes()
 
     pot_np = prepare_pots(elec_pos_list[1], names_list[1], h, pop_names, time_pts)
-    kcsd, est_pot, x, y, k = do_kcsd(pot_np, elec_pos_list[1][:, :2], -40, 40, -3500, 500)
+    kcsd_m, est_pot_m, x_m, y_m, k_m = do_kcsd(pot_np, elec_pos_list[1][:, :2], -40, 40, -3500, 500)
+    
+    lowpass = 0.5
+    highpass = 300
+    Fs = 30000
+    resamp = 12
+    tp= 760
+    
+    forfilt=np.load('npx_data.npy')
+    
+    [b,a] = butter(3, [lowpass/(Fs/2.0), highpass/(Fs/2.0)] ,btype = 'bandpass')
+    filtData = filtfilt(b,a, forfilt)
+    pots_resamp = filtData[:,::resamp]
+    pots = pots_resamp[:, :]
+    Fs=int(Fs/resamp)
+    
+    pots_for_csd = np.delete(pots, 191, axis=0)
+    ele_pos_def = npx.eles_to_coords(np.arange(384,0,-1))
+    ele_pos_for_csd = np.delete(ele_pos_def, 191, axis=0)
+    
+    k_e, est_csd_e, est_pots_e, ele_pos_e = npx.do_kcsd(ele_pos_for_csd, pots_for_csd, ele_limit = (0,320))
 
     time_pts_ds = int(time_pts/4)
-    cut = 9
+    cut1 = 9
     start_pt = 625  # 50 ms before the stimulus
     end_pt = 1375  # 250 ms after the stimulus
     # for cut in range(kcsd.shape[0]):
-    make_figure_spacetime(est_pot[:, :, start_pt:end_pt],
-                          kcsd[:, :, start_pt:end_pt], cut=cut,
+    make_figure_spacetime(est_pot_m[:, :, start_pt:end_pt],
+                          kcsd_m[:, :, start_pt:end_pt], k_e, est_pots_e, est_csd_e, tp, Fs, cut1=cut1, cut2=15,
                           titl1='Estimated LFP', titl2='Estimated CSD',
-                          fig_title=('Estimated POT and CSD in time 1stim '+ str(cut)))
+                          fig_title=('Estimated POT and CSD in time 1stim '))
+    #plot_1D_pics(k, est_csd, est_pots, tp, 15)

figures/npx/kCSD2D_reconstruction_from_npx.py

@@ -6,17 +6,19 @@
 from figure_properties import *
 plt.close('all')
 #%%
-def make_plot_spacetime(ax, xx, yy, zz, title='True CSD', cmap=cm.bwr_r, ymin=0, ymax=10000):
+def make_plot_spacetime(ax, xx, yy, zz, Fs, title='True CSD', cmap=cm.bwr_r, ymin=0, ymax=10000, ylabel=True):
     im = ax.imshow(zz,extent=[0, zz.shape[1]/Fs*1000,-3500, 500], aspect='auto',
                    vmax = 1*zz.max(),vmin = -1*zz.max(), cmap=cmap)
     ax.set_xlabel('Time (ms)')
-    ax.set_ylabel('Y ($\mu$m)')
+    if ylabel:
+        ax.set_ylabel('Y ($\mu$m)')
     if 'Pot' in title: ax.set_ylabel('Y ($\mu$m)')
     ax.set_title(title)
+    ticks = np.linspace(-zz.max(), zz.max(), 3, endpoint=True)
     if 'CSD' in title:
-        plt.colorbar(im, orientation='vertical', format='%.2f', ticks = [-0.01,0,0.01])
+        plt.colorbar(im, orientation='horizontal', format='%.2f', ticks = ticks)
     else:
-        plt.colorbar(im, orientation='vertical', format='%.1f', ticks = [-0.6,0,0.6]) 
+        plt.colorbar(im, orientation='horizontal', format='%.1f', ticks = ticks) 
     # plt.gca().invert_yaxis()
 
 def make_plot(ax, xx, yy, zz, title='True CSD', cmap=cm.bwr):
@@ -57,13 +59,13 @@ def eles_to_coords(eles):
     ys = eles_to_ycoord(eles)
     return np.array((xs, ys)).T
 
-def plot_1D_pics(k, est_csd, est_pots, tp, cut=9):
+def plot_1D_pics(k, est_csd, est_pots, tp, Fs, cut=9):
     plt.figure(figsize=(12, 8))
     # plt.suptitle('plane: '+str(k.estm_x[cut,0])+' $\mu$m '+' $\lambda$ : '+str(k.lambd)+
                  # '  R: '+ str(k.R))
     ax1 = plt.subplot(122)
     set_axis(ax1, -0.05, 1.05, letter= 'D')
-    make_plot_spacetime(ax1, k.estm_x, k.estm_y, est_csd[cut,:,:], 
+    make_plot_spacetime(ax1, k.estm_x, k.estm_y, est_csd[cut,:,:], Fs,
               title='Estimated CSD', cmap='bwr')
     for lvl, name in zip([-500,-850,-2000], ['II/III', 'IV', 'V/VI']):
         plt.axhline(lvl, ls='--', color='grey')
@@ -80,7 +82,7 @@ def plot_1D_pics(k, est_csd, est_pots, tp, cut=9):
     plt.tight_layout()
     plt.savefig('figure_1D_pics', dpi=300)
 
-def plot_2D_pics(k, est_csd, est_pots, tp, cut, save=0):
+def plot_2D_pics(k, est_csd, est_pots, tp, Fs, cut, save=0):
     plt.figure(figsize=(12, 8))
     ax1 = plt.subplot(122)
     set_axis(ax1, -0.05, 1.05, letter= 'B')
@@ -127,5 +129,5 @@ def do_kcsd(ele_pos_for_csd, pots_for_csd, ele_limit):
 
     k, est_csd, est_pots, ele_pos = do_kcsd(ele_pos_for_csd, pots_for_csd, ele_limit = (0,320))
 
-    plot_1D_pics(k, est_csd, est_pots, tp, 15) 
-    plot_2D_pics(k, est_csd, est_pots, tp=tp, cut=15)
+    plot_1D_pics(k, est_csd, est_pots, tp, Fs, 15) 
+    plot_2D_pics(k, est_csd, est_pots, tp, Fs, cut=15)