% TestPDFpoly  Test of the function PDFpoly and PolyInterPol 
%              Most functions can be approximated by piecewise polynomials
% Here a 1D function y=f(x) is approximated by polynomials
% the function is given in different ways in PDFpoly and PolyInterPol
% see help text for piecewise polynomials in Matlab, ex: help mkpp,
% and PDFpoly and PolyInterPol

%----------------------------------------------------------------------
% Copyright (c) 2002.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  karl.skretting@tn.his.no   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  18.12.2002  KS: function made
%----------------------------------------------------------------------

clear all;
Mfile='TestPDFpoly';

TestNo=4;   
% 1: test only PolyInterPol
% 2: test PDFpoly (which uses PolyInterPol)
% 3: find the "error distribution" property of a frame
% 4: find the "error distribution" property of a simple frame
            
if TestNo==1
    % approximate the function: y=e^(-x*x)
    dx=0.005;
    x=0:dx:4;
    L=length(x);
    y=exp(-x.*x);
    x=x(:);y=y(:);
    bs=[0,0.5,1,1.5,2,3,4];
    bs=bs(:);
    deg=[4,4,4,4,4,3];
    cont=[1,1,1,1,1];
    % fixval=[0,1,0;0,0,1;0.75,0.6,0];   % here an error is demaded for x=0.75
    % fixval=[0,1,0;0,0,1;0.75,exp(-0.5625),0];  % here it should be correct at x=0.75
    fixval=[0,1,0;0,0,1];   % and here only special demands for x=0
    %
    pp=PolyInterPol(x,y,bs,deg,cont,fixval);  % here the piecewise polynomial is made
    figure(1);clf;
    fnplt(pp);
    hold on;
    plot(x,y,'r-');
    grid on;
    %
    f=fnval(pp,x);
    figure(2);clf;
    plot(x,y-f);
    grid on;
    
    disp([Mfile,': 1-norm of difference is ',num2str(norm(y-f,1))]);
    disp([Mfile,': 2-norm of difference is ',num2str(norm(y-f,2))]);
    
    ppi=fnint(pp);
    disp([Mfile,': integral of f (the ppi)            is ',num2str(fnval(ppi,4),8)]);
    disp([Mfile,': trapesoid integral of f (pp)       is ',num2str((sum(f)+sum(f(2:(end-1))))*dx/2,8)]);
    disp([Mfile,': trapesoid integral of y=exp(-x.*x) is ',num2str((sum(y)+sum(y(2:(end-1))))*dx/2,8)]);
    disp([Mfile,': sqrt(pi)/2  (the true value)       is ',num2str(sqrt(pi)/2,8)]);
    % for a special value
    x1=2.0;yest=ppval(pp,x1);y1=exp(-x1*x1);e=y1-yest;disp([x1,y1,yest,e*10000]);
end

if TestNo==2
    % the normal error distribution function:  y2=(1/sqrt(2*pi))*exp(-0.5*(x.^2))
    % is approximated from some random values (from this distribution)
    data=randn(3000,1);
    bs=-4:1:4;
    deg=[1,2,2,2,2,2,2,1];
    cont=[1,2,2,2,2,2,1];
    fixval=[0,0,1; -4,0,0; 4,0,0];  % at x=0 y'=0, at x=-4 y=0, and at x=4 y=0
    N=length(bs)-1;
    L=200;
    %    figure(1);clf;  % PDFpoly also makes a figure (if Display)
    pdffun=PDFpoly(data,L,bs,deg,cont,fixval);
    %
    delta=(bs(N+1)-bs(1))/(L+1);
    x=linspace(bs(1)+delta/2,bs(N+1)-delta/2,L)';   % a column vector
    y2=(1/sqrt(2*pi))*exp(-0.5*(x.^2));  
    plot(x,y2,'r-');
    y3=fnval(pdffun,x);    title('green is histogram based line, blue is polynomial estimate, red is true pdf');
    figure(2);clf;
    plot(x,y2-y3);
    disp(['Average error: ',num2str(sum(abs(y2-y3))/L,8)]);
end

if TestNo==3
    % the error distribution function for a frame is estimated
    [r2,r2max,r2avg,r2mse]=FrameProperties('FrameEx1s20',[14,15,16,17]);
    bs=[0.0,0.4,0.5,0.6,0.8,1.0];
    deg=[3,2,2,2,3];
    cont=[1,1,1,1];
    fixval=[0,0,1; 0,0,0; 1,0,0];  % at x=0 y'=0, at x=0 y=0, and at x=1 y=0
    N=length(bs)-1;
    L=100;
    %
    figure(1);clf;  % PDFpoly also makes a figure (if Display)
    pdfr2=PDFpoly(r2,L,bs,deg,cont,fixval);
end

if TestNo==4
    % the error distribution function for a small frame is estimated
    % the frame is the points in an icosahedron, and the error is the
    % distance to the closest great circle for any pair.
    % the frame c) in Figure 7.4 in my thesis
    c=sqrt(0.5+0.1*sqrt(5));
    s=sqrt(0.5-0.1*sqrt(5));
    F=[c,c,0,0,s,-s; s,-s,c,c,0,0; 0,0,s,-s,c,c];
    [r2,r2max,r2avg,r2mse]=FrameProperties(F,[14,15,16,17]);
    bs=[0,0.05,0.09,0.12,0.14];
    deg=[3,3,3,2];
    cont=[2,2,2];
    fixval=[0.14,0,0];  % at x=0.14 y=0
    N=length(bs)-1;
    L=200;
    pdfr2=PDFpoly(r2,L,bs,deg,cont,fixval);
    title('pdf for r2');
    pause;
    %
    [r1,r1max,r1avg,r1mse]=FrameProperties(F,[10,11,12,13]);
    r1maxt=sqrt(2/3-2/15*sqrt(5)); % the teoretic value
    bs1=[0,0.525,0.55,0.61];      % pdf seems to have maximum at approx. x=0.525
    deg1=[3,3,3];
    cont1=[1,2];
    fixval1=[0,0,0;0.61,0,0];  % at x=0 y=0 and x=0.61 y=0
    N1=length(bs1)-1;
    L1=200;
    %
    pdfr1=PDFpoly(r1,L1,bs1,deg1,cont1,fixval1);
    title('pdf for r1');
    %
end