function [ boundary_envelopes ] = find_boundary( A );
%
% FIND_BOUNDARY  Find the support boundary for a particular matrix, A.
% 
% This routine finds the support boundary for a binary matrix, A, by letting the 
% support threshold, m, i.e., the number of non-zero rows that a column must contain, 
% vary, from 1 to the number of rows in A, and for each such support threshold, m,
% finding the maximum column threshold, q, i.e., the maximum number of non-zero 
% rows that a column can contain for the given row threshold.  
%
% [ boundary_envelopes ] = find_boundary( A );

% Copyright © Regents of the University of Minnesota. All rights reserved.

% Loop over all support thresholds, from 1 to the number of rows.
%   To improve efficiency we use several strategies. First, our basic 
%   underlying strategy is to look for the boundary row by row and, 
%   within each row, we perform a binary search to find the boundary.
%   Second, when we find an envelope, while doing a binary search in
%   a row, we use the actual m value of this envelope, not the specified
%   values. Since the actual m value will be at least as large this 
%   narrows the interval of search. Third, if the boundary envelope 
%   that we find 'covers' additional rows, i.e., its value of m is greater 
%   than that of the current row index, then we can skip those rows.
%   Fourth, we can bound the interval for the binary search because the 
%   maximum n value in a row is at most the n value of the previous envelope
%   minus 1. Fifth, we can also bound the interval for the binary search
%   by estimating a bound for the minimum value of n. We do this by recording
%   the n values of the last few envelopes and taking the maximum of 
%   the n values of successive envelopes. Sixth, when performing the binary
%   search, we can sometimes save one search for an envelope if realize that
%   the last good envelope.
% 
% Peform some initialization.
[ M N ] = size( A );
max_n = N;
num_envelopes = 0;
m = 1;
done = 0;

num_last_envs = 3;
estimate_failed = 0;
diffs = [];

% Loop
while ~done
    
    % If we have seen a few envelopes, try to estimate the minimum bound on n.
    if num_envelopes > num_last_envs & ~estimate_failed
        min_n = boundary_envelopes( num_envelopes ).n - max( diffs );
        min_n = max( 1, min_n );
        original_min_n = min_n;
        min_n_is_estimate = 1;
        n = min_n;
    % Otherwise, just start with a minimum bound of 1 for n.
    else
        min_n = 1;
        min_n_is_estimate = 0;
        if estimate_failed
            estimate_failed = 0;
        end
        n = ceil( ( min_n + max_n ) / 2 );
    end
    found_max_n = 0;
    last_good_envelope = [];
    
    % Do a binary search to find the non-empty envelope with the max value of n.
    while ~found_max_n
        % Find the envelope for the given value of m and n.
        env = find_envelope( A, m, n, last_good_envelope );
        %disp( n );
        
        % Set up for the next iteration
        if isempty( env )  % We have gone too far - go back.
            max_n = n - 1;
            if max_n < 1
                done = 1;
                break;
            end
            if min_n_is_estimate
                estimate_failed = 1;
                %disp('estimate failed');
                break;
            end
        else % Perhaps not far enough
            min_n = env.n;
            last_good_envelope = env;
        end
        old_n = n;
        n = ceil( ( min_n + max_n ) / 2 );
    
        % If the last envelope that we found is good, and the one we just tried
        % is not and has a value of n only one more than that envelope, the last envelope
        % is the boundary envelope.
        if isempty( env ) & ~isempty( last_good_envelope ) & max_n == last_good_envelope.n
            n = last_good_envelope.n;
            old_n = last_good_envelope.n;
            env = last_good_envelope;
        end
            
        % Quit if the old and new n are the same
        if old_n == n
            found_max_n = 1;
            max_n = n; 
            m = env.m + 1; 
            if ~isempty( env ) 
                % If the n value of the new envelope is not less than that of the 
                % previous envelope, the previous envelope is not on the boundary.
                % Thus, we will write over the previous envelope by not incrementing num_envelopes.
                if num_envelopes == 0 | ...
                   ( num_envelopes > 0 & env.n < boundary_envelopes( num_envelopes ).n )
                    num_envelopes = num_envelopes + 1;
                end
                boundary_envelopes( num_envelopes ) = env;
                if num_envelopes > 1 
                    diff = boundary_envelopes( num_envelopes - 1 ).n - boundary_envelopes( num_envelopes ).n;
                    if num_envelopes < num_last_envs + 2
                        diffs( num_envelopes - 1 ) = diff;
                    else 
                        diffs = [ diffs( 2 : end ) diff ];
                    end
                    %diffs
                end
            else
                done = 1;
            end
        end
        
        min_n_is_estimate = 0;
    end
end
