[WX,ZX,O,E]=msnn2(t,xcell,noHidCells,lrate,epoch)

function [WX,ZX,O,E]=msnn2(t,xcell,noHidCells,lrate,epoch)

% This is a msnn backpropagation algorithm with 1 hiddenlayer
% t = target
% x = input
% noHidCells = number of hidden cells
% lrate = learning rate
% epoch = number of iterations
% W = weights for hidden layer
% Z = weights for output layer
% O = output response
% E = training error
% The training set consists of K pairs
% Input is K by n_input_cells
% Target is K by n_output_cells

%[nk1,nj]=size(t);
%[nk2,ni]=size(xcell{1});
nk1=4;nj=4;
nk2=4;ni=900;
%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

% load previous weights
load wq;
load zq;
%ker=[1 1 1];
SE=strel('disk',3);



% set goal
epsilon=10^(-2);

% initialization and declaration
counter=0;
e=ones(nk1,nj);
err=zeros(1,epoch);

% construct training loop
while (sum(sum(e.^2)) > epsilon) & (counter < epoch)
   order=randperm(nk1);
   counter=counter+1;
    
% start an epoch
 for k=1:nk1
     
     % FEATURE EXTRACTION: HIT-MISS TRANSFORM
     for nn=1:4
     
     xe{nn}=double(imerode(xcell{nn},SE));
     xd{nn}=double(imdilate(xcell{nn},SE));
     xcell{nn}=abs(xe{nn}-xd{nn});
     end;
     
     for nu=1:4 
     xcelldummy{nu} = xcell{nu};
     xcelldummy{nu} = imresize(xcelldummy{nu},[30 30]);
     end;
     
     %save ht xhitmiss;
     x1=xcelldummy{1};
     x2=xcelldummy{2};
     x3=xcelldummy{3};
     x4=xcelldummy{4};
    
     v1 = x1(:); 
     v2 = x2(:); 
     v3 = x3(:); 
     v4 = x4(:); 
    
     x = [v1';v2';v3';v4'];
     
     % CLASSIFICATION: TRAINING NETWORK
     np=order(k);
     f=[x(np,:),1];
% feedforward
     h=logsigfx(f*WX); %logistic(squashing) function
     h=[h,1]; % bias=1
     o=logsigfx(h*ZX); %logistic(squashing) function
     e(np,:)=t(np,:)-o; % (t-o)=e
% correction for output layer: Z
     delta=e(np,:).*o.*(1-o);
% correction for hidden layer: W
     deltaH=h.*(1-h).*(delta*ZX');
% weight update
     ZX=ZX+(lrate*(h'*delta));
     WX=WX+(lrate*(f'*deltaH(1:noHidCells)));
  end; 
  % end of an epoch
  
  clc;
  err(counter)=sum(sum(e.^2));
  disp(['epoch=',int2str(counter),...
       ' error=',num2str(err(counter)),...
       ' lrate=',num2str(lrate)]);
end;
save zq ZX;
save wq WX;

% give final answer
  O=logsigfx([x,ones(nk1,1)]*WX);
  O=logsigfx([O,ones(nk1,1)]*ZX);
E=err;