8000 BUG: branch choices for `geomspace` on complex arguments · Issue #25644 · numpy/numpy · GitHub
[go: up one dir, main page]
Skip to content
BUG: branch choices for geomspace on complex arguments #25644
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom
Open
mcabbott opened this issue Jan 20, 2024 · 4 comments
Open

BUG: branch choices for geomspace on complex arguments #25644

mcabbott opened this issue Jan 20, 2024 · 4 comments
Labels
00 - Bug

Comments

@mcabbott
Copy link
mcabbott commented Jan 20, 2024 •
edited
Loading

Describe the issue:

The docs for np.geomspace state:

If the inputs or dtype are complex, the output will follow a logarithmic spiral in the complex plane. (There are an infinite number of spirals passing through two points; the output will follow the shortest such path.)

But the function does not appear to obey this rule. I'm not sure how it is in fact choosing branch cuts.

Searching a bit, this may be altered by #25441 but I believe that is not in any released version yet.

Reproduce the code example:

import numpy as np

np.geomspace(-1j, 0.001 + 1j, 5)   # semicircle to the right, via +1
np.geomspace(-1j, -0.001 + 1j, 5)  # semicircle to the left, via -1, as documented

x = 1.2 + 3.4j  # also x = 1+1j
delta = 0.01j
np.geomspace(x, -x + delta, 5)  # crosses real line near 3.5, middle point 3.4-1.2j
np.geomspace(x, -x - delta, 5)  # goes the same way!

Error message:

No error message. Numerical result is:

>>> x = 1.2 + 3.4j  # also x = 1+1j
>>> delta = 0.01j
>>> np.geomspace(x, -x + delta, 5)  # crosses real line near 3.5, middle point 3.4-1.2j
array([ 1.2       +3.4j       ,  3.2509223 +1.5538648j ,
        3.39499673-1.20000116j,  1.55032927-3.24738677j,
       -1.2       -3.39j      ])
>>> np.geomspace(x, -x - delta, 5)  # goes the same way!
array([ 1.2       +3.4j       ,  3.25445784+1.55740034j,
        3.40499673-1.20000115j,  1.56093588-3.25799337j,
       -1.2       -3.41j      ])

Python and NumPy Versions:

1.26.3
3.11.7 (main, Jan 16 2024, 14:42:22) [Clang 14.0.0 (clang-1400.0.29.202)]

Runtime Environment:

[{'numpy_version': '1.26.3',
  'python': '3.11.7 (main, Jan 16 2024, 14:42:22) [Clang 14.0.0 '
            '(clang-1400.0.29.202)]',
  'uname': uname_result(system='Darwin', node='ArmBook.local', release='21.6.0', version='Darwin Kernel Version 21.6.0: Sun Nov  6 23:29:57 PST 2022; root:xnu-8020.240.14~1/RELEASE_ARM64_T8101', machine='arm64')},
 {'simd_extensions': {'baseline': ['NEON', 'NEON_FP16', 'NEON_VFPV4', 'ASIMD'],
                      'found': ['ASIMDHP'],
                      'not_found': ['ASIMDFHM']}},
 {'architecture': 'armv8',
  'filepath': '/opt/homebrew/lib/python3.11/site-packages/numpy/.dylibs/libopenblas64_.0.dylib',
  'internal_api': 'openblas',
  'num_threads': 8,
  'prefix': 'libopenblas',
  'threading_layer': 'pthreads',
  'user_api': 'blas',
  'version': '0.3.23.dev'}]
None

Context for the issue:
@eendebakpt
Copy link
Contributor

I can confirm the documentation does not match the implementation. I am not sure what the use case is for a complex np.geomspace, or why the shortest path is preferred, so not sure which one needs to be changed.

A fix 8000 that modifies the code so the shortest path is used, is in main...eendebakpt:geomspace, but it is not very elegant (it ad-hoc modifies the phase of the result after the main calculation). Another option would be to deal with the radius and phase of the start and stop separately. For the radius use the normal logspace, for the phase part use a linspace (making sure to take the shortest path module 2*pi).

@endolith Any ideas on this?

@mhvk
Copy link
Contributor
mhvk commented Jan 28, 2024

I did a quick test on -dev, which includes my #25441, and now it seems more reasonable,

In [2]: np.geomspace(x, -x + delta, 5)
Out[2]: 
array([ 1.2       +3.4j       , -1.5538648 +3.2509223j ,
       -3.39499673+1.20000116j, -3.24738677-1.55032927j,
       -1.2       -3.39j      ])

In [3]: np.geomspace(x, -x - delta, 5)
Out[3]: 
array([ 1.2       +3.4j       ,  3.25445784+1.55740034j,
        3.40499673-1.20000115j,  1.56093588-3.25799337j,
       -1.2       -3.41j      ])

I remember my pleasant surprise about "... how the code in geomspace is substantially simplified by using the new definition!", but it seems it was actually more than that!

Would you be able to confirm that the issue is solved on numpy-dev?

@endolith
Copy link
Contributor

@eendebakpt I had no use case in mind for complex inputs, and wasn't sure if they even needed to be supported. I think the "shortest such path" description was more an observation than a specification. I don't care what it does with complex inputs. The main goal of geomspace was just more intuitive arguments than logspace.

@seberg seberg mentioned this issue Jan 29, 2024
Merged
@eendebakpt
Copy link
Contributor

@mhvk I created a regression test that passes: #25715, so it seems the issue has been resolved on current dev.

mhvk pushed a commit that referenced this issue Jan 29, 2024
@eendebakpt
[TEST] Add regression test for np.geomspace (#25715)
8a36b95 
See gh-25644, which showed that numpy 1.26 did not take the
shorted logarithmic spiral. Coincidentally, though, this was
already fixed by gh-25441, so this commit just adds a regression
test.
@eendebakpt eendebakpt mentioned this issue Apr 2, 2024
Closed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
00 - Bug
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
4 participants
@endolith @eendebakpt @mhvk @mcabbott
0