

%Adapting Integrate and fire from Christoff Koch chapter 14
%voltages in millivolts, currents in nanoamps, 
% resistances in megaohms
% capacitance in nanoFarads
%
%fully vectorized version

clear all
global dVictor ngrp grp dRossum tSpike

%Number of neurons
nNeuron = 20;
%timing
dt = .05; %millisecond time step
tEnd = 150; %maximum simulation time
%build the time vector
time = 0:dt:tEnd ;


%membrane parameters
tMemb = 10*ones(1,nNeuron);  %mSec time constant, was 3
rMemb = 40*ones(1,nNeuron); %memb resistance
cMemb = tMemb/rMemb *ones(1,nNeuron); %memb capacitance
vRest = 0;  %millivolts resting potential
tRef = 2*ones(1,nNeuron); %mSec: length of the absolute refactory period. v is held at zero for this time

vThresh = 16*ones(1,nNeuron); %threshold voltage

%adapting (K) cond: gAdapt
tAdapt = 20*ones(1,nNeuron); %mSec time constant of slow recovery after an AP 
Ginc = .0/dt*ones(1,nNeuron);  %.05 microSiemans: Instant increase in gAdapt when a spike occurs 

%synaptic connections
%weight matrix --
%each column represents the input weights to one neuron
%e.g. element(1,2) is the weight connecting neuron 1 to neuron 2

weight = zeros(nNeuron,nNeuron);

%weight = 3 * randn(nNeuron,nNeuron) + .0;
%weight = 2 * rand(nNeuron,nNeuron) - .25;

%synaptic current rises instantly, then decays exponentially
%decay rate of synaptic current 
%Irate = .3*ones(1,nNeuron);
%Idecay = 1 - Irate*dt;

%init state variables
vMemb = zeros(1,nNeuron);
spike = zeros(1,nNeuron);
gAdapt = zeros(1,nNeuron);
spikeTime = -inf*ones(1,nNeuron); %time of occurence of last spike

%initial input currents
%sine currents added in
Io = 0.4*ones(1,nNeuron); %0.4
%DC levels 
I1 = 1.;
I2 = 1.;
I3 = 1.;
I4 = 1.;
Is = [I1*ones(1,nNeuron/4) I2*ones(1,nNeuron/4)...
    I3*ones(1,nNeuron/4)  I4*ones(1,nNeuron/4)]; %steady current

AmpSine = [.2*ones(1,nNeuron/4) .4*ones(1,nNeuron/4)...
    .6*ones(1,nNeuron/4)  .8*ones(1,nNeuron/4)]; %sine amplitude
f = 0.1; %cycles/mSec

%init the noise state
Ir = zeros(1,nNeuron);
%step the time
tindex = 1; %pointer for output arrays

%init plotting variables
Vout = zeros(length(time),nNeuron) ;  
Aout = zeros(length(time),nNeuron) ;
Sout = zeros(length(time),nNeuron) ;


for t=time
    
    %input current to cells
    %Istate =  spike*weight + Istate .* Idecay ; %decaying synaptic currents
    Ir = 0.1*randn(size(Is)) + .90*Ir ; %rms amp approx 0.23
    %Ir = randn(size(Is));
    I =  (Is+AmpSine*sin(2*pi*f*t)) .* (t>50 & t<100) + .25*Ir + Io ;
    
    %add to plotting variables
    Vout(tindex,:) = vMemb;  
    Aout(tindex,:) = I ;
    Sout(tindex,:) = spike;
    
    
    %update state variables
    %use second order Runge-Kutta
    % dgAdapt/dt = (Ginc.*spike - gAdapt)./tAdapt
    % dv/dt = 
    
    k1g = dt*((Ginc.*spike - gAdapt./tAdapt)) ;
    
    k1v = dt*( I./cMemb - vMemb.*(1 + rMemb.*gAdapt)./tMemb) .* ...
        ((t-spikeTime) >= tRef) ;
    
    k2g = dt*(Ginc.*spike - (gAdapt + 0.5*k1g)./tAdapt) ;
    k2v = dt*( I./cMemb - (vMemb + 0.5*k1v).*(1 + rMemb.*(gAdapt + 0.5*k1g))./tMemb) .* ...
        ((t+0.5*dt-spikeTime) >= tRef) ;
    
    
    gAdapt = gAdapt + k2g;
    vMemb = vMemb + k2v;
    
    %generate the spike
    spike = (vMemb >= vThresh);
    %and its time of occurence
    spikeTime(spike) = t;
    %hit the membrane potential
    vMemb(spike) = 0.;
    
    %update time index
    Tout(tindex) = t;
    tindex = tindex+1;
    
end %main time step for-loop

figure(1)
clf
np=3;

%subplot(np,1,1)
%plot(Tout,Vout)
%ylabel('V')

subplot(np,1,1)
plot(Tout,Aout)
ylabel('Iin')

subplot(np,1,2)
plot(Tout,Sout(:,1:nNeuron/4))
hold on
plot(Tout,Sout(:,nNeuron/4:nNeuron/2)-.8)
plot(Tout,Sout(:,nNeuron/2:nNeuron*3/4)-1.6)
plot(Tout,Sout(:,nNeuron*3/4:nNeuron)-2.4)
ylabel('Spikes')
%set(gca,'ylim',[-1 1.5])
xlabel('time (mS)')
drawnow
%return 

%figure(2)
%clf
%pn=2;
%subplot(np,1,3)

%compute times of spikes in seconds into a CELL ARRAY
spiketrain = Sout(fix(50/dt):fix(100/dt),:);
for i=1:nNeuron
    tSpike{i} = .001*time(find(spiketrain(:,i)>0.5))+.05 ;
end

srange = 1:8;
for s=srange
    TimeScale(s) = .001 * 2^(s-1);
    for i=1:nNeuron
        for j=i+1:nNeuron
            %calc vanRossum dist
            %make dt depend on timescale
            tc = TimeScale(s) ; % Sec
            dt = tc/100 ;
            tk = 0:dt:4*tc;
            %convert tend to seconds
            t = 0:dt:tEnd*.001;
            
            spiketrain1 = zeros(size(t));
            spiketrain2 = zeros(size(t));
            spiketrain1(fix(tSpike{i}/dt)) = 1;
            spiketrain2(fix(tSpike{j}/dt)) = 1;
            
            ekernel = exp(-tk/tc);
            cSpike1 = conv(ekernel,spiketrain1);
            cSpike2 = conv(ekernel,spiketrain2);
            dRossum(s,i,j) =  2*dt/tc * sum((cSpike1-cSpike2).^2) ;
            dRossum(s,j,i) = dRossum(s,i,j);
        end  
    end
end


%pick a distance weighting exponent
k = -2 ;

%get bit rate for each timescale
for s=srange
    
    %define response group membership as a CELL array
    %( note curly brackets{}! )
    %so that different groups can have different numbers of traces
    ngrp = 4;
    grp{1} = 1:5  ; 
    grp{2} = 6:10 ;
    grp{3} = 11:15 ;
    grp{4} = 16:20 ;
    
    %make array of avg dist of each train to each group 
    %as an (nNeuron,ngrp) array
    %excluding the current train from the average
    for i=1:nNeuron
        for j = 1:ngrp
            traceIndex = grp{j}(find(grp{j}~=i)) ;  %excludes the current trace
            avgDist(i,j) = ( mean(dRossum(s,i,traceIndex).^k) ).^(1/k) ;
        end
    end
    
    %find minimum in each row (for each train)
    %make N matrix (see Victor)
    %with a row per stimulus group and a column per response group
    N=zeros(ngrp);
    for i=1:nNeuron
        [minv,minResponseGrp] = min(avgDist(i,:));
        %find out which stim group the trace is in
        for j=1:ngrp
            if (sum(grp{j}==i)==1)
                StimGrp = j;
            end
        end
        %form the confusion matrix
        N(StimGrp,minResponseGrp) = N(StimGrp,minResponseGrp) + 1;
    end
    
    
    %get the number of bits available--
    %max info is log2(ngrp)
    Ntotal = sum(sum(N)); 
    
    for i=1:ngrp
        for j=1:ngrp
            h(i,j) = N(i,j) * ...
                ( log2(N(i,j)+.0001) - log2(sum(N(:,j))+.0001) - ...
                log2(sum(N(i,:))+.0001) + log2(Ntotal) );     
        end
    end
    
    H(s) = sum(sum(h))/Ntotal  ;
end    

subplot(np,1,3)
semilogx(TimeScale,H)
xlabel('Time scale (sec)')
ylabel('Bits')
set(gca,'ylim',[0 2])
%text(.01,1,num2str(N))
%legend('Victor Dist','Rossum Dist')







