% June 2014, P.G. Martinsson
%
% This program compares different techniques for evaluating the rank of
% interaction between "source" and "target" regions of different
% geometries.
%
% This has not been proofed at all - use at your own risk!

function main_laplace_expansions

LOCAL_box_to_circle
%LOCAL_box_to_line
%LOCAL_box_to_box
%LOCAL_disc_to_circle_quad
%LOCAL_box_to_circle_quad
%LOCAL_box_to_line_quad
%LOCAL_curve_to_line_quad
%LOCAL_box_to_line_randomized

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_box_to_circle

rng(0)

NS     = 300;       % Number of source points.
NT     = 200;       % Number of target points.
cc     = [0;0];     % Center of expansion.
a_in   = 0.5;       % Source box is of [-a_in,a_in] x [-a_in,a_in]
k_max  = 80;        % Max rank to be investigated.

%%% Set up the source geometry
xxs          = a_in*(1 - 2*rand(2,NS));
box_geom_src = [2*a_in,cc(1),cc(2)];

%%% Set up the target geometry
tt  = linspace(0,2*pi*(1 - 1/NT),NT);
xxt = 1.5*[cos(tt);sin(tt)];

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
A       = 0.5*log(DD1.*DD1 + DD2.*DD2);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the ID
[Q,R,ind] = LOCAL_pivoted_QR3(A,k_max);
err_ID    = [ss(1),NaN*ones(1,k_max)];
for i = 2:length(kk)
  k         = kk(i);
  err_ID(i) = norm(A(:,ind) - Q(:,1:k)*R(1:k,:));
end

%%% Perform the multipole expansions.
err_mpole = [ss(1),NaN*ones(1,k_max)];
P_max     = floor(k_max/2)-1;
CC        = LOCAL_T_ofs(box_geom_src,P_max,xxs);
BB        = LOCAL_T_tfo(xxt,box_geom_src,P_max);
for p = 0:P_max
  err_mpole([2*p+2,2*p+3]) = norm(A - real(BB(:,1:(p+1))*CC(1:(p+1),:)));
end

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
han = plot(a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'MarkerSize',20)
set(han(3),'MarkerSize',20)
axis equal
axis off
title(sprintf('The geometry'))

%%% Plot the errors.
figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'x-',...
         kk,(1/ss(1))*err_ID,'x-',...
         kk,(1/ss(1))*err_mpole,'x-');
legend('svd','ID','mpole')
title(sprintf('Expansion errors'))
xlabel('k')
ylabel('||A - A_k||/||A||')

%%% Plot a skeleton
figure(3); subplot(1,1,1); hold off;
acc  = 1e-10;
k_id = max(find(err_ID > (ss(1)*acc)));
han = plot(xxs(1,:),xxs(2,:),'k.',...
           a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,ind(1:k_id)),xxs(2,ind(1:k_id)),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'MarkerSize',20)
set(han(1),'MarkerEdgeColor',0.7*[1,1,1])
set(han(2),'LineWidth',2)
set(han(3),'MarkerSize',20)
set(han(4),'MarkerSize',20)
axis equal
axis off
title(sprintf('Skeleton points: eps=%0.1e  k=%d',acc,k_id))
keyboard
    
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_box_to_line

rng(0)

NS     = 300;       % Number of source points.
NT     = 200;       % Number of target points.
cc     = [0;0];     % Center of expansion.
a_in   = 0.5;       % Source box is of [-a_in,a_in] x [-a_in,a_in]
k_max  = 80;        % Max rank to be investigated.

%%% Set up the source geometry
xxs = a_in*(1 - 2*rand(2,NS));
box_geom_src = [2*a_in,cc(1),cc(2)];

%%% Set up the target geometry
nside = ceil(NT/4);
h     = 3/nside;
tt    = linspace(-1.5,1.5-h,nside);
xxt   = [tt,1.5*ones(1,nside),-tt,-1.5*ones(1,nside);...
         -1.5*ones(1,nside),tt,1.5*ones(1,nside),-tt];

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
A       = 0.5*log(DD1.*DD1 + DD2.*DD2);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the ID
[Q,R,ind] = LOCAL_pivoted_QR3(A,k_max);
err_ID    = [ss(1),NaN*ones(1,k_max)];
for i = 2:length(kk)
  k         = kk(i);
  err_ID(i) = norm(A(:,ind) - Q(:,1:k)*R(1:k,:));
end

%%% Perform the multipole expansions.
err_mpole = [ss(1),NaN*ones(1,k_max)];
P_max     = floor(k_max/2)-1;
CC        = LOCAL_T_ofs(box_geom_src,P_max,xxs);
BB        = LOCAL_T_tfo(xxt,box_geom_src,P_max);
for p = 0:P_max
  err_mpole([2*p+2,2*p+3]) = norm(A - real(BB(:,1:(p+1))*CC(1:(p+1),:)));
end

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
han = plot(a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'MarkerSize',20)
set(han(3),'MarkerSize',20)
axis equal
axis off
title(sprintf('The geometry.'))

%%% Plot errors.
figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'x-',...
         kk,(1/ss(1))*err_ID,'x-',...
         kk,(1/ss(1))*err_mpole,'x-');
legend('svd','ID','mpole')
title(sprintf('Expansion errors.'))
xlabel('k')
ylabel('||A - A_k||/||A||')

%%% Plot a skeleton
figure(3); subplot(1,1,1); hold off;
acc  = 1e-12;
k_id = max(find(err_ID > (ss(1)*acc)));
han = plot(xxs(1,:),xxs(2,:),'k.',...
           a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,ind(1:k_id)),xxs(2,ind(1:k_id)),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'MarkerSize',20)
set(han(1),'MarkerEdgeColor',0.7*[1,1,1])
set(han(2),'LineWidth',2)
set(han(3),'MarkerSize',20)
set(han(4),'MarkerSize',20)
axis equal
axis off
title(sprintf('Skeleton points: eps=%0.1e  k=%d',acc,k_id))
keyboard
    
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_box_to_box

rng(0)

NS     = 300;       % Number of source points.
NT     = 300;       % Number of target points.
cc     = [0;0];     % Center of expansion.
a_in   = 0.5;       % Source box is of [-a_in,a_in] x [-a_in,a_in]
k_max  = 80;        % Max rank to be investigated.

%%% Set up the source geometry
xxs = a_in*(1 - 2*rand(2,NS));
box_geom_src = [2*a_in,cc(1),cc(2)];

%%% Set up the target geometry
xxt = [2.5 - rand(1,NT);0.5 - rand(1,NT)];

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
A       = 0.5*log(DD1.*DD1 + DD2.*DD2);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the ID
[Q,R,ind] = LOCAL_pivoted_QR3(A,k_max);
err_ID    = [ss(1),NaN*ones(1,k_max)];
for i = 2:length(kk)
  k         = kk(i);
  err_ID(i) = norm(A(:,ind) - Q(:,1:k)*R(1:k,:));
end

%%% Perform the multipole expansions.
err_mpole = [ss(1),NaN*ones(1,k_max)];
P_max     = floor(k_max/2)-1;
CC        = LOCAL_T_ofs(box_geom_src,P_max,xxs);
BB        = LOCAL_T_tfo(xxt,box_geom_src,P_max);
for p = 0:P_max
  err_mpole([2*p+2,2*p+3]) = norm(A - real(BB(:,1:(p+1))*CC(1:(p+1),:)));
end

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
han = plot(0.5*[-1,1,1,-1,-1],0.5*[-1,-1,1,1,-1],'k',...
           2 + 0.5*[-1,1,1,-1,-1],0.5*[-1,-1,1,1,-1],'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'LineWidth',2)
set(han(3),'MarkerSize',20)
set(han(4),'MarkerSize',20)
axis equal
axis off
title(sprintf('The geometry.'))

%%% Plot the errors
figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'x-',...
         kk,(1/ss(1))*err_ID,'x-',...
         kk,(1/ss(1))*err_mpole,'x-');
legend('svd','ID','mpole')
%title(sprintf('Expansion errors.'))
xlabel('k')
ylabel('||A - A_k||/||A||')

%%% Plot a skeleton
figure(3); subplot(1,1,1); hold off;
acc  = 1e-12;
k_id = max(find(err_ID > (ss(1)*acc)));
han = plot(xxs(1,:),xxs(2,:),'k.',...
           a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,ind(1:k_id)),xxs(2,ind(1:k_id)),'r.',...
           xxt(1,:),xxt(2,:),'b.',...
           2 + 0.5*[-1,1,1,-1,-1],0.5*[-1,-1,1,1,-1],'k');
set(han(1),'MarkerSize',20)
set(han(1),'MarkerEdgeColor',0.7*[1,1,1])
set(han(2),'LineWidth',2)
set(han(3),'MarkerSize',20)
set(han(4),'MarkerSize',20)
set(han(5),'LineWidth',2)
axis equal
axis off
title(sprintf('Skeleton points: eps=%0.1e  k=%d',acc,k_id))
keyboard
    
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_disc_to_circle_quad

rng(0)

cc     = [0;0]; % Center of geometry.
a_in   = 0.5;   % Radius of source circle.
rad    = 0.75;   % Radius of target circle.
NT     = 200;   % Number of points along the outer contour.
k_max  = 80;    % Max rank to be investigated.

%%% Set up the source geometry.
n_rad      = 20;
n_ang      = 60;
[rr,wr]    = LOCAL_clencurt(n_rad,0,a_in);
tt         = linspace(0,2*pi*(1 - 1/n_ang),n_ang);
wt         = (2*pi/n_ang)*ones(1,n_ang);
[RR,TT]    = meshgrid(rr,tt);
XX1        = RR.*cos(TT);
XX2        = RR.*sin(TT);
WW         = RR.*(wt' * wr');
NS         = numel(WW);
xxs        = [reshape(XX1,1,NS);...
              reshape(XX2,1,NS)];
wws        = reshape(WW,1,NS);
box_geom_src = [2*a_in,cc(1),cc(2)];

%%% Set up the target geometry.
tt   = linspace(0,2*pi*(1 - 1/NT),NT);
xxt  = rad*[cos(tt);sin(tt)];
wwt  = (rad*2*pi/NT)*ones(1,NT);

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
A       = 0.5*log(DD1.*DD1 + DD2.*DD2).*sqrt(wwt' * wws);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the ID
[Q,R,ind] = LOCAL_pivoted_QR3(A,k_max);
err_ID    = [ss(1),NaN*ones(1,k_max)];
for i = 2:length(kk)
  k         = kk(i);
  err_ID(i) = norm(A(:,ind) - Q(:,1:k)*R(1:k,:));
end

%%% Perform the multipole expansions.
err_mpole = [ss(1),NaN*ones(1,k_max)];
P_max     = floor(k_max/2)-1;
CC        = LOCAL_T_ofs(box_geom_src,P_max,xxs);
BB        = LOCAL_T_tfo(xxt,box_geom_src,P_max);
for p = 0:P_max
  A_approx = real(BB(:,1:(p+1))*CC(1:(p+1),:)).*sqrt(wwt' * wws);
  err_mpole([2*p+2,2*p+3]) = norm(A - A_approx);
end

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
tt = linspace(0,2*pi,500);
han = plot(a_in*cos(tt),a_in*sin(tt),'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'MarkerSize',20)
set(han(3),'MarkerSize',20)
axis equal
axis off
title(sprintf('The geometry.'))

%%% Plot the errors.
figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'x-',...
         kk,(1/ss(1))*err_ID,'x-',...
         kk,(1/ss(1))*err_mpole,'o-');
legend('svd','ID','mpole')
title(sprintf('Expansion errors.'))
xlabel('k')
ylabel('||A - A_k||/||A||')

%%% Plot a skeleton
figure(1); subplot(1,1,1); hold off;
acc  = 1e-12;
k_id = max(find(err_ID > (ss(1)*acc)));
han = plot(xxs(1,:),xxs(2,:),'k.',...
           a_in*cos(tt),a_in*sin(tt),'k',...
           xxs(1,ind(1:k_id)),xxs(2,ind(1:k_id)),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'MarkerSize',20)
set(han(1),'MarkerEdgeColor',0.7*[1,1,1])
set(han(2),'LineWidth',2)
set(han(3),'MarkerSize',20)
set(han(4),'MarkerSize',20)
axis equal
axis off
fprintf('Skeleton points: eps=%0.1e  k=%d\n',acc,k_id)

keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_box_to_circle_quad

rng(0)

cc     = [0;0]; % Center of geometry.
rad    = 1.5;   % Target circle has radius "rad".
a_in   = 0.5;   % Inner box is of size 2*a_in x 2*a_in.
NT     = 200;   % Number of points along the outer contour.
n      = 20;    % Number of Chebyshev nodes along a side of inner box.
k_max  = 80;        % Max rank to be investigated.

%%% Set up the source geometry.
[x,w]      = LOCAL_clencurt(n,-a_in,a_in);
x          = x(end:(-1):1)';
w          = w';
[XX1,XX2]  = meshgrid(x);
xxs        = [reshape(XX1,1,n*n);...
              reshape(XX2,1,n*n)];
wws        = reshape(w' * w,1,n*n);
NS         = size(xxs,2);
box_geom_src = [2*a_in,cc(1),cc(2)];

%%% Set up the target geometry.
tt   = linspace(0,2*pi*(1 - 1/NT),NT);
xxt  = rad*[cos(tt);sin(tt)];
wwt  = (rad*2*pi/NT)*ones(1,NT);

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
A       = 0.5*log(DD1.*DD1 + DD2.*DD2).*sqrt(wwt' * wws);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the ID
[Q,R,ind] = LOCAL_pivoted_QR3(A,k_max);
err_ID    = [ss(1),NaN*ones(1,k_max)];
for i = 2:length(kk)
  k         = kk(i);
  err_ID(i) = norm(A(:,ind) - Q(:,1:k)*R(1:k,:));
end

%%% Perform the multipole expansions.
err_mpole = [ss(1),NaN*ones(1,k_max)];
P_max     = floor(k_max/2)-1;
CC        = LOCAL_T_ofs(box_geom_src,P_max,xxs);
BB        = LOCAL_T_tfo(xxt,box_geom_src,P_max);
for p = 0:P_max
  A_approx = real(BB(:,1:(p+1))*CC(1:(p+1),:)).*sqrt(wwt' * wws);
  err_mpole([2*p+2,2*p+3]) = norm(A - A_approx);
end

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
han = plot(a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'MarkerSize',20)
set(han(3),'MarkerSize',20)
axis equal
axis off
title(sprintf('The geometry.'))

%%% Plot the errors.
figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'x-',...
         kk,(1/ss(1))*err_ID,'x-',...
         kk,(1/ss(1))*err_mpole,'x-');
legend('svd','ID','mpole')
title(sprintf('Expansion errors.'))
xlabel('k')
ylabel('||A - A_k||/||A||')

%%% Plot a skeleton
figure(3); subplot(1,1,1); hold off;
acc  = 1e-10;
k_id = max(find(err_ID > (ss(1)*acc)));
han = plot(xxs(1,:),xxs(2,:),'k.',...
           a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,ind(1:k_id)),xxs(2,ind(1:k_id)),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'MarkerSize',20)
set(han(1),'MarkerEdgeColor',0.7*[1,1,1])
set(han(2),'LineWidth',2)
set(han(3),'MarkerSize',20)
set(han(4),'MarkerSize',20)
axis equal
axis off
title(sprintf('Skeleton points: eps=%0.1e  k=%d',acc,k_id))
keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_box_to_line_quad

rng(0)

cc     = [0;0]; % Center of geometry.
a_in   = 0.5;   % Inner box is of size 2*a_in x 2*a_in.
a_out  = 1.5;   % Outer box is of size 2*a_out x 2*a_out.
n      = 30;    % Number of Chebyshev nodes along a side of inner box.
n_out  = 100;   % Number of Chebyshev nodes along a side of outer box.
k_max  = 80;    % Max rank to be investigated.

%%% Set up the source geometry.
[x,w]      = LOCAL_clencurt(n,-a_in,a_in);
x          = x(end:(-1):1)';
w          = w';
[XX1,XX2]  = meshgrid(x);
xxs        = [reshape(XX1,1,n*n);...
              reshape(XX2,1,n*n)];
wws        = reshape(w' * w,1,n*n);
NS         = size(xxs,2);
box_geom_src = [2*a_in,cc(1),cc(2)];

%%% Set up the target geometry.
[x,w] = LOCAL_clencurt(n_out,-a_out,a_out);
x     = x(end:(-1):1)';
w     = w';
xxt   = [x,a_out*ones(1,n_out-2),x(n_out:(-1):1),-a_out*ones(1,n_out-2);...
         -a_out*ones(1,n_out-1),x,a_out*ones(1,n_out-2),x(n_out:(-1):2)];
wwt   = [2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1))];
NT    = size(xxt,2);

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
A       = 0.5*log(DD1.*DD1 + DD2.*DD2).*sqrt(wwt' * wws);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the ID
[Q,R,ind] = LOCAL_pivoted_QR3(A,k_max);
err_ID    = [ss(1),NaN*ones(1,k_max)];
for i = 2:length(kk)
  k         = kk(i);
  err_ID(i) = norm(A(:,ind) - Q(:,1:k)*R(1:k,:));
end

%%% Perform the multipole expansions.
err_mpole = [ss(1),NaN*ones(1,k_max)];
P_max     = floor(k_max/2)-1;
CC        = LOCAL_T_ofs(box_geom_src,P_max,xxs);
BB        = LOCAL_T_tfo(xxt,box_geom_src,P_max);
for p = 0:P_max
  A_approx = real(BB(:,1:(p+1))*CC(1:(p+1),:)).*sqrt(wwt' * wws);
  err_mpole([2*p+2,2*p+3]) = norm(A - A_approx);
end

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
han = plot(a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'MarkerSize',20)
set(han(3),'MarkerSize',20)
axis equal
axis off
title(sprintf('The geometry,'))

%%% Plot the errors.
figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'x-',...
         kk,(1/ss(1))*err_ID,'x-',...
         kk,(1/ss(1))*err_mpole,'x-');
legend('svd','ID','mpole')
title(sprintf('Expansion errors,'))
xlabel('k')
ylabel('||A - A_k||/||A||')
axis([0,150,1e-15,2])

%%% Plot a skeleton
figure(1); subplot(1,1,1); hold off;
acc  = 1e-12;
k_id = max(find(err_ID > (ss(1)*acc)));
han = plot(xxs(1,:),xxs(2,:),'k.',...
           a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,ind(1:k_id)),xxs(2,ind(1:k_id)),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'MarkerSize',20)
set(han(1),'MarkerEdgeColor',0.7*[1,1,1])
set(han(2),'LineWidth',2)
set(han(3),'MarkerSize',20)
set(han(4),'MarkerSize',20)
axis equal
axis off
title(sprintf('Skeleton points: eps=%0.1e  k=%d',acc,k_id))
fprintf('Skeleton points: eps=%0.1e  k=%d\n',acc,k_id)

keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_curve_to_line_quad

rng(0)

cc     = [0;0]; % Center of geometry.
a_in   = 0.5;   % Inner box is of size 2*a_in x 2*a_in.
a_out  = 1.5;   % Outer box is of size 2*a_out x 2*a_out.
NS     = 100;   % Number of Chebyshev nodes along contour.
n_out  = 100;   % Number of Chebyshev nodes along a side of outer box.
k_max  = 40;    % Max rank to be investigated.

%%% Set up the source geometry.
[xx1,w]    = LOCAL_clencurt(NS,-a_in,a_in);
xx1        = xx1(end:(-1):1)';
w          = w';
xx2        = 0.1*exp(xx1).*sin(8*xx1);
dxx2       = 0.1*exp(xx1).*sin(8*xx1) + 0.8*exp(xx1).*cos(8*xx1);
xxs        = [xx1;xx2];
wws        = w.*sqrt(1 + dxx2.*dxx2);
box_geom_src = [2*a_in,cc(1),cc(2)];

%%% Set up the target geometry.
[x,w] = LOCAL_clencurt(n_out,-a_out,a_out);
x     = x(end:(-1):1)';
w     = w';
xxt   = [x,a_out*ones(1,n_out-2),x(n_out:(-1):1),-a_out*ones(1,n_out-2);...
         -a_out*ones(1,n_out-1),x,a_out*ones(1,n_out-2),x(n_out:(-1):2)];
wwt   = [2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1))];
NT    = size(xxt,2);

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
A       = 0.5*log(DD1.*DD1 + DD2.*DD2).*sqrt(wwt' * wws);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the ID
[Q,R,ind] = LOCAL_pivoted_QR3(A,k_max);
err_ID    = [ss(1),NaN*ones(1,k_max)];
for i = 2:length(kk)
  k         = kk(i);
  err_ID(i) = norm(A(:,ind) - Q(:,1:k)*R(1:k,:));
end

%%% Perform the multipole expansions.
err_mpole = [ss(1),NaN*ones(1,k_max)];
P_max     = floor(k_max/2)-1;
CC        = LOCAL_T_ofs(box_geom_src,P_max,xxs);
BB        = LOCAL_T_tfo(xxt,box_geom_src,P_max);
for p = 0:P_max
  A_approx = real(BB(:,1:(p+1))*CC(1:(p+1),:)).*sqrt(wwt' * wws);
  err_mpole([2*p+2,2*p+3]) = norm(A - A_approx);
end

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
han = plot(a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'MarkerSize',20)
set(han(3),'MarkerSize',20)
axis equal
axis off
title(sprintf('The geometry,'))

%%% Plot the errors.
figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'x-',...
         kk,(1/ss(1))*err_ID,'x-',...
         kk,(1/ss(1))*err_mpole,'x-');
legend('svd','ID','mpole')
title(sprintf('Expansion errors,'))
xlabel('k')
ylabel('||A - A_k||/||A||')

%%% Plot a skeleton
figure(1); subplot(1,1,1); hold off;
acc  = 1e-12;
k_id = max(find(err_ID > (ss(1)*acc)));
han = plot(xxs(1,:),xxs(2,:),'k.',...
           a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,ind(1:k_id)),xxs(2,ind(1:k_id)),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'MarkerSize',20)
set(han(1),'MarkerEdgeColor',0.7*[1,1,1])
set(han(2),'LineWidth',2)
set(han(3),'MarkerSize',20)
set(han(4),'MarkerSize',20)
axis equal
axis off
title(sprintf('Skeleton points: eps=%0.1e  k=%d',acc,k_id))
fprintf('Skeleton points: eps=%0.1e  k=%d\n',acc,k_id)

keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_box_to_line_randomized

cc     = [0;0]; % Center of geometry.
a_in   = 0.5;   % Inner box is of size 2*a_in x 2*a_in.
a_out  = 1.5;   % Outer box is of size 2*a_out x 2*a_out.
n      = 25;    % Number of Chebyshev nodes along a side of inner box.
n_out  = 60;    % Number of Chebyshev nodes along a side of outer box.
k_max  = 80;    % Max rank to be investigated.
p      = 10;    % Over-sampling parameter in randomized method.

%%% Set up the source geometry.
[x,w]      = LOCAL_clencurt(n,-a_in,a_in);
x          = x(end:(-1):1)';
w          = w';
[XX1,XX2]  = meshgrid(x);
xxs        = [reshape(XX1,1,n*n);...
              reshape(XX2,1,n*n)];
wws        = reshape(w' * w,1,n*n);
NS         = size(xxs,2);
box_geom_src = [2*a_in,cc(1),cc(2)];

%%% Set up the target geometry.
[x,w] = LOCAL_clencurt(n_out,-a_out,a_out);
x     = x(end:(-1):1)';
w     = w';
xxt   = [x,a_out*ones(1,n_out-2),x(n_out:(-1):1),-a_out*ones(1,n_out-2);...
         -a_out*ones(1,n_out-1),x,a_out*ones(1,n_out-2),x(n_out:(-1):2)];
wwt   = [2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1))];
NT    = size(xxt,2);

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
A       = 0.5*log(DD1.*DD1 + DD2.*DD2).*sqrt(wwt' * wws);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the randomized sampling.
n_exp = 20; % Number of experiments.
err_rand = [ss(1),NaN*ones(1,k_max)];
ERR = err_rand'*ones(1,n_exp);
for j = 1:size(ERR,2)
  Omega = randn(NS,k_max + p);
  Y     = A*Omega;
  for i = 2:length(kk)
    k = kk(i);
    [Q,~] = qr(Y(:,1:k),0);
    ERR(i,j) = norm(A - Q*(Q'*A));
  end
end

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
han = plot(a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'MarkerSize',20)
set(han(3),'MarkerSize',20)
axis equal
axis off
title(sprintf('The geometry,'))

%%% Plot the errors.
figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'x-',...
         kk,(1/ss(1))*ERR,'r-');
legend('svd','randomized')
title(sprintf('Expansion errors,'))
xlabel('k')
ylabel('||A - A_k||/||A||')

keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The outgoing-from-sources translation operator:
%   T(  1,j) = 1;
%   T(p+1,j) = -(1/p)*(xxj - cc)^p
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function T = LOCAL_T_ofs(box_geom_src,p_fmm,xxsrc)

nsrc = size(xxsrc,2);
DD   = ones(p_fmm,1)*((xxsrc(1,:)-box_geom_src(2)) + 1i*(xxsrc(2,:) - box_geom_src(3)));
PP   = ((1:p_fmm)')*ones(1,nsrc);
T    = [ ones(1,nsrc); -(1./PP).*(DD.^PP)];
  
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The targets-from-outgoing expansion.
%   T(i,1)   = log(xxi - cc)
%   T(i,p+1) = (1/(xxi - cc))^p
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function T = LOCAL_T_tfo(xxtrg,box_geom,p_fmm)

ntrg = size(xxtrg,2);
PP   = ones(ntrg,1)*(1:p_fmm);
dd   = ((xxtrg(1,:)-box_geom(2)) + 1i*(xxtrg(2,:)-box_geom(3))).';
T    = [log(dd), 1./((dd*ones(1,p_fmm)).^PP)];
  
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Q,R,ind] = LOCAL_pivoted_QR3(A,k)

% This function orthogonalizes the COLUMNS of A
% It uses a modified Gram-Schmidt with column pivoting

m = size(A,1);
n = size(A,2);

R = zeros(min(m,n),n);
Q = A;
ind = 1:n;

for j = 1:k
    [trash, j_max] = max(sum(Q(:,j:n).*Q(:,j:n)));
    j_max          = j_max + j - 1;
    Q(:,[j,j_max]) = Q(:,[j_max,j]);
    R(:,[j,j_max]) = R(:,[j_max,j]);
    ind([j,j_max]) = ind([j_max,j]);
    r_jj   = norm(Q(:,j));
    Q(:,j) = Q(:,j)/r_jj;
    Q(:,j) = Q(:,j) - Q(:,1:(j-1))*(Q(:,1:(j-1))'*Q(:,j));
    Q(:,j) = Q(:,j)/norm(Q(:,j));
%    Q(:,j) = Q(:,j) - Q(:,1:(j-1))*(Q(:,1:(j-1))'*Q(:,j));
%    Q(:,j) = Q(:,j)/norm(Q(:,j));
    R(j,j) = r_jj;
    rr     = Q(:,j)'*Q(:,(j+1):n);
    R(j,(j+1):n) = rr;
    Q(:,(j+1):n) = Q(:,(j+1):n) - Q(:,j)*rr;
end

Q = Q(:, 1:min(m,n));

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the Chebyshev nodes on [-1,1].
% It also computes a differentiation operator.
% It is modified from a code by Nick Trefethen.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [D,x] = LOCAL_cheb(N)

if N==0
  D=0; 
  x=1; 
  return
end
x = cos(pi*(0:N)/N)'; 
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';                  
D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal entries
D  = D - diag(sum(D'));                 % diagonal entries

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,w] = LOCAL_clencurt(N1,a,b)

% Adapted from "fclencurt.m" by Greg von Winckel.
N=N1-1; 
bma=b-a;
c=zeros(N1,2);
c(1:2:N1,1)=(2./[1 1-(2:2:N).^2 ])'; 
c(2,2)=1;
f=real(ifft([c(1:N1,:);c(N:-1:2,:)]));
w=bma*([f(1,1); 2*f(2:N,1); f(N1,1)])/2;
x=0.5*((b+a)+N*bma*f(1:N1,2));

return

