BUG: ufuncs created with np.frompyfunc do not recognize identity in reduce in numpy 1.25.X

Describe the issue:

Starting in numpy 1.25.0, it appears that the ufunc.reduce method does not recognize the identity argument for ufuncs created with np.frompyfunc

Reproduce the code example:

import numpy as np
import operator

print(np.__version__)

x = np.arange(10)
mask = x % 2 == 0
print(np.add.reduce(x, where=mask))

my_add = np.frompyfunc(operator.add, nin=2, nout=1, identity=0)
print(my_add.reduce(x, where=mask))

Error message:

Result in numpy 1.24.4:

1.24.4
20
20

Result in numpy 1.25.0:

1.25.0
20

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
[<ipython-input-7-a7d31df3c91a>](https://localhost:8080/#) in <cell line: 11>()
      9 
     10 my_add = np.frompyfunc(operator.add, nin=2, nout=1, identity=0)
---> 11 print(my_add.reduce(x, where=mask))

ValueError: reduction operation 'add (vectorized)' does not have an identity, so to use a where mask one has to specify 'initial'

Runtime information:

>>> import sys, numpy; print(numpy.__version__); print(sys.version)
1.25.0
3.10.12 (main, Jun 11 2023, 05:26:28) [GCC 11.4.0]
>>> print(numpy.show_runtime())
[{'numpy_version': '1.25.0',
  'python': '3.10.12 (main, Jun 11 2023, 05:26:28) [GCC 11.4.0]',
  'uname': uname_result(system='Linux', node='e87715ee78c7', release='5.15.109+', version='#1 SMP Fri Jun 9 10:57:30 UTC 2023', machine='x86_64')},
 {'simd_extensions': {'baseline': ['SSE', 'SSE2', 'SSE3'],
                      'found': ['SSSE3',
                                'SSE41',
                                'POPCNT',
                                'SSE42',
                                'AVX',
                                'F16C',
                                'FMA3',
                                'AVX2'],
                      'not_found': ['AVX512F',
                                    'AVX512CD',
                                    'AVX512_KNL',
                                    'AVX512_KNM',
                                    'AVX512_SKX',
                                    'AVX512_CLX',
                                    'AVX512_CNL',
                                    'AVX512_ICL']}},
 {'architecture': 'Haswell',
  'filepath': '/usr/local/lib/python3.10/dist-packages/numpy.libs/libopenblas64_p-r0-7a851222.3.23.so',
  'internal_api': 'openblas',
  'num_threads': 2,
  'prefix': 'libopenblas',
  'threading_layer': 'pthreads',
  'user_api': 'blas',
  'version': '0.3.23'}]
None

Context for the issue:

I am implementing jax.numpy.frompyfunc and had to special-case the tests in order to work around this bug: https://github.com/google/jax/pull/17275

Start by running the provided np.frompyfunc and ufunc.reduce reproduction on NumPy 1.24.4 and 1.25.0, then trace how the identity argument is handled for masked reductions. Done means the frompyfunc-created ufunc recognizes its identity with where=mask and the regression is covered by an appropriate test.