bisection(f,lb,ub,target,options)

function [x, fx, exitFlag] = bisection(f,lb,ub,target,options)
% BISECTION Vectorized root-finding method.
%   [x,fVal,ExitFlag] = BISECTION(f,LB,UB,target,options) finds x +/- TolX
%   (LB < x < UB) such that f(x) = target +/- TolFun.
% 
%   Any or all of f(scalar), f(array), LB, UB, target, TolX, or TolFun may
%   be scalar or n-dim arrays. All non-scalar arrays must be the same size.
%   All outputs will be this size.
% 
%   x = BISECTION(f,LB,UB) finds the root(s) of function f on the interval
%   [LB, UB], i.e. finds x such that f(x) = 0 where LB <= x <= UB. f will
%   never be evaluated outside of the interval specified by LB and UB. f
%   should have only one root and f(UB) and f(LB) must bound it. elements
%   of x are NaN for instances where a solution could not be found.
% 
%   x = BISECTION(f,LB,UB,target) finds x such that f(x) = target.
% 
%   x = BISECTION(f,LB,UB,target,TolX) will terminate the search when the
%   search interval is smaller than TolX (TolX must be positive).
% 
%   x = BISECTION(f,LB,UB,target,options) solves with the default
%   parameters replaced by values in the structure OPTIONS, an argument
%   created with the OPTIMSET function. Used options are TolX and TolFun.
%   Note that OPTIMSET will not allow arrays for tolerances, so set the
%   fields of the options structure manually for non-scalar TolX or TolFun.
%
%   Default values are target = 0, TolX = 1e-6, and TolFun = 0.
% 
%   [x,fVal] = BISECTION(f,...) returns the value of f evaluated at x.
%
%   [x,fVal,ExitFlag] = BISECTION(...) returns an ExitFlag that describes
%   the exit condition of BISECTION. Possible values of elements of
%   ExitFlag and the corresponding exit conditions are
%
%       1   Search interval smaller than TolX.
%       2   Function value within TolFun of target.
%       3   Search interval smaller than TolX AND function value within 
%           TolFun of target.
%      -1   No solutions found.
% 
%   Note that there is no iteration limit. This is because BISECTION (with
%   a TolX that won't introduce numerical issues) is guaranteed to converge
%   if f is a continuous function on the interval [UB, LB] and f(x) -
%   target changes sign on the interval. The absolute error is halved at
%   each step so the method converges linearly. BISECTION is a very robust
%   root-finding method. However, Brent's method (such as implemented in
%   FZERO) can converge superlinearly and is as robust. FZERO also has more
%   features and input checking, so use BISECTION in cases where either the
%   optimization toolbox is unavailable or if FZERO would have to be
%   implemented in a loop to solve multiple cases, in which case BISECTION
%   will be much faster because of vectorization.
%
%   The user should define LB, UB, target, TolX, and TolFun for each
%   specific application using great care for the following reasons:
%      -There is no iteration limit, so given an unsolvable task, such as
%       TolX = TolFun = 0, BISECTION remains in an unending loop. 
%      -Spacing between very large floating point numbers is likely to be
%       greater than TolX. 
%      -There is no initial check to make sure that f(x) - target changes
%       sign between LB and UB.
%      -Very large or very small numbers can introduce numerical issues.
%
%   Example 1: find cube root of array 'target' without using NTHROOT
%       options = optimset('TolX',1e-9);
%       target = [(-100:.1:100)' (-1000:1:1000)'];
% 
%       tic;
%       xfz = zeros(size(target));
%       for ii = 1:numel(target)
%           xfz(ii) = fzero(@(x) x.^3-target(ii),[-20 20],options);
%       end
%       fzero_time = toc
% 
%       tic;
%       xbis = bisection(@(x) x.^3,-20,20,target,options);
%       bisection_time = toc
% 
%       sprintf('Bisection took %1.1f %% the time fzero needs.',...
%                   (100*bisection_time/fzero_time))
% 
%   Example 2: find roots by varying the function coefficients
%       [A, B] = meshgrid(linspace(1,2,10),linspace(4,12,30));
%       f=@(x) A.*x.^.2+B.*x.^.87-15;
%       xstar = bisection(f,0,5);
% 
%   See also FZERO, FMINBND, OPTIMSET.
% 
%   [x,fVal,ExitFlag] = BISECTION(f,LB,UB,target,options)

%   Author: Sky Sartorius
%       http://www.mathworks.com/matlabcentral/fileexchange/authors/101715

% --- Process inputs. ---
% Set default values
tolX    = 1e-6;
tolFun  = 0;
if nargin == 5
    if isstruct(options)
        if isfield(options,'TolX') && ~isempty(options.TolX)
            tolX = options.TolX;
        end 
        if isfield(options,'TolFun') && ~isempty(options.TolFun)
            tolFun = options.TolFun;
        end
    else
        tolX = options;
    end      
end
if nargin<4 || isempty(target); target=0; end


ub_in = ub; lb_in = lb; 
f = @(x) f(x) - target;

% --- Flip UB and LB if necessary. ---
isFlipped = lb>ub;
if any(isFlipped(:))
    ub(isFlipped) = lb_in(isFlipped);
    lb(isFlipped) = ub_in(isFlipped);
    ub_in = ub; lb_in = lb;
end

% --- Make sure everything is the same size for a non-scalar problem. ---
if isscalar(lb) && isscalar(ub)
    % Test if f returns multiple outputs for scalar input.
    if ~isscalar(target)
        ub = ub*ones(size(target));
    else
        jnk = f(ub); 
        if ~isscalar(jnk)
            ub = ub*ones(size(jnk));
        end
    end
end

% Check if lb and/or ub need to be made into arrays.
if isscalar(lb) && ~isscalar(ub)    
    lb = lb*ones(size(ub));
elseif ~isscalar(lb) && isscalar(ub)    
    ub = ub*ones(size(lb));
end

testconvergence
% --- Iterate ---
while any(stillNotDone(:))
    bigger  = fx.*f(ub) > 0;
    ub(bigger)= x(bigger);
    lb(~bigger)= x(~bigger);
    
    testconvergence;
end

    function testconvergence
        x=(ub+lb)/2;
        fx=f(x);
        outsideTolFun =  abs(fx)  > tolFun;
        outsideTolX =   (ub - lb) > tolX;
        stillNotDone = outsideTolX & outsideTolFun;
    end

% --- Check that f(x+tolX) and f(x-tolX) have opposite sign. ---
fu = f(min(x+tolX,ub_in)); 
fl = f(max(x-tolX,lb_in));
unboundedRoot = (fu.*fl) > 0;

% Throw out unbounded results if not meeting TolFun convergence criteria.
x(unboundedRoot & outsideTolFun) = NaN; 

% --- Catch NaN elements of UB, LB, target, or other funky stuff. ---
x(isnan(fx)) = NaN;

% --- Characterize results. ---
fx = fx + target;
if nargout > 2 
    exitFlag                                    = +~outsideTolX;
    exitFlag(~outsideTolFun)                    =  2;
    exitFlag(~outsideTolFun & ~outsideTolX)     =  3;
    exitFlag(isnan(x))                          = -1;
end

end

% V2: July     2010
% V3: December 2012