code: def test_apply_gufunc_axes_two_kept_coredims():
      rng = da.random.default_rng()
      a = rng.normal(size=(20, 30), chunks=(10, 30))
      b = rng.normal(size=(10, 1, 40), chunks=(5, 1, 40))
  
      def outer_product(x, y):
          return np.einsum("i,j->ij", x, y)
  
      c = apply_gufunc(
          outer_product, "(i),(j)->(i,j)", a, b, vectorize=True
      )
      assert c.compute().shape == (10, 20, 30, 40)
  
imports: from __future__ import annotations
  import numpy
  import pytest
  from numpy.testing import assert_equal
  import dask.array
  from dask.array.core import Array
  from dask.array.gufunc import _parse_gufunc_signature, _validate_normalize_axes, apply_gufunc, as_gufunc, gufunc
  from dask.array.utils import assert_eq
label_desc: Use image node_id 41 for calling a global function with example usage: test_apply_gufunc_axes_two_kept_coredims() and returns: np