% eight_5.m      Problem 5 of Chapter 8 of Dayan and Abbott

% NEUR 1680

% Repeat exercise 4 for a data set with zero mean, but with subtractive
% normalization and saturation. Start with initial values for the
% weights that are chosen randomly over the full range from 0 to their
% saturation limit. When does this algorithm produce a weight vector
% aligned with the principal component axis of the input data set, and
% when does it fail to do so? Why does the weight vector sometimes
% fail to align with the principal component axis?

% Again, our solution will use the sequential (i.e., non-averaged) version of Hebb rule (equations 8.3 and 8.9)

% Last updated April 18, 2012

clear

tau_w = 2000;       % (arbitrary units)
ite = 3000;        % number of iterations (equal to the number of inputs presented to create one trajectory)
N_runs = 50;        % number of runs, i.e., trajectories starting at random point
w_max = 5;          % weight-saturation value

var1 = 1;           % variance of input 1 
var2 = .2;           % variance of input 2
cov12 = .0;        % covariance of inputs 1 and 2

C = [var1 cov12; cov12 var2];

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

for i_run = 1:N_runs 

    U = mvnrnd([0 0],C,ite)';               % input sequence is drawn from bivariate gaussian with mean [0 0]

    W       = zeros(2,ite);                 % this will be the 2-D weight vector  at times t = 1,2, ... ite
    V       = zeros(1,ite);                 % this will be the response           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
        
        u = U(:,t);                         % presynaptic activity at time t
        v = w'*u;                           % postsynaptic activity (equation 8.2)
        saturation_margin = 0.01;           % weights close enough to the upper or lower bound
                                            %    will be considered saturated
        n = (w > saturation_margin) .* (w < (w_max - saturation_margin));
        N_u = sum(n);
        w = w + v*(u - (n'*u)*n/N_u)/tau_w; % subtractive weight normalization term (equation 8.14)
                                           	% modified to take saturation into account
        w = min(w,[w_max;w_max]);           % saturation
        w = max(w,[0;0]);                   % keep weights non-negative
        W(:, t) = w;                        % store weight vector
        
    end
    
    figure(1)
    
    if i_run == 1
        plot(U(1,:),U(2,:),'g.','markersize',5)    % on first run, create scatter plot of input
        axis([-3 w_max -3 w_max])
        axis square
        hold on
        grid on
    end
    
    plot(W(1,:),W(2,:),'color','b','linewidth',1)
    plot(W(1,1),W(2,1),'b.','markersize',20)
    plot(W(1,end),W(2,end),'r.','markersize',20)

end

xlabel('w_1   and   u_1')
ylabel('w_2   and   u_2')
title('Correlation-based Hebbian learning with subtractive weight normalization')
hold off


