def _projected_sinkhorn(
        self, x: np.ndarray, x_init: np.ndarray, cost_matrix: np.ndarray, eps: np.ndarray
    ) -> np.ndarray:
        """
        The projected sinkhorn_optimizer.

        :param x: Current adversarial examples.
        :param x_init: An array with the original inputs.
        :param cost_matrix: A non-negative cost matrix.
        :param eps: Maximum perturbation that the attacker can introduce.
        :return: Adversarial examples.
        """
        # Normalize inputs
        normalization = x_init.reshape(x.shape[0], -1).sum(-1).reshape(x.shape[0], 1, 1, 1)
        x = x.copy() / normalization
        x_init = x_init.copy() / normalization

        # Dimension size for each example
        m = np.prod(x_init.shape[1:])

        # Initialize
        beta = np.log(np.ones(x.shape) / m)
        exp_beta = np.exp(-beta)

        psi = np.ones(x.shape[0])

        var_k = np.expand_dims(np.expand_dims(np.expand_dims(psi, -1), -1), -1)
        var_k = np.exp(-var_k * cost_matrix - 1)

        convergence = np.array([-np.inf])

        for _ in range(self.projected_sinkhorn_max_iter):
            # Block coordinate descent iterates
            x_init[x_init == 0.0] = EPS_LOG  # Prevent divide by zero in np.log
            alpha = np.log(self._local_transport(var_k, exp_beta, self.kernel_size)) - np.log(x_init)
            exp_alpha = np.exp(-alpha)

            beta = (
                self.regularization
                * np.exp(self.regularization * x)
                * self._local_transport(var_k, exp_alpha, self.kernel_size)
            )
            beta[beta > 1e-10] = np.real(lambertw(beta[beta > 1e-10]))
            beta -= self.regularization * x
            exp_beta = np.exp(-beta)

            # Newton step
            var_g = -eps + self._batch_dot(
                exp_alpha, self._local_transport(cost_matrix * var_k, exp_beta, self.kernel_size)
            )

            var_h = -self._batch_dot(
                exp_alpha, self._local_transport(cost_matrix * cost_matrix * var_k, exp_beta, self.kernel_size)
            )

            delta = var_g / var_h

            # Ensure psi >= 0
            tmp = np.ones(delta.shape)
            neg = psi - tmp * delta < 0

            while neg.any() and np.min(tmp) > 1e-2:
                tmp[neg] /= 2
                neg = psi - tmp * delta < 0

            psi = np.maximum(psi - tmp * delta, 0)

            # Update K
            var_k = np.expand_dims(np.expand_dims(np.expand_dims(psi, -1), -1), -1)
            var_k = np.exp(-var_k * cost_matrix - 1)

            # Check for convergence
            next_convergence = self._projected_sinkhorn_evaluation(
                x,
                x_init,
                alpha,
                exp_alpha,
                beta,
                exp_beta,
                psi,
                var_k,
                eps,
            )

            if (np.abs(convergence - next_convergence) <= 1e-4 + 1e-4 * np.abs(next_convergence)).all():
                break

            convergence = next_convergence

        result = (beta / self.regularization + x) * normalization

        return result

        The projected sinkhorn_optimizer.

        :param x: Current adversarial examples.
        :param x_init: An array with the original inputs.
        :param cost_matrix: A non-negative cost matrix.
        :param eps: Maximum perturbation that the attacker can introduce.
        :return: Adversarial examples.
        
art/attacks/evasion/wasserstein.py