% Chad Myers
% ELE 571
% File: xortrain.m
% ------------------------
% This file contains the xortrain function which, given a number of clusters,
% number of training points, number of iterations, and a flag specifying
% whether to use RBF or EBF distrubutions trains a neural network to
% discriminate between two XOR-like classes.


function [W,means,sigma,correct] = xortrain(clusters,n,iter,rbf)  % xortrain(num. of clusters, num. of training pts.,num. of training iterations,1 for RBF/0 for EBF)
n;  %number of training samples
iter; %number of training iterations
clusters;  %number of clusters
rbf; %0-EBF 1-RBF
class1_id = 1;
class2_id = -1;

W = ones(clusters,1);
if(rbf==1) 
    learnrate_w = .01;
    learnrate_c = .008;
    learnrate_sigma = .002;
else
    learnrate_w = .008;
    learnrate_sigma = .002;
    learnrate_c = .005;
end

fprintf('Training...\n');
%randomly generate training data and desired result in the window given by x E [-1,1], y E
%[-1,1]

for i=1:n,
    temp1 = 2*rand(2);
    temp2 = rand(2);
    
    if(temp1(1) > 1) data(i,1) = temp1(2);
    else data(i,1) = -temp1(2);
    end
        
    if(temp2(1) > .5) data(i,2) = temp2(2);
    else data(i,2) = -temp2(2);
    end
    
    if(abs(data(i,2)) > .5*abs(data(i,1))) desired(i,1) = class1_id;
    else desired(i,1) = class2_id;
    end
end


%use k-means and k-nearest neighbors to get initial means and variances
[temp1,means,sums] = kmeans(data,clusters,'dist','sqEuclidean');
    
%for RBF case
if(rbf == 1)
    for i=1:clusters,
    p=1;
    for j=1:n,
        if(temp1(j) == i) clusterdata(p,1:2) = data(j,1:2);
        p=p+1;
        end
    end
    
    variance = cov(clusterdata);
    sigma(i) = (variance(1,1)+variance(2,2))^.5;
    
    clear clusterdata;
    end

else %for EBF case
    for i=1:clusters,
        sum_x=0;
        sum_y=0;
        for j=1:clusters,
            sum_x = sum_x + (means(i,1)-means(j,1))^2;
            sum_y = sum_y + (means(i,2)-means(j,2))^2;
        end
        sigma(i,1) = (sum_x/clusters)^.5;
        sigma(i,2) = (sum_y/clusters)^.5;
    end
end


 %compute activation function outputs for training samples
 for j=1:n,
     for k=1:clusters,
         if(rbf==1)
             phi(k,j)= Gaussian(data(j,1:2),means(k,1:2),[sigma(k),sigma(k)]);
         else
             phi(k,j)= Gaussian(data(j,1:2),means(k,1:2),sigma(k,1:2));
         end
     end
 end
 
 %use least squares method to approximate W
 A = ones(n,clusters+1);
 A(1:n,2:(clusters+1)) = phi';
 W = inv(A'*A)*A'*desired;

figure;
%start training network
for i=1:iter    
    %recompute activation function outputs for training samples
    for j=1:n,
      for k=1:clusters,
           if(rbf==1)
               phi(k,j)= Gaussian(data(j,1:2),means(k,1:2),[sigma(k),sigma(k)]);
           else
                phi(k,j)= Gaussian(data(j,1:2),means(k,1:2),sigma(k,1:2));
            end
         end
    end
    
    %compute network output for training samples
    W_aug = W(2:(clusters+1),1);
    actual(:,i)= W(1,1) + phi'*W_aug;
    %compute error
    e(1:n,i) = (desired - actual(1:n,i));
    
    %compute new weights using previous values
    W_aug = W(2:(clusters+1),1);
    W(1,1) = W(1,1) + learnrate_w*(e(1:n,i)'*ones(n,1));
    W_aug = W_aug + learnrate_w*(phi*e(1:n,i));
    W(2:(clusters+1),1) = W_aug;
    
    %replot each iteration
    s = sign(actual(:,i));
    cla reset;
    scatter(data(:,1),data(:,2),10,s);
    hold on;
    
    x=-2:1:2;
    y=-x/2;
    plot(x,y);
    x=-2:1:2;
    y=x/2;
    plot(x,y);
    
    %plot ellipses for each cluster
    for q=1:clusters,
        if(rbf==1) dev = [sigma(q), sigma(q)];
        else dev = sigma(q,1:2);
        end
        ellipse = generate_ellipse(300,0,means(q,1:2),dev);
        plot(ellipse(:,1),ellipse(:,2),'bd','MarkerSize',1);
    end
    pause(.1);

    
    
    %adjust distributions
    if(rbf==1)  % for RBF case
        for m=1:clusters,
            s_sum=0;
            c_sum_x=0;
            c_sum_y=0;
            for p=1:n,
                s_sum = s_sum + e(p,i)*W_aug(m,1)*phi(m,p)*(norm(data(p,1:2)-means(m,1:2)))^2*sigma(m)^(-3);
                c_sum_x = c_sum_x + e(p,i)*W_aug(m,1)*phi(m,p)*(data(p,1)-means(m,1))*sigma(m)^(-2);
                c_sum_y = c_sum_y + e(p,i)*W_aug(m,1)*phi(m,p)*(data(p,2)-means(m,2))*sigma(m)^(-2);    
            end
            sigma(m) = sigma(m) + learnrate_sigma*s_sum;
            means(m,1:2) = means(m,1:2) + learnrate_c*[c_sum_x, c_sum_y];
        end
    else   % for EBF case
         for m=1:clusters,
            sum_x=0;
            sum_y=0;
            c_sum_x=0;
            c_sum_y=0;
            
            for p=1:n,
                sum_x = sum_x + e(p,i)*W_aug(m,1)*phi(m,p)*(data(p,1)-means(m,1))^2*sigma(m,1)^(-3);
                sum_y = sum_y + e(p,i)*W_aug(m,1)*phi(m,p)*(data(p,2)-means(m,2))^2*sigma(m,2)^(-3);
                c_sum_x = c_sum_x + e(p,i)*W_aug(m,1)*phi(m,p)*(data(p,1)-means(m,1))*sigma(m,1)^(-2);
                c_sum_y = c_sum_y + e(p,i)*W_aug(m,1)*phi(m,p)*(data(p,2)-means(m,2))*sigma(m,2)^(-2);  
            end
            sigma(m,1) = sigma(m,1) + learnrate_sigma*sum_x;
            sigma(m,2) = sigma(m,2) + learnrate_sigma*sum_y;
            means(m,1:2) = means(m,1:2) + learnrate_c*[c_sum_x, c_sum_y];
        end
    end
end

            
for i=1:iter,
    error(i) = norm(e(:,i));
end
figure;
plot([1:1:iter],error);

b=sign(actual(:,iter))+sign(desired);
a=abs(b(:,1));
correct=sum(a)/(2*n);
fprintf('Number of training points: %d\n',n);
fprintf('Percentage of training points correctly classified: %f\n',correct*100);
