V1_connectivity_rules_PNAS/Figure_2_Code_and_Data/panels.py

"""Create panels for figures. Helper module."""

import numpy as np
from .util import (stackmean, get_idx_match, contrast_log_scale, replace_imshow_ticklabels,
                   ori_con_2d, circ_gaussian, hyper_ratio, wrapped_gaussian, von_mises)
import matplotlib.pyplot as plt
from cycler import cycler
import datajoint as dj
from mpl_toolkits.axes_grid1.inset_locator import InsetPosition
from statsmodels.distributions.empirical_distribution import ECDF
from scipy.stats import sem


def kernels(kerns, dtopt, dtrange, ax=None, key=None, prop='grat_orientation'):
    if ax is None:
        fig, ax = plt.subplots()

    kern = np.stack(kerns['t_kernel'].groupby(prop).apply(stackmean))
    propvals = np.unique(kerns.index.get_level_values(prop))
    if prop == 'grat_orientation':
        kern = np.concatenate((kern, kern[0][np.newaxis, :]))
        propvals = np.concatenate((propvals, np.array([180]))).astype(int)
    im = ax.imshow(kern[:, :30]*1000, origin='lower', aspect='auto')
    dtidx = np.where(np.isclose(dtrange, dtopt))[0]
    ax.vlines(dtidx, -0.5, kern.shape[0]-0.5)
    ax.set_ylabel(prop[5].capitalize() + prop[6:])
    if plt.rcParams['text.usetex']:
        ax.set_xlabel(r'$\delta t$ (s)')
    else:
        ax.set_xlabel(r'dt (s)')
    replace_imshow_ticklabels(ax, dtrange[:30], propvals)
    cbar = plt.colorbar(im, ax=ax)
    cbar.set_label('Firing (Hz)', rotation=270, labelpad=7)
    return cbar


def opt_dt_distr_exc_inh(exc_t_opts, inh_t_opts, t_opts, ax=None, excitcol='darkorange',
                         inhibcol='teal'):
    if ax is None:
        fig, ax = plt.subplots()

    _, bins, _ = ax.hist(exc_t_opts, color=excitcol, bins='doane')
    ax.hist(inh_t_opts, color=inhibcol, bins=bins)
    xlim = ax.get_xlim()
    ax.set_xlim([0, xlim[1]])
    ax7a = ax.twinx()
    ecdf = ECDF(t_opts)
    ax7a.plot(ecdf.x, ecdf.y, color='k')
    ax7a.set_ylim([0, 1.01])
    if plt.rcParams['text.usetex']:
        ax.set_xlabel(r'Optimal $\delta t$ (s)')
    else:
        ax.set_xlabel('Optimal dt (s)')
    ax.set_ylabel('N')
    ax7a.set_ylabel('Proportion')
    ax7a.set_yticks([0, 0.5, 1])
    ax7a.spines['right'].set_bounds(0, 1)
    ax7a.spines['right'].set_visible(True)
    ax.spines['bottom'].set_bounds(*ax.get_xlim())


def oricon_matrix(ax, data, cbar=False):
    df = data
    m = df.pivot_table('p_spike_stim', 'grat_orientation', 'grat_contrast')
    m = np.concatenate((m.values, m.values[0][np.newaxis, :]))*1000  # convert to Hz
    im = ax.imshow(m, origin='lower')
    if cbar:
        bar = plt.colorbar(im, ax=ax)
        bar.set_label('Firing rate (Hz)')
    else:
        bar = None
    ax.set_xlabel('Contrast')
    ax.set_ylabel('Orientation ($^{\circ}$)')
    contrasts = np.unique(df.index.get_level_values('grat_contrast')).round(2)
    orientations = np.linspace(0, 180, 13, endpoint=True).astype(int)
    replace_imshow_ticklabels(ax, contrasts, orientations)
    return ax, bar


def wave_shapes(wforms, ax=None, excitcol='darkorange', inhibcol='teal', annotate=False):
    if ax is None:
        fig, ax = plt.subplots()
    wforms['wave'] = wforms.apply(
        lambda x: x['u_wave'][np.where(x['u_chans'] == x['u_maxchan'])[0][0]], axis=1)
    wforms['normwave'] = wforms.wave.apply(lambda x: x / np.abs(np.min(x)))
    wforms['ip_wave'] = wforms.normwave.apply(
        lambda x: np.interp(np.linspace(0, len(x), 120), np.arange(len(x)), x))
    wforms['al_wave'] = wforms.ip_wave.apply(lambda x: np.roll(x, 29 - np.argmin(x)))
    wtime = np.linspace(0, 120 * 0.017, 120)

    colors = [excitcol, inhibcol]
    for c, (gname, g) in zip(colors, wforms.groupby('unit_type')):
        ax.plot(wtime, np.stack(g['al_wave']).T, color=c, alpha=0.1)
    exmean = np.mean(wforms.query('unit_type=="excit"')['al_wave'], axis=0)
    inmean = np.mean(wforms.query('unit_type=="inhib"')['al_wave'], axis=0)
    ax.plot(wtime, exmean, color=excitcol)
    ax.plot(wtime, inmean, color=inhibcol)
    ax.set_xlabel('Time (ms)')
    ax.set_ylabel('Voltage')
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)
    ax.spines['left'].set_bounds(ax.dataLim.ymin, ax.dataLim.ymax)
    if annotate:
        ax.vlines(wtime[np.argmin(inmean)], np.min(inmean), 1.2*np.max(inmean),
                  linestyles='dashed', color='k', zorder=3, linewidth=1)
        ax.vlines(wtime[np.argmax(inmean)], np.max(inmean), 1.2*np.max(inmean),
                  linestyles='dashed', color='k', zorder=3, linewidth=1)
        ax.hlines(1.2*np.max(inmean), wtime[np.argmin(inmean)], wtime[np.argmax(inmean)],
                  color='k', zorder=3, linewidth=1)
        peak = np.where(inmean >= 0.5*np.max(inmean))[0]
        ax.vlines(wtime[peak[0]], 0.5*np.max(inmean), -0.5, linestyles='dashed', color='k',
                  zorder=3, linewidth=1)
        ax.vlines(wtime[peak[-1]], 0.5*np.max(inmean), -0.5, linestyles='dashed', color='k',
                  zorder=3, linewidth=1)
        ax.hlines(-0.5, wtime[peak[0]], wtime[peak[-1]], color='k', zorder=3, linewidth=1)


def wave_shapes_lgn(wforms, ax=None, color='k', annotate=False):
    if ax is None:
        fig, ax = plt.subplots()
    wforms['wave'] = wforms.apply(
        lambda x: x['u_wave'][np.where(x['u_chans'] == x['u_maxchan'])[0][0]], axis=1)
    wforms['normwave'] = wforms.wave.apply(lambda x: x / np.abs(np.min(x)))
    wforms['ip_wave'] = wforms.normwave.apply(
        lambda x: np.interp(np.linspace(0, len(x), 120), np.arange(len(x)), x))
    wforms['al_wave'] = wforms.ip_wave.apply(lambda x: np.roll(x, 29 - np.argmin(x)))
    wtime = np.linspace(0, 120 * 0.017, 120)

    ax.plot(wtime, np.stack(wforms['al_wave']).T, color=color, alpha=0.1)
    mean = np.mean(wforms['al_wave'], axis=0)
    ax.plot(wtime, mean, color=color)
    ax.set_xlabel('Time (ms)')
    ax.set_ylabel('Voltage')
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)
    ax.spines['left'].set_bounds(ax.dataLim.ymin, ax.dataLim.ymax)


def wave_props(non_rc_props, rc_props, ax=None, excitcol='darkorange', inhibcol='teal', markeredgewidth=0.2,
               **plotkwargs):
    if ax is None:
        fig, ax = plt.subplots()
    colors = [excitcol, inhibcol]
    for c, (gname, g) in zip(colors, non_rc_props.groupby('unit_type')):
        ax.plot(g['fwhm_peak'], g['t2p_time'], '.', color=c, alpha=0.25, markersize=4, **plotkwargs)
    for c, (gname, g) in zip(colors, rc_props.groupby('unit_type')):
        ax.plot(g['fwhm_peak'], g['t2p_time'], '.', color=c,
                markeredgewidth=markeredgewidth, alpha=0.5, markersize=7, **plotkwargs)
    ax.set_xlabel('Fwhm_peak (ms)')
    ax.set_ylabel('T2p_time (ms)')
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)
    ax.spines['left'].set_bounds(ax.dataLim.ymin, ax.dataLim.ymax)


def wave_props_lgn(non_rc_props, rc_props, ax=None, color='k', markeredgewidth=0.2,
                   **plotkwargs):
    if ax is None:
        fig, ax = plt.subplots()
    ax.plot(non_rc_props['fwhm_peak'], non_rc_props['t2p_time'], '.', color=color, alpha=0.25, markersize=4, **plotkwargs)
    ax.plot(rc_props['fwhm_peak'], rc_props['t2p_time'], '.', color=color,
            markeredgewidth=markeredgewidth, alpha=0.5, markersize=7, **plotkwargs)
    ax.set_xlabel('Fwhm_peak (ms)')
    ax.set_ylabel('T2p_time (ms)')
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)
    ax.spines['left'].set_bounds(ax.dataLim.ymin, ax.dataLim.ymax)


def daterrsvd(datdf, svddf, dtrange, axs=None, cbar=False, cbaraxs=None):
    if axs is None:
        fig, (ax1, ax2, ax3) = plt.subplots(1, 3)
        fig.suptitle(key)
    else:
        assert len(axs) == 3
        ax1, ax2, ax3 = axs
    # prepare data
    dat = datdf.pivot_table('p_spike_stim', 'grat_orientation', 'grat_contrast').values
    # prepare svd
    svddf['p_svd'] = svddf.apply(lambda x: x.svd_s_kernel[get_idx_match(x.dt_maxresponse,
                                                                        dtrange)], axis=1)
    svddf['p_err'] = svddf.apply(lambda x: x.err_s_kernel[get_idx_match(x.dt_maxresponse,
                                                                        dtrange)], axis=1)
    svd = svddf.pivot_table('p_svd', 'grat_orientation', 'grat_contrast').values
    svderr = svddf.pivot_table('p_err', 'grat_orientation', 'grat_contrast').values
    # wrap orientation domain
    dat = np.vstack((dat, dat[0, :]))
    svd = np.vstack((svd, svd[0, :]))
    svderr = np.vstack((svderr, svderr[0, :]))
    # prepare ticklabels
    oris = np.linspace(0, 180, 13).astype(int)
    cons = np.unique(datdf.index.get_level_values('grat_contrast'))
    # plot
    ims = []
    for ax, mat in zip([ax1, ax2, ax3], [dat, svd, svderr]):
        im = ax.imshow(mat*1000, origin='lower')  # go from probability to Hz
        ims.append(im)
        replace_imshow_ticklabels(ax, cons, oris)
    # adjust labels
    ax2.set_yticklabels([])
    ax3.set_yticklabels([])
    ax1.set_ylabel('Orientation')
    for ax in [ax1, ax2, ax3]:
        ax.set_xlabel('Contrast')
    if cbar:
        if cbaraxs is None:
            cb1 = plt.colorbar(ims[1], ax=axs[1])
            cb2 = plt.colorbar(ims[2], ax=axs[2])
        else:
            cb1 = plt.colorbar(ims[1], cbaraxs[0])
            cb2 = plt.colorbar(ims[2], cbaraxs[1])
        return cb1, cb2


def svd_vs_spatcorr(df, ax=None, exc_color='darkorange',
                    inh_color='teal'):
    if ax is None:
        fig, ax = plt.subplots()
    inv_df = df.query('g_z_opt<=1.96 or power_opt>0.95')
    dep_df = df.query('g_z_opt>1.96 and power_opt<=0.95')
    e_inv_df = inv_df.query('unit_type=="excit"')
    i_inv_df = inv_df.query('unit_type=="inhib"')
    e_dep_df = dep_df.query('unit_type=="excit"')
    i_dep_df = dep_df.query('unit_type=="inhib"')
    ax.scatter(1-e_inv_df['power_opt'], e_inv_df['g_z_opt'], color=exc_color)
    ax.scatter(1-i_inv_df['power_opt'], i_inv_df['g_z_opt'], color=inh_color)
    ax.scatter(1-e_dep_df['power_opt'], e_dep_df['g_z_opt'], color=exc_color, alpha=0.5)
    ax.scatter(1-i_dep_df['power_opt'], i_dep_df['g_z_opt'], color=inh_color, alpha=0.5)
    ax.hlines(1.96, 0, 0.95, linestyles='dashed', alpha=0.5)
    ax.vlines(0.05, 1.96, np.max(df['g_z_opt']), linestyles='dashed',
              alpha=0.5)
    ax.spines['left'].set_bounds(ax.dataLim.ymin, ax.dataLim.ymax)
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)
    ax.set_ylabel('g_z_opz')
    ax.set_xlabel('svd_power')


def svd_vs_spatcorr_lgn(df, ax=None, color='k', draw_c_inv_borders=True):
    if ax is None:
        fig, ax = plt.subplots()
    inv_df = df.query('g_z_opt<=1.96 or power_opt>0.95')
    dep_df = df.query('g_z_opt>1.96 and power_opt<=0.95')
    ax.scatter(1-inv_df['power_opt'], inv_df['g_z_opt'], color=color)
    ax.scatter(1-dep_df['power_opt'], dep_df['g_z_opt'], color=color, alpha=0.5)
    if draw_c_inv_borders:
        ax.hlines(1.96, 0.05, 0.95, linestyles='dashed', alpha=0.5)
        ax.vlines(0.05, 1.96, np.max(df['g_z_opt']), linestyles='dashed',
                  alpha=0.5)
    ax.spines['left'].set_bounds(ax.dataLim.ymin, ax.dataLim.ymax)
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)
    ax.set_ylabel('g_z_opz')
    ax.set_xlabel('svd_power')


def contrast_amp_distr(exc_data, inh_data, ax=None, force_align_bins=True,
                       excit_col='darkorange', inhib_col='teal'):
    if ax is None:
        fig, ax = plt.subplots()
    tables = [exc_data, inh_data]
    colors = [excit_col, inhib_col]
    hists = []
    for i, (data, color) in enumerate(zip(tables, colors)):
        rmax, r0 = data.rmax.values, data.r0.values
        amp = (rmax - r0)*1000
        create_bins_with_first_hist = i == 0
        if create_bins_with_first_hist:
            if force_align_bins:
                bins = np.histogram_bin_edges(amp, bins='fd')
                bins = bins - bins[np.argmin(np.abs(bins))]
                h = np.histogram(amp, bins=bins)
            else:
                h = np.histogram(amp, bins='fd')
            split = sum(h[1] <= 0) - 1
            ax.hist(h[1][:split], h[1][:split+1], weights=h[0][:split], alpha=0.5, color=color)
            ax.hist(h[1][split:-1], h[1][split:], weights=h[0][split:], alpha=1, color=color)
            hists.append(h)
        else:
            h = np.histogram(amp, bins=hists[0][1])
            split = sum(h[1] <= 0) - 1
            ax.hist(h[1][:split], h[1][:split+1], weights=h[0][:split], alpha=0.5, color=color)
            ax.hist(h[1][split:-1], h[1][split:], weights=h[0][split:], alpha=1, color=color)
            hists.append(h)
    ax.vlines(0, 0, np.max(np.stack([h[0] for h in hists])), linestyle='dashed', alpha=0.5)
    ax.set_xlabel('Response amplitude (Hz)')
    ax.set_ylabel('N')
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)
    ax.spines['left'].set_bounds(ax.dataLim.ymin, ax.dataLim.ymax)


def marginal_tuning(datdf, pars, axs=None, tuning_model='wrap_gauss'):
    if axs is None:
        gridspec_kw = {'width_ratios': [2, 1], 'height_ratios': [1, 3]}
        fig, [[ax1, delax], [ax2, ax3]] = plt.subplots(2, 2, gridspec_kw=gridspec_kw,
                                                       tight_layout=True)
        fig.delaxes(delax)
    else:
        assert len(axs) == 3
        ax1, ax2, ax3 = axs
    # data
    datmat = datdf.pivot_table('p_spike_stim', 'grat_orientation', 'grat_contrast')
    firstrow = datmat.loc[:0, :]
    firstrow = firstrow.rename(index={0: 180.0})
    datmat = datmat.append(firstrow)
    ax2.imshow(datmat.values, origin='lower', aspect='auto')
    # smooth fit
    nfitpoints = 100
    smoothoris = np.linspace(0, 180, nfitpoints)
    smoothcons = contrast_log_scale(10, nfitpoints)
    smoothstim = np.meshgrid(smoothoris, smoothcons)
    fitmat = ori_con_2d(smoothstim, *pars, orifunc=tuning_model)
    # marginals
    # contrast marginal data points
    cons = np.unique(datdf.index.get_level_values('grat_contrast'))
    ax1.plot(np.arange(len(cons)), datmat.mean(axis=0)*1000, '.', color='k')
    # contrast marginal
    con_idcs = np.linspace(0, len(cons)-1, nfitpoints)
    ax1.plot(con_idcs, fitmat.mean(axis=0)*1000, color='k')
    # orientation marginal data points
    oris = np.unique(datmat.index.get_level_values('grat_orientation'))
    ax3.plot(datmat.mean(axis=1)*1000, np.arange(len(oris)), '.', color='k')
    # orientation marginal
    oriidcs = np.linspace(0, len(oris)-1, nfitpoints)
    ax3.plot(fitmat.mean(axis=1)*1000, oriidcs, color='k')
    # fix labels
    ax1.set_xticklabels([])
    ax1.set_ylabel('Response')
    ax3.set_yticklabels([])
    ax3.set_xlabel('Response')
    replace_imshow_ticklabels(ax2, cons, oris.astype(int))
    ax2.set_ylabel('Orientation ($^{\circ}$)')
    ax2.set_xlabel('Contrast')


def goodness_of_fit(v1_r2s, ax=None, excit_col='darkorange', inhib_col='teal'):
    if ax is None:
        fig, ax = plt.subplots()
    h = np.histogram(v1_r2s.query('unit_type=="excit"')['r2'].values, bins='fd')
    split = sum(h[1] <= 0.4)
    ax.hist(h[1][:split], h[1][:split + 1], weights=h[0][:split], alpha=0.5, color=excit_col)
    ax.hist(h[1][split:-1], h[1][split:], weights=h[0][split:], alpha=1, color=excit_col)

    h2 = np.histogram(v1_r2s.query('unit_type=="inhib"')['r2'].values, bins=h[1])
    ax.hist(h2[1][:split], h2[1][:split+1], weights=h2[0][:split], alpha=0.5, color=inhib_col)
    ax.hist(h2[1][split:-1], h2[1][split:], weights=h2[0][split:], alpha=1, color=inhib_col)

    ax.vlines(0.4, 0, np.max(np.stack([h[0], h2[0]])), linestyles='dashed', alpha=0.5)
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)
    ax.spines['left'].set_bounds(ax.dataLim.ymin, ax.dataLim.ymax)
    ax.set_xlabel('Goodness of fit (r²)')
    ax.set_ylabel('N')


def goodness_of_fit_lgn(v1_r2s, ax=None, color='k'):
    if ax is None:
        fig, ax = plt.subplots()
    h = np.histogram(v1_r2s['r2'].values, bins='fd')
    split = sum(h[1] <= 0.4)
    ax.hist(h[1][:split], h[1][:split + 1], weights=h[0][:split], alpha=0.5, color=color)
    ax.hist(h[1][split:-1], h[1][split:], weights=h[0][split:], alpha=1, color=color)

    ax.vlines(0.4, 0, np.max(np.stack([h[0]])), linestyles='dashed', alpha=0.5)
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)
    ax.spines['left'].set_bounds(ax.dataLim.ymin, ax.dataLim.ymax)
    ax.set_xlabel('Goodness of fit (r²)')
    ax.set_ylabel('N')


def population_tuning(tuning_widths, ax=None, color='darkorange'):
    if ax is None:
        fig, ax = plt.subplots()

    oris = np.linspace(0, 180, 100)
    labeloris = np.linspace(-90, 90, 100)
    for tuning_width in tuning_widths:
        ax.plot(labeloris, wrapped_gaussian(oris, sigma=tuning_width), color=color, alpha=0.2)
    ax.set_xlabel('Orientation ($^{\circ}$)')
    ax.set_ylabel('Response')
    ax.spines['left'].set_bounds(ax.dataLim.ymin, ax.dataLim.ymax)
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)

def population_tuning_von_mises(tuning_widths, ax=None, color='darkorange'):
    if ax is None:
        fig, ax = plt.subplots()

    oris = np.linspace(0, 180, 100)
    labeloris = np.linspace(-90, 90, 100)
    for tuning_width in tuning_widths:
        ax.plot(labeloris, von_mises(oris, kappa=tuning_width), color=color, alpha=0.2)
    ax.set_xlabel('Orientation ($^{\circ}$)')
    ax.set_ylabel('Response')
    ax.spines['left'].set_bounds(ax.dataLim.ymin, ax.dataLim.ymax)
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)

def population_contrast(pardf, ax=None, color='k', alpha=0.5):
    if ax is None:
        fig, ax = plt.subplots()
    cons = np.linspace(0, 1, 100)
    for uname, pars in pardf.iterrows():
        conres = hyper_ratio(cons, 1, pars['c50'], pars['n'], 0, pars['sup'])
        conres = conres/np.max(conres)
        ax.plot(cons, conres, color=color, alpha=alpha)
    ax.set_xlabel('Contrast')
    ax.set_ylabel('Response')
    ax.spines['left'].set_bounds(ax.dataLim.ymin, ax.dataLim.ymax)
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)


def osi_cumulative(osidf, ax=None, unit_type='excit', region='V1',
                   excit_col='darkorange', inhib_col='teal'):
    if ax is None:
        fig, ax = plt.subplots()
    # plot OSI distribution
    osis = osidf.query(f's_region=="{region}" and unit_type=="{unit_type}"').opt_osi
    colors = {'excit': excit_col, 'inhib': inhib_col}
    ecdf = ECDF(osis)
    ax.step(ecdf.x, ecdf.y, color=colors[unit_type])
    ax.set_xlabel('Orientation selectivity index')
    ax.set_ylabel('Fraction of units')
    ax.set_xlim(0 - ax.margins()[0], 1 + ax.margins()[1])
    ax.spines['bottom'].set_bounds(0, 1)
    ax.spines['left'].set_bounds(0, 1)


def osi_cumulative_lgn(osidf, ax=None, color='k'):
    if ax is None:
        fig, ax = plt.subplots()
    # plot OSI distribution
    osis = osidf.opt_osi
    ecdf = ECDF(osis)
    ax.step(ecdf.x, ecdf.y, color=color)
    ax.set_xlabel('Orientation selectivity index')
    ax.set_ylabel('Fraction of units')
    ax.set_xlim(0 - ax.margins()[0], 1 + ax.margins()[1])
    ax.spines['bottom'].set_bounds(0, 1)
    ax.spines['left'].set_bounds(0, 1)


def osi_hist(osidf, ax=None, unit_type='excit', region='V1', keys=None,
             excit_col='darkorange', inhib_col='teal'):
    if ax is None:
        fig, ax = plt.subplots()
    # plot OSI distribution
    osis = osidf.query(f's_region=="{region}" and unit_type=="{unit_type}"').opt_osi
    colors = {'excit': excit_col, 'inhib': inhib_col}
    ax.hist(osis, bins='fd', color=colors[unit_type], density=True, histtype='step')
    ax.set_xlabel('Orientation selectivity index')
    ax.set_ylabel('Nunits')


def osi_hist_lgn(osidf, ax=None, color='k'):
    if ax is None:
        fig, ax = plt.subplots()
    # plot OSI distribution
    osis = osidf.opt_osi
    ax.hist(osis, bins='fd', color=color, density=True, histtype='step')
    ax.set_xlabel('Orientation selectivity index')
    ax.set_ylabel('Nunits')


def contrast_cumulative(pardf, ax=None, color='darkorange'):
    if ax is None:
        fig, ax = plt.subplots()
    cons = np.linspace(0, 1, 100)
    c50s = []
    for uname, pars in pardf.iterrows():
        conres = hyper_ratio(cons, 1, pars['c50'], pars['n'], 0, pars['sup'])
        conres = conres/np.max(conres)
        c50idx = np.where(conres > 0.5)[0][0]
        c50s.append(cons[c50idx])
    ecdf = ECDF(c50s)
    ax.step(ecdf.x, ecdf.y, color=color)
    ax.set_xlabel('Contrast at half-height')
    ax.set_ylabel('Fraction of units')
    ax.set_xlim(0 - ax.margins()[0], 1 + ax.margins()[1])
    ax.spines['bottom'].set_bounds(0, 1)
    ax.spines['left'].set_bounds(0, 1)


def mean_population_contrast(pardf, ax=None, color='k'):
    if ax is None:
        fig, ax = plt.subplots()
    cons = np.linspace(0, 1, 100)
    res = []
    for uname, pars in pardf.iterrows():
        conres = hyper_ratio(cons, pars['rmax'], pars['c50'], pars['n'], pars['r0'],
                             pars['sup'])
        res.append(conres)
    resmat = np.stack(res)*1000
    pop_conres = resmat.mean(axis=0)
    res_sem = sem(resmat, axis=0)
    ax.plot(cons, pop_conres, color=color)
    ax.fill_between(cons, pop_conres + res_sem, pop_conres - res_sem, color=color,
                    alpha=0.5)
    ax.set_xlabel('Contrast')
    ax.set_ylabel('Response')
    ax.spines['left'].set_bounds(0, ax.dataLim.ymax)
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)
    _, ymax = ax.get_ylim()
    ax.set_ylim([0, ymax])


def mean_norm_population_tuning(tuning_widths, ax=None, color='darkorange'):
    if ax is None:
        fig, ax = plt.subplots()
    oris = np.linspace(0, 180, 100)
    labeloris = np.linspace(-90, 90, 100)
    res = []
    for tw in tuning_widths:
        res.append(wrapped_gaussian(oris, sigma=tw))
    resmat = np.stack(res)
    mean_res = resmat.mean(axis=0)
    sem_res = sem(resmat, axis=0)
    ax.plot(labeloris, mean_res, color=color)
    ax.fill_between(labeloris, mean_res + sem_res, mean_res - sem_res,
                    color=color, alpha=0.5)
    ax.set_xlabel('Orientation ($^{\circ}$)')
    ax.set_ylabel('Response')
    ax.spines['left'].set_bounds(0, ax.dataLim.ymax)
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)
    _, ymax = ax.get_ylim()
    ax.set_ylim([0, ymax])

def mean_norm_population_tuning_von_mises(tuning_widths, ax=None, color='darkorange'):
    if ax is None:
        fig, ax = plt.subplots()
    oris = np.linspace(0, 180, 100)
    labeloris = np.linspace(-90, 90, 100)
    res = []
    for tw in tuning_widths:
        res.append(von_mises(oris, kappa=tw))
    resmat = np.stack(res)
    mean_res = resmat.mean(axis=0)
    sem_res = sem(resmat, axis=0)
    ax.plot(labeloris, mean_res, color=color)
    ax.fill_between(labeloris, mean_res + sem_res, mean_res - sem_res,
                    color=color, alpha=0.5)
    ax.set_xlabel('Orientation ($^{\circ}$)')
    ax.set_ylabel('Response')
    ax.spines['left'].set_bounds(0, ax.dataLim.ymax)
    ax.spines['bottom'].set_bounds(ax.dataLim.xmin, ax.dataLim.xmax)
    _, ymax = ax.get_ylim()
    ax.set_ylim([0, ymax])