1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042 |
- """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))
- ## 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))
|