Repository: scipy/scipy
Pull request: BUG: fft.fht: set `u.imag[-1] = 0` only when `n` is even
Closes #21661.


Problem statement:
BUG: fft.fht: should set `u.imag[-1] = 0` only when `n` is even
### Describe your issue.

In the source of scipy.fft.fht, in the function fhtcoeff(...) in  _fftlog_backend.py, one has to ensure that the array u has a u[-1] which is purely real when the data is of even size n. That requirement stems from the fact that one wants the fht transform of real data to end up real. The argument is spelled out in the original paper. 

However, the code sets u[-1].imag=0 even for odd n, which I'm pretty sure is wrong. So, I believe that on lines 94-95:

    # fix last coefficient to be real
    u.imag[-1] = 0
    
there should be a check to ensure that u.imag[-1]=0 only for even n. 


### Reproducing Code Example

```python
The issue has to do with mathematical correctness of the implemented algorithm. The code as it is doesn't crash.
```


### Error message

```shell
none
```


### SciPy/NumPy/Python version and system information

```shell
1.14.1
```

Hints from issue discussion:
cc @ntessore 
Oh, here is a code example showing the issue. Using the example in the documentation, but running it with an odd n=129 produces much larger relative errors relative to the analytical result. For even n=128 or 130, the error is much smaller. Here is the example with n=129:
```
import numpy as np
from scipy import fft
import matplotlib.pyplot as plt

mu = 0.0                     # Order mu of Bessel function
r = np.logspace(-7, 1, 129)  # Input evaluation points
dln = np.log(r[1]/r[0])      # Step size
offset = fft.fhtoffset(dln, initial=-6*np.log(10), mu=mu)
k = np.exp(offset)/r[::-1]   # Output evaluation points

def f(x, mu):
    """Analytical function: x^(mu+1) exp(-x^2/2)."""
    return x**(mu + 1)*np.exp(-x**2/2)


def f(x, mu):
    """Analytical function: x^(mu+1) exp(-x^2/2)."""
    return x**(mu + 1)*np.exp(-x**2/2)

a_r = f(r, mu)
fht = fft.fht(a_r, dln, mu=mu, offset=offset)
a_k = f(k, mu)
rel_err = abs((fht-a_k)/a_k)
figargs = {'sharex': True, 'sharey': True, 'constrained_layout': True}
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), **figargs)
ax1.set_title(r'$r^{\mu+1}\ \exp(-r^2/2)$')
ax1.loglog(r, a_r, 'k', lw=2)
ax1.set_xlabel('r')
ax2.set_title(r'$k^{\mu+1} \exp(-k^2/2)$')
ax2.loglog(k, a_k, 'k', lw=2, label='Analytical')
ax2.loglog(k, fht, 'C3--', lw=2, label='FFTLog')
ax2.set_xlabel('k')
ax2.legend(loc=3, framealpha=1)
ax2.set_ylim([1e-10, 1e1])
ax2b = ax2.twinx()
ax2b.loglog(k, rel_err, 'C0', label='Rel. Error (-)')
ax2b.set_ylabel('Rel. Error (-)', color='C0')
ax2b.tick_params(axis='y', labelcolor='C0')
ax2b.legend(loc=4, framealpha=1)
ax2b.set_ylim([1e-9, 1e-3])
plt.show()

```
FYI: Modified local copy of code to fix the u.imag[-1]=0 only for even n, and now relative error is small for odd n.

This is unrelated, but is in the same code: the documentation example at [https://docs.scipy.org/doc/scipy/reference/generated/scipy.fft.fht.html](https://docs.scipy.org/doc/scipy/reference/generated/scipy.fft.fht.html) has a bug. It should call fhtoffset as:
offset = fft.fhtoffset(dln,mu, initial=-6*np.log(10))
Currently it's being called with offset = fft.fhtoffset(dln, initial=-6*np.log(10), mu=mu), which does not produce correct result  as the argument order is wrong.
> Currently it's being called with offset = fft.fhtoffset(dln, initial=-6np.log(10), mu=mu), which does not produce correct result as the argument order is wrong.

In what way is the result unexpected? You should be able to pass keyword arguments in any order.

EDIT: indeed:

```python
In [10]: fhtoffset(dln, initial=-6*np.log(10), mu=mu)
Out[10]: np.float64(-13.879121901447917)

In [11]: fhtoffset(dln, mu=mu, initial=-6*np.log(10))
Out[11]: np.float64(-13.879121901447917)
```
Apologies! I had overwritten mu while testing. The doc string is fine.

> > Currently it's being called with offset = fft.fhtoffset(dln, initial=-6np.log(10), mu=mu), which does not produce correct result as the argument order is wrong.
> 
> In what way is the result unexpected? You should be able to pass keyword arguments in any order.
> 
> EDIT: indeed:
> 
> ```python
> In [10]: fhtoffset(dln, initial=-6*np.log(10), mu=mu)
> Out[10]: np.float64(-13.879121901447917)
> 
> In [11]: fhtoffset(dln, mu=mu, initial=-6*np.log(10))
> Out[11]: np.float64(-13.879121901447917)
> ```


Thanks for catching this! It must have been lost in translation, as my original code has a `if not n & 1:` there.
would you like to submit a fix @stassev ?

Commit message: BUG: fft.fht: should set u.imag[-1] = 0 only when n is even #21661

Patch to review:
diff --git a/scipy/fft/_fftlog_backend.py b/scipy/fft/_fftlog_backend.py
index d3d05848ac09..414f16072dc3 100644
--- a/scipy/fft/_fftlog_backend.py
+++ b/scipy/fft/_fftlog_backend.py
@@ -92,7 +92,8 @@ def fhtcoeff(n, dln, mu, offset=0.0, bias=0.0, inverse=False):
     np.exp(u, out=u)
 
     # fix last coefficient to be real
-    u.imag[-1] = 0
+    if (n%2==0):
+        u.imag[-1] = 0
 
     # deal with special cases
     if not np.isfinite(u[0]):


Code review comments:
diff_hunk: @@ -92,7 +92,8 @@ def fhtcoeff(n, dln, mu, offset=0.0, bias=0.0, inverse=False):
     np.exp(u, out=u)
 
     # fix last coefficient to be real
-    u.imag[-1] = 0
+    if (n%2==0):
path: scipy/fft/_fftlog_backend.py
text: @user1:
```suggestion
    if n % 2 == 0:
```
Small nitpick.