clear;
N=200;
LP=60;
LT=40;

ntry=1;%number of random trials. We take the average of all trials to estimate the error of regression

freq=4;%frequency of the function to regress (f)

Err1=[];Err2=[];
for i=1:ntry,
    disp(['Random trial number ',num2str(i)])
    
    phi=2*pi*rand(N,1);
    theta=rand(N,1);theta=(theta-min(theta))/(max(theta)-min(theta));theta=theta*pi-pi/2;
    x=cos(phi).*cos(theta);
    y=sin(phi).*cos(theta);
    z=sin(theta);
    G=[x y z];%The data set is generated: points scattered at random on the sphere
    
    f=sin(freq*theta);%The function to regress
    
    tP=(0:LP-1)'/(LP-1);
    tT=(0:LT-1)'/(LT-1);
    Ephi=2*pi*tP*ones(size(tT'));Ephi=Ephi';Ephi=Ephi(:);
    Etheta=pi*ones(size(tP))*tT'-pi/2;Etheta=Etheta';Etheta=Etheta(:);
    Ex=cos(Ephi).*cos(Etheta);
    Ey=sin(Ephi).*cos(Etheta);
    Ez=sin(Etheta);
    Ex=cos(Ephi).*cos(Etheta);
    Ey=sin(Ephi).*cos(Etheta);
    Ez=sin(Etheta);
    H=[Ex,Ey,Ez];%The fine grid on which one wants to predict the values of f
    
    D=1;
    
    [E,V,eigen]=extension_prepare(G,H,'Gaussian',D,100);%Gaussian geometric harmonic kernel, width=1
    c=V*f;
    Ef=E*c;%f is extended as Ef on the fine grid H. This direct interpolation is inadapted to f, since the extension kernel is wide, and f has high bandwidth
    
    pat=sin(freq*Etheta);%The ideal (true value) of f on H.
    Err1=[Err1 norm(Ef-pat)/norm(pat)];%Relative error of regression
        
    g=sin(theta);%We represent f=h(g) with g=sin(theta). g has low bandwidth
    c=V*g;
    Eg=E*c;
    Eg=(Eg-min(Eg))/(max(Eg)-min(Eg))*2-1;%g is extended and then rescale otherwise splines are forced to extrapolate -->go to hell (see below)
        
    h=spline(g,f,Eg);%h is interpolated using splines.
    Err2=[Err2 norm(h-pat)/norm(pat)];%relative error with this procedure
end
st=2;%subsampling factor for the plots as scatter3 is slow on my machine
submat=zeros(LT,LP);submat(1:st:LT,1:st:LP)=1;
submat=submat(:);sub=find(submat==1);

figure(1);scatter3(x,y,z,15,f,'filled');title('Function f','Fontsize',14);
figure(2);scatter3(Ex(sub),Ey(sub),Ez(sub),15,h(sub),'filled');title('Extension of f','Fontsize',14);
%figure(3);plot(Eg,h,'b.');hold on;plot(g,f,'r.');hold off;
figure(3);plot(g,f,'*');title('g vs h','Fontsize',14);

Err1=mean(Err1);
Err2=mean(Err2);
disp(strcat('Relative error for direct interpolation=',num2str(Err1)))
disp(strcat('Relative error for adapted interpolation=',num2str(Err2)))