Source code for allensdk.internal.model.glif.error_functions

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