function dx = projectBox(x, dx, lb, ub)
%
% projectBox is a helper function that moves x + dx to the nearest point in within the box 
% constraints, lb <= x + dx <= ub. If the bounds are already satisfied, the variable dx is unchanged.
% Nonfinite bounds lb(i) = -Inf and ub(i) = Inf are supported; however, empty bounds are not.
% 
%   Copyright 2020 The MathWorks, Inc.
%
% INPUTS:
%
%   x  - n-dimensional vector
%   lb - n-dimensional lower-bound constraints
%   ub - n-dimensional upper-bound constraints
% 
% ASSUMES:
% 
%   numel(x) == numel(lb)
%   numel(x) == numel(ub)
%   lb <= ub

%nolinear LS solvers do not support complex optimization and bounds simultaneously
if ~isreal(x)
    return
end

%Mathematically equivalent to P(x+dx) - x; however, numerically inaccurate.
%E.g. project(1,1e-16,-Inf,Inf) returns zero because 1e-16 < eps('double')
% xold = x(:);
% x(:) = max(x(:) + dx(:),lb(:));
% x(:) = min(x(:),ub(:));
% dx = x(:) - xold;
% dx = reshape(dx,size(x));

%Mathematically equivalent to P(x+dx) - x; however, skips needless
%cancellation (x+dx)-x for numerical reasons. See above. 
trialX = x(:) + dx(:);
activeLB = isfinite(lb) & trialX < lb;
activeUB = isfinite(ub) & trialX > ub;
dx(activeLB) = lb(activeLB) - x(activeLB);
dx(activeUB) = ub(activeUB) - x(activeUB);

end