%%  ART1 network for letter recognition
clear;
% Parameters & initial values
n1=12*13; n2=26; % number of input and output nodes
a1=1;  b1=1.5;  c1=1;  d1=1; L=2; rho=0.9;
wt=2*(b1-1)/d1*ones(n1,n2); wb=0.5*L/(L-1+n1)*ones(n2,n1);
x1=-b1/(1+c1)*ones(n2,1); x2=zeros(n2,1);
%% training vectors
load pattern1; rPat=reshape(pattern1', n1, n2); % matrix of inputs

%% ART update
for run=1:100;
    for pat=1:n2
        I=rPat(:,pat);
        x=randperm(n1); I(x(1:16))=1-I(x(1:16)); %adding noise
        IOR=ones(n2,1);
        
        while sum(IOR)>0;
            x1=I./(1+a1*(I+b1)+c1);  % activ. of F1 with ext input
            s1=(x1>0)*1;             % output of F1
            t=wb*s1;                 % input to F2
            [maxv,maxi]=max(IOR.*t); % winner
            u=zeros(n2,1); u(maxi)=1;% takes all
            v=wt*u;                  % top-down input to F1
            x1=(I+d1*v-b1)./(1+a1*(I+d1*v)+c1);% new activ. of F1
            s1=(x1>0)*1;             % new output of F1
            
            if sum(s1)/sum(I) > rho; %resonant state
                wt(:,maxi)=wt(:,maxi)+0.1*(s1-wt(:,maxi));
                wb=wt'*L/(L-1+sum(s1));
                break
            else
                IOR(maxi)=0; % search for other node
            end
        end
    end
end
%% Display catergory prototypes
for i=1:n2
    thisWeightArray=reshape(wt(:,i),13, 12)';
    format +; disp(thisWeightArray-0.5); format;
end