% eight_4.m      Problem 4 of Chapter 8 of Dayan and Abbott

% NEUR 1680

% Construct two-dimensional input data sets similar to those shown
% in figure 8.4 and use them to train a two-input, one-output, linear
% network using correlation- and covariance-based Hebbian learning
% rules with multiplicative normalization. Compare the final outcome
% for the weights with the principal components of the data when the
% mean of the input distribution is zero and when it is nonzero.

% 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)
alpha_coeff = .1;   % coefficient of multiplicative-normalization term (Oja term)
ite = 30000;        % 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

rule = 'corr';       % choose 'corr' for correlation-based rule 
                    %        'cov'  for covariance-based rule

m1 = 1;             % mean of input 1 (firing rate of presynaptic neuron 1)
m2 = 1;             % mean of input 2 (firing rate of presynaptic neuron 2)
var1 = 1;           % variance of input 1 
var2 = 1;           % variance of input 2
cov12 = -.8;        % covariance of inputs 1 and 2


m = [m1; m2];       % mean input will be used in covariance-based rule
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([m1 m2],C,ite)';             % input sequence is drawn from bivariate gaussian

    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  = 10*rand(2,1) - [5; 5];             % initial weight vector is picked at random in large square
    
    for t = 1:ite                           % create trajectory
        
        u = U(:,t);                         % presynaptic activity at time t
        v = w'*u;                           % postsynaptic activity (equation 8.2)
       
        switch rule
            
            case 'corr'                     % correlation rule (equation 8.3)
            
            w = w + v*u/tau_w;      
            
            case 'cov'                      % covariance rule with presynaptic threshold (equation 8.9)
                                            % mean input is used as presynaptic threshold
            w = w + v*(u - m)/tau_w;        
            
        end
        
        w = w - alpha_coeff*v^2*w/tau_w;    % subtract "multiplicative" normalization term (Oja rule, equation 8.16)
        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([-5 5 -5 5])
        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')
switch rule
    case 'corr'
        title('Correlation-based Hebbian learning with multiplicative (Oja) weight normalization')
    case 'cov'
        title('Covariance-based Hebbian learning with multiplicative (Oja) weight normalization')
end
hold off

% plot the responses for 5 of the input vectors, as well as the square amplitude of the weight vector
response_curves = W'*U(:,1:5);
figure(2)
plot(response_curves)
hold on
plot(sum(W.^2), 'linewidth',3)
hold off
xlabel('time')
legend('response curve','response curve','response curve','response curve','response curve','squared norm of weight vector')
grid on
