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 -*-

########## THIS FITTING PROGRAM IS MEANT TO ROUGHLY FOLLOW PRIYANKA'S ANALYSIS
########## This variant does not use an air kernel, rather a constant air rate.
## USAGE1: python2.6 fit_odor_pulses.py <../results/odor_pulses/pulseresult.pickle> 191.0
## USAGE2: python2.6 fit_odor_pulses.py <../results/odor_pulses/pulseresult.pickle> 191.0 <../results/odor_morphs/morphresult.pickle>

########## OBSOLETE ######### OBSOLETE ########## OBSOLETE ######### OBSOLETE ######## OBSOLETE ########
####### Use fit_odor_pulses_constair_separateodors.py instead. Priyanka fits each odor separately for each mitral.

from scipy import optimize
from scipy import special
import scipy.interpolate
from scipy import signal
import scipy.io
from pylab import *
import pickle
import sys
import math

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

from stimuliConstants import * # has SETTLETIME, PULSE_RUNTIME
from networkConstants import * # has central_glom
from data_utils import * # has axes_labels()
from analysis_utils import * # has read_morphfile() and various constants
import generate_firerates_odors as gen_frate_odors # has respiration_function()

iterationnum = 1
sigmoid_op = False
SMOOTH_KERNELS = False

## Have overridden the kernel_time=1.5s (for generating ORN frates) in stimuliConstants.py,
## as Priyanka uses 2.0s kernels.
kernel_time = 2.0

numpulses = RANDOM_PULSE_NUMS*2
## rebin the spikes as below, irrespective of previous binning
pulserebindt = fitting_dt # 50ms as Priyanka uses
pulserebins = int(PULSE_RUNTIME/pulserebindt)
bin_width_time = 2*pulserebindt ## overlapping bins with this bin width
pulsetlist = arange(pulserebindt/2.0,PULSE_RUNTIME,pulserebindt)
## very IMPORTANT to have the same dt for pulse and kernel
## when doing discrete convolution!
kerneldt = pulserebindt # same as fitting_dt above
kernel_size = int(kernel_time/kerneldt)
kerneltlist = arange(0.0,kernel_time,kerneldt)[:kernel_size]
#NUMWTS = NUMMIX-3 # remove the two pure odors and one pure air weights
#NUMPARAMS = 3*NUMBINS+2*NUMWTS+1 # two pure + one air response,
# weights of A and B for mixtures, max of output sigmoid
firstrun = True

## If you want to plot in style of Priyanka
fill_below = True

## number of the mitral to be fitted.
fitted_mitrals = [2*central_glom,2*central_glom+1]

global iternum
iternum = 0

resp_pulses = gen_frate_odors.respirationPulses
st_idx = gen_frate_odors.startidx
numt = gen_frate_odors.numt

def responseconvolve(pulses, kernel):
    ## don't use mode='same', as it takes only the valid part of the convolution.
    ## always multiply convolution_result by dt
    ## to leave response invariant on decimating stim and kernel simultaneously.
    responsecalc = convolve( pulses, kernel, mode='full' ) * pulserebindt
    ## bin0 of the simulated response is at time pulserebindt
    ## so the fitted response should also be shifted by a pulserebindt
    return array( responsecalc[1:pulserebins+1] )
    
def sigmoid_fn(valarray, max, half_point):
    ## ensure that valarray is a numpy array
    ## else where() returns empty tuple, and below exp() statement won't work.
    valarray = array(valarray)
    if sigmoid_op:
        ## exp() - logistic function - not used
        #return max/(1.0+exp(-valarray+half_point)) # numpy's exp works for arrays.
        ## thresholded erf() - used
        ## erf(0.5)~=0.5; set half_point as center of sigmoid.
        valarray[where(valarray<0)[0]] = 0
        posindices = where(valarray<0)[0]
        valarray[posindices] = max * special.erf( 0.5*valarray[posindices]/half_point )
        return valarray
    else:
        ## where value in valarray < 0, set it to zero.
        valarray[where(valarray<0)[0]]=0
        return valarray

#def inv_sigmoid_fn(valarray, max, half_point):
#    ## ensure that valarray is a numpy array
#    ## else where() returns empty tuple, and below exp() statement won't work.
#    valarray = array(valarray)
#    if sigmoid_op:
#        ## numpy has overloaded operators, also numpy ln is overloaded for arrays.
#        return half_point - log(max/valarray-1)
#    else:
#        return valarray

def SGfilter2(kernelA,kernelB):
    ## savitsky-golay filter maintains high frequency signal, unlike moving average filter
    ## window_size must be odd, savitsky_golay copied from scipy cookbook online
    window_size = 11 # same as Priyanka
    kernelA = savitzky_golay(kernelA, window_size, order=4)
    kernelB = savitzky_golay(kernelB, window_size, order=4)
    return kernelA, kernelB

def chisqfunc(params, firingbinsmeanList, firingbinserrList, bgnd, \
        pulserebinList, ret_array=True, ORNkernels=()):
    kernelA = params[0:kernel_size]
    kernelB = params[kernel_size:2*kernel_size]
    if sigmoid_op:
        sigmoid_max = params[-2]
        sigmoid_threshold = params[-1]
    else:
        sigmoid_max = None
        sigmoid_threshold = None
    if SMOOTH_KERNELS:
        kernelA, kernelB = SGfilter2(kernelA, kernelB)
    ## These below should be numpy arrays,
    ## else + acts like append in python (not element-wise addition)!!!
    responsecalcA1 = sigmoid_fn(
        responseconvolve(pulserebinList[2],kernelA) + bgnd,\
        sigmoid_max, sigmoid_threshold)
    responsecalcB1 = sigmoid_fn(
        responseconvolve(pulserebinList[3],kernelB) + bgnd,\
        sigmoid_max, sigmoid_threshold)
    responsecalcA2 = sigmoid_fn(
        responseconvolve(pulserebinList[4],kernelA) + bgnd,\
        sigmoid_max, sigmoid_threshold)
    responsecalcB2 = sigmoid_fn(
        responseconvolve(pulserebinList[5],kernelB) + bgnd,\
        sigmoid_max, sigmoid_threshold)
    ### don't use A+B to fit, just plot the result on A+B.
    #responsecalcAB = sigmoid_fn(
    #    responseconvolve(pulserebinList[6],kernelA) + responseconvolve(pulserebinList[7],kernelB) + \
    #    + bgnd, sigmoid_max, sigmoid_threshold)
    ## fitted response is also shifted by a pulserebindt, to compare index to index
    ## only compare after odor onset at PULSE_START
    starti = int(PULSE_START/pulserebindt)
    chisqarray = [ (val-firingbinsmeanList[1][i])/firingbinserrList[1][i]\
        for i,val in enumerate(responsecalcA1[starti:],start=starti) ]
    chisqarray.extend( [ (val-firingbinsmeanList[2][i])/firingbinserrList[2][i]\
        for i,val in enumerate(responsecalcB1[starti:],start=starti) ] )
    chisqarray.extend( [ (val-firingbinsmeanList[3][i])/firingbinserrList[3][i]\
        for i,val in enumerate(responsecalcA2[starti:],start=starti) ] )
    chisqarray.extend( [ (val-firingbinsmeanList[4][i])/firingbinserrList[4][i]\
        for i,val in enumerate(responsecalcB2[starti:],start=starti) ] )
    #chisqarray.extend( [ (val-firingbinsmeanList[5][i])/firingbinserrList[5][i]\
    #    for i,val in enumerate(responsecalcAB[starti:],start=starti) ] )
        
    global iternum
    iternum += 1
    ## normalize chi-sq to the number of dof
    num_dof = float(firingbinsmeanList[1][starti:].size*(numpulses-2) - params.size)
    if ret_array:
        if iternum % 10000 == 0:
            chisqarraysq = [i**2 for i in chisqarray]
            chisq = reduce(lambda x, y: x+y, chisqarraysq) / num_dof
            print chisq
        ## len(chisqarray) must be >= len(params):
        ## since leastsq() assumes it is solving M eqns in N unknowns where M>=N.
        return array(chisqarray) / sqrt(num_dof) # this must be squared and added for chi-sq
    else:
        chisq = reduce(lambda x, y: x+y**2, chisqarray) / num_dof
        if iternum % 10000 == 0:
            print chisq
        return chisq

def fit_pulses(filename,stimseed,match_morph,NOSHOW):
    ######### load in the mitral pulse responses
    f = open(filename,'r')
    #### mitral_responses_list[avgnum][pulsenum][mitralnum][binnum]
    #### mitral_responses_binned_list[avgnum][pulsenum][mitralnum][binnum]
    mitral_responses_list, mitral_responses_binned_list = pickle.load(f)
    f.close()
    ########## load in the stimuli
    seednum = stimseed
    f = open('../generators/firerates/firerates_2sgm_'+str(seednum)+'.pickle','r')
    frateOdorList,fratePulseList,randomPulseList,\
        randomPulseStepsList,randomResponseList,kernels\
        = pickle.load(f)
    #del frateOdorList, fratePulseList
    pulseList = array(randomPulseList)
    pulseStepsList = array(randomPulseStepsList)
    ## randomResponseList[glomnum][frate,...]
    ORNfrateList = array(randomResponseList)[central_glom]
    f.close()

    ##-------------------------- rebin the responses and pulses ------------------------------
    ### for overlapping pulses (uncomment in analysis_utils.py)
    #mitral_responses_binned_list = \
    #    rebin_pulses(mitral_responses_list, int(PULSE_RUNTIME/pulserebindt), PULSE_RUNTIME, 0.0, bin_width_time)
    mitral_responses_binned_list = \
        rebin_pulses(mitral_responses_list, int(PULSE_RUNTIME/pulserebindt), PULSE_RUNTIME, 0.0)
    ## don't just decimate, rather take average of each bin
    decimatefactor = int(round(pulserebindt/FIRINGFILLDT)) # ensure that these are multiples
    ## binning the pulses: average of every decimatefactor number of values
    pulserebinList = array( [ array([ pulseList[pulsenum][i:i+decimatefactor].mean() \
            for i in range(0,len(pulseList[pulsenum]),decimatefactor) ]) \
        for pulsenum in range(len(pulseList)) ] )
    ## decimating the pulses, i.e. taking one value after every decimatefactor
    pulserebinList = array( [ pulseList[pulsenum][::decimatefactor] \
        for pulsenum in range(len(pulseList)) ] )

    ## very important to convert to numpy array(),
    ## else where() below does not find any element satisfying condition.
    mitral_responses_binned_list = array(mitral_responses_binned_list)
    numtrials = len(mitral_responses_binned_list)
    mitral_responses_mean = mean(mitral_responses_binned_list, axis=0)
    mitral_responses_std = std(mitral_responses_binned_list, axis=0)
    ## since I fit the mean response, I must use standard error/deviation of the _mean_
    ## = standard deviation of a repeat / sqrt(num of repeats).
    ## NO! The model is for an individual response, not for the mean response!
    #mitral_responses_se = mitral_responses_std/sqrt(numtrials)

    fittedkernels = []
    chisqs = []
    global sigmoid_op
    ## full fitting data for both mitrals
    mdict_full = []
    fitted_responses_full = []
    if not NOSHOW:
        fig = figure(facecolor='w') # 'none' is transparent
        ax_kernelsA = plt.subplot2grid((6,3),(0,2),rowspan=3)
        ax_kernelsB = plt.subplot2grid((6,3),(3,2),rowspan=3)
    for fitted_mitral in fitted_mitrals:
        ## take the odor responses of the mitral to be fitted
        firingbinsmeanList = mitral_responses_mean[:,fitted_mitral]
        firingbinserrList = mitral_responses_std[:,fitted_mitral]
        ## slightly different from Priyanka -- not taking full air buffer, allowing to settle
        bgnd = mean( [response[int((SETTLETIME+PULSE_AIR_BUFFER/2.0)/pulserebindt):\
                int((SETTLETIME+PULSE_AIR_BUFFER)/pulserebindt)] \
            for response in firingbinsmeanList] )
        print "bgnd =",bgnd
        ### force the errors to be 1.2*sqrt(mean) which is what Adil observes.
        ### The model seems to have constant noise independent of the mean firing rate.
        ### standard error of the mean, hence /sqrt(numtrials)
        #firingbinserrList = 1.2*firingbinsmeanList/sqrt(numtrials)

        errcut = 5e-1
        #errcut = 1.0/pulserebindt/sqrt(numtrials) # errcut = 3.33Hz!
        #errcut = 10.0 #Hz
        ## put in a minimum error, else divide by zero problems.
        ## find the minimum error >= errcut
        largeerrors = firingbinserrList[where(firingbinserrList>errcut)]
        if largeerrors!=[]: errmin = largeerrors.min()
        else: errmin = errcut
        print firingbinserrList
        print 'least error = ',errmin
        errmin = errcut
        ## numpy where(), replace by errmin,
        ## all those elements in firingbinsList which are less than errmin 
        firingbinserrList = where(firingbinserrList>errcut, firingbinserrList, errmin)

        #---------------------------- fitting ------------------------------------------------

        decimatefactor = int(kerneldt/FIRINGFILLDT)
        ## kernels are only kernel_size long.
        ## numpy append() is like python's extend() - makes flat not nested lists
        ## I take only kernelsize length after decimating the ORN kernel 
        ## by taking only every int(kerneldt/FIRINGFILLDT) element
        ## since convolve_result is always multiplied by dt,
        ## the response is invariant of decimating the stimulus and the kernel simultaneously.
        ## NOTE: We are using ORN kernel as mitral kernel, typical mitral
        ## firing rates are 5 times larger. So put in factor of 5!!!
        ## If ORN kernel is too long, take only kernel_size, if too short, pad with zeroes.
        MitralORNfiringratio = 5.0
        linORNkernelA = array(kernels[1][::decimatefactor][:kernel_size])*MitralORNfiringratio
        if len(linORNkernelA)<kernel_size:
            linORNkernelA = append(linORNkernelA,zeros(kernel_size-len(linORNkernelA)))
        linORNkernelB = array(kernels[2][::decimatefactor][:kernel_size])*MitralORNfiringratio
        if len(linORNkernelB)<kernel_size:
            linORNkernelB = append(linORNkernelB,zeros(kernel_size-len(linORNkernelB)))
        ## Initial values for the parameters
        if firstrun:
            ORNkernels = (linORNkernelA, linORNkernelB)
            params = linORNkernelA # ORN kernelA is initial kernelA
            params = append(params,linORNkernelB) # ORN kernelB is initial kernelB
            if sigmoid_op:
                ## I don't have a steepness param,\
                ## else there will be redundancy in kernel scaling vs steepness.
                params = append(params, [120.0,60.0]) # 120Hz max and 60Hz threshold for sigmoid
        else:
            f = open(sys.argv[1]+'_params'+str(fitted_mitral),'r')
            sigmoid_op,chisq,params,discard_firingmeans,discard_firingsd,discard_fitresponses = pickle.load(f)
            f.close()
            ORNkernels = ()

        ## export the data for use in matlab in Priyanka's fitting function
        mdict={'firingbinsmeanList':firingbinsmeanList,
            'firingbinserrList':firingbinserrList, 'bgnd':bgnd, 'FIRINGFILLDT':FIRINGFILLDT,
            'pulseList':pulserebinList,'pulseStepsList':pulseStepsList,
            'start_kernels':ORNkernels,'pulserebindt':pulserebindt,'pulsetlist':pulsetlist}
        scipy.io.savemat('mit'+str(fitted_mitral)+'.mat', mdict=mdict)
        mdict_full.append(mdict)

        ## comment next fitting line if you just want to plot parameters.
        ## args is a tuple! if only one element write (elem, )

        #########
        ### fmin_powell is faster than fmin (simplex) but slower than leastsq
        ### However it does not enter into the very bad local minima of leastsq
        ### WARNING! But often it is worse at fitting than leastsq()
        #iternum = 0
        #params = optimize.fmin_powell( chisqfunc, params, \
        #    args=(firingbinsmeanList, firingbinserrList, pulserebinList, False),\
        #    full_output=1, maxfun=50000, maxiter=50000, ftol = 1 )
        #chisq = params[1]
        #print "The normalized chisq for mitral",fitted_mitral,"after fmin_powell =",chisq
        #params = params[0]

        #########
        ## leastsq() uses a modified version of Levenberg-Marquardt algorithm
        ## it optimizes M equations in N unknowns, where M>=N.
        ## Hence chisqfunc must return M numbers in an array
        ## which must be >= the number of params N in params0, 
        ## else: 'TypeError: Improper input parameters.'
        iternum = 0
        params = optimize.leastsq( chisqfunc, params,\
            args=(firingbinsmeanList, firingbinserrList, bgnd, pulserebinList, True, ORNkernels),\
            full_output=1, maxfev=100000, ftol=1e-100, xtol=1e-100 )
        print params[3] # print the status message
        params = params[0] # take only fitted params; leastsq returns a tuple with errmsg, etc.
        ### Calculate sum of squares of the chisqarray
        chisqarraysq = [i**2 for i in chisqfunc(params,\
            firingbinsmeanList, firingbinserrList, bgnd, pulserebinList, True, ORNkernels)]
        #chisqnumpy = array(chisqarraysq)
        #pos = where(chisqnumpy>100)[0]
        #lenf = len(firingbinsmeanList[0])
        #print "The [(pulsenum,index),...] where chisq>100 are",[ (i/lenf,i%lenf) for i in pos]
        chisq = reduce(lambda x, y: x+y, chisqarraysq)
        print "The normalized chisq for mitral",fitted_mitral,"after leastsq =",chisq
        
        ##########
        ### minimize is a wrapper for all the various fitting algorithms!
        #iternum = 0
        #fit_method = "Nelder-Mead"#"CG"#"Powell"#"Nelder-Mead"
        #params,info = optimize.minimize( chisqfunc, params,\
        #    args=(firingbinsmeanList, firingbinserrList, pulserebinList, False, ORNkernels),\
        #    method = fit_method, full_output=True,\
        #    options={'maxfev':100000, 'ftol':1e-100, 'xtol':1e-100, 'maxfun':500000, 'maxiter':500000} )
        #print info['message']
        #chisq = info['fun']
        #print "The normalized chisq for mitral",fitted_mitral,"after",fit_method,"=",chisq

        kernelA = params[0:kernel_size]
        kernelB = params[kernel_size:2*kernel_size]
        ## irrespective of smooth kernels flag, filter the output kernels for display
        kernelA, kernelB = SGfilter2(kernelA, kernelB)

        if sigmoid_op:
            sigmoid_max = params[-2]
            sigmoid_threshold = params[-1]
        else:
            sigmoid_max = None
            sigmoid_threshold = None

        #chisq = chisqfunc(params,firingbinsmeanList, firingbinserrList, pulserebinList)

        fittedkernels.append( (kernelA,kernelB) )
        chisqs.append(chisq)

        fitted_responses = []
        for pulsenum in range(1,6):
            if pulsenum in [1,3]:
                response = sigmoid_fn(responseconvolve(pulserebinList[pulsenum+1],kernelA) + bgnd,\
                    sigmoid_max, sigmoid_threshold)
            elif pulsenum in [2,4]:
                response = sigmoid_fn(responseconvolve(pulserebinList[pulsenum+1],kernelB) + bgnd,\
                    sigmoid_max, sigmoid_threshold)
            else:
                response = responseconvolve(pulserebinList[pulsenum+1],kernelA)
                response += responseconvolve(pulserebinList[pulsenum+2],kernelB)
                response += bgnd
                response = sigmoid_fn(response, sigmoid_max, sigmoid_threshold)
            fitted_responses.append(response)
        fitted_responses_full.append(fitted_responses)

        paramsfile = open(filename+'_params'+str(fitted_mitral),'w')
        pickle.dump((sigmoid_op,chisq,params,firingbinsmeanList,firingbinserrList,fitted_responses), paramsfile)
        paramsfile.close()


        ##---------------------------- kernel via the m-sequence method -----------------------------------------
        ## Doesn't work! One reason is that it doesn't take care of thresholding of neural response!

        ## only compare for PULSE_DURATION from PULSE_START
        starti = int(PULSE_START/pulserebindt)
        endi = int((PULSE_START+PULSE_DURATION)/pulserebindt)
        stimulus = pulserebinList[2][starti:endi]
        response = firingbinsmeanList[1][starti:endi] # response to odor A first pulse
        kernelA_mseq0 = circular_correlate(stimulus,response) # from data_utils.py
        stimulus = pulserebinList[4][starti:endi]
        response = firingbinsmeanList[3][starti:endi] # response to odor A first pulse
        kernelA_mseq1 = circular_correlate(stimulus,response) # from data_utils.py
        kernelA_mseq = (kernelA_mseq0+kernelA_mseq1)/2.0
        #figure()
        #pulsepoints = arange(0,PULSE_DURATION,pulserebindt)
        #plot(pulsepoints,kernelA_mseq)
        #plot(pulsepoints,stimulus)
        #plot(pulsepoints,response)


        #---------------------------- plot kernels and pulse responses -----------------------------------------

        if not NOSHOW:
            ## fig defined before the mitnum loop
            ############################## plot the kernels
            text(0.6,0.3,'arb units', fontsize=label_fontsize, rotation='vertical')

            ax_kernelsA.plot(kerneltlist, kernelA, color=(1-fitted_mitral,fitted_mitral,0), marker=',',\
                linestyle='solid',linewidth=linewidth,label='mit '+str(fitted_mitral))
            #ax.plot(kerneltlist, linORNkernelA, color=(1,0,1), marker=',',\
            #    linestyle='solid',linewidth=linewidth,label='5 * ORN kernel')
            ax_kernelsB.plot(kerneltlist, kernelB, color=(fitted_mitral,1-fitted_mitral,0), marker=',',\
                linestyle='solid',linewidth=linewidth,label='mit '+str(fitted_mitral))
            #ax.plot(kerneltlist, linORNkernelB, color=(0,1,1), marker=',',\
            #    linestyle='solid',linewidth=linewidth,label='5 * ORN kernel')
            ax_kernelsA.set_yticks([0])
            ax_kernelsB.set_yticklabels(['0'])
            ax_kernelsA.set_xlim(0,kernel_time)
            ax_kernelsB.set_xlim(0,kernel_time)
            biglegend(ax=ax_kernelsA)
            biglegend(ax=ax_kernelsB)
            ax_kernelsB.set_xticks([0,2.0])
            ax_kernelsB.set_xticklabels(['0','2.0'])
            axes_off(ax_kernelsA,x=True,y=False)
            axes_labels(ax_kernelsA,'','',adjustpos=False)
            axes_labels(ax_kernelsB,'time (s)','',adjustpos=False)

            ############################### plot the responses and the fits
            text(-0.125,0.25,'firing rate (Hz)', fontsize=label_fontsize, rotation='vertical')
            for pulseiter, pulsenum in enumerate([1,2,5]):#range(1,numpulses): # similar to Priyanka, only 3 plots
                sister_ratio = 0#(mitnum%MIT_SISTERS)/float(MIT_SISTERS)
                ## similar to Priyanka: only 3 plots, skip 2 pulse responses
                ax = plt.subplot2grid((6,3),(2*pulseiter,fitted_mitral),rowspan=2)
                ################# random pulses and ORN responses
                xpulse = int(max(firingbinsmeanList[pulsenum]+firingbinserrList[pulsenum]/sqrt(9))+1)
                randompulsetime = arange(0,PULSE_RUNTIME+1e-10,FIRINGFILLDT)
                if pulsenum in [1,2,3,4]: # odor A or B
                    if fill_below:
                        if pulsenum in [1,3]: col = (1,0,0)
                        if pulsenum in [2,4]: col = (0,1,0)
                        fill_between(randompulsetime,xpulse*pulseStepsList[pulsenum+1]+
                            pulseList[0]+ORN_BGND_MEAN,linewidth=0,color=col,alpha=0.4)
                    else:
                        plot(randompulsetime,2*pulseList[pulsenum+1]+pulseList[0]+ORN_BGND_MEAN,
                            linewidth=linewidth,label='Air+Odor waveform')
                        plot(randompulsetime,ORNfrateList[pulsenum+1]+ORNfrateList[0]+ORN_BGND_MEAN,
                            linewidth=linewidth,label='Receptors response')
                elif pulsenum in [5]: # odor A & odor B
                    if fill_below:
                        fill_between(randompulsetime,xpulse*pulseStepsList[pulsenum+1],\
                            color=(1,0,0),linewidth=0,alpha=0.4)
                        fill_between(randompulsetime,xpulse*pulseStepsList[pulsenum+2],\
                            color=(0,1,0),linewidth=0,alpha=0.4)
                    else:
                        plot(randompulsetime,pulseList[0]+ORN_BGND_MEAN,\
                            linewidth=linewidth, label='Air waveform')
                        plot(randompulsetime,2*pulseList[pulsenum+1],\
                            linewidth=linewidth, label='Odor A waveform')
                        plot(randompulsetime,2*pulseList[pulsenum+2],\
                            linewidth=linewidth, label='Odor B waveform')
                        plot(randompulsetime,ORNfrateList[pulsenum+2]+ORNfrateList[pulsenum+1]\
                            +ORNfrateList[0]+ORN_BGND_MEAN, linewidth=linewidth, label='Receptors response')
                ################### Plot the simulated responses
                ## smooth the simulated response
                ## windowsize=5 and SD=0.65 are defaults from matlab's smoothts() for gaussian smoothing
                Gwindow = signal.gaussian(5,0.65)
                ## help from http://www.scipy.org/Cookbook/SignalSmooth
                simresponse = convolve(Gwindow/Gwindow.sum(),firingbinsmeanList[pulsenum],mode='same')
                if fill_below:
                    ## numpy array, hence adds element by element
                    fill_between(pulsetlist,
                        simresponse+firingbinserrList[pulsenum]/sqrt(9),
                        simresponse-firingbinserrList[pulsenum]/sqrt(9),
                        color=(0.7,0.7,0.7))
                    plot(pulsetlist,simresponse,linewidth=linewidth,color=(0.3,0.3,0.3))
                else:
                    errorbar(pulsetlist,y=simresponse,\
                        yerr=firingbinserrList[pulsenum]/sqrt(9),\
                        color=(pulsenum/float(numpulses),1-pulsenum/float(numpulses),sister_ratio),\
                        marker='+',linestyle='solid',linewidth=linewidth,label='Mitral response')
                ################## Plot the fitted responses
                starti = int(PULSE_START/pulserebindt)
                if pulsenum in [1,3]: titlestr = 'Odor A random pulse'
                elif pulsenum in [2,4]: titlestr = 'Odor B random pulse'
                else: titlestr = 'Odor A and odor B random pulse'
                response = fitted_responses[pulsenum-1]
                if fill_below:
                    plot(pulsetlist[starti:],response[starti:],linestyle='solid',
                        linewidth=linewidth,color=(0,0,0),
                        label=['air','A','B','A','B','A&B'][pulsenum])
                    lgd = legend()
                    for k in lgd.get_lines():
                        k.set_linewidth(0)
                    lgd.draw_frame(False)
                    ltext  = lgd.get_texts()
                    for l in ltext:
                        l.set_fontsize(label_fontsize)
                    #fig.setp(ltext, fontsize=20)
                else:
                    plot(pulsetlist[starti:],response[starti:],
                        marker='o',linestyle='dashed',
                        linewidth=linewidth, label='Linear fit')
                    #title('Odor morphs w/ smooth fit',fontsize=24 )
                    #title( 'pulsenum = '+str(pulsenum)+', chisquare normalized = '+str(chisq) )
                    title(titlestr, fontsize=label_fontsize)
                    lgd = legend()
                    ltext  = lgd.get_texts()
                    for l in ltext:
                        l.set_fontsize(label_fontsize)
                    ax.set_ylim(0,xpulse)
                ax.set_yticks([xpulse])
                ax.set_yticklabels([str(xpulse)])
                for label in ax.get_yticklabels():
                    label.set_fontsize(label_fontsize)
                ax.set_xlim(0,9)
                ax.set_xticks([])
            ax.set_xticks([0,9])
            ax.set_xticklabels(['0','9'])    
            axes_labels(ax,'time (s)','',adjustpos=False)

    if not NOSHOW:
        show()
    
    return fittedkernels,chisqs,fitted_responses_full,mdict_full

if __name__ == "__main__":
    if len(sys.argv) > 2:
        ######### should we match mitral morph responses using kernel
        if '--match_morph' in sys.argv:
            match_morph = True
            mmfilenameindex = sys.argv.index('--match_morph')+1
            morph_filename = sys.argv[mmfilenameindex]
        else: match_morph = False
        if 'NOSHOW' in sys.argv: NOSHOW = True
        else: NOSHOW = False
        fit_pulses(sys.argv[1],sys.argv[2],match_morph,NOSHOW)
    else:
        print "At least specify data file containing pickled mitral responses, and ORN frate seed."
        sys.exit(1)

Loading data, please wait...