Coverage for brainbox/quality/lfp_qc.py: 0%

1""" 

2Quality control on LFP 

3Based on code by Olivier Winter: 

4https://github.com/int-brain-lab/ibllib/blob/master/examples/ibllib/ephys_qc_raw.py 

5 

6adapted by Anne Urai 

7Instructions https://docs.google.com/document/d/1lBNcssodWdBILBN0PrPWi0te4f6I8H_Bb-7ESxQdv5U/edit# 

8""" 

9 

10from pathlib import Path 

11 

12import sys 

13import glob 

14import os 

15import matplotlib.pyplot as plt 

16import numpy as np 

17import seaborn as sns 

18 

19 

20from ibllib.ephys import ephysqc 

21import one.alf.io as alfio 

22# from IPython import embed as shell 

23 

24 

25def _plot_spectra(outpath, typ, savefig=True): 

26 ''' 

27 TODO document this function 

28 ''' 

29 

30 spec = alfio.load_object(outpath, 'ephysQcFreq' + typ.upper(), namespace='spikeglx') 

31 

32 # hack to ensure a single key name 

33 if 'power.probe_00' in spec.keys(): 

34 spec['power'] = spec.pop('power.probe_00') 

35 spec['freq'] = spec.pop('freq.probe_00') 

36 elif 'power.probe_01' in spec.keys(): 

37 spec['power'] = spec.pop('power.probe_01') 

38 spec['freq'] = spec.pop('freq.probe_01') 

39 

40 # plot 

41 sns.set_style("whitegrid") 

42 plt.figure(figsize=[9, 4.5]) 

43 ax = plt.axes() 

44 ax.plot(spec['freq'], 20 * np.log10(spec['power'] + 1e-14), 

45 linewidth=0.5, color=[0.5, 0.5, 0.5]) 

46 ax.plot(spec['freq'], 20 * np.log10(np.median(spec['power'] + 1e-14, axis=1)), label='median') 

47 ax.set_xlabel(r'Frequency (Hz)') 

48 ax.set_ylabel(r'dB rel to $V^2.$Hz$^{-1}$') 

49 if typ == 'ap': 

50 ax.set_ylim([-275, -125]) 

51 elif typ == 'lf': 

52 ax.set_ylim([-260, -60]) 

53 ax.legend() 

54 ax.set_title(outpath) 

55 if savefig: 

56 plt.savefig(outpath / (typ + '_spec.png'), dpi=150) 

57 print('saved figure to %s' % (outpath / (typ + '_spec.png'))) 

58 

59 

60def _plot_rmsmap(outpath, typ, savefig=True): 

61 ''' 

62 TODO document this function 

63 ''' 

64 

65 rmsmap = alfio.load_object(outpath, 'ephysQcTime' + typ.upper(), namespace='spikeglx') 

66 

67 # hack to ensure a single key name 

68 if 'times.probe_00' in rmsmap.keys(): 

69 rmsmap['times'] = rmsmap.pop('times.probe_00') 

70 rmsmap['rms'] = rmsmap.pop('rms.probe_00') 

71 elif 'times.probe_01' in rmsmap.keys(): 

72 rmsmap['times'] = rmsmap.pop('times.probe_01') 

73 rmsmap['rms'] = rmsmap.pop('rms.probe_01') 

74 

75 plt.figure(figsize=[12, 4.5]) 

76 axim = plt.axes([0.2, 0.1, 0.7, 0.8]) 

77 axrms = plt.axes([0.05, 0.1, 0.15, 0.8]) 

78 axcb = plt.axes([0.92, 0.1, 0.02, 0.8]) 

79 

80 axrms.plot(np.median(rmsmap['rms'], axis=0)[:-1] * 1e6, np.arange(1, rmsmap['rms'].shape[1])) 

81 axrms.set_ylim(0, rmsmap['rms'].shape[1]) 

82 

83 im = axim.imshow(20 * np.log10(rmsmap['rms'].T + 1e-15), aspect='auto', origin='lower', 

84 extent=[rmsmap['times'][0], rmsmap['times'][-1], 0, rmsmap['rms'].shape[1]]) 

85 axim.set_xlabel(r'Time (s)') 

86 axrms.set_ylabel(r'Channel Number') 

87 plt.colorbar(im, cax=axcb) 

88 if typ == 'ap': 

89 im.set_clim(-110, -90) 

90 axrms.set_xlim(100, 0) 

91 elif typ == 'lf': 

92 im.set_clim(-100, -60) 

93 axrms.set_xlim(500, 0) 

94 axim.set_xlim(0, 4000) 

95 axim.set_title(outpath) 

96 if savefig: 

97 plt.savefig(outpath / (typ + '_rms.png'), dpi=150) 

98 

99 

100# ============================== ### 

101# FIND THE RIGHT FILES, RUN AS SCRIPT 

102# ============================== ### 

103 

104if __name__ == '__main__': 

105 if len(sys.argv) != 2: 

106 print("Please give the folder path as an input argument!") 

107 else: 

108 outpath = Path(sys.argv[1]) # grab from command line input 

109 fbin = glob.glob(os.path.join(outpath, '*.lf.bin')) 

110 assert len(fbin) > 0 

111 print('fbin: %s' % fbin) 

112 # make sure you send a path for the time being and not a string 

113 ephysqc.extract_rmsmap(Path(fbin[0])) 

114 _plot_spectra(outpath, 'lf') 

115 _plot_rmsmap(outpath, 'lf')