find_best_table_1de(...

function [x_f, z_f, mse_f, me_f] = find_best_table_1de(...
                                      x_0, z_0, ...
                                      mse_target, max_error, ...
                                      method)

%% find_best_table_1de
%
% This function finds the best way to reduce the size of a table to the
% so that the smaller table matches the original within specified error 
% bounds. That is, suppose you had a table of 10,000 points, but wanted to 
% find the smallest table that had no more than 1.5 units of mean squared 
% error. This function would find how many points are necessary to
% accomplish this and where the points go. "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_nde 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
%   mse_target Target mean squared error
%   max_error  Maximum allowable error anywhere
%   method     Interpolation method to use, e.g., 'linear' (see 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
%   me_f   Maximum error anywhere in final table
%
% 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.
       
    % Set some defaults.
    if nargin < 5, method = 'linear'; end;
    if nargin == 4 && ischar(max_error) % fcn(x_0, z_0, mse, method)
        method = max_error;
        max_error = inf;
    end;
    if nargin == 3 || isempty(max_error), max_error = inf; end;
    if isempty(mse_target),               mse_target = inf; end;
                                  
    % Initialize the algorithm with the simplest possible table.
    n_x = 3;
    [~, ~, mse, me] = find_best_table_1d(x_0, z_0, n_x, method);

    % Loop until we meet tolerance or have tried too many times.
    count = 0;
    max_iterations = length(x_0);
    while (mse > mse_target || me > max_error) && count < max_iterations

        % Keep count of number of iterations.
        count = count + 1;

        % Try increasing x and y independently.
        n_x = 2*n_x;
        [~, ~, mse, me] = find_best_table_1d(x_0, z_0, n_x,method);

    end

    % Ok, now we now that n_x/2 isn't good enough, but n_x is. Search 
    % between them.
    n_x_left  = n_x/2;
    n_x_right = n_x;

    % Update left and right bounds until we've corned the solution.
    while    n_x_right > n_x_left + 1 ...
          && count < max_iterations

        % Keep count of number of iterations.
        count = count + 1;

        % If we can update x, try it.
        n_x_t = round(0.5*(n_x_right - n_x_left)) + n_x_left;
        [~, ~, mse, me] = find_best_table_1d(x_0, z_0, n_x_t, method);

        % If we can update x, do so.
        if mse < mse_target && me < max_error
            n_x_right = n_x_t;
        % Otherwise, we can't update x, so move in the left bound.
        else
            n_x_left = n_x_t;
        end

    end

    % Take the (now just good enough) right bound.
    n_x = n_x_right;

    % Gather the final output.
    [x_f, z_f, mse_f, me_f] = find_best_table_1d(x_0, z_0, n_x, method);

end