%========================================================================
%FMM demo program to COMPUTE the near neighbours' list and
%to optinally illustrate interaction lists/parents' neighbours
%============================================ index_demo.m > dfmm3 > dfmm 
%Usage: out=index2(node, level, io, opts)  %e.g index2(3,2,1,'white')
%                                          % io = 0 (No), 1 (Plot), 2 (IL)
%   or  out=index2(node, level);           %e.g index2(3,2) % node=3
%   or  out=index2(node);                  %e.g index2(3)   % for level=2
%
% Stuart C. Hawkins and Ke Chen (University of Liverpool, UK)
%
% SEE also intera.m
 
% Main internal notation ___________________________________________________
%  F0        : List of starting global index addresses of each level
%  ANCESTORS : List of positions of the node's ancestors at their level
% index ---- : Convert a local index (i,j) at level j to global index k
% vert ----- : Find the local index j for (i,j) from a global k
% horiz ---- : Find the local index i for (i,j) from a global k
% neighbours : Find the neighbour's index list for any node on a level
%            : (using near_neighbours)
% fmm_draw   : Plot near neighbours' info on all levels for node/ANCESTORS
%            : (using draw_near_only)
% parent     : Find the parent's index id, given a child's index (i,j)
% children   : Find a parent's 4 children, given a children index list k
% Figure     : Automatic subplotting
% --------------------------------------------------------------------------

function out=fmm_demo(node,levs,iout,opts)
 if nargin<4, opts='color'; end
 if nargin<3, iout=0;       end  % plotting or not
 if nargin<2, levs=2;       end
 if nargin<1, help index2, return, end
 fprintf('-----index2.m------Finding near neighours for node=%d, ',node)
 fprintf('\n\t\ton a given finest level = %d\n',levs) 
  if node<1|node>4^levs, help index2, out=-1;
    fprintf('\tInput "node" out of range  /  1:%d\n',4^levs); return, end 
% get starting index of level l=levs' starting 
 out.node=node;
 for ll=1:levs
 F0(ll)=index(1,1,ll); % position = 4^0+4^1+...+4^(l-1)
 end

% update target --- node in "global" index (shift = f0-1)
 node=node-1+F0(levs);

% get neighbours of node
  near=neighbours(node,levs);
  out.near=near+1-F0(levs);out.lev=levs;
  disp(' Main task of index2.m complete ...')
  if iout<1, return, end % not to draw pictures
% ********** Main Mfile Ended Here **********************

if nargin>2 & length(findstr('white',opts))>0 
     % set colour data for b/w output
      white=[1.0,1.0,1.0];  % draw with colours specified
      fmm_draw(F0,iout,near,node,levs,white,white,white) 
else % draw with default colours
      fmm_draw(F0,iout,near,node,levs)    
end
return %*************************************************

%========================================================
function [near]=neighbours(f,lev) 
% get coordinates of f (f shifted by f0-1)
  j=vert(f,lev); i=horiz(f,j,lev); 
% get (i,j) coords of neighbours of f
  [ii,jj,ll]=near_neighbours(i,j,lev); 
% get indexes of the neighbours of f from their (i,j) coords
  near=index(ii,jj,ll);

%========================================================
function [ii,jj,ll]=near_neighbours(i,j,l) 
% initialise arrays -- neighbour (i,j) list + lev list
  ii=[]; jj=[]; ll=[]; 
% loop through coordinates----Search in +/-1 Box
  for cj=max(1,j-1):min(2^l,j+1) 
      for ci=max(1,i-1):min(2^l,i+1) 
          if( cj~=j | ci~=i )
              ll=[ll;l]; ii=[ii;ci]; jj=[jj;cj];
          end 
      end 
  end

%===================================================
function l=level(f)  %linking f=partial sums to lev
  l=floor( log(3*(f-1)+1)/log(4) ); 
% Again the storage system = 4^0+4^1+4^2+...+4^LEV 

%===================================================
function j=vert(f,l)
  j=floor( (f-1-(4.^l-1)/3)./(2.^l)+1 );

%===================================================
function i=horiz(f,j,l)
  i=f-2.^l.*(j-1)-(4.^l-1)/3;

%===================================================
function f=index(i,j,l); % (i,j) to single f (shifted)
  f=i+2.^l.*(j-1)+(4.^l-1)/3;

%=======================================================
function fmm_draw(F0,iout,near,target,l,ctarget,cnear,cfar)
%---------------------% target = shifted node on level=l
if nargin<6 % set default for ctarget  (plot colors)
    ctarget='b'; %%[1.0,0.0,0.0]; %1 0 0
end 
if nargin<7 % set default for cnear
    cnear=[1.0,0.5,0.0];
end 
if nargin<8 % set default for cfar
    cfar=[0.5,1.0,0.5]; %1 1 0
end, hold off

%---------------------%target = shifted node on level=l
%---------------------%Loop through coarse/fine levels 
for ll=l:-1:1
    Figure(ll)
    %-----------------% draw far field-----------------
    ff=1;            % initialise counter 
    for j=1:2^ll      % loop through rows 
        for i=1:2^ll  % loop through entries in row
            patch([i-1 i-1 i i],[j j-1 j-1 j],cfar); 
            t=text(i-0.7,j-0.7,num2str(ff));
        set(t,'fontsize',30/(ll+2)), set(t,'color','b')
            ff=ff+1;  % increment counter
        end 
    end 
    %-----------------% draw near field----------------------
    if ll==l % finest level from input list
       t=title(['Finest and Near-field Level ' num2str(ll)]); 
            set(t,'fontsize',15); set(t,'color','r')
       [i j]=draw_near_only(near,target,l,F0(l),...
                           cnear,ctarget);%Finest i j
    %-----------------% Find ancestors' id-------------------
           ANCESTORS=[]; % global index (adding F0)
         for lev=l-1:-1:1
           [ip jp lp]=parent(i,j,lev+1);
           ANCESTORS=[index(ip,jp,lp) ANCESTORS];
            i=ip; j=jp;
         end % lev
    else %------------% Plot near neighbours of parents------
       t=title(['Parent and Far-field Level ' num2str(ll)]); 
            set(t,'fontsize',15); set(t,'color','b')
        Target=ANCESTORS(ll);
        t=sprintf('near%d=neighbours(Target,ll);',ll);
        eval(t) %-----generate near1, near2 etc
        f0=F0(ll);
t=sprintf('draw_near_only(near%d,Target,ll,f0,cnear,ctarget);',ll);
        eval(t) %-----Plot each near1, near2 etc
    end % if
end % lev    % EXTRA below
    t=sprintf('near%d=near;',l); eval(t) %-----Make nearLEV 

if iout==2 %=========== Request to plot interaction lists=========
for ll=l-1:-1:1        %---------------------------------%
    Figure(ll+1), hold on
    t=sprintf('in_list=children(near%d,ll,F0(ll));',ll);
        eval(t) %-----Find the interaction list (Step 1) -----

    t=sprintf('IL=setdiff(in_list,near%d-F0(ll+1)+1);',ll+1);
        eval(t) %-----Find the interaction list (Step 2) -----
    jfine= vert(IL+F0(ll+1)-1,ll+1);
    ifine=horiz(IL+F0(ll+1)-1,jfine,ll+1);
    for k=1:length(ifine)
    t=plot(ifine(k)-0.5,jfine(k)-0.5,'rx'); set(t,'markersize',20/ll)
    end %k
end % for
end % iout

%===================================================
function [ii,jj,ll]=parent(i,j,l)
ll=l-1; ii=ceil(i/2); jj=ceil(j/2);

%===================================================
function listo=children(k,lev,f0)
   k = k + 1 - f0;  % shifted back to [1, 4^lev]
  two=2^lev; twop=two*2;
   i = rem(k-1,two)+1; % horizontal index
  k1 = 4*k - 2*i; % vector
  listo = [ k1-1; k1+twop-1; k1+twop; k1]; 

%=====================================================
function [i,j]=draw_near_only(near,target,l,f0,c1,c2);
        for k=1:length(near) % get coordinates of k
            l=level(near(k));
            j=vert(near(k),l);
            i=horiz(near(k),j,l); % draw patch
            patch([i-1 i-1 i i],[j j-1 j-1 j],c1); 
            t=text(i-0.7,j-0.7,num2str(near(k)-f0+1));
        set(t,'fontsize',30/(l+2)), set(t,'color','k')
        end 
        %-------------% draw target------------------
          j=vert(target,l); i=horiz(target,j,l); 
        patch([i-1 i-1 i i],[j j-1 j-1 j],c2); 
        t=text(i-0.7,j-0.7,num2str(target-f0+1));
        set(t,'fontsize',30/(l+2)), set(t,'color','k')
        t=plot(i-0.5,j-0.6,'r+'); set(t,'markersize',25/l)
        axis off
%=====================================================
function Figure(in)
   m=4; n=1;  % for m x n subplotting
   total = m*n;
   fig = floor( (in-1)/total );     % sub = subplot no
   sub = in - total*fig; fig=fig+1; % fig = figure no
   figure(fig); 
   subplot(m,n,sub); hold on
