clear all;
close all;
% Calculate power functions for polynomial interpolation
% We need the RBF parameters only for error evaluation
global RBFtype; 
global RBFpar;
global RBFscale;
global A
RBFscale=1   % 1 for unscaled Sobolev space
RBFtype='ms' % take 'ms' for Sobolev-Matern kernel
soborder=4;
RBFpar=soborder-1;  % in 2D we then are in Sobolev space RBFpar+d/2 = RBFpar+1
order=7;     % this presribes the polynomial order
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evaluation grid
he=0.02;
[xe, ye]=meshgrid(-1:he:1,-1:he:1);
Xe=[xe(:), ye(:)];
% restricted to circle 
inda=find(xe.^2+ye.^2<=1);
Xea=Xe(inda,:);
% and cut out the third quadrant 
indb=find((Xea(:,1)>=0) | (Xea(:,2)>=0));
indc=inda(indb);
Xec=Xe(indc,:);
% figure
% plot(Xec(:,1),Xec(:,2),'.')
nXec=length(Xec(:,1));
% We take the best points for a given polynomial order
pts=Xec;
polvalues=polynomials(pts,order);
[np no]=size(polvalues)
% We do a pivoted LU factorization 
% and use the points that correspond to row selection.
% This gives a Newton type basis
[L, U, P]=lu(polvalues);
PPQ=P*(1:np)';
pq=PPQ(1:no,1);
PpL=P'*L;
val=zeros(size(xe(:)));
% We now want Lagrange basis functions
Pp=P'*polvalues;
Pploc=Pp(1:no,1:no);
Lp=polvalues*inv(Pploc);
% 
figure
subplot(2,2,3)
spts=P*pts;
phi=-pi/2:0.01:pi;
rad=ones(size(phi));
[xc yc]=pol2cart(phi,rad);
plot(xc,yc)
hold on
xc=-1:0.01:0;
plot(xc,zeros(size(xc)))
hold on
xc=-1:0.01:0;
plot(zeros(size(xc)),xc)
hold on
plot(spts(1:no,1),spts(1:no,2),'ro')

axis([-1.2 1.2 -1.2 1.2])
daspect([1 1 1])
% title(sprintf('Selected %d data points',no))

A=kermat(spts(1:no,:),spts(1:no,:));
iA=inv(A);
Ae=kermat(Xec,spts(1:no,:));
uu=Ae*iA;
vv=uu.*Ae;
vvv=sum(vv'.*uu')';
vw=sum(vv')';
pow=frbf(0,0)-vw;
val(indc)=pow;
mval=min(val(:))
subplot(2,2,2)
surf(xe,ye,reshape(sqrt((val-mval)),size(xe)));
shading interp
% title('Optimal Power Function')

vv=Lp.*Ae;
vw=sum(vv')';
w=(Lp*A).*Lp;
ww=sum(w')';
gpow=frbf(0,0)-2*vw+ww;
gval=zeros(size(xe(:)));
gval(indc)=gpow;
mgval=min(gval(:))
subplot(2,2,1)
surf(xe,ye,reshape(sqrt((gval-mgval)),size(xe)));
shading interp
% title('Polynomial Power Function')

mdval=min(sqrt(gval-mgval)-sqrt(val-mval))
subplot(2,2,4)
surf(xe,ye,reshape(sqrt(gval-mgval)-sqrt(val-mval)-mdval,size(xe)));
shading interp
%title('Difference of Power Functions')
pmax=max(sqrt(val-mval))
disp('Hit Enter to continue')
pause
% We now plot Newton basis functions
% for i=1:no
%     val(indc)=PpL(:,i);
%     figure
%     surf(xe,ye,reshape(val,size(xe)))
%     % shading interp
%     % title(sprintf('Newton basis function %d',i))
%     shading interp
%     pause(1)
% end
% 

figure
subplot(2,2,1)
val(indc)=PpL(:,9);
surf(xe,ye,reshape(val,size(xe)))
% shading interp
% title(sprintf('Newton basis function %d',i))
shading interp
subplot(2,2,2)
val(indc)=PpL(:,15);
surf(xe,ye,reshape(val,size(xe)))
% shading interp
% title(sprintf('Newton basis function %d',i))
shading interp
subplot(2,2,3)
val(indc)=PpL(:,27);
surf(xe,ye,reshape(val,size(xe)))
% shading interp
% title(sprintf('Newton basis function %d',i))
shading interp
subplot(2,2,4)
val(indc)=PpL(:,28);
surf(xe,ye,reshape(val,size(xe)))
% shading interp
% title(sprintf('Newton basis function %d',i))
shading interp

disp('Hit Enter to continue')
pause
% We now plot Lagrange basis functions
% for i=1:no
%     val(indc)=Lp(:,i);
%     figure
%     surfc(xe,ye,reshape(val,size(xe)))
%     shading interp
%     title(sprintf('Lagrange basis function %d',i))
%     pause(1)
% end


