Olfactory bulb microcircuits model with dual-layer inhibition (Gilra & Bhalla 2015)

 Download zip file 
Help downloading and running models
Accession:153574
A detailed network model of the dual-layer dendro-dendritic inhibitory microcircuits in the rat olfactory bulb comprising compartmental mitral, granule and PG cells developed by Aditya Gilra, Upinder S. Bhalla (2015). All cell morphologies and network connections are in NeuroML v1.8.0. PG and granule cell channels and synapses are also in NeuroML v1.8.0. Mitral cell channels and synapses are in native python.
Reference:
1 . Gilra A, Bhalla US (2015) Bulbar microcircuit model predicts connectivity and roles of interneurons in odor coding. PLoS One 10:e0098045 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Olfactory bulb;
Cell Type(s): Olfactory bulb main mitral GLU cell; Olfactory bulb main interneuron periglomerular GABA cell; Olfactory bulb main interneuron granule MC GABA cell;
Channel(s): I A; I h; I K,Ca; I Sodium; I Calcium; I Potassium;
Gap Junctions:
Receptor(s): AMPA; NMDA; Gaba;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: Python; MOOSE/PyMOOSE;
Model Concept(s): Sensory processing; Sensory coding; Markov-type model; Olfaction;
Implementer(s): Bhalla, Upinder S [bhalla at ncbs.res.in]; Gilra, Aditya [aditya_gilra -at- yahoo -period- com];
Search NeuronDB for information about:  Olfactory bulb main mitral GLU cell; Olfactory bulb main interneuron periglomerular GABA cell; Olfactory bulb main interneuron granule MC GABA cell; AMPA; NMDA; Gaba; I A; I h; I K,Ca; I Sodium; I Calcium; I Potassium; Gaba; Glutamate;
# -*- coding: utf-8 -*-
# ALL SI UNITS
# milliMolar is same as mol/m^3

## USAGE: python2.6 average_odor_morphs.py [morphs results directory]

import os,sys,math,string
import os.path
import pickle
import subprocess
cwd = os.getcwd() # current working directory

sys.path.extend(["..","../networks","../generators","../simulations"])

from stimuliConstants import * # has RESPIRATION
from pylab import * # part of matplotlib that depends on numpy but not scipy
from scipy import stats
from sim_utils import * # has rebin() and imports data_utils.py for axes_off()
from calc_corrs import * # has calc_corrs()
## Choose whether to have linear weights and rectifier: FULLlin = True
## or just linear weights and sigmoidal non-linearity: FULLlin = False
FULLlin = True
if FULLlin:
    import fit_odor_morphs_withFULLlin as fit_om # has fit_morphs()
    linextn = 'FULLlin'
else:
    import fit_odor_morphs as fit_om # has fit_morphs()
    linextn = 'lin'
## I only need residual2noise, but average_odor_pulses imports average_odor_morphs,
## hence import <...> as <...> to avoid conflict.
import average_odor_pulses as forR2N # has function to calculate residual to noise

IN_VIVO = True
## these below are often set in the respective function arguments in __main__
DIRECTED = True ## overrides the one in networkConstants.py (came in via sim_utils.py)
FRAC_DIRECTED = 0.01 ## overrides the one in networkConstants.py (came in via sim_utils.py)
## Below two overide the variables that came in from stimuliConstantsMinimal.py via stimuliConstants.py
NONLINEAR_ORNS = False
NONLINEAR_TYPE = 'P' # P for primary glom non-linear, L for lateral gloms non-linear

fullfilename = '../results/odor_morphs/morphs_random'
if NONLINEAR_ORNS: fullfilename += 'NL'+NONLINEAR_TYPE
fullfilename += '.pickle'
#fullfile = open(fullfilename,'r')
#morphs = pickle.load(fullfile)
#fullfile.close()

PLOTRESP_NUM = 1 # whether to plot 2 respiration cycles or 1
NUMBINS = 5 ## overriding that in calc_corrs.py
BIN_WIDTH_TIME = RESPIRATION/NUMBINS
bindt = RESPIRATION/float(NUMBINS)
## I take the last PLOTRESP_NUM of respiration cycles out of NUM_RESPS simulated
responsetlist = arange( SETTLETIME+(NUM_RESPS-PLOTRESP_NUM)*RESPIRATION+bindt/2.0, \
    SETTLETIME+NUM_RESPS*RESPIRATION, RESPIRATION/NUMBINS )

## IMPORTANT!!!!!! num_gloms defined below is only for some functions.
## Typically used plot_directed() has num_gloms defined in __main__ below.
salient = False#True
if salient:
    stim_seeds = [-1,-10,-19,-28]#range(-1,-37,-1)#[-1,-2,-3,-4,-5,-6,-7,-8,-9]
    num_gloms_list = [3] # used only by some functions, see __main__ for others
    ## inh_options = [ (no_singles,no_joints,no_lat,no_PGs,varyRMP), ... ]
    ## in order, below options are:
    ## all cells; no lat; no joints, varyRMP; no PGs; no singles + no joints, only mitrals
    #inh_options = [ (0,(False,False,False,False,False)), (1,(False,False,True,False,False)), \
    #    (2,(False,True,False,False,False)), (3,(False,False,False,False,True)), \
    #    (4,(False,False,False,True,False)), (5,(True,True,False,False,False)), (6,(True,True,False,True,False))]
    inh_options = [ (0,(False,False,False,False,False)), (0,(False,False,True,False,False)) ]
else:
    #stim_seeds = append(stim_seeds,[-1,-10,-19,-28])
    
    ##################################### IMPORTANT ###################################
    ###### USE ONE OR THE OTHER FOR FIGURES
    ## I used all simulations for decorrelation FIGURES
    stim_seeds = arange(750.0,1100.0,1.0)
    ## I used only 50 seeds for fittng Adil morphs.
    #stim_seeds = arange(750.0,800.0,1.0)
    
    num_gloms_list = [3] # used only by some functions, see __main__ for others
    ## inh_options = [ (no_singles,no_joints,no_lat,no_PGs,varyRMP), ... ]
    ## in order,below options are: all cells; no lat; no joints, varyRMP; no PGs; only mitrals
    #inh_options = [ (0,(False,False,False,False,False)), (1,(False,False,True,False,False)), \
    #    (2,(False,True,False,False,False)), (3,(False,False,False,False,True)), \
    #    (6,(False,False,False,True,False)), (8,(True,True,False,True,False))]
    ###### THESE inh_options ARE USED BY ONLY SOME FUNCTIONS, SEE __main__ FOR OTHERS.
    inh_options = [ (0,(False,False,False,False,False))]#, (1,(False,False,True,False,False))]
        #(2,(False,False,False,False,True)) ]
net_seeds = [100.0,200.0]

def get_filename(netseed,stimseed,inh,numgloms,stimi,neti,inhi,
        _directed=DIRECTED,_frac_directed=FRAC_DIRECTED,\
        resultsdir='../results/odor_morphs'):
    ### read filename from the output file of automated run
    #filename = morphs[ngi,stimi,neti,inhi]
    ## construct the filename
    if inh[0]: singles_str = '_NOSINGLES'
    else: singles_str = '_SINGLES'
    if inh[1]: joints_str = '_NOJOINTS'
    else: joints_str = '_JOINTS'
    if inh[3]: pgs_str = '_NOPGS'
    else: pgs_str = '_PGS'
    if inh[2]: lat_str = '_NOLAT'
    else: lat_str = '_LAT'
    if inh[4]: varmitstr = '_VARMIT'
    else: varmitstr = '_NOVARMIT'
    ## stable enough that time tags are not needed
    filename = resultsdir+'/odormorph'+\
        '_netseed'+str(netseed)+'_stimseed'+str(stimseed)
    if NONLINEAR_ORNS: filename += '_NL'+NONLINEAR_TYPE
    filename += singles_str+joints_str+pgs_str+lat_str+varmitstr+\
        '_numgloms'+str(numgloms)
    if _directed: filename += '_directed'+str(_frac_directed)
    filename += '.pickle'
    return filename, (singles_str, joints_str, pgs_str, lat_str, varmitstr)

def read_morphfile_allmits(filename):
    f = open(filename,'r')
    (mitral_responses_list,mitral_responses_binned_list) = pickle.load(f)
    f.close()
    mitral_responses_binned_list = \
        rebin(mitral_responses_list, numbins=PLOTRESP_NUM*NUMBINS,\
            bin_width_time=BIN_WIDTH_TIME, numresps=PLOTRESP_NUM)
    numavgs = len(mitral_responses_list)
    mitral_responses_avg = mean(mitral_responses_binned_list, axis=0)
    mitral_responses_std = std(mitral_responses_binned_list, axis=0)
    return numavgs, mitral_responses_avg, mitral_responses_std

def plot_responses():
    for stimi,stimseed in enumerate(stim_seeds):
        if not salient: net_seeds_here = [stimseed]
        else: net_seeds_here = net_seeds
        numsubfigs_rows = len(net_seeds_here)*len(num_gloms_list)
        numsubfigs_cols = len(inh_options)
        fig = figure(facecolor='w')
        ax = fig.add_subplot(111)
        figtext(0.1,0.94,'A,B,air: mit0 vs mit1 : stim ='+str(stimseed),fontsize=20)
        delaxes()

        for neti,netseed in enumerate(net_seeds_here):
            for ngi,num_gloms in enumerate(num_gloms_list):
                ## inh =  (no_singles,no_joints,no_lat,no_PGs,varyRMP)
                for plotyi,(inhi,inh) in enumerate(inh_options):

                    filename, switch_strs \
                        = get_filename(netseed,stimseed,inh,num_gloms,stimi,neti,inhi)
                    ## if the result file for these seeds & tweaks doesn't exist,
                    ## then carry on to the next.
                    if not os.path.exists(filename): continue
                    switches_str = string.join(switch_strs,'')
                    print filename

                    ## read in the spiketimes/binned responses for this run
                    numavgs, mitral_responses_avg, mitral_responses_std = \
                        read_morphfile_allmits(filename)

                    ## calc phase corr-s for printing on the title
                    (air_corr,odorA_corr,odorB_corr),\
                        (tcorrlist,airxcorrgram,odorAxcorrgram,odorBxcorrgram),\
                        Dfrates = \
                        calc_corrs(filename, norm_str="overall", \
                            numbins=NUMBINS, bin_width_time=BIN_WIDTH_TIME, printinfo=False)

                    ## plot the binned firing rates highlighting VARMIT differences
                    ## add_subplot(rows,cols,fignum)
                    ## fignum fills first row to end, then next row
                    ax = fig.add_subplot(numsubfigs_rows,numsubfigs_cols,\
                        (neti+ngi)*numsubfigs_cols+plotyi+1)
                    ## air
                    ax.errorbar(x=responsetlist,y=mitral_responses_avg[6,0],\
                        yerr=mitral_responses_std[6,0],color=(0,0,0),linewidth=2)
                    ax.errorbar(x=responsetlist,y=mitral_responses_avg[6,1],\
                        yerr=mitral_responses_std[6,1],color=(0,0,0.5),linewidth=2)                    
                    ## odorA
                    ax.errorbar(x=responsetlist,y=mitral_responses_avg[5,0],\
                        yerr=mitral_responses_std[5,0],color=(1,0,0),linewidth=2)
                    ax.errorbar(x=responsetlist,y=mitral_responses_avg[5,1],\
                        yerr=mitral_responses_std[5,1],color=(1,0,0.5),linewidth=2)
                    ## odorB
                    ax.errorbar(x=responsetlist,y=mitral_responses_avg[0,0],\
                        yerr=mitral_responses_std[0,0],color=(0,1,0),linewidth=2)
                    ax.errorbar(x=responsetlist,y=mitral_responses_avg[0,1],\
                        yerr=mitral_responses_std[0,1],color=(0,1,0.5),linewidth=2)
                    ax.set_title(
                        str(netseed)+str(num_gloms)+'G'+'_%1.1f_%1.1f_%1.1f'\
                        %(air_corr,odorA_corr,odorB_corr)+\
                        switches_str,fontsize=8)

def plot_decorr_single(resultsdir,stimseed=stim_seeds[0],\
        numgloms=num_gloms_list[0],inh=inh_options[0][1],odornum=None,splax=None):
    ## inh =  (no_singles,no_joints,no_lat,no_PGs,varyRMP)
    fig = figure(facecolor='w')
    ax = fig.add_subplot(111)
    figtext(0.12,0.90,'sisters: mit0 (r) vs mit1 (g)',fontsize=34)
    netseed = stimseed
    if stimseed<0:
        stimseed = int(stimseed)
        netseed = net_seeds[0]

    filename, switch_strs \
        = get_filename(netseed,stimseed,inh,numgloms,\
            None,None,None,resultsdir=resultsdir)
    switches_str = string.join(switch_strs,'')
    print filename

    ## read in the spiketimes/binned responses for this run
    numavgs, mitral_responses_avg, mitral_responses_std = \
        read_morphfile_allmits(filename)
    ## since I plot the mean response, I must plot standard error of the mean
    ## = standard deviation of a repeat / sqrt(num of repeats).
    ## NOTE: our plot depends on number of trials.
    mitral_responses_se = mitral_responses_std/sqrt(numavgs)

    ## calc phase corr-s for printing on the title
    (air_corr,odorA_corr,odorB_corr),\
        (tcorrlist,airxcorrgram,odorAxcorrgram,odorBxcorrgram),\
        Dfrates = \
        calc_corrs(filename, norm_str="overall", \
            numbins=NUMBINS, bin_width_time=BIN_WIDTH_TIME, printinfo=False)

    odor_corr = odorA_corr
    ax.errorbar(x=responsetlist,y=mitral_responses_avg[5,0],\
        yerr=mitral_responses_std[5,0],color=(1,0,0),linewidth=4)
    ax.errorbar(x=responsetlist,y=mitral_responses_avg[5,1],\
        yerr=mitral_responses_std[5,1],color=(0,1,0),linewidth=4)
    ax.set_title('corr = %1.2f'%(odor_corr,),fontsize=30)
    #ax.set_xlim(0.75,1.25)
    #ax.set_xticks([0.75,1.25])
    ax.set_xticks([0.8,0.9,1.0,1.1,1.2])
    ax.set_xticklabels(['1','2','3','4','5'])
    ## just to scale up the ticks fontsize.
    axes_labels(ax,'phase bin','firing rate (Hz)',adjustpos=False,fontsize=34)
    
    if odornum is not None and splax is not None:
        ## mit0 (red) vs mit1 (blue)
        simresponse = mitral_responses_avg[odornum,0]
        simerr = mitral_responses_se[odornum,0]
        plottimes = arange(BIN_WIDTH_TIME/2.0,RESPIRATION,BIN_WIDTH_TIME)
        #splax.fill_between(plottimes, simresponse+simerr, simresponse-simerr,
        #    color='r',alpha=0.4)
        #splax.plot(plottimes,simresponse,color='r',linewidth=linewidth)
        splax.errorbar(x=plottimes,y=simresponse,color='r',\
            marker='s',ms=marker_size,linewidth=linewidth)
        simresponse = mitral_responses_avg[odornum,1]
        simerr = mitral_responses_se[odornum,1]
        #splax.fill_between(plottimes, simresponse+simerr, simresponse-simerr,
        #    color='b',alpha=0.4)
        #splax.plot(plottimes,simresponse,color='b',linewidth=linewidth)
        splax.errorbar(x=plottimes,y=simresponse,color='b',\
            marker='o',ms=marker_size,linewidth=linewidth,linestyle='dashed')
        splax.set_xlim(0,RESPIRATION)
        beautify_plot(splax,x0min=False,y0min=True,xticksposn='bottom',yticksposn='left')
        axes_labels(splax,'time (s)','firing rate (Hz)',adjustpos=False,fontsize=label_fontsize)
        #splax.text(-0.3,0.9,'c',fontweight='bold',fontsize=label_fontsize,transform=splax.transAxes)

def plot_responses_mits_paperfigure(resultsdir,odornum,stimseed=stim_seeds[0],\
        numgloms=num_gloms_list[0],inh=inh_options[0][1]):
    netseed = stimseed
    filename, switch_strs \
        = get_filename(netseed,stimseed,inh,numgloms,\
            None,None,None,resultsdir=resultsdir)
    ## read in the spiketimes/binned responses for this run
    numavgs, mitral_responses_avg, mitral_responses_std = \
        read_morphfile_allmits(filename)
    ## plot separate figures for responses of mitral cells
    plottimes = arange(BIN_WIDTH_TIME/2.0,RESPIRATION,BIN_WIDTH_TIME)
    for miti,mitnum in enumerate([0,1,2,4]):
        fig = figure(figsize=(columnwidth/5.0,linfig_height/4.0),\
            dpi=300,facecolor='none') # none is transparent
        ax = fig.add_subplot(111)
        simresponse = mitral_responses_avg[odornum,mitnum]
        ax.errorbar(x=plottimes,y=simresponse,color=['r','m','g','b'][miti],\
            marker='o',ms=marker_size,linewidth=linewidth,\
            linestyle=['solid','dashed','solid','dashed'][miti])
        if mitnum==4:
            add_scalebar(ax,matchx=False,matchy=False,hidex=True,hidey=True,\
                sizex=0.2,labelx='0.2 s',sizey=8,labely='8 Hz',\
                bbox_to_anchor=[0.7,-0.3],bbox_transform=ax.transAxes)
        beautify_plot(ax,x0min=False,y0min=True,\
            drawxaxis=False,drawyaxis=False,xticks=[],yticks=[])
        ax.set_ylim(0,40) # same ymax to set common scale bar
        fig.tight_layout()
        fig_clip_off(fig)
        fig.savefig('../figures/decorr/response_mit'+str(mitnum)+'.svg',\
            dpi=fig.dpi,transparent=True)
        
def plot_decorrs_special_paperfigure(resultsdir,numglomslist,_inh_options=inh_options,
        _directed=DIRECTED,_frac_directed=FRAC_DIRECTED,graph=True):
    if graph:
        fig = figure(figsize=(8, 3), dpi=150, facecolor='w')
        ax = fig.add_subplot(111)
        #figtext(0.1,0.94,'sisters: mit0 (r) vs mit1 (g): frate(Hz) vs phase bin',fontsize=34)
        delaxes()
        figair = figure(facecolor='w')
        axair = figair.add_subplot(111)
        #figtext(0.1,0.94,'Air: mit0 (r) vs mit1 (g): frate(Hz) vs phase bin',fontsize=34)
        delaxes()
    all_odor_corrs = [] # all odor corrs including nan-s
    odor_corrs = [] # not nan odor corrs
    air_corrs = []
    good_odor_corrs = []
    good_air_corrs = []
    total_frate = array([0]*NUMBINS)
    num_responses = 0
    total_airfrate = array([0]*NUMBINS)
    deltafrate_sis0 = []
    deltafrate_sis1 = []
    for stimi,stimseed in enumerate(stim_seeds):
        netseed = stimseed
        if stimseed<0:
            stimseed = int(stimseed)
            netseed = net_seeds[0]
        for ngi,num_gloms in enumerate(numglomslist):
            ## inh =  (no_singles,no_joints,no_lat,no_PGs,varyRMP)
            for plotyi,(inhi,inh) in enumerate(_inh_options):

                filename, switch_strs \
                    = get_filename(netseed,stimseed,inh,num_gloms,None,None,None,\
                        _directed,_frac_directed,resultsdir)
                ## if the result file for these seeds & tweaks doesn't exist,
                ## then carry on to the next.
                if not os.path.exists(filename): continue
                switches_str = string.join(switch_strs,'')
                #if graph: print filename

                ## read in the spiketimes/binned responses for this run
                f = open(filename,'r')
                (mitral_responses_list,mitral_responses_binned_list) = pickle.load(f)
                f.close()
                ## rebin the spikes
                mitral_responses_binned_list = \
                    rebin(mitral_responses_list, numbins=PLOTRESP_NUM*NUMBINS,\
                        bin_width_time=BIN_WIDTH_TIME, numresps=PLOTRESP_NUM)
                numavgs = len(mitral_responses_list)
                mitral_responses_avg = mean(mitral_responses_binned_list, axis=0)
                mitral_responses_std = std(mitral_responses_binned_list, axis=0)
                ## mean 1% odor firing rate
                total_frate += mitral_responses_avg[0,0]+\
                    mitral_responses_avg[0,1]+\
                    mitral_responses_avg[5,0]+\
                    mitral_responses_avg[5,1]
                num_responses += 4
                total_airfrate += mitral_responses_avg[6,0] + mitral_responses_avg[6,1]
                #if graph:
                #    print "Odors:"
                #    print mitral_responses_avg[0,0]
                #    print mitral_responses_avg[0,1]
                #    print mitral_responses_avg[5,0]
                #    print mitral_responses_avg[5,1]
                #    print "Airs:"
                #    print mitral_responses_avg[6,0]
                #    print mitral_responses_avg[6,1]
                
                ## Calc change in firing rate for each odor, for each sister
                numphasebins = float(len(mitral_responses_avg[0,0]))
                for odornum in [0,5]:
                    deltafrate0 = sum(mitral_responses_avg[odornum,0])/numphasebins \
                        - sum(mitral_responses_avg[6,0])/numphasebins
                    deltafrate1 = sum(mitral_responses_avg[odornum,1])/numphasebins \
                        - sum(mitral_responses_avg[6,1])/numphasebins
                    deltafrate_sis0.append(deltafrate0)
                    deltafrate_sis1.append(deltafrate1)

                ## calc phase corr-s for printing on the title
                (air_corr,odorA_corr,odorB_corr),\
                    (tcorrlist,airxcorrgram,odorAxcorrgram,odorBxcorrgram),\
                    Dfrates = \
                    calc_corrs(filename, norm_str="overall", \
                        numbins=NUMBINS, bin_width_time=BIN_WIDTH_TIME, printinfo=False)
                all_odor_corrs.append(odorA_corr)
                all_odor_corrs.append(odorB_corr)
                ## If one of the phasic responses is all zeroes, stats.pearsonr gives a one-time warning
                ## and returns nan-s, but nan-s are removed by below, for obtaining the distribution
                if not isnan(odorA_corr): odor_corrs.append(odorA_corr)
                if not isnan(odorB_corr): odor_corrs.append(odorB_corr)
                if not isnan(air_corr): air_corrs.append(air_corr)
                print stimseed, odorA_corr, odorB_corr, air_corr

                corr_cutoff = -0.5
                numrows = 4
                numcols = 6
                numplots = numrows*numcols
                ## nan-s always return false on comparison with non-nan, hence ignored
                ### fignum fills first row to end, then next row
                ## odorA
                if odorA_corr<corr_cutoff:
                    good_odor_corrs.append(odorA_corr)
                    if graph and len(good_odor_corrs)<=numplots:
                        ax = fig.add_subplot(numrows,numcols,len(good_odor_corrs))
                        ax.errorbar(x=responsetlist,y=mitral_responses_avg[5,0],\
                            yerr=mitral_responses_std[5,0],color=(1,0,0),linewidth=4)
                        ax.errorbar(x=responsetlist,y=mitral_responses_avg[5,1],\
                            yerr=mitral_responses_std[5,1],color=(0,1,0),linewidth=4)
                        ax.set_title('%1.2f'%(good_odor_corrs[-1]),fontsize=16)
                        ymax = ax.get_ylim()[1]
                        ax.set_yticks([0,ymax])
                        ax.set_yticklabels(['0',str(ymax)])
                        #ax.set_xlim(0.75,1.25)
                        #ax.set_xticks([0.75,1.25])
                        ax.set_xticks([])
                        #ax.set_xticklabels(['0.75','1.25'])
                        ## just to scale up the ticks fontsize.
                        axes_labels(ax,'','',adjustpos=False,fontsize=18)
                        #print filename, "corr =", good_odor_corrs[-1]
                ## odorB
                if odorB_corr<corr_cutoff:
                    good_odor_corrs.append(odorB_corr)
                    if graph and len(good_odor_corrs)<=numplots:
                        ax = fig.add_subplot(numrows,numcols,len(good_odor_corrs))
                        ax.errorbar(x=responsetlist,y=mitral_responses_avg[0,0],\
                            yerr=mitral_responses_std[0,0],color=(1,0,0),linewidth=4)
                        ax.errorbar(x=responsetlist,y=mitral_responses_avg[0,1],\
                            yerr=mitral_responses_std[0,1],color=(0,1,0),linewidth=4)
                        ax.set_title('%1.2f'%(good_odor_corrs[-1]),fontsize=16)
                        ymax = ax.get_ylim()[1]
                        ax.set_yticks([0,ymax])
                        ax.set_yticklabels(['0',str(ymax)])
                        #ax.set_xlim(0.75,1.25)
                        #ax.set_xticks([0.75,1.25])
                        ax.set_xticks([])
                        #ax.set_xticklabels(['0.75','1.25'])
                        ## just to scale up the ticks fontsize.
                        axes_labels(ax,'','',adjustpos=False,fontsize=18)
                        #print filename, "corr =", good_odor_corrs[-1]
                ## air
                if air_corr<corr_cutoff:
                    good_air_corrs.append(air_corr)
                    if graph and len(good_air_corrs)<=numplots:
                        axair = figair.add_subplot(numrows,numcols,len(good_air_corrs))
                        axair.errorbar(x=responsetlist,y=mitral_responses_avg[6,0],\
                            yerr=mitral_responses_std[6,0],color=(1,0,0),linewidth=4)
                        axair.errorbar(x=responsetlist,y=mitral_responses_avg[6,1],\
                            yerr=mitral_responses_std[6,1],color=(0,1,0),linewidth=4)
                        axair.set_title('%1.2f'%(good_air_corrs[-1]),fontsize=16)
                        ymax = ax.get_ylim()[1]
                        ax.set_yticks([0,ymax])
                        ax.set_yticklabels(['0',str(ymax)])
                        #axair.set_xlim(0.75,1.25)
                        #axair.set_xticks([0.75,1.25])
                        axair.set_xticks([])
                        #axair.set_xticklabels(['0.75','1.25'])
                        ## just to scale up the ticks fontsize.
                        axes_labels(axair,'','',adjustpos=False,fontsize=18)
                        #print filename, "corr =", good_air_corrs[-1]

    ## correlation in change in firing rate
    corr_deltafrate = \
        dot(deltafrate_sis0,deltafrate_sis1)/ \
            sqrt(dot(deltafrate_sis0,deltafrate_sis0)*dot(deltafrate_sis1,deltafrate_sis1))
    ## * is element wise multiplication.
    ## contribution of each mitral-pair--odor combo to the deltafrate corr.
    ## normalized by the normalization of the dot-product.
    corr_deltafrate_eachmit = array(deltafrate_sis0)*array(deltafrate_sis1) \
            / sqrt(dot(deltafrate_sis0,deltafrate_sis0)*dot(deltafrate_sis1,deltafrate_sis1))
                    
    if num_responses==0:
        print "No result files found ..."
        sys.exit(1)
    overall_odor_mean = sum(total_frate/float(num_responses))/float(len(total_frate))
    overall_air_mean = sum(total_airfrate/float(num_responses/2.0))/float(len(total_airfrate))
    print "Phasic odor mean firing rate =",total_frate/float(num_responses)
    print "Overall odor mean = ", overall_odor_mean
    print "Phasic air mean firing rate =",total_airfrate/float(num_responses/2.0)
    print "Overall air mean = ", overall_air_mean
    print "delta frates of sisters :",zip(deltafrate_sis0,deltafrate_sis1)
    print "Correlation between sisters' change in firing rate with odor",corr_deltafrate

    if graph:
        fig.tight_layout() # for previous figure

        ## paper figure of decorr example and correlation distribution
        fig = figure(figsize=(columnwidth,linfig_height/2.0),dpi=300,facecolor='w')
        splax = fig.add_subplot(1,2,1)
        ###### paper figure for example decorrelation.
        ## odornum should be 5 (odor A) or 0 (odor B).
        plot_decorr_single(resultsdir=resultsdir,stimseed=844.0,numgloms=3,\
            inh=(False,False,False,False,False),odornum=0,splax=splax)

        ax = fig.add_subplot(1,2,2)
        #title("PDF of correlations between non-zero sister responses",fontsize=30)
        ax.hist(air_corrs,10,range=(-1.0,1.0),normed=True,histtype='step',\
            linewidth=linewidth,color='b',ls='dotted')
        ax.hist(odor_corrs,10,range=(-1.0,1.0),normed=True,histtype='step',\
            linewidth=linewidth,color='r',ls='solid')
        #ax.text(-0.3, 0.9, 'd', fontweight='bold', fontsize=label_fontsize, transform=ax.transAxes)
        beautify_plot(ax,x0min=False,y0min=True,xticksposn='bottom',yticksposn='left')
        axes_labels(ax,'correlation','prob. density',adjustpos=False)
        ax.set_xticks([-1,0,1])
        fig.tight_layout()
        fig_clip_off(fig)
        fig.savefig('../figures/decorr/sim_decorr.svg', dpi=fig.dpi)
        fig.savefig('../figures/decorr/sim_decorr.png', dpi=fig.dpi)
    
        fig = figure(figsize=(columnwidth/3.0, linfig_height/2.0), dpi=300, facecolor='w')
        ax = fig.add_subplot(111)
        #title("PDF of correlations between non-zero sister responses",fontsize=30)
        scatter(all_odor_corrs,corr_deltafrate_eachmit,marker='o',s=marker_size)
        axes_labels(ax,'correlation','deltafrate',adjustpos=False)
        #biglegend('upper left')
        fig.tight_layout()

    return corr_deltafrate, odor_corrs, air_corrs, overall_odor_mean, overall_air_mean

def plot_xcorrgrams():
    ## plot the average xcorrgrams over all stimuli
    ## for each netseed (figure) and inhibition type (panel)
    for neti,netseed in enumerate(net_seeds):
        
        numsubfigs_rows = len(num_gloms_list)
        numsubfigs_cols = len(inh_options)
        ## fig with subpanels for spike corr
        figxcorr = figure(facecolor='w')
        figtext(0.1,0.94,'odor A,B,air xcorrs between sisters: net'+str(netseed),fontsize=20)
        delaxes() # delete the main axes to accomodate subpanels
        ## printed output header
        print 'odors A, B, air phase similarity / corr between sisters: netseed = '+str(netseed)
        ## fig with subpanels for phase corr / similarity
        figcorr = figure(facecolor='w')
        figtext(0.1,0.94,'odorA,B,air ph sim/corr between sisters: net'+str(netseed),fontsize=20)
        delaxes() # delete the main axes to accomodate subpanels
        
        for ngi,num_gloms in enumerate(num_gloms_list):

            ## inh =  (no_singles,no_joints,no_lat,no_PGs,varyRMP)
            for plotyi,(inhi,inh) in enumerate(inh_options):

                ## spike corrs for all odors and airs to this network instance
                air_xcorrgrams = []
                odorA_xcorrgrams = []
                odorB_xcorrgrams = []
                ## phase tuning corrs for all odors and airs to this network instance
                air_corrs = []
                odorA_corrs = []
                odorB_corrs = []
                ## average change in frate (compared to air) for all odors
                mit0_Dfrates = []
                mit1_Dfrates = []
                mit2_Dfrates = []
                mit3_Dfrates = []

                for stimi,stimseed in enumerate(stim_seeds):

                    filename, switch_strs = \
                        get_filename(netseed,stimseed,inh,num_gloms,stimi,neti,inhi)
                    ## if the result file for these seeds & tweaks doesn't exist,
                    ## then carry on to the next.
                    if not os.path.exists(filename): continue
                    switches_str = string.join(switch_strs,'')

                    ## calculate spike time based xcorrs:
                    (air_corr,odorA_corr,odorB_corr),\
                        (tcorrlist,airxcorrgram,odorAxcorrgram,odorBxcorrgram),\
                        Dfrates = \
                        calc_corrs(filename, norm_str="overall", \
                            numbins=NUMBINS, bin_width_time=BIN_WIDTH_TIME, printinfo=False)
                    if air_corr<0 or odorA_corr<0 or odorB_corr<0:
                        print filename
                        print "air corr =",air_corr,"odor A corr =",odorA_corr,"odor B corr =",odorB_corr
                    ## if first value in list is nan, so are all the others; exclude them
                    if not math.isnan(airxcorrgram[0]): air_xcorrgrams.append(airxcorrgram)
                    if not math.isnan(odorAxcorrgram[0]): odorA_xcorrgrams.append(odorAxcorrgram)
                    if not math.isnan(odorBxcorrgram[0]): odorB_xcorrgrams.append(odorBxcorrgram)
                    ## if the phase similarity / correlation is nan, exclude
                    if not math.isnan(air_corr): air_corrs.append(air_corr)
                    if not math.isnan(odorA_corr): odorA_corrs.append(odorA_corr)
                    if not math.isnan(odorB_corr): odorB_corrs.append(odorB_corr)
                    ## average Delta frates for the sisters 0 and 1 -- odor A and B
                    mit0_Dfrates.extend((Dfrates[0],Dfrates[2]))
                    mit1_Dfrates.extend((Dfrates[1],Dfrates[3]))
                    #mit2_Dfrates.extend((Dfrates[4],Dfrates[6]))
                    #mit3_Dfrates.extend((Dfrates[5],Dfrates[7]))
                
                #### plot various analyses for this network instance + conditions

                ## plot spike time based xcorrs:
                ax = figxcorr.add_subplot(numsubfigs_rows,numsubfigs_cols,\
                    (ngi)*numsubfigs_cols+plotyi+1)
                if len(air_xcorrgrams)>0:
                    ax.plot(tcorrlist, mean(air_xcorrgrams,axis=0), color=(0,0,0), label='air')
                if len(odorA_xcorrgrams)>0:
                    ax.plot(tcorrlist, mean(odorA_xcorrgrams,axis=0), color=(1,0,0), label='odor A')
                if len(odorB_xcorrgrams)>0:
                    ax.plot(tcorrlist, mean(odorB_xcorrgrams,axis=0), color=(0,1,0), label='odor B')
                ax.set_title(
                    str(num_gloms)+'G'+switches_str,fontsize=8)

                ## print phase similarity / corrs
                print '\t'+str(num_gloms)+'Glomeruli'+switches_str
                print "\t\tAir phase similarity =",mean(air_corrs),'SD',std(air_corrs)
                print "\t\tOdorA phase similarity =",mean(odorA_corrs),'SD',std(odorA_corrs)
                print "\t\tOdorB phase similarity =",mean(odorB_corrs),'SD',std(odorB_corrs)

                ## plot histogram of phase corr / similarity between sisters for air & odors:
                ## similar to plots in fig 7e of Ashesh Dhawale et al
                ax = figcorr.add_subplot(numsubfigs_rows,numsubfigs_cols,\
                    (ngi)*numsubfigs_cols+plotyi+1)
                ## typically air_corrs will have less data points as many are nan-s
                ## how to normalize? normed=True is not what I want.
                if len(air_corrs)>0:
                    ax.hist(air_corrs, bins=5, histtype='step', normed=False, color=(0,0,0), label='air')
                if len(odorA_corrs)>0:
                    ax.hist(odorA_corrs, bins=5, histtype='step', normed=False, color=(1,0,0), label='odor A')
                if len(odorB_corrs)>0:
                    ax.hist(odorB_corrs, bins=5, histtype='step', normed=False, color=(0,1,0), label='odor B')
                ax.set_title(str(num_gloms)+'G'+switches_str,fontsize=8)
                ax.set_xlim(-1.0,1.0)

                ## print corr between mit0 and mit1 mean Dfrates:
                #print "\t\tmit0 mean Delta firing rate change (Hz) for various odors =",mit0_Dfrates
                #print "\t\tmit1 mean Delta firing rate change (Hz) for various odors =",mit1_Dfrates
                print "\t\tCorrelation of 'mean change in firing rate' Odor Response Spectrum (Ashesh et al)"
                print "\t\t\tbetween sisters 0 and 1 =",\
                    stats.pearsonr(mit0_Dfrates,mit1_Dfrates)[0]
                #print "\t\t\tbetween non-sisters 0 and 2 =",\
                #    stats.pearsonr(mit0_Dfrates,mit2_Dfrates)[0]

def plot_directed(glomnums):
    """ Be sure to set frac_directed=0.0/0.05 above at the very beginning. """
    odor_corrs_means = []
    odor_corrs_SDs = []
    air_corrs_means = []
    air_corrs_SDs = []
    corrs_deltafrate = []
    fig = figure()
    for gni,glomnum in enumerate(glomnums):
        print "Computing phasic and deltafrate correlations for # of gloms =",glomnum
        ## Set graph=True below to plot neg corr-ed responses too.
        corr_deltafrate, odor_corrs, air_corrs, overall_odor_mean, overall_air_mean = \
            plot_decorrs_special([glomnum],graph=True)
        ax = fig.add_subplot(len(glomnums),1,gni+1)
        #hist(air_corrs,20,range=(-1.0,1.0),normed=True,histtype='step',\
        #    color='b',linewidth=2,label='air %2.1f'%overall_air_mean+'Hz')
        hist(odor_corrs,20,range=(-1.0,1.0),normed=True,histtype='step',\
            color='r',linewidth=2,label='odor %2.1f'%overall_odor_mean+'Hz')
        ax.set_xticks([])
        #ax.set_xticklabels(['0.75','1.25'])
        ## just to scale up the ticks fontsize.
        axes_labels(ax,'','',adjustpos=False,fontsize=34)

        corrs_deltafrate.append(corr_deltafrate)
        ## mean and SD of phasic correlations of odor and air
        odor_corrs_means.append(mean(odor_corrs))
        odor_corrs_SDs.append(std(odor_corrs))
        air_corrs_means.append(mean(air_corrs))
        air_corrs_SDs.append(std(air_corrs))

        ax.set_yticks([])
        #biglegend(legendlocation='upper left')
        if gni == len(glomnums)-1:
            ax.set_xticks([-1.0,0.0,1.0])
            ax.set_xticklabels(['-1','0','1'])
            axes_labels(ax,'phase correlation','',adjustpos=False,fontsize=30)
    plt.tight_layout()

    ## mean phase corr vs number of connected gloms
    fig=figure()
    ax=fig.add_subplot(111)
    #plot(glomnums,air_corrs_means,color='b',linewidth=2,label='air')
    plot(glomnums,odor_corrs_means,color='r',linewidth=2,label='odor')
    ax.set_xticks(glomnums)
    ax.set_xticklabels([str(glomnum) for glomnum in glomnums])
    axes_labels(ax,'# of connected glomeruli','phase correlation mean',\
        adjustpos=False,fontsize=30)
    #biglegend(legendlocation='lower left')
    plt.tight_layout()
    ## spread of phase corr vs number of connected gloms
    fig=figure()
    ax=fig.add_subplot(111)
    #errorbar(glomnums,air_corrs_SDs,color='b',linewidth=2,label='air')
    errorbar(glomnums,odor_corrs_SDs,color='r',linewidth=2,label='odor')
    ax.set_xticks(glomnums)
    ax.set_xticklabels([str(glomnum) for glomnum in glomnums])
    axes_labels(ax,'# of connected glomeruli','phase correlation spread',\
        adjustpos=False,fontsize=30)
    #biglegend(legendlocation='upper left')
    plt.tight_layout()
    ## delta frate corr vs number of connected gloms
    fig=figure()
    ax=fig.add_subplot(111)
    plot(glomnums,corrs_deltafrate,color='b',linewidth=2)
    ax.set_xticks(glomnums)
    ax.set_xticklabels([str(glomnum) for glomnum in glomnums])
    axes_labels(ax,'# of connected glomeruli','$\Delta$frate correlation',\
        adjustpos=False,fontsize=30)
    tight_layout()
    
def plot_across_sims_paperfigure():
    odor_corrs_means = []
    odor_corrs_SDs = []
    air_corrs_means = []
    air_corrs_SDs = []
    corrs_deltafrate = []
    fig = figure(figsize=(columnwidth/2.0,linfig_height*2.0),dpi=300,facecolor='w')
    ## [(dirextn, directed, frac_directed, inh_options, numgloms, color),...]
    ## [ ('', directed, 0.0, novarmit, [...]), ('', directed, 0.0, varmit, [...]), ('', directed, 0.05, novarmit, [...]) ]
    ## inh_options = [ (no_singles,no_joints,no_lat,no_PGs,varyRMP), ... ]
    #decorr_net_options = [ (True,0.0,(0,(False,False,False,False,False)),[2,3,6],'b'),
    #    (True,0.0,(0,(False,False,False,False,True)),[2,3],'g'),
    #    (True,0.05,(0,(False,False,False,False,False)),[2,3,4,5,6,7],'r') ]
    ####### CAUTION: EVEN THOUGH folder suffix is _0.33x_may14; the scaling is 0.5x odor, 1x air, 1x background
    decorr_net_options = [ \
        ('_0.33x_may14',False,0.0,(0,(False,False,False,False,False)),[3],'k'),\
        ('_0.33x_may14',True,0.0,(1,(False,False,False,False,False)),[3],'k'),\
        ('_0.33x_may14',True,0.0,(2,(False,False,False,False,True)),[3],'k'),\
        ('',True,0.01,(3,(False,False,False,False,False)),[3],'k'),\
        ('_2x4x_may16',True,0.01,(4,(False,False,False,False,False)),[7],'k') ]
    #decorr_net_options = [ (True,0.0,(0,(False,False,False,False,False)),[2],'b'),
    #    (True,0.05,(0,(False,False,False,False,False)),[3],'r') ]
    neti=0
    indexstrs = []
    colorslist = []
    numsubplots = sum([1 for net_option in decorr_net_options for glomnum in net_option[4]])
    for (dirextn,directed,frac_directed,inh_option,numglomslist,color) in decorr_net_options:
        for numgloms in numglomslist:
            print "Computing phasic and deltafrate correlations for # of gloms =",numgloms
            print "The inh tweak (no_singles,no_joints,no_lat,no_PGs,varyRMP) is :",inh_option
            ## Set graph=True below to plot neg corr-ed responses too.
            corr_deltafrate, odor_corrs, air_corrs, overall_odor_mean, overall_air_mean = \
                plot_decorrs_special_paperfigure('../results/odor_morphs'+dirextn, \
                    [numgloms],[inh_option],directed,frac_directed,graph=False)
            ax = fig.add_subplot(numsubplots,1,neti+1)
            #hist(air_corrs,20,range=(-1.0,1.0),normed=True,histtype='step',\
            #    color='b',linewidth=2,label='air %2.1f'%overall_air_mean+'Hz')
            counts,bins,patches = hist(odor_corrs,10,range=(-1.0,1.0),normed=True,histtype='step',\
                color=color,linewidth=linewidth,label='odor %2.1f'%overall_odor_mean+'Hz')
            bar(bins[:5],counts[:5],width=bins[1]-bins[0],facecolor='r',edgecolor='r')
            ax.set_xticks([])
            ax.set_yticks([])

            corrs_deltafrate.append(corr_deltafrate)
            ## mean and SD of phasic correlations of odor and air
            odor_corrs_mean = mean(odor_corrs)
            odor_corrs_SD = std(odor_corrs)
            odor_corrs_means.append(odor_corrs_mean)
            odor_corrs_SDs.append(odor_corrs_SD)
            air_corrs_means.append(mean(air_corrs))
            air_corrs_SDs.append(std(air_corrs))
            if neti==0:
                ymax = ax.get_ylim()[1]*1.2 # hard coded to 1.2x initial-ymax, to fit all subplots!
                ylim = (0.0,ymax)
                ax.set_ylim(ylim)
                ymid = ymax/2.0
            else:
                ax.set_ylim(ylim)
            ### The mean and SD are obvious from the figure.
            ### Drawing extra lines/bars like below are distracting -- hence commented out.
            ### draw an arrow of default 'Curve' style (only line) for the mean on the correlation distribution
            #arrmean = matplotlib.patches.FancyArrowPatch(
            #    (odor_corrs_mean, ymax*0.75), (odor_corrs_mean, ymax*0.25), linewidth=linewidth)
            #ax.add_patch(arrmean)
            ### draw an arrow of style BarAB for the SD on the correlation distribution
            #arrSD = matplotlib.patches.FancyArrowPatch(
            #    (odor_corrs_mean-odor_corrs_SD, ymid), (odor_corrs_mean+odor_corrs_SD, ymid), linewidth=linewidth )
            #arrSD.set_arrowstyle('|-|',widthA=10,angleA=90,widthB=10,angleB=90)
            #ax.add_patch(arrSD)
            beautify_plot(ax,x0min=False,y0min=True,xticksposn='bottom',yticksposn='left')
            if neti==2:
                axes_labels(ax,'','probability density',adjustpos=False,fontsize=label_fontsize)
            if neti>=numsubplots-1:
                ax.set_xticks([-1,0,1])
                ax.set_xticklabels(['-1','0','1'])
                ## scale up the ticks and axes labels fontsize.
                axes_labels(ax,'correlation','',adjustpos=False,fontsize=label_fontsize)
            else: ax.set_xticks([])
            #ax.text(0.15, 1.0, ['d','e','f','g','h','i'][neti], \
            #    fontweight='bold', fontsize=label_fontsize, transform=ax.transAxes)
            ## Print the delta-frate correlation on the plot/histogram
            ax.text(0.15, 0.75, '%1.3f'%(corr_deltafrate), fontsize=label_fontsize, transform=ax.transAxes)

            indexstr = ''
            if directed: indexstr += 'directed, '
            if frac_directed>0.0: indexstr += 'differential, '
            if inh_option[1][4]: indexstr += 'varymits, '
            indexstr += 'latgloms '+str(numgloms-1)
            indexstrs.append(indexstr)
            colorslist.append(color)
            neti+=1

    fig.tight_layout()
    fig.subplots_adjust(top=0.95)
    fig_clip_off(fig)
    fig.savefig('../figures/decorr/decorr_conns.svg',dpi=fig.dpi)
    fig.savefig('../figures/decorr/decorr_conns.png',dpi=fig.dpi)

    ## figure for plotting mean and SD of the distribution
    mainfig = figure(figsize=(8, 6), dpi=150, facecolor='w')
    ax=mainfig.add_subplot(111)
    indices = range(neti)
    ax.bar(indices, odor_corrs_means, color=colorslist,
       align='center', yerr=odor_corrs_SDs, ecolor='black', width=0.3)
    ax.set_ylim([-0.25,1.0])
    ax.set_yticks([-0.25,0,0.5,1.0])
    ax.set_yticklabels(['-0.25','0','0.5','1'])
    ax.set_xticks(indices)
    ax.set_xticklabels(indexstrs)
    mainfig.autofmt_xdate() # rotates xlabels
    axes_labels(ax,'','',\
        adjustpos=False,fontsize=16)
    mainfig.tight_layout()

def plot_peaks_tufted_vs_mitrals_paper_figure():
    fig = figure(figsize=(columnwidth/2.0,linfig_height*2.0),dpi=300,facecolor='w')
    ## [(dirextn, directed, frac_directed, inh_options, numglomslist, color),...]
    ## inh_options = [ (no_singles,no_joints,no_lat,no_PGs,varyRMP), ... ]
    ####### CAUTION: EVEN THOUGH folder suffix is _0.33x_may14; the scaling is 0.5x odor, 1x air, 1x background
    decorr_net_options = [ \
        ('',True,0.01,(1,(False,False,False,False,False)),[2],'k'),\
        ('',True,0.01,(2,(False,False,False,True,False)),[2],'k') ]
    neti=0
    indexstrs = []
    colorslist = []
    numsubplots = sum([1 for net_option in decorr_net_options for glomnum in net_option[4]])
    for (dirextn,directed,frac_directed,inh_option,numglomslist,color) in decorr_net_options:
        for numgloms in numglomslist:
            print "Computing phasic and deltafrate correlations for # of gloms =",numgloms
            print "The inh tweak (no_singles,no_joints,no_lat,no_PGs,varyRMP) is :",inh_option
            ## Set graph=True below to plot neg corr-ed responses too.
            corr_deltafrate, odor_corrs, air_corrs, overall_odor_mean, overall_air_mean = \
                plot_decorrs_special_paperfigure('../results/odor_morphs'+dirextn, \
                    [numgloms],[inh_option],directed,frac_directed,graph=False)
            ax = fig.add_subplot(numsubplots,1,neti+1)
            #hist(air_corrs,20,range=(-1.0,1.0),normed=True,histtype='step',\
            #    color='b',linewidth=2,label='air %2.1f'%overall_air_mean+'Hz')
            counts,bins,patches = hist(odor_corrs,10,range=(-1.0,1.0),normed=True,histtype='step',\
                color=color,linewidth=linewidth,label='odor %2.1f'%overall_odor_mean+'Hz')
            bar(bins[:5],counts[:5],width=bins[1]-bins[0],facecolor='r',edgecolor='r')
            ax.set_xticks([])
            ax.set_yticks([])

            corrs_deltafrate.append(corr_deltafrate)
            ## mean and SD of phasic correlations of odor and air
            odor_corrs_mean = mean(odor_corrs)
            odor_corrs_SD = std(odor_corrs)
            odor_corrs_means.append(odor_corrs_mean)
            odor_corrs_SDs.append(odor_corrs_SD)
            air_corrs_means.append(mean(air_corrs))
            air_corrs_SDs.append(std(air_corrs))
            if neti==0:
                ymax = ax.get_ylim()[1]*1.2 # hard coded to 1.2x initial-ymax, to fit all subplots!
                ylim = (0.0,ymax)
                ax.set_ylim(ylim)
                ymid = ymax/2.0
            else:
                ax.set_ylim(ylim)
            ### The mean and SD are obvious from the figure.
            ### Drawing extra lines/bars like below are distracting -- hence commented out.
            ### draw an arrow of default 'Curve' style (only line) for the mean on the correlation distribution
            #arrmean = matplotlib.patches.FancyArrowPatch(
            #    (odor_corrs_mean, ymax*0.75), (odor_corrs_mean, ymax*0.25), linewidth=linewidth)
            #ax.add_patch(arrmean)
            ### draw an arrow of style BarAB for the SD on the correlation distribution
            #arrSD = matplotlib.patches.FancyArrowPatch(
            #    (odor_corrs_mean-odor_corrs_SD, ymid), (odor_corrs_mean+odor_corrs_SD, ymid), linewidth=linewidth )
            #arrSD.set_arrowstyle('|-|',widthA=10,angleA=90,widthB=10,angleB=90)
            #ax.add_patch(arrSD)
            beautify_plot(ax,x0min=False,y0min=True,xticksposn='bottom',yticksposn='left')
            if neti==2:
                axes_labels(ax,'','probability density',adjustpos=False,fontsize=label_fontsize)
            if neti>=numsubplots-1:
                ax.set_xticks([-1,0,1])
                ax.set_xticklabels(['-1','0','1'])
                ## scale up the ticks and axes labels fontsize.
                axes_labels(ax,'correlation','',adjustpos=False,fontsize=label_fontsize)
            else: ax.set_xticks([])
            #ax.text(0.15, 1.0, ['d','e','f','g','h','i'][neti], \
            #    fontweight='bold', fontsize=label_fontsize, transform=ax.transAxes)
            ## Print the delta-frate correlation on the plot/histogram
            ax.text(0.15, 0.75, '%1.3f'%(corr_deltafrate), fontsize=label_fontsize, transform=ax.transAxes)

            indexstr = ''
            if directed: indexstr += 'directed, '
            if frac_directed>0.0: indexstr += 'differential, '
            if inh_option[1][4]: indexstr += 'varymits, '
            indexstr += 'latgloms '+str(numgloms-1)
            indexstrs.append(indexstr)
            colorslist.append(color)
            neti+=1

    fig.tight_layout()
    fig.subplots_adjust(top=0.95)
    fig_clip_off(fig)
    fig.savefig('../figures/decorr/decorr_conns.svg',dpi=fig.dpi)
    fig.savefig('../figures/decorr/decorr_conns.png',dpi=fig.dpi)

    ## figure for plotting mean and SD of the distribution
    mainfig = figure(figsize=(8, 6), dpi=150, facecolor='w')
    ax=mainfig.add_subplot(111)
    indices = range(neti)
    ax.bar(indices, odor_corrs_means, color=colorslist,
       align='center', yerr=odor_corrs_SDs, ecolor='black', width=0.3)
    ax.set_ylim([-0.25,1.0])
    ax.set_yticks([-0.25,0,0.5,1.0])
    ax.set_yticklabels(['-0.25','0','0.5','1'])
    ax.set_xticks(indices)
    ax.set_xticklabels(indexstrs)
    mainfig.autofmt_xdate() # rotates xlabels
    axes_labels(ax,'','',\
        adjustpos=False,fontsize=16)
    mainfig.tight_layout()

## Obsolete: This is called by fit_odor_morphs.py as a panel in Figure.
def plot_chisq_hist_paperfigure(ax1,ax2,resultsdir='../results/odor_morphs'):
    """ Plot chi-sq histogram of the morph fits. """
    #fig = figure(figsize=(columnwidth/3.0,columnwidth/2.0),dpi=300,facecolor='w') # 'none' is transparent
    #ax = fig.add_subplot(111,frameon=False)
    ## inh =  (no_singles,no_joints,no_lat,no_PGs,varyRMP)
    inh_options = [ (0,(False,False,False,False,False),'lat inh') ]
    for ploti,(inhi,inh,inhstr) in enumerate(inh_options):
        chisqs = []
        lin_chisqs = []
        n_accept = 0
        for stimi,stimseed in enumerate(stim_seeds):
            if not salient: net_seeds = [stimseed]
            for neti,netseed in enumerate(net_seeds):
                for ngi,num_gloms in enumerate([3]):

                    filename, switch_strs \
                        = get_filename(netseed,stimseed,inh,num_gloms,stimi,neti,inhi,resultsdir=resultsdir)
                    switches_str = string.join(switch_strs,'')
                    ## if the result file for these seeds & tweaks doesn't exist,
                    ## then carry on to the next.
                    if not os.path.exists(filename): continue
                    #print filename
                    for fitted_mitral in [0,1]:
                        ## First the weighted-linear sigmoid:
                        ## If the fitted params file does not exist, create it (them).
                        if not os.path.exists(filename+'_params'+str(fitted_mitral)):
                            print "fitting file",filename
                            ## fit the responses for this result file
                            params,chisq,inputsA,inputsB,fitted_responses,\
                                numavgs,firingbinsmeanList,firingbinserrList\
                                = fit_om.fit_morphs(filename, fitted_mitral, 'arb')
                        else:
                            f = open(filename+'_params'+str(fitted_mitral),'r')
                            params,chisq = pickle.load(f)
                            f.close()
                        chisqs.append(chisq)
                        ## linear-sigmoid [perhaps get rid of sigmoid also a la Priyanka?]
                        ## If the fitted params file does not exist, create it (them).
                        if not os.path.exists(filename+'_paramslin'+str(fitted_mitral)):
                            ## fit the responses for this result file
                            params,chisq,inputsA,inputsB,fitted_responses,\
                                numavgs,firingbinsmeanList,firingbinserrList\
                                = fit_om.fit_morphs(filename, fitted_mitral, 'lin')
                        else:
                            f = open(filename+'_paramslin'+str(fitted_mitral),'r')
                            params,chisq = pickle.load(f)
                            f.close()
                        lin_chisqs.append(chisq)
                        n_accept += 1

        ax1.hist(chisqs,20,histtype='step',linewidth=linewidth,label=inhstr,color='k')
        ax2.hist(lin_chisqs,20,histtype='step',linewidth=linewidth,label=inhstr,color='k')
        print "Number of mitral cells accepted =",n_accept
        
        ## beautify plots
        for axnum,ax in enumerate([ax1,ax2]):
            ## ax.transAxes ensures relative to axes size, rather than to data units.
            #text(0.15, 1.0, ['E','F'][axnum], fontsize=label_fontsize, transform = ax.transAxes)
            ax.get_yaxis().set_ticks_position('left')
            ax.get_xaxis().set_ticks_position('bottom')
            xmin, xmax = ax.get_xaxis().get_view_interval()
            ymin, ymax = ax.get_yaxis().get_view_interval()
            ax.set_xlim(0,xmax)
            ax.set_ylim(0,ymax)
            ax.set_xticks([0,xmax])
            ax.set_yticks([0,ymax])
            ax.add_artist(Line2D((0, 0), (0, ymax), color='black', linewidth=axes_linewidth))
            ax.add_artist(Line2D((0, xmax), (0, 0), color='black', linewidth=axes_linewidth))
            axes_labels(ax,'','',adjustpos=False) # sets font-size for tick labels also
        ax2.text(-0.4,1.3,'count',fontsize=label_fontsize, rotation='vertical', transform=ax.transAxes)
        ax2.text(0.3,-0.38,'chi-sq',fontsize=label_fontsize, transform=ax.transAxes)
    #fig.tight_layout()
    #subplots_adjust(top=0.90,wspace=1.0)
    #fig.savefig('../figures/morph_chisqs.svg', bbox_inches='tight',dpi=fig.dpi)
    #fig.savefig('../figures/morph_chisqs.png', bbox_inches='tight',dpi=fig.dpi)

def plot_onemitexample_R2N_hist_paperfigure(eg_netseed,eg_mitnum,resultsdir='../results/odor_morphs'):
    """ Plot residual to noise histogram of the morph fits. """
    fig = figure(figsize=(columnwidth,columnwidth/2.0),dpi=300,facecolor='w') # 'none' is transparent
    ax3 = fig.add_subplot(2,3,1)
    ax4 = fig.add_subplot(2,3,2)
    ax5 = fig.add_subplot(2,3,4)
    ax6 = fig.add_subplot(2,3,5)
    ax1 = fig.add_subplot(2,3,3)
    ax2 = fig.add_subplot(2,3,6)
    ## inh =  (no_singles,no_joints,no_lat,no_PGs,varyRMP)
    inh_options = [ (0,(False,False,False,False,False),'lat inh') ]
    for ploti,(inhi,inh,inhstr) in enumerate(inh_options):
        R2Ns = []
        lin_R2Ns = []
        chilist = []
        n_accept = 0
        for stimi,stimseed in enumerate(stim_seeds):
            if not salient: net_seeds = [stimseed]
            for neti,netseed in enumerate(net_seeds):
                for ngi,num_gloms in enumerate([3]):

                    filename, switch_strs \
                        = get_filename(netseed,stimseed,inh,num_gloms,stimi,neti,inhi,resultsdir=resultsdir)
                    switches_str = string.join(switch_strs,'')
                    ## if the result file for these seeds & tweaks doesn't exist,
                    ## then carry on to the next.
                    if not os.path.exists(filename): continue
                    print filename
                    for fitted_mitral in [0,1]:
                        ## First the weighted-linear sigmoid:
                        ## If the fitted params file does not exist, create it (them).
                        if not os.path.exists(filename+'_params'+str(fitted_mitral)):
                            print "fitting file",filename
                            refit = True
                        else: refit = False
                        ## read in params & responses for this result file
                        mit_fit_params = \
                            fit_om.fit_morphs(filename, fitted_mitral, 'arb', refit=refit)
                        params,chisq,inputsA,inputsB,fitted_responses,\
                                numavgs,firingbinsmeanList,firingbinserrList = mit_fit_params
                        S2N,S2R = forR2N.residual2noise(fitted_responses[-2],firingbinsmeanList[-2],\
                            firingbinserrList[-2]*sqrt(numavgs),starti=0) # odor A
                        R2N_A = S2N/S2R
                        if isnan(R2N_A): continue
                        S2N,S2R = forR2N.residual2noise(fitted_responses[0],firingbinsmeanList[0],\
                            firingbinserrList[0]*sqrt(numavgs),starti=0) # odor B
                        R2N_B = S2N/S2R
                        if isnan(R2N_B): continue
                        R2Ns.append(R2N_A)
                        R2Ns.append(R2N_B)
                        if netseed == eg_netseed and fitted_mitral == eg_mitnum:
                            fit_om.plot_example_onemit(ax3,ax4,eg_mitnum,mit_fit_params)
                        
                        ## Linear-rectifier or Linear-sigmoid depending on FULLlin variable above.
                        ## If the fitted params file does not exist, create it (them).
                        if not os.path.exists(filename+'_params'+linextn+str(fitted_mitral)):
                            print "fitting FULLlin file",filename
                            refit = True
                        else: refit = False
                        ## fit/get the params and responses for this result file
                        mit_fit_params = \
                            fit_om.fit_morphs(filename, fitted_mitral, 'lin', refit=refit)
                        params,chisq,inputsA,inputsB,fitted_responses,\
                            numavgs,firingbinsmeanList,firingbinserrList = mit_fit_params
                        S2N,S2R = forR2N.residual2noise(fitted_responses[-2],firingbinsmeanList[-2],\
                            firingbinserrList[-2]*sqrt(numavgs),starti=0) # odor A
                        R2N_A = S2N/S2R
                        if isnan(R2N_A): continue
                        S2N,S2R = forR2N.residual2noise(fitted_responses[0],firingbinsmeanList[0],\
                            firingbinserrList[0]*sqrt(numavgs),starti=0) # odor B
                        R2N_B = S2N/S2R
                        if isnan(R2N_B): continue
                        lin_R2Ns.append(R2N_A)
                        lin_R2Ns.append(R2N_B)
                        chilist.append(sqrt(chisq))
                        if netseed == eg_netseed and fitted_mitral == eg_mitnum:
                            fit_om.plot_example_onemit(ax5,ax6,eg_mitnum,mit_fit_params)

                        n_accept += 1

        R2N_max = 1.0
        ax1.hist(clip(R2Ns,0,R2N_max),20,normed=True,edgecolor='b',facecolor='b')
        _,y1 = ax1.get_ylim()
        ax2.hist(clip(lin_R2Ns,0,R2N_max),20,normed=True,edgecolor='b',facecolor='b')
        #ax2.hist(clip(chilist,0,R2N_max),20,normed=True,edgecolor='b',facecolor='b')
        _,y2 = ax2.get_ylim()
        yR2Nmax = max(y1,y2)
        print "Number of mitral cells accepted =",n_accept
        
        ## beautify plots
        for axnum,ax in enumerate([ax1,ax2]):
            xmin,xmax,ymin,ymax = \
                beautify_plot(ax,x0min=True,y0min=True,xticksposn='bottom',yticksposn='left')
            ax.set_xlim([0,R2N_max])
            ax.set_xticks([0,R2N_max])
            ax.set_ylim([0,yR2Nmax])
            ax.set_yticks([0,yR2Nmax])
        for ax in [ax1,ax3,ax4]:
            ax.set_xticklabels(['',''])
        ## axes_labels() sets sizes of tick labels too.
        axes_labels(ax1,'','prob. density',adjustpos=False,xpad=0,ypad=0)
        ax1.yaxis.set_label_coords(-0.29,-0.3)
        axes_labels(ax2,'$\sqrt{residual/noise}$','',adjustpos=False,xpad=1,ypad=0)

        axes_labels(ax3,'','firing rate (Hz)',adjustpos=False,xpad=0,ypad=0)
        ax3.yaxis.set_label_coords(-0.29,-0.3)
        axes_labels(ax5,'time (s)','',adjustpos=False,xpad=3,ypad=0)

        axes_labels(ax4,'','fitted weight',adjustpos=False,xpad=0,ypad=0)
        ax4.yaxis.set_label_coords(-0.24,-0.3)
        axes_labels(ax6,'conc (% SV)','',adjustpos=False,xpad=3,ypad=0)

        fig_clip_off(fig)
        fig.tight_layout()
        fig.subplots_adjust(hspace=0.3,wspace=0.5) # has to be after tight_layout()
        fig.savefig('../figures/morphs_R2Ns.svg',dpi=fig.dpi)
        fig.savefig('../figures/morphs_R2Ns.png',dpi=fig.dpi)

if __name__ == "__main__":
    if len(stim_seeds)<5:
        plot_responses()
        #plot_xcorrgrams()
        #plot_decorr_single()
    else:
        ## default network: set directed=True, frac_directed=0.05, varmit False (in inh_options)
        #glomnums = range(1,6)
        ## non-decorrelating networks:
        ## set directed=True, frac_directed=0.0, varmit True/False (in inh_options)
        ## OR directed=False; varmit False (in inh_options)
        #glomnums = [2,3,6]
        #plot_directed(glomnums)
        
        ## PAPER figure 7: plot the correlation distributions for various network connectivities
        ## CAUTION: SET stim_seeds AT THE TOP to (750.0,1100.0)
        #plot_across_sims_paperfigure()
        
        ###### for the default network, plot the best decorr-ed responses:
        ###### PAPER figure 6: for distribution of air and odor distributions & example.
        ## CAUTION: SET stim_seeds AT THE TOP to (750.0,1100.0)
        if len(sys.argv)>1: resultsdir = sys.argv[1]
        else: resultsdir = '../results/odor_morphs/'
        print 'Using results directory:',resultsdir
        #plot_decorrs_special_paperfigure(resultsdir,[3],[(0,(False,False,False,False,False))],\
        #    _directed=True,_frac_directed=0.01,graph=True)
        ## PAPER figure 6: example decorr: choose odornum as 5 for odorA and 0 for odor B
        ## CAUTION: SET stim_seeds AT THE TOP to (750.0,1100.0)
        #plot_responses_mits_paperfigure( resultsdir,odornum=0,stimseed=844.0,\
        #    numgloms=3,inh=(False,False,False,False,False) )
        
        ############ ODOR MORPHS -- Adil style fits
        ## OBSOLETE -- odor morphs chisq histogram
        #plot_chisq_hist()
        
        ## PAPER FIGURE supplementary figure 5 : Adil's morph fits
        ## CAUTION: SET stim_seeds AT THE TOP to (750.0,800.0)
        ## plots one mitral morphs fits example
        ## AND histogram of overall fits
        ## pass example netseed, example mitnum and resultsdir
        #plot_onemitexample_R2N_hist_paperfigure(754.0,0,resultsdir) # 752.0,0 is an alternative eg

        ## PAPER FIGURE supplementary figure 6: tufted (no PG) vs mitral (with PG) phase difference
        ## JUST STARTED, INCOMPLETE
        plot_peaks_tufted_vs_mitrals_paper_figure()

    show()

Loading data, please wait...