////  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
pattern1=read('pattern1',26*12,13);
rPat=matrix(pattern1', n1, n2); // matrix of inputs

//// ART update
for run=1:100;
    for pat=1:n2
        I=rPat(:,pat);
        x=grand(1,'prm',(1:156)'); 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,maxind]=max(IOR.*t); // winner
            u=zeros(n2,1); u(maxind)=1;// takes all
            v=wt*u;                  // top-down input to F1
            x1=(I+d1*v-b1)./(1+a1*(I+d1*v)+c1);// new activ.n of F1
            s1=(x1>0)*1;             // new output of F1
            
            if sum(s1)/sum(I) > rho; //resonant state
                wt(:,maxind)=wt(:,maxind)+0.1*(s1-wt(:,maxind));
                wb=wt'*L/(L-1+sum(s1));
                break
            else
                IOR(maxind)=0; // search for other node
            end
        end
    end
end
//// Display catergory prototypes
for i=1:n2
    thisWeightArray=matrix(wt(:,i),13, 12)';
  for j=1:12
    for k=1:13 
      thisWeightArrayDisplay(j,k)=' ';
      if (thisWeightArray(j,k)-0.5)>0
         thisWeightArrayDisplay(j,k)='+';
      end
    end
  end
  disp(thisWeightArrayDisplay)
end