Coverage for src/thunderfish/thunderlogger.py: 0%

1"""# thunderlogger 

2 

3Detect segments of interest in large data files and extract EOD waveforms. 

4""" 

5 

6import sys 

7import os 

8import glob 

9import argparse 

10import traceback 

11import datetime as dt 

12import numpy as np 

13import pandas as pd 

14from scipy.signal import butter, lfilter 

15from types import SimpleNamespace 

16import matplotlib.pyplot as plt 

17import matplotlib.dates as mdates 

18from thunderlab.configfile import ConfigFile 

19from thunderlab.dataloader import DataLoader 

20from thunderlab.tabledata import TableData, write_table_args 

21from thunderlab.eventdetection import hist_threshold 

22from .version import __version__, __year__ 

23from .eodanalysis import save_eod_waveform, save_wave_fish, save_pulse_fish 

24from .eodanalysis import save_wave_spectrum, save_pulse_spectrum, save_pulse_peaks 

25from .eodanalysis import load_eod_waveform, load_wave_fish, load_pulse_fish 

26from .eodanalysis import load_wave_spectrum, load_pulse_spectrum, load_pulse_peaks 

27from .eodanalysis import wave_similarity, pulse_similarity 

28from .thunderfish import configuration, save_configuration 

29from .thunderfish import detect_eods, remove_eod_files 

30 

31 

32def add_thunderlogger_config(cfg, detection_thresh='auto', 

33 default_thresh=0.002, thresh_fac=3.0, 

34 thresh_nbins=500): 

35 """Add parameters needed for by thunderlogger. 

36 

37 Parameters 

38 ---------- 

39 cfg: ConfigFile 

40 The configuration. 

41 detection_thresh: float or 'auto' 

42 Only data segments with standard deviation larger than this value 

43 are analyzed for EODs. If set to 'auto' a threshold is computed 

44 from all the data segments of a recording channel. 

45 default_thresh: float 

46 Threshold that is used if "detection_thresh" is set to "auto" and 

47 no data are available. 

48 thresh_fac: float 

49 The threshold for analysing data segments is set to the mean of the 

50 most-likely standard deviations plus this factor times the corresponding 

51 standard deviation. 

52 thresh_nbins: int 

53 The number of bins used to compute a histogram of the standard 

54 deviations of the data segments, from which the mean and standard 

55 deviation are estimated for automatically computing a threshold. 

56 """ 

57 cfg.add_section('Thunderlogger:') 

58 cfg.add('detectionThreshold', detection_thresh, '', 'Only analyse data segements with a standard deviation that is larger than this threshold. If set to "auto" compute threshold from all the standard deviations of a recording channel.') 

59 cfg.add('detectionThresholdDefault', default_thresh, '', 'Threshold that is used if "detectionThreshold" is set to "auto" and no data are available.') 

60 cfg.add('detectionThresholdStdFac', thresh_fac, '', 'An automatically computed threshold for analysing data segments is set to the mean of the most-likely standard deviations plus this factor times the corresponding standard deviation.') 

61 cfg.add('detectionThresholdNBins', thresh_nbins, '', 'The number of bins used to compute a histogram of the standard deviations of the data segments, from which the mean and standard deviation are estimated for automatically computing a threshold.') 

62 cfg.add('startTime', 'none', '', 'Provide a start time for the recordings overwriting the meta data of the data files (YYYY-mm-ddTHH:MM:SS format).') 

63 

64 

65def extract_eods(files, thresholds, stds_only, cfg, verbose, plot_level, 

66 thresh=0.002, max_deltaf=1.0, max_dist=0.00005, 

67 deltat_max=dt.timedelta(minutes=5), start_time=None): 

68 t0s = [] 

69 stds = None 

70 supra_thresh = None 

71 wave_fishes = None 

72 pulse_fishes = None 

73 if start_time is None: 

74 # XXX we should read this from the meta data: 

75 filename = os.path.splitext(os.path.basename(files[0]))[0] 

76 times = filename.split('-')[1] 

77 start_time = dt.datetime.strptime(times, '%Y%m%dT%H%M%S') 

78 toffs = start_time 

79 t1 = start_time 

80 unit = None 

81 for file in files: 

82 try: 

83 with DataLoader(file) as sf: 

84 # analyze: 

85 sys.stdout.write(file + ': ') 

86 unit = sf.unit 

87 if max_dist < 1.1/sf.samplerate: 

88 max_dist = 1.1/sf.samplerate 

89 window_size = cfg.value('windowSize') 

90 ndata = int(window_size * sf.samplerate) 

91 step = ndata//2 

92 b, a = butter(1, 10.0, 'hp', fs=sf.samplerate, output='ba') 

93 if stds is None: 

94 stds = [[] for c in range(sf.channels)] 

95 supra_thresh = [[] for c in range(sf.channels)] 

96 if wave_fishes is None: 

97 wave_fishes = [[] for c in range(sf.channels)] 

98 if pulse_fishes is None: 

99 pulse_fishes = [[] for c in range(sf.channels)] 

100 for k, data in enumerate(sf.blocks(ndata, step)): 

101 sys.stdout.write('.') 

102 sys.stdout.flush() 

103 t0 = toffs + dt.timedelta(seconds=k*step/sf.samplerate) 

104 t1 = t0 + dt.timedelta(seconds=ndata/sf.samplerate) 

105 t0s.append(t0) 

106 for channel in range(sf.channels): 

107 if thresholds: 

108 thresh = thresholds[channel] 

109 fdata = lfilter(b, a, data[:,channel] - np.mean(data[:ndata//20,channel])) 

110 sd = np.std(fdata) 

111 stds[channel].append(sd) 

112 supra_thresh[channel].append(1 if sd > thresh else 0) 

113 if stds_only: 

114 continue 

115 if sd > thresh: 

116 # clipping: 

117 min_clip = cfg.value('minClipAmplitude') 

118 if min_clip == 0.0: 

119 min_clip = cfg.value('minDataAmplitude') 

120 max_clip = cfg.value('maxClipAmplitude') 

121 if max_clip == 0.0: 

122 max_clip = cfg.value('maxDataAmplitude') 

123 name = file 

124 # detect EODs in the data: 

125 _, _, _, eod_props, mean_eods, spec_data, peak_data, _, _, _ = \ 

126 detect_eods(data[:,channel], sf.samplerate, 

127 min_clip, max_clip, 

128 name, verbose, plot_level, cfg) 

129 first_fish = True 

130 for props, eod, spec, peaks in zip(eod_props, mean_eods, 

131 spec_data, peak_data): 

132 fish = None 

133 fish_deltaf = 100000.0 

134 if props['type'] == 'wave': 

135 for wfish in wave_fishes[channel]: 

136 deltaf = np.abs(wfish.props['EODf'] - props['EODf']) 

137 if deltaf < fish_deltaf: 

138 fish_deltaf = deltaf 

139 fish = wfish 

140 if fish_deltaf > max_deltaf: 

141 fish = None 

142 peaks = None 

143 else: 

144 fish_dist = 10000.0 

145 for pfish in pulse_fishes[channel]: 

146 ddist = np.abs(pfish.props['P1-P1-dist'] - 

147 props['P1-P1-dist']) 

148 if ddist < fish_dist: 

149 fish_dist = ddist 

150 fish_deltaf = np.abs(pfish.props['EODf'] - 

151 props['EODf']) 

152 fish = pfish 

153 if fish_dist > max_dist or \ 

154 fish_deltaf > max_deltaf: 

155 fish = None 

156 spec = None 

157 if fish is not None and \ 

158 t0 - fish.times[-1][1] < deltat_max: 

159 if fish.times[-1][1] >= t0 and \ 

160 np.abs(fish.times[-1][2] - props['EODf']) < 0.5 and \ 

161 fish.times[-1][3] == channel and \ 

162 fish.times[-1][4] == file: 

163 fish.times[-1][1] = t1 

164 else: 

165 fish.times.append([t0, t1, props['EODf'], channel, file]) 

166 if props['p-p-amplitude'] > fish.props['p-p-amplitude']: 

167 fish.props = props 

168 fish.waveform = eod 

169 fish.spec = spec 

170 fish.peaks = peaks 

171 else: 

172 new_fish = SimpleNamespace(props=props, 

173 waveform=eod, 

174 spec=spec, 

175 peaks=peaks, 

176 times=[[t0, t1, props['EODf'], channel, file]]) 

177 if props['type'] == 'pulse': 

178 pulse_fishes[channel].append(new_fish) 

179 else: 

180 wave_fishes[channel].append(new_fish) 

181 if first_fish: 

182 sys.stdout.write('\n ') 

183 first_fish = False 

184 sys.stdout.write('%6.1fHz %5s-fish @ %s\n ' % 

185 (props['EODf'], props['type'], 

186 t0.strftime('%Y-%m-%dT%H:%M:%S'))) 

187 toffs += dt.timedelta(seconds=len(sf)/sf.samplerate) 

188 sys.stdout.write('\n') 

189 sys.stdout.flush() 

190 except EOFError as error: 

191 # XXX we need to update toffs by means of the metadata of the next file! 

192 sys.stdout.write(file + ': ' + str(error) + '\n') 

193 if pulse_fishes is not None and len(pulse_fishes) > 0: 

194 pulse_fishes = [[pulse_fishes[c][i] for i in 

195 np.argsort([fish.props['EODf'] for fish in pulse_fishes[c]])] 

196 for c in range(len(pulse_fishes))] 

197 if wave_fishes is not None and len(wave_fishes) > 0: 

198 wave_fishes = [[wave_fishes[c][i] for i in 

199 np.argsort([fish.props['EODf'] for fish in wave_fishes[c]])] 

200 for c in range(len(wave_fishes))] 

201 return pulse_fishes, wave_fishes, start_time, toffs, t0s, stds, supra_thresh, unit 

202 

203 

204def save_times(times, idx, output_basename, name, **kwargs): 

205 td = TableData() 

206 td.append('index', '', '%d', [idx] * len(times)) 

207 td.append('tstart', '', '%s', 

208 [t[0].strftime('%Y-%m-%dT%H:%M:%S') for t in times]) 

209 td.append('tend', '', '%s', 

210 [t[1].strftime('%Y-%m-%dT%H:%M:%S') for t in times]) 

211 if len(times[0]) > 2: 

212 td.append('EODf', 'Hz', '%.1f', [t[2] for t in times]) 

213 td.append('device', '', '%s', 

214 [name for t in times]) 

215 if len(times[0]) > 2: 

216 td.append('channel', '', '%d', [t[3] for t in times]) 

217 td.append('file', '', '%s', [t[4] for t in times]) 

218 fp = output_basename + '-times' 

219 if idx is not None: 

220 fp += '-%d' % idx 

221 td.write(fp, **kwargs) 

222 

223 

224def load_times(file_path): 

225 data = TableData(file_path).data_frame() 

226 data['index'] = data['index'].astype(int) 

227 data['tstart'] = pd.to_datetime(data['tstart']) 

228 data['tstart'] = pd.Series(data['tstart'].dt.to_pydatetime(), dtype=object) 

229 data['tend'] = pd.to_datetime(data['tend']) 

230 data['tend'] = pd.Series(data['tend'].dt.to_pydatetime(), dtype=object) 

231 if 'channel' in data: 

232 data['channel'] = data['channel'].astype(int) 

233 return data 

234 

235 

236def save_power(times, stds, supra_thresh, unit, output_basename, **kwargs): 

237 td = TableData() 

238 td.append('index', '', '%d', list(range(len(times)))) 

239 td.append('time', '', '%s', 

240 [t.strftime('%Y-%m-%dT%H:%M:%S') for t in times]) 

241 for c, (std, thresh) in enumerate(zip(stds, supra_thresh)): 

242 td.append('channel%d'%c, unit, '%g', std) 

243 td.append('thresh%d'%c, '', '%d', thresh) 

244 fp = output_basename + '-stdevs' 

245 td.write(fp, **kwargs) 

246 

247 

248def load_power(file_path, start_time=None): 

249 base = os.path.basename(file_path) 

250 device = base[0:base.find('-stdevs')] 

251 data = TableData(file_path) 

252 times = [] 

253 for row in range(data.rows()): 

254 times.append(dt.datetime.strptime(data[row,'time'], 

255 '%Y-%m-%dT%H:%M:%S')) 

256 if start_time is not None: 

257 deltat = start_time - times[0] 

258 for k in range(len(times)): 

259 times[k] += deltat 

260 channels = (data.columns()-2)//2 

261 stds = np.zeros((len(times), channels)) 

262 supra_thresh = np.zeros((len(times), channels), dtype=int) 

263 for c in range(channels): 

264 stds[:,c] = data[:,'channel%d'%c] 

265 supra_thresh[:,c] = data[:,'thresh%d'%c] 

266 return np.array(times), stds, supra_thresh, device 

267 

268 

269def save_data(output_folder, name, pulse_fishes, wave_fishes, 

270 tstart, tend, t0s, stds, supra_thresh, unit, cfg): 

271 output_basename = os.path.join(output_folder, name) 

272 if pulse_fishes is not None: 

273 for c in range(len(pulse_fishes)): 

274 out_path = output_basename + '-c%d' % c 

275 idx = 0 

276 # pulse fish: 

277 pulse_props = [] 

278 for fish in pulse_fishes[c]: 

279 save_eod_waveform(fish.waveform, unit, idx, out_path, 

280 **write_table_args(cfg)) 

281 if fish.peaks is not None: 

282 save_pulse_peaks(fish.peaks, unit, idx, out_path, 

283 **write_table_args(cfg)) 

284 save_times(fish.times, idx, out_path, name, 

285 **write_table_args(cfg)) 

286 pulse_props.append(fish.props) 

287 pulse_props[-1]['index'] = idx 

288 idx += 1 

289 save_pulse_fish(pulse_props, unit, out_path, 

290 **write_table_args(cfg)) 

291 # wave fish: 

292 wave_props = [] 

293 if wave_fishes is not None: 

294 for fish in wave_fishes[c]: 

295 save_eod_waveform(fish.waveform, unit, idx, out_path, 

296 **write_table_args(cfg)) 

297 if fish.spec is not None: 

298 save_wave_spectrum(fish.spec, unit, idx, out_path, 

299 **write_table_args(cfg)) 

300 save_times(fish.times, idx, out_path, name, 

301 **write_table_args(cfg)) 

302 wave_props.append(fish.props) 

303 wave_props[-1]['index'] = idx 

304 idx += 1 

305 save_wave_fish(wave_props, unit, out_path, 

306 **write_table_args(cfg)) 

307 # recording time window: 

308 save_times([(tstart, tend)], None, output_basename, name, 

309 **write_table_args(cfg)) 

310 # signal power: 

311 if stds is not None and len(stds) > 0: 

312 save_power(t0s, stds, supra_thresh, unit, output_basename, 

313 **write_table_args(cfg)) 

314 

315 

316def load_data(files, start_time=None): 

317 all_files = [] 

318 for file in files: 

319 if os.path.isdir(file): 

320 all_files.extend(glob.glob(os.path.join(file, '*fish*'))) 

321 else: 

322 all_files.append(file) 

323 pulse_fishes = [] 

324 wave_fishes = [] 

325 channels = set() 

326 for file in all_files: 

327 if 'pulse' in os.path.basename(file): 

328 pulse_props = load_pulse_fish(file) 

329 base_file, ext = os.path.splitext(file) 

330 base_file = base_file[:base_file.rfind('pulse')] 

331 count = 0 

332 for props in pulse_props: 

333 idx = props['index'] 

334 waveform, unit = \ 

335 load_eod_waveform(base_file + 'eodwaveform-%d'%idx + ext) 

336 times = load_times(base_file + 'times-%d'%idx + ext) 

337 times['index'] = idx 

338 fish = SimpleNamespace(props=props, 

339 waveform=waveform, 

340 unit=unit, 

341 times=times) 

342 for i, t in times.iterrows(): 

343 channels.add((t['device'], t['channel'])) 

344 try: 

345 peaks, unit = \ 

346 load_pulse_peaks(base_file + 'pulsepeaks-%d'%idx + ext) 

347 fish.peaks = peaks 

348 except FileNotFoundError: 

349 fish.peaks = None 

350 pulse_fishes.append(fish) 

351 count += 1 

352 #if count > 300: # XXX REMOVE 

353 # break 

354 elif 'wave' in os.path.basename(file): 

355 wave_props = load_wave_fish(file) 

356 base_file, ext = os.path.splitext(file) 

357 base_file = base_file[:base_file.rfind('wave')] 

358 count = 0 

359 for props in wave_props: 

360 idx = props['index'] 

361 waveform, unit = \ 

362 load_eod_waveform(base_file + 'eodwaveform-%d'%idx + ext) 

363 times = load_times(base_file + 'times-%d'%idx + ext) 

364 times['index'] = idx 

365 fish = SimpleNamespace(props=props, 

366 waveform=waveform, 

367 unit=unit, 

368 times=times) 

369 for i, t in times.iterrows(): 

370 channels.add((t['device'], t['channel'])) 

371 try: 

372 spec, unit = \ 

373 load_wave_spectrum(base_file + 'wavespectrum-%d'%idx + ext) 

374 fish.spec = spec 

375 except FileNotFoundError: 

376 fish.spec = None 

377 wave_fishes.append(fish) 

378 count += 1 

379 #if count > 300: # XXX REMOVE 

380 # break 

381 base_file = base_file[:base_file.rfind('-c')+1] 

382 times = load_times(base_file + 'times' + ext) 

383 tstart = times.tstart[0] 

384 tend = times.tend[0] 

385 if start_time is not None: 

386 deltat = start_time - tstart 

387 tstart += deltat 

388 tend += deltat 

389 for fish in pulse_fishes + wave_fishes: 

390 fish.times.tstart += deltat 

391 fish.times.tend += deltat 

392 return pulse_fishes, wave_fishes, tstart, tend, sorted(channels) 

393 

394 

395def plot_signal_power(times, stds, supra_threshs, devices, thresholds, 

396 title, output_folder): 

397 plt.rcParams['axes.xmargin'] = 0 

398 plt.rcParams['axes.ymargin'] = 0 

399 n = 0 

400 for s in stds: 

401 n += s.shape[1] 

402 h = n*4.0 

403 t = 0.8 if title else 0.1 

404 fig, axs = plt.subplots(n, 1, figsize=(16/2.54, h/2.54), 

405 sharex=True, sharey=True) 

406 fig.subplots_adjust(left=0.1, right=0.99, top=1-t/h, bottom=1.6/h, 

407 hspace=0) 

408 i = 0 

409 for time, cstds, threshs, device in zip(times, stds, supra_threshs, devices): 

410 for c, (std, thresh) in enumerate(zip(cstds.T, threshs.T)): 

411 ax = axs[i] 

412 ax.plot(time, std) 

413 if thresholds: 

414 ax.axhline(thresholds[i], color='k', lw=0.5) 

415 elif len(std[thresh<1]) > 0: 

416 thresh = np.max(std[thresh<1]) 

417 ax.axhline(thresh, color='k', lw=0.5) 

418 #stdm = np.ma.masked_where(thresh < 1, std) 

419 #ax.plot(time, stdm) 

420 ax.set_yscale('log') 

421 #ax.set_ylim(bottom=0) 

422 ax.set_ylabel('%s-c%d' % (device, c)) 

423 i += 1 

424 if title: 

425 axs[0].set_title(title) 

426 axs[-1].xaxis.set_major_formatter(mdates.DateFormatter('%b %d %Hh')) 

427 plt.setp(axs[-1].get_xticklabels(), ha='right', 

428 rotation=30, rotation_mode='anchor') 

429 fig.savefig(os.path.join(output_folder, 'signalpowers.pdf')) 

430 plt.show() 

431 

432 

433def merge_fish(pulse_fishes, wave_fishes, 

434 max_noise=0.1, max_deltaf=10.0, max_dist=0.0002, max_rms=0.3): 

435 pulse_eods = [] 

436 for i in np.argsort([fish.props['P2-P1-dist'] for fish in pulse_fishes]): 

437 if pulse_fishes[i].props['noise'] > max_noise: 

438 continue 

439 if pulse_fishes[i].props['P2-P1-dist'] <= 0.0: 

440 continue 

441 if len(pulse_eods) > 0 and \ 

442 np.abs(pulse_fishes[i].props['P2-P1-dist'] - pulse_eods[-1].props['P2-P1-dist']) <= max_dist and \ 

443 pulse_similarity(pulse_fishes[i].waveform, pulse_eods[-1].waveform, 4) < max_rms: 

444 pulse_eods[-1].times.append(pulse_fishes[i].times) 

445 if not hasattr(pulse_eods[-1], 'othereods'): 

446 pulse_eods[-1].othereods = [pulse_eods[-1].waveform] 

447 pulse_eods[-1].othereods.append(pulse_fishes[i].waveform) 

448 if pulse_fishes[i].props['p-p-amplitude'] > pulse_eods[-1].props['p-p-amplitude']: 

449 pulse_eods[-1].waveform = pulse_fishes[i].waveform 

450 pulse_eods[-1].props = pulse_fishes[i].props 

451 continue 

452 pulse_eods.append(pulse_fishes[i]) 

453 pulse_eods = [pulse_eods[i] for i in np.argsort([fish.props['EODf'] for fish in pulse_eods])] 

454 

455 wave_eods = [] 

456 for i in np.argsort([fish.props['EODf'] for fish in wave_fishes]): 

457 if wave_fishes[i].props['noise'] > max_noise: 

458 continue 

459 if len(wave_eods) > 0 and \ 

460 np.abs(wave_fishes[i].props['EODf'] - wave_eods[-1].props['EODf']) < max_deltaf and \ 

461 wave_similarity(wave_fishes[i].waveform, wave_eods[-1].waveform, 

462 wave_fishes[i].props['EODf'], wave_eods[-1].props['EODf']) < max_rms: 

463 wave_eods[-1].times.append(wave_fishes[i].times) 

464 if not hasattr(wave_eods[-1], 'othereods'): 

465 wave_eods[-1].othereods = [] 

466 wave_eods[-1].othereods.append(wave_fishes[i].waveform) 

467 if wave_fishes[i].props['p-p-amplitude'] > wave_eods[-1].props['p-p-amplitude']: 

468 wave_eods[-1].waveform = wave_fishes[i].waveform 

469 wave_eods[-1].props = wave_fishes[i].props 

470 continue 

471 wave_eods.append(wave_fishes[i]) 

472 return pulse_eods, wave_eods 

473 

474 

475def plot_eod_occurances(pulse_fishes, wave_fishes, tstart, tend, 

476 channels, output_folder): 

477 channel_colors = ['#2060A7', '#40A787', '#478010', '#F0D730', 

478 '#C02717', '#873770', '#008797', '#007030', 

479 '#AAB71B', '#F78017', '#D03050', '#53379B'] 

480 plt.rcParams['axes.facecolor'] = 'none' 

481 plt.rcParams['axes.xmargin'] = 0 

482 plt.rcParams['axes.ymargin'] = 0.05 

483 plt.rcParams['font.family'] = 'sans-serif' 

484 n = len(pulse_fishes) + len(wave_fishes) 

485 h = n*2.5 + 2.5 + 0.3 

486 fig, axs = plt.subplots(n, 2, squeeze=False, figsize=(16/2.54, h/2.54), 

487 gridspec_kw=dict(width_ratios=(1,2))) 

488 fig.subplots_adjust(left=0.02, right=0.97, top=1-0.3/h, bottom=2.2/h, 

489 hspace=0.2) 

490 pi = 0 

491 prev_xscale = 0.0 

492 for ax, fish in zip(axs, pulse_fishes + wave_fishes): 

493 # EOD waveform: 

494 ax[0].spines['left'].set_visible(False) 

495 ax[0].spines['right'].set_visible(False) 

496 ax[0].spines['top'].set_visible(False) 

497 ax[0].spines['bottom'].set_visible(False) 

498 ax[0].xaxis.set_visible(False) 

499 ax[0].yaxis.set_visible(False) 

500 time = 1000.0 * fish.waveform[:,0] 

501 #ax[0].plot([time[0], time[-1]], [0.0, 0.0], 

502 # zorder=-10, lw=1, color='#AAAAAA') 

503 if hasattr(fish, 'othereods'): 

504 for eod in fish.othereods: 

505 ax[0].plot(1000.0*eod[:,0], eod[:,1], 

506 zorder=5, lw=1, color='#AAAAAA') 

507 ax[0].plot(time, fish.waveform[:,1], 

508 zorder=10, lw=2, color='#C02717') 

509 ax[0].text(0.0, 1.0, '%.0f\u2009Hz' % fish.props['EODf'], 

510 transform=ax[0].transAxes, va='baseline', zorder=20) 

511 if fish.props['type'] == 'wave': 

512 lim = 750.0/fish.props['EODf'] 

513 ax[0].set_xlim([-lim, +lim]) 

514 tmax = lim 

515 else: 

516 ax[0].set_xlim(time[0], time[-1]) 

517 tmax = time[-1] 

518 xscale = 1.0 

519 if tmax < 1.0: 

520 xscale = 0.5 

521 elif tmax > 10.0: 

522 xscale = 5.0 

523 ymin = np.min(fish.waveform[:,1]) 

524 ymax = np.max(fish.waveform[:,1]) 

525 ax[0].plot((tmax-xscale, tmax), (ymin - 0.04*(ymax-ymin),)*2, 

526 'k', lw=3, clip_on=False, zorder=0) 

527 if ax[0] is axs[-1,0] or xscale != prev_xscale: 

528 if xscale < 1.0: 

529 ax[0].text(tmax-0.5*xscale, ymin - 0.1*(ymax-ymin), 

530 '%.0f\u2009\u00b5s' % (1000.0*xscale), 

531 ha='center', va='top', zorder=0) 

532 else: 

533 ax[0].text(tmax-0.5*xscale, ymin - 0.1*(ymax-ymin), 

534 '%.0f\u2009ms' % xscale, 

535 ha='center', va='top', zorder=0) 

536 prev_xscale = xscale 

537 # time bar: 

538 ax[1].xaxis.set_major_formatter(mdates.DateFormatter('%b %d %Hh')) 

539 min_eodf = 10000 

540 max_eodf = 0 

541 for index, time in fish.times.iterrows(): 

542 if time['EODf'] < min_eodf: 

543 min_eodf = time['EODf'] 

544 if time['EODf'] > max_eodf: 

545 max_eodf = time['EODf'] 

546 ax[1].plot([time['tstart'], time['tend']], 

547 [time['EODf'], time['EODf']], 

548 lw=5, color=channel_colors[channels.index((time['device'], time['channel']))%len(channel_colors)]) 

549 if max_eodf > min_eodf + 10.0: 

550 ax[1].text(0.0, 1.0, '%.0f \u2013 %.0f\u2009Hz' % (min_eodf, max_eodf), 

551 transform=ax[1].transAxes, va='baseline', zorder=20) 

552 ax[1].set_xlim(tstart, tend) 

553 ax[1].spines['left'].set_visible(False) 

554 ax[1].spines['right'].set_visible(False) 

555 ax[1].yaxis.set_visible(False) 

556 ax[1].spines['top'].set_visible(False) 

557 if ax[1] is not axs[-1,1]: 

558 ax[1].spines['bottom'].set_visible(False) 

559 ax[1].xaxis.set_visible(False) 

560 else: 

561 plt.setp(ax[1].get_xticklabels(), ha='right', 

562 rotation=30, rotation_mode='anchor') 

563 fig.savefig(os.path.join(output_folder, 'eodwaveforms.pdf')) 

564 

565 

566def main(cargs=None): 

567 # config file name: 

568 cfgfile = __package__ + '.cfg' 

569 

570 # command line arguments: 

571 if cargs is None: 

572 cargs = sys.argv[1:] 

573 parser = argparse.ArgumentParser(add_help=False, 

574 description='Extract EOD waveforms of weakly electric fish from logger data.', 

575 epilog='version %s by Benda-Lab (2015-%s)' % (__version__, __year__)) 

576 parser.add_argument('-h', '--help', action='store_true', 

577 help='show this help message and exit') 

578 parser.add_argument('--version', action='version', version=__version__) 

579 parser.add_argument('-v', action='count', dest='verbose', default=0, 

580 help='verbosity level. Increase by specifying -v multiple times, or like -vvv') 

581 parser.add_argument('-V', action='count', dest='plot_level', default=0, 

582 help='level for debugging plots. Increase by specifying -V multiple times, or like -VVV') 

583 parser.add_argument('-c', dest='save_config', action='store_true', 

584 help='save configuration to file {0} after reading all configuration files'.format(cfgfile)) 

585 parser.add_argument('-f', dest='format', default='auto', type=str, 

586 choices=TableData.formats + ['py'], 

587 help='file format used for saving analysis results, defaults to the format specified in the configuration file or "csv"') 

588 parser.add_argument('-p', dest='save_plot', action='store_true', 

589 help='plot analyzed data') 

590 parser.add_argument('-m', dest='merge', action='store_true', 

591 help='merge similar EODs before plotting') 

592 parser.add_argument('-s', dest='stds_only', action='store_true', 

593 help='analyze or plot standard deviation of data only')