def fun(self, epsilon, x0, x, b, ell, u, c, r, lambda0=None):
        """
        Computes the minimum norm necessary to reach the boundary. More precisely, we aim to solve the following
        optimization problem

            min ||delta||_2^2 s.t. lower <= x + delta <= upper AND b.dot(delta) = c

        Let's forget about the box constraints for a second, i.e.

            min ||delta||_2^2 s.t. b.dot(delta) = c

        The dual of this problem is quite straight-forward to solve,

            g(lambda, delta) = ||delta||_2^2 + lambda * (c - b.dot(delta))

        The minimum of this Lagrangian is delta^* = lambda * b / 2, and so

            inf_delta g(lambda, delta) = lambda^2 / 4 ||b||_2^2 + lambda * c

        and so the optimal lambda, which maximizes inf_delta g(lambda, delta), is given by

            lambda^* = 2c / ||b||_2^2

        which in turn yields the optimal delta:

            delta^* = c * b / ||b||_2^2

        To take into account the box-constraints we perform a binary search over lambda and apply the box
        constraint in each step.
        """
        N = x.shape[0]

        # new box constraints
        _ell, _u = self._Linf_bounds(x0, epsilon, ell, u)

        # initialize lambda
        _lambda = lambda0

        # compute delta and determine active set
        k = 0

        lambda_max, lambda_min = 1e10, -1e10

        # check whether problem is actually solvable (i.e. check whether boundary constraint can be reached)
        max_c = 0
        min_c = 0

        for n in range(N):
            if b[n] > 0:
                max_c += b[n] * (_u[n] - x[n])
                min_c += b[n] * (_ell[n] - x[n])
            else:
                max_c += b[n] * (_ell[n] - x[n])
                min_c += b[n] * (_u[n] - x[n])

        if c > max_c or c < min_c:
            return -np.inf, k, _lambda

        while True:
            k += 1
            _c = 0
            norm = 0
            _active_bnorm = 0

            for n in range(N):
                lam_step = _lambda * b[n] / 2
                if lam_step + x[n] < _ell[n]:
                    delta_step = _ell[n] - x[n]
                elif lam_step + x[n] > _u[n]:
                    delta_step = _u[n] - x[n]
                else:
                    delta_step = lam_step
                    _active_bnorm += b[n] ** 2

                _c += b[n] * delta_step
                norm += delta_step**2

            if 0.9999 * np.abs(c) - EPS < np.abs(_c) < 1.0001 * np.abs(c) + EPS:
                if norm > r**2:
                    return -np.inf, k, _lambda
                else:
                    return -epsilon, k, _lambda
            else:
                # update lambda according to active variables
                if _c > c:
                    lambda_max = _lambda
                else:
                    lambda_min = _lambda
                #
                if _active_bnorm == 0:
                    # update is stepping out of feasible region, fallback to binary search
                    _lambda = (lambda_max - lambda_min) / 2 + lambda_min
                else:
                    _lambda += 2 * (c - _c) / _active_bnorm

                dlambda = lambda_max - lambda_min
                if _lambda > lambda_max - 0.1 * dlambda or _lambda < lambda_min + 0.1 * dlambda:
                    # update is stepping out of feasible region, fallback to binary search
                    _lambda = (lambda_max - lambda_min) / 2 + lambda_min

        Computes the minimum norm necessary to reach the boundary. More precisely, we aim to solve the following
        optimization problem

            min ||delta||_2^2 s.t. lower <= x + delta <= upper AND b.dot(delta) = c

        Let's forget about the box constraints for a second, i.e.

            min ||delta||_2^2 s.t. b.dot(delta) = c

        The dual of this problem is quite straight-forward to solve,

            g(lambda, delta) = ||delta||_2^2 + lambda * (c - b.dot(delta))

        The minimum of this Lagrangian is delta^* = lambda * b / 2, and so

            inf_delta g(lambda, delta) = lambda^2 / 4 ||b||_2^2 + lambda * c

        and so the optimal lambda, which maximizes inf_delta g(lambda, delta), is given by

            lambda^* = 2c / ||b||_2^2

        which in turn yields the optimal delta:

            delta^* = c * b / ||b||_2^2

        To take into account the box-constraints we perform a binary search over lambda and apply the box
        constraint in each step.
        
art/attacks/evasion/brendel_bethge.py