find_best_table_1d(x_0, z_0, ...

function [x_f, z_f, mse_f, max_error] = find_best_table_1d(x_0, z_0, ...
                                                           n_x, ...
                                                           method)

%% find_best_table_1d
%
% 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.
% 
% Inputs:
%   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)
% 
% Outputs:
%   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)];
    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(:)));
        
    end

    % 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.');
    end
    
    % 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', ...
                                 'Display',     'notify'));

    % 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);
    
end