function [C,CM,NRCLRET,CC,CR] = wan_clus_autom(A,MAXNRPOINT,PLOT,PLOTEM,A_ORG)

A=normalise_misval(A',1,1)';

if nargin < 4 | isempty(PLOTEM) | PLOTEM==0
   PLOTEM = 0;
   AR=[];
else
   AR=randdataperm(A');
   AR=normalise(AR,1,1)';
end

if nargin < 3 | isempty(PLOT)
   PLOT = 1;
end

if nargin < 2 | isempty(MAXNRPOINT) | MAXNRPOINT < 1
   MAXNRPOINT=1;
end

samples = size(A,2);
E=size(A,1);
D=E-2;
R=sqrt(E-1);

CLTHR=0.95;
REITERTHR=0.1;

EXTTRESH=R/2;

TFPTH=5;
BPRTH=5;

cluster=0;
CM=[];
CC=[];
CR=[];
TOOFEWPOINTS=0;
C=zeros(1,samples);
NRNOCONV=1;
BPR=0;
RRES=0;

while NRNOCONV<2 & TOOFEWPOINTS < TFPTH & BPR < BPRTH & RRES < 2
   PT=find(C==0);
   SP=A(:,PT);
   if isempty(SP)
      break
   end
   [ME,CE,CONV] = wan_shr_adap(SP,EXTTRESH);
   if CONV==1
      if PLOTEM==1
         [RADNW,CONVEM,PCRET,PBRET,VARRET,EM,FIG] = exp_max(SP,ME,EXTTRESH,CLTHR,PLOTEM,AR(:,PT));
         MESS=[' - See figure ',num2str(FIG)];
      else
         [RADNW,CONVEM] = exp_max(SP,ME,EXTTRESH,CLTHR);
         MESS=[];
      end
      CONV=CONVEM;
   end
   if CONV==1
      NRNOCONV=0;
      if abs(RADNW-EXTTRESH)/EXTTRESH < REITERTHR
          Q=find(dist_misval(SP,ME,2)<RADNW); 
          BPR=0;
          RRES=0;
          if length(Q)>MAXNRPOINT
             cluster=cluster+1;
             CC=[CC ME];
             CR=[CR RADNW];
             C(PT(Q))=cluster;
             disp(['		Found cluster with ',num2str(length(Q)),' genes - Radius is ',num2str(RADNW),MESS])
             TOOFEWPOINTS=0;
             EXTTRESH=RADNW;
          else
             C(PT(Q))=-1;
             disp(['             Discarded cluster with ',num2str(length(Q)),' genes - Radius is ',num2str(RADNW),MESS])
             TOOFEWPOINTS=TOOFEWPOINTS+1;
             EXTTRESH=RADNW;
         end
      else
         EXTTRESH=RADNW;
         BPR=BPR+1;
         if BPR==BPRTH
            disp(['             Radius cannot be predicted: Algorithm stopped',MESS])
         else 
            disp(['             Badly predicted radius / Re-iteration with new radius: ',num2str(RADNW),MESS]) 
         end
      end              
   else
      NRNOCONV=NRNOCONV+1;
      if NRNOCONV==1
         EXTTRESH=R/2;
         RRES=RRES+1;
         if RRES==2
            disp(['No convergence: Algorithm stopped'])
         else
            disp(['             Failure to converge: Radius Reset to (1/2)*R / Re-iteration'])
         end
      else
         disp(['No convergence: Algorithm stopped'])
      end
   end  
   if TOOFEWPOINTS==TFPTH
      disp(['             No more significant clusters found: Algorithm stopped'])
   end
end

NC=6;

if nargin==5
   AT=A_ORG';
else
   AT=A';
end

if PLOT==1
   SP=ceil(cluster/NC);
   f=figure;
   fl=0;
   for loop=1:cluster
   MV=AT(find(C==loop),:);
      if ~isempty(MV)
         MN=mean_misval(MV',2)';
         CM=[CM MN'];
      	 hs = subplot(SP,NC,loop);
      	 set(hs,'FontSize',7);
      	 for iter=1:sum(C==loop)
	        plot(MV(iter,:),'k-')
			hold on
      	 end
         p = plot(MN,'r');
      	 set(p,'LineWidth',2.5);
         t = title(['Cluster ',num2str(loop),'  NOG=',num2str(sum(C==loop))]);
      	 set(t,'FontSize',7);
  		 set(t,'Color',[0 0 1]);
         if NC*ceil(loop/NC)>cluster & fl==0
            xl = xlabel('Experiments');
      		set(xl,'FontSize',7);
      		yl = ylabel('Normalized expr. level');
            set(yl,'FontSize',7);
            fl=1;
   		end
        hold off
     end
  end
end
   
C=C-(C>-1);
NRCLRET=cluster;




