"""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 # 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) # 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_5_S1c_coefs.csv') elif measure == 'spars': # use Steffen's LMM linregress fit, for fig5S1 fname = os.path.join('stats', 'figure_5_S1d_coefs.csv') elif measure == 'rel': # use Steffen's LMM linregress fit, for fig5S1 fname = os.path.join('stats', 'figure_5_S1e_coefs.csv') elif measure == 'snr': # use Steffen's LMM linregress fit, for fig5S1 fname = os.path.join('stats', 'figure_5_S1f_coefs.csv') elif measure == 'meanpkw': # use Steffen's LMM linregress fit, for fig5S1 fname = os.path.join('stats', 'figure_5_S1g_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))