% Demo the O(N) fine-FMM strip correction scheme for CE paper. Barnett 4/27/12
% Taken from surrogate_eval_err.m

clear; expt = 'f'; % solve
N=1100; mmax = 22; wsc = .1; % rand wiggly
o.napproxv = 3e3; rand('twister',123);  % seed
aj = wsc*(rand(1,mmax)-0.5); bj = wsc*(rand(1,mmax)-0.5); % `rough' shape
s = segment.smoothfourier(N,aj,bj,o);
d = domain(s, 1); k = 20.0; d.k = k;  % choose Laplace or Helmholtz
o = []; d.addlayerpot(s, 'D', o); % d.plot;
dx = 0.005; bb = d.boundingbox;[xx yy] = meshgrid(bb(1):dx:bb(2),bb(3):dx:bb(4));
zz = xx(:)+1i*yy(:); tic; di = d.inside(zz); ii = find(di);
fprintf('inside %d pts, %d verts = %.3g s\n',numel(zz),numel(s.approxv),toc);

% fill tau density ...
if expt=='1'
  if k>0, error('Helmholtz has no tau=1 known potential!'); end
  f = @(z) -1 + 0*z;               % known solution -1 inside
  maxtau = 1;
  co = ones(N,1);                    % unit density
elseif expt=='f'
  accsolv = 0;
  if accsolv, s.requadrature(3*N); end % solve on finer grid
  p = bvp(d);  % now solve a BVP to get a more interesting density
  d.bas{1}.quad='m'; % Kress for self-eval (if kappa known & N<1400)
  %d.bas{1}.quad='a'; d.bas{1}.ord = 16; % Alpert for self-eval (if no kappa)
  if k==0, f = @(z) real(z).*imag(z); % note: better since not zero along the ray 3*pi/4 that bump is on!
    % was: real(z).^2 - imag(z).^2;  % u(x,y) = x^2 - y^2 exact, Lap
  else, f = @(z) exp(1i*k*real(exp(1i*pi/3)'*z)); % u = plane wave
    %f = @(z) 3*besselh(0,k*abs(z-(1-1i)));  % pt src (1,1) or (.9,.8) close!
  end
  maxtau = max(abs(f(s.x)));  % obviously this *isn't* tau max in strip!
  s.setbc(-1, 'D', [], @(t) f(s.Z(t)));
  tic; p.solvecoeffs; %opts.comparefunc = f; p.showsolution(opts);
  fprintf('time fill+solve size %d lin sys = %.3g s\n',numel(p.co),toc)
  if accsolv, co = p.co(2:3:end-1); % get coincident nodes (3N->N downsample)
    s.requadrature(N); else, co = p.co; end % for plotting later
end
% native ....
fz = nan*xx; fz(ii) = f(zz(ii)); fz = reshape(fz, size(xx)); % fill fz, u...
u = nan*xx;
if N<500   % switch for dense (crude eval matrix fill) or FMM
  tic; A = d.bas{1}.eval(pointset(zz(ii)));   % do grid eval by hand
  fprintf('time for filling interior di eval matrix A = %.3g s\n',toc);
  u(ii) = A*co; 
else
  tic; u(ii) = d.bas{1}.evalFMM(co,pointset(zz(ii)));
    fprintf('time for FMM interior di eval  = %.3g s\n',toc);
end
u = reshape(u, size(xx)); unat = u;

figure; imagesc(xx(1,:),yy(:,1),log10(abs(unat-fz))); colorbar;   % fig a
v=[.52 .9 .55 1.02]; hold on;rectangle('position',[v([1 3]) v([2 4])-v([1 3])]);
caxis([-14 0]); axis equal tight; set(gca,'ydir','normal');
colormap(1-gray(256)); text(-1.2,1.25,'(a) $$\log_{10} |u-f|$$','interpreter','latex');
set(gcf,'paperposition', [0 0 4 3]);
print -depsc2 orderN1e3_native.eps

 % improved *local only* correction scheme, O(N).
  B = ceil(N/4); h = 1/(2*B); %surrogate scheme based on B rect s-plane-boxes...
  M = 18;   % M=p max degree of LE
  als = 5.0/N; % imag strip width wrt [0,1]-scaled seg, for acc to 1e-12
  v = s.Z(s.t + 1i*als); % imag-shifted N bdry pts. Now clip off loopy ones:
  vsd = abs(repmat(s.x,[1 N]) - repmat(v.',[N 1])); % find dists v to pO, O(N^2)
  [vsdm, ivsd] = min(vsd,[],1); clear vsd; % get indices of nearest pts
  oi.toz = s.t(ivsd);  % nearest on-pO initial guesses to help invertZ
  vs = s.invertZparam(v,oi); % slow since wobbly curve
  v = v(find(abs(imag(vs))>0.99*als)); % snip off too-close self-int loops
  fc = ~utils.inpolywrapper(zz(ii), v); % for ii pts, true if in strip
  ic = ii(find(fc)); zzc = zz(ic);      % pts to correct
  t0 = (1:B)/B+1i*als/2; z0 = s.Z(t0); nnf=0*z0;  % all preimages and centers
  Rmax = 2.0*max(s.speed)*als;  % max dist to search for pts near each center
  %figure; plot([s.x; s.x(1)], '-'); hold on; plot([v; v(1)], 'go-'); plot(z0, '*'); axis equal; % show geom
  Nf = 5*N; cof = quadr.interptrig(co, Nf); % fine grid: choose factor (5)
  sf = utils.copy(s); sf.requadrature(Nf);  % fine segment
  lpf = layerpot(sf, 'D'); lpf.doms = d;    % fine layerpot
  tic; uc = lpf.evalFMM(cof,pointset(zzc));    % fine FMM eval in strip
  fprintf('fine FMM in strip = %.3g s\n',toc);
  % setup a dummy local lp and seg, whose x,nx,w will be overwritten for each b:
  tic;% profile clear; profile on;
  sl = utils.copy(sf); lpl = layerpot(sl, 'D'); lpl.doms = d; % dummy lp
  nv = 10; l = (1:nv)/nv;  % polygonal approx to box
  box = [l-.1i 1-.1i+.6i*l 1+.5i-l .5i-.6i*l].'*2*h;
  for b=1:B % loop over boxes, get indices in `box', correct them...
    iR = find(abs(zzc-z0(b))<Rmax); % strip pt indices
    vb = s.Z(b/B - h + box); iz0 = utils.inpolywrapper(zzc(iR), vb);
    zzR = zzc(iR(iz0)); %hold on; plot(zzR, '.'); drawnow % check the set of pts
    bs = max(abs(s.Zp(t0(mod(b+[-3:3]-1,B)+1))))*2*h; % this box size 
    nf = find(abs(sf.x - z0(b))<bs*2.5); % distance criterion for local fine pts
    nnf(b) = numel(nf);
    sl.x = sf.x(nf); sl.nx = sf.nx(nf); sl.w = sf.w(nf); sl.speed=sf.speed(nf);
    coR = cof(nf);   % local fine coeffs
    ucorr = -lpl.eval(pointset(zzR)) * coR;    % subtract off direct eval
    LEo.besselcode='m'; LEo.real=0; LEo.rescale_rad=max(abs(zzR-z0(b)));%0.7*bs;
    LE = regfbbasis(z0(b), M, LEo); LE.doms = d;  % local exp
    S2L = layerpot.localfromDLP(sl, LE); coLE = S2L * coR; % LE coeffs
    ucorr = ucorr + LE.eval(pointset(zzR)) * coLE;   % add local-exp potential
    uc(iR(iz0)) = uc(iR(iz0)) + ucorr; % correct the fine FMM pot in local box
  end
  u(ic) = uc; % overwrite corrected strip back into main u array
  %profile off; profile viewer
  fprintf('local-exp correction on %d pts in strip = %.3g s\n',numel(uc),toc);
fprintf('max abs err = %.3g\n',max(abs(u(di)-fz(di))))

figure; imagesc(xx(1,:),yy(:,1),log10(abs(u-fz))); colorbar;   % fig c
caxis([-14 -11.2]);
axis equal tight; set(gca,'ydir','normal');
colormap(1-gray(256)); text(-1.2,1.25,'(c) $$\log_{10} |\tilde{u}-f|$$','interpreter','latex');
set(gcf,'paperposition', [0 0 4 3]);
print -depsc2 orderN1e3_localexp.eps

figure; box =  [l 1+.5i*l 1+.5i-l .5i-.5i*l].'*2*h;  % shape for plotting. fig b
b=34; % in the zoom region
iR = find(abs(zzc-z0(b))<Rmax); % strip pt indices
vb = s.Z(b/B - h + box); iz0 = utils.inpolywrapper(zzc(iR), v);
zzR = zzc(iR(iz0)); %hold on; plot(zzR, '.'); drawnow % check the set of pts
bs = max(abs(s.Zp(t0(mod(b+[-3:3]-1,B)+1))))*2*h; % this box size 
nf = find(abs(sf.x - z0(b))<bs*2.5); % distance criterion for local fine pts
fill(real(vb), imag(vb), [.8 .8 .8]); hold on;
plot(z0, 'm*'); circle(real(z0(b)), imag(z0(b)), bs*2.5, 'g-'); hold on;
plot(z0(b), 'k*'); % 'markersize', 10);
for b=20:50, vb = s.Z(b/B - h + box); plot([vb; vb(1)], 'b-'); end % show boxes
plot([s.x; s.x(1)], 'k.', 'markersize', 20, 'color', [.6 .6 .6]);
hold on; plot([sf.x; sf.x(1)], 'k.-', 'markersize', 1);
plot(s.Z(s.t + 1i*als), 'r-.');  % correction strip
plot(v, 'r-');  % correction strip after snips
plot(sf.x(nf), '.', 'markersize',7,'color',[0 0 0]);
axis equal; axis([.52 .9 .55 1.02]);
text(.55,1.04,'(b) geometry');
set(gcf,'paperposition', [0 0 3 4]);
print -depsc2 orderN1e3_geomzoom.eps
