function val=dothemall
% this plots the power function for harmonic interpolation
clear all;
close all;
global datafun      % steers the data d+function   
global domaintype   % steers the domain
global harmkertype  % steers harmonic kernel, exp and poi 
global cfact        % c factor in harmonic kernel
global Ctest        % number of test points on the boundary
global Npoints      % number of interpolation points on the boundary

% we clear the table file, if it exists 
fid = fopen('table.dat','w+');
fclose(fid);
% global controls for ALL cases
% Ctest=1000;
Npoints=20;
% now we do several cases

harmkertype='exp';
datafun='y'; 
domaintype='card';
powerfct;
return


function powerfct
global datafun
global domaintype
global harmkertype
global cfact
global xpoint        % FS point for data function 3 
global ypoint
global Ctest
global Npoints

% tp=getbndpoints(Ctest,0);
% [Ctest cd]=size(tp);
% dtp=datfct(tp);
% [tht rt]=cart2pol(tp(:,1),tp(:,2));
% cp=getbndpoints(300,0);
% [Mcoll cd]=size(cp);
% dcp=datfct(cp);
% [thc rc]=cart2pol(cp(:,1),cp(:,2));
% bval=dcp;

% hlist=[];
% errlist=[];
% nlist=[];
% clist=[];

sp=getbndpoints(Npoints,0);
h=norm(sp(1,:)-sp(2,:));
% cfact=2*h;
sp=sp(1:end-1,:);
[Npoints nd]=size(sp);
[ths rs]=cart2pol(sp(:,1),sp(:,2));

rmax=max(rs(:));
if harmkertype=='poi'
    cfact=0.85/rmax;
end
if harmkertype=='exp'
    cfact=1/rmax;
end
% 
% errlist=[];
% nlist=[];
% clist=[];
he=0.01;
[xe, ye]=meshgrid(-1.1:he:1.1,-1.1:he:1.1);
xee=xe(:);
yee=ye(:);
[the, re]=cart2pol(xee,yee);
rhoe=bndcurve(the);
ifind=find(re<=rhoe);
mask=zeros(size(xe));
msk=mask(:);
msk(ifind)=1.0;
figure
% subplot(2,1,1)
% plot(xee(ifind),ye(ifind),'.')
% hold on
plot(sp(:,1),sp(:,2),'o');
% title('Interpolation Points')
axis([-1.2 0.7 -1.2 1.2])
thee=the(ifind);
ree=re(ifind);
while 1==1
    imat=getharmker(ths,rs,ths,rs);
    si=size(imat)
    ce=condest(imat)
    emat=getharmker(thee,ree,ths,rs);
    se=size(emat)
    if ce<1.e12
        break;
    end
    cfact=cfact*1.2
    pause
end
u=imat\emat;
su=size(u)

for i=1:length(thee)
    powfct(i)=getharmker(thee(i),ree(i),thee(i),ree(i));  
end
sp=size(powfct)
pf=powfct'-2*sum(u.*emat)'+sum(u.*(imat*u))';
val=zeros(size(xee));
val(ifind)=pf;
figure
% subplot(2,1,2)
surfc(xe,ye,reshape(val,size(xe)))
axis([-1 1  -1 1 -0.5*max(abs(val)) max(abs(val))])
% title('Power Function for Harmonic Exponential Kernel')
shading interp
return

function vmat=getevalker(tha,ra,thb,rb);
RBFpar=1;
[xa ya]=pol2cart(tha, ra);
[xb yb]=pol2cart(thb, rb);
p=[xb yb];
q=[xa ya];
[np pdim]=size(p)
[nq qdim]=size(q)
dst=p*q';
dsi=size(dst)
cp=sum((p.*p)'); % squared norms of p points, as row, halved
cq=sum((q.*q)'); % squared norms of q points, as row, halved
dst=sqrt(repmat(cp',1,nq)+repmat(cq,np,1)-2*dst); 
vmat=besselk(RBFpar,dst+eps).*(dst+eps).^RBFpar;
return


function singlecase
global datafun
global domaintype
global harmkertype
global cfact
global xpoint        % FS point for data function 3 
global ypoint
global Ctest
global Npoints

tp=getbndpoints(Ctest,0);
[Ctest cd]=size(tp);
dtp=datfct(tp);
[tht rt]=cart2pol(tp(:,1),tp(:,2));
rmax=max(rt(:));
if harmkertype=='poi'
    cfact=0.85/rmax;
end
if harmkertype=='exp'
    cfact=1/rmax;
end
cp=getbndpoints(300,0);
[Mcoll cd]=size(cp);
dcp=datfct(cp);
[thc rc]=cart2pol(cp(:,1),cp(:,2));
bval=dcp;

hlist=[];
errlist=[];
nlist=[];
clist=[];

sp=getbndpoints(Npoints,0);
h=norm(sp(1,:)-sp(2,:));
% cfact=2*h;
sp=sp(1:end-1,:);
[Npoints nd]=size(sp);
[ths rs]=cart2pol(sp(:,1),sp(:,2));
errlist=[];
nlist=[];
clist=[];

mat=getharmker(ths,rs,thc,rc);
indlist=myselector(mat,bval);    
sv=svd(mat(:,indlist));

for n=1:length(indlist)   
    ce=sv(1)/sv(n);
    coeff=mat(:,indlist(1:n))\bval;
    matt=getharmker(ths(indlist(1:n),1),rs(indlist(1:n),1),tht,rt);
    res=matt*coeff;
    ner=norm(dtp-res,Inf);
    clist=[clist ce];
    errlist=[errlist ner];
    nlist=[nlist n];
end
mydisplay(domaintype, datafun, harmkertype, Npoints, n, ner)

% subplot(2,2,3)
% semilogy(nlist,errlist)
% xlabel('N')
% title('Error')
% subplot(2,2,4)
% plot(sp(:,1),sp(:,2),'.',sp(indlist,1),sp(indlist,2),'o')
% title('Points')
% subplot(2,2,2)
% plot(tht,dtp-res)
% title('Last Error')
% subplot(2,2,1)
% plot(tht,dtp)
% title('Data')
return

function points=getbndpoints(n, dist)
% construct set of n points at distance dist of boundary
points=curveval(n,dist);
return

function points=curveval(n,dist)
% gets n points on periodic curve, with double endpoint
% and addition of radial distance dist 
global bndrand       % controls random shifts on source boundary
hsc=2*pi/n-eps;
the=-pi:hsc:pi;
if (bndrand==1)
    % bndrand
    the=the+rand(size(the))*hsc;
    toolarge=find(the>pi);
    the(toolarge)=the(toolarge)-2*pi;
end
rho=bndcurve(the)+dist;
[xp yp]=pol2cart(the,rho);
points=[xp' yp'];
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rho=bndcurve(the)
% evaluates periodic curve in polar coordinates
global domaintype
switch lower(domaintype)
    case 'circ'
        rho=ones(size(the));
    case 'amoeb'
        rho=0.35*ones(size(the)).*(exp(sin(the)).*sin(2*the).^2+...
            exp(cos(the)).*cos(2*the).^2);
    case 'defo'
        rho=ones(size(the)).*(1+abs(the))/4.2;
    case 'card'
        c=0.5*sqrt(tan(pi/8));
        rho=ones(size(the)).*sqrt(c.^2+abs(sin(the/2)))/sqrt(c.^2+1);
    case 'leaf'
        rho=ones(size(the)).*sqrt(1+0.5*sin(4*the))/sqrt(2);
        
    case 'rect'
        ri=find(abs(the)<=pi/4);
        le=find(abs(the)>=3*pi/4);
        to=find(abs(the-pi/2)<pi/4);
        lo=find(abs(the+pi/2)<pi/4);
        % ss=sqrt(2);
        ss=1;
        rho(ri)=ones(size(the(ri)))./abs(cos(the(ri)))/ss;
        rho(le)=ones(size(the(le)))./abs(cos(the(le)))/ss;
        rho(to)=ones(size(the(to)))./abs(sin(the(to)))/ss;
        rho(lo)=ones(size(the(lo)))./abs(sin(the(lo)))/ss;
        return
    otherwise
        error('domain not known')
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val=datfct(points)
global datafun
global xpoint
global ypoint
switch lower(datafun)
    case 'abs'
        val=abs(points(:,2));
    case 'y'
        val=points(:,2);
    case 'trunc'
        val=max(0,points(:,2));
    case 'poly'    
        val=points(:,2).^3.*points(:,1).^2;
    case 'harm'
        val=cos(points(:,2)).*exp(points(:,1));
    case 'gauss'
        val=exp(-points(:,2).^2-points(:,1).^2);
    otherwise
        error('data function not known')
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mat=getharmker(ths,rs,thc,rc)
% gets harmonic kernel matrix from two sets of
% points in polar coordinates
global harmkertype
global cfact
Mcoll=length(thc);
Npoints=length(ths);
rcs=cfact^2*rc*rs';
tdif=repmat(ths,1,Mcoll)-repmat(thc',Npoints,1);
rc=rcs.*cos(tdif');
switch lower(harmkertype)
    case 'exp'
        mat=cos(rcs.*sin(tdif')).*exp(rcs.*cos(tdif'));
    case 'poi'
        mat=(1-rc)./(1-2*rc+rcs.*rcs);
    otherwise
        error('Unknown kernel')
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function indlist=myselector(A,bval);
disp('selector started')
[m L]=size(A);
Aold=A;
bvalold=bval;
p=zeros(1,L);  % squares of norms of columns
indlist=[];
for k=1:L
    for i=1:L      % get column norms squared, takes m*L
        p(1,i)=A(:,i)'*A(:,i);
    end
    plist=find(p>eps);
    if (length(p)==0)
        disp('all p are zero');
        break;
    end
    f=bval'*A; % scalar products of rhs with columns
    [fm ifk]=max(f(1,plist).^2./(eps+p(1,plist))); 
    ifm=plist(ifk);    
    if (length(ifm)==0)
        disp('no useful column found')
        break;
    end
    u=A(:,ifm);
    v=u/norm(u);
    A=A-v*(v'*A);
    A(:,ifm)=0;   
    bvalsq=bval'*bval;
    bval=bval-u*f(1,ifm)/p(1,ifm);
    bvalnsq=bval'*bval;
    if (bvalnsq>bvalsq-eps)
        disp('no absolute reduction')
        break;
    end
    if (bvalnsq>bvalsq*0.99999)        
        disp('no relative reduction')
        break;
    end
    indlist=[indlist ifm];
end
disp('selector ended')
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mydisplay(domaintype, datafun, harmkertype, Npoints, n, ner)
str=sprintf('%s & %s & %s & %d & %d & %e\\\\',...
    domaintype, datafun, harmkertype, Npoints, n, ner);
disp(str);
fid = fopen('table.dat','a');
fprintf(fid,'%s\n',str);
fclose(fid);
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function harmkerint
global datafun      % steers the data d+function   
global domaintype   % steers the domain
global harmkertype  % steers harmonic kernel, exp and poi 
global cfact        % c factor in harmonic kernel
global Ctest        % number of test points on the boudary
global Npoints      % number of interpolation points on the bondary

[xs ys]=meshgrid(-1:0.1:1,-1:0.1:1);
tp=[xs(:) ys(:)];
[Ctest cd]=size(tp);
dtp=datfct(tp);
[tht rt]=cart2pol(tp(:,1),tp(:,2));
rmax=max(rt(:));
if harmkertype=='poi'
    cfact=0.85/rmax;
end
if harmkertype=='exp'
    cfact=1/rmax;
end
hlist=[];
errlist=[];
nlist=[];
clist=[];

sp=2*rand(Npoints,2)-1;
[Npoints nd]=size(sp);
[ths rs]=cart2pol(sp(:,1),sp(:,2));
bval=datfct(sp);

errlist=[];
nlist=[];
clist=[];

mat=getharmker(ths,rs,ths,rs);
indlist=myselector(mat,bval);    
sv=svd(mat(:,indlist));

for n=1:length(indlist)   
    ce=sv(1)/sv(n);
    coeff=mat(:,indlist(1:n))\bval;
    matt=getharmker(ths(indlist(1:n),1),rs(indlist(1:n),1),tht,rt);
    res=matt*coeff;
    ner=norm(dtp-res,Inf);
    clist=[clist ce];
    errlist=[errlist ner];
    nlist=[nlist n];
end
subplot(2,2,1)
surf(xs,ys,reshape(dtp,size(xs)))
% axis([-1 1 -1 1 -1 1])
title('Data')
subplot(2,2,2)
surf(xs,ys,reshape(res,size(xs)))
title('Interpolant')
subplot(2,2,4)
plot(sp(indlist,1),sp(indlist,2),'ro',sp(:,1),sp(:,2),'.')
title('Data Points')
subplot(2,2,3)
surf(xs,ys,reshape(dtp-res,size(xs)))
title('Error')
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plotsing
global datafun      % steers the data d+function   
global domaintype   % steers the domain
global harmkertype  % steers harmonic kernel, exp and poi 
global cfact        % c factor in harmonic kernel
global Ctest        % number of test points on the boudary
global Npoints      % number of interpolation points on the bondary

domaintype='rect';
tp=getbndpoints(Ctest,0);
[Ctest cd]=size(tp);
[tht rt]=cart2pol(tp(:,1),tp(:,2));
rmax=max(abs(rt(:)));
cfact=0.95/rmax;
st=(1./rt(:))/cfact^2;
[xx yy]=pol2cart(tht,st);
subplot(1,2,1)
plot(tp(:,1),tp(:,2),'r.',xx,yy,'.')
daspect([1 1 1]) 
axis([-2.5 2.5 -2.5 2.5]);
title('Reflected points for rectangle')

domaintype='card';
subplot(1,2,2)
tp=getbndpoints(Ctest,0);
[Ctest cd]=size(tp);
[tht rt]=cart2pol(tp(:,1),tp(:,2));
rmax=max(abs(rt(:)));
cfact=0.95/rmax;
st=(1./rt(:))/cfact^2;
[xx yy]=pol2cart(tht,st);
plot(tp(:,1),tp(:,2),'r.',xx,yy,'.')
title('Reflected points for cardioid')
daspect([1 1 1]) 
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mat=getharmkerlocrec(ths,rs,thc,rc,c)
% gets harmonic kernel matrix from two sets of
% points in polar coordinates, with local rescaling
global harmkertype
Mcoll=length(thc);
Npoints=length(ths);
rcs=rc*(c.^2.*rs)';
tdif=repmat(ths,1,Mcoll)-repmat(thc',Npoints,1);
rc=rcs.*cos(tdif');
switch lower(harmkertype)
    case 'exp'
        mat=cos(rcs.*sin(tdif')).*exp(rcs.*cos(tdif'));
    case 'poi'
        mat=(1-rc)./(1-2*rc+rcs.*rcs);
    otherwise
        error('Unknown kernel')
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function singlecaselocrec
global datafun      % steers the data d+function   
global domaintype   % steers the domain
global harmkertype  % steers harmonic kernel, exp and poi 
global cfact        % c factor in harmonic kernel
global Ctest        % number of test points on the boudary
global Npoints      % number of interpolation points on the bondary

tp=getbndpoints(Ctest,0);
[Ctest cd]=size(tp);
dtp=datfct(tp);
[tht rt]=cart2pol(tp(:,1),tp(:,2));
rmax=max(rt(:));
cp=getbndpoints(300,0);
[Mcoll cd]=size(cp);
dcp=datfct(cp);
[thc rc]=cart2pol(cp(:,1),cp(:,2));
bval=dcp;

hlist=[];
errlist=[];
nlist=[];
clist=[];

sp=getbndpoints(Npoints,0);
h=norm(sp(1,:)-sp(2,:));
% cfact=2*h;
sp=sp(1:end-1,:);
[Npoints nd]=size(sp);
[ths rs]=cart2pol(sp(:,1),sp(:,2));
errlist=[];
nlist=[];
clist=[];
c=ones(size(rs));
reflr=ones(size(rs));
for i=1:length(rs)
    
    if harmkertype=='poi'
        if  mod(i,2)==0
            c(i)=0.95./rs(i);
        else
            c(i)=0.85/rmax;
        end
    end
    if harmkertype=='exp'
        if (mod(i,2)==0)
            c(i)=1./rs(i);
        else
            c(i)=1./rmax;
        end
    end
    reflr(i)=(1./rs(i))/c(i)^2;
end

mat=getharmkerlocrec(ths,rs,thc,rc,c);
indlist=myselector(mat,bval);    
sv=svd(mat(:,indlist));

for n=1:length(indlist)   
    ce=sv(1)/sv(n);
    coeff=mat(:,indlist(1:n))\bval;
    matt=getharmkerlocrec...
        (ths(indlist(1:n),1),rs(indlist(1:n),1),tht,rt,c(indlist(1:n)));
    res=matt*coeff;
    ner=norm(dtp-res,Inf);
    clist=[clist ce];
    errlist=[errlist ner];
    nlist=[nlist n];
end
mydisplaylocrec(domaintype, datafun, harmkertype, Npoints, n, ner)
subplot(2,1,2)
plot(tht,dtp-res)
title('Error')
subplot(2,1,1)
plot(tht,dtp)
title('Data')
figure(5)
[xr yr]=pol2cart(ths,reflr);
plot(cp(:,1),cp(:,2),'.',xr,yr,'.',xr(indlist),yr(indlist),'o')
title('Reflected points for selected interpolation points') 
% sp(indlist,1),sp(indlist,2),'o',
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mydisplaylocrec(domaintype, datafun, harmkertype, Npoints, n, ner)
str=sprintf('%s & %s & %s & %d & %d & %e locrec\\\\',...
    domaintype, datafun, harmkertype, Npoints, n, ner);
disp(str);
fid = fopen('table.dat','a');
fprintf(fid,'%s\n',str);
fclose(fid);
return
