# -*- coding: utf-8 -*-
# ALL SI UNITS
# milliMolar is same as mol/m^3

## USAGE: python2.6 average_odor_scaledpulses_noregress.py

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 scipy import signal # for gaussian smoothing
from scipy import interpolate # for interp1d
from scipy import optimize # for fsolve
from sim_utils import * # has rebin() and imports data_utils.py for axes_off()
from analysis_utils import * # has rebin_mean() and read_pulsefile()
from calc_corrs import * # has calc_corrs()
## BE CAREFUL: use _constair_separateodors below, others are obsolete.
import fit_scaledpulses_noregress as fit_sp # has fit_plot_scaledpulses and ref pulse constants
import average_odor_morphs as corr_utils # has get_filename() for morphs

IN_VIVO = True
directed = True
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
_scaledWidth = 0.2 # s # overrides that in stimuliConstantsMinimal.py
num_scalings = len(pulseList)-1

#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
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 )

min_frate_cutoff = 0.25 # Hz

salient = False#True
if salient:
    stim_seeds = [-1,-19]#range(-1,-37,-1)#[-1,-2,-3,-4,-5,-6,-7,-8,-9]
    num_gloms_list = [3]
    ## 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)), (1,(False,False,True,False,False)) ]
else:
    #stim_seeds = [157.0,160.0,190.0,191.0,212.0,441.0]
    #num_gloms_list = [5,2]
    stim_seeds = arange(750.0,800.0,1.0)#[157.0,160.0,190.0,191.0]
    num_gloms_list = [3]
    ## 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,(True,False,False,False,False)), (3,(True,True,False,False,True)), \
        (6,(False,False,False,True,False)), (8,(True,True,False,True,False))]
net_seeds = [200.0]

def get_filename(netseed,stimseed,inh,ngi,\
        nl_orns=NONLINEAR_ORNS,nl_type=NONLINEAR_TYPE,\
        resultsdir='../results/odor_pulses'):
    ### read filename from the output file of automated run
    ## 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+'/scaledpulses_width'+str(_scaledWidth)+\
        '_netseed'+str(netseed)+'_stimseed'+str(stimseed)
    if nl_orns: filename += '_NL'+nl_type
    filename += singles_str+joints_str+pgs_str+lat_str+varmitstr+\
        '_numgloms'+str(num_gloms_list[ngi])
    if directed: filename += '_directed'+str(FRAC_DIRECTED)
    filename += '.pickle'
    return filename, (singles_str, joints_str, pgs_str, lat_str, varmitstr)

def plot_lin_contribs_oldpaperfigure():
    """ plot linearity scores for all options
    """
    fig1 = figure(figsize=(columnwidth,linfig_height),dpi=300,facecolor='w')
    ax1a = plt.subplot2grid((2,5),(0,0)) # full
    ax1b = plt.subplot2grid((2,5),(1,0)) # full
    ax2a = plt.subplot2grid((2,5),(0,1)) # no lat
    ax2b = plt.subplot2grid((2,5),(1,1)) # no lat
    ax3a = plt.subplot2grid((2,5),(0,2)) # no self
    ax3b = plt.subplot2grid((2,5),(1,2)) # no self
    ax4a = plt.subplot2grid((2,5),(0,3)) # no inh
    ax4b = plt.subplot2grid((2,5),(1,3)) # no inh
    ax5a = plt.subplot2grid((2,5),(0,4)) # non-lin
    ax5b = plt.subplot2grid((2,5),(1,4))
    ## inh = (no_singles,no_joints,no_lat,no_PGs,varyRMP)
    inh_options = [
        ('',0,(False,False,False,False,False),'all',False,None,(ax1a,ax1b)), \
        ('',1,(False,False,True,False,False),'no lat-inh',False,None,(ax2a,ax2b)), \
        ('',2,(True,False,False,True,False),'no self-inh',False,None,(ax3a,ax3b)), \
        ('',3,(True,True,False,True,False),'no inh',False,None,(ax4a,ax4b)), \
        ('',4,(False,False,False,False,False),'non-lin ORNs',True,'P',(ax5a,ax5b)) ]
    Rsqs_all = [ [] for i in range(len(inh_options)) ]
    R_all = [ [] for i in range(len(inh_options)) ]
    maxy_hist = 0.0
    for ploti,(dirextn,inhi,inh,inhstr,nl_orns,nl_type,(axa,axb)) in enumerate(inh_options):
        slopes_all = []
        peaks_all = []
        avg_frates_all = []
        for stimi,stimseed in enumerate(stim_seeds):
            filename, switch_strs \
                = get_filename(stimseed,stimseed,inh,0,nl_orns,nl_type,\
                    resultsdir='../results/odor_pulses'+dirextn)
            ## if the result file for these seeds & tweaks doesn't exist,
            ## then carry on to the next.
            if not os.path.exists(filename): continue
            worker = fit_sp.fit_plot_scaledpulses(filename,stimseed,_scaledWidth)
            fits_2mits,peaks_2mits = worker.fit_pulses(dirextn,True,False,False) # (dirextn,noshow,savefig,test)
            for mitrali in [0,1]:
                _,slopes,intercepts,r_values,_,_,se_slope,avg_frates = zip(*fits_2mits[mitrali])
                ## leave out the fit for reference with itself
                r_values = delete(r_values,fit_sp.ref_response_scalenum-1)
                R_all[ploti].extend(r_values)
                Rsqs_all[ploti].extend(r_values**2)
                nan_present = False
                for slope in slopes:
                    if isnan(slope): nan_present = True
                intercepts_bad = False
                #for intercept in intercepts:
                #    if abs(intercept)>2.0: intercepts_bad = True
                if not nan_present and not intercepts_bad:
                    slopes_all.append(slopes)
                    peaks_all.append(peaks_2mits[mitrali])
                    print filename
                else:
                    print 'Left out :',filename
                avg_frates_all.append(avg_frates)
        slopes_mean = append([0],mean(slopes_all,axis=0))
        slopes_std = append([0],std(slopes_all,axis=0))
        peaks_mean = append([0],mean(peaks_all,axis=0))
        peaks_std = append([0],std(peaks_all,axis=0))
        avg_frates_mean = mean(avg_frates_all,axis=0)
        print "The average firing rates for the different scaling is",avg_frates_mean
        print "Peaks mean is",peaks_mean
        
        ## plots
        axa.hist(R_all[ploti],10,range=(0,1.0),normed=True,histtype='step',\
            linewidth=linewidth,color='k',ls='solid')
        xmin,xmax,ymin,ymax = \
            beautify_plot(axa,x0min=True,y0min=True,xticksposn='bottom',yticksposn='left',yticks=[])
        axa.set_ylim(0,ymax) # later this is set to the max over all histograms maxy_hist
        maxy_hist = max(maxy_hist,ymax)
        ## x-axis label for full first row
        if inhi==2:
            axes_labels(axa,"correlation with 1x","",xpad=2)
            #axa.xaxis.set_label_coords(1.3,-0.25)
        ### inset plots within subplots
        ### passing transform=ax.transAxes to add_axes() doesn't work, hence jugglery from
        ### http://matplotlib.1069221.n5.nabble.com/Adding-custom-axes-within-a-subplot-td20316.html
        #Bbox = matplotlib.transforms.Bbox.from_bounds(1.6, 0.4, 0.6, 0.6) 
        #trans = ax.transAxes + fig1.transFigure.inverted() 
        #l, b, w, h = matplotlib.transforms.TransformedBbox(Bbox,trans).bounds
        #axinset = fig1.add_axes([l, b, w, h])
        xconcs = array(scaledList[1:])
        ref_scale = scaledList[fit_sp.ref_response_scalenum]
        print "Slopes are =",slopes_mean*ref_scale
        #axb.plot(range(6),range(6),color=(0,0,0.7,0.5),dashes=(2.0,1.0)) # linear reference
        axb.errorbar(x=append([0],xconcs),y=peaks_mean,yerr=peaks_std,\
            color='k',linewidth=linewidth,capsize=cap_size)
        #for slopes in slopes_all:
        #    axb.plot(xconcs,array(slopes)*ref_scale,linewidth=linewidth)
        _,_,_,ymax = beautify_plot(axb,x0min=True,y0min=True,\
            xticksposn='bottom',yticksposn='left',yticks=[])
        axb.set_xlim(0,5)
        axb.set_ylim(0,180)
        axb.set_xticks([0,1,5])
        axb.set_yticks([])
        ### Draw the twin y axis (turned off always by beautify_plot)
        #for loc, spine in axb.spines.items(): # items() returns [(key,value),...]
        #    spine.set_linewidth(axes_linewidth)
        #    if loc in ['right']:
        #        spine.set_color('k') # draw spine in black

        ## x-axis label for full second row
        if inhi==2:
            axes_labels(axb,"ORN scaling","",xpad=2)
            #axb.xaxis.set_label_coords(1.3,-0.25)
        if inhi==0: # laebl twin y axes of rightmost plot
            axb.set_yticks([0,180])
            axes_labels(axb,'','peak (Hz)',ypad=-6) # to set font size for twin y ticklabels

    for i,(_,_,_,_,_,_,(axa,axb)) in enumerate(inh_options):
        axa.set_ylim(0,maxy_hist)
        if i==0: # label y axis for left-most plots
            axa.set_yticks([0,maxy_hist])
            axes_labels(axa,'','density',ypad=2)
    fig1.tight_layout()
    fig_clip_off(fig1)
    fig1.subplots_adjust(top=0.95,left=0.1,bottom=0.15,right=0.98,wspace=0.4,hspace=0.4)
    #fig1.text(0.31,0.65,'density',fontsize=label_fontsize,\
    #        rotation='vertical', transform=fig1.transFigure)
    #fig1.text(0.6,0.025,'$R^2$',fontsize=label_fontsize,transform=fig1.transFigure)
    fig1.savefig('../figures/lin_contribs_scaledpulses.svg',dpi=fig1.dpi)
    fig1.savefig('../figures/lin_contribs_scaledpulses.png',dpi=fig1.dpi)

def plot_lin_contribs_paperfigure():
    """ plot linearity scores for all options
    """
    fig1 = figure(figsize=(columnwidth,linfig_height),dpi=300,facecolor='w')
    ax1a = plt.subplot2grid((2,5),(0,0)) # full
    ax1b = plt.subplot2grid((2,5),(1,0)) # full
    ax2a = plt.subplot2grid((2,5),(0,1)) # no lat
    ax2b = plt.subplot2grid((2,5),(1,1)) # no lat
    ax3a = plt.subplot2grid((2,5),(0,2)) # no self
    ax3b = plt.subplot2grid((2,5),(1,2)) # no self
    ax4a = plt.subplot2grid((2,5),(0,3)) # no inh
    ax4b = plt.subplot2grid((2,5),(1,3)) # no inh
    ax5a = plt.subplot2grid((2,5),(0,4)) # non-lin
    ax5b = plt.subplot2grid((2,5),(1,4))
    ## inh = (no_singles,no_joints,no_lat,no_PGs,varyRMP)
    inh_options = [
        ('',0,(False,False,False,False,False),'all',False,None,(ax1a,ax1b)), \
        ('',1,(False,False,True,False,False),'no lat-inh',False,None,(ax2a,ax2b)), \
        ('',2,(True,False,False,True,False),'no self-inh',False,None,(ax3a,ax3b)), \
        ('',3,(True,True,False,True,False),'no inh',False,None,(ax4a,ax4b)), \
        ('',4,(False,False,False,False,False),'non-lin ORNs',True,'P',(ax5a,ax5b)) ]
    Rsqs_all = [ [] for i in range(len(inh_options)) ]
    R_all = [ [] for i in range(len(inh_options)) ]
    maxy_R = 0.0
    for ploti,(dirextn,inhi,inh,inhstr,nl_orns,nl_type,(axa,axb)) in enumerate(inh_options):
        slopes_all = []
        peaks_all = []
        avg_frates_all = []
        R_thisinh = []
        for stimi,stimseed in enumerate(stim_seeds):
            filename, switch_strs \
                = get_filename(stimseed,stimseed,inh,0,nl_orns,nl_type,\
                    resultsdir='../results/odor_pulses'+dirextn)
            ## if the result file for these seeds & tweaks doesn't exist,
            ## then carry on to the next.
            if not os.path.exists(filename): continue
            worker = fit_sp.fit_plot_scaledpulses(filename,stimseed,_scaledWidth)
            fits_2mits,peaks_2mits = worker.fit_pulses(dirextn,True,False,False) # (dirextn,noshow,savefig,test)
            for mitrali in [0,1]:
                _,slopes,intercepts,r_values,_,_,se_slope,avg_frates = zip(*fits_2mits[mitrali])
                R_thisinh.append(r_values)
                ## leave out the fit for reference with itself
                r_values = delete(r_values,fit_sp.ref_response_scalenum-1)
                R_all[ploti].extend(r_values)
                Rsqs_all[ploti].extend(r_values**2)
                nan_present = False
                for slope in slopes:
                    if isnan(slope): nan_present = True
                intercepts_bad = False
                #for intercept in intercepts:
                #    if abs(intercept)>2.0: intercepts_bad = True
                if not nan_present and not intercepts_bad:
                    slopes_all.append(slopes)
                    peaks_all.append(peaks_2mits[mitrali])
                    print filename
                else:
                    print 'Left out :',filename
                avg_frates_all.append(avg_frates)
        slopes_mean = append([0],mean(slopes_all,axis=0))
        slopes_std = append([0],std(slopes_all,axis=0))
        peaks_mean = append([0],mean(peaks_all,axis=0))
        peaks_std = append([0],std(peaks_all,axis=0))
        avg_frates_mean = mean(avg_frates_all,axis=0)
        print "The average firing rates for the different scaling is",avg_frates_mean
        print "Peaks mean is",peaks_mean
        
        ## plots
        R_mean = mean(R_thisinh,axis=0)
        R_std = std(R_thisinh,axis=0)
        xconcs = scaledList[1:]
        axa.errorbar(xconcs,R_mean,R_std,linewidth=linewidth,color='b',\
            ls='solid',marker='o',ms=marker_size,capsize=cap_size)
        #axa.plot(xconcs,R_std,linewidth=linewidth,color='r',ls='solid',marker='x',ms=marker_size)
        xmin,xmax,ymin,ymax = \
            beautify_plot(axa,x0min=True,y0min=True,xticksposn='bottom',yticksposn='left',yticks=[])
        axa.set_ylim(0,ymax) # later this is set to the max over all histograms maxy_hist
        axa.set_xticks([0,1,2,5])
        maxy_R = max(maxy_R,ymax)
        ## x-axis label for full first row
        if inhi==2:
            axes_labels(axa,"","",xpad=2) # ORN scaling is on the next row x-axis
            #axa.xaxis.set_label_coords(1.3,-0.25)
        ### inset plots within subplots
        ### passing transform=ax.transAxes to add_axes() doesn't work, hence jugglery from
        ### http://matplotlib.1069221.n5.nabble.com/Adding-custom-axes-within-a-subplot-td20316.html
        #Bbox = matplotlib.transforms.Bbox.from_bounds(1.6, 0.4, 0.6, 0.6) 
        #trans = ax.transAxes + fig1.transFigure.inverted() 
        #l, b, w, h = matplotlib.transforms.TransformedBbox(Bbox,trans).bounds
        #axinset = fig1.add_axes([l, b, w, h])
        xconcs = array(scaledList[1:])
        ref_scale = scaledList[fit_sp.ref_response_scalenum]
        print "Slopes are =",slopes_mean*ref_scale
        #axb.plot(range(6),range(6),color=(0,0,0.7,0.5),dashes=(2.0,1.0)) # linear reference
        axb.errorbar(x=append([0],xconcs),y=peaks_mean,yerr=peaks_std,\
            color='b',linewidth=linewidth,capsize=cap_size)
        #for slopes in slopes_all:
        #    axb.plot(xconcs,array(slopes)*ref_scale,linewidth=linewidth)
        _,_,_,ymax = beautify_plot(axb,x0min=True,y0min=True,\
            xticksposn='bottom',yticksposn='left',yticks=[])
        axb.set_xlim(0,5)
        axb.set_ylim(0,180)
        axb.set_xticks([0,1,2,5])
        axb.set_yticks([])
        ### Draw the twin y axis (turned off always by beautify_plot)
        #for loc, spine in axb.spines.items(): # items() returns [(key,value),...]
        #    spine.set_linewidth(axes_linewidth)
        #    if loc in ['right']:
        #        spine.set_color('k') # draw spine in black

        ## x-axis label for full second row
        if inhi==2:
            axes_labels(axb,"conc (% SV)","",xpad=2)
            #axb.xaxis.set_label_coords(1.3,-0.25)
        if inhi==0: # laebl twin y axes of rightmost plot
            axb.set_yticks([0,180])
            axes_labels(axb,'','peak (Hz)',ypad=-6) # to set font size for twin y ticklabels

    for i,(_,_,_,_,_,_,(axa,axb)) in enumerate(inh_options):
        axa.set_ylim(0,1)
        if i==0: # label y axis for left-most plots
            axa.set_yticks([0,1])
            axes_labels(axa,'','corr to 1%',ypad=2)
    fig1.tight_layout()
    fig_clip_off(fig1)
    fig1.subplots_adjust(top=0.95,left=0.1,bottom=0.15,right=0.98,wspace=0.4,hspace=0.4)
    #fig1.text(0.31,0.65,'density',fontsize=label_fontsize,\
    #        rotation='vertical', transform=fig1.transFigure)
    #fig1.text(0.6,0.025,'$R^2$',fontsize=label_fontsize,transform=fig1.transFigure)
    fig1.savefig('../figures/lin_contribs_scaledpulses.svg',dpi=fig1.dpi)
    fig1.savefig('../figures/lin_contribs_scaledpulses.png',dpi=fig1.dpi)

def plot_responseshapes_vs_conc():
    fig = figure(figsize=(columnwidth/3.*2,linfig_height/2.),dpi=300,facecolor='w')
    ax1 = fig.add_subplot(1,2,1)
    ax2 = fig.add_subplot(1,2,2)
    xconcs = array(scaledList[1:])
    ## inh = (no_singles,no_joints,no_lat,no_PGs,varyRMP)
    inh_option = ('',0,(False,False,False,False,False),'all',False,None)
    tpeaks_all = []
    tpeaks_normed_all = []
    durations_normed_all = []
    risetimes_all = []
    durations_all = []
    deltapeaks_all = []
    latencies_all = []
    latencies_normed_all = []
    (dirextn,inhi,inh,inhstr,nl_orns,nl_type) = inh_option
    for stimi,stimseed in enumerate(stim_seeds):
        filename, switch_strs \
            = get_filename(stimseed,stimseed,inh,0,nl_orns,nl_type,\
                resultsdir='../results/odor_pulses'+dirextn)
        ## 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
        ######### load in the mitral pulse responses
        #### mitral_responses_list[avgnum][scalenum][mitralnum][spikenum]
        #### mitral_responses_binned_list[avgnum][scalenum][mitralnum][binnum]
        mitral_responses_list, mitral_responses_binned_list = read_pulsefile(filename)
        ##-------------------------- rebin the responses and pulses ------------------------------
        ## rebin sim responses to pulserebindt=50ms, then take mean
        numtrials,mitral_responses_mean,mitral_responses_std = \
                rebin_mean(mitral_responses_list,fit_sp.pulserebindt,SCALED_RUNTIME)
        for mitrali in [0,1]:
            ## take the odor responses of the mitral to be fitted
            firingbinsmeanList = mitral_responses_mean[:,mitrali]
            
            starti = int(PULSE_START/fit_sp.pulserebindt)
            endi = int((PULSE_START+scaledWidth+kernel_time)/fit_sp.pulserebindt)
            air_bgnd = firingbinsmeanList[0]
            air_bgnd_relevant = firingbinsmeanList[0][starti:endi]
            tpeaks = []
            risetimes = []
            durations = []
            deltapeaks = []
            latencies = []
            for scalenum in [1,2,3,4,5]: ## conc scaled pulses
                ## find the mean first spike latency after odor onset
                latency = 0
                for avgi in range(len(mitral_responses_list)):
                    for spiketime in mitral_responses_list[avgi][scalenum][mitrali]:
                        if spiketime > PULSE_START:
                            latency += spiketime - PULSE_START
                            break
                latency_mean = latency / float(len(mitral_responses_list))
                latencies.append(latency_mean)
                ## find peak times and durations of binned responses
                scale = scaledList[scalenum]
                response = firingbinsmeanList[scalenum][starti:endi]-air_bgnd_relevant
                tpeak = argmax(response)
                tpeaks.append(tpeak)
                ##--------- Calculate half-max duration of the response
                ## f is an interpolated 1D function that gives
                ## down-shifted response that is zero at halfmax for fsolve/bisect below
                f = interpolate.interp1d(range(len(response)),response-max(response)/2.,\
                        bounds_error=False,fill_value=0.0)
                ### fsolve to find half-max points on either side of tpeak
                ### but often gives negative values!! so ditch this
                #halfresp_low = optimize.fsolve(f,tpeak-1)[0]
                #halfresp_high = optimize.fsolve(f,tpeak+1)[0]
                ## use bisect to find halfmax pt within an interval
                ## but bisect need opposite signs at the ends of the interval
                ## else ValueError, so increase interval, till you get it
                tmin = tpeak-1
                while True:
                    try:
                        halfresp_low = optimize.bisect(f,tpeak,tmin-1)
                        break
                    except ValueError:
                        tmin -= 1
                tmax = tpeak+1
                while True:
                    try:
                        halfresp_high = optimize.bisect(f,tpeak,tmax)
                        break
                    except ValueError:
                        tmax += 1
                risetimes.append(halfresp_low)
                duration = (halfresp_high-halfresp_low)
                durations.append(duration)
                
            ## take only the reduction in the time to peak vs concentration
            tpeaks = array(tpeaks)*fit_sp.pulserebindt*1000
            risetimes = array(risetimes)*fit_sp.pulserebindt*1000
            durations = array(durations) *fit_sp.pulserebindt*1000
            tpeaks_all.append(tpeaks)
            tpeaks_normed_all.append(tpeaks/mean(tpeaks))
            risetimes_all.append(risetimes)
            durations_all.append(durations)
            durations_normed_all.append(durations/mean(durations))
            deltapeaks = [tpeaks[-1]-tpeaks[-2],tpeaks[-2]-tpeaks[-3]]
            deltapeaks_all.extend(deltapeaks)
            latencies_all.append(latencies)
            latencies_normed_all.append(latencies/mean(latencies))
            #ax1.plot(xconcs,tpeaks)
            #ax2.plot(xconcs,durations)
    
    tpeaks_mean = mean(tpeaks_all,axis=0)
    tpeaks_std = std(tpeaks_all,axis=0)
    tpeaks_normed_mean = mean(tpeaks_normed_all,axis=0)
    tpeaks_normed_se = std(tpeaks_normed_all,axis=0) / sqrt(len(tpeaks_normed_all))
    risetimes_mean = mean(risetimes_all,axis=0)
    risetimes_std = std(risetimes_all,axis=0)
    durations_mean = mean(durations_all,axis=0)
    durations_std = std(durations_all,axis=0)
    durations_normed_mean = mean(durations_normed_all,axis=0)
    durations_normed_se = std(durations_normed_all,axis=0) / sqrt(len(durations_normed_all))
    latencies_mean = mean(latencies_all,axis=0)
    latencies_std = std(latencies_all,axis=0)
    latencies_normed_mean = mean(latencies_normed_all,axis=0)
    latencies_normed_se = std(latencies_normed_all,axis=0) / sqrt(len(latencies_normed_all))
    ## plots
    ax1.errorbar(xconcs,latencies_normed_mean,latencies_normed_se,color='b',linewidth=linewidth)
    #ax1.errorbar(xconcs,risetimes_mean,risetimes_std,color='b',linewidth=linewidth)
    ax2.errorbar(xconcs,tpeaks_normed_mean,tpeaks_normed_se,color='b',linewidth=linewidth)
    #ax2.errorbar(xconcs,durations_normed_mean,durations_normed_se,color='b',linewidth=linewidth)
    for i,ax in enumerate([ax1,ax2]):
        _,_,ymin,ymax = beautify_plot(ax,x0min=True,y0min=False,\
            xticksposn='bottom',yticksposn='left')
        ax.set_xlim(0,5)
        ax.set_xticks([0,1,2,5])
        ax.set_yticks([ymin,1,ymax])
        axes_labels(ax,'conc (% SV)',['norm-t-spike','norm-t-peak'][i])
    ### Draw the twin y axis (turned off always by beautify_plot)
    #for loc, spine in axb.spines.items(): # items() returns [(key,value),...]
    #    spine.set_linewidth(axes_linewidth)
    #    if loc in ['right']:
    #        spine.set_color('k') # draw spine in black

    fig_clip_off(fig)
    fig.tight_layout()
    fig.savefig('../figures/scaledpulses_latency.svg',dpi=fig.dpi)
    fig.savefig('../figures/scaledpulses_latency.png',dpi=fig.dpi)
    
    fig = figure()
    ax = fig.add_subplot(111)
    ax.hist(deltapeaks_all)
    axes_labels(ax,'delta peak times','#')

if __name__ == "__main__":
    ## plot linearity with parts of the network removed / tweaked
    #plot_lin_contribs_oldpaperfigure()
    ## PAPER figure. Similar to above, but plotting mean and sd of corr with 1x
    ## versus concentration instead of a combined distrib.
    plot_lin_contribs_paperfigure()
    ## Compare pulse shapes versus concentration
    plot_responseshapes_vs_conc()
    
    show()