function ezspline_time()
%EZSPLINE_TIME Compares toolbox runtimes (Ezspline vs. MATLAB Spline).
%
%   MATLAB Spline Toolbox required.
%

time2d;

time3d;

return

%
%-----------------------------------------------------------------
%
function time2d()

%
% benchmark table : columns (time, error)
%
matlabte = [];
ezte = [];
gridsizes = 10 : 10 : 100;

%
% define function & domain
%
nx = gridsizes(1);
xmin = 0; xmax = 1;
ny = gridsizes(1);
ymin = 0; ymax = 1;

x = linspace(xmin, xmax, nx);
y = linspace(ymin, ymax, ny);

[xx, yy] = ndgrid(x,y);
f = xx.^3 + 2 * yy.^3;

%
% matlab spline object
%
spl = csapi({x,y}, f);

% 
% toolbox spline object
%
[ez, ier] = ezspline_setup(f, [[0 0], [0 0]]);

if ( ier ~= 0 )
  error('EZSPLINE_TIME Error: Ezspline object failed to initialize.');
end

%
% loop over different grid sizes
% grid interp
%
for grid = 1:length(gridsizes)
        ni = gridsizes(grid);
        i = linspace(xmin, xmax, ni);

        [xi, yi] = ndgrid(i,i);
        fexact =  xi.^3 + 2 * yi.^3;

        % with MATLAB
        tic;
        fi = fnval(spl, {i, i});
        matlabte(grid, 1) = toc;
        matlabte(grid, 2) = err(fi,fexact);
        
        % with ezspline toolbox
        tic;
        [fi, ier] = ezspline_interp(ez, i, i);
        if (ier ~= 0)
           disp('ERROR: ezspline_interp failed on 2D grid.');
        end
        ezte(grid, 1) = toc;
        ezte(grid, 2) = err(fi,fexact);
        
end

%
% interp a series of points
%
cloudsizes = 10 : 10 : 50;

for grid = 1:length(cloudsizes)
        ni = cloudsizes(grid);
        i = linspace(xmin, xmax, ni);

        [xi, yi] = ndgrid(i,i);
        fexact =  xi.^3 + 2 * yi.^3;
        P = zeros(length(xi(:)),2);
        P(:,1) = xi(:);
        P(:,2) = yi(:);

        % with MATLAB
        fi = zeros(1,size(P,1));
        tic;
        for p = 1:size(P,1)
             fi(p) = fnval(spl, {P(p,1), P(p,2)});
        end
        matlabte(grid, 3) = toc;
        matlabte(grid, 4) = err(fi,fexact);
        
        % with ezspline toolbox
        tic;
        [fi, ier] = ezspline_interp(ez, P);
        ezte(grid, 3) = toc;
        ezte(grid, 4) = err(fi(:,1),fexact);
        if (ier ~= 0)
           disp('ERROR: ezspline_interp failed on 2D cloud of points.');
        end

end

plot2dresults(matlabte, ezte, gridsizes, cloudsizes.^2);

return;

%
%-----------------------------------------------------------
%
function error = err(fi,fexact)

error =  (fi(:)-fexact(:)).^2 ./ prod(size(fi)) ;
error = sum( error(:) );

return

%
%-----------------------------------------------------------
%
function plot2dresults(mte,ezte,grids,clouds)

%
% gridded interp
%
subplot(2,2,1); 
hold on;

% matlab 
[ax,h1,h2] = plotyy(grids,mte(:,1),grids,mte(:,2));
% ez
[ax,h3,h4] = plotyy(grids,ezte(:,1),grids,ezte(:,2));
% line colors
set(h3,'LineStyle', '--');
set(h4,'LineStyle', '--');
% adjust axis labels
% time
t = [min(min(ezte(:,1),mte(:,1))) max(max(ezte(:,1),mte(:,1)))];
set(ax(1), 'ylim', t);
set(ax(1), 'ytick', linspace(t(1),t(2),5) );

% error
e = [ min(min(ezte(:,2),mte(:,2))) max(max(ezte(:,2),mte(:,2))) ];
set(ax(2), 'ylim', e);
set(ax(2), 'ytick', linspace(e(1),e(2),5) );

% labels
legend([h1 h3], 'MATLAB', 'Ezspline', 2);
set(get(ax(1), 'YLabel'), 'String', 'Runtime (s)');
set(get(ax(2), 'YLabel'), 'String', 'Error');
xlabel('Grid Size (N x N)');
title('2D Grid Interpolation');

%
% cloud interp
%
subplot(2,2,2); 
hold on;

% matlab 
nc = length(clouds);
[ax,h1,h2] = plotyy(clouds,mte(1:nc,3),clouds,mte(1:nc,4));
% ez
[ax,h3,h4] = plotyy(clouds,ezte(1:nc,3),clouds,ezte(1:nc,4));
% line colors
set(h3,'LineStyle', '--');
set(h4,'LineStyle', '--');
% adjust axis labels
% time
t = [min(min(ezte(1:nc,3),mte(1:nc,3))) max(max(ezte(1:nc,3),mte(1:nc,3)))];
set(ax(1), 'ylim', t);
set(ax(1), 'ytick', linspace(t(1),t(2),5) );

% error
e = [ min(min(ezte(1:nc,4),mte(1:nc,4))) max(max(ezte(1:nc,4),mte(1:nc,4))) ];
set(ax(2), 'ylim', e);
set(ax(2), 'ytick', linspace(e(1),e(2),5) );

% labels
legend([h1 h3], 'MATLAB', 'Ezspline', 2);
set(get(ax(1), 'YLabel'), 'String', 'Runtime (s)');
set(get(ax(2), 'YLabel'), 'String', 'Error');
xlabel('Cloud Size (# points)');
title('2D Point Cloud Interpolation');

return

%
%-----------------------------------------------------------------
%
function time3d()

%
% benchmark table : columns (time, error)
%
matlabte = [];
ezte = [];
gridsizes = 10 : 10 : 100;

%
% define function & domain
%
nx = gridsizes(1);
xmin = 0; xmax = 1;
ny = gridsizes(1);
ymin = 0; ymax = 1;
nz = gridsizes(1);
zmin = 0; zmax = 1;

x = linspace(xmin, xmax, nx);
y = linspace(ymin, ymax, ny);
z = linspace(zmin, zmax, nz);

[xx, yy, zz] = ndgrid(x,y,z);
f = xx.^3 + 2 * yy.^3 + 3 * zz.^3;

%
% matlab spline object
%
spl = csapi({x,y,z}, f);

% 
% toolbox spline object
%
[ez, ier] = ezspline_setup(f, [[0 0], [0 0], [0 0]]);

if ( ier ~= 0 )
  error('EZSPLINE_TIME Error: Ezspline object failed to initialize.');
end

%
% loop over different grid sizes
% grid interp
%
for grid = 1:length(gridsizes)
        ni = gridsizes(grid);
        i = linspace(xmin, xmax, ni);

        [xi, yi, zi] = ndgrid(i,i,i);
        fexact =  xi.^3 + 2 * yi.^3 + 3 * zi.^3;

        % with MATLAB
        tic;
        fi = fnval(spl, {i, i, i});
        matlabte(grid, 1) = toc;
        matlabte(grid, 2) = err(fi,fexact);
        
        % with ezspline toolbox
        tic;
        [fi, ier] = ezspline_interp(ez, i, i, i);
        if (ier ~= 0)
           disp('ERROR: ezspline_interp failed on 3D grid.');
        end

        ezte(grid, 1) = toc;
        ezte(grid, 2) = err(fi,fexact);
        
end

%
% interp a series of points
%
cloudsizes = 10 : 10 : 50;

for grid = 1:length(cloudsizes)
        ni = cloudsizes(grid);
        i = linspace(xmin, xmax, ni);

        [xi, yi, zi] = ndgrid(i,i,i);
        fexact =  xi.^3 + 2 * yi.^3 + 3 * zi.^3;
        P = zeros(length(xi(:)),3);
        P(:,1) = xi(:);
        P(:,2) = yi(:);
        P(:,3) = zi(:);
        
        % with MATLAB
        fi = zeros(1,size(P,1));
        tic;
        for p = 1:size(P,1)
             fi(p) = fnval(spl, {P(p,1), P(p,2), P(p,3)});
        end
        matlabte(grid, 3) = toc;
        matlabte(grid, 4) = err(fi,fexact);
        
        % with ezspline toolbox
        tic;
        [fi, ier] = ezspline_interp(ez, P);
        if (ier ~= 0)
           disp('ERROR: ezspline_interp failed on 3D cloud of points.');
        end
        
        ezte(grid, 3) = toc;
        ezte(grid, 4) = err(fi(:,1),fexact);

end

plot3dresults(matlabte, ezte, gridsizes, cloudsizes.^3);

return;

%
%-----------------------------------------------------------
%
function plot3dresults(mte,ezte,grids,clouds)

%
% gridded interp
%
subplot(2,2,3); 
hold on;

% matlab 
[ax,h1,h2] = plotyy(grids,mte(:,1),grids,mte(:,2));
% ez
[ax,h3,h4] = plotyy(grids,ezte(:,1),grids,ezte(:,2));
% line colors
set(h3,'LineStyle', '--');
set(h4,'LineStyle', '--');
% adjust axis labels
% time
t = [min(min(ezte(:,1),mte(:,1))) max(max(ezte(:,1),mte(:,1)))];
set(ax(1), 'ylim', t);
set(ax(1), 'ytick', linspace(t(1),t(2),5) );

% error
e = [ min(min(ezte(:,2),mte(:,2))) max(max(ezte(:,2),mte(:,2))) ];
set(ax(2), 'ylim', e);
set(ax(2), 'ytick', linspace(e(1),e(2),5) );

% labels
legend([h1 h3], 'MATLAB', 'Ezspline', 2);
set(get(ax(1), 'YLabel'), 'String', 'Runtime (s)');
set(get(ax(2), 'YLabel'), 'String', 'Error');
xlabel('Grid Size (N x N x N)');
title('3D Grid Interpolation');

%
% cloud interp
%
subplot(2,2,4); 
hold on;

% matlab 
nc = length(clouds);
[ax,h1,h2] = plotyy(clouds,mte(1:nc,3),clouds,mte(1:nc,4));
% ez
[ax,h3,h4] = plotyy(clouds,ezte(1:nc,3),clouds,ezte(1:nc,4));
% line colors
set(h3,'LineStyle', '--');
set(h4,'LineStyle', '--');
% adjust axis labels
% time
t = [min(min(ezte(1:nc,3),mte(1:nc,3))) max(max(ezte(1:nc,3),mte(1:nc,3)))];
set(ax(1), 'ylim', t);
set(ax(1), 'ytick', linspace(t(1),t(2),5) );

% error
e = [ min(min(ezte(1:nc,4),mte(1:nc,4))) max(max(ezte(1:nc,4),mte(1:nc,4))) ];
set(ax(2), 'ylim', e);
set(ax(2), 'ytick', linspace(e(1),e(2),5) );

% labels
legend([h1 h3], 'MATLAB', 'Ezspline', 2);
set(get(ax(1), 'YLabel'), 'String', 'Runtime (s)');
set(get(ax(2), 'YLabel'), 'String', 'Error');
xlabel('Cloud Size (# points)');
title('3D Point Cloud Interpolation');

return