1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786 |
- """Collect rasters of responses for all LGN opto experiments as a function of stimulus type
- (movies, gratings), run state, and opto state. Then, calculate various measures and save
- results to DataFrames. Some variables (DataFrames and mseu string lists) can be optionally
- saved to and loaded from .pickle files to avoid long and unnecessary recalculation times.
- After connecting to server with `djd`, run this script at the IPython/Jupyter command line;
- the `-i` runs scripts in IPython's namespace instead of an empty one:
- `run -i main.py`
- `pvmvis` (default), `ntsrmvis` or `negntsrmvis` can be provided as an `--exptype=` argument
- to this script, to load and later plot PV-Cre (most of the paper), Ntsr1-Cre (Fig1S4),
- or Ntsr negative (Fig1S5) data, respectively.
- This script can also be run without a server connection, if you choose to load from pickles.
- """
- import sys
- import argparse
- from numpy import pi
- from numpy.random import choice
- import scipy.stats
- from scipy.stats import (linregress, ttest_rel, ttest_ind, ttest_1samp, wilcoxon,
- ks_2samp)
- from scipy.stats.mstats import gmean
- from scipy.optimize import least_squares
- #import statsmodels.formula.api as smf
- import matplotlib.pyplot as plt
- from matplotlib.offsetbox import AnchoredText
- from matplotlib.ticker import ScalarFormatter, FormatStrFormatter
- import seaborn as sns
- import pandas as pd
- #from pycircstat.tests import kuiper
- from collections import OrderedDict as odict
- import expio
- from djd.util import (seq2sql, intround, split_tranges, wrap_raster, issubset, key2mseustr,
- mseustr2key, key2msestr)
- from djd.plot import wintitle, simpletraster, cf, saveall, simple_plot_rfs
- from djd.signal import (raster2psth, raster2freqcomp, sparseness, reliability, snr_baden2016,
- get_psth_peaks_gac, get_trialmat, pairwisecorr, rfs_for_cdata)
- from djd.model import fit_MUA_RFs
- from djd.stats import vector_OSI, percentile_ci, rayleigh_test, niell_DSI
- from djd.model import sum_of_gaussians
- from djd.eye import i_clean
- from util import (get_exps, fitmodel, desat, get_max_snr, ms2msstr,
- msu2msustr, mseustr2msustr, mseustrs2msustrs, mseustr2mstr, mseustrs2mstrs,
- mseustr2msestr, mseustrs2msestrs, findmse,
- load, load_all, save, export2csv, axes_disable_scientific, linear_loss,
- residual_rsquared)
- #import rf # deprecated
- ## Define constants:
- EXPTYPE = 'pvmvis' # pvmvis, ntsrmvis, negntsrmvis, pgnpvmvis
- if __name__ == '__main__':
- # parse command-line arguments:
- parser = argparse.ArgumentParser(description='Calculate script for natfeedback paper')
- parser.add_argument("--exptype", help="pvmvis, ntsrmvis, negntsrmvis, pgnpvmvis")
- args = parser.parse_args()
- EXPTYPE = args.exptype or EXPTYPE
- print('*** EXPTYPE: %r' % EXPTYPE)
- PAPERPATH = '/home/mspacek/blab/natfb/paper'
- VERBOSE = False # only used for gac at the moment
- RATETHRESH = 0.01 # mean firing rate threshold for unit inclusion, across all trial types, Hz
- SNRTHRESH = 0.015
- #RSQTHRESH = 0.2
- # p-value thresh for significance of each point's deviance from diagonal in meanrate, burst
- # ratio, and reliability scatter plots:
- SCATTERPTHRESH = 0.05
- SGNF2ALPHA = {True: 1, False: 0.4}
- # visual drive unit inclusion criteria for gratings, ignores st8 and opto:
- VDNCRIT = 1 # min num points in tuning curve that are at least VDZCRIT sems away from spont rate
- VDZCRIT = 2.5 # 1.95, 2.5, 3.2
- ## TODO: these need to be optimized for mouse instead of cat:
- PSTHBASELINEMEDIANX = 2 # PSTH median multiplier to estimate baseline
- PSTHPEAKMINTHRESH = 3 # PSTH peak detection threshold relative to estimated baseline, Hz
- STIMTYPES = 'mvi', 'grt'
- stimtype2axislabel = {'grt':'grating', 'mvi':'movie'}
- # exclude 'clr' kind, since 'nat' is a superset of 'clr' anyway
- #MVIKINDS = 'nat', 'pnk', 'shf'
- MVIKINDS = 'nat',
- ALLSTIMTYPES = MVIKINDS + ('grt',)
- STIMTYPESPLUSALL = list(STIMTYPES) + ['mvi+grt'] # strings
- STIMTYPESPLUSALLI = list(STIMTYPES) + [slice(None)] # for use as indices
- # iterate over run states in this order:
- ST8S = ['run', 'sit'] # excludes 'none'
- ALLST8S = ['none', 'run', 'sit']
- ST8COMBOS = [['none'], ['run', 'sit'], ['run', 'sit']]
- OPTOCOMBOS = [[False, True], [False], [True]]
- st82crit = odict((('run', '(run_speed > 1).mean() > 0.5'),
- ('sit', '(run_speed < 0.25).mean() > 0.5'),
- ('none', 'none')))
- rungreen = '#' + ''.join([ format(val, '02x') for val in [66, 168, 117] ])
- sitorange = '#' + ''.join([ format(val, '02x') for val in [235, 115, 9] ])
- st82clr = {'run':rungreen, 'sit':sitorange, 'none':'black'}
- st82alpha = {None:1, 'run':1, 'sit':0.5} # sitting is like suppression, so desaturate colour
- st82clrmap = {'run':plt.cm.Greens, 'sit':plt.cm.Oranges, 'none':plt.cm.Greys}
- # define example mseus and their colors:
- green = '#007f00' # deep green
- violet = '#9f3fff' # pure violet (7f00ff) is a little too dark
- darkblue = '#000064' # control condition
- optoblue = '#2a7fff' # light shade of blue for opto condition
- black = '#000000'
- asterisk = 5, 2, 0 # 5 sided, asterisk style, 0 deg rotation, has rounded corners though
- DEFSZ = 50 # default scatter plot point size (open points)
- exmpli2clr = {1:violet, 2:green, 3:optoblue} # example neurons 1, 2, 3
- # use different marker shapes to help distinguish example points:
- exmpli2mrk = {1:'x', 2:'+', 3:asterisk}
- exmpli2sz = {1:60, 2:70, 3:60} # marker size
- exmpli2lw = {1:2, 2:2, 3:2} # (marker) line width
- # use large closed circles for all example points:
- #exmpli2mrk = {1:'.', 2:'.', 3:'.'}
- #EXMPLSZ = 70 # example scatter plot point size (closed points)
- #exmpli2sz = {1:EXMPLSZ, 2:EXMPLSZ, 3:EXMPLSZ}
- fig1exmpli = 1 # fig1 uses example neuron 1
- fig5exmpli = 1 # fig5 uses example neuron 1
- mvimseu2exmpli = {'PVCre_2018_0003_s03_e03_u51':1, # example neuron 1
- #'PVCre_2018_0003_s03_e03_u52',
- 'PVCre_2017_0008_s12_e06_u56':2, # example neuron 2
- #'PVCre_2017_0008_s09_e07_u22',
- }
- grtmseu2exmpli = {'PVCre_2018_0003_s03_e02_u51':1, # example neuron 1
- #'PVCre_2018_0003_s03_e02_u52',
- #'PVCre_2018_0003_s02_e02_u11',
- #'PVCre_2018_0003_s02_e02_u38':
- #'PVCre_2018_0001_s02_e06_u64':
- #'PVCre_2018_0001_s02_e06_u61':
- #'PVCre_2018_0001_s02_e06_u64':
- #'PVCre_2018_0001_s05_e02_u16':3, # old example neuron 3
- 'PVCre_2018_0003_s05_e02_u164':3, # new example neuron 3 (2019-05-31)
- }
- mvimsu2exmpli = { mseustr2msustr(k):v for k, v in mvimseu2exmpli.items() }
- grtmsu2exmpli = { mseustr2msustr(k):v for k, v in grtmseu2exmpli.items() }
- msu2exmpli = {**mvimsu2exmpli, **grtmsu2exmpli} # merge the two
- burstclr = 'red'
- OPTOS = False, True
- opto2tuni = {False:0, True:1} # for indexing into tun_params field in Tuning table
- opto2fb = {None:'', False:'feedback', True:'suppression'}
- opto2clr = {None:darkblue, False:darkblue, True:optoblue}
- opto2alpha = {None:1, False:1, True:0.5} # with opto, feedback is removed, so desaturate colour
- modmeasures = ['meanrate', 'meanrate02', 'meanrate35', 'blankmeanrate', 'blankcondmeanrate',
- 'meanburstratio', 'blankmeanburstratio', 'blankcondmeanburstratio',
- 'spars', 'rel', 'meanpkw', 'snr']
- modmeasuresnoblankcond = modmeasures.copy()
- modmeasuresnoblankcond.remove('blankcondmeanrate')
- modmeasuresnoblankcond.remove('blankcondmeanburstratio')
- modmeasuresnoblank = modmeasuresnoblankcond.copy()
- modmeasuresnoblank.remove('blankmeanrate')
- modmeasuresnoblank.remove('blankmeanburstratio')
- measure2axislabel = {'meanrate':'FR', 'meanrate02':'FR', 'meanrate35':'FR',
- 'blankmeanrate':'FR', 'blankcondmeanrate':'FR',
- 'meanburstratio':'burst ratio', 'blankmeanburstratio':'burst ratio',
- 'blankcondmeanburstratio':'burst ratio',
- 'spars':'sparseness', 'rel':'reliability',
- 'meanpkw':'mean peak width', 'snr':'SNR'}
- measure2axisunits = {'meanrate':' (spk/s)', 'meanrate02':' (spk/s)', 'meanrate35':' (spk/s)',
- 'meanpkw':' (s)'}
- short2longaxislabel = {'FR':'Firing rate'}
- # restrict FiringPattern to our usual burst definition from Lu et al, 1992:
- BURSTCRITERION = ('(fp_dtsilent BETWEEN 0.099 and 0.101) AND '
- '(fp_dtburst BETWEEN 0.0039 AND 0.0041)')
- fmodes = ['', 'nb', 'b', 'nr'] # (narrow) firing modes
- fmode2txt = {'':'all', 'nb':'nonburst', 'b':'burst', 'nr':'nonrand'}
- # Tuning table keys used for gratings:
- ivs_order = 'grat_orientation, opto1'
- tun_model = 'sum_of_gaussians'
- model = 'threshlin' # linear or threshlin
- relaverage = 'mean' # reliability averaging method: mean or median
- DEFAULTFIGURESIZE = 2.8, 2.8 # normal
- DUMBBELLFIGSIZE = 3.3, 4.4 # dumbbell plots
- RASTERSCALEX = 0.75
- OFFSETS = -0.5, 0.75 # wide raster time offsets, sec
- OFFSETS02 = 0, -3 # limits to spikes from t=0 to 2 s, assuming 5 sec movies
- OFFSETS35 = 3, 0 # limits to spikes from t=3 to 5 s
- NTRIALSMVIGRT = 120 # limit movies to first 120 trials for more direct comparison with gratings
- PSTHHEIGHT = 3
- binw, tres, kernel = 0.02, 0.001, 'gauss' # for calculating PSTHs, sec
- psthplottres = 0.010 # sec
- ssx = intround(psthplottres / tres) # subsample factor to get the desired plot tres
- MINNTRIALSAMPLETHRESH = 10
- # ON/OFF/transient classification:
- PDT = 0.050 # propagation delay time for movie onset (s)
- WINDT = 0.20 # duration separating start of pre and PDT, and PDT and end of post (s)
- EXCLTR = 0, 2*PDT # time range around PDT to treat as transient response
- assert EXCLTR[0] <= PDT <= EXCLTR[1]
- ## Query the user:
- msg = ('Constants defined. What next?\n'
- '0: calculate from scratch, including fits\n'
- '1: calculate from scratch, but load fits from pickle\n'
- '2: load from pickles and exit\n'
- '3: load from pickles and assume connected\n'
- '>> ')
- response = int(input(msg))
- # load fits from pickle? save heaps of CPU time:
- LOADFITSFROMPICKLE = False if response == 0 else True
- if response in [2, 3]:
- subfolder = EXPTYPE if EXPTYPE != 'pvmvis' else ''
- name2val = load_all(subfolder=subfolder)
- locals().update(name2val)
- if response == 2:
- sys.exit()
- ## Collect series, experiments, units from database:
- '''
- # all single movie opto experiments in LGN that have been sorted:
- nats = ((e & 'e_name LIKE "MAS_%%"' & 'e_name NOT LIKE "%%+%%"' & 'e_optowl = 470')
- & (s & 's_depth > 2500' & 's_depth < 3500') # upper limit excludes TRN series
- & u)
- # all pink noise opto experiments in LGN that have been sorted:
- pnks = ((e & 'e_name LIKE "MAS_%%PNK%%"' & 'e_optowl = 470')
- & (s & 's_depth > 2500' & 's_depth < 3500')
- & u)
- '''
- """
- PVCre mice:
- * PVCre_2018_0002_s04 - good RFs, OK tuning, excellent opto
- ** PVCre_2019_0001_s05 - excellent RFs, OK tuning, good opto
- - make sure to exclude e01 and e02 due to seizure between e02 and e03
- ** PVCre_2019_0001_s06 - excellent RFs and tuning, OK opto, started by Davide
- - sorted: single unit RFs are quite bad though, excluding
- - this mouse's pupil pointed very high under visible light illumination, in the dark
- with dilated pupil, it pointed forward normally. Another good reason to exclude
- - PVCre_2019_0001_s07 - excellent RFs and tuning, very weak opto
- * PVCre_2019_0001_s08 - excellent RFs and tuning, OK opto. stained with DiI
- * PVCre_2019_0002_s05 - excellent RFs and tuning, good opto
- * PVCre_2019_0002_s06 - excellent RFs, OK tuning, excellent opto
- * PVCre_2019_0002_s07 - excellent RFs, good tuning and opto
- ** PVCre_2019_0002_s08 - excellent RFs, tuning and opto! stained with DiI
- - sorted, single unit RFs are good
- As of 2018-12-13, both PVCre_2018_0010_s02 and s05 are also excluded via updated
- 's_region' values. Both were intended to be in dLGN, but may have somehow missed it, based on
- weird RFs
- Ntsr1cre mice:
- * Ntsr1Cre_2019_0002_s03: sorted, mov rasters look wrong (bad sorting? asleep mouse?),
- 100 trials/cond, no running, but opto effects seem consistent
- - check itracking movies, resort this one
- - Ntsr1Cre_2019_0002_s04: no mov
- * Ntsr1Cre_2019_0002_s05: sorted but only 3 units, only 2 have mov PSTH peaks, ~ consistent
- - Ntsr1Cre_2019_0003: no mov
- - Ntsr1Cre_2019_0004 has big cortical hole:
- - Ntsr1Cre_2019_0004_s02: nice mov PSTHs, but no opto effect at all for mov or grt
- - Ntsr1Cre_2019_0004_s03: ??
- ** Ntsr1Cre_2019_0007_s04 - good
- - Ntsr1Cre_2019_0007_s05 - suspicious about if it's really LGN, SU results look poor
- ** Ntsr1Cre_2019_0008_s05 - good
- - do older mice just for their grating data, do shell/core and other separation
- """
- mvis = get_exps(locals(), EXPTYPE)
- mviseries = s & mvis
- mviunits = up & mvis
- mvispikes = spk & mvis
- if EXPTYPE == 'negntsrmvis': # force inclusion of a V1 movie recording
- mvieye = EyeIxtract() & [{'m':'Ntsr1Cre_2020_0001', 's':2, 'e':11}, # V1
- {'m':'Ntsr1Cre_2020_0001', 's':4, 'e':10}] # dLGN
- else:
- mvieye = EyeIxtract() & mvis
- mvirun = rt & mvis
- # fetch unique animals included in data set:
- #np.unique(mviseries.fetch('m'))
- ## TODO: need to check which units are plausibly in dLGN, not all units on the shank will be!
- ## - this is mostly done ahead of time during spike sorting, MUA envelope RF analysis suggests
- ## some single units have maxchans that are outside the chan range of significant RFs of
- ## the population, though there are some caveats to that - maybe maxchan is too stringent,
- ## or maybe it's possible for one or two single units to have visual responses even when
- ## the population doesn't show much response to sparse noise
- # all grating oritun opto experiments in LGN that have been sorted that occurred during the
- # same series as movie experiments:
- grts = ((e & 'e_name LIKE "%%oritun%%"' & 'e_optowl IN (465, 470)')
- & mviseries
- & u)
- grtseries = mviseries & grts
- # use the same units for gratings as in movies:
- #grtunits = mviunits & grts
- # use units active during a grating experiment, not necessarily the same units as in movies:
- grtunits = up & grts
- grtspikes = spk & grts
- grttun = tun & grts
- spons = ((e & 'e_name LIKE "%%spont%%"' & 'e_optowl IN (465, 470)')
- & mviseries
- & u)
- # get all possible movie mseus, regardless of rate:
- mvimseus = []
- for mse in mvis.fetch(dj.key):
- uids = (up & mse).fetch('u')
- for uid in uids:
- mseu = mse.copy()
- mseu['u'] = uid
- mvimseus.append(mseu)
- # get all possible grating mseus, regardless of rate:
- grtmseus = []
- for mse in grts.fetch(dj.key):
- uids = (up & mse).fetch('u')
- for uid in uids:
- mseu = mse.copy()
- mseu['u'] = uid
- grtmseus.append(mseu)
- ## TODO: grtmseu filtering by visual responsiveness should be done here, but would currently
- ## result in index errors in the big grtresp loop below, so it's done piecemeal in the loop:
- #grtmseus = (grts * up).fetch(dj.key)
- #grtmseus = (grts * up) & (vd.Units() & 'n_crit=1' & 'z_crit > 2.5')
- # get all possible spontaneous mseus, regardless of rate:
- sponmseus = []
- for mse in spons.fetch(dj.key):
- uids = (up & mse).fetch('u')
- for uid in uids:
- mseu = mse.copy()
- mseu['u'] = uid
- sponmseus.append(mseu)
- msefmt = '{m}_s{s:02}_e{e:02}'
- msufmt = '{m}_s{s:02}_u{u:02}'
- mseufmt = '{m}_s{s:02}_e{e:02}_u{u:02}'
- mvimseustrs = [ mseufmt.format(**mseu) for mseu in mvimseus ]
- grtmseustrs = [ mseufmt.format(**mseu) for mseu in grtmseus ]
- sponmseustrs = [ mseufmt.format(**mseu) for mseu in sponmseus ]
- mvimsustrs = [ msufmt.format(**unit) for unit in mviunits ]
- grtmsustrs = [ msufmt.format(**unit) for unit in grtunits ]
- mvigrtmsustrs = list(np.union1d(mvimsustrs, grtmsustrs)) # superset of the two lists
- if response == 3:
- sys.exit()
- ## Calculate:
- ## collect movie rasters, meanrates, burst info, PSTHs, sparseness, reliability and PSTH
- ## peak info by mseu, movie kind, run state, and opto condition:
- levels = [mvimseustrs, MVIKINDS, ALLST8S, OPTOS]
- names = ['mseu', 'kind', 'st8', 'opto']
- mi = pd.MultiIndex.from_product(levels, names=names)
- columns = ['dt', 'optotrange', 'trialis',
- 'raster', 'braster', 'nbraster', 'wraster', 'wbraster', 'wnbraster', # rasters
- 'rates', 'rate02s', 'rate35s', 'blankrates', # single trial FRs
- 'meanrate', 'meanrate02', 'meanrate35', 'blankmeanrate', # trial-averaged FRs
- 'burstis', 'wburstis', 'burstratios', 'blankburstratios', # single trial BRs
- 'meanburstratio', 'blankmeanburstratio', # trial-averaged BRs
- 'psth', 'bpsth', 'nbpsth', 'wpsth', 'wbpsth', 'wnbpsth', 't', 'wt', 'bins', 'wbins',
- 'spars', 'rhos', 'rel', 'snr', 'pkts', 'pkws', 'meanpkw']
- mviresp = pd.DataFrame(index=mi, columns=columns)
- ntotpeaks = 0 # total number of peaks detected
- for mvi in mvis:
- msestr = msefmt.format(**mvi)
- print(msestr)
- exp = mvi['e']
- moviecond = mov & mvi
- kind2stimis = moviecond.moviekind2stimis() # movie kind to stimi mapping
- trials = trial & mvi
- alltrialis = np.sort(trials.fetch('trial_id')) # trialis corresponding to rows in tranges
- #rates = ra & mvi
- uids = (up & mvi).fetch('u')
- # all trials, in trial ID order:
- rasters, tranges, _, opto, _ = (spk & mvi).get_rasters()
- wrasters, *_ = (spk & mvi).get_rasters(offsets=OFFSETS) # w = wide, same returned tranges
- rasters02, *_ = (spk & mvi).get_rasters(offsets=OFFSETS02) # from 0 to 2 s, same tranges
- rasters35, *_ = (spk & mvi).get_rasters(offsets=OFFSETS35) # from 3 to 5 s, same tranges
- wtranges = tranges.copy() + OFFSETS
- #tranges02 = tranges.copy() + OFFSETS02
- #tranges35 = tranges.copy() + OFFSETS35
- # make sure trial duration is very consistent:
- dts = tranges.ptp(axis=1)
- dt = dts.mean() # mean trial duration
- dtstd = dts.std()
- assert dtstd / dt < 0.005
- # make sure opto tranges are consistent relative to trial tranges:
- optotrialstimt0s = tranges[opto == True][:, 0]
- evtranges = (evtimes & mvi).fetch1('ev_tranges')
- alloptotranges = evtranges - np.reshape(optotrialstimt0s, (-1, 1))
- optotrange = alloptotranges.mean(axis=0)
- optotrangestd = alloptotranges.std(axis=0)
- assert (optotrangestd / optotrange < 0.005).all()
- # get tranges and rasters for just the period before stimulus onset during which opto
- # can be on, i.e. from optotrange[0] (-0.272 s) to 0 s before each stimulus onset:
- blankdt = abs(optotrange[0]) # averages to 0.283 s across all experiments
- print('mvi blankdt:', blankdt)
- blanktranges = np.zeros_like(tranges)
- blanktranges[:, 0] = tranges[:, 0] - blankdt
- blanktranges[:, 1] = tranges[:, 0]
- blankrasters, _, _, _, _ = (spk & mvi).get_rasters(tranges=blanktranges)
- mviopto = evcond & mvi & 'ev_chan="opto1"'
- mvioptostimis, mvioptovals = mviopto.fetch('stim_id', 'ec_val')
- # calculate burst ratios for all units:
- firingpatterns = fp & mvi & BURSTCRITERION
- brus = firingpatterns.burst_ratio(tranges=tranges)
- wbrus = firingpatterns.burst_ratio(tranges=wtranges)
- blankbrus = firingpatterns.burst_ratio(tranges=blanktranges)
- for kind, kindstimis in kind2stimis.items():
- print(kind)
- if not kindstimis or kind not in MVIKINDS: # exclude empty stimis and 'clr' kind
- continue
- kindtrials = (trials & 'stim_id IN %s' % seq2sql(kindstimis))
- kindtrialis = kindtrials.fetch('trial_id')
- trialiis = np.isin(alltrialis, kindtrialis) # bool with same shape as alltrialis
- # calculate narrow and wide PSTH bins for all units:
- tmin, tmax = 0, dt
- bins = split_tranges([(tmin, tmax)], binw, tres) # narrow
- midbins = bins.mean(axis=1) # narrow
- wbins = split_tranges([(tmin+OFFSETS[0], tmax+OFFSETS[1])], binw, tres) # wide
- wmidbins = wbins.mean(axis=1) # wide
- for st8, st8crit in st82crit.items():
- if st8crit == 'none': # 'none' is not in the State.Trial table
- st8trialis = alltrialis
- else:
- st8trialis = (st8t & mvi & {'st8_crit':st8crit}).fetch('trial_id')
- st8trialis = np.sort(st8trialis)
- for opto in OPTOS:
- optois = mvioptovals == opto
- # limit stimis by opto:
- optostimis = mvioptostimis[optois]
- # also limit stimis by movie kind:
- stimis = np.intersect1d(kindstimis, optostimis)
- if len(stimis) == 0:
- continue
- optotrials = trials & 'stim_id IN %s' % seq2sql(optostimis)
- optotrialis = optotrials.fetch('trial_id')
- # intersect kind, st8 and opto trialis:
- trialis = np.intersect1d(kindtrialis, st8trialis)
- trialis = np.intersect1d(trialis, optotrialis)
- trialiis = np.isin(alltrialis, trialis) # bool with same shape as alltrialis
- if len(trialis) == 0:
- print('%s %s %s has no trials, skipping' % (kind, st8, opto))
- continue
- for uid in uids:
- if uid not in rasters:
- continue
- mseu = {'m':mvi['m'], 's':mvi['s'], 'e':mvi['e'], 'u':uid}
- mseustr = mseufmt.format(**mseu)
- raster = rasters[uid][trialiis] # subset of raster trials
- wraster = wrasters[uid][trialiis] # subset of wraster trials
- blankraster = blankrasters[uid][trialiis] # subset of blankraster trials
- raster02 = rasters02[uid][trialiis][:NTRIALSMVIGRT]
- raster35 = rasters35[uid][trialiis][:NTRIALSMVIGRT]
- rates = np.array([ len(trialspikes)/dt for trialspikes in raster ])
- rate02s = np.full(rates.shape, np.nan) # init with full shape with NaN
- rate35s = np.full(rates.shape, np.nan)
- blankrates = np.full(rates.shape, np.nan)
- # fill only the first n trials/blanks:
- rate02s[:len(raster02)] = np.array([ len(trialspikes)/2 for trialspikes
- in raster02 ])
- rate35s[:len(raster35)] = np.array([ len(trialspikes)/2 for trialspikes
- in raster35 ])
- blankrates[:len(blankraster)] = np.array([ len(bspikes)/blankdt for bspikes
- in blankraster ])
- meanrate = len(np.hstack(raster)) / dt / len(raster)
- #assert np.isclose(meanrate, rates.mean()) # sanity check
- meanrate02 = len(np.hstack(raster02)) / 2 / len(raster02)
- meanrate35 = len(np.hstack(raster35)) / 2 / len(raster35)
- blankmeanrate = len(np.hstack(blankraster)) / blankdt / len(blankraster)
- if meanrate < RATETHRESH:
- print('%s %s %s %s does not meet RATETHRESH, skipping'
- % (mseustr, kind, st8, opto))
- continue
- mviresprow = mviresp.loc[mseustr, kind, st8, opto]
- # technically, this is the wrong approach, and could result in writing
- # to a copy (an independent Series), which isn't the intention (see https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy)
- # but it seems to work because the mviresp DataFrame is being filled one
- # row at a time, and wasn't initialized explicitly by setting everything
- # to e.g. np.nan in a single vector operation
- mviresprow['dt'] = dt
- mviresprow['optotrange'] = optotrange
- mviresprow['trialis'] = trialis
- mviresprow['raster'] = raster # narrow raster
- mviresprow['wraster'] = wraster # wide raster, including OFFSETS
- mviresprow['rates'] = rates # based on narrow raster
- mviresprow['rate02s'] = rate02s
- mviresprow['rate35s'] = rate35s
- mviresprow['blankrates'] = blankrates
- mviresprow['meanrate'] = meanrate # based on narrow raster
- mviresprow['meanrate02'] = meanrate02
- mviresprow['meanrate35'] = meanrate35
- mviresprow['blankmeanrate'] = blankmeanrate
- # fill in the corresponding trial burst ratios:
- for bru, wbru, blankbru in zip(brus, wbrus, blankbrus):
- if bru['u'] == uid:
- assert bru['m'] == mseu['m']
- assert bru['s'] == mseu['s']
- assert bru['e'] == mseu['e']
- burstis = np.array(bru['burst_index'], # narrow, ragged array
- dtype=object)[trialiis]
- wburstis = np.array(wbru['burst_index'], # wide
- dtype=object)[trialiis]
- #blankburstis = np.array(blankbru['burst_index'],
- # dtype=object)[trialiis]
- burstratios = np.array(bru['burst_ratio'])[trialiis] # narrow
- blankburstratios = np.array(blankbru['burst_ratio'])[trialiis]
- mviresprow['burstis'] = burstis # narrow
- mviresprow['wburstis'] = wburstis # wide, for plotting
- mviresprow['burstratios'] = burstratios # narrow
- mviresprow['blankburstratios'] = blankburstratios
- mviresprow['meanburstratio'] = np.nanmean(burstratios)
- mviresprow['blankmeanburstratio'] = np.nanmean(blankburstratios)
- break # out of bru loop
- # collect narrow burst and non-burst spike rasters:
- braster, nbraster = [], []
- for row, bis in zip(raster, burstis): # iterate over trials in raster
- if len(bis) == 0:
- bis = [] # none
- nbis = slice(None) # all
- else:
- allis = np.arange(len(row)) # all spikes in this trial
- nbis = np.setdiff1d(allis, bis)
- braster.append(row[bis])
- nbraster.append(row[nbis])
- # collect wide burst and non-burst spike rasters:
- wbraster, wnbraster = [], []
- for row, bis in zip(wraster, wburstis): # iterate over trials in wraster
- if len(bis) == 0:
- bis = [] # none
- nbis = slice(None) # all
- else:
- allis = np.arange(len(row)) # all spikes in this trial
- nbis = np.setdiff1d(allis, bis)
- wbraster.append(row[bis])
- wnbraster.append(row[nbis])
- braster = np.array(braster, dtype=object) # ragged array
- nbraster = np.array(nbraster, dtype=object)
- wbraster = np.array(wbraster, dtype=object)
- wnbraster = np.array(wnbraster, dtype=object)
- mviresprow['braster'] = braster # narrow burst raster
- mviresprow['nbraster'] = nbraster # narrow non-burst raster
- mviresprow['wbraster'] = wbraster # wide burst raster
- mviresprow['wnbraster'] = wnbraster # wide non-burst raster
- # calculate narrow PSTHs:
- psth = raster2psth(raster, bins, binw, tres, kernel)
- bpsth = raster2psth(braster, bins, binw, tres, kernel)
- nbpsth = raster2psth(nbraster, bins, binw, tres, kernel)
- # calculate wide PSTHs:
- wpsth = raster2psth(wraster, wbins, binw, tres, kernel)
- wbpsth = raster2psth(wbraster, wbins, binw, tres, kernel)
- wnbpsth = raster2psth(nbraster, wbins, binw, tres, kernel)
- mviresprow['psth'] = psth # narrow PSTH
- mviresprow['bpsth'] = bpsth # narrow burst PSTH
- mviresprow['nbpsth'] = nbpsth # narrow non-burst PSTH
- mviresprow['wpsth'] = wpsth # wide PSTH
- mviresprow['wbpsth'] = wbpsth # wide burst PSTH
- mviresprow['wnbpsth'] = wnbpsth # wide non-burst PSTH
- mviresprow['t'] = midbins # narrow PSTH timepoints
- mviresprow['wt'] = wmidbins # wide PSTH timepoints
- mviresprow['bins'] = bins # narrow PSTH bins
- mviresprow['wbins'] = wbins # wide PSTH bins
- # calculate PSTH sparseness on narrow PSTH:
- spars = sparseness(psth)
- mviresprow['spars'] = spars
- # collect single trial signals from narrow PSTH, calc response reliability:
- trialsignals = []
- for trialspiketimes in raster:
- trialsignal = raster2psth([trialspiketimes], bins, binw, tres, kernel)
- trialsignals.append(trialsignal)
- trialsignals = np.vstack(trialsignals)
- rel, rhos = reliability(trialsignals, average=relaverage)
- if not np.isnan(rel):
- assert rel > -0.05 # can come out weakly negative
- mviresprow['rhos'] = rhos
- mviresprow['rel'] = max(rel, 0) # clamp weakly negative values to 0
- mviresprow['snr'] = snr_baden2016(trialsignals)
- # detect peaks in narrow PSTH:
- baseline = PSTHBASELINEMEDIANX * np.median(psth)
- thresh = baseline + PSTHPEAKMINTHRESH # peak detection threshold
- spikets = np.sort(np.concatenate(raster)) # flatten raster spike times
- peakis, lis, ris = get_psth_peaks_gac(spikets, midbins, psth, thresh,
- verbose=VERBOSE)
- if VERBOSE:
- print() # new line
- npeaks = len(peakis)
- if npeaks > 0:
- mviresprow['pkts'] = peakis * tres # PSTH peak times (s)
- pkws = (ris - lis) * tres # PSTH peak widths (s)
- mviresprow['pkws'] = pkws
- mviresprow['meanpkw'] = pkws.mean() # mean PSTH peak width (s)
- ## TODO: use the peakis themselves for some other kind of measure,
- ## maybe plot distributions of peaks relative to stimulus start, or
- ## interpeak intervals, e.g.
- # use absence of PSTH peaks to decide which units/conditions to exclude?
- # this would obviously only work for movies, but that's what mviresp
- # is all about...
- #if npeaks == 0:
- # nrnids[state].append(nid) # save nonresponsive nids by state
- # continue # this PSTH has no peaks, skip all subsequent measures
- ntotpeaks += npeaks
- print('Detected %d movie PSTH peaks' % ntotpeaks)
- ## collect grating meanrates, rasters, tuning, and burst info by mseu,
- ## run state, and opto condition:
- levels = [grtmseustrs, ALLST8S, OPTOS]
- names = ['mseu', 'st8', 'opto']
- mi = pd.MultiIndex.from_product(levels, names=names)
- columns = ['dt', 'optotrange', 'trialis',
- 'raster', 'braster', 'nbraster', 'wraster', # rasters
- 'rates', 'blankrates', 'blankcondrates', # single trial FRs
- 'meanrate', 'blankmeanrate', 'blankcondmeanrate', # trial-averaged FRs
- 'burstis', 'wburstis', 'burstratios', # single trial BRs
- 'blankburstratios', 'blankcondburstratios', # single trial BRs
- 'meanburstratio', 'blankmeanburstratio', 'blankcondmeanburstratio', # trial-averaged BRs
- 'stimis', 'oris', 'tfreq', 'tun', 'tunrsq']
- grtresp = pd.DataFrame(index=mi, columns=columns)
- for grt in grts:
- msestr = msefmt.format(**grt)
- print(msestr)
- exp = grt['e']
- alltrials = trial & grt # all trials, including blanks
- # all stimis and trialis, in stimis order:
- allstimis, alltrialis = alltrials.fetch('stim_id', 'trial_id')
- # map triali to stimi for all trials:
- triali2stimi = { triali:stimi for stimi, triali in zip(allstimis, alltrialis) }
- alltrialis = np.sort(alltrialis) # trialis corresponding to rows in tranges
- trials = alltrials & (grat & 'grat_contrast != 0.0') # filter out blank trials
- blankcondtrials = alltrials & (grat & 'grat_contrast = 0.0') # blank trials only
- uids = (up & grt).fetch('u')
- # all trials, in trial ID (temporal) order:
- rasters, tranges, _, opto, _ = (spk & grt).get_rasters()
- wrasters, *_ = (spk & grt).get_rasters(offsets=OFFSETS) # w = wide, same returned tranges
- wtranges = tranges.copy() + OFFSETS
- # make sure opto tranges are consistent relative to trial tranges:
- optotrialstimt0s = tranges[opto == True][:, 0]
- evtranges = (evtimes & grt).fetch1('ev_tranges')
- alloptotranges = evtranges - np.reshape(optotrialstimt0s, (-1, 1))
- optotrange = alloptotranges.mean(axis=0)
- optotrangestd = alloptotranges.std(axis=0)
- assert (optotrangestd / optotrange < 0.005).all()
- # get tranges and rasters for just the period before stimulus onset during which opto
- # can be on, i.e. from optotrange[0] to 0 s before each stimulus onset:
- blankdt = abs(optotrange[0]) # averages to 0.285 s across all experiments, as in movies
- print('grt blankdt:', blankdt)
- blanktranges = np.zeros_like(tranges)
- blanktranges[:, 0] = tranges[:, 0] - blankdt
- blanktranges[:, 1] = tranges[:, 0]
- blankrasters, _, _, _, _ = (spk & grt).get_rasters(tranges=blanktranges)
- grtopto = evcond & grt & 'ev_chan="opto1"'
- grtoptostimis, grtoptovals = grtopto.fetch('stim_id', 'ec_val')
- # calculate burst ratios for all units:
- firingpatterns = fp & grt & BURSTCRITERION
- brus = firingpatterns.burst_ratio(tranges=tranges)
- wbrus = firingpatterns.burst_ratio(tranges=wtranges)
- blankbrus = firingpatterns.burst_ratio(tranges=blanktranges)
- # make sure trial duration is very consistent:
- dts = tranges.ptp(axis=1)
- dt = dts.mean() # mean trial duration, for all trials (non-blank and blank)
- dtstd = dts.std()
- assert dtstd / dt < 0.005
- # make sure opto tranges are consistent relative to trial tranges:
- optotrialstimt0s = tranges[opto == True][:, 0]
- evtranges = (evtimes & grt).fetch1('ev_tranges')
- alloptotranges = evtranges - np.reshape(optotrialstimt0s, (-1, 1))
- optotrange = alloptotranges.mean(axis=0)
- optotrangestd = alloptotranges.std(axis=0)
- assert (optotrangestd / optotrange < 0.005).all()
- # make sure temporal freq is very consistent:
- tfreqs = (grat & grt).fetch('grat_temp_freq')
- tfreq = tfreqs.mean()
- tfreqstd = tfreqs.std()
- assert tfreqstd / tfreq < 0.005
- for st8, st8crit in st82crit.items():
- if st8crit == 'none': # 'none' is not in the State.Trial table
- st8trialis = alltrialis
- else:
- st8trialis = (st8t & grt & {'st8_crit':st8crit}).fetch('trial_id')
- st8trialis = np.sort(st8trialis)
- for opto in OPTOS:
- optois = grtoptovals == opto
- # limit stimis by opto:
- optostimis = grtoptostimis[optois]
- optotrials = trials & 'stim_id IN %s' % seq2sql(optostimis)
- blankcondoptotrials = blankcondtrials & 'stim_id IN %s' % seq2sql(optostimis)
- optotrialis = optotrials.fetch('trial_id')
- blankcondoptotrialis = blankcondoptotrials.fetch('trial_id')
- # intersect st8 and opto trialis:
- trialis = np.intersect1d(st8trialis, optotrialis)
- blankcondtrialis = np.intersect1d(st8trialis, blankcondoptotrialis)
- if len(trialis) == 0:
- print('%s %s has no trials, skipping' % (st8, opto))
- continue
- stimis = np.array([ triali2stimi[triali] for triali in trialis ]) # triali order
- oristimis, oris = (grat & grt).fetch('stim_id', 'grat_orientation')
- oris = np.hstack([ oris[oristimis == stimi] for stimi in stimis ]) # triali order
- stimsortis = stimis.argsort(kind='stable')
- trialis = trialis[stimsortis] # sorted by stimi
- stimis = stimis[stimsortis] # sorted by stimi
- oris = oris[stimsortis] # sorted by stimi
- for uid in uids:
- if uid not in rasters:
- continue
- mseu = {'m':grt['m'], 's':grt['s'], 'e':grt['e'], 'u':uid}
- mseustr = mseufmt.format(**mseu)
- ## TODO: this should really be done once, where grtmseus are defined,
- ## but I don't want to break things right now (2019-07-15)
- if not (vd.Units() & mseu & ('n_crit=%d' % VDNCRIT)
- & ('z_crit > %f' % VDZCRIT)):
- # not a visually responsive mseu (this query ignores st8 and opto)
- print('%s is not visually responsive, skipping' % mseustr)
- continue
- raster = rasters[uid][trialis] # subset of trials, sorted by stimi
- wraster = wrasters[uid][trialis] # subset of trials, sorted by stimi
- blankraster = blankrasters[uid][trialis] # subset of trials, sorted by stimi
- blankcondraster = rasters[uid][blankcondtrialis] # subset of blank condition trials
- rates = np.array([ len(trialspikes)/dt for trialspikes in raster ])
- blankrates = np.array([ len(bspikes)/blankdt for bspikes in blankraster ])
- assert len(blankrates) == len(rates) # should be same shape as rates
- blankcondrates = np.array([ len(bcspikes)/dt for bcspikes in blankcondraster ])
- meanrate = len(np.hstack(raster)) / dt / len(raster)
- #assert np.isclose(meanrate, rates.mean()) # sanity check
- blankmeanrate = len(np.hstack(blankraster)) / blankdt / len(blankraster)
- if len(blankcondraster) == 0:
- blankcondmeanrate = np.nan
- else:
- blankcondmeanrate = len(np.hstack(blankcondraster)) / dt / len(blankcondraster)
- if meanrate < RATETHRESH:
- print('%s %s %s does not meet RATETHRESH, skipping'
- % (mseustr, st8, opto))
- continue
- grtresprow = grtresp.loc[mseustr, st8, opto]
- grtresprow['dt'] = dt
- grtresprow['optotrange'] = optotrange
- # overload 'trialis' name in df - includes non-blank trials and blank condition
- # trials, in that order:
- grtresprow['trialis'] = np.concatenate([trialis, blankcondtrialis])
- grtresprow['raster'] = raster # sorted by stimi
- grtresprow['wraster'] = wraster
- grtresprow['rates'] = rates
- grtresprow['blankrates'] = blankrates
- grtresprow['blankcondrates'] = blankcondrates
- grtresprow['meanrate'] = meanrate
- grtresprow['blankmeanrate'] = blankmeanrate
- grtresprow['blankcondmeanrate'] = blankcondmeanrate
- grtresprow['stimis'] = stimis # one entry per row in raster
- grtresprow['oris'] = oris # one entry per row in raster
- grtresprow['rates'] = rates
- grtresprow['tfreq'] = tfreq
- tunrow = (tun & mseu & {'ivs_order':ivs_order, 'tun_model':tun_model,
- 'st8_crit':st8crit})
- if tunrow: # not empty
- tunparams = tunrow.fetch1('tun_pars') # 2x5 array of tuning params
- tunrsq = tunrow.fetch1('tun_rsq') # len 2 array of R2 values
- tunspon = tunrow.fetch1('tun_spon_mean') # len 2 array of spon activity
- oti = opto2tuni[opto] # opto Tuning index
- params = list(tunparams[oti]) + [(tunspon[oti])]
- grtresprow['tun'] = params # len 6 array: 5 tuning params + spon
- """[dp, rp, rn, r0, sigma, spon]:
- oripref, gauss ampl @ oripref, gauss ampl @ antipref, gauss offset,
- gauss width, spon level
- so index 1 wrt 3 and 5 is a measure of ori tuning strength"""
- grtresprow['tunrsq'] = tunrsq[oti]
- else:
- print('WARNING: empty Tuning row for %s, %s' % (mseu, st8crit))
- # fill in the corresponding trial burst ratios:
- for bru, wbru, blankbru in zip(brus, wbrus, blankbrus):
- if bru['u'] == uid:
- assert bru['m'] == mseu['m']
- assert bru['s'] == mseu['s']
- assert bru['e'] == mseu['e']
- burstis = np.array(bru['burst_index'], # narrow, ragged array
- dtype=object)[trialis]
- wburstis = np.array(wbru['burst_index'], # wide
- dtype=object)[trialis]
- burstratios = np.array(bru['burst_ratio'])[trialis] # narrow
- blankburstratios = np.array(blankbru['burst_ratio'])[trialis]
- blankcondburstratios = np.array(bru['burst_ratio'])[blankcondtrialis]
- grtresprow['burstis'] = burstis # narrow
- grtresprow['wburstis'] = wburstis # wide, for plotting
- grtresprow['burstratios'] = burstratios # narrow
- grtresprow['blankburstratios'] = blankburstratios # narrow
- grtresprow['blankcondburstratios'] = blankcondburstratios # narrow
- # nanmean can raise "Mean of empty slice" warnings if argument
- # is empty or is all nan, but that's OK, returns a nan in both cases:
- grtresprow['meanburstratio'] = np.nanmean(burstratios)
- grtresprow['blankmeanburstratio'] = np.nanmean(blankburstratios)
- if len(blankcondburstratios) == 0:
- grtresprow['blankcondmeanburstratio'] = np.nan
- else:
- grtresprow['blankcondmeanburstratio'] = np.nanmean(blankcondburstratios)
- break # out of bru loop
- # collect narrow burst and non-burst spike rasters:
- braster, nbraster = [], []
- for row, bis in zip(raster, burstis): # iterate over trials in raster
- if len(bis) == 0:
- bis = [] # none
- nbis = slice(None) # all
- else:
- allis = np.arange(len(row)) # all spikes in this trial
- nbis = np.setdiff1d(allis, bis)
- braster.append(row[bis])
- nbraster.append(row[nbis])
- braster = np.array(braster, dtype=object) # ragged array
- nbraster = np.array(nbraster, dtype=object)
- grtresprow['braster'] = braster # narrow burst raster
- grtresprow['nbraster'] = nbraster # narrow non-burst raster
- ## Generate grating ori tuning dataframe for saving to a .pickle, and for convenient
- ## offline tuning plots:
- grttundf = grttun.fetch(format='frame')
- grttuninfo = (TuningInfo() & grttun).fetch(format='frame') # contains stimulus oris
- # convert separate m, s, e, u columns in grttundf to a single mseu column in a new temporary
- # df without a MultiIndex:
- grttunresprows = []
- for rowindex, rowvalues in grttundf.iterrows():
- mid, sid, eid, ivs_order, tun_model, st8_type, st8_crit, uid = rowindex
- mean, sem, spon_mean, pars, rsq = (rowvalues['tun_mean'],
- rowvalues['tun_sem'],
- rowvalues['tun_spon_mean'],
- rowvalues['tun_pars'],
- rowvalues['tun_rsq'])
- mseustr = key2mseustr({'m':mid, 's':sid, 'e':eid, 'u':uid})
- ori = grttuninfo.loc[mid, sid, eid, ivs_order]['ti_axes'][:, 0]
- newrow = {'mseu':mseustr, 'ivs_order':ivs_order, 'tun_model':tun_model,
- 'st8_type':st8_type, 'st8_crit':st8_crit,
- 'ori':ori, 'mean':mean, 'sem':sem, 'spon_mean':spon_mean,
- 'pars':pars, 'rsq':rsq}
- grttunresprows.append(newrow)
- newindex = ['mseu', 'ivs_order', 'tun_model', 'st8_type', 'st8_crit']
- grttunresp = pd.DataFrame(grttunresprows).set_index(newindex)
- ## collect spontaneous meanrates and meanburstratios by mseu and opto cond,
- ## ignore run st8 for now:
- levels = [sponmseustrs, OPTOS]
- names = ['mseu', 'opto']
- mi = pd.MultiIndex.from_product(levels, names=names)
- columns = ['meanrate', 'meanburstratio']
- sponresp = pd.DataFrame(index=mi, columns=columns)
- for spon in spons:
- msestr = msefmt.format(**spon)
- print(msestr)
- #exp = spon['e']
- optotranges = (evtimes & spon).fetch1('ev_tranges')
- # ensure that intervals between end and start of consecutive opto pulses are at least
- # as long as the pulse length:
- optodts = np.diff(optotranges, axis=1) # opto pulse durations (sec)
- maxoptodt = optodts.max() # sec
- isitranges = np.vstack([optotranges[:-1, 0], optotranges[1:, 1]]).T # shape: npulses-1, 2
- assert np.diff(isitranges, axis=1).min() > maxoptodt # sec
- # use one opto pulse duration before each pulse as the control trange for that pulse:
- ctrltranges = optotranges - optodts
- ctrldts = np.diff(ctrltranges, axis=1) # control trange durations (sec)
- spkspon = spk & spon
- fpspon = fp & spon & BURSTCRITERION
- optorasters, _, _, optovals, _ = spkspon.get_rasters(optotranges, offsets=[0, 0])
- ctrlrasters, _, _, ctrlvals, _ = spkspon.get_rasters(ctrltranges, offsets=[0, 0])
- optobrus = fpspon.burst_ratio(tranges=optotranges)
- ctrlbrus = fpspon.burst_ratio(tranges=ctrltranges)
- assert optovals.all() == True
- assert ctrlvals.all() == False
- assert list(optorasters) == list(ctrlrasters) # should have exact same set of units
- uids = list(optorasters)
- opto2rasters = {False:ctrlrasters, True:optorasters}
- opto2brus = {False:ctrlbrus, True:optobrus}
- opto2totaldt = {False:ctrldts.sum(), True:optodts.sum()}
- for opto in OPTOS:
- rasters = opto2rasters[opto]
- brus = opto2brus[opto]
- totaldt = opto2totaldt[opto] # sec
- for uid, bru in zip(uids, brus):
- trialspikets = np.concatenate(rasters[uid]) # concatenate spikes from all "trials"
- mseu = {'m':spon['m'], 's':spon['s'], 'e':spon['e'], 'u':uid}
- mseustr = mseufmt.format(**mseu)
- sponresprow = sponresp.loc[mseustr, opto]
- nspikes = len(trialspikets)
- meanrate = nspikes / totaldt
- sponresprow['meanrate'] = meanrate
- assert bru['u'] == uid # sanity check
- burstratios = bru['burst_ratio'] # one per "trial"
- meanburstratio = np.nanmean(burstratios)
- sponresprow['meanburstratio'] = meanburstratio
- ## find best movie and grating response parameters per unit, across experiments:
- mvilevels = [mvimsustrs, MVIKINDS, ALLST8S, OPTOS]
- grtlevels = [grtmsustrs, ALLST8S, OPTOS]
- mvinames = ['msu', 'kind', 'st8', 'opto']
- grtnames = ['msu', 'st8', 'opto']
- mvimi = pd.MultiIndex.from_product(mvilevels, names=mvinames)
- grtmi = pd.MultiIndex.from_product(grtlevels, names=grtnames)
- bestmvicolumns = modmeasuresnoblankcond
- bestgrtcolumns = ['meanrate', 'blankmeanrate', 'meanburstratio', 'blankmeanburstratio']
- bestmviresp = pd.DataFrame(index=mvimi, columns=bestmvicolumns)
- bestgrtresp = pd.DataFrame(index=grtmi, columns=bestgrtcolumns)
- stimtype2resp = {'mvi':mviresp, 'grt':grtresp}
- stimtype2bestresp = {'mvi':bestmviresp, 'grt':bestgrtresp}
- for stimtype in STIMTYPES:
- resp, bestresp = stimtype2resp[stimtype], stimtype2bestresp[stimtype]
- columns = bestresp.columns
- for index, row in resp.iterrows():
- mseu = index[0]
- msu = mseustr2msustr(mseu)
- bestindex = msu, *index[1:]
- for column in columns:
- newval = row[column]
- oldbestval = bestresp[column][bestindex]
- # the biggest value isn't always the best, e.g. lower meanpkw values are better:
- if column in ['meanpkw']:
- if pd.isna(oldbestval) or (newval < oldbestval):
- bestresp[column][bestindex] = newval # overwrite
- else:
- if pd.isna(oldbestval) or (newval > oldbestval):
- bestresp[column][bestindex] = newval # overwrite
- ## collect movie model fits by kind, state and opto:
- #levels = [mvimseustrs, MVIKINDS, ALLST8S]
- levels = [mvimseustrs, ['nat'], ['none']]
- names = ['mseu', 'kind', 'st8']
- mi = pd.MultiIndex.from_product(levels, names=names)
- fitcolumns = ['m', 'th', 'rsq', 'nbm', 'nbth', 'nbrsq', 'bm', 'bth', 'brsq']
- sampfitcolumns = ['ms', 'ths', 'rsqs', 'nbms', 'nbths', 'nbrsqs', 'bms', 'bths', 'brsqs',
- 'nrms', 'nrths', 'nrrsqs']
- fits = pd.DataFrame(index=mi, columns=fitcolumns)
- sampfits = pd.DataFrame(index=mi, columns=sampfitcolumns)
- nsamples = 1000
- np.random.seed(0) # fix random seed for identical results from np.random.choice() on each run
- if LOADFITSFROMPICKLE:
- print('Loading fits')
- fits = load('fits')
- sampfits = load('sampfits')
- nsamples = len(sampfits[sampfits['ms'].notnull()]['ms'].iloc[0]) # convoluted, but works
- mseustrfititerator = []
- else:
- print('Calculating fits')
- mseustrfititerator = mvimseustrs
- for mseustr in mseustrfititerator:
- for kind in ['nat']: #MVIKINDS:
- for st8 in ['none']: #ALLST8S:
- print(mseustr, kind, st8)
- psth = mviresp.loc[mseustr, kind, st8]['psth']
- if psth.isna().any(): # test both ctrl and opto
- continue
- # init sampfits lists:
- sampfits.loc[mseustr, kind, st8] = [], [], [], [], [], [], [], [], [], [], [], []
- # iterate over (narrow) firing modes (all, non-burst, burst, non-rand):
- for fmode in fmodes:
- # first do single fit of all relevant trials, fit and test data are
- # identical (no trial resampling):
- if fmode != 'nr': # don't bother filling nr column in the fits array
- psth = mviresp.loc[mseustr, kind, st8][fmode+'psth']
- ctrlpsth, optopsth = psth[False], psth[True]
- mm, b, rsq = fitmodel(ctrlpsth, optopsth, ctrlpsth, optopsth,
- model=model)
- th = -b / mm # threshold, i.e., x intercept
- fits.loc[mseustr, kind, st8][fmode+'m'] = mm
- fits.loc[mseustr, kind, st8][fmode+'th'] = th
- fits.loc[mseustr, kind, st8][fmode+'rsq'] = rsq
- # now sample trials multiple times, separate into fit and test data:
- if fmode != 'nr':
- ctrlraster = mviresp.loc[mseustr, kind, st8, False][fmode+'raster']
- optoraster = mviresp.loc[mseustr, kind, st8, True][fmode+'raster']
- else: # fmode == 'nr'
- # find number of burst spikes per trial and condition,
- # remove that same number randomly from each trial in each condition:
- ctrlraster, optoraster = [], []
- ctrlbraster = mviresp.loc[mseustr, kind, st8, False]['braster']
- optobraster = mviresp.loc[mseustr, kind, st8, True]['braster']
- # find number of burst spikes per trial in both conditions:
- bctrlcounts = [ len(ctrlbtrial) for ctrlbtrial in ctrlbraster ]
- boptocounts = [ len(optobtrial) for optobtrial in optobraster ]
- # build nr rasters from full rasters, trial by trial:
- ctrlfullraster = mviresp.loc[mseustr, kind, st8, False]['raster']
- optofullraster = mviresp.loc[mseustr, kind, st8, True]['raster']
- for bctrlcount, ctrltrial in zip(bctrlcounts, ctrlfullraster):
- nspikes = len(ctrltrial) - bctrlcount # num spikes to sample
- resampctrltrial = choice(ctrltrial, size=nspikes, replace=False)
- ctrlraster.append(resampctrltrial)
- for boptocount, optotrial in zip(boptocounts, optofullraster):
- nspikes = len(optotrial) - boptocount # num spikes to sample
- resampoptotrial = choice(optotrial, size=nspikes, replace=False)
- optoraster.append(resampoptotrial)
- ctrlraster, optoraster = np.array(ctrlraster), np.array(optoraster)
- nctrltrials = len(ctrlraster) # typically 200 trials when st8 == 'none'
- noptotrials = len(optoraster) # typically 200 trials when st8 == 'none'
- for ntrials in [nctrltrials, noptotrials]:
- if ntrials < MINNTRIALSAMPLETHRESH:
- print('Not enough trials to sample:', mseustr, kind, st8, fmode)
- continue # don't bother sampling, leave entry in sampfits empty
- ctrltrialis = np.arange(nctrltrials)
- optotrialis = np.arange(noptotrials)
- ctrlsamplesize = intround(nctrltrials / 2) # half for fitting, half for testing
- optosamplesize = intround(noptotrials / 2) # half for fitting, half for testing
- # probably doesn't matter if use ctrl or opto bins, should be the same:
- bins = mviresp.loc[mseustr, kind, st8, False]['bins']
- for samplei in range(nsamples):
- # randomly sample half the ctrl and opto trials, without replacement:
- ctrlfitis = np.sort(choice(ctrltrialis, size=ctrlsamplesize, replace=False))
- optofitis = np.sort(choice(optotrialis, size=optosamplesize, replace=False))
- ctrltestis = np.setdiff1d(ctrltrialis, ctrlfitis) # get the complement
- optotestis = np.setdiff1d(optotrialis, optofitis) # get the complement
- ctrlfitraster, ctrltestraster = ctrlraster[ctrlfitis], ctrlraster[ctrltestis]
- optofitraster, optotestraster = optoraster[optofitis], optoraster[optotestis]
- # calculate fit and test opto PSTHs, subsampled in time:
- ctrlfitpsth = raster2psth(ctrlfitraster, bins, binw, tres, kernel)[::ssx]
- optofitpsth = raster2psth(optofitraster, bins, binw, tres, kernel)[::ssx]
- ctrltestpsth = raster2psth(ctrltestraster, bins, binw, tres, kernel)[::ssx]
- optotestpsth = raster2psth(optotestraster, bins, binw, tres, kernel)[::ssx]
- #if np.isnan(ctrlfitpsth).any() or np.isnan(optofitpsth).any():
- # continue
- #assert not (np.isnan(ctrltestpsth).any() or np.isnan(optotestpsth).any())
- mm, b, rsq = fitmodel(ctrlfitpsth, optofitpsth, ctrltestpsth, optotestpsth,
- model=model)
- #if np.isnan(rsq) or rsq < RSQTHRESH: # poor fit
- # continue # skip this sample
- th = -b / mm # threshold, i.e., x intercept
- sampfits.loc[mseustr, kind, st8][fmode+'ms'].append(mm)
- sampfits.loc[mseustr, kind, st8][fmode+'ths'].append(th)
- sampfits.loc[mseustr, kind, st8][fmode+'rsqs'].append(rsq)
- ## FMI (feedback modulation index: (feedback - suppression) / (feedback + suppression)) and
- ## RMI (run modulation index: (run - sit) / (run + sit)) for each mseu,
- ## for nat mvi and grt experiments, for all measures.
- # set up mviFMI DataFrame:
- levels = [mvimseustrs, ALLST8S]
- names = ['mseu', 'st8']
- mi = pd.MultiIndex.from_product(levels, names=names)
- mviFMI = pd.DataFrame(index=mi, columns=modmeasures)
- # set up grtFMI DataFrame:
- levels = [grtmseustrs, ALLST8S]
- names = ['mseu', 'st8']
- mi = pd.MultiIndex.from_product(levels, names=names)
- grtFMI = pd.DataFrame(index=mi, columns=modmeasures) # lots of unused columns
- # set up mviRMI DataFrame:
- levels = [mvimseustrs, OPTOS]
- names = ['mseu', 'opto']
- mi = pd.MultiIndex.from_product(levels, names=names)
- mviRMI = pd.DataFrame(index=mi, columns=modmeasures)
- # set up grtRMI DataFrame:
- levels = [grtmseustrs, OPTOS]
- names = ['mseu', 'opto']
- mi = pd.MultiIndex.from_product(levels, names=names)
- grtRMI = pd.DataFrame(index=mi, columns=modmeasures)
- stimtype2FMI = {'mvi': mviFMI, 'grt': grtFMI}
- stimtype2RMI = {'mvi': mviRMI, 'grt': grtRMI}
- for stimtype in STIMTYPES:
- if stimtype == 'mvi':
- resp = mviresp.xs('nat', level='kind') # keep only nat movie measures
- else:
- resp = grtresp
- fmidf = stimtype2FMI[stimtype]
- rmidf = stimtype2RMI[stimtype]
- for st8 in ALLST8S:
- for measure in modmeasures:
- if measure not in resp:
- continue
- measurevals = resp[measure][(slice(None), st8)] # mseustr, opto
- feedback = measurevals[(slice(None), False)] # mseustr
- suppress = measurevals[(slice(None), True)] # mseustr
- assert (feedback.dropna() >= 0).all() # negative values violate -1 <= FMI <= 1
- assert (suppress.dropna() >= 0).all()
- # filter out cases where measureval in both conditions is 0:
- keepis = (feedback != 0) | (suppress != 0)
- if stimtype == 'mvi': # also filter by SNRTHRESH:
- maxsnrs = resp['snr'][(slice(None), st8)].max(level='mseu') # max SNR per mseu
- keepis = keepis & (maxsnrs >= SNRTHRESH)
- feedback, suppress = feedback[keepis], suppress[keepis]
- fmis = (feedback - suppress) / (feedback + suppress) # indexed by mseustr
- for mseustr, fmi in zip(fmis.index.values, fmis):
- fmidf.loc[mseustr, st8][measure] = fmi # fill in the appropriate rows
- for opto in OPTOS:
- for measure in modmeasures:
- if measure not in resp:
- continue
- measurevals = resp[measure][(slice(None), slice(None), opto)] # mseustr, st8
- runvals = measurevals[(slice(None), 'run')] # mseustr
- sitvals = measurevals[(slice(None), 'sit')] # mseustr
- assert (runvals.dropna() >= 0).all() # negative values violate -1 <= RMI <= 1
- assert (sitvals.dropna() >= 0).all()
- # filter out cases where measureval in both conditions is 0:
- keepis = (runvals != 0) | (sitvals != 0)
- if stimtype == 'mvi': # also filter by SNRTHRESH:
- maxsnrs = resp['snr'][(slice(None), slice(None), opto)].max(level='mseu')
- keepis = keepis & (maxsnrs >= SNRTHRESH)
- runvals, sitvals = runvals[keepis], sitvals[keepis]
- rmis = (runvals - sitvals) / (runvals + sitvals) # indexed by mseustr
- for mseustr, rmi in zip(rmis.index.values, rmis):
- rmidf.loc[mseustr, opto][measure] = rmi # fill in the appropriate rows
- ## maxFMI and maxRMI for each unit (msu), for nat mvi and grt experiments, for all measures.
- ## For each stimtype+st8+measure+unit combination, take the abs max across experiments. Also
- ## store the argmax (i.e., the mseu). Among other things, this is for comparison
- ## of movie and grating responses.
- # set up maxFMI DataFrame:
- levels = [mvigrtmsustrs, ALLST8S, STIMTYPES]
- names = ['msu', 'st8', 'stimtype']
- mi = pd.MultiIndex.from_product(levels, names=names)
- columns = ['mseu'] + modmeasures
- maxFMI = pd.DataFrame(index=mi, columns=columns)
- # set up maxRMI DataFrame:
- levels = [mvigrtmsustrs, OPTOS, STIMTYPES]
- names = ['msu', 'opto', 'stimtype']
- mi = pd.MultiIndex.from_product(levels, names=names)
- columns = ['mseu'] + modmeasures
- maxRMI = pd.DataFrame(index=mi, columns=columns)
- for stimtype in STIMTYPES:
- fmidf = stimtype2FMI[stimtype]
- rmidf = stimtype2RMI[stimtype]
- for st8 in ALLST8S:
- for measure in modmeasures:
- mseustrs = fmidf.index.levels[0].values
- for mseustr in mseustrs:
- fmi = fmidf.loc[mseustr, st8][measure]
- if pd.isna(fmi):
- continue
- msustr = mseustr2msustr(mseustr)
- oldfmi = maxFMI.loc[msustr, st8, stimtype][measure] # check existing val
- if pd.isna(oldfmi) or (abs(fmi) > abs(oldfmi)):
- maxFMI.loc[msustr, st8, stimtype][measure] = fmi # overwrite
- maxFMI.loc[msustr, st8, stimtype]['mseu'] = mseustr # overwrite
- for opto in OPTOS:
- for measure in modmeasures:
- mseustrs = rmidf.index.levels[0].values
- for mseustr in mseustrs:
- rmi = rmidf.loc[mseustr, opto][measure]
- if pd.isna(rmi):
- continue
- msustr = mseustr2msustr(mseustr)
- oldrmi = maxRMI.loc[msustr, opto, stimtype][measure] # check existing val
- if pd.isna(oldrmi) or (abs(rmi) > abs(oldrmi)):
- maxRMI.loc[msustr, opto, stimtype][measure] = rmi # overwrite
- maxRMI.loc[msustr, opto, stimtype]['mseu'] = mseustr # overwrite
- ## calculate eye position stdev and pupil area during all 3 stimtypes,
- ## as a function of run state and opto conditions, store in ipos_st8 and ipos_opto DataFrames:
- # Requires huxley file system to be mounted. Author: Davide Crombie
- # parameters for cleaning up the pupil displacement data:
- min0t = 1 # only interpolate over stretches missing data shorter than min0t (s)
- min1t = 1 # only keep stretches of continuous data longer than min1t (s)
- min1prop = .8 # only keep stretches of data with at least 80% non-NaN
- st8_rows = [] # to be filled with dictionaries containing relevant data
- opto_rows = []
- print('Calculating ipos_st8 and ipos_opto DataFrames...')
- for stimtype, stims in zip(['mvi', 'grt', 'spon'], [mvis, grts, spons]):
- print('...for stimtype %r' % stimtype)
- stimeye = EyeIxtract() & stims
- for key in stimeye.get_keys(): # these are mse keys
- if stimtype == 'mvi':
- moviecond = mov & key
- kind2stimis = moviecond.moviekind2stimis() # movie kind to stimi mapping
- # stim_id values considered in this analysis:
- kindstimis = tuple(kind2stimis['nat'])
- if len(kindstimis) == 0:
- continue # wrong kind of movie, e.g. PNK or SHF
- msestr = "%s_s%02d_e%02d" % (key['m'], key['s'], key['e'])
- print(msestr)
- # fetch and clean the eye position data:
- irow = EyeIxtract() & key
- ts = irow.fetch1('i_frame_times')
- eyedt = np.diff(ts).mean() # eye sampling period
- # x & y pupil positions relative to baseline (median) in deg. visual angle,
- # after subtraction of corneal reflection position:
- xpos, ypos = irow.fetch_position(cr_sub=True, deg=True) # in deg
- area = irow.fetch_area(mm=False, normalize=True) # in pix
- # position timecourses after elimintation of unclean stretches, smoothing,
- # and interpolation of remaining NaN values within clean stretches:
- min0len, min1len = min0t/eyedt, min1t/eyedt # convert from s to n timepoints
- xpos_clean = i_clean(xpos, min0len=min0len, min1len=min1len, min1prop=min1prop)
- ypos_clean = i_clean(ypos, min0len=min0len, min1len=min1len, min1prop=min1prop)
- area_clean = i_clean(area, min0len=min0len, min1len=min1len, min1prop=min1prop)
- conditions = {}
- # get keys for trial table restriction based on st8:
- conditions['run'] = (st8t & key & {'st8_crit':st82crit['run']}).get_keys()
- conditions['sit'] = (st8t & key & {'st8_crit':st82crit['sit']}).get_keys()
- if len(conditions['run']) == 0: # perhaps run data wasn't acquired for this exp
- assert len(conditions['sit']) == 0 # sanity check
- del conditions['run'] # del empty list
- del conditions['sit'] # del empty list
- # get keys for trial table restriction based on opto condition:
- if stimtype == 'mvi':
- if len(kindstimis) > 1: # only include opto conditions correspond to nat movie
- conditions[False] = (evcond & key & 'ev_chan="opto1"' & 'ec_val=0' &
- 'stim_id IN {}'.format(kindstimis)).get_keys()
- conditions[True] = (evcond & key & 'ev_chan="opto1"' & 'ec_val=1' &
- 'stim_id IN {}'.format(kindstimis)).get_keys()
- else:
- conditions[False] = (evcond & key & 'ev_chan="opto1"' & 'ec_val=0').get_keys()
- conditions[True] = (evcond & key & 'ev_chan="opto1"' & 'ec_val=1').get_keys()
- # compute eye position variance for all conditions:
- for condition in conditions: # iterate over all existing keys in conditions
- row = {}
- row['mse'] = msestr
- row['stimtype'] = stimtype
- # start and stop times for trials of the current state only:
- ## TODO: spons don't have trials, just periods of opto and non-opto
- trialis, t0s, t1s = (trial & conditions[condition]).fetch('trial_id',
- 'trial_on_time',
- 'trial_off_time')
- if len(t0s) == 0 or len(t1s) == 0:
- continue
- # sort pupil displacement data into trial matrices:
- xpos_trialmat, _ = get_trialmat(xpos_clean, ts, t0s, t1s)
- ypos_trialmat, _ = get_trialmat(ypos_clean, ts, t0s, t1s)
- area_trialmat, area_trialts = get_trialmat(area_clean, ts, t0s, t1s)
- area_trialmean = np.nanmean(area_trialmat, axis=1) # mean pupil area of each trial
- # save trial IDs:
- row['trialis'] = trialis
- # save mean pupil area for each trial:
- row['area_trialmean'] = area_trialmean
- # save pupil area signal for each trial:
- row['area_trialmat'] = area_trialmat
- row['area_trialts'] = area_trialts
- # save mean pupil area across trials:
- row['mean_area'] = np.nanmean(area_trialmean)
- # save eye position variability across trials:
- row['std_xpos_cross'] = np.nanstd(xpos_trialmat, axis=0).mean()
- row['std_ypos_cross'] = np.nanstd(ypos_trialmat, axis=0).mean()
- # save eye position variability within trials:
- row['std_xpos_within'] = np.nanmean(np.nanstd(xpos_trialmat, axis=1))
- row['std_ypos_within'] = np.nanmean(np.nanstd(ypos_trialmat, axis=1))
- # add dictionary to appropriate list:
- if condition in ['run', 'sit']:
- row['st8'] = condition
- st8_rows.append(row)
- elif condition in [False, True]:
- row['opto'] = condition
- opto_rows.append(row)
- # convert list of dicts to DataFrame:
- ipos_st8 = pd.DataFrame(st8_rows)
- if len(st8_rows) > 0:
- ipos_st8.set_index(['mse', 'stimtype', 'st8'], inplace=True) # convert to MultiIndex
- ipos_opto = pd.DataFrame(opto_rows)
- if len(opto_rows) > 0:
- ipos_opto.set_index(['mse', 'stimtype', 'opto'], inplace=True)
- ## collect unit spatial position, store in upos DataFrame:
- rows = [ {'msu':msu2msustr(msu), 'x':x, 'y':y}
- for (msu, x, y) in zip(*mviunits.fetch(dj.key, 'u_xpos', 'u_ypos')) ]
- upos = pd.DataFrame(rows)
- upos = upos.set_index(['msu'])
- ## collect runspeed traces, store in runspeed DataFrame:
- runspeedrows = [] # to be filled with dictionaries containing relevant data
- for mvi in mvis.fetch(dj.key):
- msestr = "%s_s%02d_e%02d" % (mvi['m'], mvi['s'], mvi['e'])
- print(msestr)
- # fetch run data for all trials in this mse, each is an array, sometimes some trials
- # are missing from the run data:
- runtrialis, runt, runspeed = (rt & mvi).fetch('trial_id', 'run_t', 'run_speed')
- if len(runspeed) == 0:
- continue # no run data for this mse
- moviecond = mov & mvi
- kind2stimis = moviecond.moviekind2stimis() # movie kind to stimi mapping
- #trials = trial & mvi
- #alltrialis = np.sort(trials.fetch('trial_id')) # trialis corresponding to rows in tranges
- for kind in ['nat']: #MVIKINDS:
- kindstimis = kind2stimis[kind]
- if len(kindstimis) == 0:
- continue # wrong kind of movie
- # get stimis from evcond table for this mvi and kind, and restrict to opto1 chan:
- mviopto = evcond & mvi & 'stim_id IN %s' % seq2sql(kindstimis) & {'ev_chan':'opto1'}
- mvioptostimis, mvioptovals = mviopto.fetch('stim_id', 'ec_val')
- assert len(mvioptostimis) == 2
- assert len(mvioptovals) == 2
- mvioptovals = np.bool_(mvioptovals) # convert from [0, 1] to [False, True]
- opto2stimi = {False:int(mvioptostimis[mvioptovals == False]),
- True :int(mvioptostimis[mvioptovals == True])}
- for opto in OPTOS:
- row = {}
- row['mse'] = msestr
- row['kind'] = kind
- row['opto'] = opto
- # extract run trials for the stimi corresponding to this kind and opto:
- stimi = opto2stimi[opto]
- stimtrialis = np.sort((trial & mvi & {'stim_id':stimi}).fetch('trial_id'))
- # make sure that stimtrialis is a subset of runtrialis:
- assert issubset(stimtrialis, runtrialis)
- # now it's safe to ask where stimtrialis are found in runtrialis:
- trialiis = runtrialis.searchsorted(stimtrialis)
- #print(opto, len(stimtrialis), len(runtrialis))
- row['trialis'] = runtrialis[trialiis]
- row['t'] = runt[trialiis]
- row['speed'] = runspeed[trialiis]
- runspeedrows.append(row)
- runspeed = pd.DataFrame(runspeedrows)
- if len(rows) > 0:
- runspeed = runspeed.set_index(['mse', 'kind', 'opto']) # convert to MultiIndex
- ## classify units in various ways, store classification in celltype DataFrame:
- # use normal Index, not MultiIndex, since there's only a single index column, otherwise
- # assigning values with `df.loc[indexval][colname] = foo` won't work
- index = pd.Index(mvigrtmsustrs, name='msu')
- columns = ['sbczscore', 'sbc', 'dsi', 'depth', 'normdepth', 'onoff', 'trans',
- 'mvionsetpsth', 't']
- celltype = pd.DataFrame(index=index, columns=columns)
- print('Classifying SbC')
- # NOTE: seems that the zscores stored in the Condition.ZScores table are only
- # for control conditions, i.e. opto == False, which is good
- SBCZSCORETHRESH = -3 # 0.0, -1.96, -2.58, -3
- # classify by SbC, keep the most negative sbczscore of each msu across grating experiments:
- msustr2maxzsmseu = {} # map msustr to most negative sbczscore mseu key
- for grt in grts:
- zscores = zs & grt
- uids = (up & grt).fetch('u')
- for uid in uids:
- msu = {'m':grt['m'], 's':grt['s'], 'u':uid}
- mseu = {'m':grt['m'], 's':grt['s'], 'e':grt['e'], 'u':uid}
- msustr = msufmt.format(**msu)
- mseustr = mseufmt.format(**mseu)
- unitzscore = zscores & {'u':uid}
- if len(unitzscore) == 0:
- print('No Condition.Zscores entry for', mseustr)
- continue
- sbckeys = unitzscore.SbC()
- assert len(sbckeys) == 1
- sbczscore = sbckeys[0]['sbc_zscore'] # median, across ori conditions
- if np.isnan(sbczscore):
- print('WARNING: SbC() returned NaN for', mseustr)
- meanrate = grtresp.loc[mseustr, 'none', False]['meanrate']
- if np.isnan(meanrate):
- continue # skip this non visually responsive mseu, based on z_crit and n_crit
- oldsbczscore = celltype.loc[msustr]['sbczscore']
- if np.isnan(oldsbczscore) or sbczscore < oldsbczscore: # keep most -ve sbczscore
- celltype.loc[msustr]['sbczscore'] = sbczscore # update
- msustr2maxzsmseu[msustr] = mseu
- # classify units as SbC == [False, True] according to SBCZSCORETHRESH:
- for msustr in mvigrtmsustrs:
- sbczscore = celltype.loc[msustr]['sbczscore']
- if np.isnan(sbczscore):
- sbc = np.nan # this is the default anyway
- elif sbczscore <= SBCZSCORETHRESH:
- sbc = True
- else:
- sbc = False
- celltype.loc[msustr]['sbc'] = sbc
- # plot rasters, for testing SBCZSCORETHRESH, and zscores in general:
- try:
- mseu = msustr2maxzsmseu[msustr]
- except KeyError: # not a single exp with an sbczscore found for this msu
- print('Did not collect an sbczscore for', msustr)
- continue
- '''
- # plot the non-blank, opto == False trials of this msu:
- offsets = [-1, 1]
- rasters, tranges, allstimis, opto, ukeys = (spk & mseu).get_rasters(offsets=offsets)
- assert len(rasters) == 1
- uid, raster = rasters.popitem()
- stimis = allstimis[opto == False]
- raster = raster[opto == False]
- blankstimis = (grat & mseu & 'stim_id IN %s' % seq2sql(stimis)
- & 'grat_contrast = 0.0').fetch('stim_id')
- assert len(blankstimis) == 1
- blankstimi = blankstimis[0] # should be 24
- raster = raster[stimis != blankstimi] # exclude blank trials
- title = "%s, zscore=%.3f, SbC=%s" % (msustr, sbczscore, sbc)
- ax = simpletraster(raster, dt=5, offsets=offsets, scatter=True, title=title)
- ax.set_title(title)
- '''
- # classify by direction selectivity, keep biggest DSI of each msu across grating experiments:
- print('Classifying DSI')
- resp = grtresp.xs('none', level='st8') # keep only none run state
- resp = resp.xs(False, level='opto') # keep only control opto state
- for mseustr, row in resp.iterrows():
- print(mseustr)
- tunparams = row['tun']
- if np.isnan(tunparams).all():
- continue # skip this mseustr
- dp, rp, rn, r0, sigma, spon = tunparams
- dsi = niell_DSI(rp, rn, r0)
- msustr = mseustr2msustr(mseustr)
- celltyperow = celltype.loc[msustr]
- olddsi = celltyperow['dsi']
- if np.isnan(olddsi) or dsi > olddsi:
- celltyperow['dsi'] = dsi # update
- # estimate dLGN depth of each unit using sparse noise MUA RFs to help delineate
- # at what channel dLGN starts and ends:
- print('Estimating dLGN depth')
- # first plot MUA RFs for each sparse noise experiment in each series, visually inspect them,
- # and manually enter estimated start and end channel of dLGN for each series into
- # msstr2dLGNchanrange:
- for series in mviseries:
- snkeys = (Series.Experiment() & series & 'e_name LIKE "%%Noise%%"').fetch(dj.key)
- for snkey in snkeys:
- print('Sparse noise exp:', snkey)
- rfs, axes, chans = rfs_for_cdata(snkey)
- # plot MUA RFs for inspection
- simple_plot_rfs(snkey, rfs, axes, chans, interpolation=None, nrows=4)
- # Now save all plots to disk in MUA_RFs folder:
- saveall(format='pdf')
- saveall(format='png')
- # Next, manually map MOUSE_SERIES string to tuple of uppermost and lowermost clear dLGN
- # channel (end inclusive) based on qualitative estimation of start and end of MUA RF
- # progression of each series, potentially considering multiple sparse noise experiments
- # per series:
- ## TODO: switch to using LGNChans table and its ms2maxchanrange() method instead of
- ## duplicating it here:
- msstr2dLGNchanrange = {'PVCre_2017_0006_s03': (29, 8),
- 'PVCre_2017_0008_s09': (25, 7),
- 'PVCre_2017_0008_s12': (24, 2),
- 'PVCre_2017_0015_s03': (28, 5),
- 'PVCre_2017_0015_s07': (19, 1),
- 'PVCre_2018_0001_s02': (11, 1),
- 'PVCre_2018_0001_s05': (25, 6),
- 'PVCre_2018_0003_s02': (20, 3),
- 'PVCre_2018_0003_s03': (14, 2),
- 'PVCre_2018_0003_s05': (20, 9),
- 'PVCre_2019_0002_s08': (21, 1),
- 'Ntsr1Cre_2020_0001_s04': (31, 5),
- }
- # get dLGN depth of each unit in each series, add to celltype df:
- for series in mviseries:
- msstr = ms2msstr(series)
- hichan, lochan = msstr2dLGNchanrange[msstr] # hi and lo chans (in terms of depth)
- # could use coords entered into ProbeModel table, but easier to use SiteLoc dict
- # in probes.expio:
- #probemodel, chans, coords = (pm & (pr & series)).fetch1('pm', 'pm_chans', 'pm_coords')
- probemodel = (pr & series).fetch1('pm') # string
- probe = expio.probes.getprobe(probemodel)
- # y coord of upper and lowermost dLGN chan, use upper as reference and upper-lower as depth
- # normalization for all units in this series:
- y0 = probe.SiteLoc[hichan][1]
- y1 = probe.SiteLoc[lochan][1]
- uids, ys = (up & series).fetch('u', 'u_ypos')
- depths = ys - y0
- dLGNspan = y1 - y0
- normdepths = depths / dLGNspan
- for uid, depth, normdepth in zip(uids, depths, normdepths):
- msustr = msstr + '_u%02d' % uid
- if depth < 0:
- print('WARNING: %s has depth %.1f < 0, skipping entry in celltype df'
- % (msustr, depth))
- continue
- if normdepth > 1.5:
- print('WARNING: %s has normdepth %.1f > 1.5, skipping entry in celltype df'
- % (msustr, normdepth))
- continue
- celltype.loc[msustr, 'depth'] = depth
- celltype.loc[msustr, 'normdepth'] = normdepth
- '''
- # classify by layer (shell/core), pretty much deprecated, don't trust the code in rf.py
- # that does all of this...
- # NOTE: several (13/80 as of 2020-02-11) msu end up with no layer classification because
- # they fall outside of the detected range of chans with MUA envelope RFs of sufficient score
- print('Classifying shell/core')
- print('NOTE: something is wrong with this code, hangs on first series at full CPU')
- for series in mviseries:
- try:
- unitlayers = rf.get_layer(series, snexpid=-1) # list of dicts
- except RuntimeError as err:
- print(err)
- continue # no .envl.dat file?
- for k in unitlayers:
- msu = {'m':k['m'], 's':k['s'], 'u':k['u']}
- msustr = msufmt.format(**msu)
- celltype.loc[msustr]['layer'] = k['layer']
- '''
- # classify as ON/OFF/transient according to movie onset response. Movies are generally
- # darker than mid grey, so a decrease in PSTH after onset is an ON response.
- # After examining lots of raster plots, simplest and most informative may be
- # +/- 0.2 sec on either side of movie onset, after taking into account 0.05 s propagation
- # delay, i.e. from -0.15 to 0 as the spont (pre) rate, and from 0.05 to 0.25 as the evoked
- # (post) rate. Taking just a short period of time around the movie onset helps reduce the
- # impact of movie motion on the result
- print('Classifying ON/OFF/transient')
- resp = mviresp.xs('nat', level='kind') # keep only nat movie kind
- resp = resp.xs('none', level='st8') # keep only none run state
- resp = resp.xs(False, level='opto') # keep only control opto state
- for mseustr, row in resp.iterrows():
- print(mseustr)
- wpsth = row['wpsth']
- wt = row['wt']
- if np.isnan(wt).any():
- continue # skip this mseustr
- prei0, prei1 = wt.searchsorted([PDT-WINDT, EXCLTR[0]])
- posti0, posti1 = wt.searchsorted([EXCLTR[1], PDT+WINDT])
- mvionsetpsth = wpsth[prei0:posti1]
- t = wt[prei0:posti1]
- prepsth = wpsth[prei0:prei1]
- postpsth = wpsth[posti0:posti1]
- insidepsth = wpsth[prei1:posti0]
- outsidepsth = np.concatenate([prepsth, postpsth])
- pre = prepsth.mean()
- post = postpsth.mean()
- ins = insidepsth.mean()
- out = outsidepsth.mean()
- if pre == 0.0 or post == 0.0 or ins == 0.0 or out == 0.0:
- print('No pre/post/inside/outside firing rate, skipping')
- continue
- onoff = (pre - post)/(pre + post)
- trans = (ins - out)/(ins + out)
- Z = max(prepsth.max(), postpsth.max()) # normalization factor, ignore transient period
- msustr = mseustr2msustr(mseustr)
- celltyperow = celltype.loc[msustr]
- oldonoff = celltyperow['onoff']
- oldtrans = celltyperow['trans']
- if np.isnan(oldonoff) or (abs(onoff) > abs(oldonoff)):
- celltyperow['onoff'] = onoff # update
- celltyperow['mvionsetpsth'] = mvionsetpsth / Z # normalize for plotting, update
- celltyperow['t'] = t # update
- if np.isnan(oldtrans) or (abs(trans) > abs(oldtrans)):
- celltyperow['trans'] = trans # update
- ## code for avoiding annoying dataframe pointer vs copy issues:
- '''
- # overwriting/resetting values of an existing column:
- celltype['sbc'][:] = np.nan # I think the explicit slice makes the difference
- # code for adding a new column to an existing dataframe:
- celltype = celltype.copy()
- celltype['dsi'] = np.nan
- celltype['dsi'] = celltype['dsi'].astype(object)
- celltype = celltype.copy()
- celltype['onoff'] = np.nan
- celltype['trans'] = np.nan
- celltype['mvionsetpsth'] = np.nan
- celltype['t'] = np.nan
- celltype['onoff'] = celltype['onoff'].astype(object)
- celltype['trans'] = celltype['trans'].astype(object)
- celltype['mvionsetpsth'] = celltype['mvionsetpsth'].astype(object)
- celltype['t'] = celltype['t'].astype(object)
- '''
- ## collect MUA envelope RF positions for each msu, to use in calculating how well centered
- ## the screen was on the RF of each unit:
- index = pd.Index(mvigrtmsustrs, name='msu')
- columns = ['x0', 'y0']
- cellscreenpos = pd.DataFrame(index=index, columns=columns)
- print('Getting RF params from MUA envelope')
- snexpid = -1 # use last one in each series
- MUAENVLRSQTHRESH = 0.2 # consider only decent MUA envelope RF fits
- series2rffitdf = {} # collect all RF fit params for all series
- for series in mviseries:
- erows = Series.Experiment() & series & 'e_name LIKE "%%Noise%%"'
- expvals = erows.fetch(dj.key, 'e_name', 'e_displayazim', 'e_displayelev')
- snkeys, enames, azims, elevs = expvals
- snkey, ename, azim, elev = (snkeys[snexpid], enames[snexpid],
- azims[snexpid], elevs[snexpid])
- print('Sparse noise exp:', snkey, ename)
- rfs, axes, chans = rfs_for_cdata(snkey)
- #simple_plot_rfs(snkey, rfs, axes, chans, interpolation=None, nrows=4)
- fitdf = fit_MUA_RFs(snkey, rfs, chans)
- msestr = key2msestr(snkey)
- msstr = mseustr2mstr(msestr)
- series2rffitdf[msstr] = fitdf # save
- chans = fitdf.index.levels[0].to_numpy() # all the chans in the df
- chan2screenpos = {}
- for chan in chans:
- chanrow = fitdf.loc[chan, 'mix']
- if chanrow['rsq'] < MUAENVLRSQTHRESH:
- continue
- # remove azim and elev correction to get RF position wrt screen center:
- chan2screenpos[chan] = chanrow['x0']-azim, chanrow['y0']-elev
- urows = up & series
- uids, maxchans = urows.fetch('u', 'u_maxchan')
- for uid, maxchan in zip(uids, maxchans):
- if maxchan in chan2screenpos:
- screenpos = chan2screenpos[maxchan]
- msu = {'m':series['m'], 's':series['s'], 'u':uid}
- msustr = msufmt.format(**msu)
- cellscreenpos.loc[msustr] = screenpos
- ## TODO: somehow convert series2rffitdf dict to DataFrame multiindexed by msstr, chan, rftype:
- #muarfs = pd.DataFrame()
- '''
- # test reliability measure using sparse trials with randomly placed spikes:
- ntrials, nt, nspikes = 200, 1000, 50
- trials = np.zeros((ntrials, nt))
- for trial in trials:
- for spikei in range(nspikes):
- i = np.random.randint(0, nt)
- trial[i] = 1
- print('reliability:', reliability(trials)[0])
- '''
- '''
- # how to get pairs of ON and OFF values from mviresp instead of resp:
- # multiindex row slicing is confusing, need the axis=0 for some reason,
- # see https://pandas.pydata.org/pandas-docs/stable/advanced.html#using-slicers
- li = mviresp.loc(axis=0)
- sparsons = np.float64(li[:, st8, True]['spars'].values)
- sparsoffs = np.float64(li[:, st8, False]['spars'].values)
- nanis = np.isnan(sparsons)
- if not (np.isnan(sparsoffs) == nanis).all():
- raise ValueError("ON and OFF conditions have different nan indices")
- sparsons = sparsons[~nanis]
- sparsoffs = sparsoffs[~nanis]
- '''
- ## TODO: now that mseu are the row indicies instead of msu, they're all unique, and
- ## the group kwarg isn't doing any grouping, I think
- ## TODO: I think the MLM stats don't work because currently the various different MVIKINDS
- ## and ST8S are disabled, only nat and none are being used. So comment this out for now:
- '''
- print('ALL MVI MixedLM')
- mvirespr = mviresp.reset_index() # make a copy, convert mi to columns
- mvirespr = mvirespr.dropna() # drop rows with nans, but not really necessary
- mvirespr = mvirespr[mvirespr['st8'] != 'none'] # exclude 'none' state
- for col in ['meanrate', 'spars', 'rel', 'meanburstratio']:
- # each of these columns are object arrs for some reason, instead of float, which
- # confuses statsmodels and prevents printout of summary. See:
- # https://stackoverflow.com/questions/29799161/summary-not-working-for-ols-estimation
- mvirespr[col] = mvirespr[col].astype(np.float64) # convert each column to float
- formulas = ('meanrate ~ kind * opto * st8',
- 'spars ~ kind * opto * st8',
- 'rel ~ kind * opto * st8',
- 'meanburstratio ~ kind * opto * st8')
- for formula in formulas:
- # group kwarg specifies how data should be treated, prevents measures from different units
- # from being lumped together:
- mlm = smf.mixedlm(formula=formula, data=mvirespr, groups=mvirespr['mseu'])
- results = mlm.fit()
- print(formula)
- print('ngroups:', results.mlm.n_groups)
- paramp = pd.concat([results.params, results.pvalues], axis=1,
- keys=['param', 'p']) # Dataframe from 2 Series
- print(paramp)
- #print(results.summary())
- print()
- print('----------------------------------------')
- print()
- print('NAT only MixedLM')
- mvirespr = mviresp.reset_index() # make a copy, convert mi to columns
- natrespr = mvirespr[mvirespr['kind'] == 'nat']
- natrespr = natrespr.dropna() # drop rows with nans, but not really necessary
- natrespr = natrespr[natrespr['st8'] != 'none'] # exclude 'none' state
- for col in ['meanrate', 'spars', 'rel', 'meanburstratio']:
- # each of these columns are object arrs for some reason, instead of float, which
- # confuses statsmodels and prevents printout of summary. See:
- # https://stackoverflow.com/questions/29799161/summary-not-working-for-ols-estimation
- natrespr[col] = natrespr[col].astype(np.float64) # convert each column to float
- formulas = ('meanrate ~ opto * st8',
- 'spars ~ opto * st8',
- 'rel ~ opto * st8',
- 'meanburstratio ~ opto * st8')
- for formula in formulas:
- # group kwarg specifies how data should be treated, prevents measures from different units
- # from being lumped together:
- mlm = smf.mixedlm(formula=formula, data=natrespr, groups=natrespr['mseu'])
- results = mlm.fit()
- print(formula)
- print('ngroups:', results.mlm.n_groups)
- paramp = pd.concat([results.params, results.pvalues], axis=1,
- keys=['param', 'p']) # Dataframe from 2 Series
- print(paramp)
- #print(results.summary())
- print()
- print('----------------------------------------')
- print()
- print('PNK only MixedLM')
- mvirespr = mviresp.reset_index() # make a copy, convert mi to columns
- pnkrespr = mvirespr[mvirespr['kind'] == 'pnk']
- pnkrespr = pnkrespr.dropna() # drop rows with nans, but not really necessary
- pnkrespr = pnkrespr[pnkrespr['st8'] != 'none'] # exclude 'none' state
- for col in ['meanrate', 'spars', 'rel', 'meanburstratio']:
- # each of these columns are object arrs for some reason, instead of float, which
- # confuses statsmodels and prevents printout of summary. See:
- # https://stackoverflow.com/questions/29799161/summary-not-working-for-ols-estimation
- pnkrespr[col] = pnkrespr[col].astype(np.float64) # convert each column to float
- formulas = ('meanrate ~ opto * st8',
- 'spars ~ opto * st8',
- 'rel ~ opto * st8',
- 'meanburstratio ~ opto * st8')
- for formula in formulas:
- # group kwarg specifies how data should be treated, prevents measures from different units
- # from being lumped together:
- mlm = smf.mixedlm(formula=formula, data=pnkrespr, groups=pnkrespr['mseu'])
- results = mlm.fit()
- print(formula)
- print('ngroups:', results.mlm.n_groups)
- paramp = pd.concat([results.params, results.pvalues], axis=1,
- keys=['param', 'p']) # Dataframe from 2 Series
- print(paramp)
- #print(results.summary())
- print()
- print('----------------------------------------')
- print()
- print('SHF only MixedLM')
- mvirespr = mviresp.reset_index() # make a copy, convert mi to columns
- shfrespr = mvirespr[mvirespr['kind'] == 'shf']
- shfrespr = shfrespr.dropna() # drop rows with nans, but not really necessary
- shfrespr = shfrespr[shfrespr['st8'] != 'none'] # exclude 'none' state
- for col in ['meanrate', 'spars', 'rel', 'meanburstratio']:
- # each of these columns are object arrs for some reason, instead of float, which
- # confuses statsmodels and prevents printout of summary. See:
- # https://stackoverflow.com/questions/29799161/summary-not-working-for-ols-estimation
- shfrespr[col] = shfrespr[col].astype(np.float64) # convert each column to float
- formulas = ('meanrate ~ opto * st8',
- 'spars ~ opto * st8',
- 'rel ~ opto * st8',
- 'meanburstratio ~ opto * st8')
- for formula in formulas:
- # group kwarg specifies how data should be treated, prevents measures from different units
- # from being lumped together:
- mlm = smf.mixedlm(formula=formula, data=shfrespr, groups=shfrespr['mseu'])
- results = mlm.fit()
- print(formula)
- print('ngroups:', results.mlm.n_groups)
- paramp = pd.concat([results.params, results.pvalues], axis=1,
- keys=['param', 'p']) # Dataframe from 2 Series
- print(paramp)
- #print(results.summary())
- print()
- print('----------------------------------------')
- print()
- print('Grating MixedLM')
- grtrespr = grtresp.reset_index() # make a copy, convert mi to columns
- grtrespr = grtrespr.dropna() # drop rows with nans, but not really necessary
- grtrespr = grtrespr[grtrespr['st8'] != 'none'] # exclude 'none' state
- for col in ['meanrate', 'meanburstratio']:
- # each of these columns are object arrs for some reason, instead of float, which
- # confuses statsmodels and prevents printout of summary. See:
- # https://stackoverflow.com/questions/29799161/summary-not-working-for-ols-estimation
- grtrespr[col] = grtrespr[col].astype(np.float64) # convert each column to float
- formulas = ('meanrate ~ opto * st8',
- 'meanburstratio ~ opto * st8')
- for formula in formulas:
- # group kwarg specifies how data should be treated, prevents measures from different units
- # from being lumped together:
- mlm = smf.mixedlm(formula=formula, data=grtrespr, groups=grtrespr['mseu'])
- results = mlm.fit()
- print(formula)
- print('ngroups:', results.mlm.n_groups)
- paramp = pd.concat([results.params, results.pvalues], axis=1,
- keys=['param', 'p']) # Dataframe from 2 Series
- print(paramp)
- #print(results.summary())
- print()
- print('----------------------------------------')
- print()
- '''
|