BUG: numpy 1.22 incorrectly infers dimensions

[ketch]

Describe the issue:

Using f2py to wrap certain routines in Clawpack, on some systems we get a different calling signature using numpy 1.22 versus earlier versions of numpy. Some specifics are described here: clawpack/pyclaw#681

Some of the input arrays' dimensions involve two variables, so normally f2py requires that those values be passed in. Using 1.22, sometimes f2py sets up the wrapper so that both of those values are inferred, but their inferred values depend on each other. So even if one passes the values in, the routine does not work correctly.

I've found that rolling back to numpy version 1.18.5 fixes this issue. I haven't tested other versions in between.

I wonder if this is related to #20103 or #20696.

Reproduce the code example:

After installing Clawpack, run any PyClaw example.

Error message:

No response

NumPy/Python version information:

1.22.3

[charris]

ISTR another report along these lines. @HaoZeke Thoughts?

[HaoZeke]

ISTR another report along these lines. @HaoZeke Thoughts?

Yup, there have been a couple growing pains (#20721, #20709, #20853) from our new dimension solver #19805, but the fixes have been relatively quick to implement since it makes reasoning about the code much easier than before. Will look into this soon.

[HaoZeke]

This seems to be a non-starter on main.

git clone https://github.com/clawpack/clawpack.git
cd clawpack
git submodule init
git submodule update
pip install -e .
python -c "import clawpack.pyclaw.classic.classic1 as c; print(help(c))"

Gives:

Help on module clawpack.pyclaw.classic.classic1 in clawpack.pyclaw.classic:

NAME
    clawpack.pyclaw.classic.classic1

DESCRIPTION
    This module 'classic1' is auto-generated with f2py (version:1.23.0.dev0+1087.g0eaa40db3).
    Functions:
        limiter(maxm,num_ghost,mx,wave,s,mthlim,num_eqn=shape(wave, 0),num_waves=shape(wave, 1))
        philim = philim(a,b,meth)
        q,cfl = step1(num_ghost,mx,q,aux,dx,dt,method,mthlim,use_fwave,rp1,num_eqn=shape(q, 0),num_waves=shape(mthlim, 0),num_aux=shape(aux, 0),f=,wave=,s=,amdq=,apdq=,dtdx=,rp1_extra_args=())
    .

DATA
    __f2py_numpy_version__ = '1.23.0.dev0+1087.g0eaa40db3'
    limiter = <fortran object>
    philim = <fortran philim>
    step1 = <fortran object>

VERSION
    1.23.0.dev0+1087.g0eaa40db3

So it seems to be fixed on master. Will close in a day or so if there are no updates.

For reference, original bug as reported on the clawpack repository was that the dimensions were being inferred; so:

q,cfl = step1(num_ghost,mx,q,aux,dx,dt,method,mthlim,use_fwave,rp1,[num_eqn,num_waves,num_aux,f,wave,s,amdq,apdq,dtdx,rp1_extra_args])

Was instead:

q,cfl = step1(q,aux,dx,dt,method,mthlim,use_fwave,rp1,[num_eqn,num_waves,num_ghost,num_aux,mx,f,wave,s,amdq,apdq,dtdx,rp1_extra_args])

If these were cyclic dependencies (as @ketch mentions), then this was fixed by #21256.

[ketch]

In our experience the bug is not consistent. The function signature produced by f2py is different on different machines. But using numpy 1.18 always fixes it, and we've found no other way to fix it.

[HaoZeke]

Could you leave a reproducible example? Or a known machine where the incorrect results are produced? There should be no difference in the wrappers based on the underlying system. @ketch.

[ketch]

There should be no difference in the wrappers based on the underlying system.

I think we all agree on that! Some characteristics of machines with this problem are mentioned here: clawpack/pyclaw#681 (comment)

Let me know what else you need to know.

[WarrenWeckesser]

For what its worth, I experimented with this simplified version of the step1 subroutine:

foo.f90:

subroutine foo(num_eqn, num_ghost, mx, q)

    implicit double precision (a-h,o-z)
    double precision :: q(num_eqn, 1-num_ghost:mx+num_ghost)

!f2py intent(in,out) q
!f2py intent(in) num_eqn
!f2py intent(in) num_ghost
!f2py intent(in) mx

    return
end subroutine foo

I ran f2py with the command

    f2py --build-dir . -m foo foo.f90

With numpy 1.22.0 - 1.22.3, the signature of the wrapper generated by f2py is

    q = foo(q,[num_eqn,num_ghost,mx])

This is the unexpected result that is reported in this issue.

With the main development branch ('1.23.0.dev0+1087.g0eaa40db3'), and with all the older versions that I tested (1.19.0, 1.19.5, 1.20.0-1.20.3, 1.21.0-1.21.6), the signature is

    q = foo(num_ghost,mx,q,[num_eqn])

This is on a Linux machine, running Python 3.9.7 (miniconda), with all numpy versions except the main development branch installed with pip.

So for this platform and Python 3.9.7, it appears that only the 1.22 release series has the issue.

[ketch]

Okay, I can also confirm that this bug seems to be fixed in the current development version of numpy.

[HaoZeke]

Closing as fixed in 1.23... thanks @WarrenWeckesser and @ketch