natfeedback_code_eLife/fig5.py

"""Figure 5 and 5S1 plots, use run -i fig5.py"""

mi = pd.MultiIndex.from_product([mvimseustrs, ST8S, OPTOS], names=['mseu', 'st8', 'opto'])
columns = ['trialis', 'rates', 'burstratios', 'meanrate', 'meanburstratio',
           'spars', 'rel', 'meanpkw', 'snr']
fig5 = pd.DataFrame(index=mi, columns=columns) # excludes kind and opto
## NOTE: make sure that rows are not overwritten due to temporarily enabling iteration
## over kind in any of the following loops!


# tweak exactly how much horizontal space to use for raster plots and PSTHs:
PLOTOFFSETS = -0.45, 0.45

# plot movie rasters for each unit by kind, state and opto condition:
# copied and slightly modified from fig1.py:
showbursts = True
trialiis = None #np.arange(120) # plot only these trials, None plots all trials
for mseustr in mvimseu2exmpli: # mvimseustrs
    if mseustr != 'PVCre_2018_0003_s03_e03_u51': # plot only this one
        continue
    for kind in MVIKINDS:
        for st8 in ['run', 'sit']:
            for opto in OPTOS:
                dt = mviresp.loc[mseustr, kind, st8, opto]['dt'] # stimulus duration (s)
                #wdt = -PLOTOFFSETS[0] + dt + PLOTOFFSETS[1] # wide raster duration, sec
                tmin, tmax = 0 + PLOTOFFSETS[0], dt + PLOTOFFSETS[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=PLOTOFFSETS,
                                  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 + PLOTOFFSETS[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:
# copied and slightly modified from fig1.py:
for mseustr in mvimseu2exmpli: #mvimseustrs
    if mseustr != 'PVCre_2018_0003_s03_e03_u51': # plot only this one
        continue
    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 = -PLOTOFFSETS[0] + dt + PLOTOFFSETS[1] # wide raster duration, sec
            tmin, tmax = 0 + PLOTOFFSETS[0], dt + PLOTOFFSETS[1]
            #figsize = 0.872 + wdt*1.5*RASTERSCALEX, PSTHHEIGHT
            figsize = 0.95 + 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(t[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)
            a.set_ylim(-2.75, 57)
            # 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


# scatter plot movie run vs sit meanrate, during feedback condition:
figsize = DEFAULTFIGURESIZE[0]*1.05, DEFAULTFIGURESIZE[1] # tweak to make space for log units
logmin, logmax = -2, 2
logticks = np.array([-2, 0, 2])
for kind in ['nat']:#MVIKINDS:
    for opto in OPTOS:
        f, a = plt.subplots(figsize=figsize)
        wintitle('runsit meanrate movie %s %s' % (kind, opto))
        rruns, rsits, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in mvimseustrs:
            meanrate = mviresp.loc[mseustr, kind, ST8S, opto]['meanrate']
            snr = mviresp.loc[mseustr, kind, ST8S, opto]['snr']
            if meanrate.isna().any(): # missing one or both meanrates
                continue
            if (snr < SNRTHRESH).all(): # neither condition has decent SNR
                continue
            rrun, rsit = meanrate.values
            rruns.append(rrun)
            rsits.append(rsit)
            trialis = mviresp.loc[mseustr, kind, ST8S, opto]['trialis']
            rates = mviresp.loc[mseustr, kind, ST8S, opto]['rates']
            runrates = rates.xs('run', level='st8').iloc[0]
            sitrates = rates.xs('sit', level='st8').iloc[0]
            _, pval = ttest_ind(runrates, sitrates, equal_var=False)
            sgnfs.append(pval < SCATTERPTHRESH) # bool
            fig5.loc[(mseustr, ST8S, opto), 'meanrate'] = rrun, rsit # save
            fig5.loc[(mseustr, ST8S, opto), 'trialis'] = trialis.values # for trial-wise joins
            fig5.loc[(mseustr, ST8S, opto), 'rates'] = rates.values # save trial-wise values
            if mvimseu2exmpli.get(mseustr) == fig5exmpli:
                exmplis.append(keptmseui)
                exmplmseustrs.append(mseustr)
            else:
                normlis.append(keptmseui)
            keptmseui += 1 # manually increment
        rruns = np.asarray(rruns)
        rsits = np.asarray(rsits)
        sgnfs = np.asarray(sgnfs)
        # replace off-scale low values with log0min, so the points remain visible:
        #pltrruns, pltrsits = rruns.copy(), rsits.copy()
        #pltrruns[pltrruns <= 10**logmin] = 10**log0min
        #pltrsits[pltrsits <= 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(black, SGNF2ALPHA[False]) # do manual alpha mixing
        a.scatter(rsits[normlinsgnfis], rruns[normlinsgnfis], clip_on=False,
                  marker='.', c='None', edgecolor=c, s=DEFSZ)
        # plot normal (non-example) significant points:
        c = desat(black, SGNF2ALPHA[True]) # do manual alpha mixing
        a.scatter(rsits[normlsgnfis], rruns[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(rsits[exmpli], rruns[exmpli], clip_on=False,
                      marker=marker, c=c, s=sz, lw=lw)
        a.set_xlabel('Sit FR (spk/s)')
        a.set_ylabel('Run 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(rsits, rruns) # paired t-test
        #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))


# scatter plot movie run vs sit burst ratio, during feedback condition:
figsize = DEFAULTFIGURESIZE[0]*1.06, DEFAULTFIGURESIZE[1] # tweak to make space for log units
logmin, logmax = -3.55, 0
logticks = np.array([-3, -2, -1, 0])
for kind in ['nat']:#MVIKINDS:
    for opto in OPTOS:
        f, a = plt.subplots(figsize=figsize)
        wintitle('runsit burst ratio movie %s %s' % (kind, opto))
        brruns, brsits, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in mvimseustrs:
            br = mviresp.loc[mseustr, kind, ST8S, opto]['meanburstratio']
            snr = mviresp.loc[mseustr, kind, ST8S, opto]['snr']
            if br.isna().any(): # missing one or both values
                continue
            if (snr < SNRTHRESH).all(): # neither condition has decent SNR
                continue
            brrun, brsit = br.values
            brruns.append(brrun)
            brsits.append(brsit)
            burstratios = mviresp.loc[mseustr, kind, ST8S, opto]['burstratios']
            runburstratios = burstratios.xs('run', level='st8').iloc[0]
            sitburstratios = burstratios.xs('sit', level='st8').iloc[0]
            _, pval = ttest_ind(runburstratios, sitburstratios, equal_var=False)
            sgnfs.append(pval < SCATTERPTHRESH) # bool
            fig5.loc[(mseustr, ST8S, opto), 'meanburstratio'] = brrun, brsit # save
            fig5.loc[(mseustr, ST8S, opto), 'burstratios'] = burstratios.values # save trial-wise
            if mvimseu2exmpli.get(mseustr) == fig5exmpli:
                exmplis.append(keptmseui)
                exmplmseustrs.append(mseustr)
            else:
                normlis.append(keptmseui)
            keptmseui += 1 # manually increment
        brruns = np.asarray(brruns)
        brsits = np.asarray(brsits)
        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)
        # plot normal (non-example) insignificant points:
        c = desat(black, SGNF2ALPHA[False]) # do manual alpha mixing
        a.scatter(brsits[normlinsgnfis], brruns[normlinsgnfis], clip_on=True,
                  marker='.', c='None', edgecolor=c, s=DEFSZ) # clip_on=False fails
        # plot normal (non-example) significant points:
        c = desat(black, SGNF2ALPHA[True]) # do manual alpha mixing
        a.scatter(brsits[normlsgnfis], brruns[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(brsits[exmpli], brruns[exmpli], clip_on=False,
                      marker=marker, c=c, s=sz, lw=lw)
        a.set_xlabel('Sit BR')
        a.set_ylabel('Run 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(brsits, brruns) # paired t-test
        #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))


# scatter plot movie run vs sit sparseness, during feedback condition:
figsize = DEFAULTFIGURESIZE
linmin, linmax, linstep = 0, 1, 0.5
ticks = np.arange(intround(linmin), intround(linmax)+linstep, linstep)
for kind in ['nat']:#MVIKINDS:
    for opto in OPTOS:
        f, a = plt.subplots(figsize=figsize)
        wintitle('runsit sparseness movie %s %s' % (kind, opto))
        sparsruns, sparssits, exmplis, exmplmseustrs, normlis = [], [], [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in mvimseustrs:
            spars = mviresp.loc[mseustr, kind, ST8S, opto]['spars']
            snr = mviresp.loc[mseustr, kind, ST8S, opto]['snr']
            if spars.isna().any(): # missing one or both values
                continue
            if (snr < SNRTHRESH).all(): # neither condition has decent SNR
                continue
            sparsrun, sparssit = spars.values
            sparsruns.append(sparsrun)
            sparssits.append(sparssit)
            fig5.loc[(mseustr, ST8S, opto), 'spars'] = sparsrun, sparssit # save
            if mvimseu2exmpli.get(mseustr) == fig5exmpli:
                exmplis.append(keptmseui)
                exmplmseustrs.append(mseustr)
            else:
                normlis.append(keptmseui)
            keptmseui += 1 # manually increment
        sparsruns = np.asarray(sparsruns)
        sparssits = np.asarray(sparssits)
        # 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(sparssits[normlis], sparsruns[normlis], clip_on=False,
                  marker='.', c='None', edgecolor='black', 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(sparssits[exmpli], sparsruns[exmpli], clip_on=False,
                      marker=marker, c=c, s=sz, lw=lw)
        a.set_xlabel('Sit sparseness')
        a.set_ylabel('Run sparseness')
        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(sparssits, sparsruns) # paired t-test
        #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))


# scatter plot movie run vs sit reliability, during feedback condition:
figsize = DEFAULTFIGURESIZE[0]*0.98, DEFAULTFIGURESIZE[1] # tweak down to match others
linmin, linmax = 0, 0.63
ticks = 0, 0.5
for kind in ['nat']:#MVIKINDS:
    for opto in OPTOS:
        f, a = plt.subplots(figsize=figsize)
        wintitle('runsit reliability movie %s %s' % (kind, opto))
        relruns, relsits, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in mvimseustrs:
            rel = mviresp.loc[mseustr, kind, ST8S, opto]['rel']
            snr = mviresp.loc[mseustr, kind, ST8S, opto]['snr']
            if rel.isna().any(): # missing one or both values
                continue
            if (snr < SNRTHRESH).all(): # neither condition has decent SNR
                continue
            relrun, relsit = rel.values
            relruns.append(relrun)
            relsits.append(relsit)
            rhos = mviresp.loc[mseustr, kind, ST8S, opto]['rhos']
            runrhos = rhos.xs('run', level='st8').iloc[0]
            sitrhos = rhos.xs('sit', level='st8').iloc[0]
            _, pval = ttest_ind(runrhos, sitrhos, equal_var=False)
            sgnfs.append(pval < SCATTERPTHRESH) # bool
            fig5.loc[(mseustr, ST8S, opto), 'rel'] = relrun, relsit # save
            if mvimseu2exmpli.get(mseustr) == fig5exmpli:
                exmplis.append(keptmseui)
                exmplmseustrs.append(mseustr)
            else:
                normlis.append(keptmseui)
            keptmseui += 1 # manually increment
        relruns = np.asarray(relruns)
        relsits = np.asarray(relsits)
        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)
        # plot normal (non-example) insignificant points:
        c = desat(black, SGNF2ALPHA[False]) # do manual alpha mixing
        a.scatter(relsits[normlinsgnfis], relruns[normlinsgnfis], clip_on=False,
                  marker='.', c='None', edgecolor=c, s=DEFSZ)
        # plot normal (non-example) significant points:
        c = desat(black, SGNF2ALPHA[True]) # do manual alpha mixing
        a.scatter(relsits[normlsgnfis], relruns[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(relsits[exmpli], relruns[exmpli], clip_on=False,
                      marker=marker, c=c, s=sz, lw=lw)
        a.set_xlabel('Sit reliability')
        a.set_ylabel('Run reliability')
        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(relsits, relruns) # paired t-test
        #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))


# scatter plot movie sit vs. run SNR:
figsize = DEFAULTFIGURESIZE
linmin, linmax = 0, 0.65
ticks = 0, 0.5
for kind in ['nat']:#MVIKINDS:
    for opto in OPTOS:
        f, a = plt.subplots(figsize=figsize)
        wintitle('runsit SNR movie %s %s' % (kind, opto))
        snrruns, snrsits, exmplis, exmplmseustrs, normlis = [], [], [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in mvimseustrs:
            snr = mviresp.loc[mseustr, kind, ST8S, opto]['snr']
            if snr.isna().any(): # missing one or both values
                continue
            if (snr < SNRTHRESH).all(): # neither condition has decent SNR
                continue
            snrrun, snrsit = snr.values
            snrruns.append(snrrun)
            snrsits.append(snrsit)
            fig5.loc[(mseustr, ST8S, opto), 'snr'] = snrrun, snrsit # save
            if mvimseu2exmpli.get(mseustr) == fig5exmpli:
                exmplis.append(keptmseui)
                exmplmseustrs.append(mseustr)
            else:
                normlis.append(keptmseui)
            keptmseui += 1 # manually increment
        snrruns = np.asarray(snrruns)
        snrsits = np.asarray(snrsits)
        # 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(snrsits[normlis], snrruns[normlis], clip_on=False,
                  marker='.', c='None', edgecolor='black', 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(snrsits[exmpli], snrruns[exmpli], clip_on=False,
                      marker=marker, c=c, s=sz, lw=lw)
        a.set_xlabel('Sit SNR')
        a.set_ylabel('Run 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(snrsits, snrruns) # paired t-test
        #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))


# scatter plot movie run vs sit mean peak widths, during feedback condition:
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 opto in OPTOS:
        f, a = plt.subplots(figsize=figsize)
        wintitle('runsit mean peak width movie %s %s' % (kind, opto))
        meanpkwruns, meanpkwsits, exmplis, exmplmseustrs, normlis = [], [], [], [], []
        keptmseui = 0 # manually init and increment instead of using enumerate()
        for mseustr in mvimseustrs:
            meanpkw = mviresp.loc[mseustr, kind, ST8S, opto]['meanpkw']
            snr = mviresp.loc[mseustr, kind, ST8S, opto]['snr']
            if meanpkw.isna().any(): # missing one or both values
                continue
            if (snr < SNRTHRESH).all(): # neither condition has decent SNR
                continue
            meanpkwrun, meanpkwsit = meanpkw.values
            meanpkwruns.append(meanpkwrun)
            meanpkwsits.append(meanpkwsit)
            fig5.loc[(mseustr, ST8S, opto), 'meanpkw'] = meanpkwrun, meanpkwsit # save
            if mvimseu2exmpli.get(mseustr) == fig5exmpli:
                exmplis.append(keptmseui)
                exmplmseustrs.append(mseustr)
            else:
                normlis.append(keptmseui)
            keptmseui += 1 # manually increment
        meanpkwruns = np.asarray(meanpkwruns)
        meanpkwsits = np.asarray(meanpkwsits)
        # 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(meanpkwsits[normlis], meanpkwruns[normlis], clip_on=False,
                  marker='.', c='None', edgecolor='black', 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(meanpkwsits[exmpli], meanpkwruns[exmpli], clip_on=False,
                      marker=marker, c=c, s=sz, lw=lw)
        a.set_xscale('log')
        a.set_yscale('log')
        a.set_xlabel('Sit peak width (s)')
        a.set_ylabel('Run peak width (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))
        #t, p = ttest_rel(meanpkwsits, meanpkwruns) # paired t-test
        #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))


# scatter plot RMI of all movie measures vs meanrate RMI:
figsize = DEFAULTFIGURESIZE
linmin, linmax, linstep = -1, 1, 1
ticks = np.arange(linmin, linmax+linstep, linstep)
for kind in ['nat']:#MVIKINDS:
    for opto in [False]:#OPTOS:
        for measure in modmeasuresnoblank:
            if measure.startswith('meanrate'):
                continue # don't bother plotting meanrate* against meanrate
            axislabel = measure2axislabel[measure]
            rrmis, msrrmis = [], []
            exmplis, exmplmseustrs, normlis = [], [], []
            keptmseui = 0 # manually init and increment instead of using enumerate()
            for mseustr in mvimseustrs:
                rrmi = mviRMI.loc[mseustr, opto]['meanrate']
                msrrmi = mviRMI.loc[mseustr, opto][measure]
                if pd.isna(rrmi) or pd.isna(msrrmi): # missing at least one RMI value
                    continue
                rrmis.append(rrmi)
                msrrmis.append(msrrmi)
                if mvimseu2exmpli.get(mseustr) == fig1exmpli:
                    exmplis.append(keptmseui)
                    exmplmseustrs.append(mseustr)
                else:
                    normlis.append(keptmseui)
                keptmseui += 1 # manually increment
            rrmis, msrrmis = np.asarray(rrmis), np.asarray(msrrmis)
            assert len(rrmis) == len(msrrmis)
            print('%s vs meanrate: nobservations=%d' % (measure, len(rrmis)))
            ## scatter plot RMIs:
            f, a = plt.subplots(figsize=figsize)
            wintitle('opto RMI %s vs RMI rate %s %s' % (measure, kind, opto))
            # 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(rrmis[normlis], msrrmis[normlis], clip_on=False,
                      marker='.', c='None', edgecolor=opto2clr[opto], 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(rrmis[exmpli], msrrmis[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 fig5S1
                fname = os.path.join('stats', 'figure_5S1c_coefs.csv')
            elif measure == 'spars': # use Steffen's LMM linregress fit, for fig5S1
                fname = os.path.join('stats', 'figure_5S1d_coefs.csv')
            elif measure == 'rel': # use Steffen's LMM linregress fit, for fig5S1
                fname = os.path.join('stats', 'figure_5S1e_coefs.csv')
            elif measure == 'snr': # use Steffen's LMM linregress fit, for fig5S1
                fname = os.path.join('stats', 'figure_5S1f_coefs.csv')
            elif measure == 'meanpkw': # use Steffen's LMM linregress fit, for fig5S1
                fname = os.path.join('stats', 'figure_5S1g_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([rrmis.min(), rrmis.max()])
                    y = mm * x + b
                    a.plot(x, y, '-', color='red') # plot linregress fit
            a.set_xlabel('Firing rate RMI')
            if axislabel.islower():
                axislabel = axislabel.capitalize()
            ylabel = '%s RMI' % axislabel
            a.set_ylabel(ylabel)
            a.set_xlim(linmin, linmax)
            a.set_ylim(linmin, linmax)
            a.set_xticks(ticks)
            a.set_yticks(ticks)
            a.minorticks_off()
            a.set_aspect('equal')
            a.spines['left'].set_position(('outward', 4))
            a.spines['bottom'].set_position(('outward', 4))


# scatter plot movie vs grating RMI for applicable measures:
figsize = DEFAULTFIGURESIZE
linmin, linmax, linstep = -1, 1, 1
for measure in ['meanrate', 'meanburstratio']: # only 2 measures applicable to both mvi & grt
    axislabel = measure2axislabel[measure]
    for opto in OPTOS:
        f, a = plt.subplots(figsize=figsize)
        wintitle('maxRMI %s movie grating %s' % (measure, opto))
        mvimods, grtmods, exmplis, exmplmsustrs, normlis = [], [], [], [], []
        keptmsui = 0 # manually init and increment instead of using enumerate()
        for msustr in mvigrtmsustrs:
            mod = maxRMI[measure][msustr, opto]
            if mod.isna().any(): # missing one or both mod values
                continue
            mvimods.append(mod['mvi'])
            grtmods.append(mod['grt'])
            if msu2exmpli.get(msustr) == fig5exmpli:
                exmplis.append(keptmsui)
                exmplmsustrs.append(msustr)
            else:
                normlis.append(keptmsui)
            keptmsui += 1 # manually increment
        mvimods = np.asarray(mvimods)
        grtmods = np.asarray(grtmods)
        # plot y=x line:
        xyline = [linmin, linmax], [linmin, linmax]
        a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
        # plot normal (non-example) points:
        c = desat('black', opto2alpha[opto]) # do manual alpha mixing
        a.scatter(grtmods[normlis], mvimods[normlis], clip_on=False,
                  marker='.', c='None', edgecolor=c, s=DEFSZ)
        # plot example points, one at a time:
        for exmpli, msustr in zip(exmplis, exmplmsustrs):
            marker = exmpli2mrk[msu2exmpli[msustr]]
            c = exmpli2clr[msu2exmpli[msustr]]
            sz = exmpli2sz[msu2exmpli[msustr]]
            lw = exmpli2lw[msu2exmpli[msustr]]
            a.scatter(grtmods[exmpli], mvimods[exmpli], clip_on=False,
                      marker=marker, c=c, s=sz, lw=lw)
        a.set_xlabel('Grating %s RMI' % axislabel)
        a.set_ylabel('Movie %s RMI' % axislabel)
        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()) # make log 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(grtmods, mvimods) # paired t-test
        #a.add_artist(AnchoredText('$\mathregular{p=%.2g}$' % p,
        #             loc='upper left', frameon=False))