%let's calculate the rates for a set of M neurons and test them out with
%some stimuli

M=30 % number of neurons
x=10 % stimulus value
g=10 % stimulus strength (contrast)
sig=3 % tolerance of neuron
numberofevents=300 % number of detector events
stimuluschange=1 % setting this to 1 switches on an abrupt change to the stimulus.


resprate=@(x,g,sig,mu) g*exp(-((x-mu)^2)/(2*sig^2)) % returns the firing response rate for a single neuron with preferred stimulus mu, tolerance sig, input x and strength of input g.

rates=[]; % the below calculates the rates
for i=1:M
    rates=[rates resprate(x,g,sig,i)]; % this calculates all of the rates for the M neurons
end

%% crude Gillespie
who=[];
when=[];
% the below works out when the events occurr and to which neurons.
for T=1:numberofevents % moving to event-time rather than absolute time

az0=sum(rates)
r1 = rand(); 
r2 = rand(); 
tau = (1/az0)*log(1/r1)

 % locate neuron corresponding to draw
      tmp = r2*az0; 
      cumsuma = cumsum(rates);
      reaction = -1;
      for j=1:M
          if cumsuma(j) > tmp && reaction == -1
              reaction = j;
          end
      end
who=[who reaction]
when=[when tau]
end

if stimuluschange==1

%now we will simulate an abrupt change in the location of the stimulus 
x=2 % stimulus value
g=100 % stimulus strength (contrast)

rates=[];
for i=1:30
    rates=[rates resprate(x,g,sig,i)];
end

for T=1:numberofevents % moving to event-time rather than absolute time

az0=sum(rates)
r1 = rand(); 
r2 = rand(); 
tau = (1/az0)*log(1/r1)

 % locate neuron corresponding to draw
      tmp = r2*az0; 
      cumsuma = cumsum(rates);
      reaction = -1;
      for j=1:M
          if cumsuma(j) > tmp && reaction == -1
              reaction = j;
          end
      end
who=[who reaction]
when=[when tau]
end

end

when=cumsum(when)

figure; scatter(who,when)

figure; hist(who)



