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

%    LINES 36, 37 AND 38 SHOULD BE COMPLETED.

% NEUR 1680
% Last updated April 4, 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 = ;               % equation 8.5
        w = ;               % saturation at w_max
        w = ;               % 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
 

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))

