run_networked_neuron

RUN_NETWORKED_NEURON calculates the firing behavior of a single theta neuron

Description: 
Calculates the firing behavior of a single theta neuron with n input 
synapses and q input spikes. q=n for SISO in which there is exactly one 
spike per synapse. Likewise, there are r output spikes. There will be at
most one output spike after the last input spike, depending on the neuron 
phase and the sign of the baseline current, Io.

Syntax:
RUN_NETWORKED_NEURON(ThNN, ti, NeuronIndex);
RUN_NETWORKED_NEURON(..., [GradFlag], [Verbose]);
[ts, Theta, Y, Z, thm, thp] = RUN_NETWORKED_NEURON(...);

Input Parameters:
o ThNN: An object of the theta neuron network class
o ti: A qx2 array that contains the neuron indices (globally indexed to
    the network) for each spike time (column 1) and the input spike times 
    (column 2)
o NeuronIndex: The index of the current neuron within the network ThNN
o GradFlag: Optional boolean flag to indicate whether the gradient matrices Y and Z 
    should be calculated.  If the flag is 0, then -1 is returned for Y and Z.
    The default value is 1.
o Verbose: Optional flag to indicate if extra information will be
    displayed on the screen. A value of 0 displays no additional information (this is 
    the default value), while a value of 1 displays all information.  Values greater 
    than 1 display partial information. See Verbose for more details.

Output Parameters:
o ts: An rx1 array of the output firing times, returns an empty array if no
    output spikes occur
o Theta: An array of size 2+ts(end)/TimeStep. If numerical integration is used, the 
    complete trajectory of the phase is returned.  For event-driven simulation, -1 
    is returned (since the trajectory is not calculated).
o Y: An nxr matrix of the effect of each synapse n on output spike r, used
     in network training
o Z: A qxr matrix of the effect of each input spike q on output spike r,
     used in network training
o thm: A qx1 array of phase values directly before each input spike
o thp: A qx1 array of phase values directly after each input spike

Examples:
>> ThNN=theta_neuron_network;
>> ts = run_networked_neuron(ThNN, [1 3; 2 4; 1 2], 3);

>> ThNN=theta_neuron_network('Io',0.05,'SimulationMethod','Numerical');
>> [ts, Theta] = run_networked_neuron(ThNN, [1 3; 2 50], 3);
>> figure;
>> plot(ThNN.TimeStep*(1:length(Theta)),Theta,ts,zeros(1,length(ts)),'o');
>> xlabel('Time (ms)'); ylabel('Phase'); grid on;

See also theta_neuron_network, verbose

0001 function [ts, Theta, Y, Z, thm, thp] = run_networked_neuron(ThNN, ti, NeuronIndex, GradFlag, Verbose)
0002 %RUN_NETWORKED_NEURON calculates the firing behavior of a single theta neuron
0003 %
0004 %Description:
0005 %Calculates the firing behavior of a single theta neuron with n input
0006 %synapses and q input spikes. q=n for SISO in which there is exactly one
0007 %spike per synapse. Likewise, there are r output spikes. There will be at
0008 %most one output spike after the last input spike, depending on the neuron
0009 %phase and the sign of the baseline current, Io.
0010 %
0011 %Syntax:
0012 %RUN_NETWORKED_NEURON(ThNN, ti, NeuronIndex);
0013 %RUN_NETWORKED_NEURON(..., [GradFlag], [Verbose]);
0014 %[ts, Theta, Y, Z, thm, thp] = RUN_NETWORKED_NEURON(...);
0015 %
0016 %Input Parameters:
0017 %o ThNN: An object of the theta neuron network class
0018 %o ti: A qx2 array that contains the neuron indices (globally indexed to
0019 %    the network) for each spike time (column 1) and the input spike times
0020 %    (column 2)
0021 %o NeuronIndex: The index of the current neuron within the network ThNN
0022 %o GradFlag: Optional boolean flag to indicate whether the gradient matrices Y and Z
0023 %    should be calculated.  If the flag is 0, then -1 is returned for Y and Z.
0024 %    The default value is 1.
0025 %o Verbose: Optional flag to indicate if extra information will be
0026 %    displayed on the screen. A value of 0 displays no additional information (this is
0027 %    the default value), while a value of 1 displays all information.  Values greater
0028 %    than 1 display partial information. See Verbose for more details.
0029 %
0030 %Output Parameters:
0031 %o ts: An rx1 array of the output firing times, returns an empty array if no
0032 %    output spikes occur
0033 %o Theta: An array of size 2+ts(end)/TimeStep. If numerical integration is used, the
0034 %    complete trajectory of the phase is returned.  For event-driven simulation, -1
0035 %    is returned (since the trajectory is not calculated).
0036 %o Y: An nxr matrix of the effect of each synapse n on output spike r, used
0037 %     in network training
0038 %o Z: A qxr matrix of the effect of each input spike q on output spike r,
0039 %     used in network training
0040 %o thm: A qx1 array of phase values directly before each input spike
0041 %o thp: A qx1 array of phase values directly after each input spike
0042 %
0043 %Examples:
0044 %>> ThNN=theta_neuron_network;
0045 %>> ts = run_networked_neuron(ThNN, [1 3; 2 4; 1 2], 3);
0046 %
0047 %>> ThNN=theta_neuron_network('Io',0.05,'SimulationMethod','Numerical');
0048 %>> [ts, Theta] = run_networked_neuron(ThNN, [1 3; 2 50], 3);
0049 %>> figure;
0050 %>> plot(ThNN.TimeStep*(1:length(Theta)),Theta,ts,zeros(1,length(ts)),'o');
0051 %>> xlabel('Time (ms)'); ylabel('Phase'); grid on;
0052 %
0053 %See also theta_neuron_network, verbose
0054 
0055 %Copyright (C) 2008 Sam McKennoch <Samuel.McKennoch@loria.fr>
0056 
0057 
0058 %Initialize return parameters to error return values
0059 Theta=-1; ts=[]; Y=-1; Z=-1; thm=-1; thp=-1;
0060 
0061 %Error handling for input arguments
0062 if nargin<3
0063     disp('Error in run_networked_neuron: Not enough input arguements');
0064     disp(['Needed at least 3 inputs but only got ' num2str(nargin)]);
0065     return;
0066 end
0067 if (nargin < 4)
0068     GradFlag=1;
0069 end
0070 if (nargin < 5)
0071     Verbose=0;
0072 end
0073 
0074 
0075 %Extract Parameters from Network and Neuron
0076 Delays=ThNN.Delays;
0077 Weights=ThNN.Weights;
0078 SimulationMethod=ThNN.SimulationMethod;
0079 TimeStep=ThNN.TimeStep;
0080 NIPS2007Gradient=ThNN.NIPS2007Gradient;
0081 try
0082     CurrentNeuron=struct(ThNN.Neurons(NeuronIndex));
0083 catch
0084     disp('Error in run_networked_neuron: NeuronIndex is not valid');
0085     return;
0086 end
0087 Io=CurrentNeuron.Io;
0088 Alpha=CurrentNeuron.Alpha;
0089 Beta=CurrentNeuron.Beta;
0090 ThFPP=CurrentNeuron.PositiveFixedPoint;
0091 
0092 %Initialize pre and post spike phase matrices
0093 thm=zeros(1,size(ti,1));
0094 thp=thm;
0095 
0096 %Order the input spike times and thus the Weights and Delays as well
0097 for j=1:size(ti,1)
0098     ti(j,2)=max(ti(j,2)+Delays(ti(j,1),NeuronIndex),0);
0099 end
0100 if ~isempty(ti)
0101     %Inline sortrows
0102     %ti=sortrows(ti,2);
0103     [Temp1,Temp2]=sort(ti(:,2));
0104     ti=[ti(Temp2,1) Temp1];    
0105 end
0106 
0107 %Verbose Display Header
0108 if (Verbose==1)
0109     disp(' ');
0110     disp('********** Function "run" for single neurons in a theta neuron network **********');
0111     disp(['---> Using Simulation Method: ', num2str(SimulationMethod)]);
0112 end
0113 
0114 
0115 
0116 % Analytical Integration %%
0117 if strcmp(SimulationMethod,'Event-Driven')
0118     prev_time=0;
0119     prev_thp=CurrentNeuron.Phase;
0120     for j=1:size(ti,1)
0121         %Without simplification for analytical (this doesn't remove the
0122         %imaginary numbers)
0123         %Fthp=(-1/Beta)*atanh(tan(prev_thp/2)/Beta)
0124         %thm(j)=2*atan(Beta*tanh(-Beta*(ti(j)-prev_time+Fthp)))
0125 
0126         %Calculate the next Theta_minus based on next spike time
0127         c1=tan(prev_thp/2)/Beta;
0128         if Io<0
0129             thm(j)=2*atan( Beta*(exp(-2*Beta*(ti(j,2)-prev_time))*(1+c1)-(1-c1))...
0130                 /(exp(-2*Beta*(ti(j,2)-prev_time))*(1+c1)+1-c1) );
0131         else
0132             thm(j)=2*atan(Beta*((tan(Beta*(ti(j,2)-prev_time))+c1)/(1-c1*tan(Beta*(ti(j,2)-prev_time)))));
0133         end
0134         thp(j)=thm(j)+Alpha*Weights(ti(j,1),NeuronIndex)*(1+cos(thm(j)));
0135 
0136         %Check for firing on spike
0137         if (thp(j)>pi && thm(j)<=pi)
0138             ts(end+1)=ti(j,2);
0139             thp(j)=thp(j)-2*pi;
0140             if j<size(ti,1) && Verbose
0141                 disp(['---> Neuron ', num2str(NeuronIndex),' Output Spike Time ',num2str(length(ts)),': ' num2str(ts(end)), '(Fired on Non-Final Spike)']);
0142             elseif (Verbose==1)
0143                 disp(['---> Neuron ', num2str(NeuronIndex),' Output Spike Time ',num2str(length(ts)),': ' num2str(ts(end)), '(Fired on Final Input Spike)']);
0144             end
0145             %Check for firing on trajectory; If Io<0 the previous thp must be
0146             %between ThFPP and pi and thm must be between -pi and -ThFPP in order
0147             %for the neuron to have fired on trajectory, regardless of the weight.
0148         elseif (Io<0 && prev_thp<pi && prev_thp>ThFPP && thm(j)>-pi && thm(j)<=(-ThFPP+10*eps))
0149             c1=tan(prev_thp/2)/Beta;
0150             ts(end+1)=prev_time+(1/Beta)*atanh(1/c1);
0151             if (Verbose==1)
0152                 disp(['---> Neuron ', num2str(NeuronIndex),' Output Spike Time ',num2str(length(ts)),': ' num2str(ts(end)),' (Fired on Trajectory)']);
0153             end
0154             %If Io>=0, we have to check the elapsed time between input spikes
0155             %and compare that to the unforced output spike frequency.
0156         elseif (Io>0 && (thm(j)<prev_thp || (ti(j,2)-prev_time)>(pi/Beta) ))
0157             ts(end+1)=prev_time+(1/Beta)*( (pi/2) - atan(c1) );
0158             if (Verbose==1)
0159                 disp(['---> Neuron ', num2str(NeuronIndex),' Output Spike Time ',num2str(length(ts)),': ' num2str(ts(end)), ' (Fired on Trajectory)']);
0160             end
0161             %Check for case where phase completes multiple
0162             %revolutions around phase circle between input spikes
0163             if (ti(j,2)-prev_time)>(pi/Beta)
0164                 num_revs=floor((ti(j,2)-prev_time)/(pi/Beta));
0165                 for k=2:num_revs
0166                     ts(end+1)=ts(end)+pi/Beta;
0167                     if (Verbose==1)
0168                         disp(['---> Neuron ', num2str(NeuronIndex),' Output Spike Time ',num2str(length(ts)),': ' num2str(ts(end)), ' (Multiple Revolutions)']);
0169                     end
0170                 end
0171                 if (ts(end)+pi/Beta)<ti(j,2)
0172                     ts(end+1)=ts(end)+pi/Beta;
0173                     if (Verbose==1)
0174                         disp(['---> Neuron ', num2str(NeuronIndex),' Output Spike Time ',num2str(length(ts)),': ' num2str(ts(end)), ' (Final Revolution)']);
0175                     end
0176                 end
0177             end
0178         end
0179         if (Verbose==1)
0180             disp(' ');
0181             disp(['--->Input Spike Time ', num2str(j), ' (Input Neuron ', num2str(ti(j,1)), '): ', num2str(ti(j,2))]);
0182             disp(['--->Pre -Spike Phase: ' num2str(thm(j))]);
0183             disp(['--->Post-Spike Phase: ' num2str(thp(j))]);
0184             disp(['--->Weight and Delay: ' num2str(Weights(ti(j,1),NeuronIndex)), ' ', num2str(Delays(ti(j,1),NeuronIndex))]);
0185             if thp(j)<ThFPP && thp(j)>-ThFPP && Io<0;
0186                 disp(['--->Operating Below Positive Fixed Point (',num2str(ThFPP),') - Neuron may not fire again']);
0187             end
0188         end
0189         
0190         prev_time=ti(j,2);
0191         prev_thp=thp(j);
0192 
0193     end
0194     if ~isempty(thp)
0195         c1=tan(thp(end)/2)/Beta;
0196     else
0197         c1=tan(CurrentNeuron.Phase/2)/Beta;
0198     end
0199     %If Io>0, check to see if an input spike has occured since the last
0200     %output spike, and if so calculate the next output spike.
0201     if (Io>0 && (isempty(ts) || ti(end,2)>ts(end)))
0202         if ~isempty(ti)
0203             ts(end+1)=ti(end,2)+(1/Beta)*( (pi/2) - atan(c1) );
0204         else
0205             ts(end+1)=(1/Beta)*( (pi/2) - atan(c1) );
0206         end
0207         if (Verbose==1)
0208             disp(['---> Neuron ', num2str(NeuronIndex),' Output Spike Time ',num2str(length(ts)),': ' num2str(ts(end)), ' (Next Cyclic Spike)']);
0209         end
0210 
0211         %If we are in region 1 with Io<0, calculate the final ts
0212     elseif ~isempty(thp) && (thp(end)>(ThFPP+1e-5)) && Io<0
0213         ts(end+1)=ti(end,2)+(1/Beta)*atanh(1/c1);
0214         if (Verbose==1)
0215             disp(['---> Neuron ', num2str(NeuronIndex),' Output Spike Time ',num2str(length(ts)),': ' num2str(ts(end)), ' (Final Trajectory Spike)']);
0216         end
0217         %If we are in region 1 with Io<0 and no input spikes, calculate the final ts
0218     elseif isempty(thp) && (CurrentNeuron.Phase>(ThFPP+1e-5)) && Io<0
0219         ts(end+1)=(1/Beta)*atanh(1/c1);
0220         if (Verbose==1)
0221             disp(['---> Neuron ', num2str(NeuronIndex),' Output Spike Time ',num2str(length(ts)),': ' num2str(ts(end)), ' (Final Trajectory Spike, No Input Spikes)']);
0222         end
0223     end
0224     %% End Analytical Integration %%
0225 
0226 % Numerical Integration %
0227 elseif strcmp(SimulationMethod,'Numerical')
0228     f_count=0;
0229     %Integrate long enough to include effects of all input spikes, and at
0230     %least one spike for positive Io case (pi/Beta is the spike period)
0231     if Io<0
0232         t_bl=(1/Beta)*atanh(Beta/tan((ThFPP+0.0001)/2));
0233         if ~isempty(ti)
0234             int_limit=(ti(end,2)+1.5*t_bl)/TimeStep;
0235         else
0236             int_limit=1.1*t_bl/TimeStep;
0237         end
0238     else
0239         t_period=pi/Beta;
0240         int_limit=(ti(end,2)+1.5*t_period)/TimeStep;
0241     end
0242     Theta=zeros(1,ceil(int_limit));
0243     Theta(1)=ThFPP+0.0001;
0244     for j=1:int_limit
0245 
0246         %Update the unforced response of phase based on Io
0247         Theta(j+1)=Theta(j)+TimeStep*((1-cos(Theta(j)))+Alpha*Io*(1+cos(Theta(j))));
0248 
0249         %Update the forced response of phase based on Isyn
0250         for k=(f_count+1):size(ti,1)
0251             if ((ti(k,2)>=j*TimeStep) && (ti(k,2)<(j+1)*TimeStep))
0252                 if (Verbose==1)
0253                     disp(' ');
0254                     disp(['--->Input Spike Time ', num2str(k), ' (Neuron ', num2str(ti(k,1)), '): ', num2str(ti(k,2))]);
0255                 end
0256                 thm(k)=Theta(j+1);
0257                 if (Verbose==1)
0258                     disp(['--->Pre -Spike Phase: ' num2str(thm(k))]);
0259                 end
0260                 Theta(j+1)=Theta(j+1)+Weights(ti(k,1),NeuronIndex)*Alpha*(1+cos(Theta(j+1)));
0261                 thp(k)=Theta(j+1);
0262                 if (Verbose==1)
0263                     disp(['--->Post-Spike Phase: ' num2str(thp(k))]);
0264                     disp(['--->Weight and Delay: ' num2str(Weights(ti(k,1),NeuronIndex)), ' ', num2str(Delays(ti(k,1),NeuronIndex))]);
0265                     if thp(k)<ThFPP && thp(k)>-ThFPP && Io<0;
0266                         disp(['--->Operating Below Positive Fixed Point (',num2str(ThFPP),') -  Neuron may not fire again']);
0267                     end
0268                 end
0269                 f_count=f_count+1;
0270             else
0271                 break;
0272             end
0273         end
0274 
0275         %Check to see if no further spikes are possible/desired
0276         if ((f_count==size(ti,1) && Theta(j+1)<ThFPP && Io<0) || (~isempty(ts) && ti(end,2)<ts(end) && Io>0))
0277             break;
0278         end
0279 
0280         %Check for output spike
0281         if Theta(j+1)>pi && Theta(j)<=pi
0282             ts(end+1)=j*TimeStep;
0283             if (Verbose==1)
0284                 disp(['---> Neuron ', num2str(NeuronIndex),' Output Spike Time ',num2str(length(ts)),': ' num2str(ts(end))]);
0285             end
0286             Theta(j+1)=Theta(j+1)-2*pi;
0287         end
0288 
0289         %Special Case of Refractory Firing
0290         if Theta(j+1)>-pi && Theta(j)<=-pi
0291             Theta(j+1)=Theta(j+1)+2*pi;
0292         end
0293     end
0294     Theta=Theta(1:j+1);
0295     %% End Numerical Integration %%
0296     
0297 else
0298     disp('Error in running networked theta neuron: Unknown theta neuron simulation method!');
0299     if (Verbose==1)
0300         disp('******** End function "run" for single neurons in a theta neuron network ********');
0301     end
0302     return;
0303 end
0304 
0305 %Print Final Firing Times
0306 if (Verbose==1) && ~isempty(ts)
0307     disp(' ');
0308     disp(['---> Neuron ', num2str(NeuronIndex),' Output Spike Times : ' num2str(ts)]);
0309 elseif (Verbose==1)
0310     disp(' ');
0311     disp(['---> Neuron ', num2str(NeuronIndex),' does not Fire after receiving all the input spikes!']);
0312 end
0313 
0314 %For Thomas' Method, Neuron is assumed to spike at the end of the trial
0315 if NIPS2007Gradient==1 && isempty(ts)
0316     ts=20;
0317 end
0318 
0319 
0320 %Calculate MIMO Y (n x r) and Z (q x r) if desired
0321 if GradFlag
0322 
0323     %Calculate all DthpDthpm1 terms (Equation 3.1.4 in NECO Paper, Although MIMO Version is Used Here)
0324     DthpDthpm1=zeros(1,length(thp));
0325     DthpDthpm1(1)=1;
0326     for m=2:length(thp)
0327         if Io<0
0328             c1=tan(thp(m-1)/2)/Beta;
0329             c2=tanh(-Beta*(ti(m,2)-ti(m-1,2)));
0330             DthpDthpm1(m)=(1-c2^2)*(1+(Beta*c1)^2)/((1+c1*c2)^2+(Beta*(c1+c2))^2);
0331 
0332             if (Verbose==1)
0333                 disp(' ');
0334                 disp(['Test Numerical Gradient for Effect of Spike ', num2str(m-1), ' on Phase at Spike ', num2str(m)]);
0335                 c1=tan(thp(m-1)/2)/Beta;
0336                 thm_test(2)=(1-Alpha*Weights(ti(m,1),NeuronIndex)*sin(thm(m)))*...
0337                     2*atan( Beta*(exp(-2*Beta*(ti(m,2)-ti(m-1,2)))*(1+c1)-(1-c1))...
0338                     /(exp(-2*Beta*(ti(m,2)-ti(m-1,2)))*(1+c1)+1-c1) );
0339                 c1=tan((thp(m-1)-0.001)/2)/Beta;
0340                 thm_test(1)=(1-Alpha*Weights(ti(m,1),NeuronIndex)*sin(thm(m)))*...
0341                     2*atan( Beta*(exp(-2*Beta*(ti(m,2)-ti(m-1,2)))*(1+c1)-(1-c1))...
0342                     /(exp(-2*Beta*(ti(m,2)-ti(m-1,2)))*(1+c1)+1-c1) );
0343                 c1=tan((thp(m-1)+0.001)/2)/Beta;
0344                 thm_test(3)=(1-Alpha*Weights(ti(m,1),NeuronIndex)*sin(thm(m)))*...
0345                     2*atan( Beta*(exp(-2*Beta*(ti(m,2)-ti(m-1,2)))*(1+c1)-(1-c1))...
0346                     /(exp(-2*Beta*(ti(m,2)-ti(m-1,2)))*(1+c1)+1-c1) );
0347                 df=gradient(thm_test,0.001);
0348                 disp('Numerical and Analytical Should Agree');
0349                 disp(['Numerical Value: ', num2str(df(2))]);
0350                 disp(['Analytical Value: ', num2str(DthpDthpm1(m)*(1-Alpha*Weights(ti(m,1),NeuronIndex)*sin(thm(m))))]);
0351             end
0352         else
0353             c1=tan(thp(m-1)/2)/Beta;
0354             c2=tan(Beta*(ti(m,2)-ti(m-1,2)));
0355             DthpDthpm1(m)=(1+c2^2)*(1+(Beta*c1)^2)/((1-c1*c2)^2+(Beta*(c1+c2))^2);
0356 
0357             if (Verbose==1)
0358                 disp(' ');
0359                 disp(['Test Numerical Gradient for Effect of Spike ', num2str(m-1), ' on Phase at Spike ', num2str(m)]);
0360                 c1=tan((thp(m-1)-0.001)/2)/Beta;
0361                 thm_test(1)=(1-Alpha*Weights(ti(m,1),NeuronIndex)*sin(thm(m)))*...
0362                     2*atan(Beta*((tan(Beta*(ti(m,2)-ti(m-1,2)))+c1)/(1-c1*tan(Beta*(ti(m,2)-ti(m-1,2))))));
0363                 c1=tan(thp(m-1)/2)/Beta;
0364                 thm_test(2)=(1-Alpha*Weights(ti(m,1),NeuronIndex)*sin(thm(m)))*...
0365                     2*atan(Beta*((tan(Beta*(ti(m,2)-ti(m-1,2)))+c1)/(1-c1*tan(Beta*(ti(m,2)-ti(m-1,2))))));
0366                 c1=tan((thp(m-1)+0.001)/2)/Beta;
0367                 thm_test(3)=(1-Alpha*Weights(ti(m,1),NeuronIndex)*sin(thm(m)))*...
0368                     2*atan(Beta*((tan(Beta*(ti(m,2)-ti(m-1,2)))+c1)/(1-c1*tan(Beta*(ti(m,2)-ti(m-1,2))))));
0369                 df=gradient(thm_test,0.001);
0370                 disp('Numerical and Analytical For DthpDthpm1 Should Agree');
0371                 disp(['Numerical Value: ', num2str(df(2))]);
0372                 disp(['Analytical Value: ', num2str(DthpDthpm1(m)*(1-Alpha*Weights(ti(m,1),NeuronIndex)*sin(thm(m))))]);
0373             end
0374 
0375         end
0376         DthpDthpm1(m)=DthpDthpm1(m)*(1-Alpha*Weights(ti(m,1),NeuronIndex)*sin(thm(m)));
0377     end
0378     
0379     %Initialize gradient training matrices
0380     %Y: an nxr matrix of the effect of each synapse n on output spike r
0381     %Z: a qxr matrix of the effect of each input spike q on output spike r
0382     %Y=zeros(max(ti(:,1)), length(ts)); %Will it be a problem if there are two synapses, but the one with connection index one does not fire?
0383     Y=zeros(size(Weights,1), length(ts)); %Will it be a problem if there are two synapses, but the one with connection index one does not fire?
0384     Z=zeros(size(ti,1), length(ts));
0385 
0386     %Cycle through each output spike, r
0387     for r=1:length(ts)
0388         if isempty(ti) break; end %For case where no input spikes are provided?
0389         last_fire_index=find(ti(:,2)<ts(r),1,'last');
0390         yD=(1-cos(thp(last_fire_index)) + Alpha*Io* (1+cos(thp(last_fire_index))));
0391 
0392         %Cycle through each input spike, q
0393         SpikeCount=ones(max(ti(:,1)),1);
0394         for q=1:last_fire_index
0395             n=ti(q,1);
0396 
0397             %Including the old gradient calculation method for comparison
0398             %purposes.  Note that it is only valid for networks with no
0399             %hidden layers, and probably on SISO data sets.
0400             if NIPS2007Gradient==1 %NIPS 2006 method
0401                 y_temp=((-Alpha*(1+cos(thm(q)))) / ((1-cos(thp(q)))+ Alpha*Io*(1+cos(thp(q)))));
0402             else %Regular method
0403                 y_temp=((-Alpha*(1+cos(thm(q)))) / yD) * prod(DthpDthpm1((q+1):last_fire_index));
0404             end
0405 
0406             %Verify first that there is a connection between n and NeuronIndex, otherwise there is no connection and y reduces to 0.
0407             if Weights(n,NeuronIndex)
0408                 Y(n,r)=Y(n,r)+y_temp;
0409             end
0410             SpikeCount(n)=SpikeCount(n)+1;
0411 
0412 
0413             gamma=(1/(Alpha*(1+cos(thm(q))))) * (-((1-cos(thp(q)))+Alpha*Io*(1+cos(thp(q)))) + ...
0414                 ((1-Alpha*Weights(n,NeuronIndex)*sin(thm(q))) * ((1-cos(thm(q)))+Alpha*Io*(1+cos(thm(q))))) );
0415             Z(q,r)=gamma*y_temp;
0416 
0417         end
0418         SpikeCount=max(SpikeCount-1,1);
0419 %        for sc=1:length(SpikeCount)
0420 %            SpikeCount(sc)=max(SpikeCount(sc)-1,1);
0421 %        end
0422         %Y(:,r)=Y(:,r)./SpikeCount;
0423     end
0424 
0425     %NIPS Credit Bound
0426     if NIPS2007Gradient==1 && min(min(Y))<-1000
0427         Y(find(Y<-1000))=-1000;
0428     end
0429 end
0430 
0431 %Verbose Footer
0432 if (Verbose==1)
0433     disp('******** End function "run" for single neurons in a theta neuron network ********');
0434 end
0435 return;