clear all;
close all;
% Program for P-greedy method using Newton basis.
% No need for solving equations.
% No need for storing large matrices.
global RBFtype; 
global RBFpar;
global RBFscale;
RBFscale=3
RBFtype='mq'
RBFpar=-2
dx=2 % space dimension
dy=dx;
he=0.01; % 1 D stepsize for evaluation grid. 
% The grid can consist of very many points.
% he=0.01 gives Np=40401 points to work on
[xe, ye]=meshgrid(-1:he:1,-1:he:1);
xee=xe(:);
yee=ye(:);
ind=find((xee.^2+yee.^2<=1)&((xee>=0)|(yee>=0)) );
msk=zeros(size(xee));
msk(ind)=1.0;
figure
plot(xee(ind),yee(ind),'.')
title('Discretized domain')
disp('Hit enter to proceed')
pause
Xe=[xee(ind), yee(ind)];
Np=length(xee(ind));
Gtol=1.0e-12
% [v, mi, Lmatrix, qval, hval] =PFGreedyNewtonMovie(Xe, 200, Gtol);
maxpoints=40
[Np dx]=size(Xe);
Lmatrix=zeros(maxpoints,maxpoints); % the A=L*L' decomposition matrix
iter=1 % iteration counter
% We start with the power function for no points. 
% It is the constant K(x,x)=phi(0)
res=kermat(zeros(1,dx),zeros(1,dx))*ones(Np,1); 
% This will be the matrix containing the Newton basis, columnwise,
% as functions on the evaluation grid.
% Since we store L=Lmatrix, we can easily evaluate the 
% Newton basis on other points as well.
v=zeros(Np,maxpoints);
% Now the main iteration. Each step gives a new point 
% and a new basis function.
hv(1)=sqrt(8);
qv(1)=hv(1);
for iter=1:maxpoints
    iter
    %         figure
    %         val=zeros(size(xe));
    %         val(ind)=res;
    %         surf(xe,ye,reshape(val,size(xe)))
    %         shading interp
    %         title(sprintf('Power function for %d points',iter-1))
    %         F(iter)=getframe;
    [maxerr, maxind]=max(abs(res)); % maximum of power function 
    mi(iter,1)=maxind; % the point index where the max. is attained
    p(iter,:)=Xe(maxind,:); % we put the new point z into a list         
    u=kermat(Xe,p(iter,:)); % We need the values K(x,z) of the kernel,
    % where x runs over all evaluation points.
    % This is a VECTOR, not a matrix.
    uu=u(:); % we rearrange the values into a long vector
    % The following is the update formula for A=L*L' 
    % implemented directly.
    lnormsq=0; % this will be the squared norm of the lower L row
    % except for the diagonal. It is zero on startup.
    if iter>1
        hqvec=distsqh(Xe,p(iter,:))';
        hminvec=min(hqvec,hminvec);
        qvec=hqvec(1,mi(1:iter-1,1));
        [qval, iqmin]=min(qvec);
        qv(iter)=sqrt(qval/2);
        hv(iter)=sqrt(2*max(hminvec));
        %         if iqmin<q2h
        %             q2h=qmin;
        %             if q2h<qold/2
        %                 RBFscale=RBFscale/2
        %                 qold=qold/2;
        %             end
        %         end
        uv=uu(mi(1:iter-1,1),1); % select the kernel values for old points
        Lv=backsubL(Lmatrix(1:iter-1,1:iter-1),uv); % solve with old L backwards
        % Note that this does NOT need a solver, just an evaluator.
        % To avoid the backward substitution routine backsubL, here is
        % Lv=Lmatrix(1:iter-1,1:iter-1)\uv; % solve with old L backwards, primitively
        Lmatrix(iter,1:iter-1)=Lv'; % this gives the new lower row of L
        % except for the diagonal element
        lnormsq=Lv'*Lv; % the norm of that row vector
    else
        hminvec=distsqh(Xe,p(1,:))';
    end   
    % whch finally allows us to get the new diagonal element,
    % to finish the L update process:
    Lmatrix(iter,iter)=sqrt(uu(maxind)-lnormsq);
    Linvnorm(iter)=norm(inv(Lmatrix(1:iter,1:iter)));
    % We now calculate the new Newton basis function, using the updated
    % Lmatrix
    v(:,iter)=uu;
    for i=1:iter-1
        v(:,iter)=v(:,iter)-Lmatrix(iter,i)*v(:,i);
    end
    v(:,iter)=v(:,iter)/Lmatrix(iter,iter);
    %     val=zeros(size(xe));
    %     val(ind)=v(:,iter);
    %     figure
    %     surf(xe,ye,reshape(val,size(xe)))
    %     shading interp
    %     % title(sprintf('Newton basis function %d, rescaled',iter))
    %     G(iter)=getframe;
    
    % We now update the power function
    res=res-v(:,iter).^2;
    % and do some statistics
    resnorm(iter)=norm(res,Inf);
    %     resplot(iter)=resnorm;
    % and stop the iteration for error and for condition
    if resnorm(iter)<Gtol
        break;
    end
    ce(iter)=condest(Lmatrix(1:iter,1:iter));
    if ce(iter)>1.0e14
        break;
    end
    % pause(1)
end
% now some selected single plots
figure
subplot(2,2,1)
val=zeros(size(xe));
val(ind)=v(:,21);
surf(xe,ye,reshape(val,size(xe)))
shading interp
% title(sprintf('Newton basis function %d, rescaled',iter))
subplot(2,2,2)
val(ind)=v(:,33);
surf(xe,ye,reshape(val,size(xe)))
shading interp
% title(sprintf('Newton basis function %d, rescaled',iter))
subplot(2,2,3)
val(ind)=v(:,37);
surf(xe,ye,reshape(val,size(xe)))
shading interp
% title(sprintf('Newton basis function %d, rescaled',iter))






% qval=sqrt(2*qval);
% hval=sqrt(2*max(hminvec'));
% figure
% semilogy(1:length(qv),qv,1:length(hv),hv)
% legend('q','h')
% figure
% plot(1:length(hv)-1,hv(2:end)./hv(1:end-1))
% legend('h jumps')
% 

iter=length(mi(:,1))
subplot(2,2,4)

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
axis([-1.2 1.2 -1.2 1.2])
daspect([1 1 1])
plot(Xe(mi(1:iter,1),1), Xe(mi(1:iter,1),2),'ro')
% title(sprintf('Selected %d points',iter))
disp('Hit enter to proceed')
pause

res=kermat(zeros(1,dx),zeros(1,dx))*ones(Np,1); 
res=res-sum(v'.*v')';
val=zeros(size(xe));
val(ind)=res;
figure
surf(xe,ye,reshape(val,size(xe)))
title('Final Power Function')
shading interp
% figure
% movie(F,1,1)
% movie2avi(F,'powermovie','fps',1)
% figure 
% movie(G,1,1)
% movie2avi(G,'newtonmovie','fps',1)
disp('Hit enter')
pause
figure
val=zeros(size(xe));
val(ind)=sum(abs(v'));
surf(xe,ye,reshape(val,size(xe)))
title('Lebesgue Function')
shading interp
figure
semilogy(Linvnorm);
title('Norm of inverse of L')
figure
semilogy(resnorm);
title('Norm of Power Function Squared')
disp('Hit enter to close figures')
pause

close all;