diff --git a/figures/fig_invariance_log_hp.pdf b/figures/fig_invariance_log_hp.pdf index 6779cd9..07dc478 100644 Binary files a/figures/fig_invariance_log_hp.pdf and b/figures/fig_invariance_log_hp.pdf differ diff --git a/figures/fig_invariance_thresh_lp.pdf b/figures/fig_invariance_thresh_lp.pdf index 651d1d1..d8a4cda 100644 Binary files a/figures/fig_invariance_thresh_lp.pdf and b/figures/fig_invariance_thresh_lp.pdf differ diff --git a/python/fig_invariance_full.py b/python/fig_invariance_full.py new file mode 100644 index 0000000..474da97 --- /dev/null +++ b/python/fig_invariance_full.py @@ -0,0 +1,272 @@ +import plotstyle_plt +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.transforms import BboxTransformTo +from itertools import product +from thunderhopper.filetools import search_files +from thunderhopper.modeltools import load_data +from color_functions import load_colors +from plot_functions import hide_axis, ylimits, xlabel, ylabel, plot_line, plot_barcode, strip_zeros +from IPython import embed + +def time_bar(ax, dur, y0=0.9, y1=0.95, xshift=0.5, parent=None, transform=None, **kwargs): + t0, t1 = ax.get_xlim() + offset = (t1 - t0 - dur) * xshift + x0 = t0 + offset + x1 = x0 + dur + if parent is None: + parent = ax + if transform is None: + transform = BboxTransformTo(parent.bbox) + if transform is not ax.transData: + trans = ax.transData + transform.inverted() + x0 = trans.transform((x0, 0))[0] + x1 = trans.transform((x1, 0))[0] + parent.add_artist(plt.Rectangle((x0, y0), x1 - x0, y1 - y0, + transform=transform, **kwargs)) + return None + +def add_snip_axes(fig, grid_kwargs): + grid = fig.add_gridspec(**grid_kwargs) + axes = np.zeros((grid.nrows, grid.ncols), dtype=object) + for i, j in product(range(grid.nrows), range(grid.ncols)): + axes[i, j] = fig.add_subplot(grid[i, j]) + [hide_axis(ax, 'left') for ax in axes.flatten()] + [hide_axis(ax, 'bottom') for ax in axes.flatten()] + return axes + +def plot_snippets(axes, time, snippets, ymin=None, ymax=None, **kwargs): + ymin, ymax = ylimits(snippets, minval=ymin, maxval=ymax, pad=0.05) + for ax, snippet in zip(axes, snippets.T): + plot_line(ax, time, snippet, ymin=ymin, ymax=ymax, **kwargs) + return None + +def plot_bi_snippets(axes, time, binary, **kwargs): + for ax, binary in zip(axes, binary.T): + plot_barcode(ax, time, binary[:, None], **kwargs) + return None + + +# GENERAL SETTINGS: +target = 'Omocestus_rufipes' +data_paths = search_files(target, excl='noise', dir='../data/inv/full/') +stages = ['filt', 'env', 'log', 'inv', 'conv', 'bi', 'feat'] +load_kwargs = dict( + files=stages, + keywords=['scales', 'measure', 'spread'] +) +save_path = None#'../figures/fig_invariance_full.pdf' + +# GRAPH SETTINGS: +fig_kwargs = dict( + figsize=(32/2.54, 16/2.54), +) +super_grid_kwargs = dict( + nrows=2, + ncols=2, + wspace=0, + hspace=0, + left=0, + right=1, + bottom=0, + top=1 +) +subfig_specs = dict( + pure=(0, 0), + noise=(1, 0), + analysis=(slice(None), 1) +) +pure_grid_kwargs = dict( + nrows=len(stages), + ncols=None, + wspace=0.05, + hspace=0.1, + left=0.07, + right=0.95, + bottom=0.15, + top=0.9 +) +noise_grid_kwargs = dict( + nrows=len(stages), + ncols=None, + wspace=0.05, + hspace=0.1, + left=0.07, + right=0.95, + bottom=0.15, + top=0.9 +) +analysis_grid_kwargs = dict( + nrows=1, + ncols=1, + wspace=0, + hspace=0, + left=0.15, + right=0.96, + bottom=0.1, + top=0.95 +) +snip_specs = dict( + conv=(0, slice(None)), + bi=(1, slice(None)), + feat=(2, slice(None)) +) + +# PLOT SETTINGS: +colors = load_colors('../data/stage_colors.npz') +lw_snippets = dict( + conv=0.5, + feat=2 +) +lw_analysis = 3 +xlabels = dict( + analysis='scale $\\alpha$', +) +xlab_analysis_kwargs = dict( + y=0.01, + fontsize=16, + ha='center', + va='bottom', +) +ylabels = dict( + conv='$c_i$', + bi='$b_i$', + feat='$f_i$', + analysis='ratio $\\text{SD}_{\\alpha}\\,/\\,\\text{SD}_{\\min[\\alpha]}$', + # analysis='ratio $\\sigma_{\\alpha}\\,/\\,\\sigma_{\\min[\\alpha]}$', +) +ylab_snip_kwargs = dict( + x=0.01, + fontsize=20, + rotation=0, + ha='left', + va='center', +) +ylab_analysis_kwargs = dict( + x=0.02, + fontsize=16, + ha='center', + va='top', +) +xloc = dict( + analysis=10, +) +letter_snip_kwargs = dict( + x=0.02, + y=1, + ha='left', + va='top', + fontsize=22, + fontweight='bold' +) +letter_analysis_kwargs = dict( + x=0, + y=1, + ha='left', + va='top', + fontsize=22, + fontweight='bold' +) +bar_time = 5 +bar_kwargs = dict( + y0=0.7, + y1=0.8, + color='k', + lw=0, +) +spread_kwargs = dict( + alpha=0.3, + lw=0, + zorder=0 +) +kernel_ind = 0 + +# EXECUTION: +for data_path in data_paths: + print(f'Processing {data_path}') + + # Load invariance data: + data, config = load_data(data_path, **load_kwargs) + t_full = np.arange(data['conv'].shape[0]) / config['env_rate'] + + # Reduce snippet data to kernel subset: + data['conv'] = data['conv'][:, kernel_ind] + data['bi'] = data['bi'][:, kernel_ind] + data['feat'] = data['feat'][:, kernel_ind] + + # Prepare overall graph: + fig = plt.figure(**fig_kwargs) + super_grid = fig.add_gridspec(**super_grid_kwargs) + + # Prepare pure-song snippet axes: + pure_subfig = fig.add_subfigure(super_grid[subfig_specs['pure']]) + pure_grid_kwargs['nrows' if pure_grid_kwargs['nrows'] is None else 'ncols'] = data['example_scales'].size + pure_axes = add_snip_axes(pure_subfig, pure_grid_kwargs) + for ax, stage in zip(pure_axes[:, 0], stages): + ylabel(ax, ylabels[stage], **ylab_snip_kwargs, + transform=pure_subfig.transSubfigure) + for ax, scale in zip(pure_axes[snip_specs['conv']], data['example_scales']): + ax.set_title(f'$\\alpha={strip_zeros(scale)}$') + pure_subfig.text(s='a', **letter_snip_kwargs) + + # Prepare analysis axis: + analysis_subfig = fig.add_subfigure(super_grid[subfig_specs['analysis']]) + analysis_grid = analysis_subfig.add_gridspec(**analysis_grid_kwargs) + analysis_ax = analysis_subfig.add_subplot(analysis_grid[0, 0]) + analysis_ax.set_xlim(data['scales'].min(), data['scales'].max()) + analysis_ax.xaxis.set_major_locator(plt.MultipleLocator(xloc['analysis'])) + xlabel(analysis_ax, xlabels['analysis'], **xlab_analysis_kwargs, + transform=analysis_subfig.transSubfigure) + # analysis_ax.set_yscale('log') + ylabel(analysis_ax, ylabels['analysis'], **ylab_analysis_kwargs, + transform=analysis_subfig.transSubfigure) + analysis_subfig.text(s='c', **letter_analysis_kwargs) + + # Plot pure-song kernel response snippets: + plot_snippets(pure_axes[snip_specs['conv']], t_full, data['conv'], + c=colors['conv'], lw=lw_snippets['conv']) + + # Plot pure-song binary snippets: + plot_bi_snippets(pure_axes[snip_specs['bi']], t_full, data['bi'], + color=colors['bi'], lw=0) + + # Plot pure-song feature snippets: + plot_snippets(pure_axes[snip_specs['feat']], t_full, data['feat'], + ymin=0, ymax=1, c=colors['feat'], lw=lw_snippets['feat']) + + # Indicate time scale: + time_bar(pure_axes[snip_specs['conv']][0], bar_time, **bar_kwargs) + + # # Plot noise-song kernel response snippets: + # plot_snippets(noise_axes[snip_specs['conv']], t_full, noise_data['conv'], + # c=colors['conv'], lw=lw_snippets['conv']) + + # # Plot noise-song binary snippets: + # plot_bi_snippets(noise_axes[snip_specs['bi']], t_full, noise_data['bi'], + # color=colors['bi'], lw=0) + + # # Plot noise-song feature snippets: + # plot_snippets(noise_axes[snip_specs['feat']], t_full, noise_data['feat'], + # ymin=0, ymax=1, c=colors['feat'], lw=lw_snippets['feat']) + + # # Indicate time scale: + # time_bar(noise_axes[snip_specs['conv']][0], bar_time, **bar_kwargs) + + # Plot noise-song SD ratios (limited): + analysis_ax.plot(data['scales'], data['measure_conv'], + c=colors['conv'], lw=lw_analysis) + lower, upper = data['spread_conv'] + analysis_ax.fill_between(data['scales'], lower, upper, + color=colors['conv'], **spread_kwargs) + analysis_ax.plot(data['scales'], data['measure_feat'], + c=colors['feat'], lw=lw_analysis) + lower, upper = data['spread_feat'] + analysis_ax.fill_between(data['scales'], lower, upper, + color=colors['feat'], **spread_kwargs) + + if save_path is not None: + fig.savefig(save_path) + plt.show() + +print('Done.') +embed() diff --git a/python/fig_invariance_log-hp.py b/python/fig_invariance_log-hp.py index 7d45678..7bc7906 100644 --- a/python/fig_invariance_log-hp.py +++ b/python/fig_invariance_log-hp.py @@ -1,31 +1,14 @@ import plotstyle_plt import numpy as np import matplotlib.pyplot as plt -from matplotlib.transforms import BboxTransformTo from itertools import product from thunderhopper.filetools import search_files from thunderhopper.modeltools import load_data from color_functions import load_colors -from plot_functions import hide_axis, ylimits, xlabel, ylabel, plot_line, strip_zeros +from plot_functions import hide_axis, ylimits, xlabel, ylabel,\ + plot_line, strip_zeros, time_bar from IPython import embed -def time_bar(ax, dur, y0=0.9, y1=0.95, xshift=0.5, parent=None, transform=None, **kwargs): - t0, t1 = ax.get_xlim() - offset = (t1 - t0 - dur) * xshift - x0 = t0 + offset - x1 = x0 + dur - if parent is None: - parent = ax - if transform is None: - transform = BboxTransformTo(parent.bbox) - if transform is not ax.transData: - trans = ax.transData + transform.inverted() - x0 = trans.transform((x0, 0))[0] - x1 = trans.transform((x1, 0))[0] - parent.add_artist(plt.Rectangle((x0, y0), x1 - x0, y1 - y0, - transform=transform, **kwargs)) - return None - def add_snip_axes(fig, grid_kwargs): grid = fig.add_gridspec(**grid_kwargs) axes = np.zeros((grid.nrows, grid.ncols), dtype=object) @@ -46,8 +29,10 @@ def plot_snippets(axes, time, snippets, ymin=None, ymax=None, **kwargs): target = 'Omocestus_rufipes' data_paths = search_files(target, excl='noise', dir='../data/inv/log_hp/') stages = ['env', 'log', 'inv'] -files = stages + ['scales', 'example_scales', 'limit', - 'measure_env', 'measure_log', 'measure_inv'] +load_kwargs = dict( + files=stages, + keywords=['scales', 'measure'] +) save_path = '../figures/fig_invariance_log_hp.pdf' # GRAPH SETTINGS: @@ -177,8 +162,8 @@ for data_path in data_paths: print(f'Processing {data_path}') # Load invariance data: - pure_data, config = load_data(data_path, files) - noise_data, _ = load_data(data_path.replace('.npz', '_noise.npz'), files) + pure_data, config = load_data(data_path, **load_kwargs) + noise_data, _ = load_data(data_path.replace('.npz', '_noise.npz'), **load_kwargs) t_full = np.arange(pure_data['env'].shape[0]) / config['env_rate'] # Prepare overall graph: diff --git a/python/fig_invariance_thresh-lp.py b/python/fig_invariance_thresh-lp.py index 1ca102e..b1857e2 100644 --- a/python/fig_invariance_thresh-lp.py +++ b/python/fig_invariance_thresh-lp.py @@ -1,33 +1,14 @@ -from pyparsing import alphanums - import plotstyle_plt import numpy as np import matplotlib.pyplot as plt -from matplotlib.transforms import BboxTransformTo from itertools import product from thunderhopper.filetools import search_files from thunderhopper.modeltools import load_data from color_functions import load_colors -from plot_functions import hide_axis, ylimits, xlabel, ylabel, plot_line, plot_barcode, strip_zeros +from plot_functions import hide_axis, ylimits, xlabel, ylabel,\ + plot_line, plot_barcode, strip_zeros, time_bar from IPython import embed -def time_bar(ax, dur, y0=0.9, y1=0.95, xshift=0.5, parent=None, transform=None, **kwargs): - t0, t1 = ax.get_xlim() - offset = (t1 - t0 - dur) * xshift - x0 = t0 + offset - x1 = x0 + dur - if parent is None: - parent = ax - if transform is None: - transform = BboxTransformTo(parent.bbox) - if transform is not ax.transData: - trans = ax.transData + transform.inverted() - x0 = trans.transform((x0, 0))[0] - x1 = trans.transform((x1, 0))[0] - parent.add_artist(plt.Rectangle((x0, y0), x1 - x0, y1 - y0, - transform=transform, **kwargs)) - return None - def add_snip_axes(fig, grid_kwargs): grid = fig.add_gridspec(**grid_kwargs) axes = np.zeros((grid.nrows, grid.ncols), dtype=object) @@ -53,8 +34,10 @@ def plot_bi_snippets(axes, time, binary, **kwargs): target = 'Omocestus_rufipes' data_paths = search_files(target, excl='noise', dir='../data/inv/thresh_lp/') stages = ['conv', 'bi', 'feat'] -files = stages + ['scales', 'example_scales', 'measure_conv', 'spread_conv', - 'measure_feat', 'spread_feat'] +load_kwargs = dict( + files=stages, + keywords=['scales', 'measure', 'spread'] +) save_path = '../figures/fig_invariance_thresh_lp.pdf' # GRAPH SETTINGS: @@ -81,7 +64,7 @@ pure_grid_kwargs = dict( ncols=None, wspace=0.05, hspace=0.1, - left=0.13, + left=0.07, right=0.95, bottom=0.15, top=0.9 @@ -91,7 +74,7 @@ noise_grid_kwargs = dict( ncols=None, wspace=0.05, hspace=0.1, - left=0.13, + left=0.07, right=0.95, bottom=0.15, top=0.9 @@ -114,7 +97,10 @@ snip_specs = dict( # PLOT SETTINGS: colors = load_colors('../data/stage_colors.npz') -lw_snippets = 0.5 +lw_snippets = dict( + conv=0.5, + feat=2 +) lw_analysis = 3 xlabels = dict( analysis='scale $\\alpha$', @@ -166,10 +152,10 @@ letter_analysis_kwargs = dict( ) bar_time = 5 bar_kwargs = dict( - y0=0.5, - y1=0.6, + y0=0.7, + y1=0.8, color='k', - lw = 0, + lw=0, ) spread_kwargs = dict( alpha=0.3, @@ -183,8 +169,8 @@ for data_path in data_paths: print(f'Processing {data_path}') # Load invariance data: - pure_data, config = load_data(data_path, files) - noise_data, _ = load_data(data_path.replace('.npz', '_noise.npz'), files) + pure_data, config = load_data(data_path, **load_kwargs) + noise_data, _ = load_data(data_path.replace('.npz', '_noise.npz'), **load_kwargs) t_full = np.arange(pure_data['conv'].shape[0]) / config['env_rate'] # Reduce snippet data to kernel subset: @@ -237,7 +223,7 @@ for data_path in data_paths: # Plot pure-song kernel response snippets: plot_snippets(pure_axes[snip_specs['conv']], t_full, pure_data['conv'], - ymin=0, c=colors['conv'], lw=lw_snippets) + c=colors['conv'], lw=lw_snippets['conv']) # Plot pure-song binary snippets: plot_bi_snippets(pure_axes[snip_specs['bi']], t_full, pure_data['bi'], @@ -245,14 +231,14 @@ for data_path in data_paths: # Plot pure-song feature snippets: plot_snippets(pure_axes[snip_specs['feat']], t_full, pure_data['feat'], - c=colors['feat'], lw=lw_snippets) + ymin=0, ymax=1, c=colors['feat'], lw=lw_snippets['feat']) # Indicate time scale: time_bar(pure_axes[snip_specs['conv']][0], bar_time, **bar_kwargs) # Plot noise-song kernel response snippets: plot_snippets(noise_axes[snip_specs['conv']], t_full, noise_data['conv'], - ymin=0, c=colors['conv'], lw=lw_snippets) + c=colors['conv'], lw=lw_snippets['conv']) # Plot noise-song binary snippets: plot_bi_snippets(noise_axes[snip_specs['bi']], t_full, noise_data['bi'], @@ -260,7 +246,7 @@ for data_path in data_paths: # Plot noise-song feature snippets: plot_snippets(noise_axes[snip_specs['feat']], t_full, noise_data['feat'], - c=colors['feat'], lw=lw_snippets) + ymin=0, ymax=1, c=colors['feat'], lw=lw_snippets['feat']) # Indicate time scale: time_bar(noise_axes[snip_specs['conv']][0], bar_time, **bar_kwargs) diff --git a/python/fig_invariance_thresh-lp_single.py b/python/fig_invariance_thresh-lp_single.py new file mode 100644 index 0000000..d6f85b2 --- /dev/null +++ b/python/fig_invariance_thresh-lp_single.py @@ -0,0 +1,323 @@ +import plotstyle_plt +import numpy as np +import matplotlib.pyplot as plt +from itertools import product +from thunderhopper.filetools import search_files +from thunderhopper.modeltools import load_data +from thunderhopper.filtertools import find_kern_specs +from color_functions import load_colors +from plot_functions import hide_axis, xlimits, ylimits, xlabel, ylabel, super_ylabel,\ + plot_line, plot_barcode, strip_zeros, time_bar +from IPython import embed + +def add_snip_axes(fig, grid_kwargs): + grid = fig.add_gridspec(**grid_kwargs) + axes = np.zeros((grid.nrows, grid.ncols), dtype=object) + for i, j in product(range(grid.nrows), range(grid.ncols)): + axes[i, j] = fig.add_subplot(grid[i, j]) + [hide_axis(ax, 'left') for ax in axes.flatten()] + [hide_axis(ax, 'bottom') for ax in axes.flatten()] + return axes + +def plot_snippets(axes, time, snippets, ymin=None, ymax=None, **kwargs): + ymin, ymax = ylimits(snippets, minval=ymin, maxval=ymax, pad=0.05) + for ax, snippet in zip(axes, snippets.T): + plot_line(ax, time, snippet, ymin=ymin, ymax=ymax, **kwargs) + return None + +def plot_bi_snippets(axes, time, binary, **kwargs): + for ax, binary in zip(axes, binary.T): + plot_barcode(ax, time, binary[:, None], **kwargs) + return None + +def side_distributions(axes, snippets, inset_bounds, thresh, + ymin=None, ymax=None): + bins = np.linspace(snippets.min(), snippets.max(), 50) + centers = bins[:-1] + (bins[1] - bins[0]) / 2 + for ax, snippet in zip(axes, snippets.T): + inset = ax.inset_axes(inset_bounds) + inset.axis('off') + pdf, _ = np.histogram(snippet, bins, density=True) + inset.plot(pdf, centers, c='k', lw=1) + inset.fill_betweenx(centers, pdf.min(), pdf, where=(centers > thresh), + color=colors['bi'], lw=0) + inset.set_xlim(0, pdf.max()) + ylimits(centers, inset, minval=ymin, maxval=ymax, pad=0) + return None + + +# GENERAL SETTINGS: +with_noise = True +target = 'Omocestus_rufipes' +search_kwargs = dict( + incl='subset' if not with_noise else 'subset_noise', + dir='../data/inv/thresh_lp/' +) +data_paths = search_files(target, **search_kwargs) +stages = ['conv', 'bi', 'feat'] +load_kwargs = dict( + files=stages, + keywords=['scales', 'snip', 'measure', 'thresh'] +) +save_path = None#'../figures/fig_invariance_thresh_lp_single' +if with_noise and save_path is not None: + save_path += '_noise' + +# GRAPH SETTINGS: +fig_kwargs = dict( + figsize=(32/2.54, 16/2.54), +) +super_grid_kwargs = dict( + nrows=None, + ncols=2, + wspace=0, + hspace=0, + left=0, + right=1, + bottom=0, + top=1 +) +snip_grid_kwargs = dict( + nrows=len(stages), + ncols=None, + wspace=0.11, + hspace=0.1, + left=0.1, + right=0.95, + bottom=0.01, + top=0.85 +) +big_grid_kwargs = dict( + nrows=1, + ncols=1, + wspace=0, + hspace=0, + left=0.15, + right=0.96, + bottom=0.1, + top=0.99 +) +inset_bounds = [1, 0, 0.1, 1] + +# PLOT SETTINGS: +colors = load_colors('../data/stage_colors.npz') +# lw_snippets = dict( +# conv=0.5, +# feat=2 +# ) +# lw_analysis = 3 +xlabels = dict( + big='scale $\\alpha$', +) +xlab_big_kwargs = dict( + y=0.01, + fontsize=16, + ha='center', + va='bottom', +) +ylabels = dict( + conv='$c_i$', + bi='$b_i$', + feat='$f_i$', + big='$\\mu_f$', +) +ylab_snip_kwargs = dict( + x=0.08, + fontsize=20, + rotation=0, + ha='right', + va='center', +) +ylab_super_kwargs = dict( + x=0.005, + fontsize=16, + ha='left', + va='center', +) +ylab_big_kwargs = dict( + x=0.02, + fontsize=20, + ha='center', + va='top', +) +# xloc = dict( +# analysis=10, +# ) +# letter_snip_kwargs = dict( +# x=0.02, +# y=1, +# ha='left', +# va='top', +# fontsize=22, +# fontweight='bold' +# ) +# letter_analysis_kwargs = dict( +# x=0, +# y=1, +# ha='left', +# va='top', +# fontsize=22, +# fontweight='bold' +# ) +# bar_time = 5 +# bar_kwargs = dict( +# y0=0.7, +# y1=0.8, +# color='k', +# lw=0, +# ) +kernel = np.array([ + [2, 0.008], + [4, 0.008], +])[:1] +zoom_rel = np.array([0.5, 0.55]) + + +# EXECUTION: +for data_path in data_paths: + print(f'Processing {data_path}') + + # Load invariance data: + data, config = load_data(data_path, **load_kwargs) + t_full = np.arange(data['snip_conv'].shape[0]) / config['env_rate'] + zoom_abs = zoom_rel * t_full[-1] + zoom_inds = (t_full >= zoom_abs[0]) & (t_full <= zoom_abs[1]) + kern_ind = find_kern_specs(config['k_specs'], kerns=kernel)[0] + + # Reduce to kernel subset and crop time to zoom frame: + data['snip_conv'] = data['snip_conv'][zoom_inds, kern_ind, ...] + data['snip_bi'] = data['snip_bi'][zoom_inds, kern_ind, ...] + data['snip_feat'] = data['snip_feat'][zoom_inds, kern_ind, ...] + data['measure_conv'] = data['measure_conv'][:, kern_ind, :] + data['measure_feat'] = data['measure_feat'][:, kern_ind, :] + data['threshs'] = data['threshs'][:, kern_ind] + t_full = np.arange(data['snip_conv'].shape[0]) / config['env_rate'] + + # Adjust grid parameters: + super_grid_kwargs['nrows'] = data['thresh_perc'].size + snip_grid_kwargs['ncols'] = data['example_scales'].size + + # Prepare overall graph: + fig = plt.figure(**fig_kwargs) + super_grid = fig.add_gridspec(**super_grid_kwargs) + + # Prepare analysis axis: + big_subfig = fig.add_subfigure(super_grid[slice(None), 1]) + big_grid = big_subfig.add_gridspec(**big_grid_kwargs) + big_ax = big_subfig.add_subplot(big_grid[0, 0]) + xlabel(big_ax, xlabels['big'], **xlab_big_kwargs, + transform=big_subfig.transSubfigure) + ylabel(big_ax, ylabels['big'], **ylab_big_kwargs, + transform=big_subfig.transSubfigure) + big_ax.set_xlim(data['scales'].min(), data['scales'].max()) + ylimits(data['measure_feat'], big_ax, minval=0, pad=0.05) + big_ax.set_xscale('symlog', linthresh=data['scales'][1], linscale=0.5) + + # Prepare snippet axes: + snip_axes = {} + for i in range(data['thresh_perc'].size): + snip_subfig = fig.add_subfigure(super_grid[i, 0]) + axes = add_snip_axes(snip_subfig, snip_grid_kwargs) + snip_axes[snip_subfig] = axes + super_ylabel(f'{data["thresh_perc"][i]}%', snip_subfig, + axes[0, 0], axes[-1, 0], **ylab_super_kwargs) + for ax, stage in zip(axes[:, 0], stages): + ylabel(ax, ylabels[stage], **ylab_snip_kwargs, + transform=snip_subfig.transSubfigure) + if i == 0: + for ax, scale in zip(axes[0, :], data['example_scales']): + ax.set_title(f'$\\alpha={strip_zeros(scale)}$') + + # Plot representation snippets per threshold: + for i, (subfig, axes) in enumerate(snip_axes.items()): + # Plot kernel response snippets: + plot_snippets(axes[0, :], t_full, data['snip_conv'][:, :, i], + c=colors['conv'], lw=0.5) + # Plot binary snippets: + plot_bi_snippets(axes[1, :], t_full, data['snip_bi'][:, :, i], + color=colors['bi'], lw=0) + # Plot feature snippets: + plot_snippets(axes[2, :], t_full, data['snip_feat'][:, :, i], + ymin=0, ymax=1, c=colors['feat'], lw=2) + + # Plot kernel response distribution: + side_distributions(axes[0, :], data['snip_conv'][:, :, i], + inset_bounds, data['threshs'][i]) + + # Plot analysis results: + big_ax.plot(data['scales'], data['measure_feat'], + c=colors['feat'], lw=3) + + + # # Prepare pure-song snippet axes: + # pure_subfig = fig.add_subfigure(super_grid[subfig_specs['pure']]) + # pure_grid_kwargs['nrows' if pure_grid_kwargs['nrows'] is None else 'ncols'] = pure_data['example_scales'].size + # pure_axes = add_snip_axes(pure_subfig, pure_grid_kwargs) + # for ax, stage in zip(pure_axes[:, 0], stages): + # ylabel(ax, ylabels[stage], **ylab_snip_kwargs, + # transform=pure_subfig.transSubfigure) + # for ax, scale in zip(pure_axes[snip_specs['conv']], pure_data['example_scales']): + # ax.set_title(f'$\\alpha={strip_zeros(scale)}$') + # pure_subfig.text(s='a', **letter_snip_kwargs) + + # # Prepare analysis axis: + # analysis_subfig = fig.add_subfigure(super_grid[subfig_specs['analysis']]) + # analysis_grid = analysis_subfig.add_gridspec(**analysis_grid_kwargs) + # analysis_ax = analysis_subfig.add_subplot(analysis_grid[0, 0]) + # analysis_ax.set_xlim(noise_data['scales'].min(), noise_data['scales'].max()) + # analysis_ax.xaxis.set_major_locator(plt.MultipleLocator(xloc['analysis'])) + # xlabel(analysis_ax, xlabels['analysis'], **xlab_analysis_kwargs, + # transform=analysis_subfig.transSubfigure) + # # analysis_ax.set_yscale('log') + # ylabel(analysis_ax, ylabels['analysis'], **ylab_analysis_kwargs, + # transform=analysis_subfig.transSubfigure) + # analysis_subfig.text(s='c', **letter_analysis_kwargs) + + # # Plot pure-song kernel response snippets: + # plot_snippets(pure_axes[snip_specs['conv']], t_full, pure_data['conv'], + # c=colors['conv'], lw=lw_snippets['conv']) + + # # Plot pure-song binary snippets: + # plot_bi_snippets(pure_axes[snip_specs['bi']], t_full, pure_data['bi'], + # color=colors['bi'], lw=0) + + # # Plot pure-song feature snippets: + # plot_snippets(pure_axes[snip_specs['feat']], t_full, pure_data['feat'], + # ymin=0, ymax=1, c=colors['feat'], lw=lw_snippets['feat']) + + # # Indicate time scale: + # time_bar(pure_axes[snip_specs['conv']][0], bar_time, **bar_kwargs) + + # # Plot noise-song kernel response snippets: + # plot_snippets(noise_axes[snip_specs['conv']], t_full, noise_data['conv'], + # c=colors['conv'], lw=lw_snippets['conv']) + + # # Plot noise-song binary snippets: + # plot_bi_snippets(noise_axes[snip_specs['bi']], t_full, noise_data['bi'], + # color=colors['bi'], lw=0) + + # # Plot noise-song feature snippets: + # plot_snippets(noise_axes[snip_specs['feat']], t_full, noise_data['feat'], + # ymin=0, ymax=1, c=colors['feat'], lw=lw_snippets['feat']) + + # # Indicate time scale: + # time_bar(noise_axes[snip_specs['conv']][0], bar_time, **bar_kwargs) + + # # Plot noise-song SD ratios (limited): + # analysis_ax.plot(noise_data['scales'], noise_data['measure_conv'], + # c=colors['conv'], lw=lw_analysis) + # lower, upper = noise_data['spread_conv'] + # analysis_ax.fill_between(noise_data['scales'], lower, upper, + # color=colors['conv'], **spread_kwargs) + # analysis_ax.plot(noise_data['scales'], noise_data['measure_feat'], + # c=colors['feat'], lw=lw_analysis) + # lower, upper = noise_data['spread_feat'] + # analysis_ax.fill_between(noise_data['scales'], lower, upper, + # color=colors['feat'], **spread_kwargs) + + if save_path is not None: + fig.savefig(save_path) + plt.show() + +print('Done.') +embed() diff --git a/python/plot_functions.py b/python/plot_functions.py index 34c407f..4c995ff 100644 --- a/python/plot_functions.py +++ b/python/plot_functions.py @@ -1,6 +1,7 @@ import string import numpy as np import matplotlib.pyplot as plt +from matplotlib.transforms import BboxTransformTo from itertools import product def prepare_fig(nrows, ncols, width=8, height=None, rheight=2, unit=1/2.54, @@ -154,3 +155,19 @@ def strip_zeros(num, right_digits=5): return f'{left}.{right}' return left +def time_bar(ax, dur, y0=0.9, y1=0.95, xshift=0.5, parent=None, transform=None, **kwargs): + t0, t1 = ax.get_xlim() + offset = (t1 - t0 - dur) * xshift + x0 = t0 + offset + x1 = x0 + dur + if parent is None: + parent = ax + if transform is None: + transform = BboxTransformTo(parent.bbox) + if transform is not ax.transData: + trans = ax.transData + transform.inverted() + x0 = trans.transform((x0, 0))[0] + x1 = trans.transform((x1, 0))[0] + parent.add_artist(plt.Rectangle((x0, y0), x1 - x0, y1 - y0, + transform=transform, **kwargs)) + return None diff --git a/python/save_inv_data_full.py b/python/save_inv_data_full.py new file mode 100644 index 0000000..1771712 --- /dev/null +++ b/python/save_inv_data_full.py @@ -0,0 +1,115 @@ +import glob +import numpy as np +import matplotlib.pyplot as plt +from thunderhopper.modeltools import load_data, save_data +from thunderhopper.filetools import crop_paths +from thunderhopper.model import process_signal +from IPython import embed + +# GENERAL SETTINGS: +target = 'Omocestus_rufipes' +data_paths = glob.glob(f'../data/processed/{target}*.npz') +stages = ['raw', 'filt', 'env', 'log', 'inv', 'conv', 'bi', 'feat'] +save_path = '../data/inv/full/' + +# ANALYSIS SETTINGS: +example_scales = np.array([0, 0.5, 1, 5, 10]) +scales = np.linspace(0, 10, 100) +scales = np.unique(np.concatenate((scales, example_scales))) + +# EXECUTION: +for data_path, name in zip(data_paths, crop_paths(data_paths)): + print(f'Processing {name}') + + # Get normalized song recording: + data, config = load_data(data_path, files='raw') + song, rate = data['raw'], config['rate'] + song /= song.std(axis=0) + + # Get normalized noise: + rng = np.random.default_rng() + noise = rng.normal(size=song.shape[0]) + noise /= noise.std() + + # Get song segment to be analyzed: + time = np.arange(song.shape[0]) / rate + start, end = data['songs_0'].ravel() + segment = (time >= start) & (time <= end) + + # Prepare snippet storage: + shape_low = (song.shape[0], example_scales.size) + shape_high = (song.shape[0], config['k_specs'].shape[0], example_scales.size) + snippets = dict( + raw=np.zeros(shape_low, dtype=float), + filt=np.zeros(shape_low, dtype=float), + env=np.zeros(shape_low, dtype=float), + log=np.zeros(shape_low, dtype=float), + inv=np.zeros(shape_low, dtype=float), + conv=np.zeros(shape_high, dtype=float), + bi=np.zeros(shape_high, dtype=float), + feat=np.zeros(shape_high, dtype=float) + ) + + # Prepare measure storage: + shape_low = (scales.size,) + shape_high = (scales.size, config['k_specs'].shape[0]) + measures = dict( + measure_raw=np.zeros(shape_low, dtype=float), + measure_filt=np.zeros(shape_low, dtype=float), + measure_env=np.zeros(shape_low, dtype=float), + measure_log=np.zeros(shape_low, dtype=float), + measure_inv=np.zeros(shape_low, dtype=float), + measure_conv=np.zeros(shape_high, dtype=float), + measure_feat=np.zeros(shape_high, dtype=float) + ) + + # Execute piecewise: + for i, scale in enumerate(scales): + print('Simulating scale ', scale) + + # Rescale song and add noise: + scaled = song * scale + noise + + # Process mixture: + signals, rates = process_signal(config, returns=stages, + signal=scaled, rate=rate) + # Store results: + for stage in stages: + key = f'measure_{stage}' + + # Log snippet data: + if scale in example_scales: + scale_ind = np.nonzero(example_scales == scale)[0][0] + snippets[stage][:, ..., scale_ind] = signals[stage] + + # Log "intensity measure" per stage: + if stage in ['raw', 'filt', 'env', 'log', 'inv', 'conv']: + measures[key][i] = signals[stage][segment, ...].std(axis=0) + elif stage == 'feat': + measures[key][i] = signals[stage][segment, :].mean(axis=0) / signals[stage][segment, :].std(axis=0) + + # Relate to smallest scale: + base_ind = np.argmin(scales) + for stage in stages: + if stage == 'bi': + continue + key = f'measure_{stage}' + measures[key] /= measures[key][base_ind, ...] + if stage in ['conv', 'feat']: + spread = np.zeros((2, scales.size)) + spread[0] = np.percentile(measures[key], 25, axis=1) + spread[1] = np.percentile(measures[key], 75, axis=1) + measures[f'spread_{stage}'] = spread + measures[key] = np.median(measures[key], axis=1) + + # Save analysis results: + if save_path is not None: + data = dict( + scales=scales, + example_scales=example_scales, + ) + data.update(snippets) + data.update(measures) + save_data(save_path + name, data, config, overwrite=True) +print('Done.') +embed() diff --git a/python/save_inv_data_log-hp.py b/python/save_inv_data_log-hp.py index 0fc42a5..cfef28a 100644 --- a/python/save_inv_data_log-hp.py +++ b/python/save_inv_data_log-hp.py @@ -12,25 +12,29 @@ data_paths = glob.glob(f'../data/processed/{target}*.npz') save_path = '../data/inv/log_hp/' # ANALYSIS SETTINGS: -add_noise = False +add_noise = True single_db_ref = True -# find_saturation = add_noise and False example_scales = np.array([0, 0.1, 1, 10, 100, 200]) scales = np.geomspace(0.1, 1000, 100) if not add_noise: example_scales = example_scales[example_scales > 0] scales = np.unique(np.concatenate((scales, example_scales))) -# if find_saturation: -# scales = np.append(scales, 10e10) # EXECUTION: for data_path, name in zip(data_paths, crop_paths(data_paths)): print(f'Processing {name}') - # Get normalized song envelope: + # Get song envelope: data, config = load_data(data_path, files='env') song, rate = data['env'], config['env_rate'] - song /= song.std() + + # Get song segment to be analyzed: + time = np.arange(song.shape[0]) / rate + start, end = data['songs_0'].ravel() + segment = (time >= start) & (time <= end) + + # Normalize song component: + song /= song[segment].std() # Rescale song component: mix = song[:, None] * scales[None, :] @@ -40,7 +44,7 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)): rng = np.random.default_rng() noise = rng.normal(size=song.shape) noise = extract_env(noise, rate, config=config) - noise /= noise.std() + noise /= noise[segment].std() mix += noise[:, None] # Process mixture: @@ -49,17 +53,9 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)): padtype='constant', padlen=config['padlen']) # Get "intensity measure" per stage: - measure_env = mix.std(axis=0) - measure_log = mix_log.std(axis=0) - measure_inv = mix_inv.std(axis=0) - - # # Find saturation level: - # if find_saturation: - # limit = measure_inv[-1] - # scales = scales[:-1] - # measure_env = measure_env[:-1] - # measure_log = measure_log[:-1] - # measure_inv = measure_inv[:-1] + measure_env = mix[segment, :].std(axis=0) + measure_log = mix_log[segment, :].std(axis=0) + measure_inv = mix_inv[segment, :].std(axis=0) # Save analysis results: save_inds = np.nonzero(np.isin(scales, example_scales))[0] @@ -74,8 +70,6 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)): measure_log=measure_log, measure_inv=measure_inv, ) - # if find_saturation: - # data['limit'] = limit file_name = save_path + name if add_noise: file_name += '_noise' diff --git a/python/save_inv_data_thresh-lp.py b/python/save_inv_data_thresh-lp.py index 7c8d83f..f4eea6f 100644 --- a/python/save_inv_data_thresh-lp.py +++ b/python/save_inv_data_thresh-lp.py @@ -1,5 +1,6 @@ import glob import numpy as np +import matplotlib.pyplot as plt from thunderhopper.modeltools import load_data, save_data from thunderhopper.filetools import crop_paths from thunderhopper.filters import sosfilter @@ -12,30 +13,40 @@ save_path = '../data/inv/thresh_lp/' # ANALYSIS SETTINGS: add_noise = False -threshold = 0.5 -example_scales = np.array([threshold, 0.6, 1, 10, 50, 100]) -scales = np.linspace(threshold + 0.1, 100, 100) -if not add_noise: - example_scales = example_scales[example_scales > threshold] +thresh_percent = 90 +example_scales = np.array([0, 0.5, 1, 10, 50]) +scales = np.geomspace(0.01, 50, 100) scales = np.unique(np.concatenate((scales, example_scales))) +plot_results = True # EXECUTION: for data_path, name in zip(data_paths, crop_paths(data_paths)): print(f'Processing {name}') + save_name = save_path + name - # Get normalized pure-song kernel responses: + # Get pure-song kernel responses: data, config = load_data(data_path, files='conv') song, rate = data['conv'], data['conv_rate'] - song /= song.std(axis=0) - # Prepare kernel-specific thresholds: - threshold *= song.max(axis=0, keepdims=True) + # Get song segment to be analyzed: + time = np.arange(song.shape[0]) / rate + start, end = data['songs_0'].ravel() + segment = (time >= start) & (time <= end) + + # Normalize song component: + song /= song[segment, :].std(axis=0) if add_noise: # Get normalized noise: rng = np.random.default_rng() noise = rng.normal(size=(song.shape[0], 1)) - noise /= noise.std() + noise /= noise[segment].std() + + # Prepare noise-bound threshold: + threshold = np.percentile(noise, thresh_percent, axis=0) + else: + # Reuse threshold from previous noise run: + threshold = np.load(save_name + '_noise.npz')['thresh'] # Prepare snippet storage: shape = song.shape + (example_scales.size,) @@ -71,23 +82,32 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)): feat[:, :, scale_ind] = scaled_feat # Get "intensity measure" per stage: - measure_conv[i] = scaled_conv.std(axis=0) - measure_feat[i] = scaled_feat.mean(axis=0) + measure_conv[i] = scaled_conv[segment, :].std(axis=0) + measure_feat[i] = scaled_feat[segment, :].mean(axis=0) - # Relate to smallest scale: - base_ind = np.argmin(scales) - measure_conv /= measure_conv[base_ind, :] - measure_feat /= measure_feat[base_ind, :] + # # Relate to smallest scale: + # base_ind = np.argmin(scales) + # measure_conv /= measure_conv[base_ind, :] + + if plot_results: + fig, (ax1, ax2) = plt.subplots(2, 1) + ax1.plot(scales, measure_conv) + ax1.plot(scales, measure_conv.mean(axis=1), c='k') + ax1.plot(scales, np.median(measure_conv, axis=1), c='k', ls='--') + ax2.plot(scales, measure_feat) + ax2.plot(scales, np.nanmean(measure_feat, axis=1), c='k') + ax2.plot(scales, np.nanmedian(measure_feat, axis=1), c='k', ls='--') + plt.show() # Condense measures across kernels: spread_conv = np.zeros((2, scales.size)) - spread_conv[0] = np.percentile(measure_conv, 25, axis=1) - spread_conv[1] = np.percentile(measure_conv, 75, axis=1) - measure_conv = np.median(measure_conv, axis=1) + spread_conv[0] = np.nanpercentile(measure_conv, 25, axis=1) + spread_conv[1] = np.nanpercentile(measure_conv, 75, axis=1) + measure_conv = np.nanmedian(measure_conv, axis=1) spread_feat = np.zeros((2, scales.size)) - spread_feat[0] = np.percentile(measure_feat, 25, axis=1) - spread_feat[1] = np.percentile(measure_feat, 75, axis=1) - measure_feat = np.median(measure_feat, axis=1) + spread_feat[0] = np.nanpercentile(measure_feat, 25, axis=1) + spread_feat[1] = np.nanpercentile(measure_feat, 75, axis=1) + measure_feat = np.nanmedian(measure_feat, axis=1) # Save analysis results: if save_path is not None: @@ -101,10 +121,11 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)): spread_conv=spread_conv, measure_feat=measure_feat, spread_feat=spread_feat, + thresh=threshold, + thresh_perc=thresh_percent, ) - file_name = save_path + name if add_noise: - file_name += '_noise' - save_data(file_name, data, config, overwrite=True) + save_name += '_noise' + save_data(save_name, data, config, overwrite=True) print('Done.') embed() diff --git a/python/save_inv_data_thresh-lp_subset.py b/python/save_inv_data_thresh-lp_subset.py new file mode 100644 index 0000000..e88d937 --- /dev/null +++ b/python/save_inv_data_thresh-lp_subset.py @@ -0,0 +1,134 @@ +import glob +import numpy as np +import matplotlib.pyplot as plt +from thunderhopper.modeltools import load_data, save_data +from thunderhopper.filetools import crop_paths +from thunderhopper.filters import sosfilter +from thunderhopper.filtertools import find_kern_specs +from IPython import embed + +# GENERAL SETTINGS: +target = 'Omocestus_rufipes' +data_paths = glob.glob(f'../data/processed/{target}*.npz') +save_path = '../data/inv/thresh_lp/' + +# ANALYSIS SETTINGS: +add_noise = True +thresh_percent = np.array([50, 75, 100]) +example_scales = np.array([0, 1, 10, 50]) +scales = np.geomspace(0.01, 100, 100) +scales = np.unique(np.concatenate((scales, example_scales))) +plot_results = False +kernels = np.array([ + [2, 0.008], + [4, 0.008], +]) + +# EXECUTION: +for data_path, name in zip(data_paths, crop_paths(data_paths)): + print(f'Processing {name}') + save_name = save_path + name + '_subset' + + # Get pure-song kernel responses: + data, config = load_data(data_path, files='conv') + conv, rate = data['conv'], data['conv_rate'] + + # Get song segment to be analyzed: + time = np.arange(conv.shape[0]) / rate + start, end = data['songs_0'].ravel() + segment = (time >= start) & (time <= end) + + # Reduce to kernel subset: + kern_inds = find_kern_specs(config['k_specs'], kerns=kernels) + 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] + conv = conv[:, kern_inds] + + # Normalize kernel responses: + conv /= conv[segment, :].std(axis=0) + + if add_noise: + # Get normalized noise: + rng = np.random.default_rng() + noise = rng.normal(size=(conv.shape[0], 1)) + noise /= noise[segment].std() + + # Prepare snippet storage: + shape = conv.shape + (example_scales.size, thresh_percent.size) + snip_conv = np.zeros(shape, dtype=float) + snip_bi = np.zeros(shape, dtype=float) + snip_feat = np.zeros(shape, dtype=float) + + # Prepare measure storage: + shape = (scales.size, conv.shape[1], thresh_percent.size) + measure_conv = np.zeros(shape, dtype=float) + measure_feat = np.zeros(shape, dtype=float) + + # Compute kernel-specific thresholds (thresholds, kernels): + thresholds = np.percentile(conv, thresh_percent, axis=0) + + # Execute piecewise analysis: + for i, threshs in enumerate(thresholds): + print('\nSimulating threshold ', thresh_percent[i]) + + for j, scale in enumerate(scales): + print('Simulating scale ', scale) + + # Rescale conv component: + scaled_conv = conv * scale + if add_noise: + # Add noise: + scaled_conv += noise + + # Process mixture: + scaled_bi = (scaled_conv > threshs).astype(float) + scaled_feat = sosfilter(scaled_bi, rate, config['feat_fcut'], 'lp', + padtype='fixed', padlen=config['padlen']) + + # Log snippet data: + if scale in example_scales: + scale_ind = np.nonzero(example_scales == scale)[0][0] + snip_conv[:, :, scale_ind, i] = scaled_conv + snip_bi[:, :, scale_ind, i] = scaled_bi + snip_feat[:, :, scale_ind, i] = scaled_feat + + # Get intensity measure per stage: + measure_conv[j, :, i] = scaled_conv[segment, :].std(axis=0) + measure_feat[j, :, i] = scaled_feat[segment, :].mean(axis=0) + + if plot_results: + fig, axes = plt.subplots(thresh_percent.size, kernels.shape[0], + figsize=(16, 9), layout='constrained', + sharex=True, sharey=True, squeeze=True) + axes[0, 0].set_xscale('symlog', linthresh=scales[scales>0].min(), + linscale=0.25) + + for i, thresh in enumerate(thresh_percent): + for j, kernel in enumerate(kernels): + ax = axes[i, j] + ax.plot(scales, measure_feat[:, j, i], 'k') + if i == 0: + ax.set_title(f'Kernel {kernel}') + if j == 0: + ax.set_ylabel(f'{thresh}%') + plt.show() + + # Save analysis results: + if save_path is not None: + data = dict( + scales=scales, + example_scales=example_scales, + snip_conv=snip_conv, + snip_bi=snip_bi, + snip_feat=snip_feat, + measure_conv=measure_conv, + measure_feat=measure_feat, + thresh_perc=thresh_percent, + threshs=thresholds, + ) + if add_noise: + save_name += '_noise' + save_data(save_name, data, config, overwrite=True) +print('Done.') +embed() diff --git a/python/save_snippet_data.py b/python/save_snippet_data.py index b9bbfa9..6d0e55a 100644 --- a/python/save_snippet_data.py +++ b/python/save_snippet_data.py @@ -7,11 +7,10 @@ from IPython import embed ## SETTINGS: # General: -overwrite = True input_folder = '../data/raw/' output_folder = '../data/processed/' stages = ['raw', 'filt', 'env', 'log', 'inv', 'conv', 'bi', 'feat', 'norm'] -if True: +if False: # Overwrites edited: stages.append('songs') @@ -51,8 +50,7 @@ for path, name in zip(input_paths, path_names): # Fetch and store representations: save = None if output_folder is None else output_folder + f'{name}.npz' - process_signal(config, stages, path, save=save, - label_edit=gui, overwrite=overwrite) + process_signal(config, stages, path, save=save, label_edit=gui) # Cross-control: if reload_saved: