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 

14import matplotlib.pyplot as plt 

15import matplotlib.dates as mdates 

16 

17from scipy.signal import butter, lfilter 

18from types import SimpleNamespace 

19from thunderlab.configfile import ConfigFile 

20from thunderlab.dataloader import DataLoader 

21from thunderlab.tabledata import TableData, write_table_args 

22from thunderlab.eventdetection import hist_threshold 

23 

24from .version import __version__, __year__ 

25from .eodanalysis import save_eod_waveform, save_wave_fish, save_pulse_fish 

26from .eodanalysis import save_wave_spectrum, save_pulse_spectrum, save_pulse_phases 

27from .eodanalysis import load_eod_waveform, load_wave_fish, load_pulse_fish 

28from .eodanalysis import load_wave_spectrum, load_pulse_spectrum, load_pulse_phases 

29from .eodanalysis import wave_similarity, pulse_similarity 

30from .thunderfish import configuration, save_configuration 

31from .thunderfish import detect_eods, remove_eod_files 

32 

33 

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

35 default_thresh=0.002, thresh_fac=3.0, 

36 thresh_nbins=500): 

37 """Add parameters needed for by thunderlogger. 

38 

39 Parameters 

40 ---------- 

41 cfg: ConfigFile 

42 The configuration. 

43 detection_thresh: float or 'auto' 

44 Only data segments with standard deviation larger than this value 

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

46 from all the data segments of a recording channel. 

47 default_thresh: float 

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

49 no data are available. 

50 thresh_fac: float 

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

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

53 standard deviation. 

54 thresh_nbins: int 

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

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

57 deviation are estimated for automatically computing a threshold. 

58 """ 

59 cfg.add_section('Thunderlogger:') 

60 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.') 

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

62 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.') 

63 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.') 

64 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).') 

65 

66 

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

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

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

70 t0s = [] 

71 stds = None 

72 supra_thresh = None 

73 wave_fishes = None 

74 pulse_fishes = None 

75 if start_time is None: 

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

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

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

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

80 toffs = start_time 

81 t1 = start_time 

82 unit = None 

83 for file in files: 

84 try: 

85 with DataLoader(file) as sf: 

86 # analyze: 

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

88 unit = sf.unit 

89 if max_dist < 1.1/sf.rate: 

90 max_dist = 1.1/sf.rate 

91 window_size = cfg.value('windowSize') 

92 ndata = int(window_size * sf.rate) 

93 step = ndata//2 

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

95 if stds is None: 

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

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

98 if wave_fishes is None: 

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

100 if pulse_fishes is None: 

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

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

103 sys.stdout.write('.') 

104 sys.stdout.flush() 

105 t0 = toffs + dt.timedelta(seconds=k*step/sf.rate) 

106 t1 = t0 + dt.timedelta(seconds=ndata/sf.rate) 

107 t0s.append(t0) 

108 for channel in range(sf.channels): 

109 if thresholds: 

110 thresh = thresholds[channel] 

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

112 sd = np.std(fdata) 

113 stds[channel].append(sd) 

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

115 if stds_only: 

116 continue 

117 if sd > thresh: 

118 # clipping: 

119 min_clip = cfg.value('minClipAmplitude') 

120 if min_clip == 0.0: 

121 min_clip = cfg.value('minDataAmplitude') 

122 max_clip = cfg.value('maxClipAmplitude') 

123 if max_clip == 0.0: 

124 max_clip = cfg.value('maxDataAmplitude') 

125 name = file 

126 # detect EODs in the data: 

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

128 detect_eods(data[:,channel], sf.rate, 

129 min_clip, max_clip, 

130 name, verbose, plot_level, cfg) 

131 first_fish = True 

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

133 spec_data, peak_data): 

134 fish = None 

135 fish_deltaf = 100000.0 

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

137 for wfish in wave_fishes[channel]: 

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

139 if deltaf < fish_deltaf: 

140 fish_deltaf = deltaf 

141 fish = wfish 

142 if fish_deltaf > max_deltaf: 

143 fish = None 

144 peaks = None 

145 else: 

146 fish_dist = 10000.0 

147 for pfish in pulse_fishes[channel]: 

148 ddist = np.abs(pfish.props['peak-dist'] - 

149 props['peak-dist']) 

150 if ddist < fish_dist: 

151 fish_dist = ddist 

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

153 props['EODf']) 

154 fish = pfish 

155 if fish_dist > max_dist or \ 

156 fish_deltaf > max_deltaf: 

157 fish = None 

158 spec = None 

159 if fish is not None and \ 

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

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

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

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

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

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

166 else: 

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

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

169 fish.props = props 

170 fish.waveform = eod 

171 fish.spec = spec 

172 fish.peaks = peaks 

173 else: 

174 new_fish = SimpleNamespace(props=props, 

175 waveform=eod, 

176 spec=spec, 

177 peaks=peaks, 

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

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

180 pulse_fishes[channel].append(new_fish) 

181 else: 

182 wave_fishes[channel].append(new_fish) 

183 if first_fish: 

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

185 first_fish = False 

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

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

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

189 toffs += dt.timedelta(seconds=len(sf)/sf.rate) 

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

191 sys.stdout.flush() 

192 except EOFError as error: 

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

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

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

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

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

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

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

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

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

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

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

204 

205 

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

207 td = TableData() 

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

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

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

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

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

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

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

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

216 [name for t in times]) 

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

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

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

220 fp = output_basename + '-times' 

221 if idx is not None: 

222 fp += '-%d' % idx 

223 td.write(fp, **kwargs) 

224 

225 

226def load_times(file_path): 

227 data = TableData(file_path).data_frame() 

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

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

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

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

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

233 if 'channel' in data: 

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

235 return data 

236 

237 

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

239 td = TableData() 

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

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

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

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

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

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

246 fp = output_basename + '-stdevs' 

247 td.write(fp, **kwargs) 

248 

249 

250def load_power(file_path, start_time=None): 

251 base = os.path.basename(file_path) 

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

253 data = TableData(file_path) 

254 times = [] 

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

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

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

258 if start_time is not None: 

259 deltat = start_time - times[0] 

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

261 times[k] += deltat 

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

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

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

265 for c in range(channels): 

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

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

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

269 

270 

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

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

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

274 if pulse_fishes is not None: 

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

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

277 idx = 0 

278 # pulse fish: 

279 pulse_props = [] 

280 for fish in pulse_fishes[c]: 

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

282 **write_table_args(cfg)) 

283 if fish.peaks is not None: 

284 save_pulse_phases(fish.peaks, unit, idx, out_path, 

285 **write_table_args(cfg)) 

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

287 **write_table_args(cfg)) 

288 pulse_props.append(fish.props) 

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

290 idx += 1 

291 save_pulse_fish(pulse_props, unit, out_path, 

292 **write_table_args(cfg)) 

293 # wave fish: 

294 wave_props = [] 

295 if wave_fishes is not None: 

296 for fish in wave_fishes[c]: 

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

298 **write_table_args(cfg)) 

299 if fish.spec is not None: 

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

301 **write_table_args(cfg)) 

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

303 **write_table_args(cfg)) 

304 wave_props.append(fish.props) 

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

306 idx += 1 

307 save_wave_fish(wave_props, unit, out_path, 

308 **write_table_args(cfg)) 

309 # recording time window: 

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

311 **write_table_args(cfg)) 

312 # signal power: 

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

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

315 **write_table_args(cfg)) 

316 

317 

318def load_data(files, start_time=None): 

319 all_files = [] 

320 for file in files: 

321 if os.path.isdir(file): 

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

323 else: 

324 all_files.append(file) 

325 pulse_fishes = [] 

326 wave_fishes = [] 

327 channels = set() 

328 for file in all_files: 

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

330 pulse_props = load_pulse_fish(file) 

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

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

333 count = 0 

334 for props in pulse_props: 

335 idx = props['index'] 

336 waveform, unit = \ 

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

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

339 times['index'] = idx 

340 fish = SimpleNamespace(props=props, 

341 waveform=waveform, 

342 unit=unit, 

343 times=times) 

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

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

346 try: 

347 peaks, unit = \ 

348 load_pulse_phases(base_file + 'pulsephases-%d'%idx + ext) 

349 fish.peaks = peaks 

350 except FileNotFoundError: 

351 fish.peaks = None 

352 pulse_fishes.append(fish) 

353 count += 1 

354 #if count > 300: # XXX REMOVE 

355 # break 

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

357 wave_props = load_wave_fish(file) 

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

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

360 count = 0 

361 for props in wave_props: 

362 idx = props['index'] 

363 waveform, unit = \ 

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

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

366 times['index'] = idx 

367 fish = SimpleNamespace(props=props, 

368 waveform=waveform, 

369 unit=unit, 

370 times=times) 

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

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

373 try: 

374 spec, unit = \ 

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

376 fish.spec = spec 

377 except FileNotFoundError: 

378 fish.spec = None 

379 wave_fishes.append(fish) 

380 count += 1 

381 #if count > 300: # XXX REMOVE 

382 # break 

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

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

385 tstart = times.tstart[0] 

386 tend = times.tend[0] 

387 if start_time is not None: 

388 deltat = start_time - tstart 

389 tstart += deltat 

390 tend += deltat 

391 for fish in pulse_fishes + wave_fishes: 

392 fish.times.tstart += deltat 

393 fish.times.tend += deltat 

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

395 

396 

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

398 title, output_folder): 

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

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

401 n = 0 

402 for s in stds: 

403 n += s.shape[1] 

404 h = n*4.0 

405 t = 0.8 if title else 0.1 

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

407 sharex=True, sharey=True) 

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

409 hspace=0) 

410 i = 0 

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

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

413 ax = axs[i] 

414 ax.plot(time, std) 

415 if thresholds: 

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

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

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

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

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

421 #ax.plot(time, stdm) 

422 ax.set_yscale('log') 

423 #ax.set_ylim(bottom=0) 

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

425 i += 1 

426 if title: 

427 axs[0].set_title(title) 

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

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

430 rotation=30, rotation_mode='anchor') 

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

432 plt.show() 

433 

434 

435def merge_fish(pulse_fishes, wave_fishes, 

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

437 pulse_eods = [] 

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

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

440 continue 

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

442 continue 

443 if len(pulse_eods) > 0 and \ 

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

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

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

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

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

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

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

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

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

453 continue 

454 pulse_eods.append(pulse_fishes[i]) 

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

456 

457 wave_eods = [] 

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

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

460 continue 

461 if len(wave_eods) > 0 and \ 

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

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

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

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

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

467 wave_eods[-1].othereods = [] 

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

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

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

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

472 continue 

473 wave_eods.append(wave_fishes[i]) 

474 return pulse_eods, wave_eods 

475 

476 

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

478 channels, output_folder): 

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

480 '#C02717', '#873770', '#008797', '#007030', 

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

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

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

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

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

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

487 h = n*2.5 + 2.5 + 0.3 

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

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

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

491 hspace=0.2) 

492 pi = 0 

493 prev_xscale = 0.0 

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

495 # EOD waveform: 

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

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

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

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

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

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

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

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

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

505 if hasattr(fish, 'othereods'): 

506 for eod in fish.othereods: 

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

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

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

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

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

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

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

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

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

516 tmax = lim 

517 else: 

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

519 tmax = time[-1] 

520 xscale = 1.0 

521 if tmax < 1.0: 

522 xscale = 0.5 

523 elif tmax > 10.0: 

524 xscale = 5.0 

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

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

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

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

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

530 if xscale < 1.0: 

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

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

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

534 else: 

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

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

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

538 prev_xscale = xscale 

539 # time bar: 

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

541 min_eodf = 10000 

542 max_eodf = 0 

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

544 if time['EODf'] < min_eodf: 

545 min_eodf = time['EODf'] 

546 if time['EODf'] > max_eodf: 

547 max_eodf = time['EODf'] 

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

549 [time['EODf'], time['EODf']], 

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

551 if max_eodf > min_eodf + 10.0: 

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

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

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

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

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

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

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

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

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

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

562 else: 

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

564 rotation=30, rotation_mode='anchor') 

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

566 

567 

568def main(cargs=None): 

569 # config file name: 

570 cfgfile = __package__ + '.cfg' 

571 

572 # command line arguments: 

573 if cargs is None: 

574 cargs = sys.argv[1:] 

575 parser = argparse.ArgumentParser(add_help=False, 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

591 help='plot analyzed data') 

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

593 help='merge similar EODs before plotting') 

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