%% Sparse Grid Demos
%
% Author(s): Paul R. Miles & Ralph C. Smith
%
% Date Created: April 10, 2019
%
% Exampled created using MATLAB2018a.  Certain features may not be
% compatible with previous versions of MATLAB.
%
% The sparsegrid toolbox is available for download on Andreas Klimke's
% ResearchGate page:
% <https://www.researchgate.net/publication/330956911_Sparse_Grid_Interpolation_Toolbox_v511>

% Setup workspace
close all; clear; clc;

% Add path to sparse grid toolbox
% Note: I assume the uncompressed folder has been placed in the current
% working directory.
addpath('spinterp_v5.1.1');

%% Example 1
% $f(x,y) = \sin(x) + \cos(y)$
% Domain: $[0,\pi]\times[0,\pi]$

% define function
f = @(x,y) sin(x) + cos(y);
% construct interpolant
spi = spvals(f, 2, [0 pi; 0 pi]);
% Compute interpolated values - random points in domain
x1 = pi*rand(1, 5); x2 = pi*rand(1, 5);
y = spinterp(spi, x1, x2);
error = max(abs(y - f(x1,x2)));
% display comparison between interpolated and exact
fprintf('(%4s,%4s) -> %6s, %6s\n', 'x1', 'x2', 'y(intrp)', 'f(x1,x2)');
for ii = 1:5
    fprintf('(%4.2f,%4.2f) -> %+6.4f, %+6.4f\n', x1(ii), x2(ii), y(ii), f(x1(ii),x2(ii)));
end
fprintf('max(abs(y-f(x1,x2))) = %4.2e\n', error);
%% Visualize grid
figure(1);
h = plotgrid(spi.maxLevel,spi.d);
pcol = distribute_color_spectrum(length(h), jet);
for ii = 1:length(h)
    set(h(ii), 'Marker','o');
    set(h(ii), 'MarkerSize', 10);
    set(h(ii), 'LineWidth', 2);
    set(h(ii), 'Color', pcol(ii,:));
end
xlabel('X')
ylabel('Y')
set(gca, 'FontSize', 22);
% set_figure_dimensions(8,6)

%% Visualize the interpolant
figure(2)
subplot(1,2,1);
ezmesh(f,[0,pi]);
title('f(x,y)=sin(x) + cos(y)');
subplot(1,2,2);
ezmesh(@(x,y) spinterp(spi,x,y), [0 pi]);
title('Sparse grid interpolant');
set_figure_dimensions(10,6);

%% Compute integral

% using spquad
ns = 100; % perform integration ns times
tic
for ii = 1:ns
    Qsg = spquad(spi);
end
sgtime = toc;
% using MATLAB's integral2
tic
for ii = 1:ns
    Q2 = integral2(f, 0, pi, 0, pi);
end
i2time = toc;
% display comparison
fprintf('Qsg = %6.4f, Q2 = %6.4f\n', Qsg, Q2);
fprintf('Avg. Comp. Time (sec):\n\tQsg -> %4.3e\n\tQ2 -> %4.3e\n', sgtime/ns, i2time/ns);

%% Example 2
% peaks($x, y$)
% Domain: $[0,2]\times[0,2]$

% define function
f = @(x,y) peaks(x,y);
% define sparse grid options
options = spset('DimensionAdaptive', 'on', 'RelTol', 1e-4, 'GridType', 'Chebyshev');
spi = spvals(f, 2, [0,2; 0,2], options);
%% Visualize grid
figure(3);
h = plotgrid(spi.maxLevel,spi.d);
tmpcol = jet;
pcol = distribute_color_spectrum(2*length(h), tmpcol(1:end,:));
for ii = 1:length(h)
    set(h(ii), 'Marker','o');
    set(h(ii), 'MarkerSize', 6);
    set(h(ii), 'LineWidth', 2);
    set(h(ii), 'Color', pcol(2*ii,:));
end
xlabel('X')
ylabel('Y')
cbar=colorbar;
caxis([1, length(h)])
set(cbar, 'YTick', 1:1:length(h))
set(gca, 'FontSize', 22);
set_figure_dimensions(8,6)

%% Compute integral

% using spquad
ns = 100; % perform integration ns times
tic
for ii = 1:ns
    Qsg = spquad(spi);
end
sgtime = toc;
% using MATLAB's integral2
tic
for ii = 1:ns
    Q2 = integral2(f, 0, 2, 0, 2);
end
i2time = toc;
% display comparison
fprintf('Qsg = %6.4f, Q2 = %6.4f\n', Qsg, Q2);
fprintf('Avg. Comp. Time (sec):\n\tQsg -> %4.3e\n\tQ2 -> %4.3e\n',sgtime/ns, i2time/ns);

%% Example 3
% $f = \sin(\pi x)\cdot\sin(\pi y)\cdot2z$
% Domain: $[0,1]\times[0,1]\times[0,1]$

% define function
f = @(x,y,z) sin(pi*x).*sin(pi*y).*2.*z;
% construct interpolant
options = spset('DimensionAdaptive', 'on', 'RelTol', 1e-4, 'GridType', 'Chebyshev');
spi = spvals(f, 3, [0,1; 0,1; 0,1], options);

%% Visualize grid
figure(4);
h = plotgrid(spi.maxLevel,spi.d);
tmpcol = jet;
pcol = distribute_color_spectrum(2*length(h), tmpcol(1:end,:));
for ii = 1:length(h)
    set(h(ii), 'Marker','o');
    set(h(ii), 'MarkerSize', 6);
    set(h(ii), 'LineWidth', 2);
    set(h(ii), 'Color', pcol(2*ii,:));
end
xlabel('X')
ylabel('Y')
zlabel('Z')
cbar=colorbar;
caxis([1, length(h)])
set(cbar, 'YTick', 1:1:length(h))
set(gca, 'FontSize', 22);
set_figure_dimensions(8,6)

%% Compute integral

% using spquad
ns = 100; % perform integration ns times
tic
for ii = 1:ns
    Qsg = spquad(spi);
end
sgtime = toc;

% using MATLAB's integral3
tic
for ii = 1:ns
    Q3 = integral3(f, 0, 1, 0, 1, 0, 1);
end
i3time = toc;
% display comparison
fprintf('Qsg = %6.4f, Q3 = %6.4f\n', Qsg, Q3);
fprintf('Avg. Comp. Time (sec):\n\tQsg -> %4.3e\n\tQ3 -> %4.3e\n',sgtime/ns, i3time/ns);

%% Support Functions
function [psty] = distribute_color_spectrum(nds, cmap)
% evenly distribution color palette across nds data sets
if nargin < 2 || isempty(cmap)
    cmap = jet; % default color map
end
cmp = colormap(cmap);
[mc,~] = size(cmp);
nc = linspace(1,mc,nds)';
psty = zeros(nds,3);
for ii = 1:nds
    psty(ii,:) = interp1((1:1:mc)', cmp, nc(ii));
end
end

function set_figure_dimensions(width, height)
% adjust figure dimensions
set(gcf,'units','inches');
pos = get(gca, 'pos');
set(gcf, 'Position', [pos(1) pos(2) width height]);
end