%Usage: [pr,pc,no] = detect(A,rat,iout); % for row index "pr" and col index "pc"
%      A =  Original matrix              % with no = no of swaps in pr pc
%    rat =  ratio of large entry:average entry (say 1.5)
%   iout =  1/0 for plottting or not
% e.g. n=16; G=numgrid('A',n); A=full(delsq(G)); m=floor(size(A,1)/2); close all
%      A(2:3,m-1:m+1)=5;A(4:7,m:m+4)=5; [row col number]=detect(A,1.1,1);
%      subplot(221),spy(A),subplot(222),mesh(A);axis tight
%      subplot(223),spy(A(row,col)),subplot(224),mesh(A(row,col)),axis tight

%      spy(A), figure;spy(A(row,col))

function [pr,pc,no] = detect(A,rat,iout);
    if nargin<1, help detect, return; end
    if nargin<3, iout=0; end
  [m,n]=size(A); [dr,dc,avg_r,avg_c] = gradien(A,iout,rat);
        pr=1:n; pc=1:n; 
   i = find(dr>=avg_r); pr = [setdiff(pr,i) i]; nr=length(i);
if iout==1,subplot(211), plot(i,ones(size(i))*avg_r,'g^'), end
   i = find(dc>=avg_c); pc = [setdiff(pc,i) i]; no=length(i);
if iout==1,subplot(212), plot(i,ones(size(i))*avg_c,'g^'), figure, end
   no=nr+no; % total perms
  fprintf('%s:  permutations identified = %d\n', mfilename,no)

function [dr,dc,avg_r,avg_c] = gradien(A,iout,rat);
[m,n]=size(A);
%%dr=sum(abs([A(2:m,:);A(1,:)] - A)'); % by ROWS -- ONLY good for boundaries!!
%%dr=max(dr,sum(abs([A(m,:);A(1:m-1,:)] - A)') ); % by ROWS
   dr=max(abs(A)'); % by col info
%%dc=sum(abs([A(:,2:n),A(:,1)] - A) ); % by COLS
%%dc=max(dc,sum(abs([A(:,n),A(:,1:n-1)] - A) )); % by COLS
   dc=max(abs(A)); % by row value info
%dr=dr+dc; dc=dr; % Joint for symmetric perm !!
  avg_r=mean(dr)*rat; avg_c=mean(dc)*rat;
if iout==1
subplot(211), plot(dr,'-b+'),hold on ,plot(1:m,avg_r*ones(m,1),'r-'),axis tight
  title('detection by rows')
subplot(212), plot(dc,'-b+'),hold on ,plot(1:n,avg_c*ones(n,1),'r-'),axis tight
  title('detection by columns')
end
return
