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 detect_eods 

26from .eodanalysis import save_eod_waveform, save_wave_fish, save_pulse_fish 

27from .eodanalysis import save_wave_spectrum, save_pulse_spectrum, save_pulse_phases 

28from .eodanalysis import load_eod_waveform, load_wave_fish, load_pulse_fish 

29from .eodanalysis import load_wave_spectrum, load_pulse_spectrum, load_pulse_phases 

30from .eodanalysis import wave_similarity, pulse_similarity 

31from .thunderfish import configuration, save_configuration 

32from .thunderfish import remove_eod_files 

33 

34 

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

36 default_thresh=0.002, thresh_fac=3.0, 

37 thresh_nbins=500): 

38 """Add parameters needed for by thunderlogger. 

39 

40 Parameters 

41 ---------- 

42 cfg: ConfigFile 

43 The configuration. 

44 detection_thresh: float or 'auto' 

45 Only data segments with standard deviation larger than this value 

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

47 from all the data segments of a recording channel. 

48 default_thresh: float 

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

50 no data are available. 

51 thresh_fac: float 

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

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

54 standard deviation. 

55 thresh_nbins: int 

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

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

58 deviation are estimated for automatically computing a threshold. 

59 """ 

60 cfg.add_section('Thunderlogger:') 

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

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

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

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

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

66 

67 

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

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

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

71 t0s = [] 

72 stds = None 

73 supra_thresh = None 

74 wave_fishes = None 

75 pulse_fishes = None 

76 if start_time is None: 

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

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

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

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

81 toffs = start_time 

82 t1 = start_time 

83 unit = None 

84 for file in files: 

85 try: 

86 with DataLoader(file) as sf: 

87 # analyze: 

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

89 unit = sf.unit 

90 if max_dist < 1.1/sf.rate: 

91 max_dist = 1.1/sf.rate 

92 window_size = cfg.value('windowSize') 

93 ndata = int(window_size * sf.rate) 

94 step = ndata//2 

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

96 if stds is None: 

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

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

99 if wave_fishes is None: 

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

101 if pulse_fishes is None: 

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

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

104 sys.stdout.write('.') 

105 sys.stdout.flush() 

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

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

108 t0s.append(t0) 

109 for channel in range(sf.channels): 

110 if thresholds: 

111 thresh = thresholds[channel] 

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

113 sd = np.std(fdata) 

114 stds[channel].append(sd) 

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

116 if stds_only: 

117 continue 

118 if sd > thresh: 

119 # clipping: 

120 min_clip = cfg.value('minClipAmplitude') 

121 if min_clip == 0.0: 

122 min_clip = cfg.value('minDataAmplitude') 

123 max_clip = cfg.value('maxClipAmplitude') 

124 if max_clip == 0.0: 

125 max_clip = cfg.value('maxDataAmplitude') 

126 name = file 

127 # detect EODs in the data: 

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

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

130 None, None, None, 

131 min_clip, max_clip, 

132 name, verbose, plot_level, 

133 cfg) 

134 first_fish = True 

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

136 spec_data, peak_data): 

137 fish = None 

138 fish_deltaf = 100000.0 

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

140 for wfish in wave_fishes[channel]: 

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

142 if deltaf < fish_deltaf: 

143 fish_deltaf = deltaf 

144 fish = wfish 

145 if fish_deltaf > max_deltaf: 

146 fish = None 

147 peaks = None 

148 else: 

149 fish_dist = 10000.0 

150 for pfish in pulse_fishes[channel]: 

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

152 props['peak-dist']) 

153 if ddist < fish_dist: 

154 fish_dist = ddist 

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

156 props['EODf']) 

157 fish = pfish 

158 if fish_dist > max_dist or \ 

159 fish_deltaf > max_deltaf: 

160 fish = None 

161 spec = None 

162 if fish is not None and \ 

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

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

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

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

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

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

169 else: 

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

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

172 fish.props = props 

173 fish.waveform = eod 

174 fish.spec = spec 

175 fish.peaks = peaks 

176 else: 

177 new_fish = SimpleNamespace(props=props, 

178 waveform=eod, 

179 spec=spec, 

180 peaks=peaks, 

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

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

183 pulse_fishes[channel].append(new_fish) 

184 else: 

185 wave_fishes[channel].append(new_fish) 

186 if first_fish: 

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

188 first_fish = False 

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

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

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

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

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

194 sys.stdout.flush() 

195 except EOFError as error: 

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

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

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

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

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

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

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

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

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

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

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

207 

208 

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

210 td = TableData() 

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

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

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

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

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

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

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

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

219 [name for t in times]) 

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

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

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

223 fp = output_basename + '-times' 

224 if idx is not None: 

225 fp += '-%d' % idx 

226 td.write(fp, **kwargs) 

227 

228 

229def load_times(file_path): 

230 data = TableData(file_path).data_frame() 

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

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

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

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

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

236 if 'channel' in data: 

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

238 return data 

239 

240 

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

242 td = TableData() 

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

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

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

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

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

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

249 fp = output_basename + '-stdevs' 

250 td.write(fp, **kwargs) 

251 

252 

253def load_power(file_path, start_time=None): 

254 base = os.path.basename(file_path) 

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

256 data = TableData(file_path) 

257 times = [] 

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

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

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

261 if start_time is not None: 

262 deltat = start_time - times[0] 

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

264 times[k] += deltat 

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

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

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

268 for c in range(channels): 

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

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

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

272 

273 

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

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

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

277 if pulse_fishes is not None: 

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

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

280 idx = 0 

281 # pulse fish: 

282 pulse_props = [] 

283 for fish in pulse_fishes[c]: 

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

285 **write_table_args(cfg)) 

286 if fish.peaks is not None: 

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

288 **write_table_args(cfg)) 

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

290 **write_table_args(cfg)) 

291 pulse_props.append(fish.props) 

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

293 idx += 1 

294 save_pulse_fish(pulse_props, unit, out_path, 

295 **write_table_args(cfg)) 

296 # wave fish: 

297 wave_props = [] 

298 if wave_fishes is not None: 

299 for fish in wave_fishes[c]: 

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

301 **write_table_args(cfg)) 

302 if fish.spec is not None: 

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

304 **write_table_args(cfg)) 

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

306 **write_table_args(cfg)) 

307 wave_props.append(fish.props) 

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

309 idx += 1 

310 save_wave_fish(wave_props, unit, out_path, 

311 **write_table_args(cfg)) 

312 # recording time window: 

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

314 **write_table_args(cfg)) 

315 # signal power: 

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

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

318 **write_table_args(cfg)) 

319 

320 

321def load_data(files, start_time=None): 

322 all_files = [] 

323 for file in files: 

324 if os.path.isdir(file): 

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

326 else: 

327 all_files.append(file) 

328 pulse_fishes = [] 

329 wave_fishes = [] 

330 channels = set() 

331 for file in all_files: 

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

333 pulse_props = load_pulse_fish(file) 

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

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

336 count = 0 

337 for props in pulse_props: 

338 idx = props['index'] 

339 waveform, unit = \ 

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

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

342 times['index'] = idx 

343 fish = SimpleNamespace(props=props, 

344 waveform=waveform, 

345 unit=unit, 

346 times=times) 

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

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

349 try: 

350 peaks, unit = \ 

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

352 fish.peaks = peaks 

353 except FileNotFoundError: 

354 fish.peaks = None 

355 pulse_fishes.append(fish) 

356 count += 1 

357 #if count > 300: # XXX REMOVE 

358 # break 

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

360 wave_props = load_wave_fish(file) 

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

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

363 count = 0 

364 for props in wave_props: 

365 idx = props['index'] 

366 waveform, unit = \ 

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

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

369 times['index'] = idx 

370 fish = SimpleNamespace(props=props, 

371 waveform=waveform, 

372 unit=unit, 

373 times=times) 

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

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

376 try: 

377 spec, unit = \ 

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

379 fish.spec = spec 

380 except FileNotFoundError: 

381 fish.spec = None 

382 wave_fishes.append(fish) 

383 count += 1 

384 #if count > 300: # XXX REMOVE 

385 # break 

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

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

388 tstart = times.tstart[0] 

389 tend = times.tend[0] 

390 if start_time is not None: 

391 deltat = start_time - tstart 

392 tstart += deltat 

393 tend += deltat 

394 for fish in pulse_fishes + wave_fishes: 

395 fish.times.tstart += deltat 

396 fish.times.tend += deltat 

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

398 

399 

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

401 title, output_folder): 

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

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

404 n = 0 

405 for s in stds: 

406 n += s.shape[1] 

407 h = n*4.0 

408 t = 0.8 if title else 0.1 

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

410 sharex=True, sharey=True) 

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

412 hspace=0) 

413 i = 0 

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

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

416 ax = axs[i] 

417 ax.plot(time, std) 

418 if thresholds: 

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

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

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

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

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

424 #ax.plot(time, stdm) 

425 ax.set_yscale('log') 

426 #ax.set_ylim(bottom=0) 

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

428 i += 1 

429 if title: 

430 axs[0].set_title(title) 

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

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

433 rotation=30, rotation_mode='anchor') 

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

435 plt.show() 

436 

437 

438def merge_fish(pulse_fishes, wave_fishes, 

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

440 pulse_eods = [] 

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

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

443 continue 

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

445 continue 

446 if len(pulse_eods) > 0 and \ 

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

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

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

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

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

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

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

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

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

456 continue 

457 pulse_eods.append(pulse_fishes[i]) 

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

459 

460 wave_eods = [] 

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

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

463 continue 

464 if len(wave_eods) > 0 and \ 

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

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

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

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

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

470 wave_eods[-1].othereods = [] 

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

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

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

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

475 continue 

476 wave_eods.append(wave_fishes[i]) 

477 return pulse_eods, wave_eods 

478 

479 

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

481 channels, output_folder): 

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

483 '#C02717', '#873770', '#008797', '#007030', 

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

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

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

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

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

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

490 h = n*2.5 + 2.5 + 0.3 

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

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

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

494 hspace=0.2) 

495 pi = 0 

496 prev_xscale = 0.0 

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

498 # EOD waveform: 

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

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

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

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

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

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

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

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

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

508 if hasattr(fish, 'othereods'): 

509 for eod in fish.othereods: 

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

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

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

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

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

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

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

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

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

519 tmax = lim 

520 else: 

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

522 tmax = time[-1] 

523 xscale = 1.0 

524 if tmax < 1.0: 

525 xscale = 0.5 

526 elif tmax > 10.0: 

527 xscale = 5.0 

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

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

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

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

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

533 if xscale < 1.0: 

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

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

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

537 else: 

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

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

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

541 prev_xscale = xscale 

542 # time bar: 

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

544 min_eodf = 10000 

545 max_eodf = 0 

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

547 if time['EODf'] < min_eodf: 

548 min_eodf = time['EODf'] 

549 if time['EODf'] > max_eodf: 

550 max_eodf = time['EODf'] 

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

552 [time['EODf'], time['EODf']], 

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

554 if max_eodf > min_eodf + 10.0: 

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

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

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

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

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

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

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

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

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

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

565 else: 

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

567 rotation=30, rotation_mode='anchor') 

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

569 

570 

571def main(cargs=None): 

572 # config file name: 

573 cfgfile = __package__ + '.cfg' 

574 

575 # command line arguments: 

576 if cargs is None: 

577 cargs = sys.argv[1:] 

578 parser = argparse.ArgumentParser(add_help=False, 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

594 help='plot analyzed data') 

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