Complex2

[v923z]

This is the re-based version of #366, which will be closed now.

[v923z]

@wired8 Could you, please, check, whether you still have the error that you mentioned in #366 (comment) ?

[wired8]

@wired8 Could you, please, check, whether you still have the error that you mentioned in #366 (comment) ?

Looks good, no more kernal panic.

Doing some testing with the following code:

>>> N = 5
>>> z = numpy.array([])
>>> m = numpy.arange(-N+1, N, 2)
>>> p = -numpy.exp(1j * np.pi * m / (2 * N))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'complex' object isn't iterable

Can we add support for iterables over complex types?

[v923z]

@wired8

Doing some testing with the following code:

>>> N = 5
>>> z = numpy.array([])
>>> m = numpy.arange(-N+1, N, 2)
>>> p = -numpy.exp(1j * np.pi * m / (2 * N))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'complex' object isn't iterable

Can we add support for iterables over complex types?

The problem is that the exponential function is not defined for complexes at the moment. I could implement that, but the question really is, what would be useful in this regard. There are some 20 universal functions, and I am not sure, whether implementing everything for complex types makes sense. The exponential function should definitely be included. What else? Here is a rudimentary list.

• unary operators
• binary operators
  • addition
  • multiplication
  • division
  • subtraction
  • equal
  • not equal
• slicing
• iterators
• all
• any
• exp function
• sqrt function
• imag function
• conjugate function
• real function
• compress function
• sort_complex function
• extend convolve for the complex case
• FFT in a user-configurable way, so that complete numpy-compatibility can be achieved.

EDIT: actually, I think that the error that you encountered comes from here 1j * np.pi * m: multiplication (and binary operators in general) is not implemented for complexes.
And one more comment, though, it is somewhat unrelated: the expression np.pi * m / (2 * N) is very expensive, because you create two extra arrays, one for each multiplication. It is much better to use something like

m = np.linspace(-np.pi*(N+1)/(2*N), np.pi*N/(2*N), num=N)

or something similar.

[wired8]

The list you provided is a great start. Implementing binary, exp would give us partial butterworth filter support.

[wired8]

np.exp(-3.141592653589793j)
(-1+8.742278e-08j)

np.exp([-3.141592653589793j])
Traceback (most recent call last):
File "", line 1, in
TypeError: can't convert complex to float

[v923z]

np.exp(-3.141592653589793j)
(-1+8.742278e-08j)

np.exp([-3.141592653589793j])
Traceback (most recent call last):
File "", line 1, in
TypeError: can't convert complex to float

I think this is OK that way, at least, intentional. If you need the exponential of a complex list, you can cast it first as

a = np.array([3.14j], dtype=np.complex)
np.exp(a)

The reason for this is that, when applying a function to an object, first you have to figure out,
whether the output is going to be complex, or floating point. But you can learn that only by inspecting
each element, so you can do that in a double pass only. I am not sure it's worth the trouble.

Floats are simpler in this regard, because applying a mathematical function to a number always
produces a float, so you can automatically assume that the output is of float dtype.

[wired8]

numpy.convolve

We can definitely extend convolve for the complex case.

numpy.conjugate

This is trivial, although, you can already do this:

a = np.array(..., dtype=np.complex)
b = real(a) - imag(a) * 1j

numpy.compress

Is it this one https://numpy.org/doc/stable/reference/generated/numpy.compress.html?

Yes

numpy.sort_complex

This could be added. As far as I see, this works with 1D arrays only, so the implementation is almost trivial.

The following is the guts of the zpk2tf method, so if we can implement the aforementioned we should be able to support this function which is the last part of the iirfilter.

 if issubclass(b.dtype.type, numpy.complexfloating):
        # if complex roots are all complex conjugates, the roots are real.
        roots = numpy.asarray(z, complex)
        pos_roots = numpy.compress(roots.imag > 0, roots)
        neg_roots = numpy.conjugate(numpy.compress(roots.imag < 0, roots))
        if len(pos_roots) == len(neg_roots):
            if numpy.all(numpy.sort_complex(neg_roots) ==
                         numpy.sort_complex(pos_roots)):
                b = b.real.copy()

[v923z]

 if issubclass(b.dtype.type, numpy.complexfloating):
        # if complex roots are all complex conjugates, the roots are real.
        roots = numpy.asarray(z, complex)
        pos_roots = numpy.compress(roots.imag > 0, roots)
        neg_roots = numpy.conjugate(numpy.compress(roots.imag < 0, roots))
        if len(pos_roots) == len(neg_roots):
            if numpy.all(numpy.sort_complex(neg_roots) ==
                         numpy.sort_complex(pos_roots)):
                b = b.real.copy()

Should real/imag be properties of the ndarray? Otherwise, you have to do

b = real(b)

That also copies the content automatically.

[wired8]

 if issubclass(b.dtype.type, numpy.complexfloating):
        # if complex roots are all complex conjugates, the roots are real.
        roots = numpy.asarray(z, complex)
        pos_roots = numpy.compress(roots.imag > 0, roots)
        neg_roots = numpy.conjugate(numpy.compress(roots.imag < 0, roots))
        if len(pos_roots) == len(neg_roots):
            if numpy.all(numpy.sort_complex(neg_roots) ==
                         numpy.sort_complex(pos_roots)):
                b = b.real.copy()

Should real/imag be properties of the ndarray? Otherwise, you have to do

b = real(b)

That also copies the content automatically.

Keeps it to spec with numpy, so why not!

[v923z]

Should real/imag be properties of the ndarray? Otherwise, you have to do

b = real(b)

That also copies the content automatically.

Keeps it to spec with numpy, so why not!

Well, because if you want to have numpy, then use numpy! Implementing extra features (i.e., those whose functionality can be gotten by other means) adds to the firmware size, eats RAM, and consumes development resources. By the way, I think you don't need to copy in b.real.copy(), because b.real should return a new copy, and not a view. At least, if I read the specs correctly: https://numpy.org/doc/stable/reference/generated/numpy.ndarray.real.html.

But to address your comment, b.real seems to be simply an alias for real(b), so we could add it with little overhead. b.imag, too.

[wired8]

@v923z awesome work, I am now able to get the correct output for a butter filter with one small caveat that the all function doesn't support complex yet. I do have a work around for now. https://micropython-ulab.readthedocs.io/en/latest/numpy-functions.html#all

[v923z]

@v923z awesome work, I am now able to get the correct output for a butter filter with one small caveat that the all function doesn't support complex yet.

@wired8 There, you have it: d1b3d40

[wired8]

Possible bug with np.all

p = np.array([1j,2j,3j,4j], dtype=np.complex)
q = np.array([1j,2j,3j,4j], dtype=np.complex)
print(np.all(p==q))

>>> False

Expect True

[wired8]

Possible bug with sorting and conjugate

roots = np.array([0.98-0.01j,0.99+0.01j,0.98+0.01j,0.99-0.01j], dtype=np.complex)
p = np.sort_complex(roots)
q = np.sort_complex(np.conjugate(roots))
print(p)
print(q)

>>> array([0.9800000000000001-0.01j, 0.9800000000000001+0.01j, 0.99-0.01j, 0.99+0.01j], dtype=complex)
>>> array([0.9800000000000001+0.01j, 0.9800000000000001-0.01j, 0.99+0.01j, 0.99-0.01j], dtype=complex)

Expect:

[0.98-0.01j 0.98+0.01j 0.99-0.01j 0.99+0.01j]
[0.98-0.01j 0.98+0.01j 0.99-0.01j 0.99+0.01j]

[v923z]

Possible bug with np.all

p = np.array([1j,2j,3j,4j], dtype=np.complex)
q = np.array([1j,2j,3j,4j], dtype=np.complex)
print(np.all(p==q))

>>> False

Expect True

I think the problem is that == is not implemented for complex. I could add that. That is not so expensive, given that the operator is commutative.

I should probably also raise an exception for <=, <, and similar binary operators, if one of the operands is a complex array.

[v923z]

Possible bug with sorting and conjugate

roots = np.array([0.98-0.01j,0.99+0.01j,0.98+0.01j,0.99-0.01j], dtype=np.complex)
p = np.sort_complex(roots)
q = np.sort_complex(np.conjugate(roots))
print(p)
print(q)

>>> array([0.9800000000000001-0.01j, 0.9800000000000001+0.01j, 0.99-0.01j, 0.99+0.01j], dtype=complex)
>>> array([0.9800000000000001+0.01j, 0.9800000000000001-0.01j, 0.99+0.01j, 0.99-0.01j], dtype=complex)

Expect:

[0.98-0.01j 0.98+0.01j 0.99-0.01j 0.99+0.01j] [0.98-0.01j 0.98+0.01j 0.99-0.01j 0.99+0.01j]

Thanks for reporting the issue! I'll sort it out. (Pun intended.)

[wired8]

Hey @v923z your fixes work and I now get the same expected output from a butter filter in ulab as I do in numpy/scipy. Thank you so much for all you hard work. I want to contribute my partial implementation of these libraries back into ulab, I'm guessing I add them as snippets?

I was thinking of adding a few libraries under snippets/scipy, snippets/numpy and some tests to boot. Please advise on how to proceed.

[v923z]

Hey @v923z your fixes work and I now get the same expected output from a butter filter in ulab as I do in numpy/scipy. Thank you so much for all you hard work.

Not at all. I am glad that it works now. There were a couple of smaller glitches that I discovered after you had reported the problem with sort_complex, so your tests were doubly useful. Thanks for sticking with it till the very end!

I want to contribute my partial implementation of these libraries back into ulab, I'm guessing I add them as snippets?

If the implementation is in python, then I would suggest the snippets. If it is in C, then it should be part of code. I don't know, whether we can gain anything in speed or RAM, if we implement at least chunks in C. As far as I understand numpy, it has helper functions written in C that you don't usually call from python directly, only via the python implementation of a particular function. We could call them low-level functions that don't have a pythonic interface. I would also be OK with that. If you can point me to the relevant pieces of the scipy code, we could discuss what is most sensible path.

One argument against the C implementation (beyond the obvious) is that it requires flash space, even if you don't use the function. So, snippets are somewhat more flexible.

I was thinking of adding a few libraries under snippets/scipy, snippets/numpy and some tests to boot. Please advise on how to proceed.

You have just described what to do. I would like to merge this branch, and release 4.0 very soon, but that doesn't really interfere with your stuff. You can PR against this, or the master branch, whenever you are ready. I am quite certain that contributions of this kind would be very useful for the wider community.

Actually, if you have any comments on the release draft, please, let me know!

Also, I think #465 is obsolete now. Can we close that?