%Implements a simple network and plots it dynamically (useful when you
%don't understand what's going on)

clear; %<= Never forget this one!

%0. Basic Setup

xsize = 10; %Size of input
minV  =  1; %Minimum input value
maxV  =  5; %Maximum input value
A = 0.5;
w = rand(10);
w = w / sum(sum(w));    
input = [minV*ones(xsize/2,1); maxV*ones(xsize/2,1)];

x = zeros(xsize,1);     %First layer
y = zeros(xsize,1);     %Second layer

dt    = 0.05; %step for Euler's method
tmax  = 2; %Simulation time, in milliseconds
index = 2; %Current position in array

%1. Set up Figure
f = figure(1);
set(f,'position',[1 31 1280 920]);

s(1) = subplot(3,1,2);
s(2) = subplot(3,1,1);

%Ylims
xMax = maxV/A;      
yMax = xMax/7/A;
set(s(1),'ylim',[0 xMax]);
set(s(2),'ylim',[0 yMax]);

s(3) = subplot(3,1,3);
plot(input,'linewidth',2);
xlabel('Unit','fontsize',25);
ylabel('Luminance (cd/m^2)','fontsize',15);

axes(s(1));
title('x','fontsize',25);
ylabel(['Firing rate (Hz)       ';'Membrane potential (mV)'],'fontsize',15);

axes(s(2));
title('y','fontsize',25);
ylabel('Activation (?)','fontsize',15);

set(s,'fontsize',30,'xlim',[0.5 xsize + 0.5]);

%2. Run simulation
for i=dt:dt:tmax

    %Integration step
    x(:,index) = x(:,index-1) + dt * ( -A*x(:,index-1) + input );	%x receives one-to-one input
    y(:,index) = y(:,index-1) + dt * ( -A*y(:,index-1) + w*x(:,index-1) ); %y receives weighted summation of x

    %Update plots
    subplot(3,1,2); hold on;
    cla;
    p(1) = plot(x(:,index),'o-');

    subplot(3,1,1); hold on;
    cla;    
    p(2) = plot(y(:,index),'o-');

    set(p,'linewidth',2);

    pause(0.1);
  
    %Advance index
    index = index+1;
    
end

%3. Alternative way to plot things in time
f2 = figure(2);
s2(1) = subplot(2,1,2);
mesh(x);
title('y','fontsize',25);
ylabel('Unit','fontsize',15);
xlabel('Time (ms)','fontsize',15);
zlabel('Activation (?)','fontsize',15);

s2(2) = subplot(2,1,1);
mesh(y);
title('x','fontsize',25);
set(s2,'fontsize',15);