natfeedback_code_eLife/fig3.py

"""Figure 3 plots, use run -i fig3.py"""

mi = pd.MultiIndex.from_product([grtmseustrs, OPTOS], names=['mseu', 'opto'])
columns = ['trialis', 'rates', 'blankrates', 'blankcondrates', # single trials
           'burstratios', 'blankburstratios', 'blankcondburstratios', # single trials
           'meanrate', 'blankmeanrate', 'blankcondmeanrate',
           'meanburstratio', 'blankmeanburstratio', 'blankcondmeanburstratio',
           'osi', 'oripref', 'f1f0',
           'phi', 'bphi', 'nbphi', 'dphi', 'dbphi', 'dnbphi']
fig3 = pd.DataFrame(index=mi, columns=columns) # excludes kind and st8
## NOTE: make sure that rows are not overwritten due to temporarily enabling iteration
## over st8 in any of the following loops!


# plot grating rasters for each unit by state and opto condition:
showbursts = True
## TODO: fix hard-coding:
if EXPTYPE == 'negntsrmvis':
    NTRIALSPERCONDITION = 8 # per stimi and opto combination
    inchespertrial = 1/20*RASTERSCALEX
else:
    NTRIALSPERCONDITION = 10 # per stimi and opto combination
    inchespertrial = 1/25*RASTERSCALEX
#outliers = ['PVCre_2018_0001_s02_e02_u36', 'PVCre_2018_0001_s02_e06_u36',
#            'PVCre_2018_0001_s02_e02_u64', 'PVCre_2018_0001_s02_e06_u64']
for mseustr in []:#grtmseu2exmpli: #grtmseustrs: #outliers:
    #['Ntsr1Cre_2019_0008_s06_e02_u46']
    #['Ntsr1Cre_2020_0001_s04_e04_u12'] # negntsr fig
    for st8 in ['none']:#ALLST8S:
        for opto in OPTOS:
            dt = grtresp.loc[mseustr, st8, opto]['dt']
            wdt = -OFFSETS[0] + dt + OFFSETS[1] # wide raster duration, sec
            tmin, tmax = 0 + OFFSETS[0], dt + OFFSETS[1]
            optotrange = grtresp.loc[mseustr, st8, opto]['optotrange']
            wraster = grtresp.loc[mseustr, st8, opto]['wraster']
            if np.any(pd.isna(wraster)): # no raster for this condition
                continue
            ntrials = len(wraster)
            title = '%s grating %s %s raster' % (mseustr, st8, opto)
            burstis = None
            if showbursts:
                wburstis = grtresp.loc[mseustr, st8, opto]['wburstis']
            a = simpletraster(raster=wraster, alpha=0.5, s=1, dt=dt, offsets=OFFSETS,
                              scatter=True, # False: low quality, but fast
                              title=title, inchespersec=1.5*RASTERSCALEX,
                              inchespertrial=inchespertrial,
                              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.02
            #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)
            #a.set_xticks(range(0, intround(dt)+1)) # force integer 1 sec tick marks
            # force ticks on trial axis separating ori conditions:
            a.set_yticks(range(0, ntrials+1, NTRIALSPERCONDITION))


# Plot grating tuning curves for example units. Requires a database connection.
# Uses effectively same grey as 0.4 alpha for opto True condition:
try:
    tun
except:
    print("WARNING: Can't plot tuning curves without database connection")
    tun = None
if tun is not None:
    for mseustr in grtmseu2exmpli: #grtmseustrs:#['Ntsr1Cre_2020_0001_s04_e04_u55']
        tunrow = (tun & {**mseustr2key(mseustr), 'st8_crit':'none'})
        if tunrow:
            tunrow.plot(cs=['k', optoblue])


## TODO: plot f1/f0 vs. ori for some example units, e.g. Ntsr1Cre_2020_0001_s04_e04_u55
## Some of the time, or maybe even much of the time, f1/f0 vs. ori might show better
## orientation tuning than firing rate vs. ori.


## calculate cycle averages (i.e. grating PSTHs), F1/F0, and phi for all units:
## TODO: fit sinusoids and use fit error to evaluate how linear the responses are?
PLOTCYCPSTH = False
CYCPSTHTHRESH = 10 # minimum cycle average peak, Hz
F1F0THRESH = 0.25 # minimum grating modulation index
# minimum ratio of suppression burst cycle average peak to suppression of all spikes
# cycle average peak, at maxori:
BURSTRATIOTHRESH = 0.1
TESTPEAKS = 'both' # test peak amplitude of 'either' or 'both' cycle PSTH opto conditions
ncyclesoffset = 0.25 # fraction of previous and next cycle to plot
figsize = DEFAULTFIGURESIZE
for st8 in ['none']:#ALLST8S:
    color = st82clr[st8]
    f1f0offs, f1f0ons = [], []
    phioffs, phions, bphioffs, bphions, nbphioffs, nbphions = [], [], [], [], [], []
    dphis, dbphis, dnbphis = [], [], []
    exmplis, exmplmseustrs, normlis = [], [], []
    keptmseui = 0 # manually init and increment instead of using enumerate()
    for mseustr in grtmseustrs: #['Ntsr1Cre_2020_0001_s04_e04_u12']
        cycpsths, cycbpsths, cycts = {}, {}, {} # indexed by ori and opto
        f1f0s, phis = {}, {} # indexed by ori and opto
        rasters = grtresp.loc[mseustr, st8]['raster']
        brasters = grtresp.loc[mseustr, st8]['braster']
        if rasters.isna().any():
            print('%s: missing one or both raster opto conditions, skipping' % mseustr)
            continue
        if brasters.isna().any():
            print('%s: missing one or both braster opto conditions, skipping' % mseustr)
            continue
        orisopto = grtresp.loc[mseustr, st8]['oris'] # sorted, correspond to rows in raster
        ## TODO: maybe instead of enforcing this, just pull out all existing oris in both
        ## opto conditions, and do unique on those, and if any don't exist then handle that
        ## in the innermost loop?
        assert np.equal(*orisopto).all() # make sure oris are identical across opto
        # get oris from arbitrary condition, convert to int for use as dict keys:
        uoris = np.int64(np.unique(orisopto[False]))
        for ori in uoris:
            cycpsths[ori], cycbpsths[ori], cycts[ori] = {}, {}, {}
            f1f0s[ori], phis[ori] = {}, {}
            for opto in OPTOS:
                raster = rasters[opto]
                braster = brasters[opto]
                oris = orisopto[opto]
                dt = grtresp.loc[mseustr, st8, opto]['dt']
                tfreq = grtresp.loc[mseustr, st8, opto]['tfreq']
                cycdt = 1 / tfreq # cycle duration, s
                t0 = 0 + cycdt # exclude first cycle...
                t1 = dt # ...include up to end of last cycle
                offsets = np.array([-cycdt, cycdt]) * ncyclesoffset
                # to avoid edge effects of Guassian smoothing in raster2psth(), give it a wide
                # set of bins, and then slice out a slightly narrower (one bin width less at
                # each end) set of bins:
                wtrange = 0+offsets[0], cycdt+offsets[1] # wide (full) trange
                ntrange = wtrange[0]+binw, wtrange[1]-binw # narrow trange
                wcycbins = split_tranges([wtrange], binw, tres)
                ncycbins = split_tranges([ntrange], binw, tres)
                # find which rows in wcycbins correpond to those in ncycbins, based on the
                # start time of each bin:
                nbinis = wcycbins[:, 0].searchsorted(ncycbins[:, 0])
                orirowis = oris == ori
                rast = raster[orirowis] # raster restricted to rows corrsponding to ori
                brast = braster[orirowis]
                cycraster = wrap_raster(rast, t0, t1, cycdt, offsets=offsets)
                cycbraster = wrap_raster(brast, t0, t1, cycdt, offsets=offsets)
                f0, _ = raster2freqcomp(rast, dt, 0, mean='vector')
                f1, phi = raster2freqcomp(rast, dt, tfreq, mean='vector')
                if f0 == 0.0: # no spikes at all for this ori and opto combination
                    f1f0s[ori][opto] = np.nan
                    phis[ori][opto] = np.nan
                else:
                    f1f0s[ori][opto] = f1 / f0
                    phis[ori][opto] = np.angle(np.exp(1j*phi)) + pi # wrapped +ve angle, rad
                cycts[ori][opto] = ncycbins.mean(axis=1) # save narrow raster timepoints
                wrast = raster2psth(cycraster, wcycbins, binw, tres, kernel) # wide raster
                wbrast = raster2psth(cycbraster, wcycbins, binw, tres, kernel) # wide burst raster
                cycpsths[ori][opto] = wrast[nbinis] # slice down to narrow raster
                cycbpsths[ori][opto] = wbrast[nbinis]
        # find ori with max cycle average peak, across opto conditions:
        maxcycpsthvals = []
        for ori in uoris:
            maxcycpsthvals.append(max(cycpsths[ori][False].max(), cycpsths[ori][True].max()))
        maxcycpsthvali = np.argmax(maxcycpsthvals)
        maxori = uoris[maxcycpsthvali]
        if TESTPEAKS == 'either':
            # test for minimum PSTH peak height in at least one opto condition:
            maxcycpsthval = maxcycpsthvals[maxcycpsthvali]
            if maxcycpsthval < CYCPSTHTHRESH:
                continue # skip this mseu
        elif TESTPEAKS == 'both':
            # test for minimum PSTH peak height in both opto conditions:
            maxcycpsthoff = cycpsths[maxori][False].max()
            maxcycpsthon = cycpsths[maxori][True].max()
            if (np.array([maxcycpsthoff, maxcycpsthon]) < CYCPSTHTHRESH).any():
                continue # skip this mseu
        # test for minimum f1/f0 in at least one opto condition:
        maxf1f0atmaxori = max(f1f0s[maxori][False], f1f0s[maxori][True])
        if maxf1f0atmaxori < F1F0THRESH:
            continue # skip this mseu
        # collect f1f0 and phi values at maxori:
        f1f0off, f1f0on = f1f0s[maxori][False], f1f0s[maxori][True]
        f1f0offs.append(f1f0off)
        f1f0ons.append(f1f0on)
        phioff, phion = phis[maxori][False], phis[maxori][True]
        if EXPTYPE == 'pvmvis':
            # wrap some points:
            if (phioff < 60*pi/180) & (phion > 240*pi/180):
                phioff += 360*pi/180 # wrap point up over the y=x line
            elif (phioff > 340*pi/180) & (phion < 15*pi/180):
                phioff -= 360*pi/180 # wrap point down over the y=x line
        elif EXPTYPE == 'ntsrmvis':
            # wrap one point:
            if (phioff < 30*pi/180) & (phion > 340*pi/180):
                phioff += 360*pi/180 # wrap point up over the y=x line
        phioffs.append(phioff)
        phions.append(phion)
        fig3.loc[mseustr, 'phi'] = phioff*180/pi, phion*180/pi # save
        dphi = np.angle(np.exp(1j*phioff) / np.exp(1j*phion)) # phioff - phion
        dphis.append(dphi)
        if dphi < (-75*pi/180):
            print('OUTLIER: %s dphi=%g' % (mseustr, dphi))
        fig3.loc[mseustr, 'dphi'] = dphi*180/pi, dphi*180/pi # save, ignore opto field
        # collect phi of mseus that meet minimum ratio of suppression burst cycle average peak
        # to suppression all spikes cycle average peak, at maxori:
        cycbpsthmax, cycpsthmax = cycbpsths[maxori][True].max(), cycpsths[maxori][True].max()
        if cycpsthmax == 0.0:
            if cycbpsthmax == 0.0: # numerator and denominator are both 0
                suppburstratioatmaxori = 0.0
            else:
                suppburstratioatmaxori = np.inf
        else:
            suppburstratioatmaxori = cycbpsthmax / cycpsthmax
        # if suppression burst cycle average peak is big enough, treat is as a burst mseu,
        # otherwise treat it as a nonburst mseu:
        if suppburstratioatmaxori >= BURSTRATIOTHRESH:
            bphioffs.append(phioff)
            bphions.append(phion)
            fig3.loc[mseustr, 'bphi'] = phioff*180/pi, phion*180/pi # save
            dbphi = np.angle(np.exp(1j*phioff) / np.exp(1j*phion)) # bphioff - bphion
            dbphis.append(dbphi)
            fig3.loc[mseustr, 'dbphi'] = dbphi*180/pi, dbphi*180/pi # save, ignore opto field
        else:
            nbphioffs.append(phioff)
            nbphions.append(phion)
            fig3.loc[mseustr, 'nbphi'] = phioff*180/pi, phion*180/pi # save
            dnbphi = np.angle(np.exp(1j*phioff) / np.exp(1j*phion)) # nbphioff - nbphion
            dnbphis.append(dnbphi)
            fig3.loc[mseustr, 'dnbphi'] = dnbphi*180/pi, dnbphi*180/pi # save, ignore opto field
        ## plot maxori cycle average:
        if PLOTCYCPSTH:# and mseustr in grtmseu2exmpli:#['PVCre_2018_0003_s03_e02_u51']
            f, a = plt.subplots(figsize=(2, figsize[1]))
            title = ('%s %s PSTH maxori=%d offset=%g'
                     % (mseustr, st8, maxori, ncyclesoffset))
            wintitle(title)
            for opto in OPTOS:
                c = desat(color, opto2alpha[opto]) # do manual alpha mixing
                bc = desat(burstclr, opto2alpha[opto]) # do manual alpha mixing
                fb = opto2fb[opto].title()
                if st8 == 'none':
                    label = fb
                else:
                    label = ', '.join([fb, st8.title()])
                cyct = cycts[maxori][opto]
                cycpsth = cycpsths[maxori][opto]
                cycbpsth = cycbpsths[maxori][opto]
                a.plot(cyct, cycpsth, '-', color=c, label=label)
                a.plot(cyct, cycbpsth, '-', color=bc, label=label+', Burst')
                a.axvline(x=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
                a.axvline(x=cycdt, ls='--', marker='', color='lightgray', zorder=-np.inf)
                #l = a.legend(frameon=False)
                #l.set_draggable(True)
                a.set_xlabel('Time (s)')
                a.set_ylabel('Firing rate (spk/s)')
                a.set_xlim(cyct[0], cyct[-1])
        if mseustr in grtmseu2exmpli:
            exmplis.append(keptmseui)
            exmplmseustrs.append(mseustr)
        else:
            normlis.append(keptmseui)
        keptmseui += 1 # manually increment
    f1f0offs, f1f0ons = np.asarray(f1f0offs), np.asarray(f1f0ons)
    phioffs, phions = np.asarray(phioffs), np.asarray(phions)
    bphioffs, bphions = np.asarray(bphioffs), np.asarray(bphions)
    nbphioffs, nbphions = np.asarray(nbphioffs), np.asarray(nbphions)
    dphis, dbphis, dnbphis = np.asarray(dphis), np.asarray(dbphis), np.asarray(dnbphis)
    ## scatter plot collected maxori f1f0s:
    figsize = DEFAULTFIGURESIZE
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto f1f0 maxori grating %s CYCPSTHTHRESH=%g F1F0THRESH=%g '
             'TESTPEAKS=%s' % (st8, CYCPSTHTHRESH, F1F0THRESH, TESTPEAKS))
    # plot y=x line:
    linmin, linmax = 0, 2
    xyline = [linmin, linmax], [linmin, linmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot normal (non-example) points:
    a.scatter(f1f0ons[normlis], f1f0offs[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[grtmseu2exmpli[mseustr]]
        c = exmpli2clr[grtmseu2exmpli[mseustr]]
        sz = exmpli2sz[grtmseu2exmpli[mseustr]]
        lw = exmpli2lw[grtmseu2exmpli[mseustr]]
        a.scatter(f1f0ons[exmpli], f1f0offs[exmpli], clip_on=False,
                  marker=marker, c=c, s=sz, lw=lw)
    a.set_xlabel('Suppression F1/F0')
    a.set_ylabel('Feedback F1/F0')
    #a.set_xscale('log', basex=2)
    #a.set_yscale('log', basey=2)
    #a.set_xlim(0.125, 2)
    #a.set_ylim(0.125, 2)
    #a.set_xticks([0.125, 0.25, 0.5, 1, 2])
    #a.set_yticks(a.get_xticks()) # make scale y ticks the same as x ticks
    #axes_disable_scientific(a)
    a.set_xlim(linmin, linmax)
    a.set_ylim(linmin, linmax)
    a.set_xticks(np.arange(linmin, linmax+1))
    a.set_yticks(a.get_xticks()) # make scale y ticks the same as x ticks
    a.minorticks_off()
    a.set_aspect('equal')
    a.spines['left'].set_position(('outward', 4))
    a.spines['bottom'].set_position(('outward', 4))
    t, p = ttest_rel(f1f0ons, f1f0offs) # paired t-test
    a.add_artist(AnchoredText('p$=$%.1g' % p, loc='upper left', frameon=False))
    ## scatter plot all collected maxori phis, regardless of burst classification:
    figsize = DEFAULTFIGURESIZE[0]*1.02, DEFAULTFIGURESIZE[1] # tweak to make space for labels
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto phi maxori grating %s CYCPSTHTHRESH=%g F1F0THRESH=%g '
             'TESTPEAKS=%s' % (st8, CYCPSTHTHRESH, F1F0THRESH, TESTPEAKS))
    # plot y=x line:
    #linmin, linmax, linstep = -180, 180, 180
    linmin, linmax, linstep = 0, 360, 180
    xyline = [linmin, linmax], [linmin, linmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot normal (non-example) points:
    a.scatter(phions[normlis]*180/pi, phioffs[normlis]*180/pi, 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[grtmseu2exmpli[mseustr]]
        c = exmpli2clr[grtmseu2exmpli[mseustr]]
        sz = exmpli2sz[grtmseu2exmpli[mseustr]]
        lw = exmpli2lw[grtmseu2exmpli[mseustr]]
        a.scatter(phions[exmpli]*180/pi, phioffs[exmpli]*180/pi, clip_on=False,
                  marker=marker, c=c, s=sz, lw=lw)
    a.set_xlabel('Suppression')# $\phi$ ($\degree$)')
    a.set_ylabel('Feedback')# $\phi$ ($\degree$)')
    a.set_xlim(linmin, linmax)
    a.set_ylim(linmin-10.5, linmax+52)
    a.set_xticks(np.arange(linmin, linmax+linstep, linstep))
    a.set_yticks(a.get_xticks()) # make y ticks the same as x ticks
    a.minorticks_off()
    a.set_aspect('equal')
    a.spines['left'].set_position(('outward', 4))
    a.spines['bottom'].set_position(('outward', 4))
    #z, p = rayleigh_test(dphis) # test if dphis is uniform, not sure that's useful here though
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
    ## plot maxori delta phi PDF:
    figsize = DEFAULTFIGURESIZE
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto delta phi grating pdf %s BURSTRATIOTHRESH=%g' % (st8, BURSTRATIOTHRESH))
    phibins = np.arange(-180, 180+15, 15)
    a.hist(dnbphis*180/pi, bins=phibins, color='black', histtype='step') # non burst mseus
    a.hist(dbphis*180/pi, bins=phibins, color='red', histtype='step') # burst mseus
    #a.axvline(x=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
    a.set_xlabel('$\Delta\phi=\phi_{\mathregular{fb}}-\phi_{\mathregular{sup}}$ ($\degree$)')
    a.set_ylabel('Unit count')
    a.set_xlim(-180, 180)
    a.set_xticks([-180, 0, 180])
    #p, k = kuiper(dnbphis, dbphis, axis=0) # test if dnbphis and dbphis are different
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
    ## plot maxori delta phi CDF:
    figsize = DEFAULTFIGURESIZE
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto delta phi grating cdf %s BURSTRATIOTHRESH=%g' % (st8, BURSTRATIOTHRESH))
    # add extra max bin to prevent vertical line:
    dnbphibins = list(np.unique(dnbphis*180/pi)) + [180]
    a.hist(dnbphis*180/pi, bins=dnbphibins, density=True, cumulative=True, histtype='step',
           color='black') # non burst mseus
    dbphibins = list(np.unique(dbphis*180/pi)) + [180]
    a.hist(dbphis*180/pi, bins=dbphibins, density=True, cumulative=True, histtype='step',
           color='red') # non burst mseus
    a.axvline(x=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
    a.set_xlabel('$\Delta\phi=\phi_{\mathregular{fb}}-\phi_{\mathregular{sup}}$ ($\degree$)')
    a.set_ylabel('Cumulative probability')
    a.set_xlim(-110, 110)
    #a.set_xticks([-180, 0, 180])
    a.set_yticks([0, 0.5, 1])
    nnbadv = (dnbphis > 0).sum() # non-bursting units whose phase advanced
    nnbtotal = len(dnbphis) # total number of non-bursting units
    pnb = scipy.stats.binom_test(nnbadv, nnbtotal) # test if nnbadv and nnbtotal are equal
    print('non-burst units that advanced: %d/%d, p_binom=%g' % (nnbadv, nnbtotal, pnb))
    nbadv = (dbphis > 0).sum() # bursting units whose phase advanced
    nbtotal = len(dbphis) # total number of bursting units
    pb = scipy.stats.binom_test(nbadv, nbtotal) # test if nbadv and nbtotal are equal
    print('burst units that advanced: %d/%d, p_binom=%g' % (nbadv, nbtotal, pb))


# scatter plot grating meanrates:
figsize = DEFAULTFIGURESIZE[0]*1.02, DEFAULTFIGURESIZE[1] # tweak to make space for log units
logmin, logmax = -1.1, 2
if EXPTYPE == 'ntsrmvis':
    logmin, logmax = -1.5, 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' % st8)
    rons, roffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], []
    keptmseui = 0 # manually init and increment instead of using enumerate()
    for mseustr in grtmseustrs:
        trialis = grtresp.loc[mseustr, st8]['trialis'] # non-blank & blank trialis
        if trialis.isna().any(): # missing for at least one opto condition
            continue
        ntrials = { opto:len(trialis[opto]) for opto in OPTOS } # should be equal
        ratesfull = pd.Series({ opto:np.full(ntrials[opto], np.nan) # pad
                                for opto in [False, True] })
        meanrate = grtresp.loc[mseustr, st8]['meanrate']
        if meanrate.isna().any(): # missing for at least one opto condition, can't plot
            fig3.loc[mseustr][blnksname] = ratesfull # save nans for all trials
            continue
        rons.append(meanrate[True])
        roffs.append(meanrate[False])
        rates = grtresp.loc[mseustr, st8]['rates'] # trial-wise non-blank rates
        nnblnktrials = { opto:len(rates[opto]) for opto in OPTOS } # num non-blank trials
        for opto in OPTOS:
            ratesfull[opto][:nnblnktrials[opto]] = rates[opto]
        _, pval = ttest_ind(rates[False], rates[True], equal_var=False)
        sgnfs.append(pval < SCATTERPTHRESH) # bool
        fig3.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
        fig3.loc[mseustr]['trialis'] = trialis # save, for joining DataFrames trial-wise
        fig3.loc[mseustr]['rates'] = ratesfull # save padded trial-wise values
        if mseustr in grtmseu2exmpli:
            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)
    nsgnabovediag = (roffs[sgnfis] > rons[sgnfis]).sum()
    nsgnbelowdiag = (roffs[sgnfis] < rons[sgnfis]).sum()
    print('npoints=%d, nsgnf=%d, ninsgnf=%d, nsgnabovediag=%d, nsgnbelowdiag=%d'
          % (len(sgnfs), len(sgnfis), len(insgnfis), nsgnabovediag, nsgnbelowdiag))
    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[grtmseu2exmpli[mseustr]]
        c = exmpli2clr[grtmseu2exmpli[mseustr]]
        sz = exmpli2sz[grtmseu2exmpli[mseustr]]
        lw = exmpli2lw[grtmseu2exmpli[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))


# scatter plot grating 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 st8 in ['none']:#ALLST8S:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto burst ratio grating %s' % st8)
    brons, broffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], []
    keptmseui = 0 # manually init and increment instead of using enumerate()
    for mseustr in grtmseustrs:
        trialis = grtresp.loc[mseustr, st8]['trialis'] # non-blank & blank trialis
        if trialis.isna().any(): # missing for at least one opto condition
            continue
        ntrials = { opto:len(trialis[opto]) for opto in OPTOS } # should be equal
        burstratiosfull = pd.Series({ opto:np.full(ntrials[opto], np.nan) # pad
                                      for opto in [False, True] })
        br = grtresp.loc[mseustr, st8]['meanburstratio']
        if br.isna().any(): # missing for at least one opto condition, can't plot
            fig3.loc[mseustr]['burstratios'] = burstratiosfull # save nans for all trials
            continue
        brons.append(br[True])
        broffs.append(br[False])
        burstratios = grtresp.loc[mseustr, st8]['burstratios'] # non-blank burst ratios
        nnblnktrials = { opto:len(burstratios[opto]) for opto in OPTOS } # num non-blank trials
        for opto in OPTOS:
            burstratiosfull[opto][:nnblnktrials[opto]] = burstratios[opto]
        _, pval = ttest_ind(burstratios[False], burstratios[True], equal_var=False)
        sgnfs.append(pval < SCATTERPTHRESH) # bool
        fig3.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
        fig3.loc[mseustr]['burstratios'] = burstratiosfull # save padded trial-wise values
        if mseustr in grtmseu2exmpli:
            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[grtmseu2exmpli[mseustr]]
        c = exmpli2clr[grtmseu2exmpli[mseustr]]
        sz = exmpli2sz[grtmseu2exmpli[mseustr]]
        lw = exmpli2lw[grtmseu2exmpli[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 grating opto OSI, R2, and orientation preference:
oris = np.arange(360) # get nice smooth model output in 1 deg steps
OSITHRESH = 0.02 # minimum OSI required in both conditions to include in oripref scatter
R2THRESH = 0.5 # minimum R2 tuning curve fit threshold to include in oripref scatter
for st8 in ['none']:#ALLST8S:
    color = st82clr[st8]
    osioffs, osions, oproffs, oprons, R2offs, R2ons = [], [], [], [], [], []
    exmplis, exmplmseustrs, normlis, mseustrs = [], [], [], []
    keptmseui = 0 # manually init and increment instead of using enumerate()
    for mseustr in grtmseustrs:
        tunparams = grtresp.loc[mseustr, st8]['tun']
        tunrsqs = grtresp.loc[mseustr, st8]['tunrsq']
        if tunparams.isna().any(): # missing one or both tunparams
            print('%s: missing one or both opto conditions, skipping' % mseustr)
            continue
        offparams, onparams = tunparams[False], tunparams[True]
        # don't subtract spont (params[5]) - negative rates lead to OSI exceeding 1 a lot more:
        offrates = sum_of_gaussians(oris, *offparams[:5])# - offparams[5] # model output
        onrates = sum_of_gaussians(oris, *onparams[:5])# - onparams[5] # model output
        osioff, osion = vector_OSI(oris, offrates), vector_OSI(oris, onrates)
        osioffs.append(osioff)
        osions.append(osion)
        fig3.loc[mseustr, 'osi'] = osioff, osion # save
        if mseustr in grtmseu2exmpli: # print out OSI values for example units
            print('Example OSI %s, opto_on: %.3f, opto_off: %.3f' % (mseustr, osion, osioff))
        oproff, opron = offparams[0] % 180, onparams[0] % 180 # direction pref mod 180
        oproffs.append(oproff)
        oprons.append(opron)
        R2offs.append(tunrsqs[False])
        R2ons.append(tunrsqs[True])
        if mseustr in grtmseu2exmpli:
            exmplis.append(keptmseui)
            exmplmseustrs.append(mseustr)
        else:
            normlis.append(keptmseui)
        mseustrs.append(mseustr)
        keptmseui += 1 # manually increment
    osioffs, osions = np.asarray(osioffs), np.asarray(osions)
    oproffs, oprons = np.asarray(oproffs), np.asarray(oprons)
    R2offs, R2ons = np.asarray(R2offs), np.asarray(R2ons)
    mseustrs = np.asarray(mseustrs)
    ## scatter plot OSI:
    figsize = DEFAULTFIGURESIZE[0]*1.02, DEFAULTFIGURESIZE[1] # tweak for log units
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto osi grating %s' % st8)
    logmin, logmax, logstep = -2.2, 0, 1
    logticks = np.array([-2, -1, 0])
    # replace off-scale low values with log0min, so the points remain visible:
    #pltosions, pltosioffs = osions.copy(), osioffs.copy()
    #pltosions[pltosions <= 10**logmin] = 10**log0min
    #pltosioffs[pltosioffs <= 10**logmin] = 10**log0min
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    #xyline = [linmin, linmax], [linmin, linmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot normal (non-example) points:
    a.scatter(osions[normlis], osioffs[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[grtmseu2exmpli[mseustr]]
        c = exmpli2clr[grtmseu2exmpli[mseustr]]
        sz = exmpli2sz[grtmseu2exmpli[mseustr]]
        lw = exmpli2lw[grtmseu2exmpli[mseustr]]
        a.scatter(osions[exmpli], osioffs[exmpli], clip_on=False,
                  marker=marker, c=c, s=sz, lw=lw)
    a.set_xlabel('Suppression')# OSI')
    a.set_ylabel('Feedback')# OSI')
    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 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(osions, osioffs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
    ## scatter plot R2:
    figsize = DEFAULTFIGURESIZE
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto R2 grating %s' % st8)
    linmin, linmax, linstep = 0, 1, 0.5
    xyline = [linmin, linmax], [linmin, linmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot normal (non-example) points, coloring by state:
    a.scatter(R2ons[normlis], R2offs[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[grtmseu2exmpli[mseustr]]
        c = exmpli2clr[grtmseu2exmpli[mseustr]]
        sz = exmpli2sz[grtmseu2exmpli[mseustr]]
        lw = exmpli2lw[grtmseu2exmpli[mseustr]]
        a.scatter(R2ons[exmpli], R2offs[exmpli], clip_on=False,
                  marker=marker, c=c, s=sz, lw=lw)
    a.set_xlabel(r'Suppression $\mathregular{R^2}$')
    a.set_ylabel(r'Feedback $\mathregular{R^2}$')
    a.set_xlim(linmin, linmax)
    a.set_ylim(linmin, linmax)
    a.set_xticks(np.arange(linmin, linmax+linstep, linstep))
    a.set_yticks(a.get_xticks())
    a.minorticks_off()
    a.set_aspect('equal')
    a.spines['left'].set_position(('outward', 4))
    a.spines['bottom'].set_position(('outward', 4))
    t, p = ttest_rel(R2ons, R2offs) # paired t-test
    a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
    ## scatter plot oripref:
    figsize = DEFAULTFIGURESIZE[0]*1.02, DEFAULTFIGURESIZE[1] # tweak tick labels
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto oripref grating %s OSITHRESH=%g R2THRESH=%g' % (st8, OSITHRESH, R2THRESH))
    linmin, linmax, linstep = 0, 180, 90
    xyline = [linmin, linmax], [linmin, linmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    if EXPTYPE == 'pvmvis':
        # wrap some points:
        wrapis = (oproffs < 15) & (oprons > 175)
        assert wrapis.sum() == 2
        oproffs[wrapis] = oproffs[wrapis] + 180 # wrap 2 points up over the y=x line
    # plot only units whose OSI satisfies OSITHRESH and R2 satisfies R2THRESH,
    # get indices of units that are well-tuned and well-fit in both conditions:
    tunis, = np.where((osioffs >= OSITHRESH) & (osions >= OSITHRESH) &
                      (R2offs >= R2THRESH) & (R2ons >= R2THRESH))
    normlis = np.intersect1d(tunis, normlis)
    exmplis = np.intersect1d(tunis, exmplis)
    for mseustr, oproff, opron in zip(mseustrs[tunis], oproffs[tunis], oprons[tunis]):
        fig3.loc[mseustr, 'oripref'] = oproff, opron # save
    # plot normal (non-example) points, coloring by state:
    a.scatter(oprons[normlis], oproffs[normlis], clip_on=False,
              marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
    # plot normal (non-example) points, coloring by min OSI:
    '''
    minosi = np.vstack([osions[normlis], osioffs[normlis]]).min(axis=0)
    scatt = a.scatter(oproffs[normlis], oprons[normlis],
                      marker='.', c=minosi, vmax=0.1, s=100)
    f.colorbar(scatt)
    '''
    # plot example points, one at a time:
    for exmpli, mseustr in zip(exmplis, exmplmseustrs):
        marker = exmpli2mrk[grtmseu2exmpli[mseustr]]
        c = exmpli2clr[grtmseu2exmpli[mseustr]]
        sz = exmpli2sz[grtmseu2exmpli[mseustr]]
        lw = exmpli2lw[grtmseu2exmpli[mseustr]]
        a.scatter(oprons[exmpli], oproffs[exmpli], clip_on=False,
                  marker=marker, c=c, s=sz, lw=lw)
    a.set_xlabel(r'Suppression')# $\theta$ ($\degree$)')
    a.set_ylabel(r'Feedback')# $\theta$ ($\degree$)')
    a.set_xlim(linmin, linmax)
    a.set_ylim(ymin=linmin)
    a.set_xticks(np.arange(linmin, linmax+linstep, linstep))
    a.set_yticks(a.get_xticks())
    a.minorticks_off()
    a.set_aspect('equal')
    a.spines['left'].set_position(('outward', 4))
    a.spines['bottom'].set_position(('outward', 4))
    #t, p = ttest_rel(oprons[tunis], oproffs[tunis]) # paired t-test
    #oriprefm, oriprefb, oriprefr, _, _ = linregress(oproffs[tunis], oprons[tunis])
    #oriprefrsq = oriprefr * oriprefr
    #txt = ('$\mathregular{R^2=}$%.2g\n'
    #       '$\mathregular{p=}$%.2g' % (oriprefrsq, p))
    #a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))
    ## plot OSI distributions, for choosing OSITHRESH:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto osi grating pdf %s' % st8)
    osibins = np.logspace(-2, 0, 15)
    offclr, onclr = desat(color, opto2alpha[False]), desat(color, opto2alpha[True])
    a.hist(osioffs, bins=osibins, color=offclr, histtype='step')
    a.hist(osions, bins=osibins, color=onclr, histtype='step')
    a.axvline(x=OSITHRESH, ls='--', marker='', color='lightgray', zorder=-np.inf)
    a.set_xlabel('OSI')
    a.set_ylabel('Unit count')
    a.set_xscale('log')
    ## plot R2 distributions, for choosing R2THRESH:
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto R2 grating pdf %s' % st8)
    R2bins = np.linspace(0, 1, 11)
    a.hist(R2offs, bins=R2bins, color=offclr, histtype='step')
    a.hist(R2ons, bins=R2bins, color=onclr, histtype='step')
    a.axvline(x=R2THRESH, ls='--', marker='', color='lightgray', zorder=-np.inf)
    a.set_xlabel('Fit $\mathregular{R^2}$')
    a.set_ylabel('Unit count')


# scatter plot grating opto mean f1 and f1/f0:
logmin, logmax = -1, 2
logticks = np.array([-1, 0, 1, 2])
for st8 in ['none']:#ALLST8S:
    f1offs, f1ons, f1f0offs, f1f0ons, exmplis, exmplmseustrs, normlis = [], [], [], [], [], [], []
    keptmseui = 0 # manually init and increment instead of using enumerate()
    for mseustr in grtmseustrs:
        rasters = grtresp.loc[mseustr, st8]['raster']
        dt = grtresp.loc[mseustr, st8, False]['dt']
        tfreq = grtresp.loc[mseustr, st8, False]['tfreq']
        if rasters.isna().any(): # missing one or both rasters
            print('%s: missing one or both opto conditions, skipping' % mseustr)
            continue
        offraster, onraster = rasters[False], rasters[True]
        ## maybe consider std as well, to gauge significance of each f1f0 measure...
        f0off, _ = raster2freqcomp(offraster, dt, 0, mean='scalar')
        f1off, _ = raster2freqcomp(offraster, dt, tfreq, mean='scalar')
        f0on, _ = raster2freqcomp(onraster, dt, 0, mean='scalar')
        f1on, _ = raster2freqcomp(onraster, dt, tfreq, mean='scalar')
        f1f0off = f1off / f0off
        f1f0on = f1on / f0on
        f1offs.append(f1off)
        f1ons.append(f1on)
        f1f0offs.append(f1f0off)
        f1f0ons.append(f1f0on)
        fig3.loc[mseustr, 'f1f0'] = f1f0off, f1f0on # save
        if mseustr in grtmseu2exmpli:
            exmplis.append(keptmseui)
            exmplmseustrs.append(mseustr)
        else:
            normlis.append(keptmseui)
        keptmseui += 1 # manually increment
    f1offs = np.asarray(f1offs)
    f1ons = np.asarray(f1ons)
    f1f0offs = np.asarray(f1f0offs)
    f1f0ons = np.asarray(f1f0ons)
    ## scatter plot mean f1:
    figsize = DEFAULTFIGURESIZE
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto mean f1 grating %s' % st8)
    # plot y=x line:
    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
    #linmin, linmax = 0, 75
    #xyline = [linmin, linmax], [linmin, linmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot normal (non-example) points:
    a.scatter(f1ons[normlis], f1offs[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[grtmseu2exmpli[mseustr]]
        c = exmpli2clr[grtmseu2exmpli[mseustr]]
        sz = exmpli2sz[grtmseu2exmpli[mseustr]]
        lw = exmpli2lw[grtmseu2exmpli[mseustr]]
        a.scatter(f1ons[exmpli], f1offs[exmpli], marker=marker, clip_on=False,
                  c=c, s=sz, lw=lw)
    a.set_xlabel('Suppression F1')
    a.set_ylabel('Feedback F1')
    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 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(f1ons, f1offs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
    ## scatter plot mean f1/f0
    figsize = DEFAULTFIGURESIZE
    f, a = plt.subplots(figsize=figsize)
    wintitle('opto mean f1f0 grating %s' % st8)
    # plot y=x line:
    linmin, linmax = 0, 2
    xyline = [linmin, linmax], [linmin, linmax]
    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
    # plot normal (non-example) points:
    a.scatter(f1f0ons[normlis], f1f0offs[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[grtmseu2exmpli[mseustr]]
        c = exmpli2clr[grtmseu2exmpli[mseustr]]
        sz = exmpli2sz[grtmseu2exmpli[mseustr]]
        lw = exmpli2lw[grtmseu2exmpli[mseustr]]
        a.scatter(f1f0ons[exmpli], f1f0offs[exmpli], clip_on=False,
                  marker=marker, c=c, s=sz, lw=lw)
    a.set_xlabel('Suppression')# F1/F0')
    a.set_ylabel('Feedback')# F1/F0')
    #a.set_xscale('log', basex=2)
    #a.set_yscale('log', basey=2)
    #a.set_xlim(0.125, 2)
    #a.set_ylim(0.125, 2)
    #a.set_xticks([0.125, 0.25, 0.5, 1, 2])
    #a.set_yticks(a.get_xticks()) # make scale y ticks the same as x ticks
    #axes_disable_scientific(a)
    a.set_xlim(linmin, linmax)
    a.set_ylim(linmin, linmax)
    a.set_xticks(np.arange(linmin, linmax+1))
    a.set_yticks(a.get_xticks()) # make scale y ticks the same as x ticks
    a.minorticks_off()
    a.set_aspect('equal')
    a.spines['left'].set_position(('outward', 4))
    a.spines['bottom'].set_position(('outward', 4))
    #t, p = ttest_rel(f1f0ons, f1f0offs) # paired t-test
    #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))

'''
# scatter plot grating opto local/global max f1/f0 and max f1:
figsize = DEFAULTFIGURESIZE
logmin, logmax = -1.5, 2
logticks = np.array([-1, 0, 1, 2])
for maxmethod in ['local', 'global']:
    for st8 in ['none']:#ALLST8S:
        maxf1f0s, maxf1s, exmplis, exmplmseustrs, normlis = [], [], [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in grtmseustrs:
            f0s, f1s = [[], []], [[], []]
            rasters = grtresp.loc[mseustr, st8]['raster']
            dt = grtresp.loc[mseustr, st8, False]['dt']
            tfreq = grtresp.loc[mseustr, st8, False]['tfreq']
            if rasters.isna().any(): # missing one or both rasters
                print('%s: missing one or both opto conditions, skipping' % mseustr)
                continue
            offraster, onraster = rasters[False], rasters[True]
            for optoi, opto in enumerate(OPTOS):
                stimis = grtresp.loc[mseustr, st8, opto]['stimis']
                ustimis = np.unique(stimis)
                for stimi in ustimis:
                    trialis = stimis == stimi # boolean array
                    stimiraster = rasters[opto][trialis] # raster including only stimi's trials
                    f0, _ = raster2freqcomp(stimiraster, dt, 0, mean='vector')
                    f1, _ = raster2freqcomp(stimiraster, dt, tfreq, mean='vector')
                    f0s[optoi].append(f0)
                    f1s[optoi].append(f1)
            f0s = np.array(f0s) # 2D array, row0:False, row1:True
            f0s[f0s == 0] = np.nan # prevent div by 0, replace 0s with nans
            f1s = np.array(f1s)
            f1f0s = f1s / f0s # 2D array, row0:False, row1:True
            if maxmethod == 'local':
                # get max for each opto condition separately:
                #maxf1f0 = np.nanmax(f1f0s, axis=1)
                maxrowis = np.array([0, 1])
                maxcolis = np.nanargmax(f1f0s, axis=1) # 2 values, one per row
                maxf1f0 = f1f0s[maxrowis, maxcolis]
                maxf1 = f1s[maxrowis, maxcolis]
            elif maxmethod == 'global':
                # find global max, use it to pick paired value of other opto condition:
                flatmaxi = np.nanargmax(f1f0s)
                maxrowi, maxcoli = np.unravel_index(flatmaxi, f1f0s.shape)
                maxf1f0 = f1f0s[:, maxcoli]
                maxf1 = f1s[:, maxcoli]
            maxf1f0s.append(maxf1f0)
            maxf1s.append(maxf1)
            if mseustr in grtmseu2exmpli:
                exmplis.append(keptmseui)
                exmplmseustrs.append(mseustr)
            else:
                normlis.append(keptmseui)
            keptmseui += 1 # manually increment
        maxf1f0s = np.asarray(maxf1f0s).T # 2D array: optoi, mseui
        maxf1s = np.asarray(maxf1s).T # 2D array: optoi, mseui
        ## scatter plot max f1/f0:
        f, a = plt.subplots(figsize=figsize)
        wintitle('opto %s max f1f0 grating %s' % (maxmethod, st8))
        # plot y=x line:
        #xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
        linmin, linmax = 0, 2
        xyline = [linmin, linmax], [linmin, linmax]
        a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
        # plot normal (non-example) points:
        a.scatter(maxf1f0s[1, normlis], maxf1f0s[0, normlis],
                  marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
        # plot example points, one at a time:
        for exmpli, mseustr in zip(exmplis, exmplmseustrs):
            marker = exmpli2mrk[grtmseu2exmpli[mseustr]]
            c = exmpli2clr[grtmseu2exmpli[mseustr]]
            sz = exmpli2sz[grtmseu2exmpli[mseustr]]
            a.scatter(maxf1f0s[1, exmpli], maxf1f0s[0, exmpli], marker=marker, c=c, s=sz)
        a.set_xlabel('Suppression max F1/F0')
        a.set_ylabel('Feedback max F1/F0')
        a.set_xlim(linmin, linmax+0.1)
        a.set_ylim(linmin, linmax+0.1)
        a.set_xticks(np.arange(linmin, linmax+1))
        #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 scale y ticks the same as x ticks
        a.minorticks_off()
        a.set_aspect('equal')
        t, p = ttest_rel(maxf1f0s[1], maxf1f0s[0]) # paired t-test
        a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
        ## scatter plot max f1:
        f, a = plt.subplots(figsize=figsize)
        wintitle('opto %s max f1 grating %s' % (maxmethod, st8))
        # plot y=x line:
        xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
        #linmin, linmax = 0, 2
        #xyline = [linmin, linmax], [linmin, linmax]
        a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
        # plot normal (non-example) points:
        a.scatter(maxf1s[1, normlis], maxf1s[0, normlis],
                  marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
        # plot example points, one at a time:
        for exmpli, mseustr in zip(exmplis, exmplmseustrs):
            marker = exmpli2mrk[grtmseu2exmpli[mseustr]]
            c = exmpli2clr[grtmseu2exmpli[mseustr]]
            sz = exmpli2sz[grtmseu2exmpli[mseustr]]
            a.scatter(maxf1s[1, exmpli], maxf1s[0, exmpli], marker=marker, c=c, s=sz)
        a.set_xlabel('Suppression max F1')
        a.set_ylabel('Feedback max F1')
        #a.set_xlim(linmin, linmax+0.1)
        #a.set_ylim(linmin, linmax+0.1)
        #a.set_xticks(np.arange(linmin, linmax+1))
        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 scale y ticks the same as x ticks
        a.minorticks_off()
        axes_disable_scientific(a)
        a.set_aspect('equal')
        t, p = ttest_rel(maxf1s[1], maxf1s[0]) # paired t-test
        a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
'''

# fig4S1a: scatter plot grating meanburstratio FMI vs meanrate FMI:
figsize = DEFAULTFIGURESIZE[0]*0.95, DEFAULTFIGURESIZE[1] # tweak to match others in 4S1
linmin, linmax, linstep = -1, 1, 1
ticks = np.arange(linmin, linmax+linstep, linstep)
for st8 in ['none']:#ALLST8S:
    for measure in ['meanburstratio']:
        if measure.startswith('meanrate'):
            continue # don't bother plotting meanrate* against meanrate
        axislabel = measure2axislabel[measure]
        rfmis, msrfmis = [], []
        exmplis, exmplmseustrs, normlis = [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in grtmseustrs:
            rfmi = grtFMI.loc[mseustr, st8]['meanrate']
            msrfmi = grtFMI.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 mseustr in grtmseu2exmpli:
                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' % (measure, 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 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(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[grtmseu2exmpli[mseustr]]
            c = exmpli2clr[grtmseu2exmpli[mseustr]]
            sz = exmpli2sz[grtmseu2exmpli[mseustr]]
            lw = exmpli2lw[grtmseu2exmpli[mseustr]]
            a.scatter(rfmis[exmpli], msrfmis[exmpli], marker=marker, c=c, s=sz, lw=lw)
        #mm, b, rr, p, stderr = linregress(rfmis, msrfmis)
        #rsq = rr * rr
        #x = np.array([rfmis.min(), rfmis.max()])
        #y = mm * x + b
        #txt = ('$\mathregular{R^{2}=%.2g}$' % rsq)
        #a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))
        #a.plot(x, y, '-', color='red') # plot linregress fit
        a.set_xlabel('Firing rate 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))