% eight_1.m      Problem 1 for Chapter 8 of Dayan and Abbott

% NEUR 1680

% Simulate the course of Hebbian learning for the case of figure 8.3.
% Find the ranges of initial weight values, (w_1, w_2), that lead to saturation at (1,1).
% Can you predict the result analytically?
% If the off-diagonal term in the correlation matrix is -2 instead of -0.4 and
% there are no saturation boundaries, what happens to the sum of the
% weights? Could this be used as a way of normalizing the weights?

% Last updated April 18, 2012

clear

% We will use equation 8.5 from the book, with the following parameters:


tau_w = 100;        % (arbitrary units)
ite = 1000;         % number of iterations 
N_runs = 100;       % number of runs, i.e., trajectories starting at random point
w_max = 1;          % saturation value         

q_s = 1;        % diagonal element of correlation matrix
q_d = -0.4;     % off-diagonal element of correlation matrix

Q = [q_s q_d; q_d q_s]; % correlation matrix

figure(1)
axis([0 w_max 0 w_max])
axis square
grid on
xlabel('W_1')
ylabel('W_2')

    
 for i_run = 1:N_runs
    
    W = zeros(2,ite);                   % this will be the 2-D weight vector  at times t = 1,2, ... ite
    w = w_max*rand(2,1);                % initial weight vector is picked at random in square

    for t = 1:ite                       % create trajectory
        
        w = w + Q*w/tau_w;              % equation 8.5
        w = min(w,[w_max;w_max]);       % saturation at w_max
        w = max(w,[0;0]);               % keep weights non-negative
        W(:,t) = w;                     % store weight vector
        
    end
    
    plot(W(1,:),W(2,:),'color','m','linewidth',1)       % plot trajectory in magenta
    hold on
    plot(W(1,end),W(2,end),'m.','markersize',20)        % mark end point in magenta

    % mark with a red point the initial states of trajectories ending near (w_max, w_max):
    if      (norm(W(:,end) - [w_max; w_max])< .1)
        plot(W(1,1),W(2,1),'r.','markersize',20)
        axis([0 w_max 0 w_max])
        axis square
        grid on
        drawnow
    end
    
    % mark with a blue point the initial states of trajectories ending near (0, w_max) or (w_max, 0):
    if  (norm(W(:,end )- [w_max; 0])< .1) || (norm(W(:,end )- [0; w_max])< .1)
        plot(W(1,1),W(2,1),'b.','markersize',20)
        drawnow
    end
    
 end
 
 
% To derive the curve that separates the attraction basin of (1, 0) from that of (1, 1),
% we first note that all trajectories are of the form:
%   y_1 = C_1 exp(.6t)
%   y_2 = C_2 exp(1.4t)
% where:
%   y1 = w1 + w2
%   y2 = w1 - w2.
% To find C_1 and C-2, we note that we want w2' = 0 when w1 = 1.
% Putting it all together, we derive the following implicit form for the trajectory:
%   ((w1-w2)/.6)^.6 - ((w1+w2)/1.4)^1.4 = 0


% To plot this trajectory, we use the command ezplot:
h = ezplot('((w1-w2)/.6)^.6 - ((w1+w2)/1.4)^1.4',[0 w_max 0 w_max]);
set(h, 'linewidth',3, 'color','k')

% To plot the curve that separates the attraction basin of (0, 1) from that of (1, 1)
%    we switch the roles of w1 and w2
h = ezplot('((w2-w1)/.6)^.6 - ((w2+w1)/1.4)^1.4',[0 w_max 0 w_max]);
set(h, 'linewidth',3, 'color','k')

title(['BASIC HEBB WITH SATURATION AT  0  AND  ' num2str(w_max)])
hold off

[V,D] = eig(Q);         % eigenvalues and eigenvectors of correlation matrix
disp(' ')
disp('eigenvalues:')
disp(num2str(D(logical(eye(2)))'))
disp('eigenvectors:')
disp(num2str(V))

