Browse Source

Clean up irrelevant comments & code blocks, fix .csv filenames

Martin Spacek 2 years ago
parent
commit
08505455c2
15 changed files with 46 additions and 2372 deletions
  1. 0 14
      README.md
  2. 5 20
      fig1.py
  3. 20 1362
      fig1S33S1.py
  4. 0 298
      fig1S33S1_ntsr.py
  5. 0 13
      fig1S6.py
  6. 1 8
      fig2.py
  7. 0 152
      fig3.py
  8. 1 302
      fig4.py
  9. 0 8
      fig4S1.py
  10. 6 12
      fig5.py
  11. 0 2
      fig5S2.py
  12. 1 123
      fig6.py
  13. 2 28
      ipos.py
  14. 10 11
      main.py
  15. 0 19
      util.py

+ 0 - 14
README.md

@@ -36,20 +36,6 @@ run -i fig4S1.py
 run -i ipos.py
 ```
 
-Here is the mapping between the old and new supplemental figure names under
-eLife's parent/child naming scheme. This applies to both the names in the code and the paper:
-
-| Old     | New |
-| ------- | --- |
-| S1      | 1S1 |
-| S2      | 1S2 |
-| S3      | 4S1 |
-| S4a-j   | 1S3 |
-| S4k-t   | 3S1 |
-| S5      | 1S4 |
-| negntsr | 1S5 |
-| S6      | 5S1 |
-
 Note that depending on your system, you may have to call `plt.show()` to actually display
 figure windows after running one or more of the above scripts. Use `cf()` to quickly close
 all figures.

+ 5 - 20
fig1.py

@@ -5,8 +5,6 @@ mi = pd.MultiIndex.from_product([mvimseustrs, OPTOS], names=['mseu', 'opto'])
 columns = ['trialis', 'rates', 'rate02s', 'rate35s', 'burstratios',
            'blankrates', 'blankburstratios'] + modmeasuresnoblankcond
 fig1 = pd.DataFrame(index=mi, columns=columns) # excludes kind and st8
-## NOTE: make sure that rows are not overwritten due to temporarily enabling iteration
-## over kind and/or st8 in any of the following loops!
 
 
 # plot movie rasters for each unit by kind, state and opto condition:
@@ -14,8 +12,6 @@ showbursts = True
 SCATTER = True # True: low quality, but fast; False: high quality, but slow
 trialiis = None #np.arange(120) # plot only these trials, None plots all trials
 for mseustr in []:#mvimseu2exmpli: #mvimseustrs:
-    #['Ntsr1Cre_2019_0008_s03_e04_u18']
-    #['Ntsr1Cre_2020_0001_s04_e10_u12']
     for kind in ['nat']: #MVIKINDS
         for st8 in ['none']:
             for opto in OPTOS:
@@ -311,10 +307,6 @@ for kind in ['nat']:#MVIKINDS:
         rons = np.asarray(rons)
         roffs = np.asarray(roffs)
         sgnfs = np.asarray(sgnfs)
-        # replace off-scale low values with log0min, so the points remain visible:
-        #pltrons, pltroffs = rons.copy(), roffs.copy()
-        #pltrons[pltrons <= 10**logmin] = 10**log0min
-        #pltroffs[pltroffs <= 10**logmin] = 10**log0min
         # plot y=x line:
         xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
         a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
@@ -682,9 +674,6 @@ for kind in ['nat']:#MVIKINDS:
 
 
 # plot movie PSTH peak width distributions:
-## TODO: these width measures should maybe be normalized separately for each unit,
-## so that the effect of feedback suppression isn't drowned out by the variability of
-## PSTH peak widths across units
 normalizepeakwidths = True
 figsize = DEFAULTFIGURESIZE
 for kind in ['nat']:#MVIKINDS:
@@ -837,15 +826,15 @@ for kind in ['nat']:#MVIKINDS:
                           marker=marker, c=c, s=sz, lw=lw)
             # get fname of appropriate LMM .cvs file:
             if measure == 'meanburstratio': # use Steffen's LMM linregress fit, for fig1S2
-                fname = os.path.join('stats', 'figure_1S2c_coefs.csv')
+                fname = os.path.join('stats', 'figure_1_S2c_coefs.csv')
             elif measure == 'spars': # use Steffen's LMM linregress fit, for fig1S2
-                fname = os.path.join('stats', 'figure_1S2d_coefs.csv')
+                fname = os.path.join('stats', 'figure_1_S2d_coefs.csv')
             elif measure == 'rel': # use Steffen's LMM linregress fit, for fig1S2
-                fname = os.path.join('stats', 'figure_1S2e_coefs.csv')
+                fname = os.path.join('stats', 'figure_1_S2e_coefs.csv')
             elif measure == 'snr': # use Steffen's LMM linregress fit, for fig1S2
-                fname = os.path.join('stats', 'figure_1S2f_coefs.csv')
+                fname = os.path.join('stats', 'figure_1_S2f_coefs.csv')
             elif measure == 'meanpkw': # use Steffen's LMM linregress fit, for fig1S2
-                fname = os.path.join('stats', 'figure_1S2g_coefs.csv')
+                fname = os.path.join('stats', 'figure_1_S2g_coefs.csv')
             else:
                 print("WARNING: No LMM stats for %s measure" % measure)
                 fname = None
@@ -916,10 +905,6 @@ for meanratestr, rangestr in zip(['meanrate02', 'meanrate35'], ['0-2', '3-5']):
     rons = np.asarray(rons)
     roffs = np.asarray(roffs)
     sgnfs = np.asarray(sgnfs)
-    # replace off-scale low values with log0min, so the points remain visible:
-    #pltrons, pltroffs = rons.copy(), roffs.copy()
-    #pltrons[pltrons <= 10**logmin] = 10**log0min
-    #pltroffs[pltroffs <= 10**logmin] = 10**log0min
     # plot y=x line:
     xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
     a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)

File diff suppressed because it is too large
+ 20 - 1362
fig1S33S1.py


+ 0 - 298
fig1S33S1_ntsr.py

@@ -1,298 +0,0 @@
-"""Figure 1S33S1 Ntsr versions of fig1S33S1.py, use run -i fig1S33S1_ntsr.py. Currently used
-for rebuttal figure R2"""
-
-assert EXPTYPE == 'ntsrmvis'
-
-from matplotlib.patches import Rectangle
-
-
-
-chirptype = pd.read_csv('csv/ntsr_unit_celltypes_yannik.csv')
-chirptype.set_index(['msu'], inplace=True)
-
-
-## strip plot max FMI of movie and grating meanrate and meanburstratio,
-## each msu classified as chirp ON/OFF/transient/mixed by Yannik:
-np.random.seed(0) # to get identical horizontal jitter in strip plots on every run
-figsize = DEFAULTFIGURESIZE
-chirp2clr = {'On':green, 'Off':'red', 'Transient':'black', 'Mixed':'gray'}
-for stimtype in STIMTYPES:
-    stimtypelabel = stimtype2axislabel[stimtype]
-    for measure in ['meanrate', 'meanburstratio']:
-        axislabel = measure2axislabel[measure]
-        if axislabel.islower():
-            axislabel = axislabel.capitalize()
-        fmis, onoffs = [], []
-        for msustr in mvigrtmsustrs:
-            fmi = maxFMI.loc[msustr, 'none', stimtype][measure]
-            if pd.isna(fmi):
-                continue
-            try:
-                onoff = chirptype.loc[msustr]['clu_label']
-            except KeyError:
-                continue # no chirp cell type for this msu
-            fmis.append(fmi)
-            onoffs.append(onoff)
-        fmis = np.asarray(fmis)
-        onoffs = np.asarray(onoffs)
-        onfmis = fmis[onoffs == 'ON-sust.']
-        offfmis = fmis[onoffs == 'OFF-sust.']
-        onofffmis = fmis[onoffs == 'ON-OFF-trans.']
-        mixedfmis = fmis[onoffs == 'mixed']
-        f, a = plt.subplots(figsize=figsize)
-        wintitle('maxFMI chirp %s %s stripplot' % (stimtypelabel, measure))
-        dfdict = {'On':onfmis, 'Off':offfmis, 'Transient':onofffmis, 'Mixed':mixedfmis}
-        data = pd.DataFrame.from_dict(dfdict, orient='index').transpose()
-        sns.stripplot(ax=a, data=data, palette=chirp2clr, clip_on=False,
-                      marker='.', size=np.sqrt(70))
-        a.set_ylabel('FMI')
-        a.set_xlim(0, 3)
-        a.set_ylim(-1, 1)
-        a.set_yticks([-1, 0, 1])
-        a.spines['left'].set_position(('outward', 10))
-        a.spines['bottom'].set_position(('outward', 5))
-        #a.spines['bottom'].set_visible(False)
-        #a.tick_params(bottom=False)
-        plt.xticks(rotation=45)
-        print(stimtype, measure)
-        print(data)
-
-
-## scatter plot grating meanrates, coloured by chirp ON sustained/OFF sustained/transient/mixed
-## as classified by Yannik:
-figsize = DEFAULTFIGURESIZE
-logmin, logmax = -1.2, 2
-logticks = np.array([-1, 0, 1, 2])
-#log0min = logmin + 0.05
-for st8 in ['none']:#ALLST8S:
-    f, a = plt.subplots(figsize=figsize)
-    wintitle('opto meanrate grating %s chirp' % st8)
-    rons, roffs, onoffs = [], [], []
-    for mseustr in grtmseustrs:
-        meanrate = grtresp.loc[mseustr, st8]['meanrate']
-        msustr = mseustr2msustr(mseustr)
-        try:
-            onoff = chirptype.loc[msustr]['clu_label']
-        except KeyError:
-            continue # no chirp experiment for cell typing
-        if meanrate.isna().any(): # missing one or both meanrates
-            #print('%s: missing one or both opto conditions, skipping' % mseustr)
-            continue
-        rons.append(meanrate[True])
-        roffs.append(meanrate[False])
-        onoffs.append(onoff)
-        #fig3.loc[mseustr, 'meanrate'] = meanrate[False], meanrate[True] # save
-    rons = np.asarray(rons)
-    roffs = np.asarray(roffs)
-    onoffs = np.asarray(onoffs)
-    clrs = []
-    chirptype2clr = {'ON-sust.':green, 'OFF-sust.':'red', 'ON-OFF-trans.':'black',
-                     'mixed':'gray'}
-    for onoff in onoffs:
-        clr = chirptype2clr[onoff]
-        clrs.append(clr)
-    # plot y=x line:
-    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
-    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
-    # plot all points:
-    a.scatter(rons, roffs, marker='.', c='None', edgecolor=clrs, s=DEFSZ)
-    a.set_ylabel('Feedback FR (spk/s)')
-    a.set_xlabel('Suppression FR (spk/s)')
-    a.set_xscale('log')
-    a.set_yscale('log')
-    a.set_xlim(10**logmin, 10**logmax)
-    a.set_ylim(10**logmin, 10**logmax)
-    a.set_xticks(10.0**logticks)
-    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
-    axes_disable_scientific(a)
-    a.minorticks_off()
-    a.set_aspect('equal')
-
-
-## scatter plot movie meanrates, coloured by chirp ON sustained/OFF sustained/transient/mixed
-## as classified by Yannik:
-figsize = DEFAULTFIGURESIZE
-logmin, logmax = -1.2, 2
-logticks = np.array([-1, 0, 1, 2])
-#log0min = logmin + 0.05
-for st8 in ['none']:#ALLST8S:
-    f, a = plt.subplots(figsize=figsize)
-    wintitle('opto meanrate movie nat %s chirp' % st8)
-    rons, roffs, onoffs = [], [], []
-    for mseustr in mvimseustrs:
-        meanrate = mviresp.loc[mseustr, 'nat', st8]['meanrate']
-        msustr = mseustr2msustr(mseustr)
-        try:
-            onoff = chirptype.loc[msustr]['clu_label']
-        except KeyError:
-            continue # no chirp experiment for cell typing
-        if meanrate.isna().any(): # missing one or both meanrates
-            #print('%s: missing one or both opto conditions, skipping' % mseustr)
-            continue
-        print(mseustr)
-        rons.append(meanrate[True])
-        roffs.append(meanrate[False])
-        onoffs.append(onoff)
-    rons = np.asarray(rons)
-    roffs = np.asarray(roffs)
-    clrs = []
-    chirptype2clr = {'ON-sust.':green, 'OFF-sust.':'red', 'ON-OFF-trans.':'black',
-                     'mixed':'gray'}
-    for onoff in onoffs:
-        clr = chirptype2clr[onoff]
-        clrs.append(clr)
-    # plot y=x line:
-    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
-    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
-    # plot all points:
-    a.scatter(rons, roffs, marker='.', c='None', edgecolor=clrs, s=DEFSZ)
-    a.set_ylabel('Feedback FR (spk/s)')
-    a.set_xlabel('Suppression FR (spk/s)')
-    a.set_xscale('log')
-    a.set_yscale('log')
-    a.set_xlim(10**logmin, 10**logmax)
-    a.set_ylim(10**logmin, 10**logmax)
-    a.set_xticks(10.0**logticks)
-    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
-    axes_disable_scientific(a)
-    a.minorticks_off()
-    a.set_aspect('equal')
-
-
-## scatter plot grating burst ratio, coloured by chirp ON sustained/OFF sustained/transient/mixed
-## as classified by Yannik:
-figsize = DEFAULTFIGURESIZE
-logmin, logmax = -3.55, 0
-logticks = np.array([-3, -2, -1, 0])
-for st8 in ['none']:#ALLST8S:
-    f, a = plt.subplots(figsize=figsize)
-    wintitle('opto burst ratio grating %s chirp' % st8)
-    brons, broffs, onoffs = [], [], []
-    for mseustr in grtmseustrs:
-        br = grtresp.loc[mseustr, st8]['meanburstratio']
-        msustr = mseustr2msustr(mseustr)
-        try:
-            onoff = chirptype.loc[msustr]['clu_label']
-        except KeyError:
-            continue # no chirp experiment for cell typing
-        if br.isna().any(): # missing for at least one opto condition
-            continue
-        brons.append(br[True])
-        broffs.append(br[False])
-        onoffs.append(onoff)
-        #fig3.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
-    brons, broffs, onoffs = np.asarray(brons), np.asarray(broffs), np.asarray(onoffs)
-    clrs = []
-    chirptype2clr = {'ON-sust.':green, 'OFF-sust.':'red', 'ON-OFF-trans.':'black',
-                     'mixed':'gray'}
-    for onoff in onoffs:
-        clr = chirptype2clr[onoff]
-        clrs.append(clr)
-    # plot y=x line:
-    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
-    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
-    # plot all points:
-    a.scatter(brons, broffs, marker='.', c='None', edgecolor=clrs, s=DEFSZ)
-    a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
-    a.set_ylabel('Feedback BR')
-    a.set_xscale('log')
-    a.set_yscale('log')
-    a.set_xlim(10**logmin, 10**logmax)
-    a.set_ylim(10**logmin, 10**logmax)
-    a.set_xticks(10.0**logticks)
-    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
-    a.minorticks_off()
-    axes_disable_scientific(a)
-    a.set_aspect('equal')
-
-
-## scatter plot movie burst ratio, coloured by chirp ON sustained/OFF sustained/transient/mixed
-## as classified by Yannik:
-figsize = DEFAULTFIGURESIZE
-logmin, logmax = -3.55, 0
-logticks = np.array([-3, -2, -1, 0])
-for st8 in ['none']:#ALLST8S:
-    f, a = plt.subplots(figsize=figsize)
-    wintitle('opto burst ratio movie nat %s chirp' % st8)
-    brons, broffs, onoffs = [], [], []
-    for mseustr in mvimseustrs:
-        br = mviresp.loc[mseustr, 'nat', st8]['meanburstratio']
-        msustr = mseustr2msustr(mseustr)
-        msustr = mseustr2msustr(mseustr)
-        try:
-            onoff = chirptype.loc[msustr]['clu_label']
-        except KeyError:
-            continue # no chirp experiment for cell typing
-        if br.isna().any(): # missing for at least one opto condition
-            continue
-        brons.append(br[True])
-        broffs.append(br[False])
-        onoffs.append(onoff)
-        #fig3.loc[mseustr, 'meanburstratio'] = br[False], br[True] # save
-    brons, broffs, onoffs = np.asarray(brons), np.asarray(broffs), np.asarray(onoffs)
-    clrs = []
-    chirptype2clr = {'ON-sust.':green, 'OFF-sust.':'red', 'ON-OFF-trans.':'black',
-                     'mixed':'gray'}
-    for onoff in onoffs:
-        clr = chirptype2clr[onoff]
-        clrs.append(clr)
-    # plot y=x line:
-    xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
-    a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
-    # plot all points:
-    a.scatter(brons, broffs, marker='.', c='None', edgecolor=clrs, s=DEFSZ)
-    a.set_xlabel('Suppression BR') # keep it short to maximize space for axes
-    a.set_ylabel('Feedback BR')
-    a.set_xscale('log')
-    a.set_yscale('log')
-    a.set_xlim(10**logmin, 10**logmax)
-    a.set_ylim(10**logmin, 10**logmax)
-    a.set_xticks(10.0**logticks)
-    a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
-    a.minorticks_off()
-    axes_disable_scientific(a)
-    a.set_aspect('equal')
-
-
-#TRANSTHRESH = 0.2
-#ONOFFTHRESH = 0.2
-
-
-## generate celltype vs chirptype on/off/transient/other matrix, requires that celltype code
-## from fig1S33S1.py has been run with EXPTYPE == 'ntsrmvis' set.
-## The two classification methods correspond very poorly:
-# join on msu index, those rows that don't overlap have chirp type fields left as NaN
-celltypechirptype = celltype.join(chirptype)
-mat = np.zeros((4, 4), dtype=int)
-chrptypestr2i = {'ON-sust.':0, 'OFF-sust.':1, 'ON-OFF-trans.':2, 'mixed':3}
-mvitypestr2i  = {'On':0, 'Off':1, 'Trans':2, '':3}
-for msustr, row in celltypechirptype.iterrows():
-    chrptypestr = row['clu_label']
-    if pd.isna(chrptypestr):
-        continue # mvi unit doesn't have a chirp experiment?
-    onoff = row['onoff']
-    trans = row['trans']
-    mvitypestr = ''
-    if trans >= TRANSTHRESH: # first try and classify by trans
-        mvitypestr = 'Trans'
-    elif onoff >= ONOFFTHRESH: # try and classify by on/off
-        mvitypestr = 'On'
-    elif onoff < -ONOFFTHRESH:
-        mvitypestr = 'Off'
-    chrptypei = chrptypestr2i[chrptypestr]
-    mvitypei = mvitypestr2i[mvitypestr]
-    mat[mvitypei, chrptypei] += 1
-figsize = 4, 3
-f, a = plt.subplots(figsize=figsize)
-wintitle('celltype vs chirptype matrix ONOFFTHRESH=%.2f TRANSTHRESH=%.2f'
-         % (ONOFFTHRESH, TRANSTHRESH))
-#a.imshow(mat, origin='upper', aspect='equal', cmap='gray')
-img = a.matshow(mat)
-#a.set_xticks(np.arange(4))
-#a.set_yticks(np.arange(4))
-#a.set_xticklabels(list(chrptypestr2i.keys()))
-#a.set_yticklabels(list(mvitypestr2i.keys()))
-f.colorbar(img, ticks=[mat.min(), mat.max()], label='Unit count')
-a.set_xlabel('Chirp type')
-a.set_ylabel('Movie onset type')
-f.tight_layout(pad=1)

+ 0 - 13
fig1S6.py

@@ -25,10 +25,6 @@ plt.show()
 print("Now save the first 4 scatter plot figures, but don't save any dataframes")
 
 
-## TODO: for reviewer's interest, also calculate average normalized pupil area trace
-## triggered off opto onset, separately for Ntsr and PV
-
-
 ## scatter plot pupil area during feedback vs suppression, one point (mean over trials) per mse:
 for stimtype, stimtypei in zip(STIMTYPESPLUSALL, STIMTYPESPLUSALLI):
     linmin, linmax = 0, 0.6
@@ -101,7 +97,6 @@ for stimtype, stimtypei in zip(STIMTYPESPLUSALL, STIMTYPESPLUSALLI):
         areafmi = (feedback - suppress) / (feedback + suppress)
         mse2areafmi[umse] = areafmi
     # now collect neural FMIs for each mseu:
-    ## TODO: calculating neural FMIs over all trials, but could also choose only run or sit trials
     FMI = stimtype2FMI[stimtype]
     st8 = 'none'
     mseus = np.unique(FMI.reset_index()['mseu'].values)
@@ -152,10 +147,6 @@ for stimtype, stimtypei in zip(STIMTYPESPLUSALL, STIMTYPESPLUSALLI):
     a.set_yticks([ymin, 0, ymax])
     a.spines['left'].set_position(('outward', 4))
     a.spines['bottom'].set_position(('outward', 4))
-    #txt = ('$\mathregular{p=%.2g}$\n' # prob that slope is 0
-    #       '$\mathregular{R^{2}=%.2g}$\n'
-    #       '$\mathregular{slope=%.2g}$' % (p, rsq, mm))
-    #a.add_artist(AnchoredText(txt, loc='lower left', frameon=False))
     ## scatter plot meanburstratio FMI vs pupil area FMI:
     f, a = plt.subplots(figsize=figsize)
     wintitle('meanburstratio %s FMI vs pupil area FMI' % stimtype)
@@ -192,10 +183,6 @@ for stimtype, stimtypei in zip(STIMTYPESPLUSALL, STIMTYPESPLUSALLI):
     a.set_yticks([ymin, 0, ymax])
     a.spines['left'].set_position(('outward', 4))
     a.spines['bottom'].set_position(('outward', 4))
-    #txt = ('$\mathregular{p=%.2g}$\n' # prob that slope is 0
-    #       '$\mathregular{R^{2}=%.2g}$\n'
-    #       '$\mathregular{slope=%.2g}$' % (p, rsq, mm))
-    #a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))
     # save to iposmi dataframe:
     if stimtype == 'mvi+grt':
         continue # don't overwrite stimtype column in iposmi, other values should be identical

+ 1 - 8
fig2.py

@@ -12,8 +12,6 @@ linmin = 0
 
 
 # plot single-sample V1 model-fit PSTH overtop of test feedback and suppression PSTHs:
-## TODO: correlate burst PSTHs with (model-empirical)**2, scatter plot them against each other,
-## or just plot distribution of corrcoefs
 dt = 5 # stimulus duration sec
 #tmin, tmax = 0 + OFFSETS[0], dt + OFFSETS[1]
 tmin, tmax = 0, dt
@@ -185,7 +183,7 @@ for kind in ['nat']:#MVIKINDS:
             for i, row in rows.iterrows(): # for each mseu
                 mseustr = row['mseu']
                 ms, ths, rsqs = row[fmode+'ms'], row[fmode+'ths'], row[fmode+'rsqs']
-                if np.isnan([ms, ths, rsqs]).any(): ## TODO: does this throw out too many units?
+                if np.isnan([ms, ths, rsqs]).any():
                     continue # skip this mseu
                 # skip mseus that had no meaningful signal to fit in either condition:
                 maxsnr = get_max_snr(mvirespr, mseustr, kind, st8)
@@ -194,8 +192,6 @@ for kind in ['nat']:#MVIKINDS:
                 mmed, ml, mh = percentile_ci(ms)
                 thmed, thl, thh = percentile_ci(ths)
                 rsqmed = np.median(rsqs)
-                #if rsqmed < RSQTHRESH: # skip mseus with really bad median fits
-                #    continue # skip
                 if rsqmed < 0: # catastrophic fit failure
                     print('catastrophic fit failure:', mseustr)
                     continue # skip
@@ -392,7 +388,6 @@ for kind in ['nat']:#MVIKINDS:
         for i, row in rows.iterrows(): # for each mseu
             mseustr = row['mseu']
             allrsqs, nbrsqs, nrrsqs = row['rsqs'], row['nbrsqs'], row['nrrsqs']
-            ## TODO: does this throw out too many units?:
             if np.isnan([allrsqs, nbrsqs, nrrsqs]).any():
                 continue # skip this mseu
             # skip mseus that had no meaningful signal to fit in either condition:
@@ -402,8 +397,6 @@ for kind in ['nat']:#MVIKINDS:
             allrsqmed = np.median(allrsqs)
             nbrsqmed = np.median(nbrsqs)
             nrrsqmed = np.median(nrrsqs)
-            #if rsqmed < RSQTHRESH: # skip mseus with really bad median fits
-            #    continue # skip
             if allrsqmed < 0 or nbrsqmed < 0 or nrrsqmed < 0: # catastrophic fit failure
                 print('catastrophic fit failure:', mseustr)
                 continue # skip

+ 0 - 152
fig3.py

@@ -8,13 +8,10 @@ columns = ['trialis', 'rates', 'blankrates', 'blankcondrates', # single trials
            'osi', 'oripref', 'f1f0',
            'phi', 'bphi', 'nbphi', 'dphi', 'dbphi', 'dnbphi']
 fig3 = pd.DataFrame(index=mi, columns=columns) # excludes kind and st8
-## NOTE: make sure that rows are not overwritten due to temporarily enabling iteration
-## over st8 in any of the following loops!
 
 
 # plot grating rasters for each unit by state and opto condition:
 showbursts = True
-## TODO: fix hard-coding:
 if EXPTYPE == 'negntsrmvis':
     NTRIALSPERCONDITION = 8 # per stimi and opto combination
     inchespertrial = 1/20*RASTERSCALEX
@@ -24,8 +21,6 @@ else:
 #outliers = ['PVCre_2018_0001_s02_e02_u36', 'PVCre_2018_0001_s02_e06_u36',
 #            'PVCre_2018_0001_s02_e02_u64', 'PVCre_2018_0001_s02_e06_u64']
 for mseustr in []:#grtmseu2exmpli: #grtmseustrs: #outliers:
-    #['Ntsr1Cre_2019_0008_s06_e02_u46']
-    #['Ntsr1Cre_2020_0001_s04_e04_u12'] # negntsr fig
     for st8 in ['none']:#ALLST8S:
         for opto in OPTOS:
             dt = grtresp.loc[mseustr, st8, opto]['dt']
@@ -79,13 +74,7 @@ if tun is not None:
             tunrow.plot(cs=['k', optoblue])
 
 
-## TODO: plot f1/f0 vs. ori for some example units, e.g. Ntsr1Cre_2020_0001_s04_e04_u55
-## Some of the time, or maybe even much of the time, f1/f0 vs. ori might show better
-## orientation tuning than firing rate vs. ori.
-
-
 ## calculate cycle averages (i.e. grating PSTHs), F1/F0, and phi for all units:
-## TODO: fit sinusoids and use fit error to evaluate how linear the responses are?
 PLOTCYCPSTH = False
 CYCPSTHTHRESH = 10 # minimum cycle average peak, Hz
 F1F0THRESH = 0.25 # minimum grating modulation index
@@ -114,9 +103,6 @@ for st8 in ['none']:#ALLST8S:
             print('%s: missing one or both braster opto conditions, skipping' % mseustr)
             continue
         orisopto = grtresp.loc[mseustr, st8]['oris'] # sorted, correspond to rows in raster
-        ## TODO: maybe instead of enforcing this, just pull out all existing oris in both
-        ## opto conditions, and do unique on those, and if any don't exist then handle that
-        ## in the innermost loop?
         assert np.equal(*orisopto).all() # make sure oris are identical across opto
         # get oris from arbitrary condition, convert to int for use as dict keys:
         uoris = np.int64(np.unique(orisopto[False]))
@@ -422,10 +408,6 @@ for st8 in ['none']:#ALLST8S:
     rons = np.asarray(rons)
     roffs = np.asarray(roffs)
     sgnfs = np.asarray(sgnfs)
-    # replace off-scale low values with log0min, so the points remain visible:
-    #pltrons, pltroffs = rons.copy(), roffs.copy()
-    #pltrons[pltrons <= 10**logmin] = 10**log0min
-    #pltroffs[pltroffs <= 10**logmin] = 10**log0min
     # plot y=x line:
     xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
     a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
@@ -604,10 +586,6 @@ for st8 in ['none']:#ALLST8S:
     wintitle('opto osi grating %s' % st8)
     logmin, logmax, logstep = -2.2, 0, 1
     logticks = np.array([-2, -1, 0])
-    # replace off-scale low values with log0min, so the points remain visible:
-    #pltosions, pltosioffs = osions.copy(), osioffs.copy()
-    #pltosions[pltosions <= 10**logmin] = 10**log0min
-    #pltosioffs[pltosioffs <= 10**logmin] = 10**log0min
     # plot y=x line:
     xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
     #xyline = [linmin, linmax], [linmin, linmax]
@@ -691,13 +669,6 @@ for st8 in ['none']:#ALLST8S:
     # plot normal (non-example) points, coloring by state:
     a.scatter(oprons[normlis], oproffs[normlis], clip_on=False,
               marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
-    # plot normal (non-example) points, coloring by min OSI:
-    '''
-    minosi = np.vstack([osions[normlis], osioffs[normlis]]).min(axis=0)
-    scatt = a.scatter(oproffs[normlis], oprons[normlis],
-                      marker='.', c=minosi, vmax=0.1, s=100)
-    f.colorbar(scatt)
-    '''
     # plot example points, one at a time:
     for exmpli, mseustr in zip(exmplis, exmplmseustrs):
         marker = exmpli2mrk[grtmseu2exmpli[mseustr]]
@@ -786,8 +757,6 @@ for st8 in ['none']:#ALLST8S:
     wintitle('opto mean f1 grating %s' % st8)
     # plot y=x line:
     xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
-    #linmin, linmax = 0, 75
-    #xyline = [linmin, linmax], [linmin, linmax]
     a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
     # plot normal (non-example) points:
     a.scatter(f1ons[normlis], f1offs[normlis], clip_on=False,
@@ -854,127 +823,6 @@ for st8 in ['none']:#ALLST8S:
     #t, p = ttest_rel(f1f0ons, f1f0offs) # paired t-test
     #a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
 
-'''
-# scatter plot grating opto local/global max f1/f0 and max f1:
-figsize = DEFAULTFIGURESIZE
-logmin, logmax = -1.5, 2
-logticks = np.array([-1, 0, 1, 2])
-for maxmethod in ['local', 'global']:
-    for st8 in ['none']:#ALLST8S:
-        maxf1f0s, maxf1s, exmplis, exmplmseustrs, normlis = [], [], [], [], []
-        keptmseui = 0 # manually init and increment instead of using enumerate()
-        for mseustr in grtmseustrs:
-            f0s, f1s = [[], []], [[], []]
-            rasters = grtresp.loc[mseustr, st8]['raster']
-            dt = grtresp.loc[mseustr, st8, False]['dt']
-            tfreq = grtresp.loc[mseustr, st8, False]['tfreq']
-            if rasters.isna().any(): # missing one or both rasters
-                print('%s: missing one or both opto conditions, skipping' % mseustr)
-                continue
-            offraster, onraster = rasters[False], rasters[True]
-            for optoi, opto in enumerate(OPTOS):
-                stimis = grtresp.loc[mseustr, st8, opto]['stimis']
-                ustimis = np.unique(stimis)
-                for stimi in ustimis:
-                    trialis = stimis == stimi # boolean array
-                    stimiraster = rasters[opto][trialis] # raster including only stimi's trials
-                    f0, _ = raster2freqcomp(stimiraster, dt, 0, mean='vector')
-                    f1, _ = raster2freqcomp(stimiraster, dt, tfreq, mean='vector')
-                    f0s[optoi].append(f0)
-                    f1s[optoi].append(f1)
-            f0s = np.array(f0s) # 2D array, row0:False, row1:True
-            f0s[f0s == 0] = np.nan # prevent div by 0, replace 0s with nans
-            f1s = np.array(f1s)
-            f1f0s = f1s / f0s # 2D array, row0:False, row1:True
-            if maxmethod == 'local':
-                # get max for each opto condition separately:
-                #maxf1f0 = np.nanmax(f1f0s, axis=1)
-                maxrowis = np.array([0, 1])
-                maxcolis = np.nanargmax(f1f0s, axis=1) # 2 values, one per row
-                maxf1f0 = f1f0s[maxrowis, maxcolis]
-                maxf1 = f1s[maxrowis, maxcolis]
-            elif maxmethod == 'global':
-                # find global max, use it to pick paired value of other opto condition:
-                flatmaxi = np.nanargmax(f1f0s)
-                maxrowi, maxcoli = np.unravel_index(flatmaxi, f1f0s.shape)
-                maxf1f0 = f1f0s[:, maxcoli]
-                maxf1 = f1s[:, maxcoli]
-            maxf1f0s.append(maxf1f0)
-            maxf1s.append(maxf1)
-            if mseustr in grtmseu2exmpli:
-                exmplis.append(keptmseui)
-                exmplmseustrs.append(mseustr)
-            else:
-                normlis.append(keptmseui)
-            keptmseui += 1 # manually increment
-        maxf1f0s = np.asarray(maxf1f0s).T # 2D array: optoi, mseui
-        maxf1s = np.asarray(maxf1s).T # 2D array: optoi, mseui
-        ## scatter plot max f1/f0:
-        f, a = plt.subplots(figsize=figsize)
-        wintitle('opto %s max f1f0 grating %s' % (maxmethod, st8))
-        # plot y=x line:
-        #xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
-        linmin, linmax = 0, 2
-        xyline = [linmin, linmax], [linmin, linmax]
-        a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
-        # plot normal (non-example) points:
-        a.scatter(maxf1f0s[1, normlis], maxf1f0s[0, normlis],
-                  marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
-        # plot example points, one at a time:
-        for exmpli, mseustr in zip(exmplis, exmplmseustrs):
-            marker = exmpli2mrk[grtmseu2exmpli[mseustr]]
-            c = exmpli2clr[grtmseu2exmpli[mseustr]]
-            sz = exmpli2sz[grtmseu2exmpli[mseustr]]
-            a.scatter(maxf1f0s[1, exmpli], maxf1f0s[0, exmpli], marker=marker, c=c, s=sz)
-        a.set_xlabel('Suppression max F1/F0')
-        a.set_ylabel('Feedback max F1/F0')
-        a.set_xlim(linmin, linmax+0.1)
-        a.set_ylim(linmin, linmax+0.1)
-        a.set_xticks(np.arange(linmin, linmax+1))
-        #a.set_xscale('log')
-        #a.set_yscale('log')
-        #a.set_xlim(10**logmin, 10**logmax)
-        #a.set_ylim(10**logmin, 10**logmax)
-        #a.set_xticks(10.0**logticks)
-        a.set_yticks(a.get_xticks()) # make scale y ticks the same as x ticks
-        a.minorticks_off()
-        a.set_aspect('equal')
-        t, p = ttest_rel(maxf1f0s[1], maxf1f0s[0]) # paired t-test
-        a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
-        ## scatter plot max f1:
-        f, a = plt.subplots(figsize=figsize)
-        wintitle('opto %s max f1 grating %s' % (maxmethod, st8))
-        # plot y=x line:
-        xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
-        #linmin, linmax = 0, 2
-        #xyline = [linmin, linmax], [linmin, linmax]
-        a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
-        # plot normal (non-example) points:
-        a.scatter(maxf1s[1, normlis], maxf1s[0, normlis],
-                  marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
-        # plot example points, one at a time:
-        for exmpli, mseustr in zip(exmplis, exmplmseustrs):
-            marker = exmpli2mrk[grtmseu2exmpli[mseustr]]
-            c = exmpli2clr[grtmseu2exmpli[mseustr]]
-            sz = exmpli2sz[grtmseu2exmpli[mseustr]]
-            a.scatter(maxf1s[1, exmpli], maxf1s[0, exmpli], marker=marker, c=c, s=sz)
-        a.set_xlabel('Suppression max F1')
-        a.set_ylabel('Feedback max F1')
-        #a.set_xlim(linmin, linmax+0.1)
-        #a.set_ylim(linmin, linmax+0.1)
-        #a.set_xticks(np.arange(linmin, linmax+1))
-        a.set_xscale('log')
-        a.set_yscale('log')
-        a.set_xlim(10**logmin, 10**logmax)
-        a.set_ylim(10**logmin, 10**logmax)
-        a.set_xticks(10.0**logticks)
-        a.set_yticks(a.get_xticks()) # make scale y ticks the same as x ticks
-        a.minorticks_off()
-        axes_disable_scientific(a)
-        a.set_aspect('equal')
-        t, p = ttest_rel(maxf1s[1], maxf1s[0]) # paired t-test
-        a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper left', frameon=False))
-'''
 
 # fig4S1a: scatter plot grating meanburstratio FMI vs meanrate FMI:
 figsize = DEFAULTFIGURESIZE[0]*0.95, DEFAULTFIGURESIZE[1] # tweak to match others in 4S1

+ 1 - 302
fig4.py

@@ -7,11 +7,9 @@ grtmeasures = ['meanrate', 'meanrate',   'meanrate',   'blankmeanrate',
 mi = pd.MultiIndex.from_product([mvigrtmsustrs, STIMTYPES, [False, 'prestim', 'cond']],
                                 names=['msu', 'stimtype', 'blank'])
 fig4 = pd.DataFrame(index=mi, columns=['meanrate', 'meanburstratio'])
-## NOTE: make sure that rows are not overwritten due to temporarily enabling iteration
-## over st8 in any of the following loops!
 
 
-# strip, scatter and dumbbell plot movie vs grating FMI for applicable measures:
+# strip plot movie vs grating FMI for applicable measures:
 np.random.seed(0) # to get identical horizontal jitter in strip plots on every run
 figsize = DEFAULTFIGURESIZE
 linmin, linmax, linstep = -1, 1, 1
@@ -28,36 +26,14 @@ for mvimeasure, grtmeasure in zip(mvimeasures, grtmeasures):
     else:
         Axislabel = axislabel
     for st8 in ['none']:#ALLST8S:
-        #mvifmis, grtfmis = [], []
         mvimaxfmis, grtmaxfmis, grtmaxfmi2s = [], [], []
         exmplis, exmplmsustrs, normlis = [], [], []
         # collect all movie FMIs:
-        '''
-        for mseustr in mvimseustrs:
-            mvifmi = mviFMI[mvimeasure][mseustr, st8]
-            if pd.isna(mvifmi):
-                continue
-            mvifmis.append(mvifmi)
-        # collect all grating FMIs:
-        for mseustr in grtmseustrs:
-            grtfmi = grtFMI[grtmeasure][mseustr, st8]
-            if pd.isna(grtfmi):
-                continue
-            grtfmis.append(grtfmi)
-        '''
         # collect paired movie and grating max FMIs:
         keptmsui = 0
         for msustr in mvigrtmsustrs:
             mvimaxfmi = maxFMI[mvimeasure][msustr, st8, 'mvi']
             grtmaxfmi = maxFMI[grtmeasure][msustr, st8, 'grt']
-            '''
-            ## hack to replicate Steffen's code in R:
-            if mvimeasure == 'meanrate':
-                othermvimaxfmi = maxFMI['meanburstratio'][msustr, st8, 'mvi']
-                othergrtmaxfmi = maxFMI['meanburstratio'][msustr, st8, 'grt']
-                if pd.isna(othermvimaxfmi) or pd.isna(othergrtmaxfmi):
-                    continue
-            '''
             if pd.isna(mvimaxfmi) or pd.isna(grtmaxfmi): # missing one or both maxFMI values
                 continue
             if isblank: # add a second grtmaxfmi measure for blank cond
@@ -78,51 +54,17 @@ for mvimeasure, grtmeasure in zip(mvimeasures, grtmeasures):
                 else:
                     fig4.loc[msustr, 'mvi', False][measure] = mvimaxfmi # save
                     fig4.loc[msustr, 'grt', False][measure] = grtmaxfmi # save
-            #if grtmeasure == 'meanburstratio' and grtmaxfmi > 0.8:
-            #    print('OUTLIER: %s grtmaxfmi=%g' % (msustr, grtmaxfmi))
             if msustr in msu2exmpli:
                 exmplis.append(keptmsui)
                 exmplmsustrs.append(msustr)
             else:
                 normlis.append(keptmsui)
             keptmsui += 1 # manually increment
-        #mvifmis = np.asarray(mvifmis)
-        #grtfmis = np.asarray(grtfmis)
         mvimaxfmis = np.asarray(mvimaxfmis)
         grtmaxfmis = np.asarray(grtmaxfmis)
         grtmaxfmi2s = np.asarray(grtmaxfmi2s)
         assert len(mvimaxfmis) == len(grtmaxfmis)
         npairs = len(mvimaxfmis)
-        '''
-        ## strip plot movie and grating FMIs, for all mseus:
-        f, a = plt.subplots(figsize=figsize)
-        wintitle('FMI %s nat movie grating stripplot %s' % (mvimeasure, st8))
-        # plot y=0 line:
-        a.axhline(y=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
-        data = pd.DataFrame(data={'stimtype':['Movie']*len(mvifmis)+['Grating']*len(grtfmis),
-                                  'FMI':np.concatenate([mvifmis, grtfmis])})
-        sns.stripplot(x='stimtype', y='FMI', data=data, jitter=True, clip_on=False,
-                      marker='.', color='None', edgecolor=st82clr[st8], size=np.sqrt(50))
-        # exclude extreme +/- 1 FMI values from mean and ttest:
-        #mvifmisnon1 = mvifmis[abs(mvifmis) != 1]
-        #grtfmisnon1 = grtfmis[abs(grtfmis) != 1]
-        # plot mean with short horizontal lines:
-        #meanmvifmi, meangrtfmi = mvifmisnon1.mean(), grtfmisnon1.mean()
-        #a.plot([-0.25, 0.25], [meanmvifmi, meanmvifmi], '-', lw=2, c='red', zorder=np.inf)
-        #a.plot([0.75, 1.25], [meangrtfmi, meangrtfmi], '-', lw=2, c='red', zorder=np.inf)
-        a.set_xlabel('')
-        a.set_ylabel('%s FMI' % Axislabel)
-        a.set_ylim(-1, 1)
-        t, p = ttest_ind(mvifmisnon1, grtfmisnon1, equal_var=False) # non-paired t-test
-        a.add_artist(AnchoredText('p$=$%.2g' % p, loc='upper center', frameon=False))
-        # connect the dots:
-        x = np.array([[0]*npairs, [1]*npairs])
-        y = np.array([mvimaxfmis, grtmaxfmis])
-        a.plot(x, y, '-', c='k', alpha=0.2, lw=1)
-        # due to jitter, dots don't perfectly connect. Can get actual data using:
-        #mvixy, grtxy = a.collections[0].get_offsets(), a.collections[1].get_offsets()
-        # but some of those points aren't paired, which makes it complicated
-        '''
         ## strip plot paired movie and grating max FMIs:
         if mvimeasure == 'meanrate': # 2 column wide strip plot
             f, a = plt.subplots(figsize=(figsize[0]*1.35, figsize[1]*1.5))
@@ -176,10 +118,6 @@ for mvimeasure, grtmeasure in zip(mvimeasures, grtmeasures):
                     a.plot([1.75, 2.25], [meangrtmaxfmi2, meangrtmaxfmi2], '-', lw=2, c='red',
                            zorder=np.inf)
         a.set_ylabel('%s FMI' % Axislabel)
-        #if measure.startswith('meanrate'):
-        #    a.set_ylim(-0.6, 1)
-        #    a.set_yticks([-0.5, 0, 0.5, 1])
-        #else:
         a.set_ylim(-1, 1)
         a.set_yticks([-1, -0.5, 0, 0.5, 1])
         a.spines['bottom'].set_position(('outward', 5))
@@ -203,242 +141,3 @@ for mvimeasure, grtmeasure in zip(mvimeasures, grtmeasures):
         #a.plot(x[:, signchangeis], y[:, signchangeis], '-', c='k', alpha=1.0, lw=1)
         # due to jitter, dots don't perfectly connect. Can get actual data using:
         #mvixy, grtxy = a.collections[0].get_offsets(), a.collections[1].get_offsets()
-        ## scatter plot movie vs. grating max FMIs:
-        '''
-        f, a = plt.subplots(figsize=figsize)
-        wintitle('maxFMI %s movie grating scatter' % mvimeasure)
-        # plot x=0 and y=0 lines:
-        a.axhline(y=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
-        a.axvline(x=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
-        # plot y=x line:
-        xyline = [linmin, linmax], [linmin, linmax]
-        a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
-        # plot normal (non-example) points:
-        a.scatter(grtmaxfmis[normlis], mvimaxfmis[normlis], clip_on=False,
-                  marker='.', c='None', edgecolor=st82clr[st8], s=DEFSZ)
-        # plot example points, one at a time:
-        for exmpli, msustr in zip(exmplis, exmplmsustrs):
-            marker = exmpli2mrk[msu2exmpli[msustr]]
-            c = exmpli2clr[msu2exmpli[msustr]]
-            sz = exmpli2sz[msu2exmpli[msustr]]
-            lw = exmpli2lw[msu2exmpli[msustr]]
-            a.scatter(grtmaxfmis[exmpli], mvimaxfmis[exmpli], marker=marker, c=c, s=sz, lw=lw)
-        # exclude extreme +/- 1 FMI values from mean and ttest:
-        non1is = (abs(grtmaxfmis) != 1) & (abs(mvimaxfmis) != 1)
-        grtmaxfmisnon1 = grtmaxfmis[non1is]
-        mvimaxfmisnon1 = mvimaxfmis[non1is]
-        # plot mean:
-        a.scatter(np.mean(grtmaxfmisnon1), np.mean(mvimaxfmisnon1),
-                  c='red', edgecolor='red', s=50, marker='^')
-        a.set_xlabel('Grating %s FMI' % axislabel)
-        a.set_ylabel('Movie %s FMI' % axislabel)
-        a.set_xlim(linmin, linmax)
-        a.set_ylim(linmin, linmax)
-        a.set_xticks(ticks)
-        a.set_yticks(a.get_xticks()) # make log scale y ticks the same as x ticks
-        a.minorticks_off()
-        a.set_aspect('equal')
-        a.spines['left'].set_position(('outward', 4))
-        a.spines['bottom'].set_position(('outward', 4))
-        t, p = ttest_rel(grtmaxfmisnon1, mvimaxfmisnon1) # paired t-test
-        a.add_artist(AnchoredText('p$=$%.2g' % p, loc='lower right', frameon=False))
-        '''
-        ## dumbbell plot:
-        '''
-        f, a = plt.subplots(figsize=DUMBBELLFIGSIZE)
-        wintitle('maxFMI %s movie grating dumbbell' % mvimeasure)
-        # plot x=0 line:
-        a.axvline(x=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
-        # plot horizontal lines:
-        mvimaxfmisortis = mvimaxfmis.argsort() # y vals for hlines
-        npairs = len(mvimaxfmis)
-        ys = np.arange(npairs)
-        xmvis, xgrts = mvimaxfmis[mvimaxfmisortis], grtmaxfmis[mvimaxfmisortis]
-        a.hlines(ys, xmvis, xgrts, zorder=-100)
-        sz = 100
-        # plot grt points:
-        a.scatter(xgrts, ys, clip_on=False, marker='.', c='white',
-                  edgecolor=st82clr[st8], s=sz)
-        # plot mvi points:
-        a.scatter(xmvis, ys, clip_on=False, marker='.', c=st82clr[st8], s=sz)
-        a.set_xlabel('%s FMI' % Axislabel)
-        a.set_ylabel('Neurons')
-        a.set_xlim(linmin, linmax)
-        a.set_ylim(-1, npairs)
-        a.set_xticks(ticks)
-        a.set_yticks([])
-        a.spines['left'].set_position(('outward', 4))
-        a.spines['left'].set_visible(False)
-        '''
-
-'''
-# plot CDFs of spatial suppression index for blankscreen movies/gratings wrt full screen
-# movies/gratings, for meanrate and meanburstratio, for control and opto conditions:
-figsize = DEFAULTFIGURESIZE
-stimtype2resp = {'mvi':mviresp, 'grt':grtresp}
-for stimtype in ['mvi', 'grt']:
-    resp = stimtype2resp[stimtype]
-    if stimtype == 'mvi':
-        resp = resp.xs('nat', level='kind') # dereference movie 'kind' index level
-    for measure in ['meanrate', 'meanburstratio']:
-        f, a = plt.subplots(figsize=figsize)
-        wintitle('SSI %s %s blank vs fullscreen CDF' % (stimtype, measure))
-        opto2SSI = {}
-        for opto in OPTOS:
-            rows = resp.xs(['none', opto], level=['st8', 'opto']) # only mseu index left
-            fullscreen = rows[measure]
-            blank = rows['blank'+measure]
-            SSI = (blank - fullscreen) / (blank + fullscreen)
-            SSI = np.float64(SSI[SSI.notna()].values) # pull float array out of object Series
-            opto2SSI[opto] = SSI
-            SSIbins = list(np.unique(SSI)) + [10]
-            fb = opto2fb[opto].capitalize()
-            a.hist(SSI, bins=SSIbins, density=True, histtype='step',
-                   cumulative=True, clip_on=True, lw=1.5,
-                   label=fb, color=opto2clr[opto], alpha=1)
-        a.set_xlabel('Spatial suppression index')
-        a.set_ylabel('Cumulative probability')
-        #a.legend(frameon=False)
-        a.set_xlim(-1, 1)
-        a.set_ylim(0, 1.01)
-        a.set_yticks([0, 0.5, 1])
-        #a.spines['left'].set_position(('outward', 4))
-        _, KS_p = ks_2samp(opto2SSI[False], opto2SSI[True])
-        txt = '$\mathregular{p_{KS}=%.2g}$' % KS_p # prob that distribs are the same
-        a.add_artist(AnchoredText(txt, loc='lower right', frameon=False))
-
-
-# stripplot spatial suppression index for blankscreen movies/gratings wrt full screen
-# movies/gratings, for meanrate and meanburstratio, for both opto conditions, for all mseu:
-figsize = DEFAULTFIGURESIZE
-stimtype2resp = {'mvi':mviresp, 'grt':grtresp}
-for measure in ['meanrate', 'meanburstratio']:
-    opto2SSI = {}
-    for stimtype in ['mvi', 'grt']:
-        resp = stimtype2resp[stimtype]
-        if stimtype == 'mvi':
-            resp = resp.xs('nat', level='kind') # dereference movie 'kind' index level
-        for opto in OPTOS:
-            rows = resp.xs(['none', opto], level=['st8', 'opto']) # only mseu index left
-            fullscreen = rows[measure]
-            blank = rows['blank'+measure]
-            SSI = (blank - fullscreen) / (blank + fullscreen)
-            #SSI = (blank - fullscreen) / blank
-            SSI = np.float64(SSI[SSI.notna()].values) # pull float array out of object Series
-            opto2SSI[opto] = SSI
-        # stripplot:
-        SSIons, SSIoffs = opto2SSI[True], opto2SSI[False]
-        f, a = plt.subplots(figsize=figsize)
-        wintitle('SSI %s %s blank vs fullscreen stripplot' % (stimtype, measure))
-        data = pd.DataFrame.from_dict({'Feedback':SSIoffs, 'Suppression':SSIons},
-                                      orient='index').transpose()
-        sns.stripplot(ax=a, data=data, clip_on=False, marker='.',
-                      color='None', edgecolor='black', size=np.sqrt(50))
-        # exclude extreme +/- 1 SSI values from mean:
-        SSIonsnon1 = SSIons[abs(SSIons) != 1]
-        SSIoffsnon1 = SSIoffs[abs(SSIoffs) != 1]
-        # plot mean with short horizontal lines:
-        meanSSIon, meanSSIoff = SSIonsnon1.mean(), SSIoffsnon1.mean()
-        a.plot([-0.25, 0.25], [meanSSIoff, meanSSIoff], '-', lw=2, c='red', zorder=np.inf)
-        a.plot([0.75, 1.25], [meanSSIon, meanSSIon], '-', lw=2, c='red', zorder=np.inf)
-        a.set_ylabel('Spatial suppression index')
-        a.set_ylim(-1, 1)
-        a.spines['bottom'].set_position(('outward', 5))
-        a.spines['bottom'].set_visible(False)
-        a.tick_params(bottom=False)
-
-
-# stripplot FMI spatial suppression index for blankscreen movies/gratings wrt full screen
-# movies/gratings, for meanrate and meanburstratio, for all mseu:
-figsize = DEFAULTFIGURESIZE
-stimtype2FMI = {'mvi':mviFMI, 'grt':grtFMI}
-for measure in ['meanrate', 'meanburstratio']:
-    stimtype2SSI = {}
-    for stimtype in ['mvi', 'grt']:
-        FMI = stimtype2FMI[stimtype]
-        rows = FMI.xs('none', level='st8') # only mseu index left
-        fullscreen = rows[measure]
-        blank = rows['blank'+measure]
-        SSI = (blank - fullscreen) / (abs(blank) + abs(fullscreen))
-        #SSI = (blank - fullscreen) / blank
-        SSI = np.float64(SSI[SSI.notna()].values) # pull float array out of object Series
-        stimtype2SSI[stimtype] = SSI
-        # plot CDF:
-        f, a = plt.subplots(figsize=figsize)
-        wintitle('SSI %s %s FMI blank vs fullscreen' % (stimtype, measure))
-        SSIbins = list(np.unique(SSI)) + [10]
-        a.hist(SSI, bins=SSIbins, density=True, histtype='step',
-               cumulative=True, clip_on=True, lw=1.5, color='black', alpha=1)
-        a.set_xlabel('FMI spatial suppression index')
-        a.set_ylabel('Cumulative probability')
-        #a.legend(frameon=False)
-        a.set_xlim(-1, 1)
-        a.set_ylim(0, 1.01)
-        a.set_yticks([0, 0.5, 1])
-        #a.spines['left'].set_position(('outward', 4))
-    # stripplot:
-    mviSSI, grtSSI = stimtype2SSI['mvi'], stimtype2SSI['grt']
-    f, a = plt.subplots(figsize=figsize)
-    wintitle('SSI %s FMI blank vs fullscreen stripplot' % measure)
-    data = pd.DataFrame.from_dict({'Movie':mviSSI, 'Grating':grtSSI},
-                                  orient='index').transpose()
-    sns.stripplot(ax=a, data=data, clip_on=False, marker='.',
-                  color='None', edgecolor='black', size=np.sqrt(50))
-    # exclude extreme +/- 1 SSI values from mean:
-    mviSSInon1 = mviSSI[abs(mviSSI) != 1]
-    grtSSInon1 = grtSSI[abs(grtSSI) != 1]
-    # plot mean with short horizontal lines:
-    meanmvifmi, meangrtfmi = mviSSInon1.mean(), grtSSInon1.mean()
-    a.plot([-0.25, 0.25], [meanmvifmi, meanmvifmi], '-', lw=2, c='red', zorder=np.inf)
-    a.plot([0.75, 1.25], [meangrtfmi, meangrtfmi], '-', lw=2, c='red', zorder=np.inf)
-    a.set_ylabel('FMI spatial suppression index')
-    a.set_ylim(-1, 1)
-    a.spines['bottom'].set_position(('outward', 5))
-    a.spines['bottom'].set_visible(False)
-    a.tick_params(bottom=False)
-
-
-# dumbbell plot FMI spatial suppression index for blankscreen movies/gratings wrt full screen
-# movies/gratings, for meanrate and meanburstratio, for all mseu:
-figsize = DEFAULTFIGURESIZE
-for measure in ['meanrate', 'meanburstratio']:
-    stimtype2SSI = {}
-    for stimtype in ['mvi', 'grt']:
-        rows = maxFMI.xs(['none', stimtype], level=['st8', 'stimtype']) # only msu index left
-        fullscreen = rows[measure]
-        blank = rows['blank'+measure]
-        SSI = (blank - fullscreen) / (abs(blank) + abs(fullscreen))
-        #SSI = (blank - fullscreen) / blank
-        SSI = np.float64(SSI.values) # pull float array out of object Series
-        stimtype2SSI[stimtype] = SSI
-
-    mviSSI, grtSSI = stimtype2SSI['mvi'], stimtype2SSI['grt']
-    keepis = pd.notna(mviSSI) & pd.notna(grtSSI) # remove any unit that has NaN in either stim
-    mviSSI, grtSSI = mviSSI[keepis], grtSSI[keepis]
-    f, a = plt.subplots(figsize=figsize)
-    wintitle('SSI %s maxFMI blank vs fullscreen stripplot dumbbell' % measure)
-    # plot x=0 line:
-    a.axvline(x=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
-    # plot horizontal lines:
-    mviSSIsortis = mviSSI.argsort() # y vals for hlines
-    npairs = len(mviSSI)
-    ys = np.arange(npairs)
-    xmvis, xgrts = mviSSI[mviSSIsortis], grtSSI[mviSSIsortis]
-    a.hlines(ys, xmvis, xgrts, zorder=-100)
-    sz = 100
-    # plot grt points:
-    a.scatter(xgrts, ys, marker='.', c='white', edgecolor=st82clr[st8], s=sz)
-    # plot mvi points:
-    a.scatter(xmvis, ys, marker='.', c=st82clr[st8], s=sz)
-    axislabel = measure2axislabel[measure]
-    if axislabel[0].islower():
-        Axislabel = axislabel.capitalize()
-    else:
-        Axislabel = axislabel
-    a.set_xlabel('%s FMI SSI' % Axislabel)
-    a.set_ylabel('Neurons')
-    #a.set_xlim(linmin-eps, linmax+eps)
-    a.set_ylim(-1, npairs)
-    #a.set_xticks(ticks)
-    a.set_yticks([])
-'''

+ 0 - 8
fig4S1.py

@@ -154,10 +154,6 @@ for kind in ['nat']:#MVIKINDS:
             keptmseui += 1 # manually increment
         rons = np.asarray(rons)
         roffs = np.asarray(roffs)
-        # replace off-scale low values with log0min, so the points remain visible:
-        #pltrons, pltroffs = rons.copy(), roffs.copy()
-        #pltrons[pltrons <= 10**logmin] = 10**log0min
-        #pltroffs[pltroffs <= 10**logmin] = 10**log0min
         # plot y=x line:
         xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
         a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
@@ -298,10 +294,6 @@ for st8 in ['none']:#ALLST8S:
             keptmseui += 1 # manually increment
         rons = np.asarray(rons)
         roffs = np.asarray(roffs)
-        # replace off-scale low values with log0min, so the points remain visible:
-        #pltrons, pltroffs = rons.copy(), roffs.copy()
-        #pltrons[pltrons <= 10**logmin] = 10**log0min
-        #pltroffs[pltroffs <= 10**logmin] = 10**log0min
         # plot y=x line:
         xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
         a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)

+ 6 - 12
fig5.py

@@ -4,13 +4,11 @@ mi = pd.MultiIndex.from_product([mvimseustrs, ST8S, OPTOS], names=['mseu', 'st8'
 columns = ['trialis', 'rates', 'burstratios', 'meanrate', 'meanburstratio',
            'spars', 'rel', 'meanpkw', 'snr']
 fig5 = pd.DataFrame(index=mi, columns=columns) # excludes kind and opto
-## NOTE: make sure that rows are not overwritten due to temporarily enabling iteration
-## over kind in any of the following loops!
-
 
 # tweak exactly how much horizontal space to use for raster plots and PSTHs:
 PLOTOFFSETS = -0.45, 0.45
 
+
 # plot movie rasters for each unit by kind, state and opto condition:
 # copied and slightly modified from fig1.py:
 showbursts = True
@@ -167,10 +165,6 @@ for kind in ['nat']:#MVIKINDS:
         rruns = np.asarray(rruns)
         rsits = np.asarray(rsits)
         sgnfs = np.asarray(sgnfs)
-        # replace off-scale low values with log0min, so the points remain visible:
-        #pltrruns, pltrsits = rruns.copy(), rsits.copy()
-        #pltrruns[pltrruns <= 10**logmin] = 10**log0min
-        #pltrsits[pltrsits <= 10**logmin] = 10**log0min
         # plot y=x line:
         xyline = [10**logmin, 10**logmax], [10**logmin, 10**logmax]
         a.plot(xyline[0], xyline[1], '--', color='gray', zorder=-1)
@@ -577,15 +571,15 @@ for kind in ['nat']:#MVIKINDS:
                           marker=marker, c=c, s=sz, lw=lw)
             # get fname of appropriate LMM .cvs file:
             if measure == 'meanburstratio': # use Steffen's LMM linregress fit, for fig5S1
-                fname = os.path.join('stats', 'figure_5S1c_coefs.csv')
+                fname = os.path.join('stats', 'figure_5_S1c_coefs.csv')
             elif measure == 'spars': # use Steffen's LMM linregress fit, for fig5S1
-                fname = os.path.join('stats', 'figure_5S1d_coefs.csv')
+                fname = os.path.join('stats', 'figure_5_S1d_coefs.csv')
             elif measure == 'rel': # use Steffen's LMM linregress fit, for fig5S1
-                fname = os.path.join('stats', 'figure_5S1e_coefs.csv')
+                fname = os.path.join('stats', 'figure_5_S1e_coefs.csv')
             elif measure == 'snr': # use Steffen's LMM linregress fit, for fig5S1
-                fname = os.path.join('stats', 'figure_5S1f_coefs.csv')
+                fname = os.path.join('stats', 'figure_5_S1f_coefs.csv')
             elif measure == 'meanpkw': # use Steffen's LMM linregress fit, for fig5S1
-                fname = os.path.join('stats', 'figure_5S1g_coefs.csv')
+                fname = os.path.join('stats', 'figure_5_S1g_coefs.csv')
             else:
                 print("WARNING: No LMM stats for %s measure" % measure)
                 fname = None

+ 0 - 2
fig5S2.py

@@ -35,7 +35,6 @@ for index, row in ipos_opto.iterrows():
         f = scipy.interpolate.interp1d(run_ts[i], speed[i], axis=-1, kind='linear',
                                        fill_value='extrapolate')
         rs_speed = f(arearate_ts) # runspeed signal resampled to arearate_ts
-        ## TODO: maybe apply continuous_runs() here to further smooth out running periods
         # mask out arearate during running periods:
         runmask = rs_speed > 0.25 # cm/s
         arearate[i, runmask] = np.nan # these points shouldn't influence the percentile calc
@@ -336,7 +335,6 @@ for opto in OPTOS:
     f, a = plt.subplots(figsize=figsize)
     wintitle('top bottom qrt scatter meanrate opto=%s' % opto)
     # replace any FRs < MINFRVAL with MINFRVAL, and any > MAXFRVAL with MAXFRVAL, for plotting:
-    #logmin, logmax = -0.5, 1.7
     logmin, logmax = -0.9, 1.9
     logticks = np.array([0, 1])
     MINFRVAL, MAXFRVAL = 10**logmin, 10**logmax

+ 1 - 123
fig6.py

@@ -7,8 +7,6 @@ mi = pd.MultiIndex.from_product([allmseustrs, STIMTYPES, modis],
 columns = ['meanrate', 'meanburstratio', 'spars', 'rel', 'meanpkw', 'snr']
 fig6 = pd.DataFrame(index=mi, columns=columns)
 
-# redefine these here, in case saved data was loaded from pickles and these definitions
-# later in calc.py were therefore skipped:
 stimtype2FMI = {'mvi': mviFMI, 'grt': grtFMI}
 stimtype2RMI = {'mvi': mviRMI, 'grt': grtRMI}
 
@@ -16,7 +14,7 @@ loss = 'soft_l1' # loss function, less sensitive to outliers than ordinary least
 f_scale = 0.25 # residual value at which to start treating points as outliers
 
 
-# scatter and dumbbell plot run modulation index for feedback vs. suppression trials,
+# scatter plot run modulation index for feedback vs. suppression trials,
 # for movies and gratings, for all measures:
 figsize = DEFAULTFIGURESIZE
 linmin, linmax, linstep = -1, 1, 1
@@ -98,23 +96,6 @@ for measure in modmeasuresnoblank:
                 x = np.array([np.nanmin(supprmis), np.nanmax(supprmis)])
                 y = mm * x + b
                 a.plot(x, y, '-', color='red') # plot linregress fit
-        '''
-        else:
-            # calc/plot linregress line:
-            # exclude extreme +/- 1 FMI and RMI values from linregress:
-            non1is = (abs(supprmis) != 1) & (abs(fbrmis) != 1)
-            supprmisnon1 = supprmis[non1is]
-            fbrmisnon1 = fbrmis[non1is]
-            # robust least squares:
-            result = least_squares(linear_loss, (0, 0), args=(supprmisnon1, fbrmisnon1),
-                                   loss=loss, f_scale=f_scale)
-            mm, b = result['x']
-            residuals = result['fun']
-            #rsq = residual_rsquared(fbrmisnon1, residuals)
-            #mm, b, rr, p, stderr = linregress(rmisnon1, fmisnon1)
-            #rsq = rr * rr
-            x = np.array([supprmisnon1.min(), supprmisnon1.max()])
-        '''
         #a.set_title('%s RMI' % axislabel) # typeface too big by default, decreases panel size
         a.set_xlabel('Suppression')
         a.set_ylabel('Feedback')
@@ -126,34 +107,6 @@ for measure in modmeasuresnoblank:
         a.set_aspect('equal')
         a.spines['left'].set_position(('outward', 4))
         a.spines['bottom'].set_position(('outward', 4))
-        #txt = ('$\mathregular{R^{2}=%.2g}$' % rsq)
-        #.add_artist(AnchoredText(txt, loc='lower right', frameon=False))
-        #txt = ('$\mathregular{p=%.2g}$\n'
-        #       '$\mathregular{R^{2}=%.2g}$' % (p, rsq)) # prob that slope is 0
-        #a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))
-        ## dumbbell plot feedback vs suppression RMI
-        '''
-        f, a = plt.subplots(figsize=DUMBBELLFIGSIZE)
-        wintitle('fbsupp RMI dumbbell %s %s' % (measure, stimtypelabel))
-        # plot x=0 line:
-        a.axvline(x=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
-        # plot horizontal lines:
-        sortis = fbrmis.argsort() # y vals for hlines
-        npairs = len(fbrmis)
-        ys = np.arange(npairs)
-        srtdfbrmis, srtdsupprmis = fbrmis[sortis], supprmis[sortis]
-        a.hlines(ys, srtdfbrmis, srtdsupprmis, zorder=-100)
-        sz = 50
-        # plot supprmis points:
-        a.scatter(srtdsupprmis, ys, marker='.', c='white', edgecolor='black', s=sz)
-        # plot fbrmis points:
-        a.scatter(srtdfbrmis, ys, marker='.', c='black', s=sz)
-        a.set_xlabel('%s RMI' % axislabel)
-        a.set_ylabel('Neurons')
-        a.set_xlim(linmin, linmax)
-        a.set_xticks(ticks)
-        a.set_yticks([])
-        '''
 
 
 # scatter plot feedback modulation index for run vs. sit trials,
@@ -244,23 +197,6 @@ for measure in modmeasuresnoblank:
                 x = np.array([np.nanmin(sitfmis), np.nanmax(sitfmis)])
                 y = mm * x + b
                 a.plot(x, y, '-', color='red') # plot linregress fit
-        '''
-        else:
-            # calc/plot linregress line:
-            # exclude extreme +/- 1 FMI values from linregress:
-            non1is = (abs(sitfmis) != 1) & (abs(runfmis) != 1)
-            runfmisnon1 = runfmis[non1is]
-            sitfmisnon1 = sitfmis[non1is]
-            # robust least squares:
-            result = least_squares(linear_loss, (0, 0), args=(sitfmisnon1, runfmisnon1),
-                                   loss=loss, f_scale=f_scale)
-            mm, b = result['x']
-            residuals = result['fun']
-            #rsq = residual_rsquared(runfmisnon1, residuals)
-            #mm, b, rr, p, stderr = linregress(sitfmisnon1, runfmisnon1)
-            #rsq = rr * rr
-            x = np.array([sitfmisnon1.min(), sitfmisnon1.max()])
-        '''
         #a.set_title('%s FMI' % axislabel) # typeface too big by default, decreases panel size
         a.set_xlabel('Sit')
         a.set_ylabel('Run')
@@ -272,34 +208,6 @@ for measure in modmeasuresnoblank:
         a.set_aspect('equal')
         a.spines['left'].set_position(('outward', 4))
         a.spines['bottom'].set_position(('outward', 4))
-        #txt = ('$\mathregular{R^{2}=%.2g}$' % rsq)
-        #a.add_artist(AnchoredText(txt, loc='lower right', frameon=False))
-        #txt = ('$\mathregular{p=%.2g}$\n'
-        #       '$\mathregular{R^{2}=%.2g}$' % (p, rsq)) # prob that slope is 0
-        #a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))
-        ## dumbbell plot run vs sit FMI
-        '''
-        f, a = plt.subplots(figsize=DUMBBELLFIGSIZE)
-        wintitle('runsit FMI dumbbell %s %s' % (measure, stimtypelabel))
-        # plot x=0 line:
-        a.axvline(x=0, ls='--', marker='', color='lightgray', zorder=-np.inf)
-        # plot horizontal lines:
-        sortis = runfmis.argsort() # y vals for hlines
-        npairs = len(runfmis)
-        ys = np.arange(npairs)
-        srtdrunfmis, srtdsitfmis = runfmis[sortis], sitfmis[sortis]
-        a.hlines(ys, srtdrunfmis, srtdsitfmis, zorder=-100)
-        sz = 50
-        # plot sitfmis points:
-        a.scatter(srtdsitfmis, ys, marker='.', c='white', edgecolor='black', s=sz)
-        # plot runfmis points:
-        a.scatter(srtdrunfmis, ys, marker='.', c='black', s=sz)
-        a.set_xlabel('%s FMI' % short2longaxislabel.get(axislabel, axislabel))
-        a.set_ylabel('Neurons')
-        a.set_xlim(linmin, linmax)
-        a.set_xticks(ticks)
-        a.set_yticks([])
-        '''
 
 
 # scatter plot feedback modulation index vs. run modulation index,
@@ -321,14 +229,6 @@ for measure in modmeasuresnoblank:
         keptmseui = 0 # manually init and increment instead of using enumerate()
         for mseustr in mseustrs:
             fmi = FMI[measure][mseustr, 'none'] # ignore run condition for FMI
-            """
-            TODO: plot RMI calc'd across both opto conditions, but before doing that, you
-            have to normalize by the relative mean measure (say rate) of the two opto
-            conditions, and weight them appropriately, says Jackson Smith from Oxford. Should
-            really do the equivalent for FMI, especially cuz there are not equal numbers of
-            trials for run vs sit, but the analysis that plots FMI as a function of run vs sit
-            trials more or less shows that you can pool across run condition to calculate FMI
-            """
             rmi = RMI[measure][mseustr, False] # use only feedback condition for RMI
             if np.isnan(fmi) or np.isnan(rmi): # missing one or both mod values
                 continue
@@ -399,23 +299,6 @@ for measure in modmeasuresnoblank:
                 x = np.array([np.nanmin(rmis), np.nanmax(rmis)])
                 y = mm * x + b
                 a.plot(x, y, '-', color='red') # plot linregress fit
-        '''
-        else:
-            # calc/plot linregress line:
-            # exclude extreme +/- 1 FMI and RMI values from linregress:
-            non1is = (abs(rmis) != 1) & (abs(fmis) != 1)
-            rmisnon1 = rmis[non1is]
-            fmisnon1 = fmis[non1is]
-            # robust least squares:
-            result = least_squares(linear_loss, (0, 0), args=(rmisnon1, fmisnon1),
-                                   loss=loss, f_scale=f_scale)
-            mm, b = result['x']
-            residuals = result['fun']
-            rsq = residual_rsquared(fmisnon1, residuals)
-            mm, b, rr, p, stderr = linregress(rmisnon1, fmisnon1)
-            rsq = rr * rr
-            x = np.array([rmisnon1.min(), rmisnon1.max()])
-        '''
         #a.set_title('%s' % axislabel) # typeface too big by default, decreases panel size
         a.set_xlabel('RMI')
         a.set_ylabel('FMI')
@@ -427,8 +310,3 @@ for measure in modmeasuresnoblank:
         a.set_aspect('equal')
         a.spines['left'].set_position(('outward', 4))
         a.spines['bottom'].set_position(('outward', 4))
-        #txt = ('$\mathregular{R^{2}=%.2g}$' % rsq)
-        #a.add_artist(AnchoredText(txt, loc='lower right', frameon=False))
-        #txt = ('$\mathregular{p=%.2g}$\n'
-        #       '$\mathregular{R^{2}=%.2g}$' % (p, rsq)) # prob that slope is 0
-        #a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))

+ 2 - 28
ipos.py

@@ -45,13 +45,11 @@ for stimtype, stimtypei in zip(STIMTYPESPLUSALL, STIMTYPESPLUSALLI):
     a.set_yticks([0, 0.5, 1])
     #a.legend(frameon=False)
     _, KS_p = ks_2samp(stdctrl, stdopto)
-    #_, MW_p = scipy.stats.mannwhitneyu(stdrun, stdsit)
     txt = '$\mathregular{p_{KS}=%.2g}$' % KS_p # prob that distribs are the same
     a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))
 
 
 ## scatter plot movie reliability FMI vs eye position stdev FMI:
-## TODO: this correlation is done with FMI calculated from run trials. Is this optimal?
 # run reliability FMI for each movie mseu:
 mseus, relfmis = zip(*mviFMI['rel'][:, 'run'].items()) # unpack indices and relfmi values
 mseus, relfmis = np.array(mseus), np.array(relfmis)
@@ -76,7 +74,7 @@ wintitle('rel FMI vs ipos stdev FMI')
 a.scatter(iposfmis[valid], relfmis[valid], clip_on=False,
           marker='.', c='None', edgecolor='black', s=DEFSZ)
 # plot regression:
-fname = os.path.join('stats', 'figure_1S2i_coefs.csv')
+fname = os.path.join('stats', 'figure_1_S2i_coefs.csv')
 try:
     df = pd.read_csv(fname)
     foundregression = True
@@ -99,10 +97,6 @@ a.set_xticks([xmin, 0, xmax])
 a.set_yticks([ymin, 0, ymax])
 a.spines['left'].set_position(('outward', 4))
 a.spines['bottom'].set_position(('outward', 4))
-#txt = ('$\mathregular{p=%.2g}$\n' # prob that slope is 0
-#       '$\mathregular{R^{2}=%.2g}$\n'
-#       '$\mathregular{slope=%.2g}$' % (p, rsq, mm))
-#a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))
 # save to iposmi dataframe:
 for mseu, iposfmi, relfmi in zip(mseus, iposfmis, relfmis):
     iposmi.loc[mseu, 'stimtype'] = 'mvi'
@@ -153,21 +147,6 @@ for stimtype, stimtypei in zip(STIMTYPESPLUSALL, STIMTYPESPLUSALLI):
     txt = '$\mathregular{p_{KS}=%.2g}$' % KS_p # prob that distribs are the same
     a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))
 
-'''
-## plot within trial stdev PDFs for run and sit trials:
-f, a = plt.subplots(figsize=DEFAULTFIGURESIZE)
-wintitle('ipos_st8 stdev within trial PDF')
-stdrun = ipos_st8['std_xpos_within'][:, 'run']
-stdsit = ipos_st8['std_xpos_within'][:, 'sit']
-a.hist(stdrun, bins=iposedges, histtype='step', lw=1.5,
-       color=st82clr['run'], alpha=1, label='Run')
-a.hist(stdsit, bins=iposedges, histtype='step', lw=1.5,
-       color=st82clr['sit'], alpha=1, label='Sit')
-a.set_xlabel('Eye position $\sigma$ ($\mathregular{\degree}$)')
-a.set_ylabel('Experiment count')
-#a.legend(frameon=False)
-print('st8 std_within:', wilcoxon(stdrun, stdsit))
-'''
 
 ## scatter plot movie reliability RMI vs eye position stdev FMI:
 # feedback reliability RMI for each movie mseu:
@@ -194,7 +173,7 @@ wintitle('rel RMI vs ipos stdev RMI')
 a.scatter(iposrmis[valid], relrmis[valid], clip_on=False,
           marker='.', c='None', edgecolor='black', s=DEFSZ)
 # plot regression:
-fname = os.path.join('stats', 'figure_5S1i_coefs.csv')
+fname = os.path.join('stats', 'figure_5_S1i_coefs.csv')
 try:
     df = pd.read_csv(fname)
     foundregression = True
@@ -216,13 +195,8 @@ a.set_xticks([0, 0.5])
 a.set_yticks([ymin, 0, ymax])
 a.spines['left'].set_position(('outward', 4))
 a.spines['bottom'].set_position(('outward', 4))
-#txt = ('$\mathregular{p=%.2g}$\n' # prob that slope is 0
-#       '$\mathregular{R^{2}=%.2g}$\n'
-#       '$\mathregular{slope=%.2g}$' % (p, rsq, mm))
-#a.add_artist(AnchoredText(txt, loc='upper left', frameon=False))
 # save to iposmi dataframe:
 for mseu, iposrmi, relrmi in zip(mseus, iposrmis, relrmis):
     iposmi.loc[mseu, 'stimtype'] = 'mvi' # automatically adds a new mseu row if needed
     iposmi.loc[mseu, 'iposrmi'] = iposrmi
     iposmi.loc[mseu, 'relrmi'] = relrmi
-

+ 10 - 11
main.py

@@ -102,20 +102,20 @@ fig5exmpli = 1 # fig5 uses example neuron 1
 
 # designate movie example units:
 pvmvimseu2exmpli = {'PVCre_2018_0003_s03_e03_u51':1, # example neuron 1
-                    'PVCre_2017_0008_s12_e06_u56':2, # example neuron 2
-                    }
-ntsrmvimseu2exmpli = {'Ntsr1Cre_2019_0008_s03_e04_u18':1,
-                     }
-exptype2mviexmplis = {'pvmvis':pvmvimseu2exmpli, 'ntsrmvis':ntsrmvimseu2exmpli}
+                    'PVCre_2017_0008_s12_e06_u56':2} # example neuron 2
+ntsrmvimseu2exmpli = {'Ntsr1Cre_2019_0008_s03_e04_u18':1}
+negntsrmvimseu2exmpli = {}
+exptype2mviexmplis = {'pvmvis':pvmvimseu2exmpli, 'ntsrmvis':ntsrmvimseu2exmpli,
+                      'negntsrmvis':negntsrmvimseu2exmpli}
 mvimseu2exmpli = exptype2mviexmplis[EXPTYPE]
 
 # designate grating example units:
 pvgrtmseu2exmpli = {'PVCre_2018_0003_s03_e02_u51':1,  # example neuron 1
-                    'PVCre_2018_0003_s05_e02_u164':3, # example neuron 3
-                   }
-ntsrgrtmseu2exmpli = {'Ntsr1Cre_2019_0008_s06_e02_u46':1,
-                     }
-exptype2grtexmplis = {'pvmvis':pvgrtmseu2exmpli, 'ntsrmvis':ntsrgrtmseu2exmpli}
+                    'PVCre_2018_0003_s05_e02_u164':3} # example neuron 3
+ntsrgrtmseu2exmpli = {'Ntsr1Cre_2019_0008_s06_e02_u46':1}
+ntsrgrtmseu2exmpli = {}
+exptype2grtexmplis = {'pvmvis':pvgrtmseu2exmpli, 'ntsrmvis':ntsrgrtmseu2exmpli,
+                      'negntsrmvis':ntsrgrtmseu2exmpli}
 grtmseu2exmpli = exptype2grtexmplis[EXPTYPE]
 
 mvimsu2exmpli = { mseustr2msustr(k):v for k, v in mvimseu2exmpli.items() }
@@ -160,7 +160,6 @@ 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

+ 0 - 19
util.py

@@ -571,25 +571,6 @@ def raster2freqcomp(raster, dt, f, mean='scalar'):
         # discard rare spikes in raster that for some reason (screen vsyncs?) fall outside
         # the expected trial duration:
         spikes = spikes[spikes <= dt]
-        ## TODO: getHarmResps.m only ever used an integer number of cycles, but why? Can't we
-        ## instead use the exact fractional number of cycles that the spikes occupy?
-        '''
-        ## UNTESTED:
-        if f != 0:
-            period = 1 / f
-            nperiods = int(np.floor(dt/period))
-            if nperiods == 0:
-                dt = period
-            else:
-                dt = nperiods*period
-            n = dt / f # number of cycles in this trial's spike train
-            nint = int(np.floor(n)) # number of full cycles before the end
-            nfrac = n % 1 # fractional number of cycles at the end
-            # need to split off spikes from last fractional cycle:
-            dtint, dtfrac = nint*f, nfrac*f # (s)
-            fracspikei = spikes.searchsorted(dtint)
-            fullspikes, fracspikes = spikes[:fracspikei], spikes[frackspikei:]
-        '''
         omega = 2 * np.pi * f # angular frequency (rad/s)
         res[triali] = (np.cos(omega*spikes)).sum() / dt
         ims[triali] = (np.sin(omega*spikes)).sum() / dt