"""Figure 1S3 and 3S1 unit classification plots, use run -i fig1S33S1.py"""

index = pd.Index(mvigrtmsustrs, name='msu')
columns = ['mvi_meanrate', 'mvi_meanburstratio', 'grt_meanrate', 'grt_meanburstratio',
           'sbc', 'depth', 'normdepth', 'dsi', 'rfdist',
           'mvi_meanrate_raw', 'mvi_meanburstratio_raw',
           'grt_meanrate_raw', 'grt_meanburstratio_raw']
fig1S33S1 = pd.DataFrame(index=index, columns=columns)

from scipy.stats import gmean

import matplotlib as mpl
from matplotlib.patches import Rectangle

TRANSTHRESH = 0.2
ONOFFTHRESH = 0.2


# stripplot FMI by SbC/non-SbC, mvi and grt, meanrate and meanburstratio:
np.random.seed(0) # to get identical horizontal jitter in strip plots on every run
figsize = DEFAULTFIGURESIZE
for stimtype in STIMTYPES:
    stimtypelabel = stimtype2axislabel[stimtype]
    for measure in ['meanrate', 'meanburstratio']:
        axislabel = measure2axislabel[measure]
        axislabel = short2longaxislabel.get(axislabel, axislabel)
        if axislabel.islower():
            axislabel = axislabel.capitalize()
        fmis, sbcs = [], []
        colname = '_'.join([stimtype, measure]) # e.g. 'mvi_meanrate'
        for msustr in mvigrtmsustrs:
            # save sbc to fig1S33S1 df, regardless of FMI val, will be overwritten multiple
            # times with the same value, but that's OK:
            sbc = celltype.loc[msustr]['sbc'] # bool
            fig1S33S1.loc[msustr]['sbc'] = sbc
            fmi = maxFMI.loc[msustr, 'none', stimtype][measure] # ignore run condition for FMI
            if np.isnan(fmi):
                continue
            fmis.append(fmi)
            sbcs.append(sbc)
            fig1S33S1.loc[msustr][colname] = fmi
        fmis = np.asarray(fmis)
        sbcs = np.asarray(sbcs)
        sbcfmis = fmis[sbcs == True]
        nonsbcfmis = fmis[sbcs == False]
        f, a = plt.subplots(figsize=figsize)
        wintitle('FMI SbC %s %s strip' % (stimtypelabel, measure))
        # plot y=0 line:
        a.axhline(y=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
        data = pd.DataFrame.from_dict({'SbC':sbcfmis, 'Non-SbC':nonsbcfmis},
                                      orient='index').transpose()
        sns.stripplot(ax=a, data=data, clip_on=False, marker='.',
                      color='None', edgecolor='black', size=np.sqrt(50))
        # get fname of appropriate LMM .cvs file:
        fname = None # sanity check: clear from previous loop
        if stimtype == 'mvi':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_1S3a_pred_means.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_1S3f_pred_means.csv')
        elif stimtype == 'grt':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_3S1a_pred_means.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_3S1f_pred_means.csv')
        # fetch LMM means from .csv:
        df = pd.read_csv(fname)
        meannonsbcfmi = df['non_sbc'][0]
        meansbcfmi = df['sbc'][0]
        # plot mean with short horizontal lines:
        a.plot([-0.25, 0.25], [meansbcfmi, meansbcfmi], '-', lw=2, c='red', zorder=np.inf)
        a.plot([0.75, 1.25], [meannonsbcfmi, meannonsbcfmi], '-', lw=2, c='red', zorder=np.inf)
        a.set_ylabel('%s FMI' % axislabel)
        a.set_ylim(-1, 1)
        a.set_yticks([-1, 0, 1])
        a.tick_params(bottom=False)
        a.spines['bottom'].set_position(('outward', 5))
        a.spines['bottom'].set_visible(False)

'''
# stripplot FMI by shell/core, mvi and grt, meanrate and meanburstratio:
np.random.seed(0) # to get identical horizontal jitter in strip plots on every run
figsize = DEFAULTFIGURESIZE
for stimtype in STIMTYPES:
    stimtypelabel = stimtype2axislabel[stimtype]
    for measure in ['meanrate', 'meanburstratio']:
        axislabel = measure2axislabel[measure]
        axislabel = short2longaxislabel.get(axislabel, axislabel)
        if axislabel.islower():
            axislabel = axislabel.capitalize()
        fmis, layers = [], []
        colname = '_'.join([stimtype, measure]) # e.g. 'mvi_meanrate'
        for msustr in mvigrtmsustrs:
            # save layer to fig1S33S1 df, regardless of FMI val, will be overwritten multiple
            # times with the same value, but that's OK:
            layer = celltype.loc[msustr]['layer'] # str
            fig1S33S1.loc[msustr]['layer'] = layer
            fmi = maxFMI.loc[msustr, 'none', stimtype][measure] # ignore run condition for FMI
            if np.isnan(fmi):
                continue
            fmis.append(fmi)
            layers.append(layer)
            fig1S33S1.loc[msustr][colname] = fmi # might overwrite w/ identical values
        fmis = np.asarray(fmis)
        layers = np.asarray(layers)
        shellfmis = fmis[layers == 'shell']
        corefmis = fmis[layers == 'core']
        f, a = plt.subplots(figsize=figsize)
        wintitle('FMI layer %s %s strip' % (stimtypelabel, measure))
        # plot y=0 line:
        a.axhline(y=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
        data = pd.DataFrame.from_dict({'Shell':shellfmis, 'Core':corefmis},
                                      orient='index').transpose()
        sns.stripplot(ax=a, data=data, clip_on=False, marker='.',
                      color='None', edgecolor='black', size=np.sqrt(50))
        # get fname of appropriate LMM .cvs file:
        fname = None # sanity check: clear from previous loop
        if stimtype == 'mvi':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_1S3b_pred_means.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_1S3g_pred_means.csv')
        elif stimtype == 'grt':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_3S1b_pred_means.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_3S1g_pred_means.csv')
        # fetch LMM means from .csv:
        df = pd.read_csv(fname)
        meancorefmi = df['core'][0]
        meanshellfmi = df['shell'][0]
        # plot mean with short horizontal lines:
        a.plot([-0.25, 0.25], [meanshellfmi, meanshellfmi], '-', lw=2, c='red', zorder=np.inf)
        a.plot([0.75, 1.25], [meancorefmi, meancorefmi], '-', lw=2, c='red', zorder=np.inf)
        a.set_ylabel('%s FMI' % axislabel)
        a.set_ylim(-1, 1)
        a.set_yticks([-1, 0, 1])
        a.tick_params(bottom=False)
        a.spines['bottom'].set_position(('outward', 5))
        a.spines['bottom'].set_visible(False)
'''


# scatter plot FMI vs. depth by mvi and grt, meanrate and meanburstratio:
figsize = DEFAULTFIGURESIZE
for stimtype in STIMTYPES:
    stimtypelabel = stimtype2axislabel[stimtype]
    for measure in ['meanrate', 'meanburstratio']:
        axislabel = measure2axislabel[measure]
        axislabel = short2longaxislabel.get(axislabel, axislabel)
        if axislabel.islower():
            axislabel = axislabel.capitalize()
        fmis, depths = [], []
        colname = '_'.join([stimtype, measure]) # e.g. 'mvi_meanrate'
        for msustr in mvigrtmsustrs:
            # save depth to fig1S33S1 df, regardless of FMI val, will be overwritten multiple
            # times with the same value, but that's OK:
            depth = celltype.loc[msustr]['depth'] # float
            fig1S33S1.loc[msustr]['depth'] = depth
            fmi = maxFMI.loc[msustr, 'none', stimtype][measure] # ignore run condition for FMI
            if np.isnan(fmi):
                continue
            fmis.append(fmi)
            depths.append(depth)
            fig1S33S1.loc[msustr][colname] = fmi # might overwrite w/ identical values
        fmis = np.asarray(fmis)
        depths = np.asarray(depths)
        f, a = plt.subplots(figsize=figsize)
        wintitle('FMI depth %s %s' % (stimtypelabel, measure))
        # plot y=0 line:
        a.axhline(y=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
        a.scatter(depths, fmis, clip_on=False, marker='.', c='None', edgecolor='k', s=DEFSZ)
        # get fname of appropriate LMM .cvs file:
        fname = None # sanity check: clear from previous loop
        if stimtype == 'mvi':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_1S3b_coefs.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_1S3g_coefs.csv')
        elif stimtype == 'grt':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_3S1b_coefs.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_3S1g_coefs.csv')
        # fetch LMM linregress fit params from .csv:
        df = pd.read_csv(fname)
        mm = df['slope'][0]
        b = df['intercept'][0]
        x = np.array([np.nanmin(depths), np.nanmax(depths)])
        y = mm * x + b
        a.plot(x, y, '-', color='red') # plot linregress fit
        a.set_xlabel('Depth ($\mathregular{\mu}$m)')
        a.set_ylabel('%s FMI' % axislabel)
        a.set_xlim(0, 500)
        a.set_ylim(-1, 1)
        a.set_yticks([-1, 0, 1])
        a.spines['left'].set_position(('outward', 4))
        a.spines['bottom'].set_position(('outward', 4))


# scatter plot FMI vs. normalized depth by mvi and grt, meanrate and meanburstratio:
figsize = DEFAULTFIGURESIZE
for stimtype in STIMTYPES:
    stimtypelabel = stimtype2axislabel[stimtype]
    for measure in ['meanrate', 'meanburstratio']:
        axislabel = measure2axislabel[measure]
        axislabel = short2longaxislabel.get(axislabel, axislabel)
        if axislabel.islower():
            axislabel = axislabel.capitalize()
        fmis, normdepths = [], []
        colname = '_'.join([stimtype, measure]) # e.g. 'mvi_meanrate'
        for msustr in mvigrtmsustrs:
            # save normdepth to fig1S33S1 df, regardless of FMI val, will be overwritten
            # multiple times with the same value, but that's OK:
            normdepth = celltype.loc[msustr]['normdepth'] # float
            fig1S33S1.loc[msustr]['normdepth'] = normdepth
            fmi = maxFMI.loc[msustr, 'none', stimtype][measure] # ignore run condition for FMI
            if np.isnan(fmi):
                continue
            fmis.append(fmi)
            normdepths.append(normdepth)
            fig1S33S1.loc[msustr][colname] = fmi # might overwrite w/ identical values
        fmis = np.asarray(fmis)
        normdepths = np.asarray(normdepths)
        f, a = plt.subplots(figsize=figsize)
        wintitle('FMI normdepth %s %s' % (stimtypelabel, measure))
        # plot y=0 line:
        a.axhline(y=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
        a.scatter(normdepths, fmis, clip_on=False, marker='.', c='None', edgecolor='k', s=DEFSZ)
        '''
        # get fname of appropriate LMM .cvs file:
        fname = None # sanity check: clear from previous loop
        if stimtype == 'mvi':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_1S3b_coefs.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_1S3g_coefs.csv')
        elif stimtype == 'grt':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_3S1b_coefs.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_3S1g_coefs.csv')
        # fetch LMM linregress fit params from .csv:
        df = pd.read_csv(fname)
        mm = df['slope'][0]
        b = df['intercept'][0]
        x = np.array([np.nanmin(normdepths), np.nanmax(normdepths)])
        y = mm * x + b
        a.plot(x, y, '-', color='red') # plot linregress fit
        '''
        a.set_xlabel('Normalized depth')
        a.set_ylabel('%s FMI' % axislabel)
        #a.set_xlim(0, 2)
        a.set_ylim(-1, 1)
        a.set_yticks([-1, 0, 1])
        a.spines['left'].set_position(('outward', 4))
        a.spines['bottom'].set_position(('outward', 4))


# scatter plot FMI vs. DSI by mvi and grt, meanrate and meanburstratio:
figsize = DEFAULTFIGURESIZE
for stimtype in STIMTYPES:
    stimtypelabel = stimtype2axislabel[stimtype]
    for measure in ['meanrate', 'meanburstratio']:
        axislabel = measure2axislabel[measure]
        axislabel = short2longaxislabel.get(axislabel, axislabel)
        if axislabel.islower():
            axislabel = axislabel.capitalize()
        fmis, dsis = [], []
        colname = '_'.join([stimtype, measure]) # e.g. 'mvi_meanrate'
        for msustr in mvigrtmsustrs:
            # save DSI to fig1S33S1 df, regardless of FMI val, will be overwritten multiple
            # times with the same value, but that's OK:
            dsi = celltype.loc[msustr]['dsi'] # float
            fig1S33S1.loc[msustr]['dsi'] = dsi
            fmi = maxFMI.loc[msustr, 'none', stimtype][measure] # ignore run condition for FMI
            if np.isnan(fmi):
                continue
            fmis.append(fmi)
            dsis.append(dsi)
            fig1S33S1.loc[msustr][colname] = fmi # might overwrite w/ identical values
        fmis = np.asarray(fmis)
        dsis = np.asarray(dsis)
        f, a = plt.subplots(figsize=figsize)
        wintitle('FMI DSI %s %s' % (stimtypelabel, measure))
        # plot y=0 line:
        a.axhline(y=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
        a.scatter(dsis, fmis, clip_on=False, marker='.', c='None', edgecolor='k', s=DEFSZ)
        # get fname of appropriate LMM .cvs file:
        fname = None # sanity check: clear from previous loop
        if stimtype == 'mvi':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_1S3c_coefs.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_1S3h_coefs.csv')
        elif stimtype == 'grt':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_3S1c_coefs.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_3S1h_coefs.csv')
        # fetch LMM linregress fit params from .csv:
        df = pd.read_csv(fname)
        mm = df['slope'][0]
        b = df['intercept'][0]
        x = np.array([np.nanmin(dsis), np.nanmax(dsis)])
        y = mm * x + b
        a.plot(x, y, '-', color='red') # plot linregress fit
        a.set_xlabel('DSI')
        a.set_ylabel('%s FMI' % axislabel)
        a.set_xlim(0, 1)
        a.set_ylim(-1, 1)
        a.set_yticks([-1, 0, 1])
        a.spines['left'].set_position(('outward', 4))
        a.spines['bottom'].set_position(('outward', 4))


# scatter plot FMI vs. distance of MUA envl RF from screen center, by mvi and grt,
# meanrate and meanburstratio:
figsize = DEFAULTFIGURESIZE
for stimtype in STIMTYPES:
    stimtypelabel = stimtype2axislabel[stimtype]
    for measure in ['meanrate', 'meanburstratio']:
        axislabel = measure2axislabel[measure]
        axislabel = short2longaxislabel.get(axislabel, axislabel)
        if axislabel.islower():
            axislabel = axislabel.capitalize()
        fmis, ds = [], []
        colname = '_'.join([stimtype, measure]) # e.g. 'mvi_meanrate'
        for msustr in mvigrtmsustrs:
            # save rfdist to fig1S33S1 df, regardless of FMI val, will be overwritten multiple
            # times with the same value, but that's OK:
            x0, y0 = cellscreenpos.loc[msustr]
            d = np.sqrt(x0**2 + y0**2) # distance from screen center, deg
            fig1S33S1.loc[msustr]['rfdist'] = d
            fmi = maxFMI.loc[msustr, 'none', stimtype][measure] # ignore run condition for FMI
            if np.isnan(fmi):
                continue
            fmis.append(fmi)
            ds.append(d)
            fig1S33S1.loc[msustr][colname] = fmi # might overwrite w/ identical values
        fmis = np.asarray(fmis)
        ds = np.asarray(ds)
        f, a = plt.subplots(figsize=figsize)
        wintitle('FMI rfdist %s %s' % (stimtypelabel, measure))
        # plot y=0 line:
        a.axhline(y=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
        a.scatter(ds, fmis, clip_on=False, marker='.', c='None', edgecolor='k', s=DEFSZ)
        # get fname of appropriate LMM .cvs file:
        fname = None # sanity check: clear from previous loop
        if stimtype == 'mvi':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_1S3d_coefs.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_1S3i_coefs.csv')
        elif stimtype == 'grt':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_3S1d_coefs.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_3S1i_coefs.csv')
        # fetch LMM linregress fit params from .csv:
        df = pd.read_csv(fname)
        mm = df['slope'][0]
        b = df['intercept'][0]
        x = np.array([np.nanmin(ds), np.nanmax(ds)])
        y = mm * x + b
        a.plot(x, y, '-', color='red') # plot linregress fit
        a.set_xlabel('RF dist. from center ($\degree$)')
        a.set_ylabel('%s FMI' % axislabel)
        a.set_xlim(0, 40)
        a.set_ylim(-1, 1)
        a.set_yticks([-1, 0, 1])
        a.spines['left'].set_position(('outward', 4))
        a.spines['bottom'].set_position(('outward', 4))


# scatter plot FMI vs raw measure, for rate and burst ratio during control condition,
# by mvi and grt:
stimtype2resp = {'mvi':mviresp, 'grt':grtresp}
figsize = DEFAULTFIGURESIZE
for stimtype in STIMTYPES:
    stimtypelabel = stimtype2axislabel[stimtype]
    resp = stimtype2resp[stimtype]
    if stimtype == 'mvi':
        resp = resp.xs('nat', level='kind') # dereference movie 'kind' index level
    for measure in ['meanrate', 'meanburstratio']:
        axislabel = measure2axislabel[measure]
        axislabel = short2longaxislabel.get(axislabel, axislabel)
        axisunits = measure2axisunits.get(measure, '')
        if axislabel.islower():
            axislabel = axislabel.capitalize()
        fmis, msrs = [], []
        fmicolname = '_'.join([stimtype, measure]) # e.g. 'mvi_meanrate'
        msrcolname = '_'.join([stimtype, measure, 'raw']) # e.g. 'mvi_meanrate_raw'
        for msustr in mvigrtmsustrs:
            mseustr, fmi = maxFMI.loc[msustr, 'none', stimtype][['mseu', measure]]
            if pd.isna(mseustr) or pd.isna(fmi):
                continue
            msr = resp.loc[mseustr, 'none', False][measure]
            fmis.append(fmi)
            msrs.append(msr)
            fig1S33S1.loc[msustr][fmicolname] = fmi # might overwrite w/ identical values
            fig1S33S1.loc[msustr][msrcolname] = msr
        fmis, msrs = np.asarray(fmis), np.asarray(msrs)
        ## scatter plot FMI vs raw measure:
        f, a = plt.subplots(figsize=figsize)
        wintitle('FMI raw %s %s' % (stimtypelabel, measure))
        # plot y=0 line:
        a.axhline(y=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
        a.scatter(msrs, fmis, clip_on=False, marker='.', c='None', edgecolor='k', s=DEFSZ)
        # get fname of appropriate LMM .cvs file:
        fname = None # sanity check: clear from previous loop
        if stimtype == 'mvi':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_1S3e_coefs.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_1S3j_coefs.csv')
        elif stimtype == 'grt':
            if measure == 'meanrate':
                fname = os.path.join('stats', 'figure_3S1e_coefs.csv')
            elif measure == 'meanburstratio':
                fname = os.path.join('stats', 'figure_3S1j_coefs.csv')
        # fetch LMM linregress fit params from .csv:
        df = pd.read_csv(fname)
        mm = df['slope'][0]
        b = df['intercept'][0]
        x = np.array([np.nanmin(msrs), np.nanmax(msrs)])
        y = mm * x + b
        a.plot(x, y, '-', color='red') # plot linregress fit
        a.set_xlabel('%s' % axislabel+axisunits)
        a.set_ylabel('%s FMI' % axislabel)
        a.set_xlim(xmin=0)
        a.set_ylim(-1, 1)
        #a.set_xticks(ticks)
        a.set_yticks([-1, 0, 1])
        a.spines['left'].set_position(('outward', 4))
        a.spines['bottom'].set_position(('outward', 4))


#### old opto scatter plots:
'''
# scatter plot grating meanrates, coloured by sbc:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -1.2, 2
logticks = np.array([-1, 0, 1, 2])
#log0min = logmin + 0.05
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto meanrate grating %s sbc' % st8)
    rons, roffs, sbcs = [], [], []
    for mseustr in grtmseustrs:
        meanrate = grtresp.loc[mseustr, st8]['meanrate']
        msustr = mseustr2msustr(mseustr)
        sbc = celltype.loc[msustr]['sbc']
        if np.isnan(sbc):
            sbc = None # works properly as sbcclrs dict key
        if meanrate.isna().any(): # missing one or both meanrates
            #print('%s: missing one or both opto conditions, skipping' % mseustr)
            continue
        rons.append(meanrate[True])
        roffs.append(meanrate[False])
        sbcs.append(sbc)
        #fig3.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
    rons = np.asarray(rons)
    roffs = np.asarray(roffs)
    sbcs = np.asarray(sbcs)
    sbcclrs = [ {True:'r', False:green, None:'black'}[sbc] for sbc in sbcs ]
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(rons, roffs, marker='.', c='None', edgecolor=sbcclrs, s=DEFSZ)
    # plot mean of points grouped by sbc:
    sbclogmean = gmean(rons[sbcs == True]), gmean(roffs[sbcs == True])
    notsbclogmean = gmean(rons[sbcs == False]), gmean(roffs[sbcs == False])
    a.scatter(sbclogmean[0], sbclogmean[1], marker='+', c='r', s=100)
    a.scatter(notsbclogmean[0], notsbclogmean[1], marker='+', c=green, s=100)
    a.set_ylabel('Feedback FR (spk/s)')
    a.set_xlabel('Suppression FR (spk/s)')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    axes_disable_scientific(a)
    a.minorticks_off()
    a.set_aspect('equal')
    #t, p = ttest_rel(rons, roffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
    patch = Rectangle((0, 0), 0, 0)
    l = a.legend([patch, patch], ['SbC', 'Non-SbC'], loc='upper left',
                 handlelength=0, handletextpad=0, frameon=False)
    for t, c in zip(l.get_texts(), ['r', green]):
        t.set_color(c)


# scatter plot movie meanrates, coloured by sbc:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -1.2, 2
logticks = np.array([-1, 0, 1, 2])
#log0min = logmin + 0.05
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto meanrate movie nat %s sbc' % st8)
    rons, roffs, sbcs = [], [], []
    for mseustr in mvimseustrs:
        meanrate = mviresp.loc[mseustr, 'nat', st8]['meanrate']
        msustr = mseustr2msustr(mseustr)
        sbc = celltype.loc[msustr]['sbc']
        if np.isnan(sbc):
            sbc = None # works properly as sbcclrs dict key
        if meanrate.isna().any(): # missing one or both meanrates
            #print('%s: missing one or both opto conditions, skipping' % mseustr)
            continue
        rons.append(meanrate[True])
        roffs.append(meanrate[False])
        sbcs.append(sbc)
        #fig3.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
    rons = np.asarray(rons)
    roffs = np.asarray(roffs)
    sbcs = np.asarray(sbcs)
    sbcclrs = [ {True:'r', False:green, None:'black'}[sbc] for sbc in sbcs ]
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(rons, roffs, marker='.', c='None', edgecolor=sbcclrs, s=DEFSZ)
    # plot mean of points grouped by sbc:
    sbclogmean = gmean(rons[sbcs == True]), gmean(roffs[sbcs == True])
    notsbclogmean = gmean(rons[sbcs == False]), gmean(roffs[sbcs == False])
    a.scatter(sbclogmean[0], sbclogmean[1], marker='+', c='r', s=100)
    a.scatter(notsbclogmean[0], notsbclogmean[1], marker='+', c=green, s=100)
    a.set_ylabel('Feedback FR (spk/s)')
    a.set_xlabel('Suppression FR (spk/s)')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    axes_disable_scientific(a)
    a.minorticks_off()
    a.set_aspect('equal')
    #t, p = ttest_rel(rons, roffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
    patch = Rectangle((0, 0), 0, 0)
    l = a.legend([patch, patch], ['SbC', 'Non-SbC'], loc='upper left',
                 handlelength=0, handletextpad=0, frameon=False)
    for t, c in zip(l.get_texts(), ['r', green]):
        t.set_color(c)


# scatter plot grating burst ratio, coloured by sbc:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -3.55, 0
logticks = np.array([-3, -2, -1, 0])
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto burst ratio grating %s sbc' % st8)
    brons, broffs, sbcs = [], [], []
    for mseustr in grtmseustrs:
        br = grtresp.loc[mseustr, st8]['meanburstratio']
        msustr = mseustr2msustr(mseustr)
        sbc = celltype.loc[msustr]['sbc']
        if np.isnan(sbc):
            sbc = None # works properly as sbcclrs dict key
        if br.isna().any(): # missing for at least one opto condition
            continue
        brons.append(br[True])
        broffs.append(br[False])
        sbcs.append(sbc)
        #fig3.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
    brons, broffs, sbcs = np.asarray(brons), np.asarray(broffs), np.asarray(sbcs)
    sbcclrs = [ {True:'r', False:green, None:'black'}[sbc] for sbc in sbcs ]
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(brons, broffs, marker='.', c='None', edgecolor=sbcclrs, s=DEFSZ)
    # plot mean of points grouped by sbc:
    # filter out 0 value burst ratios before calculating geometric mean:
    sbcbrons = brons[sbcs == True][brons[sbcs == True] > 0]
    sbcbroffs = broffs[sbcs == True][broffs[sbcs == True] > 0]
    notsbcbrons = brons[sbcs == False][brons[sbcs == False] > 0]
    notsbcbroffs = broffs[sbcs == False][broffs[sbcs == False] > 0]
    sbclogmean = gmean(sbcbrons), gmean(sbcbroffs)
    notsbclogmean = gmean(notsbcbrons), gmean(notsbcbroffs)
    a.scatter(sbclogmean[0], sbclogmean[1], marker='+', c='r', s=100)
    a.scatter(notsbclogmean[0], notsbclogmean[1], marker='+', c=green, s=100)
    a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
    a.set_ylabel('Feedback BR')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    a.minorticks_off()
    axes_disable_scientific(a)
    a.set_aspect('equal')
    #t, p = ttest_rel(brons, broffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
    patch = Rectangle((0, 0), 0, 0)
    l = a.legend([patch, patch], ['SbC', 'Non-SbC'], loc='upper left',
                 handlelength=0, handletextpad=0, frameon=False)
    for t, c in zip(l.get_texts(), ['r', green]):
        t.set_color(c)


# scatter plot movie burst ratio, coloured by sbc:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -3.55, 0
logticks = np.array([-3, -2, -1, 0])
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto burst ratio movie nat %s sbc' % st8)
    brons, broffs, sbcs = [], [], []
    for mseustr in mvimseustrs:
        br = mviresp.loc[mseustr, 'nat', st8]['meanburstratio']
        msustr = mseustr2msustr(mseustr)
        sbc = celltype.loc[msustr]['sbc']
        if np.isnan(sbc):
            sbc = None
        if br.isna().any(): # missing for at least one opto condition
            continue
        brons.append(br[True])
        broffs.append(br[False])
        sbcs.append(sbc)
        #fig3.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
    brons, broffs, sbcs = np.asarray(brons), np.asarray(broffs), np.asarray(sbcs)
    sbcclrs = [ {True:'r', False:green, None:'black'}[sbc] for sbc in sbcs ]
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(brons, broffs, marker='.', c='None', edgecolor=sbcclrs, s=DEFSZ)
    # plot mean of points grouped by sbc:
    # filter out 0 value burst ratios before calculating geometric mean:
    sbcbrons = brons[sbcs == True][brons[sbcs == True] > 0]
    sbcbroffs = broffs[sbcs == True][broffs[sbcs == True] > 0]
    notsbcbrons = brons[sbcs == False][brons[sbcs == False] > 0]
    notsbcbroffs = broffs[sbcs == False][broffs[sbcs == False] > 0]
    sbclogmean = gmean(sbcbrons), gmean(sbcbroffs)
    notsbclogmean = gmean(notsbcbrons), gmean(notsbcbroffs)
    a.scatter(sbclogmean[0], sbclogmean[1], marker='+', c='r', s=100)
    a.scatter(notsbclogmean[0], notsbclogmean[1], marker='+', c=green, s=100)
    a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
    a.set_ylabel('Feedback BR')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    a.minorticks_off()
    axes_disable_scientific(a)
    a.set_aspect('equal')
    #t, p = ttest_rel(brons, broffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
    patch = Rectangle((0, 0), 0, 0)
    l = a.legend([patch, patch], ['SbC', 'Non-SbC'], loc='upper left',
                 handlelength=0, handletextpad=0, frameon=False)
    for t, c in zip(l.get_texts(), ['r', green]):
        t.set_color(c)


# scatter plot grating meanrates, coloured by DSI:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -1.2, 2
logticks = np.array([-1, 0, 1, 2])
#log0min = logmin + 0.05
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto meanrate grating %s dsi' % st8)
    rons, roffs, dsis = [], [], []
    for mseustr in grtmseustrs:
        meanrate = grtresp.loc[mseustr, st8]['meanrate']
        msustr = mseustr2msustr(mseustr)
        dsi = celltype.loc[msustr]['dsi']
        if meanrate.isna().any(): # missing one or both meanrates
            #print('%s: missing one or both opto conditions, skipping' % mseustr)
            continue
        rons.append(meanrate[True])
        roffs.append(meanrate[False])
        dsis.append(dsi)
        #fig3.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
    rons = np.asarray(rons)
    roffs = np.asarray(roffs)
    dsis = np.asarray(dsis)
    cmap = plt.get_cmap() # use default, normalized from 0 to 1
    dsiclrs = []
    for dsi in dsis:
        if np.isnan(dsi):
            dsiclrs.append('gray')
        else:
            dsiclrs.append(cmap(dsi))
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(rons, roffs, marker='.', c='None', edgecolor=dsiclrs, s=DEFSZ)
    a.set_ylabel('Feedback FR (spk/s)')
    a.set_xlabel('Suppression FR (spk/s)')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    axes_disable_scientific(a)
    a.minorticks_off()
    a.set_aspect('equal')
    #t, p = ttest_rel(rons, roffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
    #scalarmappable = mpl.cm.ScalarMappable(norm=None, cmap=cmap)
    #f.colorbar(scalarmappable, label='DSI')


# scatter plot movie meanrates, coloured by dsi:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -1.2, 2
logticks = np.array([-1, 0, 1, 2])
#log0min = logmin + 0.05
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto meanrate movie nat %s dsi' % st8)
    rons, roffs, dsis = [], [], []
    for mseustr in mvimseustrs:
        meanrate = mviresp.loc[mseustr, 'nat', st8]['meanrate']
        msustr = mseustr2msustr(mseustr)
        dsi = celltype.loc[msustr]['dsi']
        if meanrate.isna().any(): # missing one or both meanrates
            #print('%s: missing one or both opto conditions, skipping' % mseustr)
            continue
        rons.append(meanrate[True])
        roffs.append(meanrate[False])
        dsis.append(dsi)
        #fig3.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
    rons = np.asarray(rons)
    roffs = np.asarray(roffs)
    dsis = np.asarray(dsis)
    cmap = plt.get_cmap() # use default, normalized from 0 to 1
    dsiclrs = []
    for dsi in dsis:
        if np.isnan(dsi):
            dsiclrs.append('gray')
        else:
            dsiclrs.append(cmap(dsi))
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(rons, roffs, marker='.', c='None', edgecolor=dsiclrs, s=DEFSZ)
    a.set_ylabel('Feedback FR (spk/s)')
    a.set_xlabel('Suppression FR (spk/s)')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    axes_disable_scientific(a)
    a.minorticks_off()
    a.set_aspect('equal')
    #t, p = ttest_rel(rons, roffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
    #scalarmappable = mpl.cm.ScalarMappable(norm=None, cmap=cmap)
    #f.colorbar(scalarmappable, label='DSI')


# scatter plot grating burst ratio, coloured by dsi:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -3.55, 0
logticks = np.array([-3, -2, -1, 0])
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto burst ratio grating %s dsi' % st8)
    brons, broffs, dsis = [], [], []
    for mseustr in grtmseustrs:
        br = grtresp.loc[mseustr, st8]['meanburstratio']
        msustr = mseustr2msustr(mseustr)
        dsi = celltype.loc[msustr]['dsi']
        if br.isna().any(): # missing for at least one opto condition
            continue
        brons.append(br[True])
        broffs.append(br[False])
        dsis.append(dsi)
        #fig3.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
    brons, broffs, dsis = np.asarray(brons), np.asarray(broffs), np.asarray(dsis)
    cmap = plt.get_cmap() # use default, normalized from 0 to 1
    dsiclrs = []
    for dsi in dsis:
        if np.isnan(dsi):
            dsiclrs.append('gray')
        else:
            dsiclrs.append(cmap(dsi))
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(brons, broffs, marker='.', c='None', edgecolor=dsiclrs, s=DEFSZ)
    a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
    a.set_ylabel('Feedback BR')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    a.minorticks_off()
    axes_disable_scientific(a)
    a.set_aspect('equal')
    #t, p = ttest_rel(brons, broffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
    #scalarmappable = mpl.cm.ScalarMappable(norm=None, cmap=cmap)
    #f.colorbar(scalarmappable, label='DSI')


# scatter plot movie burst ratio, coloured by dsi:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -3.55, 0
logticks = np.array([-3, -2, -1, 0])
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto burst ratio movie nat %s dsi' % st8)
    brons, broffs, dsis = [], [], []
    for mseustr in mvimseustrs:
        br = mviresp.loc[mseustr, 'nat', st8]['meanburstratio']
        msustr = mseustr2msustr(mseustr)
        dsi = celltype.loc[msustr]['dsi']
        if br.isna().any(): # missing for at least one opto condition
            continue
        brons.append(br[True])
        broffs.append(br[False])
        dsis.append(dsi)
        #fig3.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
    brons, broffs, dsis = np.asarray(brons), np.asarray(broffs), np.asarray(dsis)
    cmap = plt.get_cmap() # use default, normalized from 0 to 1
    dsiclrs = []
    for dsi in dsis:
        if np.isnan(dsi):
            dsiclrs.append('gray')
        else:
            dsiclrs.append(cmap(dsi))
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(brons, broffs, marker='.', c='None', edgecolor=dsiclrs, s=DEFSZ)
    a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
    a.set_ylabel('Feedback BR')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    a.minorticks_off()
    axes_disable_scientific(a)
    a.set_aspect('equal')
    #t, p = ttest_rel(brons, broffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
    #scalarmappable = mpl.cm.ScalarMappable(norm=None, cmap=cmap)
    #f.colorbar(scalarmappable, label='DSI')


# scatter plot grating meanrates, coloured by dLGN layer (shell/core):
figsize = DEFAULTFIGURESIZE
logmin, logmax = -1.2, 2
logticks = np.array([-1, 0, 1, 2])
#log0min = logmin + 0.05
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto meanrate grating %s layer' % st8)
    rons, roffs, layers = [], [], []
    for mseustr in grtmseustrs:
        meanrate = grtresp.loc[mseustr, st8]['meanrate']
        msustr = mseustr2msustr(mseustr)
        layer = celltype.loc[msustr]['layer']
        if meanrate.isna().any(): # missing one or both meanrates
            #print('%s: missing one or both opto conditions, skipping' % mseustr)
            continue
        rons.append(meanrate[True])
        roffs.append(meanrate[False])
        layers.append(layer)
        #fig3.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
    rons = np.asarray(rons)
    roffs = np.asarray(roffs)
    layers = np.asarray(layers)
    layerclrs = [ {'shell':violet, 'core':green, 'nan':'black'}[layer] for layer in layers ]
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(rons, roffs, marker='.', c='None', edgecolor=layerclrs, s=DEFSZ)
    # plot mean of points grouped by layer:
    shelllogmean = gmean(rons[layers == 'shell']), gmean(roffs[layers == 'shell'])
    corelogmean = gmean(rons[layers == 'core']), gmean(roffs[layers == 'core'])
    a.scatter(shelllogmean[0], shelllogmean[1], marker='+', c=violet, s=100)
    a.scatter(corelogmean[0], corelogmean[1], marker='+', c=green, s=100)
    a.set_ylabel('Feedback FR (spk/s)')
    a.set_xlabel('Suppression FR (spk/s)')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    axes_disable_scientific(a)
    a.minorticks_off()
    a.set_aspect('equal')
    #t, p = ttest_rel(rons, roffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
    patch = Rectangle((0, 0), 0, 0)
    l = a.legend([patch, patch], ['Shell', 'Core'], loc='upper left',
                 handlelength=0, handletextpad=0, frameon=False)
    for t, c in zip(l.get_texts(), [violet, green]):
        t.set_color(c)


# scatter plot movie meanrates, coloured by dLGN layer (shell/core):
figsize = DEFAULTFIGURESIZE
logmin, logmax = -1.2, 2
logticks = np.array([-1, 0, 1, 2])
#log0min = logmin + 0.05
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto meanrate movie nat %s layer' % st8)
    rons, roffs, layers = [], [], []
    for mseustr in mvimseustrs:
        meanrate = mviresp.loc[mseustr, 'nat', st8]['meanrate']
        msustr = mseustr2msustr(mseustr)
        layer = celltype.loc[msustr]['layer']
        if meanrate.isna().any(): # missing one or both meanrates
            #print('%s: missing one or both opto conditions, skipping' % mseustr)
            continue
        rons.append(meanrate[True])
        roffs.append(meanrate[False])
        layers.append(layer)
        #fig3.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
    rons = np.asarray(rons)
    roffs = np.asarray(roffs)
    layers = np.asarray(layers)
    layerclrs = [ {'shell':violet, 'core':green, 'nan':'black'}[layer] for layer in layers ]
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(rons, roffs, marker='.', c='None', edgecolor=layerclrs, s=DEFSZ)
    # plot mean of points grouped by layer:
    shelllogmean = gmean(rons[layers == 'shell']), gmean(roffs[layers == 'shell'])
    corelogmean = gmean(rons[layers == 'core']), gmean(roffs[layers == 'core'])
    a.scatter(shelllogmean[0], shelllogmean[1], marker='+', c=violet, s=100)
    a.scatter(corelogmean[0], corelogmean[1], marker='+', c=green, s=100)
    a.set_ylabel('Feedback FR (spk/s)')
    a.set_xlabel('Suppression FR (spk/s)')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    axes_disable_scientific(a)
    a.minorticks_off()
    a.set_aspect('equal')
    #t, p = ttest_rel(rons, roffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
    patch = Rectangle((0, 0), 0, 0)
    l = a.legend([patch, patch], ['Shell', 'Core'], loc='upper left',
                 handlelength=0, handletextpad=0, frameon=False)
    for t, c in zip(l.get_texts(), [violet, green]):
        t.set_color(c)


# scatter plot grating burst ratio, coloured by dLGN layer (shell/core):
figsize = DEFAULTFIGURESIZE
logmin, logmax = -3.55, 0
logticks = np.array([-3, -2, -1, 0])
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto burst ratio grating %s layer' % st8)
    brons, broffs, layers = [], [], []
    for mseustr in grtmseustrs:
        br = grtresp.loc[mseustr, st8]['meanburstratio']
        msustr = mseustr2msustr(mseustr)
        layer = celltype.loc[msustr]['layer']
        if br.isna().any(): # missing for at least one opto condition
            continue
        brons.append(br[True])
        broffs.append(br[False])
        layers.append(layer)
        #fig3.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
    brons, broffs, layers = np.asarray(brons), np.asarray(broffs), np.asarray(layers)
    layerclrs = [ {'shell':violet, 'core':green, 'nan':'black'}[layer] for layer in layers ]
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(brons, broffs, marker='.', c='None', edgecolor=layerclrs, s=DEFSZ)
    # plot mean of points grouped by layer:
    # filter out 0 vals before calculating geometric mean:
    shellis = layers == 'shell'
    coreis = layers == 'core'
    shellbrons = brons[shellis][brons[shellis] > 0]
    shellbroffs = broffs[shellis][broffs[shellis] > 0]
    corebrons = brons[coreis][brons[coreis] > 0]
    corebroffs = broffs[coreis][broffs[coreis] > 0]
    shelllogmean = gmean(shellbrons), gmean(shellbroffs)
    corelogmean = gmean(corebrons), gmean(corebroffs)
    a.scatter(shelllogmean[0], shelllogmean[1], marker='+', c=violet, s=100)
    a.scatter(corelogmean[0], corelogmean[1], marker='+', c=green, s=100)
    a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
    a.set_ylabel('Feedback BR')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    a.minorticks_off()
    axes_disable_scientific(a)
    a.set_aspect('equal')
    #t, p = ttest_rel(brons, broffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
    patch = Rectangle((0, 0), 0, 0)
    l = a.legend([patch, patch], ['Shell', 'Core'], loc='upper left',
                 handlelength=0, handletextpad=0, frameon=False)
    for t, c in zip(l.get_texts(), [violet, green]):
        t.set_color(c)


# scatter plot movie burst ratio, coloured by dLGN layer (shell/core):
figsize = DEFAULTFIGURESIZE
logmin, logmax = -3.55, 0
logticks = np.array([-3, -2, -1, 0])
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto burst ratio movie nat %s layer' % st8)
    brons, broffs, layers = [], [], []
    for mseustr in mvimseustrs:
        br = mviresp.loc[mseustr, 'nat', st8]['meanburstratio']
        msustr = mseustr2msustr(mseustr)
        layer = celltype.loc[msustr]['layer']
        if br.isna().any(): # missing for at least one opto condition
            continue
        brons.append(br[True])
        broffs.append(br[False])
        layers.append(layer)
        #fig3.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
    brons, broffs, layers = np.asarray(brons), np.asarray(broffs), np.asarray(layers)
    layerclrs = [ {'shell':violet, 'core':green, 'nan':'black'}[layer] for layer in layers ]
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(brons, broffs, marker='.', c='None', edgecolor=layerclrs, s=DEFSZ)
    # plot mean of points grouped by layer:
    # filter out 0 vals before calculating geometric mean:
    shellis = layers == 'shell'
    coreis = layers == 'core'
    shellbrons = brons[shellis][brons[shellis] > 0]
    shellbroffs = broffs[shellis][broffs[shellis] > 0]
    corebrons = brons[coreis][brons[coreis] > 0]
    corebroffs = broffs[coreis][broffs[coreis] > 0]
    shelllogmean = gmean(shellbrons), gmean(shellbroffs)
    corelogmean = gmean(corebrons), gmean(corebroffs)
    a.scatter(shelllogmean[0], shelllogmean[1], marker='+', c=violet, s=100)
    a.scatter(corelogmean[0], corelogmean[1], marker='+', c=green, s=100)
    a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
    a.set_ylabel('Feedback BR')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    a.minorticks_off()
    axes_disable_scientific(a)
    a.set_aspect('equal')
    #t, p = ttest_rel(brons, broffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
    patch = Rectangle((0, 0), 0, 0)
    l = a.legend([patch, patch], ['Shell', 'Core'], loc='upper left',
                 handlelength=0, handletextpad=0, frameon=False)
    for t, c in zip(l.get_texts(), [violet, green]):
        t.set_color(c)
'''
'''
# scatter plot spontaneous meanrates:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -1.2, 2
logticks = np.array([-1, 0, 1, 2])
#log0min = logmin + 0.05
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto meanrate spontaneous %s' % st8)
    rons, roffs = [], []
    for mseustr in sponmseustrs:
        meanrate = sponresp.loc[mseustr]['meanrate']
        if meanrate.isna().any(): # missing one or both meanrates
            #print('%s: missing one or both opto conditions, skipping' % mseustr)
            continue
        rons.append(meanrate[True])
        roffs.append(meanrate[False])
        #fig3.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
    rons = np.asarray(rons)
    roffs = np.asarray(roffs)
    # replace off-scale low values with log0min, so the points remain visible:
    #pltrons, pltroffs = rons.copy(), roffs.copy()
    #pltrons[pltrons <= 10**logmin] = 10**log0min
    #pltroffs[pltroffs <= 10**logmin] = 10**log0min
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    a.scatter(rons, roffs, marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
    a.set_xlabel('Suppression FR (spk/s)')
    a.set_ylabel('Feedback FR (spk/s)')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    a.minorticks_off()
    axes_disable_scientific(a)
    a.set_aspect('equal')
    t, p = ttest_rel(rons, roffs) # paired t-test
    a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
    mu = rons.mean(), roffs.mean()
    txt = '$\mathregular{\mu=%.1f, %.1f}$' % mu
    a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))


# scatter plot spontaneous burst ratio:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -3.55, 0
logticks = np.array([-3, -2, -1, 0])
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto burst ratio spontaneous %s' % st8)
    brons, broffs = [], []
    for mseustr in sponmseustrs:
        br = sponresp.loc[mseustr]['meanburstratio']
        if br.isna().any(): # missing for at least one opto condition
            continue
        brons.append(br[True])
        broffs.append(br[False])
        #fig3.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
    brons, broffs = np.asarray(brons), np.asarray(broffs)
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    a.scatter(brons, broffs, marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
    a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
    a.set_ylabel('Feedback BR')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    a.minorticks_off()
    axes_disable_scientific(a)
    a.set_aspect('equal')
    t, p = ttest_rel(brons, broffs) # paired t-test
    a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
'''

## movie onset ON/OFF/transient classification, also doesn't cluster, but doesn't seem
## to match classical sense of ON/OFF/transient from luminance step function in chirp
## stimulus (at least in Ntsr mice which had both movie and chirp stim, see fig1S33S1_ntsr.py),
## so leaving this out for now:
'''
# plot normalized PSTHs wrt movie onset coloured by ON/OFF/transient classification
figsize = 10, 5
#logmin, logmax = -1.2, 2
#logticks = np.array([-1, 0, 1, 2])
#log0min = logmin + 0.05
f, a = plt.subplots(figsize=figsize)
wintitle('mvionsetpsth nat ONOFFTHRESH=%.2f TRANSTHRESH=%.2f' % (ONOFFTHRESH, TRANSTHRESH))
clr2zorder = {'gray': -3, 'black':-2, green:-1, 'red':0}
for msustr, row in celltype.iterrows():
    onoff = row['onoff']
    trans = row['trans']
    mvionsetpsth = row['mvionsetpsth']
    t = row['t']
    clr = 'gray'
    if trans >= TRANSTHRESH: # first try and classify by trans
        clr = 'black'
    elif onoff >= ONOFFTHRESH: # try and classify by on/off
        clr = green
    elif onoff < -ONOFFTHRESH:
        clr = 'red'
    zorder = clr2zorder[clr]
    a.plot(t, mvionsetpsth, '-', color=clr, alpha=1, zorder=zorder)
a.axvline(x=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
a.axvline(x=0.1, ls='--', marker='', color='lightgray', zorder=-np.inf)
a.set_ylabel('Normalized firing rate')
a.set_xlabel('Time from movie onset (s)')
dt = np.diff(t).mean()
a.set_xlim(t[0], t[-1]+dt)
#a.set_ylim(-0.1, 2)
patch = Rectangle((0, 0), 0, 0)
l = a.legend([patch, patch, patch], ['On', 'Off', 'Transient'], loc='upper left',
             handlelength=0, handletextpad=0, frameon=False)
for t, c in zip(l.get_texts(), [green, 'red', 'black']):
    t.set_color(c)


# scatter plot grating meanrates, coloured by ON/OFF/transient:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -1.2, 2
logticks = np.array([-1, 0, 1, 2])
#log0min = logmin + 0.05
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto meanrate grating %s ONOFFTHRESH=%.2f TRANSTHRESH=%.2f'
             % (st8, ONOFFTHRESH, TRANSTHRESH))
    rons, roffs, onoffs, transs = [], [], [], []
    for mseustr in grtmseustrs:
        meanrate = grtresp.loc[mseustr, st8]['meanrate']
        msustr = mseustr2msustr(mseustr)
        onoff = celltype.loc[msustr]['onoff']
        trans = celltype.loc[msustr]['trans']
        if meanrate.isna().any(): # missing one or both meanrates
            #print('%s: missing one or both opto conditions, skipping' % mseustr)
            continue
        rons.append(meanrate[True])
        roffs.append(meanrate[False])
        onoffs.append(onoff)
        transs.append(trans)
        #fig3.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
    rons = np.asarray(rons)
    roffs = np.asarray(roffs)
    clrs = []
    for onoff, trans in zip(onoffs, transs):
        clr = 'gray'
        if trans >= TRANSTHRESH: # first try and classify by trans
            clr = 'black'
        elif onoff >= ONOFFTHRESH: # try and classify by on/off
            clr = green
        elif onoff < -ONOFFTHRESH:
            clr = 'red'
        clrs.append(clr)
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(rons, roffs, marker='.', c='None', edgecolor=clrs, s=DEFSZ)
    a.set_ylabel('Feedback FR (spk/s)')
    a.set_xlabel('Suppression FR (spk/s)')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    axes_disable_scientific(a)
    a.minorticks_off()
    a.set_aspect('equal')
    #t, p = ttest_rel(rons, roffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
    patch = Rectangle((0, 0), 0, 0)
    l = a.legend([patch, patch, patch], ['On', 'Off', 'Transient'], loc='upper left',
                 handlelength=0, handletextpad=0, frameon=False)
    for t, c in zip(l.get_texts(), [green, 'red', 'black']):
        t.set_color(c)


# scatter plot movie meanrates, coloured by ON/OFF/transient:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -1.2, 2
logticks = np.array([-1, 0, 1, 2])
#log0min = logmin + 0.05
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto meanrate movie nat %s ONOFFTHRESH=%.2f TRANSTHRESH=%.2f'
             % (st8, ONOFFTHRESH, TRANSTHRESH))
    rons, roffs, onoffs, trans = [], [], [], []
    for mseustr in mvimseustrs:
        meanrate = mviresp.loc[mseustr, 'nat', st8]['meanrate']
        msustr = mseustr2msustr(mseustr)
        onoff = celltype.loc[msustr]['onoff']
        trans = celltype.loc[msustr]['trans']
        if meanrate.isna().any(): # missing one or both meanrates
            #print('%s: missing one or both opto conditions, skipping' % mseustr)
            continue
        rons.append(meanrate[True])
        roffs.append(meanrate[False])
        onoffs.append(onoff)
        transs.append(trans)
        #fig3.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
    rons = np.asarray(rons)
    roffs = np.asarray(roffs)
    clrs = []
    for onoff, trans in zip(onoffs, transs):
        clr = 'gray'
        if trans >= TRANSTHRESH: # first try and classify by trans
            clr = 'black'
        elif onoff >= ONOFFTHRESH: # try and classify by on/off
            clr = green
        elif onoff < -ONOFFTHRESH:
            clr = 'red'
        clrs.append(clr)
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(rons, roffs, marker='.', c='None', edgecolor=clrs, s=DEFSZ)
    a.set_ylabel('Feedback FR (spk/s)')
    a.set_xlabel('Suppression FR (spk/s)')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    axes_disable_scientific(a)
    a.minorticks_off()
    a.set_aspect('equal')
    #t, p = ttest_rel(rons, roffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
    patch = Rectangle((0, 0), 0, 0)
    l = a.legend([patch, patch, patch], ['On', 'Off', 'Transient'], loc='upper left',
                 handlelength=0, handletextpad=0, frameon=False)
    for t, c in zip(l.get_texts(), [green, 'red', 'black']):
        t.set_color(c)


# scatter plot grating burst ratio, coloured by ON/OFF/transient:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -3.55, 0
logticks = np.array([-3, -2, -1, 0])
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto burst ratio grating %s ONOFFTHRESH=%.2f TRANSTHRESH=%.2f'
             % (st8, ONOFFTHRESH, TRANSTHRESH))
    brons, broffs, onoffs, trans = [], [], [], []
    for mseustr in grtmseustrs:
        br = grtresp.loc[mseustr, st8]['meanburstratio']
        msustr = mseustr2msustr(mseustr)
        onoff = celltype.loc[msustr]['onoff']
        trans = celltype.loc[msustr]['trans']
        if br.isna().any(): # missing for at least one opto condition
            continue
        brons.append(br[True])
        broffs.append(br[False])
        onoffs.append(onoff)
        transs.append(trans)
    brons, broffs = np.asarray(brons), np.asarray(broffs)
    clrs = []
    for onoff, trans in zip(onoffs, transs):
        clr = 'gray'
        if trans >= TRANSTHRESH: # first try and classify by trans
            clr = 'black'
        elif onoff >= ONOFFTHRESH: # try and classify by on/off
            clr = green
        elif onoff < -ONOFFTHRESH:
            clr = 'red'
        clrs.append(clr)
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(brons, broffs, marker='.', c='None', edgecolor=clrs, s=DEFSZ)
    a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
    a.set_ylabel('Feedback BR')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    a.minorticks_off()
    axes_disable_scientific(a)
    a.set_aspect('equal')
    #t, p = ttest_rel(brons, broffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
    patch = Rectangle((0, 0), 0, 0)
    l = a.legend([patch, patch, patch], ['On', 'Off', 'Transient'], loc='upper left',
                 handlelength=0, handletextpad=0, frameon=False)
    for t, c in zip(l.get_texts(), [green, 'red', 'black']):
        t.set_color(c)


# scatter plot movie burst ratio, coloured by ON/OFF/transient:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -3.55, 0
logticks = np.array([-3, -2, -1, 0])
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto burst ratio movie nat %s ONOFFTHRESH=%.2f TRANSTHRESH=%.2f'
             % (st8, ONOFFTHRESH, TRANSTHRESH))
    brons, broffs, onoffs, trans = [], [], [], []
    for mseustr in mvimseustrs:
        br = mviresp.loc[mseustr, 'nat', st8]['meanburstratio']
        msustr = mseustr2msustr(mseustr)
        onoff = celltype.loc[msustr]['onoff']
        trans = celltype.loc[msustr]['trans']
        if br.isna().any(): # missing for at least one opto condition
            continue
        brons.append(br[True])
        broffs.append(br[False])
        onoffs.append(onoff)
        transs.append(trans)
        #fig3.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
    brons, broffs = np.asarray(brons), np.asarray(broffs)
    for onoff, trans in zip(onoffs, transs):
        clr = 'gray'
        if trans >= TRANSTHRESH: # first try and classify by trans
            clr = 'black'
        elif onoff >= ONOFFTHRESH: # try and classify by on/off
            clr = green
        elif onoff < -ONOFFTHRESH:
            clr = 'red'
        clrs.append(clr)
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(brons, broffs, marker='.', c='None', edgecolor=clrs, s=DEFSZ)
    a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
    a.set_ylabel('Feedback BR')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    a.minorticks_off()
    axes_disable_scientific(a)
    a.set_aspect('equal')
    #t, p = ttest_rel(brons, broffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
    patch = Rectangle((0, 0), 0, 0)
    l = a.legend([patch, patch, patch], ['On', 'Off', 'Transient'], loc='upper left',
                 handlelength=0, handletextpad=0, frameon=False)
    for t, c in zip(l.get_texts(), [green, 'red', 'black']):
        t.set_color(c)


# scatter plot grating meanrates, coloured by distance of MUA envl RF from center of screen:
figsize = 3.5, DEFAULTFIGURESIZE[1] #DEFAULTFIGURESIZE
logmin, logmax = -1.2, 2
logticks = np.array([-1, 0, 1, 2])
#log0min = logmin + 0.05
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto meanrate grating %s rfdist' % st8)
    rons, roffs, ds = [], [], []
    for mseustr in grtmseustrs:
        meanrate = grtresp.loc[mseustr, st8]['meanrate']
        msustr = mseustr2msustr(mseustr)
        x0, y0 = cellscreenpos.loc[msustr]
        d = np.sqrt(x0**2 + y0**2) # distance from screen center, deg
        if meanrate.isna().any(): # missing one or both meanrates
            #print('%s: missing one or both opto conditions, skipping' % mseustr)
            continue
        rons.append(meanrate[True])
        roffs.append(meanrate[False])
        ds.append(d)
        #fig3.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
    rons = np.asarray(rons)
    roffs = np.asarray(roffs)
    ds = np.asarray(ds)
    maxds = intround(np.nanmax(ds))
    normds = ds / maxds
    cmap = plt.get_cmap() # use default, normalized from 0 to 1
    dclrs = []
    for normd in normds:
        if np.isnan(normd):
            dclrs.append('gray')
        else:
            dclrs.append(cmap(normd))
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(rons, roffs, marker='.', c='none', edgecolor=dclrs, s=DEFSZ, cmap='viridis')
    a.set_ylabel('Feedback FR (spk/s)')
    a.set_xlabel('Suppression FR (spk/s)')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    axes_disable_scientific(a)
    a.minorticks_off()
    a.set_aspect('equal')
    #t, p = ttest_rel(rons, roffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
    scalarmappable = mpl.cm.ScalarMappable(norm=None, cmap=cmap)
    scalarmappable.set_clim(0, maxds)
    f.colorbar(scalarmappable, ticks=[0, maxds], label='Dist. from center ($\degree$)')


# scatter plot movie meanrates, coloured by distance of MUA envl RF from center of screen:
figsize = 3.5, DEFAULTFIGURESIZE[1] #DEFAULTFIGURESIZE
logmin, logmax = -1.2, 2
logticks = np.array([-1, 0, 1, 2])
#log0min = logmin + 0.05
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto meanrate movie nat %s rfdist' % st8)
    rons, roffs, ds = [], [], []
    for mseustr in mvimseustrs:
        meanrate = mviresp.loc[mseustr, 'nat', st8]['meanrate']
        msustr = mseustr2msustr(mseustr)
        x0, y0 = cellscreenpos.loc[msustr]
        d = np.sqrt(x0**2 + y0**2) # distance from screen center, deg
        if meanrate.isna().any(): # missing one or both meanrates
            #print('%s: missing one or both opto conditions, skipping' % mseustr)
            continue
        rons.append(meanrate[True])
        roffs.append(meanrate[False])
        ds.append(d)
        #fig3.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
    rons = np.asarray(rons)
    roffs = np.asarray(roffs)
    ds = np.asarray(ds)
    maxds = intround(np.nanmax(ds))
    normds = ds / maxds
    cmap = plt.get_cmap() # use default, normalized from 0 to 1
    dclrs = []
    for normd in normds:
        if np.isnan(normd):
            dclrs.append('gray')
        else:
            dclrs.append(cmap(normd))
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(rons, roffs, marker='.', c='None', edgecolor=dclrs, s=DEFSZ)
    a.set_ylabel('Feedback FR (spk/s)')
    a.set_xlabel('Suppression FR (spk/s)')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    axes_disable_scientific(a)
    a.minorticks_off()
    a.set_aspect('equal')
    #t, p = ttest_rel(rons, roffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
    scalarmappable = mpl.cm.ScalarMappable(norm=None, cmap=cmap)
    scalarmappable.set_clim(0, maxds)
    f.colorbar(scalarmappable, ticks=[0, maxds], label='Dist. from center ($\degree$)')


# scatter plot grating burst ratio, coloured by distance of MUA envl RF from center of screen:
figsize = 3.6, DEFAULTFIGURESIZE[1] #DEFAULTFIGURESIZE
logmin, logmax = -3.55, 0
logticks = np.array([-3, -2, -1, 0])
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto burst ratio grating %s rfdist' % st8)
    brons, broffs, ds = [], [], []
    for mseustr in grtmseustrs:
        br = grtresp.loc[mseustr, st8]['meanburstratio']
        msustr = mseustr2msustr(mseustr)
        x0, y0 = cellscreenpos.loc[msustr]
        d = np.sqrt(x0**2 + y0**2) # distance from screen center, deg
        if br.isna().any(): # missing for at least one opto condition
            continue
        brons.append(br[True])
        broffs.append(br[False])
        ds.append(d)
        #fig3.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
    brons, broffs, ds = np.asarray(brons), np.asarray(broffs), np.asarray(ds)
    maxds = intround(np.nanmax(ds))
    normds = ds / maxds
    cmap = plt.get_cmap() # use default, normalized from 0 to 1
    dclrs = []
    for normd in normds:
        if np.isnan(normd):
            dclrs.append('gray')
        else:
            dclrs.append(cmap(normd))
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(brons, broffs, marker='.', c='None', edgecolor=dclrs, s=DEFSZ)
    a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
    a.set_ylabel('Feedback BR')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    a.minorticks_off()
    axes_disable_scientific(a)
    a.set_aspect('equal')
    #t, p = ttest_rel(brons, broffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
    scalarmappable = mpl.cm.ScalarMappable(norm=None, cmap=cmap)
    scalarmappable.set_clim(0, maxds)
    f.colorbar(scalarmappable, ticks=[0, maxds], label='Dist. from center ($\degree$)')


# scatter plot movie burst ratio, coloured by distance of MUA envl RF from center of screen:
figsize = 3.6, DEFAULTFIGURESIZE[1] #DEFAULTFIGURESIZE
logmin, logmax = -3.55, 0
logticks = np.array([-3, -2, -1, 0])
for st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto burst ratio movie nat %s rfdist' % st8)
    brons, broffs, ds = [], [], []
    for mseustr in mvimseustrs:
        br = mviresp.loc[mseustr, 'nat', st8]['meanburstratio']
        msustr = mseustr2msustr(mseustr)
        x0, y0 = cellscreenpos.loc[msustr]
        d = np.sqrt(x0**2 + y0**2) # distance from screen center, deg
        if br.isna().any(): # missing for at least one opto condition
            continue
        brons.append(br[True])
        broffs.append(br[False])
        ds.append(d)
        #fig3.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
    brons, broffs, ds = np.asarray(brons), np.asarray(broffs), np.asarray(ds)
    maxds = intround(np.nanmax(ds))
    normds = ds / maxds
    cmap = plt.get_cmap() # use default, normalized from 0 to 1
    dclrs = []
    for normd in normds:
        if np.isnan(normd):
            dclrs.append('gray')
        else:
            dclrs.append(cmap(normd))
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot all points:
    a.scatter(brons, broffs, marker='.', c='None', edgecolor=dclrs, s=DEFSZ)
    a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
    a.set_ylabel('Feedback BR')
    a.set_xscale('log')
    a.set_yscale('log')
    a.set_xlim(10**logmin, 10**logmax)
    a.set_ylim(10**logmin, 10**logmax)
    a.set_xticks(10.0**logticks)
    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
    a.minorticks_off()
    axes_disable_scientific(a)
    a.set_aspect('equal')
    #t, p = ttest_rel(brons, broffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
    scalarmappable = mpl.cm.ScalarMappable(norm=None, cmap=cmap)
    scalarmappable.set_clim(0, maxds)
    f.colorbar(scalarmappable, ticks=[0, maxds], label='Dist. from center ($\degree$)')
'''