natfeedback_code_eLife/fig1.py

"""Figure 1 and 1S2 plots, and first and last 2 sec scatter plots for 4S1, use run -i fig1.py"""

mi = pd.MultiIndex.from_product([mvimseustrs, OPTOS], names=['mseu', 'opto'])
# include single trial and trial-averaged columns:
columns = ['trialis', 'rates', 'rate02s', 'rate35s', 'burstratios',
           'blankrates', 'blankburstratios'] + modmeasuresnoblankcond
fig1 = pd.DataFrame(index=mi, columns=columns) # excludes kind and st8
## NOTE: make sure that rows are not overwritten due to temporarily enabling iteration
## over kind and/or st8 in any of the following loops!


# plot movie rasters for each unit by kind, state and opto condition:
showbursts = True
SCATTER = True # True: low quality, but fast; False: high quality, but slow
trialiis = None #np.arange(120) # plot only these trials, None plots all trials
for mseustr in []:#mvimseu2exmpli: #mvimseustrs:
    #['Ntsr1Cre_2019_0008_s03_e04_u18']
    #['Ntsr1Cre_2020_0001_s04_e10_u12']
    for kind in ['nat']: #MVIKINDS
        for st8 in ['none']:
            for opto in OPTOS:
                dt = mviresp.loc[mseustr, kind, st8, opto]['dt'] # stimulus duration (s)
                wdt = -OFFSETS[0] + dt + OFFSETS[1] # wide raster duration, sec
                tmin, tmax = 0 + OFFSETS[0], dt + OFFSETS[1]
                optotrange = mviresp.loc[mseustr, kind, st8, opto]['optotrange']
                raster = mviresp.loc[mseustr, kind, st8, opto]['wraster']
                if np.any(pd.isna(raster)): # no raster for this condition
                    continue
                title = '%s movie %s %s %s raster' % (mseustr, kind, st8, opto)
                if trialiis is not None:
                    raster = raster[trialiis]
                    title += ' %dtrials' % len(trialiis)
                burstis = None
                if showbursts:
                    wburstis = mviresp.loc[mseustr, kind, st8, opto]['wburstis']
                a = simpletraster(raster=raster, alpha=0.5, s=1, dt=dt, offsets=OFFSETS,
                                  scatter=SCATTER,
                                  title=title, inchespersec=1.5*RASTERSCALEX,
                                  inchespertrial=1/25*RASTERSCALEX,
                                  xaxis=True, burstis=wburstis)
                if opto:
                    # plot horizontal line signifying opto ON period, just above axes:
                    ton, toff = optotrange
                    a.hlines(-5, ton, toff, colors=optoblue, lw=4, clip_on=False,
                             in_layout=False)
                # plot horizontal line signifying stimulus period, just above axes:
                a.hlines(-2, 0, dt, colors='k', lw=4, clip_on=False, in_layout=False)
                # plot vertical coloured line just right of axes:
                #vlinexpos = dt + OFFSETS[1] + wdt*0.01
                #y0, y1 = a.get_ylim()
                #clr = st82clr[st8]
                #alpha = opto2alpha[opto]
                #a.vlines(vlinexpos, y0, y1, colors=clr, lw=4, alpha=alpha, clip_on=False,
                #         in_layout=False)
                a.set_xlim(tmin, tmax)


# plot movie PSTHs for each unit by kind, overplot state and opto combinations:
for mseustr in []: #mvimseu2exmpli:#['Ntsr1Cre_2020_0001_s04_e10_u12']: #mvimseustrs
    for kind in MVIKINDS:
        wt = mviresp.loc[mseustr, kind]['wt'] # DataFrame
        wpsth = mviresp.loc[mseustr, kind]['wpsth'] # DataFrame
        wbpsth = mviresp.loc[mseustr, kind]['wbpsth'] # DataFrame
        #t = mviresp.loc[mseustr, kind]['t'] # DataFrame
        psth = mviresp.loc[mseustr, kind]['psth'] # DataFrame
        if np.isnan(np.hstack(wpsth)).any(): # no PSTH for at least one condition for this mseu
            continue
        optorho, optorho_p = scipy.stats.pearsonr(psth['none', False], psth['none', True])
        runsitrho, runsitrho_p = scipy.stats.pearsonr(psth['run', False], psth['sit', False])
        skewoff = scipy.stats.skew(psth['none', False])
        skewon = scipy.stats.skew(psth['none', True])
        skewrun = scipy.stats.skew(psth['run', False])
        skewsit = scipy.stats.skew(psth['sit', False])
        _, optoKS_p = ks_2samp(psth['none', False], psth['none', True])
        _, runsitKS_p = ks_2samp(psth['run', False], psth['sit', False])
        print('%s example PSTH:\n'
              'optorho=%g optorho_p=%g skewoff=%g skewon=%g optoKS_p=%g\n'
              'runsitrho=%g runsitrho_p=%g skewrun=%g skewsit=%g runsitKS_p=%g'
              % (mseustr,
                 optorho, optorho_p, skewoff, skewon, optoKS_p,
                 runsitrho, runsitrho_p, skewrun, skewsit, runsitKS_p))
        for st8combo, optocombo in zip(ST8COMBOS, OPTOCOMBOS):
            st8combostr = ' '.join(st8combo)
            optocombostr = ' '.join([str(opto) for opto in optocombo])
            dt = mviresp.loc[mseustr, kind, 'none']['dt'].mean() # mean stimulus duration (s)
            wdt = -OFFSETS[0] + dt + OFFSETS[1] # wide raster duration, sec
            tmin, tmax = 0 + OFFSETS[0], dt + OFFSETS[1]
            figsize = 0.872 + wdt*1.5*RASTERSCALEX, PSTHHEIGHT
            f, a = plt.subplots(figsize=figsize)
            title = '%s %s %s %s PSTH' % (mseustr, kind, st8combostr, optocombostr)
            wintitle(title)
            # subsample as in PSTH scatter plots, and plot all 4 st8/opto combinations:
            for st8 in st8combo:
                color = st82clr[st8]
                for opto in optocombo:
                    if st8combo == ['run', 'sit'] and optocombo == [True]:
                        optoalpha = 1
                    else:
                        optoalpha = opto2alpha[opto]
                    if optocombo == [False]: # desaturate bursts by st8, not opto
                        burstalpha = st82alpha[st8]
                    else:
                        burstalpha = opto2alpha[opto]
                    c = desat(color, optoalpha) # do manual alpha mixing
                    bc = desat(burstclr, burstalpha) # do manual alpha mixing
                    fb = opto2fb[opto].title()
                    if st8 == 'none':
                        label = fb
                    else:
                        label = ', '.join([fb, st8.title()])
                    a.plot(wt[st8, opto][::ssx], wpsth[st8, opto][::ssx],
                           '-', color=c, label=label)
                    a.plot(wt[st8, opto][::ssx], wbpsth[st8, opto][::ssx],
                           '-', color=bc, label=label+', Burst')
            l = a.legend(frameon=False)
            l.set_draggable(True)
            a.set_xlabel('Time (s)')
            a.set_ylabel('Firing rate (spk/s)')
            a.set_xlim(tmin, tmax)
            # plot horizontal line signifying stimulus period, just below data:
            ymin, ymax = a.get_ylim()
            a.hlines(-ymax*0.025, 0, dt, colors='k', lw=4, clip_on=False, in_layout=False)
            a.set_ylim(ymin, ymax) # restore y limits from before plotting the hline


# plot PDFs, CDFs and QQ of PSTHs for each unit by kind, overplot state and opto combinations:
percentiles = np.arange(0, 100, 2)
figsize = DEFAULTFIGURESIZE
for mseustr in mvimseu2exmpli: #mvimseustrs
    if EXPTYPE != 'pvmvis':
        break # skip these plots
    for kind in MVIKINDS:
        psths = mviresp.loc[mseustr, kind]['psth']
        if np.isnan(np.hstack(psths)).any(): # no PSTH for at least one condition for this mseu
            continue
        for st8combo, optocombo in zip(ST8COMBOS, OPTOCOMBOS):
            st8combostr = ' '.join(st8combo)
            optocombostr = ' '.join([str(opto) for opto in optocombo])
            pdff, pdfa = plt.subplots(figsize=figsize) # PDF
            title = '%s %s %s %s PSTH PDF' % (mseustr, kind, st8combostr, optocombostr)
            wintitle(title)
            hpdff, hpdfa = plt.subplots(figsize=(figsize[0], 3)) # horizontal, match PSTH height
            title = '%s %s %s %s horizontal PSTH PDF' % (mseustr, kind, st8combostr, optocombostr)
            wintitle(title)
            npdff, npdfa = plt.subplots(figsize=figsize) # normalized PDF
            title = '%s %s %s %s PSTH nPDF' % (mseustr, kind, st8combostr, optocombostr)
            wintitle(title)
            cdff, cdfa = plt.subplots(figsize=figsize)
            title = '%s %s %s %s PSTH CDF' % (mseustr, kind, st8combostr, optocombostr)
            wintitle(title)
            qqf, qqa = plt.subplots(figsize=figsize)
            title = '%s %s %s %s PSTH QQ' % (mseustr, kind, st8combostr, optocombostr)
            wintitle(title)
            # plot PDFs, CDFs and QQ of all st8/opto combinations:
            for st8 in st8combo:
                c = st82clr[st8]
                qq = {}
                for opto in optocombo:
                    psth = psths[st8, opto][::ssx]
                    npsth = psth / psth.max() # normalize
                    uvals = np.unique(npsth)
                    qq[opto] = np.percentile(psth, percentiles)
                    fb = opto2fb[opto].title()
                    if st8 == 'none':
                        label = fb
                    else:
                        label = ', '.join([fb, st8.title()])
                    alpha = opto2alpha[opto]
                    nbins = intround(np.sqrt(len(psth)))
                    pdfa.hist(psth, bins=nbins, density=True, histtype='step',
                              color=c, alpha=alpha, label=label)
                    hpdfa.hist(psth, bins=nbins, density=True, histtype='step',
                               color=c, alpha=alpha, label=label, orientation='horizontal')
                    npdfa.hist(npsth, bins=nbins, density=True, histtype='step',
                               color=c, alpha=alpha, label=label)
                    cdfa.hist(npsth, bins=uvals, density=True, cumulative=True, histtype='step',
                              color=c, alpha=alpha, label=label)
                if len(optocombo) == 2:
                    qqa.scatter(qq[False], qq[True], marker='.', c=c)
            #l = a.legend(frameon=False)
            #l.set_draggable(True)
            pdfa.set_xlabel('Firing rate (spk/s)')
            pdfa.set_ylabel('Probability density')
            #pdfa.set_yticks([0, 0.5, 1])
            pdff.show()
            hpdfa.set_xlabel('Probability density')
            hpdfa.set_ylabel('Firing rate (spk/s)')
            hpdff.show()
            npdfa.set_xlabel('Normalized rate')
            npdfa.set_ylabel('Probability density')
            #pdfa.set_yticks([0, 0.5, 1])
            npdff.show()
            cdfa.set_xlabel('Normalized rate')
            cdfa.set_ylabel('Cumulative probability')
            cdfa.set_yticks([0, 0.5, 1])
            cdff.show()
            if len(optocombo) == 2:
                qqmax = np.ceil(np.array(qqa.dataLim).max())
                qqxyline = [0, qqmax], [0, qqmax]
                qqa.plot(qqxyline[0], qqxyline[1], '--', color='gray', zorder=-1)
                qqa.set_xlabel('Feedback rate (spk/s)')
                qqa.set_ylabel('Suppression rate (spk/s)')
                qqa.set_xlim([0, qqmax])
                qqa.set_ylim([0, qqmax])
                qqa.set_aspect('equal')
                #qqa.set_xticks([0, 0.5, 1])
                #qqa.set_yticks([0, 0.5, 1])
                qqf.show()


# plot one PSTH PDF, CDF and QQ, for each combination of kind, state and opto,
# each pooled across all mseu:
percentiles = np.arange(0, 100, 2)
figsize = DEFAULTFIGURESIZE
for kind in MVIKINDS:
    if EXPTYPE != 'pvmvis':
        break # skip these plots
    for st8combo in ST8COMBOS:
        pdff, pdfa = plt.subplots(figsize=figsize)
        title = 'mean %s %s PSTH PDF' % (kind, ' '.join(st8combo))
        wintitle(title)
        cdff, cdfa = plt.subplots(figsize=figsize)
        title = 'mean %s %s PSTH CDF' % (kind, ' '.join(st8combo))
        wintitle(title)
        qqf, qqa = plt.subplots(figsize=figsize)
        title = 'mean %s %s PSTH QQ' % (kind, ' '.join(st8combo))
        wintitle(title)
        for st8 in st8combo:
            c = st82clr[st8]
            qq = {}
            for opto in OPTOS:
                psths, npsths = [], []
                for mseustr in mvimseustrs: #mvimseu2exmpli
                    psth = mviresp.loc[mseustr, kind, st8, opto]['psth']
                    if np.isnan(psth).any(): # no PSTH for this mseu
                        continue
                    #psth = psth[::ssx] # no need to subsample, only calculating a few CDFs
                    psths.append(psth)
                    npsths.append(psth / psth.max()) # normalize
                psths, npsths = np.concatenate(psths), np.concatenate(npsths)
                uvals = np.unique(npsths)
                qq[opto] = np.percentile(psths, percentiles)
                fb = opto2fb[opto].title()
                if st8 == 'none':
                    label = fb
                else:
                    label = ', '.join([fb, st8.title()])
                alpha = opto2alpha[opto]
                pdfa.hist(npsths, bins=50, density=True, histtype='step',
                          color=c, alpha=alpha, label=label)
                cdfa.hist(npsths, uvals, density=True, cumulative=True, histtype='step',
                          color=c, alpha=alpha, label=label)
            qqa.scatter(qq[False], qq[True], marker='.', c=c)
        #l = a.legend(frameon=False)
        #l.set_draggable(True)
        pdfa.set_xlabel('Normalized rate')
        pdfa.set_ylabel('Probability density')
        #pdfa.set_yticks([0, 0.5, 1])
        pdff.show()
        cdfa.set_xlabel('Normalized rate')
        cdfa.set_ylabel('Cumulative probability')
        cdfa.set_yticks([0, 0.5, 1])
        cdff.show()
        qqmax = np.ceil(np.array(qqa.dataLim).max())
        qqxyline = [0, qqmax], [0, qqmax]
        qqa.plot(qqxyline[0], qqxyline[1], '--', color='gray', zorder=-1)
        qqa.set_xlabel('Feedback rate (spk/s)')
        qqa.set_ylabel('Suppression rate (spk/s)')
        #qqa.set_aspect('equal', 'datalim')
        qqa.set_xlim([0, qqmax])
        qqa.set_ylim([0, qqmax])
        #qqa.set_xticks([0, 0.5, 1])
        #print('%s xticks: %r' % (title, qqa.get_xticks()))
        #qqa.set_yticks(qqa.get_xticks())
        qqf.show()


# scatter plot movie meanrates:
figsize = DEFAULTFIGURESIZE[0]*1.02, DEFAULTFIGURESIZE[1] # tweak to make space for log units
logmin, logmax = -1, 2
logticks = np.array([-1, 0, 1, 2])
#log0min = logmin + 0.05
for kind in ['nat']:#MVIKINDS:
    for st8 in ['none']:#ALLST8S:
        f, a = plt.subplots(figsize=figsize)
        wintitle('opto meanrate movie %s %s' % (kind, st8))
        rons, roffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in mvimseustrs:
            meanrate = mviresp.loc[mseustr, kind, st8]['meanrate']
            snr = mviresp.loc[mseustr, kind, st8]['snr']
            if meanrate.isna().any(): # missing one or both meanrates
                continue
            if (snr < SNRTHRESH).all(): # neither condition has decent SNR
                continue
            rons.append(meanrate[True])
            roffs.append(meanrate[False])
            trialis = mviresp.loc[mseustr, kind, st8]['trialis']
            rates = mviresp.loc[mseustr, kind, st8]['rates']
            _, pval = ttest_ind(rates[False], rates[True], equal_var=False)
            sgnfs.append(pval < SCATTERPTHRESH) # bool
            fig1.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
            fig1.loc[mseustr]['trialis'] = trialis # save, for joining DataFrames trial-wise
            fig1.loc[mseustr]['rates'] = rates # save trial-wise values
            if mvimseu2exmpli.get(mseustr) == fig1exmpli:
                exmplis.append(keptmseui)
                exmplmseustrs.append(mseustr)
            else:
                normlis.append(keptmseui)
            keptmseui += 1 # manually increment
        rons = np.asarray(rons)
        roffs = np.asarray(roffs)
        sgnfs = np.asarray(sgnfs)
        # 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)
        insgnfis, = np.where(sgnfs == False)
        sgnfis, = np.where(sgnfs == True)
        normlinsgnfis = np.intersect1d(normlis, insgnfis)
        normlsgnfis = np.intersect1d(normlis, sgnfis)
        if len(normlinsgnfis) > 0:
            # plot normal (non-example) insignificant points:
            c = desat(st82clr[st8], SGNF2ALPHA[False]) # do manual alpha mixing
            a.scatter(rons[normlinsgnfis], roffs[normlinsgnfis], clip_on=False,
                      marker='.', c='None', edgecolor=c, s=DEFSZ)
        if len(normlsgnfis) > 0:
            # plot normal (non-example) significant points:
            c = desat(st82clr[st8], SGNF2ALPHA[True]) # do manual alpha mixing
            a.scatter(rons[normlsgnfis], roffs[normlsgnfis], clip_on=False,
                      marker='.', c='None', edgecolor=c, s=DEFSZ)
        # plot example points, one at a time:
        for exmpli, mseustr in zip(exmplis, exmplmseustrs):
            marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
            c = exmpli2clr[mvimseu2exmpli[mseustr]]
            sz = exmpli2sz[mvimseu2exmpli[mseustr]]
            lw = exmpli2lw[mvimseu2exmpli[mseustr]]
            a.scatter(rons[exmpli], roffs[exmpli], clip_on=False,
                      marker=marker, c=c, s=sz, lw=lw)
        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')
        a.spines['left'].set_position(('outward', 4))
        a.spines['bottom'].set_position(('outward', 4))
        #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))
        # plot mean firing rate PDF:
        f, a = plt.subplots(figsize=figsize)
        wintitle('opto meanrate PDF movie %s %s' % (kind, st8))
        bins = np.logspace(-1, 2, num=25, base=10) # end inclusive
        a.hist(roffs, bins=bins, histtype='step', density=False, color='k',
               label='Feedback', clip_on=False)
        a.hist(rons, bins=bins, histtype='step', density=False, color=optoblue,
               label='V1 Suppression', clip_on=False)
        a.set_xlabel('Mean firing rate (spk/s)')
        a.set_ylabel('Cell count')
        a.set_xscale('log')
        a.set_xlim(10**logmin, 10**logmax)
        a.set_xticks(10.0**logticks)
        a.minorticks_off()
        axes_disable_scientific(a)
        #l = a.legend(frameon=False)
        #l.set_draggable(True)



# scatter plot movie burst ratio:
figsize = DEFAULTFIGURESIZE[0]*1.04, DEFAULTFIGURESIZE[1] # tweak to make space for log units
logmin, logmax = -3.55, 0
if EXPTYPE == 'ntsrmvis':
    logmin, logmax = -3.83, 0
logticks = np.array([-3, -2, -1, 0])
for kind in ['nat']:#MVIKINDS:
    for st8 in ['none']:#ALLST8S:
        f, a = plt.subplots(figsize=figsize)
        wintitle('opto burst ratio movie %s %s' % (kind, st8))
        brons, broffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in mvimseustrs:
            br = mviresp.loc[mseustr, kind, st8]['meanburstratio']
            snr = mviresp.loc[mseustr, kind, st8]['snr']
            if br.isna().any(): # missing for at least one opto condition
                continue
            if (snr < SNRTHRESH).all(): # neither condition has decent SNR
                continue
            brons.append(br[True])
            broffs.append(br[False])
            burstratios = mviresp.loc[mseustr, kind, st8]['burstratios']
            _, pval = ttest_ind(burstratios[False], burstratios[True], equal_var=False)
            sgnfs.append(pval < SCATTERPTHRESH) # bool
            fig1.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
            fig1.loc[mseustr]['burstratios'] = burstratios # save trial-wise values
            if mvimseu2exmpli.get(mseustr) == fig1exmpli:
                exmplis.append(keptmseui)
                exmplmseustrs.append(mseustr)
            else:
                normlis.append(keptmseui)
            keptmseui += 1 # manually increment
        brons, broffs = np.asarray(brons), np.asarray(broffs)
        sgnfs = np.asarray(sgnfs)
        # plot y=x line:
        xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
        a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
        insgnfis, = np.where(sgnfs == False)
        sgnfis, = np.where(sgnfs == True)
        normlinsgnfis = np.intersect1d(normlis, insgnfis)
        normlsgnfis = np.intersect1d(normlis, sgnfis)
        if len(normlinsgnfis) > 0:
            # plot normal (non-example) insignificant points:
            c = desat(st82clr[st8], SGNF2ALPHA[False]) # do manual alpha mixing
            a.scatter(brons[normlinsgnfis], broffs[normlinsgnfis], clip_on=True,
                      marker='.', c='None', edgecolor=c, s=DEFSZ) # clip_on=False fails
        if len(normlsgnfis) > 0:
            # plot normal (non-example) significant points:
            c = desat(st82clr[st8], SGNF2ALPHA[True]) # do manual alpha mixing
            a.scatter(brons[normlsgnfis], broffs[normlsgnfis], clip_on=True,
                      marker='.', c='None', edgecolor=c, s=DEFSZ) # clip_on=False fails
        # plot example points, one at a time:
        for exmpli, mseustr in zip(exmplis, exmplmseustrs):
            marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
            c = exmpli2clr[mvimseu2exmpli[mseustr]]
            sz = exmpli2sz[mvimseu2exmpli[mseustr]]
            lw = exmpli2lw[mvimseu2exmpli[mseustr]]
            a.scatter(brons[exmpli], broffs[exmpli], clip_on=False,
                      marker=marker, c=c, s=sz, lw=lw)
        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')
        a.spines['left'].set_position(('outward', 4))
        a.spines['bottom'].set_position(('outward', 4))
        #t, p = ttest_rel(brons, broffs) # paired t-test
        #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))


# scatter plot movie PSTH sparseness:
figsize = DEFAULTFIGURESIZE
linmin, linmax, linstep = 0, 1, 0.5
ticks = np.arange(intround(linmin), intround(linmax)+linstep, linstep)
for kind in ['nat']:#MVIKINDS:
    for st8 in ['none']:#ALLST8S:
        f, a = plt.subplots(figsize=figsize)
        wintitle('opto sparseness %s %s' % (kind, st8))
        sparsons, sparsoffs, exmplis, exmplmseustrs, normlis = [], [], [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in mvimseustrs:
            spars = mviresp.loc[mseustr, kind, st8]['spars']
            snr = mviresp.loc[mseustr, kind, st8]['snr']
            if spars.isna().any(): # missing for at least one opto condition
                continue
            if (snr < SNRTHRESH).all(): # neither condition has decent SNR
                continue
            sparsons.append(spars[True])
            sparsoffs.append(spars[False])
            fig1.loc[mseustr, 'spars'] = spars[False], spars[True] # save
            if mvimseu2exmpli.get(mseustr) == fig1exmpli:
                exmplis.append(keptmseui)
                exmplmseustrs.append(mseustr)
            else:
                normlis.append(keptmseui)
            keptmseui += 1 # manually increment
        sparsons, sparsoffs = np.asarray(sparsons), np.asarray(sparsoffs)
        # plot y=x line:
        xyline = [linmin, linmax], [linmin, linmax]
        a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
        # plot normal (non-example) points:
        a.scatter(sparsons[normlis], sparsoffs[normlis], clip_on=False,
                  marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
        # plot example points, one at a time:
        for exmpli, mseustr in zip(exmplis, exmplmseustrs):
            marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
            c = exmpli2clr[mvimseu2exmpli[mseustr]]
            sz = exmpli2sz[mvimseu2exmpli[mseustr]]
            lw = exmpli2lw[mvimseu2exmpli[mseustr]]
            a.scatter(sparsons[exmpli], sparsoffs[exmpli], clip_on=False,
                      marker=marker, c=c, s=sz, lw=lw)
        a.set_xlabel('Suppression spars') # keep it short to maximize space for axes
        a.set_ylabel('Feedback spars')
        a.set_xlim(linmin, linmax)
        a.set_ylim(linmin, linmax)
        a.set_xticks(ticks)
        a.set_yticks(ticks)
        a.set_aspect('equal')
        a.spines['left'].set_position(('outward', 4))
        a.spines['bottom'].set_position(('outward', 4))
        #t, p = ttest_rel(sparsons, sparsoffs) # paired t-test
        #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))


# scatter plot movie PSTH reliability:
figsize = DEFAULTFIGURESIZE[0]*0.98, DEFAULTFIGURESIZE[1] # tweak to match others
linmin, linmax = 0, 0.52
ticks = 0, 0.5
for kind in ['nat']:#MVIKINDS:
    for st8 in ['none']:#ALLST8S:
        f, a = plt.subplots(figsize=figsize)
        wintitle('opto reliability %s %s' % (kind, st8))
        relons, reloffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in mvimseustrs:
            rel = mviresp.loc[mseustr, kind, st8]['rel']
            snr = mviresp.loc[mseustr, kind, st8]['snr']
            if rel.isna().any(): # missing for at least one opto condition
                continue
            if (snr < SNRTHRESH).all(): # neither condition has decent SNR
                continue
            relons.append(rel[True])
            reloffs.append(rel[False])
            # collect rhos from all trial pairs, reduce to single rho per trial, use
            # to estimate significance of individual units:
            rhos = mviresp.loc[mseustr, kind, st8]['rhos']
            fullrhos, meanrhos = {}, {}
            for opto in OPTOS:
                npairs = len(rhos[opto])
                # from quadratic formula solution to npairs = ntrials*(ntrials-1) / 2
                ntrials = (1 + np.sqrt(1 + 8*npairs)) / 2
                assert np.isclose(ntrials, int(ntrials)) # make sure this float is int
                ntrials = int(ntrials)
                uti = np.triu_indices(ntrials, k=1) # upper indices, excludes diagonal
                lti = np.tril_indices(ntrials, k=-1) # lower indices, excludes diagonal
                # empty full corr matrix, nans will remain on diagonal and will be ignored:
                fullrhos[opto] = np.tile(np.nan, (ntrials, ntrials))
                fullrhos[opto][uti] = rhos[opto] # fill the upper triangle
                # copy upper to lower triangle,
                # see https://stackoverflow.com/a/42209263/2020363:
                fullrhos[opto][lti] = fullrhos[opto].T[lti]
                meanrhos[opto] = np.nanmean(fullrhos[opto], axis=1) # along either axis is fine
            _, pval = ttest_ind(meanrhos[False], meanrhos[True], equal_var=False)
            sgnfs.append(pval < SCATTERPTHRESH) # bool
            fig1.loc[mseustr, 'rel'] = rel[False], rel[True] # save
            if mvimseu2exmpli.get(mseustr) == fig1exmpli:
                exmplis.append(keptmseui)
                exmplmseustrs.append(mseustr)
            else:
                normlis.append(keptmseui)
            keptmseui += 1 # manually increment
        relons, reloffs = np.asarray(relons), np.asarray(reloffs)
        sgnfs = np.asarray(sgnfs)
        # plot y=x line:
        xyline = [linmin, linmax], [linmin, linmax]
        a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
        insgnfis, = np.where(sgnfs == False)
        sgnfis, = np.where(sgnfs == True)
        normlinsgnfis = np.intersect1d(normlis, insgnfis)
        normlsgnfis = np.intersect1d(normlis, sgnfis)
        if len(normlinsgnfis) > 0:
            # plot normal (non-example) insignificant points:
            c = desat(st82clr[st8], SGNF2ALPHA[False]) # do manual alpha mixing
            a.scatter(relons[normlinsgnfis], reloffs[normlinsgnfis], clip_on=False,
                      marker='.', c='None', edgecolor=c, s=DEFSZ)
        if len(normlsgnfis) > 0:
            # plot normal (non-example) significant points:
            c = desat(st82clr[st8], SGNF2ALPHA[True]) # do manual alpha mixing
            a.scatter(relons[normlsgnfis], reloffs[normlsgnfis], clip_on=False,
                      marker='.', c='None', edgecolor=c, s=DEFSZ)
        # plot example points, one at a time:
        for exmpli, mseustr in zip(exmplis, exmplmseustrs):
            marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
            c = exmpli2clr[mvimseu2exmpli[mseustr]]
            sz = exmpli2sz[mvimseu2exmpli[mseustr]]
            lw = exmpli2lw[mvimseu2exmpli[mseustr]]
            a.scatter(relons[exmpli], reloffs[exmpli], clip_on=False,
                      marker=marker, c=c, s=sz, lw=lw)
        a.set_xlabel('Suppression rel') # keep it short to maximize space for axes
        a.set_ylabel('Feedback rel')
        a.set_xlim(linmin, linmax)
        a.set_ylim(linmin, linmax)
        a.set_xticks(ticks)
        a.set_yticks(ticks)
        a.set_aspect('equal')
        a.spines['left'].set_position(('outward', 4))
        a.spines['bottom'].set_position(('outward', 4))
        #t, p = ttest_rel(relons, reloffs, nan_policy='omit') # paired t-test
        #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))


# plot movie SNR PDFs, show SNRTHRESH:
figsize = DEFAULTFIGURESIZE
bins = np.arange(0, 0.5, 0.005)
#bins = np.logspace(-2, -0.2, 50, base=10)
for kind in ['nat']:#MVIKINDS:
    for st8 in ['none']:#ALLST8S:
        for opto in [None, False, True]:
            if opto is not None:
                snrs = mviresp['snr'][(slice(None), kind, st8, opto)].dropna() # one val per mseu
            else:
                snrs = mviresp['snr'][(slice(None), kind, st8)].dropna() # two vals per mseu?
            snrs = np.float64(snrs) # convert from object array
            f, a = plt.subplots(figsize=figsize)
            titlestr = 'SNR PDF %s %s' % (kind, st8)
            if opto is not None:
                titlestr += ' %s' % opto2fb[opto]
            wintitle(titlestr, f)
            c = st82clr[st8]
            a.hist(snrs, bins=bins, density=True, color=c,
                   alpha=opto2alpha[opto], label=opto2fb[opto])
            a.axvline(x=SNRTHRESH, ls='-', marker='', c='gray')#, zorder=-np.inf)
            txt = 'n=%d' % len(snrs)
            a.add_artist(AnchoredText(txt, loc='upper right', frameon=False,
                                      prop=dict(color=c, alpha=opto2alpha[opto])))
            #a.set_xscale('log')
            #l = a.legend(frameon=False)
            #l.set_draggable(True)
            a.set_xlabel('SNR')
            if opto is not None:
                a.set_xlabel('%s SNR' % opto2fb[opto].title())
            a.set_ylabel('Probability density')


# scatter plot movie SNR:
figsize = DEFAULTFIGURESIZE
linmin, linmax = 0, 0.6
ticks = 0, 0.5
for kind in ['nat']:#MVIKINDS:
    for st8 in ['none']:#ALLST8S:
        f, a = plt.subplots(figsize=figsize)
        wintitle('opto SNR %s %s' % (kind, st8))
        snrons, snroffs, exmplis, exmplmseustrs, normlis = [], [], [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in mvimseustrs:
            snr = mviresp.loc[mseustr, kind, st8]['snr']
            if snr.isna().any(): # missing for at least one opto condition
                continue
            if (snr < SNRTHRESH).all(): # neither condition has decent SNR
                continue
            snrons.append(snr[True])
            snroffs.append(snr[False])
            fig1.loc[mseustr, 'snr'] = snr[False], snr[True] # save
            if mvimseu2exmpli.get(mseustr) == fig1exmpli:
                exmplis.append(keptmseui)
                exmplmseustrs.append(mseustr)
            else:
                normlis.append(keptmseui)
            keptmseui += 1 # manually increment
        snrons, snroffs = np.asarray(snrons), np.asarray(snroffs)
        # plot y=x line:
        xyline = [linmin, linmax], [linmin, linmax]
        a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
        # plot normal (non-example) points:
        a.scatter(snrons[normlis], snroffs[normlis], clip_on=False,
                  marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
        # plot example points, one at a time:
        for exmpli, mseustr in zip(exmplis, exmplmseustrs):
            marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
            c = exmpli2clr[mvimseu2exmpli[mseustr]]
            sz = exmpli2sz[mvimseu2exmpli[mseustr]]
            lw = exmpli2lw[mvimseu2exmpli[mseustr]]
            a.scatter(snrons[exmpli], snroffs[exmpli], clip_on=False,
                      marker=marker, c=c, s=sz, lw=lw)
        a.set_xlabel('Suppression SNR')
        a.set_ylabel('Feedback SNR')
        a.set_xlim(linmin, linmax)
        a.set_ylim(linmin, linmax)
        a.set_xticks(ticks)
        a.set_yticks(ticks)
        a.set_aspect('equal')
        a.spines['left'].set_position(('outward', 4))
        a.spines['bottom'].set_position(('outward', 4))
        #t, p = ttest_rel(snrons, snroffs, nan_policy='omit') # paired t-test
        #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))


# plot movie PSTH peak width distributions:
## TODO: these width measures should maybe be normalized separately for each unit,
## so that the effect of feedback suppression isn't drowned out by the variability of
## PSTH peak widths across units
normalizepeakwidths = True
figsize = DEFAULTFIGURESIZE
for kind in ['nat']:#MVIKINDS:
    for st8 in ['none']:#ALLST8S:
        f, a = plt.subplots(figsize=figsize)
        titlestr = ('opto peak width PDF %s %s PSTHBASELINEMEDIANX=%d PSTHPEAKMINTHRESH=%d'
                    % (kind, st8, PSTHBASELINEMEDIANX, PSTHPEAKMINTHRESH))
        if normalizepeakwidths:
            titlestr += " normed"
        wintitle(titlestr, f)
        pkwons, pkwoffs = [], []
        for mseustr in mvimseustrs:
            pkws = mviresp.loc[mseustr, kind, st8]['pkws']
            snr = mviresp.loc[mseustr, kind, st8]['snr']
            if pkws.isna().any():
                continue # skip this mseu
            if (snr < SNRTHRESH).all(): # neither condition has decent SNR
                continue
            if normalizepeakwidths:
                maxpkw = np.concatenate(pkws.values).max() # max peak width for this mseu
                pkws = pkws / maxpkw # normalize all peak widths by maxpkw
            pkwons.append(pkws[True])
            pkwoffs.append(pkws[False])
        c = st82clr[st8]
        pkwons, pkwoffs = np.concatenate(pkwons), np.concatenate(pkwoffs)
        for opto, pkws, loc in zip(OPTOS, [pkwoffs, pkwons], ['upper right', 'upper left']):
            a.hist(pkws, bins=50, density=True, histtype='step', color=c,
                   alpha=opto2alpha[opto], label=opto2fb[opto])
            txt = 'n=%d' % len(pkws)
            a.add_artist(AnchoredText(txt, loc=loc, frameon=False,
                                      prop=dict(color=c, alpha=opto2alpha[opto])))
        a.set_xscale('log')
        #a.set_yscale('log')
        axes_disable_scientific(a)
        #l = a.legend(frameon=False)
        #l.set_draggable(True)
        a.set_xlabel('PSTH peak width (s)')
        a.set_ylabel('Probability density')
        #t, p = ttest_rel(pkwoffs, pkwons, nan_policy='omit') # paired t-test
        #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))


# scatter plot movie PSTH mean peak widths:
figsize = DEFAULTFIGURESIZE[0]*1.05, DEFAULTFIGURESIZE[1] # tweak to make space for log units
linmin, linmax = 0.025, 0.2
ticks = 0.05, 0.1, 0.2
for kind in ['nat']:#MVIKINDS:
    for st8 in ['none']:#ALLST8S:
        f, a = plt.subplots(figsize=figsize)
        wintitle('opto mean peak width %s %s' % (kind, st8))
        meanpkwons, meanpkwoffs, exmplis, exmplmseustrs, normlis = [], [], [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in mvimseustrs:
            meanpkw = mviresp.loc[mseustr, kind, st8]['meanpkw']
            snr = mviresp.loc[mseustr, kind, st8]['snr']
            if meanpkw.isna().any(): # missing for at least one opto condition
                continue
            if (snr < SNRTHRESH).all(): # neither condition has decent SNR
                continue
            meanpkwons.append(meanpkw[True])
            meanpkwoffs.append(meanpkw[False])
            fig1.loc[mseustr, 'meanpkw'] = meanpkw[False], meanpkw[True] # save
            if mvimseu2exmpli.get(mseustr) == fig1exmpli:
                exmplis.append(keptmseui)
                exmplmseustrs.append(mseustr)
            else:
                normlis.append(keptmseui)
            keptmseui += 1 # manually increment
        meanpkwons, meanpkwoffs = np.asarray(meanpkwons), np.asarray(meanpkwoffs)
        # plot y=x line:
        xyline = [linmin, linmax], [linmin, linmax]
        a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
        # plot normal (non-example) points:
        a.scatter(meanpkwons[normlis], meanpkwoffs[normlis], clip_on=False,
                  marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
        # plot example points, one at a time:
        for exmpli, mseustr in zip(exmplis, exmplmseustrs):
            marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
            c = exmpli2clr[mvimseu2exmpli[mseustr]]
            sz = exmpli2sz[mvimseu2exmpli[mseustr]]
            lw = exmpli2lw[mvimseu2exmpli[mseustr]]
            a.scatter(meanpkwons[exmpli], meanpkwoffs[exmpli], clip_on=False,
                      marker=marker, c=c, s=sz, lw=lw)
        a.set_xscale('log', basex=2)
        a.set_yscale('log', basey=2)
        a.set_xlabel('Suppression pkw (s)') # keep it short to maximize space for axes
        a.set_ylabel('Feedback pkw (s)')
        a.set_xlim(linmin, linmax)
        a.set_ylim(linmin, linmax)
        a.set_xticks(ticks)
        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')
        a.spines['left'].set_position(('outward', 4))
        a.spines['bottom'].set_position(('outward', 4))
        a.xaxis.set_major_formatter(ScalarFormatter())
        a.yaxis.set_major_formatter(ScalarFormatter())
        a.xaxis.set_major_formatter(FormatStrFormatter("%.2g"))
        a.yaxis.set_major_formatter(FormatStrFormatter("%.2g"))
        #t, p = ttest_rel(meanpkwoffs, meanpkwons, nan_policy='omit') # paired t-test
        #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))


"""Figure 1S2 FMI plots"""

# scatter plot FMI of all movie measures vs meanrate FMI:
figsize = DEFAULTFIGURESIZE
linmin, linmax, linstep = -1, 1, 1
ticks = np.arange(linmin, linmax+linstep, linstep)
for kind in ['nat']:#MVIKINDS:
    for st8 in ['none']:#ALLST8S:
        for measure in modmeasuresnoblank:
            if measure.startswith('meanrate'):
                continue # don't bother plotting meanrate* against meanrate
            axislabel = measure2axislabel[measure]
            rfmis, msrfmis = [], [] # rate FMIs, measure FMIs
            exmplis, exmplmseustrs, normlis = [], [], []
            keptmseui = 0 # manually init and increment instead of using enumerate()
            for mseustr in mvimseustrs:
                rfmi = mviFMI.loc[mseustr, st8]['meanrate']
                msrfmi = mviFMI.loc[mseustr, st8][measure]
                if pd.isna(rfmi) or pd.isna(msrfmi): # missing at least one FMI value
                    continue
                rfmis.append(rfmi)
                msrfmis.append(msrfmi)
                if mvimseu2exmpli.get(mseustr) == fig1exmpli:
                    exmplis.append(keptmseui)
                    exmplmseustrs.append(mseustr)
                else:
                    normlis.append(keptmseui)
                keptmseui += 1 # manually increment
            rfmis, msrfmis = np.asarray(rfmis), np.asarray(msrfmis)
            ## scatter plot FMIs:
            f, a = plt.subplots(figsize=figsize)
            wintitle('opto FMI %s vs FMI rate %s %s' % (measure, kind, st8))
            # plot x=0 and y=0 lines:
            a.axhline(y=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
            a.axvline(x=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
            # plot normal (non-example) points:
            a.scatter(rfmis[normlis], msrfmis[normlis], clip_on=False,
                      marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
            # plot example points, one at a time:
            for exmpli, mseustr in zip(exmplis, exmplmseustrs):
                marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
                c = exmpli2clr[mvimseu2exmpli[mseustr]]
                sz = exmpli2sz[mvimseu2exmpli[mseustr]]
                lw = exmpli2lw[mvimseu2exmpli[mseustr]]
                a.scatter(rfmis[exmpli], msrfmis[exmpli], clip_on=False,
                          marker=marker, c=c, s=sz, lw=lw)
            # get fname of appropriate LMM .cvs file:
            if measure == 'meanburstratio': # use Steffen's LMM linregress fit, for fig1S2
                fname = os.path.join(PAPERPATH, 'stats', 'figure_1S2c_coefs.csv')
            elif measure == 'spars': # use Steffen's LMM linregress fit, for fig1S2
                fname = os.path.join(PAPERPATH, 'stats', 'figure_1S2d_coefs.csv')
            elif measure == 'rel': # use Steffen's LMM linregress fit, for fig1S2
                fname = os.path.join(PAPERPATH, 'stats', 'figure_1S2e_coefs.csv')
            elif measure == 'snr': # use Steffen's LMM linregress fit, for fig1S2
                fname = os.path.join(PAPERPATH, 'stats', 'figure_1S2f_coefs.csv')
            elif measure == 'meanpkw': # use Steffen's LMM linregress fit, for fig1S2
                fname = os.path.join(PAPERPATH, 'stats', 'figure_1S2g_coefs.csv')
            else:
                print("WARNING: No LMM stats for %s measure" % measure)
                fname = None
            # fetch LMM linregress fit params from .csv:
            if fname:
                try:
                    df = pd.read_csv(fname)
                    foundregression = True
                except FileNotFoundError:
                    print('Missing file: %s' % fname)
                    foundregression = False
                if foundregression:
                    mm = df['slope'][0]
                    b = df['intercept'][0]
                    x = np.array([rfmis.min(), rfmis.max()])
                    y = mm * x + b
                    a.plot(x, y, '-', color='red') # plot linregress fit
            a.set_xlabel('FR FMI')
            if axislabel.islower():
                axislabel = axislabel.capitalize()
            ylabel = '%s FMI' % axislabel
            a.set_ylabel(ylabel)
            a.set_xlim(linmin, linmax)
            a.set_ylim(linmin, linmax)
            a.set_xticks(ticks)
            a.set_yticks(ticks)
            a.set_aspect('equal')
            a.spines['left'].set_position(('outward', 4))
            a.spines['bottom'].set_position(('outward', 4))


"""Some figure 4S1 plots"""

# scatter plot meanrates of first 2 and last 2 s of movies, limited to first 120 trials
# per condition, as per standard grating experiments:
figsize = DEFAULTFIGURESIZE[0]*1.04, DEFAULTFIGURESIZE[1] # tweak to make space for log units
logmin, logmax = -2, 2
logticks = np.array([-2, 0, 2])
#log0min = logmin + 0.05
kind = 'nat'
st8 = 'none'
for meanratestr, rangestr in zip(['meanrate02', 'meanrate35'], ['0-2', '3-5']):
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto meanrate movie t=%s s %d trials %s %s' % (rangestr, NTRIALSMVIGRT, kind, st8))
    rons, roffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], []
    keptmseui = 0 # manually init and increment instead of using enumerate()
    for mseustr in mvimseustrs:
        meanrate = mviresp.loc[mseustr, kind, st8][meanratestr]
        snr = mviresp.loc[mseustr, kind, st8]['snr']
        if meanrate.isna().any(): # missing one or both meanrates
            continue
        if (snr < SNRTHRESH).all(): # neither condition has decent SNR on the full range
            continue
        rons.append(meanrate[True])
        roffs.append(meanrate[False])
        ratesstr = meanratestr.split('mean')[1]+'s'
        rates = mviresp.loc[mseustr, kind, st8][ratesstr]
        _, pval = ttest_ind(rates[False], rates[True], equal_var=False)
        sgnfs.append(pval < SCATTERPTHRESH) # bool
        fig1.loc[mseustr, meanratestr] = meanrate[False], meanrate[True] # save
        fig1.loc[mseustr][ratesstr] = rates # save trial-wise values
        if mvimseu2exmpli.get(mseustr) == fig1exmpli:
            exmplis.append(keptmseui)
            exmplmseustrs.append(mseustr)
        else:
            normlis.append(keptmseui)
        keptmseui += 1 # manually increment
    rons = np.asarray(rons)
    roffs = np.asarray(roffs)
    sgnfs = np.asarray(sgnfs)
    # 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)
    insgnfis, = np.where(sgnfs == False)
    sgnfis, = np.where(sgnfs == True)
    normlinsgnfis = np.intersect1d(normlis, insgnfis)
    normlsgnfis = np.intersect1d(normlis, sgnfis)
    # plot normal (non-example) insignificant points:
    c = desat(st82clr[st8], SGNF2ALPHA[False]) # do manual alpha mixing
    a.scatter(rons[normlinsgnfis], roffs[normlinsgnfis], clip_on=False,
              marker='.', c='None', edgecolor=c, s=DEFSZ)
    # plot normal (non-example) significant points:
    c = desat(st82clr[st8], SGNF2ALPHA[True]) # do manual alpha mixing
    a.scatter(rons[normlsgnfis], roffs[normlsgnfis], clip_on=False,
              marker='.', c='None', edgecolor=c, s=DEFSZ)
    # plot normal (non-example) points:
    #a.scatter(rons[normlis], roffs[normlis], clip_on=False,
    #          marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
    # plot example points, one at a time:
    for exmpli, mseustr in zip(exmplis, exmplmseustrs):
        marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
        c = exmpli2clr[mvimseu2exmpli[mseustr]]
        sz = exmpli2sz[mvimseu2exmpli[mseustr]]
        lw = exmpli2lw[mvimseu2exmpli[mseustr]]
        a.scatter(rons[exmpli], roffs[exmpli], marker=marker, clip_on=False,
                  c=c, s=sz, lw=lw)
    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')
    a.spines['left'].set_position(('outward', 4))
    a.spines['bottom'].set_position(('outward', 4))
    #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))