close all
clear all
%--------------------------------------------------------------------------
% parameters of the experiment
%--------------------------------------------------------------------------

degvalues=22:4:38;
method='weightedbergman'%'kernel';
name='diskprofile238WK2near.fig';
namepdf='diskprofile238WK2near.pdf';
name2='diskcontour238WK2near.fig';
name2pdf='diskcontour238WK2near.pdf';
alpha=2;
xm=-2;
xM=2;
ym=-2;
yM=2;
npts=100;
toll=10^-5;
N=length(degvalues);
AccMeth='row';
AccSeqMeth='originalsequence';
AccSeq=[];
col=3;
cutofferror=10^-5;



%--------------------------------------------------------------------------
% Experiment Body 1: mesh/grid generation, func. eval./reconstruction
%--------------------------------------------------------------------------
% boundary of the regular unit m-agon


% computational grid (evaluation points)
x=linspace(xm,xM,npts);
y=linspace(ym,yM,npts);
[X,Y]=meshgrid(x,y);
pts=[X(:),Y(:)];
% true solution
Vtrue=SiciakRealDisk(X,Y);
fid=fopen('seriesresults.txt','a');
fprintf(fid,'\n');
fprintf(fid,'settype: unit ball\t alpha %f\n',alpha);
Vold=zeros(size(X));
columnindex=1;
tableofV=zeros(npts^2,length(degvalues));
shiftind=1;
for deg=degvalues
    deg
%alpha=max(1,log(deg));
    degwam=round(alpha*deg);
% mesh creation
    [Xw,Yw,C]=CircleAM(degwam,4);
    wam=[Xw(:),Yw(:)];
    %wam = approxfek(degwam,wam,1,wam,[]);
    ind=X.^2+Y.^2-1<=0;
% approximated solution
    V = SZEFw(deg,wam,pts,method);
    V=reshape(V,size(Vtrue));
    m=max(V(ind));
% compute errors and subsitutes values on the set
    ErrL=abs(V-Vtrue);
    ErrL(ind)=cutofferror;
    Step(shiftind)=sum(sum(abs(V-Vold)))/(npts^2);
    AL1(shiftind)=sum(sum(ErrL))/(npts^2);
    RL1(shiftind)=sum(sum(ErrL))/sum(sum(V));
    AM(shiftind)=max(max(ErrL));
% prints results on file    
    fprintf(fid,...
        'deg %2d\t wamdeg %d\t absolute L1 error %f\t relative L1 error %f\t absolute max error %f\n',...
        deg,degwam,AL1(shiftind),RL1(shiftind),AM(shiftind))
    Vold=V;
    tableofV(:,columnindex)=V(:);
    columnindex=columnindex+1;
    shiftind=shiftind+1;
end
%%
% accelerated algorithm and computations of errors
[tableofV1,TI]=scaledvectorrho(tableofV,degvalues,AccSeqMeth,toll,AccSeq);

switch AccMeth
    case 'diagonal'
degvaluesAcc=diag(TI);
degvaluesAcc=degvaluesAcc(degvaluesAcc>0);
jmax=length(degvaluesAcc);
for j=1:jmax
    V1temp=reshape(tableofV1(j,:,j)',size(Vtrue));
% compute errors and subsitutes values on the set
    ErrLAcc=abs(V1temp-Vtrue);
    ErrLAcc(ind)=cutofferror;
    AL1Acc(j)=sum(sum(ErrLAcc))/(npts^2);
    RL1Acc(j)=sum(sum(ErrLAcc))/sum(sum(V1temp));
    AMAcc(j)=max(max(ErrLAcc));
end
    case 'column'
degvaluesAcc=TI(:,col);
degvaluesAcc=degvaluesAcc(degvaluesAcc>0);
jmax=length(degvaluesAcc);
for j=1:jmax
    V1temp=reshape(tableofV1(j,:,col),size(Vtrue));
% compute errors and subsitutes values on the set
    ErrLAcc=abs(V1temp-Vtrue);
    ErrLAcc(ind)=cutofferror;
    AL1Acc(j)=sum(sum(ErrLAcc))/(npts^2);
    RL1Acc(j)=sum(sum(ErrLAcc))/sum(sum(V1temp));
    AMAcc(j)=max(max(ErrLAcc));
end
case 'row'
degvaluesAcc=TI(col,:);
degvaluesAcc=degvaluesAcc(degvaluesAcc>0);
jmax=length(degvaluesAcc);
for j=1:jmax
    V1temp=reshape(tableofV1(col,:,j),size(Vtrue));
% compute errors and subsitutes values on the set
    ErrLAcc=abs(V1temp-Vtrue);
    ErrLAcc(ind)=cutofferror;
    AL1Acc(j)=sum(sum(ErrLAcc))/(npts^2);
    RL1Acc(j)=sum(sum(ErrLAcc))/sum(sum(V1temp));
    AMAcc(j)=max(max(ErrLAcc));
end
end
%%
% save('AL1','AL1')
% save('RL1','RL1')
% save('AM','AM')
% save('Vold','Vold')
% save('V1','V1');

%%
fclose(fid)
%%
logal1=log10(AL1);
logdegvalues=log10(degvalues);
cal1=polyfit(logdegvalues,logal1,1);
AL1fit=10.^(polyval(cal1,logdegvalues));

logrl1=log10(RL1);
logdegvalues=log10(degvalues);
crl1=polyfit(logdegvalues,logrl1,1);
RL1fit=10.^(polyval(crl1,logdegvalues));

logam=log10(AM);
logdegvalues=log10(degvalues);
cam=polyfit(logdegvalues,logam,1);
AMfit=10.^(polyval(cam,logdegvalues));
f1=figure(1)
loglog(degvalues,AL1,'.')
hold on;
loglog(degvalues,AL1fit)
loglog(degvalues,RL1,'g.')
loglog(degvalues,RL1fit,'g')
loglog(degvalues,AM,'k.')
loglog(degvalues,AMfit,'k')
loglog(degvaluesAcc,AL1Acc,'r')
loglog(degvaluesAcc,RL1Acc,'-.r')
loglog(degvaluesAcc,AMAcc,'-or')
loglog(degvalues,log(degvalues)./degvalues,'-.m')
legend('absolute L^1 error','lin fitting','relative L-1 error','lin fitting'...
    ,'absolute max error','lin fitting','Acc. absolute L-1 error',...
    'Acc. relative L^1 error','Acc. absolute max error','(log k)/k',...
    'Location','SouthWest');
saveas(gcf,name)
saveas(f1,namepdf)
f2=figure(2);
subplot(1,2,1),contourf(X,Y,log10(ErrL))
axis equal
colorbar('South');
subplot(1,2,2),contourf(X,Y,log10(ErrLAcc))
axis equal
colorbar('South');

saveas(gcf,name2)
print(f2,name2pdf)

% clear degvaluesAcc
% clear AL1Acc
% clear RL1Acc
% clear AMAcc