j_hartling/paper_2025
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Syncing to home.

Browse Source
This commit is contained in:
j-hartling
2026-05-08 18:21:47 +02:00
parent 4b4a04ab2a
commit f14de13823
16 changed files with 578 additions and 395 deletions
Download Patch File Download Diff File Expand all files Collapse all files

452
python/fig_invariance_field.py
View File

@@ -5,11 +5,10 @@ from itertools import product
from thunderhopper.filetools import search_files
from thunderhopper.modeltools import load_data
from thunderhopper.filtertools import find_kern_specs
from misc_functions import get_saturation, reduce_kernel_set, exclude_zero_scale,\
divide_by_zero, x_dist, y_dist
from misc_functions import reduce_kernel_set, divide_by_zero, y_dist
from color_functions import load_colors
from plot_functions import hide_axis, reorder_by_sd, ylimits, super_xlabel, ylabel, title_subplot,\
plot_line, strip_zeros, time_bar, assign_colors,\
plot_line, xlabel, time_bar, assign_colors,\
letter_subplot, letter_subplots, hide_ticks
from IPython import embed
@@ -21,12 +20,12 @@ def plot_snippets(axes, time, snippets, ymin=None, ymax=None, **kwargs):
ymin=ymin, ymax=ymax, **kwargs))
return handles
def plot_curves(ax, scales, measures, **kwargs):
def plot_curves(ax, distances, measures, **kwargs):
if measures.ndim == 1:
handles = ax.plot(scales, measures, **kwargs)
handles = ax.plot(distances, measures, **kwargs)
return handles, measures
median_measure = np.nanmedian(measures, axis=1)
line_handle = ax.plot(scales, median_measure, **kwargs)[0]
line_handle = ax.plot(distances, median_measure, **kwargs)[0]
return line_handle, median_measure
def crop_noise_snippets(snippets, nin, nout, stages=['filt', 'env', 'log', 'inv', 'conv', 'feat']):
@@ -43,17 +42,17 @@ search_target = 'Pseudochorthippus_parallelus'
stages = ['filt', 'env', 'log', 'inv', 'conv', 'feat']
song_example = 'Pseudochorthippus_parallelus_micarray-short_JJ_20240815T160355-20240815T160755-1m10s690ms-1m13s614ms'
noise_example = 'merged_noise'
song_path = '../data/inv/field/song/'
noise_path = '../data/inv/field/noise/'
raw_path = search_files(search_target, incl='unnormed', dir=song_path + 'condensed/')[0]
base_path = search_files(search_target, incl='base', dir=song_path + 'condensed/')[0]
range_path = search_files(search_target, incl='range', dir=song_path + 'condensed/')[0]
song_snip_path = search_files(song_example, dir=song_path)[0]
noise_snip_path = search_files(noise_example, dir=noise_path)[0]
song_path = search_files(song_example, dir='../data/inv/field/song/')[0]
noise_path = search_files(noise_example, dir='../data/inv/field/noise/')[0]
save_path = '../figures/fig_invariance_field.pdf'
# ANALYSIS SETTINGS:
offset_distance = 10 # centimeter
thresh_rel = np.array([0, 0.5, 1, 1.5, 2, 2.5, 3])[4]
kern_subset_kwargs = dict(
combis=[['measure', 'snip'], ['conv', 'feat']],
keys=['thresh_abs'],
)
# SUBSET SETTINGS:
types = np.array([1, -1, 2, -2, 3, -3, 4, -4])
@@ -69,39 +68,53 @@ fig_kwargs = dict(
)
super_grid_kwargs = dict(
nrows=2,
ncols=1,
ncols=2,
wspace=0,
hspace=0,
left=0,
right=1,
bottom=0,
top=1,
height_ratios=[3, 2]
height_ratios=[1, 1]
)
subfig_specs = dict(
snip=(0, 0),
big=(1, 0),
snip=(0, slice(None)),
raw=(1, 0),
base=(1, 1),
)
snip_grid_kwargs = dict(
nrows=len(stages),
ncols=None,
wspace=0.1,
hspace=0.4,
left=0.11,
left=0.13,
right=0.98,
bottom=0.08,
bottom=0.05,
top=0.95
)
big_grid_kwargs = dict(
nrows=1,
ncols=3,
wspace=0.4,
hspace=0,
left=snip_grid_kwargs['left'],
right=snip_grid_kwargs['right'],
bottom=0.13,
top=0.98
raw_grid_kwargs = dict(
nrows=2,
ncols=1,
wspace=0,
hspace=0.15,
left=0.14,
right=0.9,
bottom=0.25,
top=0.95,
height_ratios=[0.8, 0.2]
)
base_grid_kwargs = dict(
nrows=3,
ncols=1,
wspace=0,
hspace=0.25,
left=raw_grid_kwargs['left'],
right=raw_grid_kwargs['right'],
bottom=raw_grid_kwargs['bottom'],
top=raw_grid_kwargs['top'],
)
inset_dist_bounds = [1.01, 0, 0.95, 1]
inset_ax_bounds = [raw_grid_kwargs['left'], 0.1, raw_grid_kwargs['right'] - raw_grid_kwargs['left'], 0.01]
# PLOT SETTINGS:
fs = dict(
@@ -112,9 +125,11 @@ fs = dict(
tit_tex=20,
bar=16,
)
colors = load_colors('../data/stage_colors.npz')
conv_colors = load_colors('../data/conv_colors_all.npz')
feat_colors = load_colors('../data/feat_colors_all.npz')
stage_colors = load_colors('../data/stage_colors.npz')
kern_colors = dict(
conv=load_colors('../data/conv_colors_subset.npz'),
feat=load_colors('../data/feat_colors_subset.npz')
)
lw = dict(
filt=0.25,
env=0.25,
@@ -122,11 +137,14 @@ lw = dict(
inv=0.25,
conv=0.25,
feat=1,
big=3,
plateau=1.5,
single=3,
swarm=1,
legend=5,
dist=1
)
xlabels = dict(
big='distance [cm]',
high='$1\\,/\\,d\\,\\sim\\,\\alpha$ [cm$^{-1}$]',
low='distance $d$ [cm]',
)
ylabels = dict(
filt='$x_{\\text{filt}}$\n$[\\text{a.u.}]$',
@@ -135,33 +153,41 @@ ylabels = dict(
inv='$x_{\\text{adapt}}$\n$[\\text{dB}]$',
conv='$c_i$\n$[\\text{dB}]$',
feat='$f_i$',
big=['measure', 'rel. measure', 'norm. measure']
raw=['$m$', '$\\mu_{f_i}$'],
base=['$m\\,/\\,m_{\\eta}$', '$\\sigma_{c_i}\\,/\\,\\sigma_{\\eta_i}$', '$\\mu_{f_i}\\,/\\,\\mu_{\\eta_i}$']
)
xlab_big_kwargs = dict(
xlab_high_kwargs = dict(
y=0.15,
fontsize=fs['lab_norm'],
ha='center',
va='bottom',
)
xlab_low_kwargs = dict(
y=0,
fontsize=fs['lab_norm'],
ha='center',
va='bottom',
)
ylab_snip_kwargs = dict(
x=0,
x=0.03,
fontsize=fs['lab_tex'],
rotation=0,
ha='left',
va='center'
ha='center',
va='center',
ma='center'
)
ylab_big_kwargs = dict(
x=-0.2,
x=0,
fontsize=fs['lab_norm'],
ha='center',
va='bottom',
va='top',
)
yloc = dict(
filt=0.03,
env=0.01,
filt=300,
env=100,
log=50,
inv=20,
conv=1,
conv=0.5,
feat=1,
)
title_kwargs = dict(
@@ -179,7 +205,7 @@ letter_snip_kwargs = dict(
fontsize=fs['letter'],
)
letter_big_kwargs = dict(
x=0,
xref=0,
y=1,
ha='left',
va='bottom',
@@ -206,49 +232,64 @@ noise_bar_time = 0.5
noise_bar_kwargs = song_bar_kwargs.copy()
noise_bar_kwargs['dur'] = noise_bar_time
noise_bar_kwargs['text_str'] = f'${int(1000 * noise_bar_time)}\\,\\text{{ms}}$'
plateau_settings = dict(
low=0.05,
high=0.95,
first=True,
last=True,
condense=None,
leg_labels = dict(
filt='$x_{\\text{filt}}$',
env='$x_{\\text{env}}$',
log='$x_{\\text{log}}$',
inv='$x_{\\text{adapt}}$',
conv='$c_i$',
feat='$f_i$'
)
plateau_line_kwargs = dict(
lw=lw['plateau'],
ls='--',
zorder=1,
leg_kwargs = dict(
ncols=3,
loc='upper left',
bbox_to_anchor=(0.025, 0.9, 0.95, 0.1),
frameon=False,
prop=dict(
size=20,
),
borderpad=0,
borderaxespad=0,
handlelength=1,
columnspacing=1,
handletextpad=0.5,
labelspacing=0.1
)
plateau_dot_kwargs = dict(
marker='o',
markersize=8,
markeredgewidth=1,
clip_on=False,
dist_line_kwargs = dict(
lw=lw['dist'],
)
dist_fill_kwargs = dict(
lw=lw['dist'],
)
# EXECUTION:
# Load raw (unnormed) invariance data:
data, config = load_data(raw_path, files='distances', keywords='mean')
dists = data['distances'] + offset_distance
# Load snippet data:
song_snip, _ = load_data(song_snip_path, keywords='snip')
t_song = np.arange(song_snip['snip_filt'].shape[0]) / config['rate']
noise_snip, _ = load_data(noise_snip_path, keywords='snip')
noise_snip = crop_noise_snippets(noise_snip, noise_snip['snip_filt'].shape[0], t_song.size)
t_noise = np.arange(noise_snip['snip_filt'].shape[0]) / config['rate']
# Load song invariance data:
song_data, config = load_data(song_path, files='distances', keywords=['measure', 'snip', 'thresh'])
t_song = np.arange(song_data['snip_filt'].shape[0]) / config['rate']
dists = song_data['distances'] + offset_distance
scales = 1 / dists
snip_dists = ['noise'] + [f'{int(d)}$\\,$cm' for d in dists]
# Optional kernel subset:
reduce_kernels = False
if any(var is not None for var in [kernels, types, sigmas]):
# Load noise invariance data:
noise_data, _ = load_data(noise_path, keywords=['measure', 'snip', 'thresh'])
noise_data = crop_noise_snippets(noise_data, noise_data['snip_filt'].shape[0], t_song.size)
t_noise = np.arange(noise_data['snip_filt'].shape[0]) / config['rate']
# Reduce kernels:
if reduce_kernels:
kern_inds = find_kern_specs(config['k_specs'], kernels, types, sigmas)
data = reduce_kernel_set(data, kern_inds, keyword='mean')
song_snip = reduce_kernel_set(song_snip, kern_inds, keyword='snip')
noise_snip = reduce_kernel_set(noise_snip, kern_inds, keyword='snip')
config['k_specs'] = config['k_specs'][kern_inds, :]
config['kernels'] = config['kernels'][:, kern_inds]
reduce_kernels = True
song_data = reduce_kernel_set(song_data, kern_inds, **kern_subset_kwargs)
noise_data = reduce_kernel_set(noise_data, kern_inds, **kern_subset_kwargs)
# Reduce thresholds:
thresh_ind = np.nonzero(song_data['thresh_rel'] == thresh_rel)[0][0]
song_data['measure_feat'] = song_data['measure_feat'][:, :, thresh_ind]
song_data['snip_feat'] = song_data['snip_feat'][:, :, :, thresh_ind]
noise_data['measure_feat'] = noise_data['measure_feat'][:, :, thresh_ind]
noise_data['snip_feat'] = noise_data['snip_feat'][:, :, :, thresh_ind]
# Adjust grid parameters:
snip_grid_kwargs['ncols'] = len(snip_dists)
@@ -278,136 +319,173 @@ time_bar(snip_axes[-1, -1], **song_bar_kwargs)
# time_bar(snip_axes[-1, 0], **noise_bar_kwargs)
letter_subplot(snip_subfig, 'a', ref=title, **letter_snip_kwargs)
# Prepare analysis axes:
big_subfig = fig.add_subfigure(super_grid[subfig_specs['big']])
big_grid = big_subfig.add_gridspec(**big_grid_kwargs)
big_axes = np.zeros((big_grid.ncols,), dtype=object)
for i in range(big_grid.ncols):
ax = big_subfig.add_subplot(big_grid[0, i])
ax.set_xlim(dists[0], 0)
# ax.set_xscale('symlog', linthresh=offset_distance, linscale=0.5)
ax.set_yscale('symlog', linthresh=0.01, linscale=0.1)
ylabel(ax, ylabels['big'][i], **ylab_big_kwargs)
# if i < (big_grid.ncols - 1):
# ax.set_ylim(scales[0], scales[-1])
# else:
# ax.set_ylim(0, 1)
big_axes[i] = ax
super_xlabel(xlabels['big'], big_subfig, big_axes[0], big_axes[-1], **xlab_big_kwargs)
letter_subplots(big_axes, 'bcd', **letter_big_kwargs)
# Prepare raw analysis axes:
raw_subfig = fig.add_subfigure(super_grid[subfig_specs['raw']])
raw_grid = raw_subfig.add_gridspec(**raw_grid_kwargs)
raw_axes = np.zeros((raw_grid.nrows,), dtype=object)
for i in range(raw_grid.nrows):
ax = raw_subfig.add_subplot(raw_grid[i, 0])
ax.set_xlim(scales[0], scales[-1])
ax.set_xscale('log')
ylabel(ax, ylabels['raw'][i], transform=raw_subfig.transSubfigure, **ylab_big_kwargs)
if i == 0:
ax.set_yscale('symlog', linthresh=0.00001, linscale=0.1)
hide_ticks(ax, 'bottom')
else:
transform = raw_subfig.transSubfigure + ax.transAxes.inverted()
inset_x1 = transform.transform((inset_dist_bounds[2], 0))[0]
inset_dist_bounds[2] = inset_x1 - inset_dist_bounds[0]
raw_inset = ax.inset_axes(inset_dist_bounds)
raw_inset.axis('off')
raw_axes[i] = ax
letter_subplots(raw_axes, 'bc', ref=raw_subfig, **letter_big_kwargs)
xlabel(raw_axes[-1], xlabels['high'], transform=raw_subfig.transSubfigure, **xlab_high_kwargs)
# Prepare base analysis axes:
base_subfig = fig.add_subfigure(super_grid[subfig_specs['base']])
base_grid = base_subfig.add_gridspec(**base_grid_kwargs)
base_axes = np.zeros((base_grid.nrows,), dtype=object)
base_insets = np.zeros((base_grid.nrows - 1,), dtype=object)
for i in range(base_grid.nrows):
ax = base_subfig.add_subplot(base_grid[i, 0])
ax.set_xlim(scales[0], scales[-1])
ax.set_xscale('log')
ax.set_yscale('log')
ylabel(ax, ylabels['base'][i], transform=base_subfig.transSubfigure, **ylab_big_kwargs)
if i < base_grid_kwargs['nrows'] - 1:
hide_ticks(ax, 'bottom')
if i > 0:
inset = ax.inset_axes(inset_dist_bounds)
inset.set_yscale('log')
inset.axis('off')
base_insets[i - 1] = inset
base_axes[i] = ax
letter_subplots(base_axes, 'def', ref=base_subfig, **letter_big_kwargs)
xlabel(base_axes[-1], xlabels['high'], transform=base_subfig.transSubfigure, **xlab_high_kwargs)
if True:
# Plot filtered snippets:
plot_snippets(snip_axes[0, 1:], t_song, song_snip['snip_filt'],
c=colors['filt'], lw=lw['filt'])
plot_line(snip_axes[0, 0], t_noise, noise_snip['snip_filt'][:, 0],
*snip_axes[0, 1].get_ylim(), c=colors['filt'], lw=lw['filt'])
plot_snippets(snip_axes[0, 1:], t_song, song_data['snip_filt'],
c=stage_colors['filt'], lw=lw['filt'])
plot_line(snip_axes[0, 0], t_noise, noise_data['snip_filt'][:, 0],
*snip_axes[0, 1].get_ylim(), c=stage_colors['filt'], lw=lw['filt'])
# Plot envelope snippets:
plot_snippets(snip_axes[1, 1:], t_song, song_snip['snip_env'],
ymin=0, c=colors['env'], lw=lw['env'])
plot_line(snip_axes[1, 0], t_noise, noise_snip['snip_env'][:, 0],
*snip_axes[1, 1].get_ylim(), c=colors['env'], lw=lw['env'])
plot_snippets(snip_axes[1, 1:], t_song, song_data['snip_env'],
ymin=0, c=stage_colors['env'], lw=lw['env'])
plot_line(snip_axes[1, 0], t_noise, noise_data['snip_env'][:, 0],
*snip_axes[1, 1].get_ylim(), c=stage_colors['env'], lw=lw['env'])
# Plot logarithmic snippets:
plot_snippets(snip_axes[2, 1:], t_song, song_snip['snip_log'],
c=colors['log'], lw=lw['log'])
plot_line(snip_axes[2, 0], t_noise, noise_snip['snip_log'][:, 0],
*snip_axes[2, 1].get_ylim(), c=colors['log'], lw=lw['log'])
plot_snippets(snip_axes[2, 1:], t_song, song_data['snip_log'],
c=stage_colors['log'], lw=lw['log'])
plot_line(snip_axes[2, 0], t_noise, noise_data['snip_log'][:, 0],
*snip_axes[2, 1].get_ylim(), c=stage_colors['log'], lw=lw['log'])
# Plot invariant snippets:
plot_snippets(snip_axes[3, 1:], t_song, song_snip['snip_inv'],
c=colors['inv'], lw=lw['inv'])
plot_line(snip_axes[3, 0], t_noise, noise_snip['snip_inv'][:, 0],
*snip_axes[3, 1].get_ylim(), c=colors['inv'], lw=lw['inv'])
plot_snippets(snip_axes[3, 1:], t_song, song_data['snip_inv'],
c=stage_colors['inv'], lw=lw['inv'])
plot_line(snip_axes[3, 0], t_noise, noise_data['snip_inv'][:, 0],
*snip_axes[3, 1].get_ylim(), c=stage_colors['inv'], lw=lw['inv'])
# Plot kernel response snippets:
all_handles = plot_snippets(snip_axes[4, 1:], t_song, song_snip['snip_conv'],
c=colors['conv'], lw=lw['conv'])
all_handles = plot_snippets(snip_axes[4, 1:], t_song, song_data['snip_conv'],
c=stage_colors['conv'], lw=lw['conv'])
for i, handles in enumerate(all_handles):
assign_colors(handles, config['k_specs'][:, 0], conv_colors)
reorder_by_sd(handles, song_snip['snip_conv'][..., i])
handles = plot_line(snip_axes[4, 0], t_noise, noise_snip['snip_conv'][:, 0],
*snip_axes[4, 1].get_ylim(), c=colors['conv'], lw=lw['conv'])
assign_colors(handles, config['k_specs'][:, 0], conv_colors)
reorder_by_sd(handles, noise_snip['snip_conv'][:, 0])
assign_colors(handles, config['k_specs'][:, 0], kern_colors['conv'])
reorder_by_sd(handles, song_data['snip_conv'][..., i])
handles = plot_line(snip_axes[4, 0], t_noise, noise_data['snip_conv'][:, 0],
*snip_axes[4, 1].get_ylim(), c=stage_colors['conv'], lw=lw['conv'])
assign_colors(handles, config['k_specs'][:, 0], kern_colors['conv'])
reorder_by_sd(handles, noise_data['snip_conv'][:, 0])
# Plot feature snippets:
all_handles = plot_snippets(snip_axes[5, 1:], t_song, song_snip['snip_feat'],
ymin=0, ymax=1, c=colors['feat'], lw=lw['feat'])
all_handles = plot_snippets(snip_axes[5, 1:], t_song, song_data['snip_feat'],
ymin=0, ymax=1, c=stage_colors['feat'], lw=lw['feat'])
for i, handles in enumerate(all_handles):
assign_colors(handles, config['k_specs'][:, 0], feat_colors)
reorder_by_sd(handles, song_snip['snip_feat'][..., i])
handles = plot_line(snip_axes[5, 0], t_noise, noise_snip['snip_feat'][:, 0],
ymin=0, ymax=1, c=colors['feat'], lw=lw['feat'])
assign_colors(handles, config['k_specs'][:, 0], feat_colors)
reorder_by_sd(handles, noise_snip['snip_feat'][:, 0])
del song_snip, noise_snip
assign_colors(handles, config['k_specs'][:, 0], kern_colors['feat'])
reorder_by_sd(handles, song_data['snip_feat'][..., i])
handles = plot_line(snip_axes[5, 0], t_noise, noise_data['snip_feat'][:, 0],
ymin=0, ymax=1, c=stage_colors['feat'], lw=lw['feat'])
assign_colors(handles, config['k_specs'][:, 0], kern_colors['feat'])
reorder_by_sd(handles, noise_data['snip_feat'][:, 0])
# Remember saturation points:
crit_inds, crit_dists = {}, {}
# Unnormed measures:
# Plot analysis results:
leg_handles = []
for stage in stages:
# Plot average intensity measure across recordings:
curve = plot_curves(big_axes[0], dists, data[f'mean_{stage}'],
c=colors[stage], lw=lw['big'],
fill_kwargs=dict(color=colors[stage], alpha=0.25))
# # Indicate saturation point:
# if stage in ['log', 'inv', 'conv', 'feat']:
# ind = get_saturation(curve, **plateau_settings)[1]
# dist = dists[ind]
# big_axes[0].plot(dist, 0, c='w', alpha=1, zorder=5.5, **plateau_dot_kwargs,
# transform=big_axes[0].get_xaxis_transform())
# big_axes[0].plot(dist, 0, mfc=colors[stage], mec='k', alpha=0.75, zorder=6, **plateau_dot_kwargs,
# transform=big_axes[0].get_xaxis_transform())
# big_axes[0].vlines(dist, big_axes[0].get_ylim()[0], curve[ind],
# color=colors[stage], **plateau_line_kwargs)
# # Log saturation point:
# crit_inds[stage] = ind
# crit_dists[stage] = dist
del data
mkey = f'measure_{stage}'
measure = song_data[mkey]
color = stage_colors[stage]
# Noise baseline-related measures:
data, _ = load_data(base_path, files='scales', keywords='mean')
if reduce_kernels:
data = reduce_kernel_set(data, kern_inds, keyword='mean')
for stage in stages:
# Plot average intensity measure across recordings:
curve = plot_curves(big_axes[1], dists, data[f'mean_{stage}'],
c=colors[stage], lw=lw['big'],
fill_kwargs=dict(color=colors[stage], alpha=0.25))
# Indicate saturation point:
# if stage in ['log', 'inv', 'conv', 'feat']:
# ind, dist = crit_inds[stage], crit_dists[stage]
# big_axes[1].plot(dist, 0, c='w', alpha=1, zorder=5.5, **plateau_dot_kwargs,
# transform=big_axes[1].get_xaxis_transform())
# big_axes[1].plot(dist, 0, mfc=colors[stage], mec='k', alpha=0.75, zorder=6, **plateau_dot_kwargs,
# transform=big_axes[1].get_xaxis_transform())
# big_axes[1].vlines(dist, big_axes[1].get_ylim()[0], curve[ind],
# color=colors[stage], **plateau_line_kwargs)
del data
## UNNORMALIZED MEASURE:
# Min-max normalized measures:
data, _ = load_data(range_path, files='scales', keywords='mean')
if reduce_kernels:
data = reduce_kernel_set(data, kern_inds, keyword='mean')
for stage in stages:
# Plot average intensity measure across recordings:
curve = plot_curves(big_axes[2], dists, data[f'mean_{stage}'],
c=colors[stage], lw=lw['big'],
fill_kwargs=dict(color=colors[stage], alpha=0.25))
# Plot single raw intensity curve (median where necessary):
handles, curve = plot_curves(raw_axes[0], scales, measure, c=color, lw=lw['single'])
# # Indicate saturation point:
# if stage in ['log', 'inv', 'conv', 'feat']:
# ind, dist = crit_inds[stage], crit_dists[stage]
# big_axes[2].plot(dist, 0, c='w', alpha=1, zorder=5.5, **plateau_dot_kwargs,
# transform=big_axes[2].get_xaxis_transform())
# big_axes[2].plot(dist, 0, mfc=colors[stage], mec='k', alpha=0.75, zorder=6, **plateau_dot_kwargs,
# transform=big_axes[2].get_xaxis_transform())
# big_axes[2].vlines(dist, big_axes[2].get_ylim()[0], curve[ind],
# color=colors[stage], **plateau_line_kwargs)
del data
# Add stage-specific proxy legend artist:
leg_handles.append(raw_axes[0].plot([], [], c=color, label=leg_labels[stage])[0])
# Plot curve swarm:
if stage == 'feat':
# Sync y-limits:
ylimits(measure, raw_axes[1], minval=0, pad=0.05)
raw_inset.set_ylim(raw_axes[1].get_ylim())
# Plot swarm:
handles = raw_axes[1].plot(scales, measure, lw=lw['swarm'])
assign_colors(handles, config['k_specs'][:, 0], kern_colors[stage])
reorder_by_sd(handles, measure)
# Plot distribution of saturation levels:
line_kwargs = dist_line_kwargs | dict(c=color)
fill_kwargs = dist_fill_kwargs | dict(color=color)
y_dist(raw_inset, measure[-1], nbins=75, log=False,
line_kwargs=line_kwargs, fill_kwargs=fill_kwargs)
## NORMALIZED MEASURE:
# Relate to noise baseline:
measure = divide_by_zero(song_data[mkey], noise_data[mkey].mean(axis=0))
# Plot single baseline-normalized intensity curve (median where necessary):
handles, curve = plot_curves(base_axes[0], scales, measure, c=color, lw=lw['single'])
# Plot curve swarm:
if stage in ['conv', 'feat']:
i0, i1 = (1, 0) if stage == 'conv' else (2, 1)
# Sync y-limits:
ylimits(measure, base_axes[i0], minval=0.9, pad=0.05)
base_insets[i1].set_ylim(base_axes[i0].get_ylim())
# Plot swarm:
handles = base_axes[i0].plot(scales, measure, lw=lw['swarm'])
assign_colors(handles, config['k_specs'][:, 0], kern_colors[stage])
reorder_by_sd(handles, measure)
# Plot distribution of saturation levels:
line_kwargs = dist_line_kwargs | dict(c=color)
fill_kwargs = dist_fill_kwargs | dict(color=color)
y_dist(base_insets[i1], measure[-1], nbins=100, log=True,
line_kwargs=line_kwargs, fill_kwargs=fill_kwargs)
# Posthoc adjustments:
raw_axes[0].set_ylim(top=100)
base_axes[0].set_ylim(1, 100)
base_axes[1].set_ylim(bottom=1)
base_insets[0].set_ylim(bottom=1)
base_axes[2].set_ylim(bottom=1)
base_insets[1].set_ylim(bottom=1)
# Add secondary x-axes:
for subfig in [raw_subfig, base_subfig]:
dual_ax = subfig.add_subplot(inset_ax_bounds)
dual_ax.set_xlim(scales[0], scales[-1])
dual_ax.set_xscale('log')
dual_ax.tick_params(axis='x', which='minor', bottom=False)
dual_ax.tick_params(axis='x', which='major', labelrotation=45)
dual_ax.set_xticks(scales, dists)
hide_axis(dual_ax, 'left')
xlabel(dual_ax, xlabels['low'], transform=subfig.transSubfigure, **xlab_low_kwargs)
# Add legend to first analysis axis:
legend = raw_axes[0].legend(handles=leg_handles, **leg_kwargs)
[handle.set_lw(lw['legend']) for handle in legend.get_lines()]
# Save graph:
if save_path is not None:

1
python/fig_invariance_full.py
View File

@@ -56,7 +56,6 @@ save_path = '../figures/fig_invariance_full.pdf'
# ANALYSIS SETTINGS:
exclude_zero = True
thresh_rel = np.array([0, 0.5, 1, 1.5, 2, 2.5, 3])[4]
percentiles = np.array([0, 100])
scale_subset_kwargs = dict(
combis=[['measure'], stages],
)

4
python/save_field_data.py
View File

@@ -8,11 +8,11 @@ from IPython import embed
# General:
search_target = '*'
mode = ['song', 'noise'][1]
mode = ['song', 'noise'][0]
input_folder = f'../data/field/raw/{mode}/'
output_folder = f'../data/field/processed/{mode}/'
stages = ['raw', 'norm']
if True:
if False:
# Overwrites edited:
stages.append('songs')

112
python/save_inv_data_field.py
View File

@@ -1,7 +1,9 @@
import numpy as np
import matplotlib.pyplot as plt
from thunderhopper.modeltools import load_data, save_data
from thunderhopper.filetools import search_files, crop_paths
from thunderhopper.filtertools import find_kern_specs
from thunderhopper.filters import sosfilter
from thunderhopper.model import process_signal
from IPython import embed
@@ -13,31 +15,24 @@ example_file = dict(
)[mode]
search_path = f'../data/field/processed/{mode}/'
data_paths = search_files('*', ext='npz', dir=search_path)
ref_path = '../data/inv/field/ref_measures.npz'
thresh_path = '../data/inv/field/thresholds.npz'
stages = ['raw', 'filt', 'env', 'log', 'inv', 'conv', 'feat']
pre_stages = stages[:-1]
save_path = f'../data/inv/field/{mode}/'
# ANALYSIS SETTINGS:
distances = np.load('../data/field/recording_distances.npy')[::-1]
thresh_rel = 0.5
thresh_rel = np.array([0, 0.5, 1, 1.5, 2, 2.5, 3])
init_scale = 10000
# SUBSET SETTINGS:
kernels = np.array([
[1, 0.002],
[-1, 0.002],
[2, 0.004],
[-2, 0.004],
[3, 0.032],
[-3, 0.032]
])
kernels = None
types = None#np.array([-1])
sigmas = None#np.array([0.001, 0.002, 0.004, 0.008, 0.016, 0.032])
types = None
sigmas = None
# PREPARATION:
if thresh_rel is not None:
# Get threshold values from pure-noise response SD:
thresh_abs = np.load(ref_path)['conv'] * thresh_rel
thresh_data = dict(np.load(thresh_path))
thresh_abs = thresh_rel[:, None] * thresh_data['sds'][None, :]
# EXECUTION:
for data_path, name in zip(data_paths, crop_paths(data_paths)):
@@ -48,9 +43,8 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)):
data, config = load_data(data_path, files='raw')
song, rate = data['raw'], config['rate']
if thresh_rel is not None:
# Set kernel-specific thresholds:
config['feat_thresh'] = thresh_abs
# Sort max to min distance:
song = song[:, ::-1] * init_scale
# Reduce to kernel subset:
if any(var is not None for var in [kernels, types, sigmas]):
@@ -58,7 +52,7 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)):
config['kernels'] = config['kernels'][:, kern_inds]
config['k_specs'] = config['k_specs'][kern_inds, :]
config['k_props'] = [config['k_props'][i] for i in kern_inds]
config['feat_thresh'] = config['feat_thresh'][kern_inds]
thresh_abs = thresh_abs[:, kern_inds]
# Get song segment to be analyzed:
time = np.arange(song.shape[0]) / rate
@@ -66,37 +60,81 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)):
segment = (time >= start) & (time <= end)
# Prepare storage:
measures = {}
shape = (distances.size, config['k_specs'].shape[0], thresh_rel.size)
measures = dict(measure_feat=np.zeros(shape, dtype=float))
if save_detailed:
snippets = {}
shape = (song.shape[0], config['k_specs'].shape[0], distances.size, thresh_rel.size)
snippets = dict(snip_feat=np.zeros(shape, dtype=float))
# Process snippet:
signals, rates = process_signal(config, returns=stages, signal=song, rate=rate)
for stage in stages:
# Sort largest to smallest distance:
signals[stage] = signals[stage][..., ::-1]
# Process snippet (excluding features):
signals, rates = process_signal(config, returns=pre_stages, signal=song, rate=rate)
# Store results:
for stage in stages:
# Log intensity measures:
# Store non-feature results:
for stage in pre_stages:
mkey = f'measure_{stage}'
if stage == 'feat':
measures[mkey] = signals[stage][segment, ...].mean(axis=0)
else:
measures[mkey] = signals[stage][segment, ...].std(axis=0)
if measures[mkey].ndim == 2:
# Make shape (distances, kernels):
# Log intensity measures:
measures[mkey] = signals[stage][segment, ...].std(axis=0)
if stage == 'conv':
# Make shape (distances, kernels) for consistency:
measures[mkey] = np.moveaxis(measures[mkey], 1, 0)
# Log optional snippet data:
if save_detailed:
# Log optional snippet data:
snippets[f'snip_{stage}'] = signals[stage]
conv = signals['conv']
# Execute piecewise per threshold:
for i, thresholds in enumerate(thresh_abs):
# Execute piecewise per distance:
for j in range(conv.shape[-1]):
feat = sosfilter((conv[:, :, j] > thresholds).astype(float),
rate, config['feat_fcut'], 'lp',
padtype='fixed', padlen=config['padlen'])
# Log intensity measure:
measure = feat[segment, ...].mean(axis=0)
measures['measure_feat'][j, :, i] = measure
if save_detailed:
# Log optional snippet data:
snippets['snip_feat'][:, :, j, i] = feat
# # Log intensity measure, ensuring shape (distances, kernels, thresholds):
# measures['measure_feat'][:, :, i] = np.moveaxis(feat[segment, ...].mean(axis=0), 1, 0)
# if save_detailed:
# # Log optional snippet data:
# snippets['snip_feat'][:, :, :, i] = feat
# thresholds = thresholds[None, :, None]
# embed()
# # Finalize processing:
# feat = sosfilter((signals['conv'] > thresholds).astype(float),
# rate, config['feat_fcut'], 'lp',
# padtype='fixed', padlen=config['padlen'])
# if i == thresholds.shape[0] - 1:
# fig, axes = plt.subplots(1, 8, sharex=True, sharey=True, figsize=(16, 9))
# for j, ax in enumerate(axes):
# ax.plot(time, feat[..., j])
# plt.show()
# embed()
# # Log intensity measure, ensuring shape (distances, kernels, thresholds):
# measures['measure_feat'][:, :, i] = np.moveaxis(feat[segment, ...].mean(axis=0), 1, 0)
# if save_detailed:
# # Log optional snippet data:
# snippets['snip_feat'][:, :, :, i] = feat
# Save analysis results:
if save_path is not None:
data = dict(
distances=distances,
thresh_rel=thresh_rel,
thresh_abs=thresh_abs,
)
data.update(measures)
if save_detailed:

4
python/save_inv_data_short.py
View File

@@ -12,8 +12,8 @@ from IPython import embed
target_species = [
# 'Chorthippus_biguttulus',
'Chorthippus_mollis',
'Chrysochraon_dispar',
'Euchorthippus_declivus',
# 'Chrysochraon_dispar',
# 'Euchorthippus_declivus',
'Gomphocerippus_rufus',
# 'Omocestus_rufipes',
# 'Pseudochorthippus_parallelus',

8
python/save_thresholds.py
View File

@@ -1,4 +1,5 @@
import numpy as np
import matplotlib.pyplot as plt
from thunderhopper.filters import sosfilter
from thunderhopper.model import convolve_kernels, process_signal
from thunderhopper.modeltools import load_data
@@ -54,7 +55,12 @@ elif mode == 'short':
conv = convolve_kernels(inv, config['kernels'], config['k_specs'])
elif mode == 'field':
starter = starter[:, channels].ravel(order='F')
conv = process_signal(config, 'conv', signal=starter, rate=config['rate'])[0]['conv']
conv = process_signal(config, 'conv', signal=starter, rate=config['rate'])[0]['conv']
# fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
# ax1.plot(starter)
# ax2.plot(conv)
# plt.show()
# embed()
# Get baseline kernel response SDs:
sds = conv[segment, :].std(axis=0)

28
python/temp_BS.py Normal file
View File

@@ -0,0 +1,28 @@
import numpy as np
import matplotlib.pyplot as plt
from thunderhopper.modeltools import load_data
from thunderhopper.filetools import search_files
from thunderhopper.model import process_signal
paths = search_files('Pseudochorthippus_parallelus_micarray-short_JJ_20240815T160355-20240815T160755-1m10s690ms-1m13s614ms', dir='../data/field/processed/song/')
thresh_rel = np.array([0, 0.5, 1, 1.5, 2, 2.5, 3])[-1]
thresh_abs = np.load('../data/inv/field/thresholds.npz')['sds'] * thresh_rel
for path in paths:
print(f'Processing {path}')
data, config = load_data(path, files='raw')
config['feat_thresh'] = thresh_abs
song, rate = data['raw'], config['rate']
time = np.arange(song.shape[0]) / rate
start, end = data['songs_0'].ravel()
segment = (time >= start) & (time <= end)
signals, rates = process_signal(config, 'feat', signal=song, rate=rate)
feat = signals['feat']
fig, axes = plt.subplots(1, 8, sharex=True, sharey=True, figsize=(16, 9))
for i, ax in enumerate(axes):
ax.plot(feat[..., i])
plt.show()