implemented electrode extraction from metadata and added kCSD plot

Author: danek8317

figures/npx/dan_kCSD_from_npx.py

@@ -0,0 +1,259 @@
+import numpy as np
+from kcsd import KCSD2D
+from pathlib import Path
+from openpyxl import load_workbook
+import DemoReadSGLXData.readSGLX as readSGLX
+# from DemoReadSGLXData.readSGLX import readMeta, SampRate, makeMemMapRaw,
+#  GainCorrectIM, GainCorrectNI, ExtractDigital
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+# from matplotlib import gridspec
+
+
+def make_plot(xx, yy, zz, title='True CSD', cmap=cm.bwr):
+    fig = plt.figure(figsize=(7, 7))
+    ax = plt.subplot(111)
+    ax.set_aspect('equal')
+    t_max = np.max(np.abs(zz))
+    levels = np.linspace(-1 * t_max, t_max, 32)
+    im = ax.contourf(xx, yy, zz, levels=levels, cmap=cmap)
+    ax.set_xlabel('X (mm)')
+    ax.set_ylabel('Y (mm)')
+    ax.set_title(title)
+    ticks = np.linspace(-1 * t_max, t_max, 3, endpoint=True)
+    plt.colorbar(im, orientation='horizontal', format='%.2f', ticks=ticks)
+    return ax
+
+
+# Specific to Ewas experimental setup
+def load_chann_map():
+    book = load_workbook('NP_do_map.xlsx')
+    sheet = book.get_sheet_by_name('sov12 sorted')
+    eleid = sheet['C3':'C386']
+    chanid = sheet['J3':'J386']
+    chan_ele_dict = {}
+    ele_chan_dict = {}
+    for e,c in zip(eleid, chanid):
+        chan_ele_dict[int(c[0].value)] = int(e[0].value)
+        ele_chan_dict[int(e[0].value)] = int(c[0].value)
+    return ele_chan_dict, chan_ele_dict
+
+
+def dan_fetch_electrodes(meta):
+    imroList = meta['imroTbl'].split(sep=')')
+    
+    # One entry for each channel plus header entry,
+    # plus a final empty entry following the last ')'
+    nChan = len(imroList) - 2
+    electrode = np.zeros(nChan, dtype=int)        # default type = float
+    channel = np.zeros(nChan, dtype=int)
+    bank = np.zeros(nChan, dtype=int)
+    for i in range(0, nChan):
+        currList = imroList[i+1].split(sep=' ')
+        channel[i] = int(currList[0][1:])
+        bank[i] = int(currList[1])
+        # reference_electrode[i] = currList[2] 
+    
+    # Channel N => Electrode (1+N+384*A), where N = 0:383, A=0:2
+    electrode = 1 + channel + 384 * bank
+    return(electrode, channel)
+    
+
+def fetch_channels(eles):
+    chans = []
+    exist_ele = []
+    for ii in eles:
+        try:
+            chans.append(ele_chan_dict[ii])
+            exist_ele.append(ii)
+        except KeyError:
+            print('Not recording from ele', ii)
+    return chans, exist_ele
+
+def eles_to_rows(eles):
+    rows = []
+    for ele in eles:
+        rows.append(np.int(np.ceil(ele/2)))
+    return rows
+
+def eles_to_ycoord(eles):
+    rows = eles_to_rows(eles)
+    y_coords = []
+    for ii in rows:
+        y_coords.append(int((480 - ii)*20))
+    return y_coords
+
+def eles_to_xcoord(eles):
+    x_coords = []
+    for ele in eles:
+        off = ele%4
+        if off == 1:
+            x_coords.append(-24)
+        elif off == 2:
+            x_coords.append(8)
+        elif off == 3:
+            x_coords.append(-8)
+        else:
+            x_coords.append(24)
+    return x_coords
+
+def eles_to_coords(eles):
+    xs = eles_to_xcoord(eles)
+    ys = eles_to_ycoord(eles)
+    return np.array((xs, ys)).T
+
+
+# File with the data
+# old
+# binFullPath = Path('./data/08_refGND_APx500_LFPx125_ApfiltON_corr_banks_stim50V_g0_t0.imec0.lf.bin')
+# Daniel
+binFullPath = Path('/mnt/zasoby/data/neuropixel/Neuropixel data from Ewa Kublik/SOV_12/data/08_refGND_APx500_LFPx125_ApfiltON_corr_banks_stim50V_g0_t0.imec0.lf.bin')
+# Chaitanya
+# binFullPath = Path('/home/chaitanya/LFP/SOV_12/data/08_refGND_APx500_LFPx125_ApfiltON_corr_banks_stim50V_g0_t0.imec0.lf.bin')
+
+meta = readSGLX.readMeta(binFullPath)
+sRate = readSGLX.SampRate(meta)
+
+tStart, tEnd = 500., 600.  # 0., 1.        # in seconds
+
+firstSamp = int(sRate*tStart)
+lastSamp = int(sRate*tEnd)
+
+
+# Return array of original channel IDs. As an example, suppose we want the
+# imec gain for the ith channel stored in the binary data. A gain array
+# can be obtained using ChanGainsIM(), but we need an original channel
+# index to do the lookup. Because you can selectively save channels, the
+# ith channel in the file isn't necessarily the ith acquired channel.
+# Use this function to convert from ith stored to original index.
+# Note that the SpikeGLX channels are 0 based.
+#
+# chans = readSGLX.OriginalChans(meta)
+
+electrodes, channels = dan_fetch_electrodes(meta)
+# for Ewa's initial file, channel 768 is SY
+# and it hould be removed - this has not been done yet
+# DANIEL
+
+
+
+
+
+# =============================================================================
+# # chanList = [0, 6, 9, 383]
+# # eleList = np.arange(769, 860)
+# eleList = np.arange(0, 959)
+# 
+# ele_chan_dict, chan_ele_dict = load_chann_map()
+# # print(ele_dict)
+# chanList, eleList = fetch_channels(eleList)
+# 
+# =============================================================================
+
+
+# Which digital word to read. 
+# For imec, there is only 1 digital word, dw = 0.
+# For NI, digital lines 0-15 are in word 0, lines 16-31 are in word 1, etc.
+dw = 0
+# Which lines within the digital word, zero-based
+# Note that the SYNC line for PXI 3B is stored in line 6.
+dLineList = [6]
+
+rawData = readSGLX.makeMemMapRaw(binFullPath, meta)
+selectData = rawData[channels, firstSamp:lastSamp+1]
+# digArray = readSGLX.ExtractDigital(rawData, firstSamp, lastSamp, dw, dLineList, meta)
+
+# convData is the potential in uV or mV
+if meta['typeThis'] == 'imec':
+    # apply gain correction and convert to uV
+    convData = 1e6*readSGLX.GainCorrectIM(selectData, channels, meta)
+else:
+    # apply gain correction and convert to mV
+    convData = 1e3*readSGLX.GainCorrectNI(selectData, channels, meta)
+
+tDat = np.arange(firstSamp, lastSamp+1)
+tDat = 1000*tDat/sRate      # plot time axis in msec
+
+
+
+ele_pos = eles_to_coords(electrodes)
+print(ele_pos)
+csd_at_time = 30.
+pots = []
+for ii, chann in enumerate(channels):
+    print(ii, chann)
+    pots.append(convData[ii, int(sRate*csd_at_time)])
+
+pots = np.array(pots)
+print(pots.shape)
+
+
+pots = pots.reshape((len(channels), 1))
+R_init = 5. # 0.3
+h = 20.   # 50
+sigma = 0.3
+k = KCSD2D(ele_pos, pots, h=h, sigma=sigma,
+           xmin=-35, xmax=35,
+           ymin=1100, ymax=2000,
+           # ymin=1000, ymax=10000,
+           gdx=10, gdy=10,
+           R_init=R_init, n_src_init=1000,
+           src_type='gauss')   # rest of the parameters are set at default
+k.cross_validate(Rs=np.logspace(-1., 1., 10), lambdas=None)
+# k.cross_validate(Rs=np.linspace(0.1, 1.001, 2), lambdas=None)
+# 2 -> 20
+
+
+est_csd = k.values('CSD')
+est_csd = est_csd.reshape(7, 90)
+est_pots = k.values('POT')
+est_pots = est_pots.reshape(7, 90)
+
+make_plot(k.estm_x, k.estm_y, est_csd[:, :],
+          title='Estimated CSD without CV', cmap=cm.bwr)
+
+make_plot(k.estm_x, k.estm_y, est_pots[:, :],
+          title='Estimated POT without CV', cmap=cm.PRGn)
+
+
+# # ax = plt.subplot(121)
+# # for ii, chan in enumerate(chanList):
+# #     ax.plot(tDat, convData[ii, :], label=str(chan)+' Ele'+str(chan_dict[chan]))
+# # plt.legend()
+# # ax = plt.subplot(122)
+# # for i in range(0, len(dLineList)):
+# #     ax.plot(tDat, digArray[i, :])
+
+# rowList = eles_to_rows(eleList)
+# num_rows = max(rowList) - min(rowList) + 1
+# print(num_rows)
+# fig = plt.figure(figsize=(4, num_rows))
+# gs = gridspec.GridSpec(nrows=num_rows, ncols=4, wspace=0, hspace=0)
+# all_maxy = -100
+# axs = []
+# for ii, chann in enumerate(chanList):
+#     ee = chan_ele_dict[chann]
+#     rr = eles_to_rows([ee])[0] - min(rowList) # last row first
+#     rr = num_rows - rr - 1
+#     print(rr, ee, num_rows-rr)
+#     off = ee%4
+#     if off == 0:
+#         ax = fig.add_subplot(gs[rr, 3])
+#     elif off == 1:
+#         ax = fig.add_subplot(gs[rr, 0])
+#     elif off == 2:
+#         ax = fig.add_subplot(gs[rr, 2])
+#     else:
+#         ax = fig.add_subplot(gs[rr, 1])
+#     ax.plot(tDat, convData[ii, :])
+#     all_maxy = max(all_maxy, max(convData[ii, :]))
+#     ax.spines['right'].set_visible(False)
+#     ax.spines['top'].set_visible(False)
+#     # ax.spines['left'].set_visible(False)
+#     # ax.set_yticklabels([])
+#     # ax.set_yticks([])
+#     ax.set_title('E('+str(ee)+')')
+#     axs.append(ax)
+# print(all_maxy)
+# plt.show()
+