import logging
import sys
import os
from scipy.stats import norm
import numpy as np
from allensdk.internal.model.glif.glif_optimizer_neuron import GlifNeuronException
from allensdk.internal.model.glif.glif_optimizer_neuron import GlifBadInitializationException
from allensdk.model.glif.glif_neuron import GlifBadResetException
# TODO: clean up
# TODO: license
[docs]def MLIN_list_error(param_guess, experiment, input_data):
#TODO: binning is now done in preprocessor so perhaps should take it out of here.
voltage_variance = input_data['subthreshold_long_square_voltage_variance']
# voltage_distribution=norm(loc=0, scale=np.sqrt(voltage_variance)*10.)
sv=input_data['sv_for_expsymm'] #used in the expsymm function
# tau_4AC=experiment.neuron.R_input*experiment.neuron.C
tau_from_AC=input_data['tau_from_AC']
spike_length=int(experiment.neuron.spike_cut_length)
noSpike_bin_size_ind=int(tau_from_AC/experiment.neuron.dt)
spike_bin_size_time=.005 #TODO: PLAY WITH THIS VALUE. 1 MS, 2, 4, 8
spike_bin_size_ind=int(spike_bin_size_time/experiment.neuron.dt)
logging.info('running parameter guess: %s' % param_guess)
dataSpikeTimes=experiment.grid_spike_times
MLIN_list = []
try:
run_data = experiment.run(param_guess)
except GlifNeuronException as e:
out=e.data
raise e
modelSpikeISI=[e.data['interpolated_ISI']]
except GlifBadInitializationException as e:
logging.error('voltage STARTS above threshold: setting error to be large.Difference between thresh and voltage is: %f' % e.dv)
raise Exception()
except GlifBadResetException as e:
logging.error('THIS REALLY SHOULDNT HAPPEN WITH NEW INITIALIZATION EXCEPTION: voltage is above threshold at reset: setting error to be large. Difference between thresh and voltage is: %f' % e.dv)
raise Exception()
v_model_list=[]
th_model_list=[]
non_spike_bin_ind_edges_list=[]
spike_edges_list_list=[]
spike_bins_list=[]
noSpike_bins_list=[]
spike_prob_list=[]
noSpike_prob_list=[]
for stim_list_index in range(0,len(experiment.stim_list)):
#TODO: the following line is a hack to take care of the case when there are no spikes in a sweep
if len(experiment.spike_time_steps[stim_list_index])==0:
MLIN=0
raise Exception('there are no spikes in the sweep')
else:
bio_spike_ind=experiment.spike_time_steps[stim_list_index]
v_model=run_data['voltage'][stim_list_index]
th_model=run_data['threshold'][stim_list_index]
#------------------------------------------------------------------------------------
#---------------make all the bins----------------------------------------------------
#------------------------------------------------------------------------------------
#--for every spike define a region for making non spiking bins
between_spike_edges_list=[]
between_spike_edges_list.append([0, bio_spike_ind[0]-spike_bin_size_ind]) #edges to first spike
for ii in range(0, len(bio_spike_ind)-1):
between_spike_edges_list.append([bio_spike_ind[ii]+spike_length, bio_spike_ind[ii+1]-spike_bin_size_ind])
#--define no spike bin edges
non_spike_bin_ind_edges_in_ISI=[]
for btw_spike_edges in between_spike_edges_list:
temp=range(btw_spike_edges[0], btw_spike_edges[1], noSpike_bin_size_ind)
temp.append(btw_spike_edges[1])
non_spike_bin_ind_edges_in_ISI.append(temp)
#--define spike bin edges
spike_edges_list=[]
for spike in bio_spike_ind:
spike_edges_list.append([spike-spike_bin_size_ind, spike])
#--nonspike bin edges need to be arranged correctly because they are just all the edges for a given ISI
non_spike_bin_ind_edges=[]
for edges_in_1_spike in non_spike_bin_ind_edges_in_ISI:
for ii in range(0,len(edges_in_1_spike)-1):
non_spike_bin_ind_edges.append([edges_in_1_spike[ii], edges_in_1_spike[ii+1]])
#-----------------------------------------------------------------------------
#-------------finding values in bins------------------------------------------
#-----------------------------------------------------------------------------
spike_bins={}
# spike_bins['vmax']=[]
spike_bins['th']=[]
spike_bins['v']=[]
spike_bins['v_th_diff']=[]
# spike_bins['ind_of_max_v']=[]
spike_bins['ind_of_max_diff_btw_v_th']=[]
for bin_edges in spike_edges_list:
bin_ind=range(bin_edges[0], bin_edges[1])
#vmax_in_bin=max(v_model[bin_ind])
diff_vector=th_model[bin_ind]-v_model[bin_ind]
#the_ind=bin_ind[np.where(v_model[bin_ind]==vmax_in_bin)[0]] this is used to use in where v is at a max (as opposed to the difference between v and th)
# print("***********************************************************")
# print('th_model[bin_ind]', th_model[bin_ind])
# print('v_model[bin_ind]',v_model[bin_ind])
# print('diff_vector', diff_vector)
# print('np.where(diff_vector==min(diff_vector))[0]', np.where(diff_vector==min(diff_vector))[0])
the_ind=bin_ind[np.where(diff_vector==min(diff_vector[~np.isnan(diff_vector)]))[0][0]]
th_in_bin=th_model[the_ind]
v_in_bin=v_model[the_ind]
# diffV=th_in_bin-vmax_in_bin
diffV=th_model[the_ind]-v_model[the_ind]
spike_bins['v'].append(v_in_bin)
spike_bins['th'].append(th_in_bin)
spike_bins['v_th_diff'].append(diffV)
spike_bins['ind_of_max_diff_btw_v_th'].append(the_ind)
noSpike_bins={}
# noSpike_bins['vmax']=[]
noSpike_bins['th']=[]
noSpike_bins['v']=[]
noSpike_bins['v_th_diff']=[]
noSpike_bins['ind_of_max_diff_btw_v_th']=[]
for bin_edges in non_spike_bin_ind_edges:
bin_ind=range(bin_edges[0], bin_edges[1])
# vmax_in_bin=max(v_model[bin_ind])
diff_vector=th_model[bin_ind]-v_model[bin_ind]
# max_indicies=np.where(v_model[bin_ind]==vmax_in_bin)[0]
min_indicies=np.where(diff_vector==min(diff_vector))[0]
# if len(max_indicies)>1:
# print('there is more than one maximum indicie in a bin at', max_indicies, 'choosing last value for computation')
# print('all voltages in the bin are', v_model[bin_ind])
the_ind=bin_ind[min_indicies[-1]] #this is here just incase there is more than one value at max voltage in a bin
th_in_bin=th_model[the_ind]
v_in_bin=v_model[the_ind]
diffV=th_model[the_ind]-v_model[the_ind]
noSpike_bins['v'].append(v_in_bin)
noSpike_bins['th'].append(th_in_bin)
noSpike_bins['v_th_diff'].append(diffV)
noSpike_bins['ind_of_max_diff_btw_v_th'].append(the_ind)
#-----------------------------------------------------------------------------
#-------------calculate MLIN--------------------------------------------------
#-----------------------------------------------------------------------------
# this was the version with the normal distribution
# noSpike_prob=np.log(np.spacing(1)+voltage_distribution.cdf(noSpike_bins['v_th_diff']))
# spike_prob=np.log(1+np.spacing(1)-voltage_distribution.cdf(spike_bins['v_th_diff']))
#version with the expsymm function
#---OPTION ONE--------
# noSpike_prob=np.log(np.spacing(1)+expsymm_cdf(noSpike_bins['v_th_diff'], sv))
# spike_prob=np.log(1+np.spacing(1)-expsymm_cdf(spike_bins['v_th_diff'], sv))
# #---OPTION TWO--------
# N_spike=np.float(len(spike_bins['v_th_diff']))
# N_noSpike=np.float(len(noSpike_bins['v_th_diff']))
# noSpike_prob=(N_spike/(N_noSpike+N_spike))*np.log(np.spacing(1)+expsymm_cdf(noSpike_bins['v_th_diff'], sv))
# spike_prob=(N_noSpike/(N_noSpike+N_spike))*np.log(1+np.spacing(1)-expsymm_cdf(spike_bins['v_th_diff'], sv))
#---OPTION THREE
noSpike_negDiff=(-np.log(2.)+np.array(noSpike_bins['v_th_diff'])/sv)[np.array(noSpike_bins['v_th_diff'])<=0.0]
noSpike_posDiff=(np.log(1.-0.5*np.exp(-np.array(noSpike_bins['v_th_diff'])/sv)))[np.array(noSpike_bins['v_th_diff'])>0.0]
noSpike_prob=np.append(noSpike_negDiff, noSpike_posDiff) #!!NOTE: this may not line up correctly in outputs of MLIN HACK
spike_negDiff=(np.log(1.-0.5*np.exp(np.array(spike_bins['v_th_diff'])/sv)))[np.array(spike_bins['v_th_diff'])<=0.0]
spike_posDiff=(-np.log(2.)-np.array(spike_bins['v_th_diff'])/sv)[np.array(spike_bins['v_th_diff'])>0.0]
spike_prob=np.append(spike_negDiff, spike_posDiff)
MLIN=-(sum(noSpike_prob)+sum(spike_prob))
logging.info('MLIN: %f', MLIN)
MLIN_list.append([MLIN])
v_model_list.append(v_model)
th_model_list.append(th_model)
non_spike_bin_ind_edges_list.append(non_spike_bin_ind_edges)
spike_edges_list_list.append(spike_edges_list)
spike_bins_list.append(spike_bins)
noSpike_bins_list.append(noSpike_bins)
spike_prob_list.append(spike_prob)
noSpike_prob_list.append(noSpike_prob)
concatenateMLINList=np.concatenate(MLIN_list)
experiment.spike_errors.append(concatenateMLINList)
# print('param Guess', param_guess, 'TRD', np.mean(concatenateTRDList))
out =np.mean(concatenateMLINList)
logging.info('MLIN: %f', np.mean(concatenateMLINList))
#-------------------------------------------------------------------
#--------------------------plotting------------------------------------
#---------------------------------------------------------------------
time=np.arange(0, len(v_model))*experiment.neuron.dt
# plt.subplot(2,1,1)
# plt.title('Model', fontsize=16)
# plt.plot(time, v_model, 'b-', label='voltage')
# plt.plot(time, th_model, 'b--', label='threshold')
# plt.plot(np.concatenate(non_spike_bin_ind_edges)*experiment.neuron.dt, v_model[np.concatenate(non_spike_bin_ind_edges)], 'k|', ms=16)
# plt.plot(np.concatenate(spike_edges_list)*experiment.neuron.dt, v_model[np.concatenate(spike_edges_list)], 'r|', ms=16)
# plt.plot(np.array(spike_bins['ind_of_max_diff_btw_v_th'])*experiment.neuron.dt, v_model[spike_bins['ind_of_max_diff_btw_v_th']], 'r.', ms=6)
# plt.plot(np.array(spike_bins['ind_of_max_diff_btw_v_th'])*experiment.neuron.dt, th_model[spike_bins['ind_of_max_diff_btw_v_th']], 'r.', ms=6)
# plt.plot(np.array(noSpike_bins['ind_of_max_diff_btw_v_th'])*experiment.neuron.dt, v_model[noSpike_bins['ind_of_max_diff_btw_v_th']], 'k.', ms=6)
# plt.plot(np.array(noSpike_bins['ind_of_max_diff_btw_v_th'])*experiment.neuron.dt, th_model[noSpike_bins['ind_of_max_diff_btw_v_th']], '.', ms=6)
# plt.title(' MLIN='+ str(MLIN)+' : sv='+str(sv)+' : ac_tau='+str(tau_from_AC)+' : spike bin size='+str(spike_bin_size_time)+'!!!!!! !!!!!', fontsize=20)
# plt.xlim([0, time[-1]])
#
# plt.subplot(2,1,2)
# plt.plot(np.array(spike_bins['ind_of_max_diff_btw_v_th'])*experiment.neuron.dt, spike_prob, 'r.', ms=16, label='spike probablility')
# plt.plot(np.array(noSpike_bins['ind_of_max_diff_btw_v_th'])*experiment.neuron.dt, noSpike_prob, 'b.', ms=16, label='no spike probablility')
# plt.legend()
#
# print("coming out of function", out)
#
# plt.show()
#converted to list 4_11_15
experiment.MLIN_HACK={
'v_model': v_model_list,
'th_model': th_model_list,
'non_spike_bin_ind_edges': non_spike_bin_ind_edges_list,
'spike_edges_list': spike_edges_list_list, #TODO Corinne: why was this originally a list but nothing else a list
'spike_bins': spike_bins_list,
'noSpike_bins': noSpike_bins_list,
'tau_from_AC': tau_from_AC,
'spike_prob': spike_prob_list,
'noSpike_prob': noSpike_prob_list,
'spike_bin_size_time' : spike_bin_size_time,
'sv': sv
}
#
return out