%% IZFORCE CLASSIFIER WITH ADAPTATION 
clear all
clc 

T = 1000000; %Total time 
dt = 0.04; %integration time step 
nt = round(T/dt);
N =  2000;  %number of neurons 
%% Izhikevich Parameters
C = 250;
vr = -60; 
b = 0; 
k = 2.5; 
vpeak = 30; 
vreset = -65;
vt = vr+40-(b/k); %threshold 
Er = 0; %Reversal Potential 
u = zeros(N,1); 
a = 0.01;
d = 200; 
tr = 2; 
td = 20; 
p = 0.1; %sparsity 
G =6*10^3; %scale weight matrix  
BIAS = 1000; %Bias current
OMEGA =  G*(randn(N,N)).*(rand(N,N)<p)/(p*sqrt(N));
%% Initialize currents, FORCE method, other parameters
IPSC = zeros(N,1); %post synaptic current 
h = zeros(N,1);
r = zeros(N,1);
hr = zeros(N,1);
JD = zeros(N,1);
BPhi = zeros(N,1);

%-----Initialization---------------------------------------------
v = vr+(vpeak-vr)*rand(N,1); %initial distribution 
v_ = v; %These are just used for Euler integration, previous time step storage
IPSC = zeros(N,1);
Q = 5*10^3; %Scale feedback term, Q in paper
E = (2*rand(N,1)-1)*Q; %scale feedback term
WE2 = 5*10^2; %scale input weights
Psi = 2*pi*rand(N,1); 
Ein = [cos(Psi),sin(Psi)]*WE2;
load lineardata.mat;  
%% lineardata.mat corresponds to linear boundaries while nonlineardata.mat corresponds to nonlinear boundaries 

%Load the data set for classification.  P contains class, z contains the points.  
inputfreq = 4; %Present a data point; 
j =0;
Xin = zeros(nt,2); 
zx = zeros(nt,1); 
nx = round(1000/(dt*inputfreq)); 
zx = 0.5*abs(sin(2*pi*(1:1:nt)*dt*inputfreq/2000)); 


%% Create supervisor and inputs.  
j = 1; 
k2 = 0;
for i =1:1:nt 
 zx(i) = zx(i)*(2*P(j)-1)*mod(k2,2);
 Xin(i,:) = [z(j,1),z(j,2)]*mod(k2,2);
if mod(i,nx)==1 
    k2 = k2 + 1; 
j = ceil(rand*2000);
end
end
clear z P

%% initialization 
z = 0; 
tspike = zeros(nt,2);
ns = 0;
%% RLS parameters.  
 Pinv = eye(N)*30;
 step = 20;
 imin = round(1000/dt);
 icrit = round(400000/dt);
 current = zeros(nt,1); 
 RECB = zeros(nt,5);
 REC = zeros(nt,10);
 i=1;
%% SIMULATION
tic
ilast = i ;
for i = ilast:1:nt; 
%% EULER INTEGRATE
I = IPSC + E*z + Ein*(Xin(i,:)') + BIAS; 
v = v + dt*(( k.*(v-vr).*(v-vt) - u + I))/C ; % v(t) = v(t-1)+dt*v'(t-1)
u = u + dt*(a*(b*(v_-vr)-u)); %same with u, the v_ term makes it so that the integration of u uses v(t-1), instead of the updated v(t)

%% 
index = find(v>=vpeak);
if length(index)>0
JD = sum(OMEGA(:,index),2); %compute the increase in current due to spiking  
%tspike(ns+1:ns+length(index),:) = [index,0*index+dt*i];
ns = ns + length(index); 
end


if tr == 0 
    IPSC = IPSC*exp(-dt/td)+   JD*(length(index)>0)/(td);
    r = r *exp(-dt/td) + (v>=vpeak)/td;
else
IPSC = IPSC*exp(-dt/tr) + h*dt;
h = h*exp(-dt/td) + JD*(length(index)>0)/(tr*td);  %Integrate the current

r = r*exp(-dt/tr) + hr*dt; 
hr = hr*exp(-dt/td) + (v>=vpeak)/(tr*td);
end


%% Apply RLS 
 z = BPhi'*r;
 err = z - zx(i);
if mod(i,step)==1
if i > imin 
 if i < icrit 
   cd = Pinv*r;
   BPhi = BPhi - (cd*err');
   Pinv = Pinv -((cd)*(cd'))/( 1 + (r')*(cd));
 end 
end 
end




%% COMPUTE S, APPLY RESETS
u = u + d*(v>=vpeak);  %implements set u to u+d if v>vpeak, component by component. 
v = v+(vreset-v).*(v>=vpeak); %implements v = c if v>vpeak add 0 if false, add c-v if true, v+c-v = c
v_ = v;  % sets v(t-1) = v for the next itteration of loop
REC(i,:) = [v(1:5)',u(1:5)']; 
current(i,:) = z; 
RECB(i,:)=BPhi(1:5);
if mod(i,round(100/dt))==1 
drawnow

figure(2)
plot(dt*(1:1:i),current(1:1:i,1)), hold on 
plot(dt*(1:1:i),zx(1:1:i),'k--'), hold off
ylim([-0.6,0.6])
xlim([dt*i-1000,dt*i])
xlabel('Time')
ylabel('Network Response')
legend('Network Output','Target Signal')


end   
end
%% 