find_best_table_nde(...

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

%% find_best_table_nde
%
% 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 100-by-100 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 
% along each axis to accomplish this and where those 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.
% 
% Inputs:
%   x_0    Cell array of values for the first independent variables (e.g., 
%          {x1, x2, ... xn}.
%   z_0    Matrix of values for the dependent variable for all x_0 and
%          y_0 using ndgrid-style syntax (e.g., z_0(2, 3) corresponds
%          to x{1}(2) and x{2}(3)). Note that this is different from what
%          is produced by meshgrid.
%   mse_target Target mean squared error
%   method     Interpolation method to use, e.g., 'linear' (help interp2)
% 
% Outputs:
%   x_f    Cell array of best divisions of independent variables found
%   z_f    Resulting fit of z_0 on x_f
%   mse_f  Mean square error of 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.
           
    % Default to linear interpolation.
    if nargin < 5, method = 'linear'; end;
    if nargin == 4 && ischar(max_error)
        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;

    % Number of dimensions
    n_d = length(x_0);
    
    % Initialize the algorithm with the simplest possible table.
    n = 3 * ones(1, n_d);

    % See how well simplest possible table does.
    [~, ~, mse, me] = find_best_table_nd(x_0, z_0, n, method);

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

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

        % Try increasing the number of points in each dimension by 2 
        % separately.
        mse_i = zeros(1, n_d);
        me_i  = zeros(1, n_d);
        for k = 1:n_d
            n_i = n;
            n_i(k) = 2*n(k);
            [~, ~, mse_i(k), me_i(k)] = find_best_table_nd(x_0, z_0, ...
                                                           n_i, method);
        end
        
        % Find which made the most improvement and use it alone.
        [~, best_index] = min(mse_i);
        n(best_index)   = 2*n(best_index);
        mse             = mse_i(best_index);
        me              = me_i(best_index);
        
    end

    % Ok, now we now that (n_x/2, n_y) and (n_x, n_y/2) aren't good enough,
    % but (n_x, n_y) is. Search between them.
    n_left  = max(floor(n/2), 3);
    n_right = n;

    % Try to bring in any of the right bounds if it's possible, starting
    % with what's currently the smallest number of indices. If it's not
    % possible to move in any of the right bounds, move in the left bounds.
    count = 0;
    while    any(n_right > n_left + 1) ...
          && count < max_iterations

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

        % Sort the indices; we want to decrease the number of breakpoints
        % for the biggest dimension first.
        [~, sort_indices] = sort(n_right, 2, 'descend');

        % We're going to look for the first dimension we can update. If we
        % can't update anything, we'll allow multiple dimensions to grow.
        updated = false;
        for k = sort_indices
            
            % If the bounds are already pressed together, skip to the next
            % dimension.
            if n_right(k) <= n_left(k) + 1
                continue;
            end

            % Update this dimension as the midpoint between left and right
            % bounds.
            n_t = n_right;
            n_t(k) = round(0.5*(n_right(k) - n_left(k))) + n_left(k);
            
            % See how well it does.
            [~,~, mse_i, me_i] = find_best_table_nd(x_0, z_0, n_t, method);
            
            % If we're still good, move the right bound in.
            if mse_i < mse_target && me_i < max_error
                n_right = n_t;
                updated = true;
                break;
            end

        end
        
        % Otherwise, we can't update any dimension, so move in the left 
        % bounds.
        if ~updated
            n_left = n_t;
        end
        
    end

    % Take the (now just good enough) right bounds.
    n = n_right;

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

end