# -*- 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_scaledpulses.py <../results/odor_pulses/scaledpulses_ ... .pickle> <stimseed> [SAVEFIG]
from scipy import optimize
from scipy import special
import scipy.interpolate
from scipy import signal
import scipy.io
from scipy.stats import linregress
from pylab import *
import pickle
import sys
import math
import copy as cp
sys.path.extend(["..","../networks","../generators","../simulations"])
from stimuliConstants import * # has SETTLETIME, SCALED_RUNTIME
from networkConstants import * # has central_glom
from data_utils import * # has axes_labels()
from analysis_utils import * # has load_morph...(), predict_respresp() and set_min_err()
## index of the reference response
ref_response_scalenum = 3
## if peak detect is True, peak scaling is plotted,
## rather than the scaling wrt reference response above
peak_detect = False#True
## rebin the spikes as below, irrespective of previous binning
pulserebindt = 50e-3#fitting_dt # 50ms as Priyanka uses
pulserebins = int(SCALED_RUNTIME/pulserebindt)
#bin_width_time = 2*pulserebindt ## overlapping bins with this bin width
pulsetlist = arange(pulserebindt/2.0,SCALED_RUNTIME,pulserebindt)
NOISE_ANALYSIS = False
## number of the mitral to be fitted.
fitted_mitrals = [2*central_glom,2*central_glom+1]
class fit_plot_scaledpulses():
def __init__(self,filename,stimseed,scaledWidth=scaledWidth):
self.filename = filename
self.stimseed = stimseed
self.scaledWidth = scaledWidth
def fit_pulses(self,dirextn='',_noshow=True,_savefig=False, _test=False):
if _test:
########## load in the stimuli
## scaledPulseList[glomnum][scalenum][odornum][binnum]
scaledPulsesList,genORNkernels \
= read_scaledpulses_stimuli(self.stimseed,self.scaledWidth)
## decimate pulseList by taking every convolutiondt/FIRINGFILLDT bins
## and decimate ORNkernels by taking every kerneldt/FIRINGFILLDT=50/1=50 bins
## pulserebinList[scalenum,binnum] numpy array
pulserebinList,linORNkernelA,linORNkernelB = \
decimate(scaledPulsesList[0,:,0],pulserebindt,genORNkernels,kerneldt,kernel_size)
numtrials = 1
else:
######### 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(self.filename)
if NOISE_ANALYSIS:
##---------- Print/plot for best bin size, but this bin size is not used presently
## From Neural Computation 19, 1503–1527 (2007) Shimazaki & Shinomoto, pg 6
DeltaList = arange(25.0e-3,500.0e-3,25.0e-3)
SSvals = []
SSvals_extended = []
for Delta in DeltaList:
numtrials,mitral_responses_mean1,mitral_responses_std1 = \
rebin_mean(mitral_responses_list,Delta,SCALED_RUNTIME)
fitted_mitral = 0
firingbinsmeanList1 = mitral_responses_mean1[:,fitted_mitral]
firingbinserrList1 = mitral_responses_std1[:,fitted_mitral]
## each val in firingbinsmeanList[scalenum][binnum] is k_i/(n*Delta)
kbar = mean(firingbinsmeanList1)*numtrials*Delta
kvariance = (std(firingbinsmeanList1)*numtrials*Delta)**2
shimazaki_shinomoto_val = (2*kbar-kvariance)/((numtrials*Delta)**2)
SSvals.append(shimazaki_shinomoto_val)
print 'Delta t =',Delta,"shimazaki_shinomoto_val =",shimazaki_shinomoto_val
## Cost function i.e. shimazaki_shinomoto_val for different # of trials
Priyanka_numtrials = 12
ssval_extended = (1.0/Priyanka_numtrials - 1.0/numtrials)*kbar/numtrials/Delta**2 + \
shimazaki_shinomoto_val
SSvals_extended.append(ssval_extended)
if not _noshow:
fig = figure()
ax = fig.add_subplot(111)
ax.plot(DeltaList,SSvals,label='#='+str(numtrials))
ax.plot(DeltaList,SSvals_extended,label='#='+str(Priyanka_numtrials))
ax.legend()
axes_labels(ax,'Delta t (s)','SS Cost function (Hz^2)',fontsize=20)
##-------------------------- 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,pulserebindt,SCALED_RUNTIME)
## full fitting data for both mitrals
fits_2mits = []
peak_scales_2mits = []
for mit_i,fitted_mitral in enumerate(fitted_mitrals):
if _test:
firingbinsmeanList = pulserebinList
firingbinsmeanList += uniform(-0,0,shape(firingbinsmeanList))
firingbinserrList = zeros(shape(firingbinsmeanList))
else:
## take the odor responses of the mitral to be fitted
firingbinsmeanList = mitral_responses_mean[:,fitted_mitral]
## The model predicts the individual response not the mean.
## Hence below fitting uses standard deviation, not standard error of the mean.
firingbinserrList = mitral_responses_std[:,fitted_mitral]
starti = int(PULSE_START/pulserebindt)
endi = int((PULSE_START+scaledWidth+kernel_time)/pulserebindt)
air_bgnd = firingbinsmeanList[0]
air_bgnd_relevant = firingbinsmeanList[0][starti:endi]
##---------------------------- fit scaled pulse responses ------------------------------------------
## define the reference response / scaling
ref_scale = scaledList[ref_response_scalenum]
ref_response = firingbinsmeanList[ref_response_scalenum][starti:endi]-air_bgnd_relevant
fits = []
peak_scales = []
for scalenum in [1,2,3,4,5]: ## conc scaled pulses
scale = scaledList[scalenum]
response = firingbinsmeanList[scalenum][starti:endi]-air_bgnd_relevant
## http://www.jerrydallal.com/LHSP/slrout.htm for defn of std error of the estimate
## SEE is std error of data about the regression line
slope, intercept, r_value, p_value, see = linregress(ref_response,response)
## SEE is called \hat{\sigma}_\epsilon i.e. sqrt(MSE) here:
## http://en.wikipedia.org/wiki/Regression_analysis
## formula for std error of slope is also from above Wikipedia article
se_slope = see/std(ref_response)
avgfrate = sum(firingbinsmeanList[scalenum][starti:endi])/float(endi-starti)
fits.append((scale,slope,intercept,r_value,p_value,see,se_slope,avgfrate))
## peak scaling
peak_scales.append(max(response))
peak_scales = array(peak_scales)#/peak_scales[ref_response_scalenum-1]*scaledList[ref_response_scalenum]
peak_scales_2mits.append(peak_scales)
fits_2mits.append(fits)
##---------------------------- plot scaled pulse responses -----------------------------------------
if not _noshow:
if peak_detect: print "BEWARE. Using peak scaling"
else: print "BEWARE. Using fitted scaling"
############################### plot the responses and the fits
fig = figure(figsize=(columnwidth,linfig_height/2),dpi=300,facecolor='w') # 'none' is transparent
## conc scaled pulses, leave the 0th pulse which is air_bgnd
ax = plt.subplot2grid((1,3),(0,0),rowspan=1)
for scaleiter,scale in enumerate(scaledList[1:]): ## conc scaled pulses
sister_ratio = (fitted_mitral%MIT_SISTERS)/float(MIT_SISTERS)
scaledpulsetime = array(pulsetlist[starti:endi]) - pulsetlist[starti] # start from t=0
################### 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[scaleiter+1]-air_bgnd,mode='same')
## ditch the smoothing above for scaled pulses
simresponse = firingbinsmeanList[scaleiter+1][starti:endi]-air_bgnd_relevant
## numpy array, hence adds element by element
scale_color = ['r','b','g','m','c'][scaleiter]
scale_label = ['1/3x','2/3x','1x','2x','5x'][scaleiter]+\
' %1.2f'%(fits_2mits[-1][scaleiter][3]) # R-value = corr
print scale_label
fill_between(scaledpulsetime,
simresponse+firingbinserrList[scaleiter+1][starti:endi]/sqrt(numtrials),
simresponse-firingbinserrList[scaleiter+1][starti:endi]/sqrt(numtrials),
color=scale_color,alpha=0.3)
plot(scaledpulsetime,simresponse,linewidth=plot_linewidth,color=scale_color,label=scale_label)
xmin,xmax,ymin,ymax = beautify_plot(ax,x0min=True,y0min=False,xticksposn='bottom',yticksposn='left')
#biglegend(fontsize=label_fontsize-2, labelspacing=0., handletextpad=0.)
ax.set_xticks([0,xmax])
if ymin<10: ax.set_yticks([ymin,0,ymax])
else: ax.set_yticks([ymin,ymax])
plot([0.,scaledWidth],[ymin+2,ymin+2],linewidth=plot_linewidth*3,color='r')
axes_labels(ax,'time (s)','rate (Hz)',adjustpos=False,xpad=0,ypad=-4)
#ax.yaxis.set_label_coords(-0.4,1.2)
## plot corr vs conc
ax = plt.subplot2grid((1,3),(0,1))
concratios,slopevsconc,_,_,_,_,se_slope,avg_frates = zip(*fits)
ax.plot(scaledList[1:],array(fits_2mits[-1])[:,3],color='b',linewidth=linewidth,\
marker='o',ms=marker_size)
beautify_plot(ax,x0min=False,y0min=False,xticksposn='bottom',yticksposn='left')
ax.set_xlim(0,5)
ax.set_xticks([0,1,2,5])
ax.set_ylim(0,1)
ax.set_yticks([0,1])
axes_labels(ax,'conc (% SV)','corr to 1%',adjustpos=False,xpad=1,ypad=2)
## plot response scaling vs conc scaling
ax = plt.subplot2grid((1,3),(0,2))
concratios,slopevsconc,_,_,_,_,se_slope,avg_frates = zip(*fits)
print "Average firing rates for mitral",fitted_mitral,\
"for different scales is",avg_frates
print "Peak firing rates for mitral",fitted_mitral,\
"for different scales is",peak_scales
ax.plot(scaledList,append([0],peak_scales),color='b',linewidth=linewidth,\
marker='o',ms=marker_size)
beautify_plot(ax,x0min=False,y0min=True,xticksposn='bottom',yticksposn='left')
ax.set_xlim(0,5)
ax.set_xticks([0,1,2,5])
axes_labels(ax,'conc (% SV)','peak (Hz)',adjustpos=False,xpad=1,ypad=-3)
#ax.yaxis.set_label_coords(-0.3,1.2)
fig.tight_layout()
fig_clip_off(fig)
#fig.subplots_adjust(top=0.94,left=0.1,right=0.99,hspace=0.4,wspace=0.5)
if _savefig:
fig.savefig('../figures/scalelinearity_example_'+str(self.stimseed)+\
'_mit'+str(mit_i)+'.svg',dpi=fig.dpi)
fig.savefig('../figures/scalelinearity_example_'+str(self.stimseed)+\
'_mit'+str(mit_i)+'.png',dpi=fig.dpi)
if NOISE_ANALYSIS:
## plot the variance vs firing rate mean for each mitral
## variance = mean/bintime of firng rate for Poisson process
fig2 = figure()
ax2 = fig2.add_subplot(111)
for scaleiter,scale in enumerate(scaledList): ## conc scaled pulses including air
ax2.scatter( firingbinsmeanList[scaleiter], \
firingbinserrList[scaleiter]**2, color='r' )
beautify_plot(ax2)
axes_labels(ax2,'mean rate (Hz)','variance (Hz^2)',fontsize=14)
## plot individual trials for a given response
fig3 = figure()
ax3 = fig3.add_subplot(111)
for trialspikelist in mitral_responses_list:
plot( plotBins( trialspikelist[0][fitted_mitral],\
pulserebins, SCALED_RUNTIME, 0.0) )
return fits_2mits, peak_scales_2mits
if __name__ == "__main__":
NOSHOW = False
if 'SAVEFIG' in sys.argv: SAVEFIG = True
else: SAVEFIG = False
if len(sys.argv) > 2:
filename = sys.argv[1]
stimseed = sys.argv[2]
worker = fit_plot_scaledpulses(filename,stimseed)
post_pulses = filename.split('odor_pulses')[1]
dirextn = post_pulses.split('/')[0]
print 'directory extension is',dirextn
if 'TEST' in sys.argv: TEST=True
else: TEST=False
worker.fit_pulses(dirextn,NOSHOW,SAVEFIG,TEST)
show()
else:
print "At least specify data file containing pickled mitral responses, and ORN frate seed."
sys.exit(1)
