function [x_f, z_f, mse_f, max_error] = find_best_table_1d(x_0, z_0, ...
% This function finds the best way to reduce the size of a table to the
% specified number of points. That is, suppose you had a table of 10,000
% points, but wanted to find the closest 10-point table. This function
% would find what 10 points best represent the 10,000 point table.
% This is accomplished using traditional constrained optimization
% techniques to best fit a smaller table with the specified dimensions.
% "Goodness" will be measured by the difference between the original table
% and the new table sampled at the original data points. See
% find_best_table_demo.m for examples.
% Note: Using the find_best_table_nd function is actually faster. This
% implementation remains because it's easier to read for those interested
% in learning about how the algorithm works and for backwards compatibility
% with the original release.
% x_0 Array of values for the first independent variable
% z_0 Matrix of values for the dependent variable for all x_0
% n_x Number of divisions of the first independent variable in output
% method Interpolation method to use, e.g., 'linear' (see help interp1)
% x_f Best division of first independent variable found
% z_f Resulting fit of z_0 on x_f
% mse_f Mean square error of final table
% max_error Maximum error anywhere
% This builds off of work by Richard Willey, Stuart Kozola, Eric Johnson,
% Sarah Zaranek, and Peter Maloney at The MathWorks, Inc.
% Requires the Optimization Toolbox (TM).
% Supports the Parallel Computing Toolbox (TM).
% Tucker McClure @ The MathWorks
% Copyright 2013 The MathWorks, Inc.
% Convert from the spacing between the points to the grid points.
function [x_i] = spacing_to_x(spacing_i)
x_i = [[0; cumsum(spacing_i)] + x_0(1); x_0(end)];
% Determine how well the current spacing fits the original data.
function [mse, max_error, z_i] = evaluate_spacing(spacing_i)
% Convert to grid points.
x_i = spacing_to_x(spacing_i);
% Evaluate it and find the difference between the new mesh and z_0,
% evaluated at x_0. Get the mean square error from these.
z_i = interp1(x_0, z_0, x_i, method);
se = (z_0 - interp1(x_i, z_i, x_0, method)).^2;
mse = mean(se(:));
max_error = sqrt(max(se(:)));
% Do some suick error checking and set defaults for arguments.
if nargin == 3
method = 'linear';
elseif nargin > 4 || nargin < 3
error('Incorrect number of arguments.');
% Our state will be the spacing between the xs.
dx_i = (x_0(end) - x_0(1))/(n_x - 1);
spacing_0 = dx_i * ones(1, n_x - 2)';
% Set up constraints. We want all the xs and ys to add up to the table
% spans, and we can't let any dx be <= 0.
% Set the lower bound to be very small, 1e-3 * the difference in
% uniformly spaced x values.
lb = 1e-3 * dx_i * ones(1, n_x - 2);
ub = ;
% There are no inequality constraints.
A = ones(1, n_x - 2);
b = x_0(end) - x_0(1) - lb(1);
% Total of distances between points needs to equal the span.
Aeq = ;
beq = ;
% Call an optimization routine.
spacing_f = fmincon(@evaluate_spacing, spacing_0, ...
A, b, Aeq, beq, lb, ub, , ...
optimset('Algorithm', 'active-set', ...
'MaxFunEvals', 10000, ...
'UseParallel', 'always', ...
% Convert back to grid points, and we're done.
x_f = spacing_to_x(spacing_f);
% Output the final table.
[mse_f, max_error, z_f] = evaluate_spacing(spacing_f);