Complex casting warning

[theXYZT]

I am trying to use SPGL1 and noticed that line 204 of spgl1.py throws a ComplexWarning over dtype casting when running complex-valued problems: https://github.com/drrelyea/spgl1/blob/master/spgl1/spgl1.py#L204

Here's a toy example (a simple noisy DFT):

import numpy as np
from spgl1 import spgl1

def noise(shape):
    f = np.random.default_rng(seed=42).standard_normal
    return (f(shape) + 1j * f(shape)) / np.sqrt(2)

t = np.arange(128)
f = np.linspace(-0.5, 0.5, 32, endpoint=False)
A = np.exp(2j * np.pi * t[:, None] * f)  # DFT Matrix

y = (10 + 5j) * np.exp(2j * np.pi * t * 0.125)  # Complex Exponential
y += noise(len(t))  # Add complex noise with variance = 1.

x = spgl1(A, y, sigma=1, iscomplex=True)[0]
print(f[20], x[20])

This returns the output:

spgl1\spgl1.py:204: ComplexWarning: Casting complex values to real discards the imaginary part
  x *= s.astype(x.dtype)
0.125 (10.06103959036724+5.000239476329228j)

The printed output reveals that the problem was solved correctly (I get the expected coefficient at the right spot). But this also returns a ComplexWarning that might be worth dealing with in case it leads to unintended effects down the road.

[mrava87]

Hi @theXYZT, thanks for reporting this.

I think you may have found something a bit more subtle that just a warning. It looks like that we use np.all(np.isreal(xx)): to check if an array has complex numbers and in that case compute its abs value before passing it to oneprojector. However np.isreal(np.array(np.ones(4)+1j*np.zeros(4))) gives True which matematically makes sense but it doesnt do what we want. I changed all not np.all(np.isreal(xx)) into np.iscomplexobj(x) and the Warning is now gone :)

Basically in your problem the data is complex valued so when the code enters oneprojector the first thing that this routines does is:

s = np.sign(b)
b = np.abs(b)

Now b is complex valued but the sign of a complex-value number is not uniquely defined. Numpy seems to use this convention: For complex inputs, thesignfunction returns``sign(x.real) + 0j if x.real != 0 else sign(x.imag) + 0j``.. But we should have never got to the point if we stopped complex array to pass through :)