coarserefine2d(f,intrv,theta,n,topc)

function [x,param,values,lambda,flagbdry] = coarserefine2d(f,intrv,theta,n,topc)
% 2-D adaptive residual subsampling method for radial basis function 
% interpolation
%
% f: function defined on interval [a b] x [a b].
% intrv: interval [a b] where f is defined.
% theta: threshold interval [thetar thetac] where thetac < thetar.
%        thetar: refinement threshold.
%        thetac: coarsening threshold.
% n: will be used to generate NxN midpoint centers.
% alpha: global parameters for multiquadric rbfs
%
% Example 1 :  f = @(x) exp(-sum(x.^2,2))
%              coarserefine(f,[-1 1],[1e-3 1e-6],8,10)
%
% Example 2 :  coarserefine(@(x)myfun(x,5),[-1 1],[1e-3 1e-6],8,10)
%              %----------------------%
%              function y = myfun2(x,c)
%              y = exp(-c^2*sum(x.^2,2));
%              %----------------------% 
%
% Note that x is a 2-column vector containing absis and coordinate of
% centers on its first and second column respectively. 
%  
% For reference, see:
% Adaptive residual subsampling methods for radial basis function
% interpolation and collocation problems. submitted to Computers Math. Appl
% Tobin A. Driscoll and Alfa R.H. Heryudono    05/14/2006
%  
% MATLAB 7 is recommended.

% Adaptive residual subsampling parameters
thetar = theta(1); thetac = theta(2);
a = intrv(1); b=intrv(2);

% Corner centers (never coarsened)
xc = [b b;a b;a a;b a];
cc = topc*ones(4,1); % corner shape parameters
u0c = feval(f,xc);

% Perform initial coarse discretization
ibox = initboxes2d(n,[a b]);

N = length(ibox);
% Create first generation of "residual boxes." of interior points
[ibox,newres] = dorefine2d(ibox,true(N,1),1);

center = false(size(ibox));  center(1:N) = true;
residual = false(size(ibox));  residual(newres) = true;

% Begin main adaptive interpolation iteration
flag = true;
while any(flag)
  N = sum(center);  M = sum(residual);
  centerbox = ibox(center); residbox = ibox(residual);

  x = cat(1,centerbox.midpt);
  c = topc*2.^(cat(1,centerbox.depth));
  w = cat(1,residbox.midpt);

  x = [x;xc];
  c = [c;cc];

  [A,B] = deal(zeros(N+4),zeros(M,N+4));
  for j=1:N+4
      A(:,j) = mq2d(x,x(j,:),c(j));
      B(:,j) = mq2d(w,x(j,:),c(j));
  end

  lambda = A\feval(f,x);
  resid = B*lambda - feval(f,w);

  %boxarea = cat(1,residbox.dimen);
  %zerolength = boxarea == 0;
  %boxarea(zerolength) = 1;
  %boxarea = prod(boxarea,2);

  normerr = -ones(size(ibox));
  % Error (residual) indicator
  normerr(residual) = abs(resid); % Pointwise error Indicator
  %normerr(residual) = abs(resid).*boxarea; % Area-based error Indicator

  % Look for coarsening opportunities in the region.
  count = 0;
  for i = find(center(:)).'
    check = abs(normerr(ibox(i).child));
    if sum(check < thetac)==ibox(i).Nchild
      center(i) = false;
      residual(i) = true;
      residual(ibox(i).child) = false;
      count = count+1;
    end
  end
  fprintf('Removing %i centers and ',count)

  % Boxes needing refinement in the region.
  flag = normerr > thetar;
  newcen = find(flag);
  center(newcen) = true;
  residual(newcen) = false;
  [ibox,newbox] = dorefine2d(ibox,flag,0);
  center(newbox) = false;
  residual(newbox) = true;

  fprintf('Adding %i centers.\n',sum(flag))
  %fprintf('------Press Enter------\n')
end

x = cat(1,ibox(center).midpt); 
flagbdry = cat(1,ibox(center).isboundary)==true;
c = topc*2.^cat(1,ibox(center).depth);
x = [xc;x];
c = [cc;c];
flagbdry = [false(4,1);flagbdry];
N = length(x);

if nargout > 1
    param = c;
    if nargout > 2
        values = feval(f,x);
        if nargout > 3
          A = zeros(N);
          for j=1:N
          A(:,j) = mq2d(x,x(j,:),c(j));
          end
          lambda = A\feval(f,x); 
        end
    end
end