fig1.py 48 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966
  1. """Figure 1 and 1S2 plots, and first and last 2 sec scatter plots for 4S1, use run -i fig1.py"""
  2. mi = pd.MultiIndex.from_product([mvimseustrs, OPTOS], names=['mseu', 'opto'])
  3. # include single trial and trial-averaged columns:
  4. columns = ['trialis', 'rates', 'rate02s', 'rate35s', 'burstratios',
  5. 'blankrates', 'blankburstratios'] + modmeasuresnoblankcond
  6. fig1 = pd.DataFrame(index=mi, columns=columns) # excludes kind and st8
  7. ## NOTE: make sure that rows are not overwritten due to temporarily enabling iteration
  8. ## over kind and/or st8 in any of the following loops!
  9. # plot movie rasters for each unit by kind, state and opto condition:
  10. showbursts = True
  11. SCATTER = True # True: low quality, but fast; False: high quality, but slow
  12. trialiis = None #np.arange(120) # plot only these trials, None plots all trials
  13. for mseustr in []:#mvimseu2exmpli: #mvimseustrs:
  14. #['Ntsr1Cre_2019_0008_s03_e04_u18']
  15. #['Ntsr1Cre_2020_0001_s04_e10_u12']
  16. for kind in ['nat']: #MVIKINDS
  17. for st8 in ['none']:
  18. for opto in OPTOS:
  19. dt = mviresp.loc[mseustr, kind, st8, opto]['dt'] # stimulus duration (s)
  20. wdt = -OFFSETS[0] + dt + OFFSETS[1] # wide raster duration, sec
  21. tmin, tmax = 0 + OFFSETS[0], dt + OFFSETS[1]
  22. optotrange = mviresp.loc[mseustr, kind, st8, opto]['optotrange']
  23. raster = mviresp.loc[mseustr, kind, st8, opto]['wraster']
  24. if np.any(pd.isna(raster)): # no raster for this condition
  25. continue
  26. title = '%s movie %s %s %s raster' % (mseustr, kind, st8, opto)
  27. if trialiis is not None:
  28. raster = raster[trialiis]
  29. title += ' %dtrials' % len(trialiis)
  30. burstis = None
  31. if showbursts:
  32. wburstis = mviresp.loc[mseustr, kind, st8, opto]['wburstis']
  33. a = simpletraster(raster=raster, alpha=0.5, s=1, dt=dt, offsets=OFFSETS,
  34. scatter=SCATTER,
  35. title=title, inchespersec=1.5*RASTERSCALEX,
  36. inchespertrial=1/25*RASTERSCALEX,
  37. xaxis=True, burstis=wburstis)
  38. if opto:
  39. # plot horizontal line signifying opto ON period, just above axes:
  40. ton, toff = optotrange
  41. a.hlines(-5, ton, toff, colors=optoblue, lw=4, clip_on=False,
  42. in_layout=False)
  43. # plot horizontal line signifying stimulus period, just above axes:
  44. a.hlines(-2, 0, dt, colors='k', lw=4, clip_on=False, in_layout=False)
  45. # plot vertical coloured line just right of axes:
  46. #vlinexpos = dt + OFFSETS[1] + wdt*0.01
  47. #y0, y1 = a.get_ylim()
  48. #clr = st82clr[st8]
  49. #alpha = opto2alpha[opto]
  50. #a.vlines(vlinexpos, y0, y1, colors=clr, lw=4, alpha=alpha, clip_on=False,
  51. # in_layout=False)
  52. a.set_xlim(tmin, tmax)
  53. # plot movie PSTHs for each unit by kind, overplot state and opto combinations:
  54. for mseustr in []: #mvimseu2exmpli:#['Ntsr1Cre_2020_0001_s04_e10_u12']: #mvimseustrs
  55. for kind in MVIKINDS:
  56. wt = mviresp.loc[mseustr, kind]['wt'] # DataFrame
  57. wpsth = mviresp.loc[mseustr, kind]['wpsth'] # DataFrame
  58. wbpsth = mviresp.loc[mseustr, kind]['wbpsth'] # DataFrame
  59. #t = mviresp.loc[mseustr, kind]['t'] # DataFrame
  60. psth = mviresp.loc[mseustr, kind]['psth'] # DataFrame
  61. if np.isnan(np.hstack(wpsth)).any(): # no PSTH for at least one condition for this mseu
  62. continue
  63. optorho, optorho_p = scipy.stats.pearsonr(psth['none', False], psth['none', True])
  64. runsitrho, runsitrho_p = scipy.stats.pearsonr(psth['run', False], psth['sit', False])
  65. skewoff = scipy.stats.skew(psth['none', False])
  66. skewon = scipy.stats.skew(psth['none', True])
  67. skewrun = scipy.stats.skew(psth['run', False])
  68. skewsit = scipy.stats.skew(psth['sit', False])
  69. _, optoKS_p = ks_2samp(psth['none', False], psth['none', True])
  70. _, runsitKS_p = ks_2samp(psth['run', False], psth['sit', False])
  71. print('%s example PSTH:\n'
  72. 'optorho=%g optorho_p=%g skewoff=%g skewon=%g optoKS_p=%g\n'
  73. 'runsitrho=%g runsitrho_p=%g skewrun=%g skewsit=%g runsitKS_p=%g'
  74. % (mseustr,
  75. optorho, optorho_p, skewoff, skewon, optoKS_p,
  76. runsitrho, runsitrho_p, skewrun, skewsit, runsitKS_p))
  77. for st8combo, optocombo in zip(ST8COMBOS, OPTOCOMBOS):
  78. st8combostr = ' '.join(st8combo)
  79. optocombostr = ' '.join([str(opto) for opto in optocombo])
  80. dt = mviresp.loc[mseustr, kind, 'none']['dt'].mean() # mean stimulus duration (s)
  81. wdt = -OFFSETS[0] + dt + OFFSETS[1] # wide raster duration, sec
  82. tmin, tmax = 0 + OFFSETS[0], dt + OFFSETS[1]
  83. figsize = 0.872 + wdt*1.5*RASTERSCALEX, PSTHHEIGHT
  84. f, a = plt.subplots(figsize=figsize)
  85. title = '%s %s %s %s PSTH' % (mseustr, kind, st8combostr, optocombostr)
  86. wintitle(title)
  87. # subsample as in PSTH scatter plots, and plot all 4 st8/opto combinations:
  88. for st8 in st8combo:
  89. color = st82clr[st8]
  90. for opto in optocombo:
  91. if st8combo == ['run', 'sit'] and optocombo == [True]:
  92. optoalpha = 1
  93. else:
  94. optoalpha = opto2alpha[opto]
  95. if optocombo == [False]: # desaturate bursts by st8, not opto
  96. burstalpha = st82alpha[st8]
  97. else:
  98. burstalpha = opto2alpha[opto]
  99. c = desat(color, optoalpha) # do manual alpha mixing
  100. bc = desat(burstclr, burstalpha) # do manual alpha mixing
  101. fb = opto2fb[opto].title()
  102. if st8 == 'none':
  103. label = fb
  104. else:
  105. label = ', '.join([fb, st8.title()])
  106. a.plot(wt[st8, opto][::ssx], wpsth[st8, opto][::ssx],
  107. '-', color=c, label=label)
  108. a.plot(wt[st8, opto][::ssx], wbpsth[st8, opto][::ssx],
  109. '-', color=bc, label=label+', Burst')
  110. l = a.legend(frameon=False)
  111. l.set_draggable(True)
  112. a.set_xlabel('Time (s)')
  113. a.set_ylabel('Firing rate (spk/s)')
  114. a.set_xlim(tmin, tmax)
  115. # plot horizontal line signifying stimulus period, just below data:
  116. ymin, ymax = a.get_ylim()
  117. a.hlines(-ymax*0.025, 0, dt, colors='k', lw=4, clip_on=False, in_layout=False)
  118. a.set_ylim(ymin, ymax) # restore y limits from before plotting the hline
  119. # plot PDFs, CDFs and QQ of PSTHs for each unit by kind, overplot state and opto combinations:
  120. percentiles = np.arange(0, 100, 2)
  121. figsize = DEFAULTFIGURESIZE
  122. for mseustr in mvimseu2exmpli: #mvimseustrs
  123. if EXPTYPE != 'pvmvis':
  124. break # skip these plots
  125. for kind in MVIKINDS:
  126. psths = mviresp.loc[mseustr, kind]['psth']
  127. if np.isnan(np.hstack(psths)).any(): # no PSTH for at least one condition for this mseu
  128. continue
  129. for st8combo, optocombo in zip(ST8COMBOS, OPTOCOMBOS):
  130. st8combostr = ' '.join(st8combo)
  131. optocombostr = ' '.join([str(opto) for opto in optocombo])
  132. pdff, pdfa = plt.subplots(figsize=figsize) # PDF
  133. title = '%s %s %s %s PSTH PDF' % (mseustr, kind, st8combostr, optocombostr)
  134. wintitle(title)
  135. hpdff, hpdfa = plt.subplots(figsize=(figsize[0], 3)) # horizontal, match PSTH height
  136. title = '%s %s %s %s horizontal PSTH PDF' % (mseustr, kind, st8combostr, optocombostr)
  137. wintitle(title)
  138. npdff, npdfa = plt.subplots(figsize=figsize) # normalized PDF
  139. title = '%s %s %s %s PSTH nPDF' % (mseustr, kind, st8combostr, optocombostr)
  140. wintitle(title)
  141. cdff, cdfa = plt.subplots(figsize=figsize)
  142. title = '%s %s %s %s PSTH CDF' % (mseustr, kind, st8combostr, optocombostr)
  143. wintitle(title)
  144. qqf, qqa = plt.subplots(figsize=figsize)
  145. title = '%s %s %s %s PSTH QQ' % (mseustr, kind, st8combostr, optocombostr)
  146. wintitle(title)
  147. # plot PDFs, CDFs and QQ of all st8/opto combinations:
  148. for st8 in st8combo:
  149. c = st82clr[st8]
  150. qq = {}
  151. for opto in optocombo:
  152. psth = psths[st8, opto][::ssx]
  153. npsth = psth / psth.max() # normalize
  154. uvals = np.unique(npsth)
  155. qq[opto] = np.percentile(psth, percentiles)
  156. fb = opto2fb[opto].title()
  157. if st8 == 'none':
  158. label = fb
  159. else:
  160. label = ', '.join([fb, st8.title()])
  161. alpha = opto2alpha[opto]
  162. nbins = intround(np.sqrt(len(psth)))
  163. pdfa.hist(psth, bins=nbins, density=True, histtype='step',
  164. color=c, alpha=alpha, label=label)
  165. hpdfa.hist(psth, bins=nbins, density=True, histtype='step',
  166. color=c, alpha=alpha, label=label, orientation='horizontal')
  167. npdfa.hist(npsth, bins=nbins, density=True, histtype='step',
  168. color=c, alpha=alpha, label=label)
  169. cdfa.hist(npsth, bins=uvals, density=True, cumulative=True, histtype='step',
  170. color=c, alpha=alpha, label=label)
  171. if len(optocombo) == 2:
  172. qqa.scatter(qq[False], qq[True], marker='.', c=c)
  173. #l = a.legend(frameon=False)
  174. #l.set_draggable(True)
  175. pdfa.set_xlabel('Firing rate (spk/s)')
  176. pdfa.set_ylabel('Probability density')
  177. #pdfa.set_yticks([0, 0.5, 1])
  178. pdff.show()
  179. hpdfa.set_xlabel('Probability density')
  180. hpdfa.set_ylabel('Firing rate (spk/s)')
  181. hpdff.show()
  182. npdfa.set_xlabel('Normalized rate')
  183. npdfa.set_ylabel('Probability density')
  184. #pdfa.set_yticks([0, 0.5, 1])
  185. npdff.show()
  186. cdfa.set_xlabel('Normalized rate')
  187. cdfa.set_ylabel('Cumulative probability')
  188. cdfa.set_yticks([0, 0.5, 1])
  189. cdff.show()
  190. if len(optocombo) == 2:
  191. qqmax = np.ceil(np.array(qqa.dataLim).max())
  192. qqxyline = [0, qqmax], [0, qqmax]
  193. qqa.plot(qqxyline[0], qqxyline[1], '--', color='gray', zorder=-1)
  194. qqa.set_xlabel('Feedback rate (spk/s)')
  195. qqa.set_ylabel('Suppression rate (spk/s)')
  196. qqa.set_xlim([0, qqmax])
  197. qqa.set_ylim([0, qqmax])
  198. qqa.set_aspect('equal')
  199. #qqa.set_xticks([0, 0.5, 1])
  200. #qqa.set_yticks([0, 0.5, 1])
  201. qqf.show()
  202. # plot one PSTH PDF, CDF and QQ, for each combination of kind, state and opto,
  203. # each pooled across all mseu:
  204. percentiles = np.arange(0, 100, 2)
  205. figsize = DEFAULTFIGURESIZE
  206. for kind in MVIKINDS:
  207. if EXPTYPE != 'pvmvis':
  208. break # skip these plots
  209. for st8combo in ST8COMBOS:
  210. pdff, pdfa = plt.subplots(figsize=figsize)
  211. title = 'mean %s %s PSTH PDF' % (kind, ' '.join(st8combo))
  212. wintitle(title)
  213. cdff, cdfa = plt.subplots(figsize=figsize)
  214. title = 'mean %s %s PSTH CDF' % (kind, ' '.join(st8combo))
  215. wintitle(title)
  216. qqf, qqa = plt.subplots(figsize=figsize)
  217. title = 'mean %s %s PSTH QQ' % (kind, ' '.join(st8combo))
  218. wintitle(title)
  219. for st8 in st8combo:
  220. c = st82clr[st8]
  221. qq = {}
  222. for opto in OPTOS:
  223. psths, npsths = [], []
  224. for mseustr in mvimseustrs: #mvimseu2exmpli
  225. psth = mviresp.loc[mseustr, kind, st8, opto]['psth']
  226. if np.isnan(psth).any(): # no PSTH for this mseu
  227. continue
  228. #psth = psth[::ssx] # no need to subsample, only calculating a few CDFs
  229. psths.append(psth)
  230. npsths.append(psth / psth.max()) # normalize
  231. psths, npsths = np.concatenate(psths), np.concatenate(npsths)
  232. uvals = np.unique(npsths)
  233. qq[opto] = np.percentile(psths, percentiles)
  234. fb = opto2fb[opto].title()
  235. if st8 == 'none':
  236. label = fb
  237. else:
  238. label = ', '.join([fb, st8.title()])
  239. alpha = opto2alpha[opto]
  240. pdfa.hist(npsths, bins=50, density=True, histtype='step',
  241. color=c, alpha=alpha, label=label)
  242. cdfa.hist(npsths, uvals, density=True, cumulative=True, histtype='step',
  243. color=c, alpha=alpha, label=label)
  244. qqa.scatter(qq[False], qq[True], marker='.', c=c)
  245. #l = a.legend(frameon=False)
  246. #l.set_draggable(True)
  247. pdfa.set_xlabel('Normalized rate')
  248. pdfa.set_ylabel('Probability density')
  249. #pdfa.set_yticks([0, 0.5, 1])
  250. pdff.show()
  251. cdfa.set_xlabel('Normalized rate')
  252. cdfa.set_ylabel('Cumulative probability')
  253. cdfa.set_yticks([0, 0.5, 1])
  254. cdff.show()
  255. qqmax = np.ceil(np.array(qqa.dataLim).max())
  256. qqxyline = [0, qqmax], [0, qqmax]
  257. qqa.plot(qqxyline[0], qqxyline[1], '--', color='gray', zorder=-1)
  258. qqa.set_xlabel('Feedback rate (spk/s)')
  259. qqa.set_ylabel('Suppression rate (spk/s)')
  260. #qqa.set_aspect('equal', 'datalim')
  261. qqa.set_xlim([0, qqmax])
  262. qqa.set_ylim([0, qqmax])
  263. #qqa.set_xticks([0, 0.5, 1])
  264. #print('%s xticks: %r' % (title, qqa.get_xticks()))
  265. #qqa.set_yticks(qqa.get_xticks())
  266. qqf.show()
  267. # scatter plot movie meanrates:
  268. figsize = DEFAULTFIGURESIZE[0]*1.02, DEFAULTFIGURESIZE[1] # tweak to make space for log units
  269. logmin, logmax = -1, 2
  270. logticks = np.array([-1, 0, 1, 2])
  271. #log0min = logmin + 0.05
  272. for kind in ['nat']:#MVIKINDS:
  273. for st8 in ['none']:#ALLST8S:
  274. f, a = plt.subplots(figsize=figsize)
  275. wintitle('opto meanrate movie %s %s' % (kind, st8))
  276. rons, roffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], []
  277. keptmseui = 0 # manually init and increment instead of using enumerate()
  278. for mseustr in mvimseustrs:
  279. meanrate = mviresp.loc[mseustr, kind, st8]['meanrate']
  280. snr = mviresp.loc[mseustr, kind, st8]['snr']
  281. if meanrate.isna().any(): # missing one or both meanrates
  282. continue
  283. if (snr < SNRTHRESH).all(): # neither condition has decent SNR
  284. continue
  285. rons.append(meanrate[True])
  286. roffs.append(meanrate[False])
  287. trialis = mviresp.loc[mseustr, kind, st8]['trialis']
  288. rates = mviresp.loc[mseustr, kind, st8]['rates']
  289. _, pval = ttest_ind(rates[False], rates[True], equal_var=False)
  290. sgnfs.append(pval < SCATTERPTHRESH) # bool
  291. fig1.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
  292. fig1.loc[mseustr]['trialis'] = trialis # save, for joining DataFrames trial-wise
  293. fig1.loc[mseustr]['rates'] = rates # save trial-wise values
  294. if mvimseu2exmpli.get(mseustr) == fig1exmpli:
  295. exmplis.append(keptmseui)
  296. exmplmseustrs.append(mseustr)
  297. else:
  298. normlis.append(keptmseui)
  299. keptmseui += 1 # manually increment
  300. rons = np.asarray(rons)
  301. roffs = np.asarray(roffs)
  302. sgnfs = np.asarray(sgnfs)
  303. # replace off-scale low values with log0min, so the points remain visible:
  304. #pltrons, pltroffs = rons.copy(), roffs.copy()
  305. #pltrons[pltrons <= 10**logmin] = 10**log0min
  306. #pltroffs[pltroffs <= 10**logmin] = 10**log0min
  307. # plot y=x line:
  308. xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
  309. a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
  310. insgnfis, = np.where(sgnfs == False)
  311. sgnfis, = np.where(sgnfs == True)
  312. normlinsgnfis = np.intersect1d(normlis, insgnfis)
  313. normlsgnfis = np.intersect1d(normlis, sgnfis)
  314. if len(normlinsgnfis) > 0:
  315. # plot normal (non-example) insignificant points:
  316. c = desat(st82clr[st8], SGNF2ALPHA[False]) # do manual alpha mixing
  317. a.scatter(rons[normlinsgnfis], roffs[normlinsgnfis], clip_on=False,
  318. marker='.', c='None', edgecolor=c, s=DEFSZ)
  319. if len(normlsgnfis) > 0:
  320. # plot normal (non-example) significant points:
  321. c = desat(st82clr[st8], SGNF2ALPHA[True]) # do manual alpha mixing
  322. a.scatter(rons[normlsgnfis], roffs[normlsgnfis], clip_on=False,
  323. marker='.', c='None', edgecolor=c, s=DEFSZ)
  324. # plot example points, one at a time:
  325. for exmpli, mseustr in zip(exmplis, exmplmseustrs):
  326. marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
  327. c = exmpli2clr[mvimseu2exmpli[mseustr]]
  328. sz = exmpli2sz[mvimseu2exmpli[mseustr]]
  329. lw = exmpli2lw[mvimseu2exmpli[mseustr]]
  330. a.scatter(rons[exmpli], roffs[exmpli], clip_on=False,
  331. marker=marker, c=c, s=sz, lw=lw)
  332. a.set_xlabel('Suppression FR (spk/s)')
  333. a.set_ylabel('Feedback FR (spk/s)')
  334. a.set_xscale('log')
  335. a.set_yscale('log')
  336. a.set_xlim(10**logmin, 10**logmax)
  337. a.set_ylim(10**logmin, 10**logmax)
  338. a.set_xticks(10.0**logticks)
  339. a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
  340. a.minorticks_off()
  341. axes_disable_scientific(a)
  342. a.set_aspect('equal')
  343. a.spines['left'].set_position(('outward', 4))
  344. a.spines['bottom'].set_position(('outward', 4))
  345. #t, p = ttest_rel(rons, roffs) # paired t-test
  346. #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
  347. #mu = rons.mean(), roffs.mean()
  348. #txt = '$\mathregular{\mu=%.1f, %.1f}$' % mu
  349. #a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))
  350. # plot mean firing rate PDF:
  351. f, a = plt.subplots(figsize=figsize)
  352. wintitle('opto meanrate PDF movie %s %s' % (kind, st8))
  353. bins = np.logspace(-1, 2, num=25, base=10) # end inclusive
  354. a.hist(roffs, bins=bins, histtype='step', density=False, color='k',
  355. label='Feedback', clip_on=False)
  356. a.hist(rons, bins=bins, histtype='step', density=False, color=optoblue,
  357. label='V1 Suppression', clip_on=False)
  358. a.set_xlabel('Mean firing rate (spk/s)')
  359. a.set_ylabel('Cell count')
  360. a.set_xscale('log')
  361. a.set_xlim(10**logmin, 10**logmax)
  362. a.set_xticks(10.0**logticks)
  363. a.minorticks_off()
  364. axes_disable_scientific(a)
  365. #l = a.legend(frameon=False)
  366. #l.set_draggable(True)
  367. # scatter plot movie burst ratio:
  368. figsize = DEFAULTFIGURESIZE[0]*1.04, DEFAULTFIGURESIZE[1] # tweak to make space for log units
  369. logmin, logmax = -3.55, 0
  370. if EXPTYPE == 'ntsrmvis':
  371. logmin, logmax = -3.83, 0
  372. logticks = np.array([-3, -2, -1, 0])
  373. for kind in ['nat']:#MVIKINDS:
  374. for st8 in ['none']:#ALLST8S:
  375. f, a = plt.subplots(figsize=figsize)
  376. wintitle('opto burst ratio movie %s %s' % (kind, st8))
  377. brons, broffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], []
  378. keptmseui = 0 # manually init and increment instead of using enumerate()
  379. for mseustr in mvimseustrs:
  380. br = mviresp.loc[mseustr, kind, st8]['meanburstratio']
  381. snr = mviresp.loc[mseustr, kind, st8]['snr']
  382. if br.isna().any(): # missing for at least one opto condition
  383. continue
  384. if (snr < SNRTHRESH).all(): # neither condition has decent SNR
  385. continue
  386. brons.append(br[True])
  387. broffs.append(br[False])
  388. burstratios = mviresp.loc[mseustr, kind, st8]['burstratios']
  389. _, pval = ttest_ind(burstratios[False], burstratios[True], equal_var=False)
  390. sgnfs.append(pval < SCATTERPTHRESH) # bool
  391. fig1.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
  392. fig1.loc[mseustr]['burstratios'] = burstratios # save trial-wise values
  393. if mvimseu2exmpli.get(mseustr) == fig1exmpli:
  394. exmplis.append(keptmseui)
  395. exmplmseustrs.append(mseustr)
  396. else:
  397. normlis.append(keptmseui)
  398. keptmseui += 1 # manually increment
  399. brons, broffs = np.asarray(brons), np.asarray(broffs)
  400. sgnfs = np.asarray(sgnfs)
  401. # plot y=x line:
  402. xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
  403. a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
  404. insgnfis, = np.where(sgnfs == False)
  405. sgnfis, = np.where(sgnfs == True)
  406. normlinsgnfis = np.intersect1d(normlis, insgnfis)
  407. normlsgnfis = np.intersect1d(normlis, sgnfis)
  408. if len(normlinsgnfis) > 0:
  409. # plot normal (non-example) insignificant points:
  410. c = desat(st82clr[st8], SGNF2ALPHA[False]) # do manual alpha mixing
  411. a.scatter(brons[normlinsgnfis], broffs[normlinsgnfis], clip_on=True,
  412. marker='.', c='None', edgecolor=c, s=DEFSZ) # clip_on=False fails
  413. if len(normlsgnfis) > 0:
  414. # plot normal (non-example) significant points:
  415. c = desat(st82clr[st8], SGNF2ALPHA[True]) # do manual alpha mixing
  416. a.scatter(brons[normlsgnfis], broffs[normlsgnfis], clip_on=True,
  417. marker='.', c='None', edgecolor=c, s=DEFSZ) # clip_on=False fails
  418. # plot example points, one at a time:
  419. for exmpli, mseustr in zip(exmplis, exmplmseustrs):
  420. marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
  421. c = exmpli2clr[mvimseu2exmpli[mseustr]]
  422. sz = exmpli2sz[mvimseu2exmpli[mseustr]]
  423. lw = exmpli2lw[mvimseu2exmpli[mseustr]]
  424. a.scatter(brons[exmpli], broffs[exmpli], clip_on=False,
  425. marker=marker, c=c, s=sz, lw=lw)
  426. a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
  427. a.set_ylabel('Feedback BR')
  428. a.set_xscale('log')
  429. a.set_yscale('log')
  430. a.set_xlim(10**logmin, 10**logmax)
  431. a.set_ylim(10**logmin, 10**logmax)
  432. a.set_xticks(10.0**logticks)
  433. a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
  434. a.minorticks_off()
  435. axes_disable_scientific(a)
  436. a.set_aspect('equal')
  437. a.spines['left'].set_position(('outward', 4))
  438. a.spines['bottom'].set_position(('outward', 4))
  439. #t, p = ttest_rel(brons, broffs) # paired t-test
  440. #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
  441. # scatter plot movie PSTH sparseness:
  442. figsize = DEFAULTFIGURESIZE
  443. linmin, linmax, linstep = 0, 1, 0.5
  444. ticks = np.arange(intround(linmin), intround(linmax)+linstep, linstep)
  445. for kind in ['nat']:#MVIKINDS:
  446. for st8 in ['none']:#ALLST8S:
  447. f, a = plt.subplots(figsize=figsize)
  448. wintitle('opto sparseness %s %s' % (kind, st8))
  449. sparsons, sparsoffs, exmplis, exmplmseustrs, normlis = [], [], [], [], []
  450. keptmseui = 0 # manually init and increment instead of using enumerate()
  451. for mseustr in mvimseustrs:
  452. spars = mviresp.loc[mseustr, kind, st8]['spars']
  453. snr = mviresp.loc[mseustr, kind, st8]['snr']
  454. if spars.isna().any(): # missing for at least one opto condition
  455. continue
  456. if (snr < SNRTHRESH).all(): # neither condition has decent SNR
  457. continue
  458. sparsons.append(spars[True])
  459. sparsoffs.append(spars[False])
  460. fig1.loc[mseustr, 'spars'] = spars[False], spars[True] # save
  461. if mvimseu2exmpli.get(mseustr) == fig1exmpli:
  462. exmplis.append(keptmseui)
  463. exmplmseustrs.append(mseustr)
  464. else:
  465. normlis.append(keptmseui)
  466. keptmseui += 1 # manually increment
  467. sparsons, sparsoffs = np.asarray(sparsons), np.asarray(sparsoffs)
  468. # plot y=x line:
  469. xyline = [linmin, linmax], [linmin, linmax]
  470. a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
  471. # plot normal (non-example) points:
  472. a.scatter(sparsons[normlis], sparsoffs[normlis], clip_on=False,
  473. marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
  474. # plot example points, one at a time:
  475. for exmpli, mseustr in zip(exmplis, exmplmseustrs):
  476. marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
  477. c = exmpli2clr[mvimseu2exmpli[mseustr]]
  478. sz = exmpli2sz[mvimseu2exmpli[mseustr]]
  479. lw = exmpli2lw[mvimseu2exmpli[mseustr]]
  480. a.scatter(sparsons[exmpli], sparsoffs[exmpli], clip_on=False,
  481. marker=marker, c=c, s=sz, lw=lw)
  482. a.set_xlabel('Suppression spars') # keep it short to maximize space for axes
  483. a.set_ylabel('Feedback spars')
  484. a.set_xlim(linmin, linmax)
  485. a.set_ylim(linmin, linmax)
  486. a.set_xticks(ticks)
  487. a.set_yticks(ticks)
  488. a.set_aspect('equal')
  489. a.spines['left'].set_position(('outward', 4))
  490. a.spines['bottom'].set_position(('outward', 4))
  491. #t, p = ttest_rel(sparsons, sparsoffs) # paired t-test
  492. #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
  493. # scatter plot movie PSTH reliability:
  494. figsize = DEFAULTFIGURESIZE[0]*0.98, DEFAULTFIGURESIZE[1] # tweak to match others
  495. linmin, linmax = 0, 0.52
  496. ticks = 0, 0.5
  497. for kind in ['nat']:#MVIKINDS:
  498. for st8 in ['none']:#ALLST8S:
  499. f, a = plt.subplots(figsize=figsize)
  500. wintitle('opto reliability %s %s' % (kind, st8))
  501. relons, reloffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], []
  502. keptmseui = 0 # manually init and increment instead of using enumerate()
  503. for mseustr in mvimseustrs:
  504. rel = mviresp.loc[mseustr, kind, st8]['rel']
  505. snr = mviresp.loc[mseustr, kind, st8]['snr']
  506. if rel.isna().any(): # missing for at least one opto condition
  507. continue
  508. if (snr < SNRTHRESH).all(): # neither condition has decent SNR
  509. continue
  510. relons.append(rel[True])
  511. reloffs.append(rel[False])
  512. # collect rhos from all trial pairs, reduce to single rho per trial, use
  513. # to estimate significance of individual units:
  514. rhos = mviresp.loc[mseustr, kind, st8]['rhos']
  515. fullrhos, meanrhos = {}, {}
  516. for opto in OPTOS:
  517. npairs = len(rhos[opto])
  518. # from quadratic formula solution to npairs = ntrials*(ntrials-1) / 2
  519. ntrials = (1 + np.sqrt(1 + 8*npairs)) / 2
  520. assert np.isclose(ntrials, int(ntrials)) # make sure this float is int
  521. ntrials = int(ntrials)
  522. uti = np.triu_indices(ntrials, k=1) # upper indices, excludes diagonal
  523. lti = np.tril_indices(ntrials, k=-1) # lower indices, excludes diagonal
  524. # empty full corr matrix, nans will remain on diagonal and will be ignored:
  525. fullrhos[opto] = np.tile(np.nan, (ntrials, ntrials))
  526. fullrhos[opto][uti] = rhos[opto] # fill the upper triangle
  527. # copy upper to lower triangle,
  528. # see https://stackoverflow.com/a/42209263/2020363:
  529. fullrhos[opto][lti] = fullrhos[opto].T[lti]
  530. meanrhos[opto] = np.nanmean(fullrhos[opto], axis=1) # along either axis is fine
  531. _, pval = ttest_ind(meanrhos[False], meanrhos[True], equal_var=False)
  532. sgnfs.append(pval < SCATTERPTHRESH) # bool
  533. fig1.loc[mseustr, 'rel'] = rel[False], rel[True] # save
  534. if mvimseu2exmpli.get(mseustr) == fig1exmpli:
  535. exmplis.append(keptmseui)
  536. exmplmseustrs.append(mseustr)
  537. else:
  538. normlis.append(keptmseui)
  539. keptmseui += 1 # manually increment
  540. relons, reloffs = np.asarray(relons), np.asarray(reloffs)
  541. sgnfs = np.asarray(sgnfs)
  542. # plot y=x line:
  543. xyline = [linmin, linmax], [linmin, linmax]
  544. a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
  545. insgnfis, = np.where(sgnfs == False)
  546. sgnfis, = np.where(sgnfs == True)
  547. normlinsgnfis = np.intersect1d(normlis, insgnfis)
  548. normlsgnfis = np.intersect1d(normlis, sgnfis)
  549. if len(normlinsgnfis) > 0:
  550. # plot normal (non-example) insignificant points:
  551. c = desat(st82clr[st8], SGNF2ALPHA[False]) # do manual alpha mixing
  552. a.scatter(relons[normlinsgnfis], reloffs[normlinsgnfis], clip_on=False,
  553. marker='.', c='None', edgecolor=c, s=DEFSZ)
  554. if len(normlsgnfis) > 0:
  555. # plot normal (non-example) significant points:
  556. c = desat(st82clr[st8], SGNF2ALPHA[True]) # do manual alpha mixing
  557. a.scatter(relons[normlsgnfis], reloffs[normlsgnfis], clip_on=False,
  558. marker='.', c='None', edgecolor=c, s=DEFSZ)
  559. # plot example points, one at a time:
  560. for exmpli, mseustr in zip(exmplis, exmplmseustrs):
  561. marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
  562. c = exmpli2clr[mvimseu2exmpli[mseustr]]
  563. sz = exmpli2sz[mvimseu2exmpli[mseustr]]
  564. lw = exmpli2lw[mvimseu2exmpli[mseustr]]
  565. a.scatter(relons[exmpli], reloffs[exmpli], clip_on=False,
  566. marker=marker, c=c, s=sz, lw=lw)
  567. a.set_xlabel('Suppression rel') # keep it short to maximize space for axes
  568. a.set_ylabel('Feedback rel')
  569. a.set_xlim(linmin, linmax)
  570. a.set_ylim(linmin, linmax)
  571. a.set_xticks(ticks)
  572. a.set_yticks(ticks)
  573. a.set_aspect('equal')
  574. a.spines['left'].set_position(('outward', 4))
  575. a.spines['bottom'].set_position(('outward', 4))
  576. #t, p = ttest_rel(relons, reloffs, nan_policy='omit') # paired t-test
  577. #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
  578. # plot movie SNR PDFs, show SNRTHRESH:
  579. figsize = DEFAULTFIGURESIZE
  580. bins = np.arange(0, 0.5, 0.005)
  581. #bins = np.logspace(-2, -0.2, 50, base=10)
  582. for kind in ['nat']:#MVIKINDS:
  583. for st8 in ['none']:#ALLST8S:
  584. for opto in [None, False, True]:
  585. if opto is not None:
  586. snrs = mviresp['snr'][(slice(None), kind, st8, opto)].dropna() # one val per mseu
  587. else:
  588. snrs = mviresp['snr'][(slice(None), kind, st8)].dropna() # two vals per mseu?
  589. snrs = np.float64(snrs) # convert from object array
  590. f, a = plt.subplots(figsize=figsize)
  591. titlestr = 'SNR PDF %s %s' % (kind, st8)
  592. if opto is not None:
  593. titlestr += ' %s' % opto2fb[opto]
  594. wintitle(titlestr, f)
  595. c = st82clr[st8]
  596. a.hist(snrs, bins=bins, density=True, color=c,
  597. alpha=opto2alpha[opto], label=opto2fb[opto])
  598. a.axvline(x=SNRTHRESH, ls='-', marker='', c='gray')#, zorder=-np.inf)
  599. txt = 'n=%d' % len(snrs)
  600. a.add_artist(AnchoredText(txt, loc='upper right', frameon=False,
  601. prop=dict(color=c, alpha=opto2alpha[opto])))
  602. #a.set_xscale('log')
  603. #l = a.legend(frameon=False)
  604. #l.set_draggable(True)
  605. a.set_xlabel('SNR')
  606. if opto is not None:
  607. a.set_xlabel('%s SNR' % opto2fb[opto].title())
  608. a.set_ylabel('Probability density')
  609. # scatter plot movie SNR:
  610. figsize = DEFAULTFIGURESIZE
  611. linmin, linmax = 0, 0.6
  612. ticks = 0, 0.5
  613. for kind in ['nat']:#MVIKINDS:
  614. for st8 in ['none']:#ALLST8S:
  615. f, a = plt.subplots(figsize=figsize)
  616. wintitle('opto SNR %s %s' % (kind, st8))
  617. snrons, snroffs, exmplis, exmplmseustrs, normlis = [], [], [], [], []
  618. keptmseui = 0 # manually init and increment instead of using enumerate()
  619. for mseustr in mvimseustrs:
  620. snr = mviresp.loc[mseustr, kind, st8]['snr']
  621. if snr.isna().any(): # missing for at least one opto condition
  622. continue
  623. if (snr < SNRTHRESH).all(): # neither condition has decent SNR
  624. continue
  625. snrons.append(snr[True])
  626. snroffs.append(snr[False])
  627. fig1.loc[mseustr, 'snr'] = snr[False], snr[True] # save
  628. if mvimseu2exmpli.get(mseustr) == fig1exmpli:
  629. exmplis.append(keptmseui)
  630. exmplmseustrs.append(mseustr)
  631. else:
  632. normlis.append(keptmseui)
  633. keptmseui += 1 # manually increment
  634. snrons, snroffs = np.asarray(snrons), np.asarray(snroffs)
  635. # plot y=x line:
  636. xyline = [linmin, linmax], [linmin, linmax]
  637. a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
  638. # plot normal (non-example) points:
  639. a.scatter(snrons[normlis], snroffs[normlis], clip_on=False,
  640. marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
  641. # plot example points, one at a time:
  642. for exmpli, mseustr in zip(exmplis, exmplmseustrs):
  643. marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
  644. c = exmpli2clr[mvimseu2exmpli[mseustr]]
  645. sz = exmpli2sz[mvimseu2exmpli[mseustr]]
  646. lw = exmpli2lw[mvimseu2exmpli[mseustr]]
  647. a.scatter(snrons[exmpli], snroffs[exmpli], clip_on=False,
  648. marker=marker, c=c, s=sz, lw=lw)
  649. a.set_xlabel('Suppression SNR')
  650. a.set_ylabel('Feedback SNR')
  651. a.set_xlim(linmin, linmax)
  652. a.set_ylim(linmin, linmax)
  653. a.set_xticks(ticks)
  654. a.set_yticks(ticks)
  655. a.set_aspect('equal')
  656. a.spines['left'].set_position(('outward', 4))
  657. a.spines['bottom'].set_position(('outward', 4))
  658. #t, p = ttest_rel(snrons, snroffs, nan_policy='omit') # paired t-test
  659. #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
  660. # plot movie PSTH peak width distributions:
  661. ## TODO: these width measures should maybe be normalized separately for each unit,
  662. ## so that the effect of feedback suppression isn't drowned out by the variability of
  663. ## PSTH peak widths across units
  664. normalizepeakwidths = True
  665. figsize = DEFAULTFIGURESIZE
  666. for kind in ['nat']:#MVIKINDS:
  667. for st8 in ['none']:#ALLST8S:
  668. f, a = plt.subplots(figsize=figsize)
  669. titlestr = ('opto peak width PDF %s %s PSTHBASELINEMEDIANX=%d PSTHPEAKMINTHRESH=%d'
  670. % (kind, st8, PSTHBASELINEMEDIANX, PSTHPEAKMINTHRESH))
  671. if normalizepeakwidths:
  672. titlestr += " normed"
  673. wintitle(titlestr, f)
  674. pkwons, pkwoffs = [], []
  675. for mseustr in mvimseustrs:
  676. pkws = mviresp.loc[mseustr, kind, st8]['pkws']
  677. snr = mviresp.loc[mseustr, kind, st8]['snr']
  678. if pkws.isna().any():
  679. continue # skip this mseu
  680. if (snr < SNRTHRESH).all(): # neither condition has decent SNR
  681. continue
  682. if normalizepeakwidths:
  683. maxpkw = np.concatenate(pkws.values).max() # max peak width for this mseu
  684. pkws = pkws / maxpkw # normalize all peak widths by maxpkw
  685. pkwons.append(pkws[True])
  686. pkwoffs.append(pkws[False])
  687. c = st82clr[st8]
  688. pkwons, pkwoffs = np.concatenate(pkwons), np.concatenate(pkwoffs)
  689. for opto, pkws, loc in zip(OPTOS, [pkwoffs, pkwons], ['upper right', 'upper left']):
  690. a.hist(pkws, bins=50, density=True, histtype='step', color=c,
  691. alpha=opto2alpha[opto], label=opto2fb[opto])
  692. txt = 'n=%d' % len(pkws)
  693. a.add_artist(AnchoredText(txt, loc=loc, frameon=False,
  694. prop=dict(color=c, alpha=opto2alpha[opto])))
  695. a.set_xscale('log')
  696. #a.set_yscale('log')
  697. axes_disable_scientific(a)
  698. #l = a.legend(frameon=False)
  699. #l.set_draggable(True)
  700. a.set_xlabel('PSTH peak width (s)')
  701. a.set_ylabel('Probability density')
  702. #t, p = ttest_rel(pkwoffs, pkwons, nan_policy='omit') # paired t-test
  703. #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
  704. # scatter plot movie PSTH mean peak widths:
  705. figsize = DEFAULTFIGURESIZE[0]*1.05, DEFAULTFIGURESIZE[1] # tweak to make space for log units
  706. linmin, linmax = 0.025, 0.2
  707. ticks = 0.05, 0.1, 0.2
  708. for kind in ['nat']:#MVIKINDS:
  709. for st8 in ['none']:#ALLST8S:
  710. f, a = plt.subplots(figsize=figsize)
  711. wintitle('opto mean peak width %s %s' % (kind, st8))
  712. meanpkwons, meanpkwoffs, exmplis, exmplmseustrs, normlis = [], [], [], [], []
  713. keptmseui = 0 # manually init and increment instead of using enumerate()
  714. for mseustr in mvimseustrs:
  715. meanpkw = mviresp.loc[mseustr, kind, st8]['meanpkw']
  716. snr = mviresp.loc[mseustr, kind, st8]['snr']
  717. if meanpkw.isna().any(): # missing for at least one opto condition
  718. continue
  719. if (snr < SNRTHRESH).all(): # neither condition has decent SNR
  720. continue
  721. meanpkwons.append(meanpkw[True])
  722. meanpkwoffs.append(meanpkw[False])
  723. fig1.loc[mseustr, 'meanpkw'] = meanpkw[False], meanpkw[True] # save
  724. if mvimseu2exmpli.get(mseustr) == fig1exmpli:
  725. exmplis.append(keptmseui)
  726. exmplmseustrs.append(mseustr)
  727. else:
  728. normlis.append(keptmseui)
  729. keptmseui += 1 # manually increment
  730. meanpkwons, meanpkwoffs = np.asarray(meanpkwons), np.asarray(meanpkwoffs)
  731. # plot y=x line:
  732. xyline = [linmin, linmax], [linmin, linmax]
  733. a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
  734. # plot normal (non-example) points:
  735. a.scatter(meanpkwons[normlis], meanpkwoffs[normlis], clip_on=False,
  736. marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
  737. # plot example points, one at a time:
  738. for exmpli, mseustr in zip(exmplis, exmplmseustrs):
  739. marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
  740. c = exmpli2clr[mvimseu2exmpli[mseustr]]
  741. sz = exmpli2sz[mvimseu2exmpli[mseustr]]
  742. lw = exmpli2lw[mvimseu2exmpli[mseustr]]
  743. a.scatter(meanpkwons[exmpli], meanpkwoffs[exmpli], clip_on=False,
  744. marker=marker, c=c, s=sz, lw=lw)
  745. a.set_xscale('log', basex=2)
  746. a.set_yscale('log', basey=2)
  747. a.set_xlabel('Suppression pkw (s)') # keep it short to maximize space for axes
  748. a.set_ylabel('Feedback pkw (s)')
  749. a.set_xlim(linmin, linmax)
  750. a.set_ylim(linmin, linmax)
  751. a.set_xticks(ticks)
  752. a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
  753. a.minorticks_off()
  754. axes_disable_scientific(a)
  755. a.set_aspect('equal')
  756. a.spines['left'].set_position(('outward', 4))
  757. a.spines['bottom'].set_position(('outward', 4))
  758. a.xaxis.set_major_formatter(ScalarFormatter())
  759. a.yaxis.set_major_formatter(ScalarFormatter())
  760. a.xaxis.set_major_formatter(FormatStrFormatter("%.2g"))
  761. a.yaxis.set_major_formatter(FormatStrFormatter("%.2g"))
  762. #t, p = ttest_rel(meanpkwoffs, meanpkwons, nan_policy='omit') # paired t-test
  763. #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
  764. """Figure 1S2 FMI plots"""
  765. # scatter plot FMI of all movie measures vs meanrate FMI:
  766. figsize = DEFAULTFIGURESIZE
  767. linmin, linmax, linstep = -1, 1, 1
  768. ticks = np.arange(linmin, linmax+linstep, linstep)
  769. for kind in ['nat']:#MVIKINDS:
  770. for st8 in ['none']:#ALLST8S:
  771. for measure in modmeasuresnoblank:
  772. if measure.startswith('meanrate'):
  773. continue # don't bother plotting meanrate* against meanrate
  774. axislabel = measure2axislabel[measure]
  775. rfmis, msrfmis = [], [] # rate FMIs, measure FMIs
  776. exmplis, exmplmseustrs, normlis = [], [], []
  777. keptmseui = 0 # manually init and increment instead of using enumerate()
  778. for mseustr in mvimseustrs:
  779. rfmi = mviFMI.loc[mseustr, st8]['meanrate']
  780. msrfmi = mviFMI.loc[mseustr, st8][measure]
  781. if pd.isna(rfmi) or pd.isna(msrfmi): # missing at least one FMI value
  782. continue
  783. rfmis.append(rfmi)
  784. msrfmis.append(msrfmi)
  785. if mvimseu2exmpli.get(mseustr) == fig1exmpli:
  786. exmplis.append(keptmseui)
  787. exmplmseustrs.append(mseustr)
  788. else:
  789. normlis.append(keptmseui)
  790. keptmseui += 1 # manually increment
  791. rfmis, msrfmis = np.asarray(rfmis), np.asarray(msrfmis)
  792. ## scatter plot FMIs:
  793. f, a = plt.subplots(figsize=figsize)
  794. wintitle('opto FMI %s vs FMI rate %s %s' % (measure, kind, st8))
  795. # plot x=0 and y=0 lines:
  796. a.axhline(y=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
  797. a.axvline(x=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
  798. # plot normal (non-example) points:
  799. a.scatter(rfmis[normlis], msrfmis[normlis], clip_on=False,
  800. marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
  801. # plot example points, one at a time:
  802. for exmpli, mseustr in zip(exmplis, exmplmseustrs):
  803. marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
  804. c = exmpli2clr[mvimseu2exmpli[mseustr]]
  805. sz = exmpli2sz[mvimseu2exmpli[mseustr]]
  806. lw = exmpli2lw[mvimseu2exmpli[mseustr]]
  807. a.scatter(rfmis[exmpli], msrfmis[exmpli], clip_on=False,
  808. marker=marker, c=c, s=sz, lw=lw)
  809. # get fname of appropriate LMM .cvs file:
  810. if measure == 'meanburstratio': # use Steffen's LMM linregress fit, for fig1S2
  811. fname = os.path.join(PAPERPATH, 'stats', 'figure_1S2c_coefs.csv')
  812. elif measure == 'spars': # use Steffen's LMM linregress fit, for fig1S2
  813. fname = os.path.join(PAPERPATH, 'stats', 'figure_1S2d_coefs.csv')
  814. elif measure == 'rel': # use Steffen's LMM linregress fit, for fig1S2
  815. fname = os.path.join(PAPERPATH, 'stats', 'figure_1S2e_coefs.csv')
  816. elif measure == 'snr': # use Steffen's LMM linregress fit, for fig1S2
  817. fname = os.path.join(PAPERPATH, 'stats', 'figure_1S2f_coefs.csv')
  818. elif measure == 'meanpkw': # use Steffen's LMM linregress fit, for fig1S2
  819. fname = os.path.join(PAPERPATH, 'stats', 'figure_1S2g_coefs.csv')
  820. else:
  821. print("WARNING: No LMM stats for %s measure" % measure)
  822. fname = None
  823. # fetch LMM linregress fit params from .csv:
  824. if fname:
  825. try:
  826. df = pd.read_csv(fname)
  827. foundregression = True
  828. except FileNotFoundError:
  829. print('Missing file: %s' % fname)
  830. foundregression = False
  831. if foundregression:
  832. mm = df['slope'][0]
  833. b = df['intercept'][0]
  834. x = np.array([rfmis.min(), rfmis.max()])
  835. y = mm * x + b
  836. a.plot(x, y, '-', color='red') # plot linregress fit
  837. a.set_xlabel('FR FMI')
  838. if axislabel.islower():
  839. axislabel = axislabel.capitalize()
  840. ylabel = '%s FMI' % axislabel
  841. a.set_ylabel(ylabel)
  842. a.set_xlim(linmin, linmax)
  843. a.set_ylim(linmin, linmax)
  844. a.set_xticks(ticks)
  845. a.set_yticks(ticks)
  846. a.set_aspect('equal')
  847. a.spines['left'].set_position(('outward', 4))
  848. a.spines['bottom'].set_position(('outward', 4))
  849. """Some figure 4S1 plots"""
  850. # scatter plot meanrates of first 2 and last 2 s of movies, limited to first 120 trials
  851. # per condition, as per standard grating experiments:
  852. figsize = DEFAULTFIGURESIZE[0]*1.04, DEFAULTFIGURESIZE[1] # tweak to make space for log units
  853. logmin, logmax = -2, 2
  854. logticks = np.array([-2, 0, 2])
  855. #log0min = logmin + 0.05
  856. kind = 'nat'
  857. st8 = 'none'
  858. for meanratestr, rangestr in zip(['meanrate02', 'meanrate35'], ['0-2', '3-5']):
  859. f, a = plt.subplots(figsize=figsize)
  860. wintitle('opto meanrate movie t=%s s %d trials %s %s' % (rangestr, NTRIALSMVIGRT, kind, st8))
  861. rons, roffs, sgnfs, exmplis, exmplmseustrs, normlis = [], [], [], [], [], []
  862. keptmseui = 0 # manually init and increment instead of using enumerate()
  863. for mseustr in mvimseustrs:
  864. meanrate = mviresp.loc[mseustr, kind, st8][meanratestr]
  865. snr = mviresp.loc[mseustr, kind, st8]['snr']
  866. if meanrate.isna().any(): # missing one or both meanrates
  867. continue
  868. if (snr < SNRTHRESH).all(): # neither condition has decent SNR on the full range
  869. continue
  870. rons.append(meanrate[True])
  871. roffs.append(meanrate[False])
  872. ratesstr = meanratestr.split('mean')[1]+'s'
  873. rates = mviresp.loc[mseustr, kind, st8][ratesstr]
  874. _, pval = ttest_ind(rates[False], rates[True], equal_var=False)
  875. sgnfs.append(pval < SCATTERPTHRESH) # bool
  876. fig1.loc[mseustr, meanratestr] = meanrate[False], meanrate[True] # save
  877. fig1.loc[mseustr][ratesstr] = rates # save trial-wise values
  878. if mvimseu2exmpli.get(mseustr) == fig1exmpli:
  879. exmplis.append(keptmseui)
  880. exmplmseustrs.append(mseustr)
  881. else:
  882. normlis.append(keptmseui)
  883. keptmseui += 1 # manually increment
  884. rons = np.asarray(rons)
  885. roffs = np.asarray(roffs)
  886. sgnfs = np.asarray(sgnfs)
  887. # replace off-scale low values with log0min, so the points remain visible:
  888. #pltrons, pltroffs = rons.copy(), roffs.copy()
  889. #pltrons[pltrons <= 10**logmin] = 10**log0min
  890. #pltroffs[pltroffs <= 10**logmin] = 10**log0min
  891. # plot y=x line:
  892. xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
  893. a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
  894. insgnfis, = np.where(sgnfs == False)
  895. sgnfis, = np.where(sgnfs == True)
  896. normlinsgnfis = np.intersect1d(normlis, insgnfis)
  897. normlsgnfis = np.intersect1d(normlis, sgnfis)
  898. # plot normal (non-example) insignificant points:
  899. c = desat(st82clr[st8], SGNF2ALPHA[False]) # do manual alpha mixing
  900. a.scatter(rons[normlinsgnfis], roffs[normlinsgnfis], clip_on=False,
  901. marker='.', c='None', edgecolor=c, s=DEFSZ)
  902. # plot normal (non-example) significant points:
  903. c = desat(st82clr[st8], SGNF2ALPHA[True]) # do manual alpha mixing
  904. a.scatter(rons[normlsgnfis], roffs[normlsgnfis], clip_on=False,
  905. marker='.', c='None', edgecolor=c, s=DEFSZ)
  906. # plot normal (non-example) points:
  907. #a.scatter(rons[normlis], roffs[normlis], clip_on=False,
  908. # marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
  909. # plot example points, one at a time:
  910. for exmpli, mseustr in zip(exmplis, exmplmseustrs):
  911. marker = exmpli2mrk[mvimseu2exmpli[mseustr]]
  912. c = exmpli2clr[mvimseu2exmpli[mseustr]]
  913. sz = exmpli2sz[mvimseu2exmpli[mseustr]]
  914. lw = exmpli2lw[mvimseu2exmpli[mseustr]]
  915. a.scatter(rons[exmpli], roffs[exmpli], marker=marker, clip_on=False,
  916. c=c, s=sz, lw=lw)
  917. a.set_xlabel('Suppression FR (spk/s)')
  918. a.set_ylabel('Feedback FR (spk/s)')
  919. a.set_xscale('log')
  920. a.set_yscale('log')
  921. a.set_xlim(10**logmin, 10**logmax)
  922. a.set_ylim(10**logmin, 10**logmax)
  923. a.set_xticks(10.0**logticks)
  924. a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
  925. a.minorticks_off()
  926. axes_disable_scientific(a)
  927. a.set_aspect('equal')
  928. a.spines['left'].set_position(('outward', 4))
  929. a.spines['bottom'].set_position(('outward', 4))
  930. #t, p = ttest_rel(rons, roffs) # paired t-test
  931. #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
  932. #mu = rons.mean(), roffs.mean()
  933. #txt = '$\mathregular{\mu=%.1f, %.1f}$' % mu
  934. #a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))