% Check pointwise evaluation error of native & close-eval quadrature scheme in
% BIE Laplace & Helmholtz
% Barnett 3/15/12

clear; expt='f';   % choose expt and N = base # bdry pts, eg 70
%N = 100; s = segment.smoothstar(N, 0.3, 3, 0.2); % ......... choose shape........
%N = 100; s = segment.smoothnonsym(N, 0.3, 0.3, 3);  % standard
%N=3000; o.napproxv = 1e3; s = segment.smoothnonsym(N, 0.3, 0.3, 20, o);  % wiggly

N=1000; mmax = 22; wsc = .1; % rand wiggly, 1e-10 acc
%N=2000; mmax = 20; wsc = .15; % rand wiggly, typ fuzz sc (20, .15)
%N=2000; mmax = 30; wsc = .11; % rand wiggly, 1e-10
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);

%a=.3;b=.3;w=3; a1=0.1; a2=0.1; t0=3*pi/4; % smoothnonsym w/ gaussian bump @ t0
%s = segment.radialfunc(N, {@(q) 1 + a*cos(w*(q+b*cos(q))) - a1*exp(-0.5*((q-t0)/a2).^2), @(q) -a*sin(w*(q+b*cos(q))).*w.*(1-b*sin(q)) + (a1/a2^2)*(q-t0).*exp(-0.5*((q-t0)/a2).^2) }); %don't need curvatue yet since Alpert acc ok given bump
% .......................................................................

%s.qtype = 'g'; s.requadrature(N);  % debugging (should be periodic)
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);
% getting slow! 3 secs for dx=.005 (2.5e5 pts), napproxv=3e3 (2e8 pts*verts/sec)

% 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));
figure; subplot(1,2,1); imagesc(xx(1,:),yy(:,1),log10(abs(u-fz))); colorbar;
caxis([-16 0]); axis equal tight; set(gca,'ydir','normal'); title('native');
drawnow;

% z0= 0; % do all through Jfilter method of layerpot. choose origin z0
%o = []; o.Jfilter.M = 50; o.Jfilter.rescale_rad = 1.0; o.Jfilter.origin = z0;
%Ab = lpf.eval(pointset(zz(ii)), o);

subplot(1,2,2); view(2); colorbar; caxis([-16 0]); axis equal tight; set(gca,'ydir','normal'); title('surrogate');

% now do surrogate box evals (plot the boxes must be as you go)
if 0   % naive O(N^2) LE coeff eval, replacing boxes values...
%B = 25; h = pi/B; % surrogate scheme based on B square s-plane-boxes...
B = 400; h = pi/B; % surrogate scheme based on B square s-plane-boxes...
Nf = 4*N; cof = quadr.interptrig(co, Nf); % fine grid
sf = utils.copy(s); sf.requadrature(Nf);  % fine segment
lpf = layerpot(sf, 'D'); lpf.doms = d;    % fine layerpot
alpha0 = h; M = 10;            % M=p, max degree of local exp
tic; nv = 5; l = h*(1:nv)/nv; % pts to build little polygon, real from [0 h]
%profile clear; profile on;
for b=1:B, t = 2*h*b; z0 = s.Z((t+1i*alpha0)/(2*pi)); % loop over boxes
  %v = s.Z([t-h+l, t+l, t+h+1i*l, t+1i*h+h-l, t+1i*h-l, t-h+1i*h-1i*l].'/(2*pi));  % half-height 2x1 boxes
  % square boxes:
  v = s.Z([t-h+l, t+l, t+h+2i*l, t+2i*h+h-l, t+2i*h-l, t-h+2i*h-2i*l].'/(2*pi));
  ib = utils.inpolywrapper(zz, v);
  plot3(real(v), imag(v), 1+0*v, 'k-'); hold on;
  plot3(real(z0),imag(z0),1,'k.');
  LEo.besselcode = 'm'; LEo.real = 0; LEo.rescale_rad = s.dist(z0);
  LE = regfbbasis(z0, M, LEo); LE.doms = d;
  S2L = layerpot.localfromDLP(sf, LE); coLE = S2L * cof;
  %Ab = LE.eval(pointset(zz(ib))) * S2L; % matrix-matrix mult.
  %coLE = S2L*cof % debug in expt='1' k=0 case
  if ~isempty(find(ib)), u(ib) = LE.eval(pointset(zz(ib))) * coLE; end
  %fprintf('b=%d, max err in box = %.3g\n', b, max(abs(u(ib)-fz(ib))))
  if mod(b,10)==0, hold on; imagesc(xx(1,:),yy(:,1),log10(abs(u-fz))); view(2); drawnow; end
end
fprintf('time for strip eval = %.3g s\n',toc); %profile off; profile viewer

else   % 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; s.plot; hold on; plot([v; v(1)], 'go-'); plot(z0, '*'); % 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
    % select local indices via being closest to this center:
    %  z0R = z0(find(abs(z0-z0(b))<Rmax)); % nearby centers
    %  iz0 = find(min(abs(repmat(zzc(iR),[1 numel(z0R)])-repmat(z0R,[numel(iR) 1])),[],2) >= abs(zzc(iR)-z0(b)));  % indices in strip which have this z0 closest
    % OR, use polygon-approximated box as above: (1.05 is spill factor)
    v = 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
    %bnei = 8; % how far away in box sizes to include in local neighborhood?
    %nf = mod((floor((b-bnei)/B*Nf):ceil((b+bnei)/B*Nf))-1,Nf)+1; % fine indices
    % found the above not good for curves which approach themselves; use dist.
    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);
end
fprintf('max abs err = %.3g\n',max(abs(u(:)-fz(:))))

% usual color image:  %subplot(1,2,2);
view(2); imagesc(xx(1,:),yy(:,1),log10(abs(u-fz))); axis equal tight; set(gca,'ydir','normal'); title('surrogate'); colorbar; caxis([-16 0]); %hold on; plot(s.x, 'w.'); plot(z0, 'k*'); 
%figure; hist(nnf); title('distribution in number of local correction pts');


if 0 % contour comparison:
figure; errs = 10.^(-16:2:0); %10.^(-16:2:0);
contourf(xx,yy,log10(abs(u-fz)), log10(errs)); caxis([-16 0]);
colorbar; axis equal; hold on; plot([s.x; s.x(1)], 'k.-');
taus = log(maxtau./errs)/N; %[0:.05:.3]; % use complex map imag shifts
o = []; o.levels = taus/(2*pi);
if exist('t0') % hack to get bump to plot imag dist ok
  n=10; no=10; o.to = [((1:n)-.5)/n, (t0+2*a2*(-no:2:no)/no)/(2*pi)] + 0.05i;
end
[c h mis gx gy] = d.showimagparam(o);
title(sprintf('log10 eval errors of native quadr, and predicted, N=%d',N));
xlabel x; ylabel y; axis tight; %axis(1.2*[-1 1 -1 1]);
%print -depsc2 native_eval_tau1_nonsym_N60.eps
end
