[FEATURE REQUEST] Add numpy.polynomial.polynomial.Polynomial.roots()

[jrmoserbaltimore]

Describe the solution you'd like ulab already has polyfit(), but not Polynomial.

I would like to see Polynomial implementing Polynomial.fit() and Polynomial.roots(), as well as Polynomial.coef as in numpy.

Additional context For certain sensors, particularly Hall effect sensors, the reading is a function of the square of the distance to an object, or some other polynomial relationship. With calibration data, a curve can be generated with Polynomial.fit() or polyfit(); to get the real-world reading, subtract the sensor reading from the final coefficient (e.g. ax**2+bx+c, subtract the raw sensor reading from c) and find the roots. For a Hall sensor, the smallest root is then the real-world distance.

ulab doesn't implement roots(), so this is not possible today. A partial implementation of Polynomial would be good for compatibility with the modern numpy library.

[v923z]

@jrmoserbaltimore Just to make sure that we are on the same page, you mean some of the functions from this package: https://numpy.org/doc/stable/reference/routines.polynomials.package.html

ulab doesn't implement roots(), so this is not possible today. A partial implementation of Polynomial would be good for compatibility with the modern numpy library.

As far as I understand, you want to fit a quadratic polynomial, right? But you can find the roots of that via a closed formula. I wonder, whether the following snippet would actually address the issue

epsilon = 1e-6

def root2(p):
    if p[2] < epsilon:
        raise ValueError('input is a linear polynomial')
    if p[1] ** 2 < 4 * p[0] * p[2]:
        raise ValueError('polynomial has no real roots')

    sqrt_root = sqrt(p[1] ** 2 - 4 * p[0] * p[2])
    return (-p[1] + sqrt_root) / (2 * p[2]),  (-p[1] - sqrt_root) / (2 * p[2]),  

p = np.polyfit(data, 2)
data = p[0] - data

From your description, I don't actually see, where you need the roots, but you would get it with the root2 function.

I am not against adding new functions, but in this case, I have the feeling that the functionality can easily be implemented in python itself. On the other hand, finding the roots of a general polynomial of order n is not trivial, in fact, outright impossible, if n > 4.

[jrmoserbaltimore]

True. Part of this is just looking at cross-platform compatibility, writing code to function in both without writing wrappers around ulab to replace numpy functions that are called on CPython. polyfit is deprecated in numpy and all.

If you did implement it, it would be viable to just raise an exception if finding the root with order n>2 or whatever you think is too complex for the poor microcontroller trying to work it out. Part of a compatible strict subset is reasonably that not everything is implemented, or implemented fully.

[v923z]

@jrmoserbaltimore

True. Part of this is just looking at cross-platform compatibility, writing code to function in both without writing wrappers around ulab to replace numpy functions that are called on CPython. polyfit is deprecated in numpy and all.

OK, I haven't known about this, and this is definitely a valid point. But this would simply mean that we would have to re-label the python-facing functions, and we can keep the implementation. So this would not be too much of an overhead. Adding a module as in numpy would have the benefit that extra functions (like differentiation, integration etc.) could be added later. The bottom line is, this could definitely be useful, and I could pull this off in a couple of days.

If you did implement it, it would be viable to just raise an exception if finding the root with order n>2 or whatever you think is too complex for the poor microcontroller trying to work it out.

That's fine, but I see another issue here, namely, what should happen, if the roots are complex? ulab doesn't support complex arrays for the simple reason that the overhead is really too much: first, all functions would have to handle complexes, second, and this is the more problematic part, if you add, subtract, multiplies etc. two arrays, you have to be able to resolve the C types of both arrays, which means that the cost of all binary operations scales with the square of data types. At the moment we have five and a half (Booleans are the half), meaning 25 possible combinations. It is possible to take shortcuts here and there, but you can't really save too much. Now, if we wanted to support complexes, then all binary operators would take 50% more flash space.

If we wanted to keep the number of supported dtypes at 5, what should we do with functions that produce complex outputs?

Finally, could you spell out the differences between Polynomial.roots, and polyroots? It seems the me that Polynomial is actually a class with a number of methods, while polynomial is a module with a number of methods. Would it be OK, if we tried to stick with the module, and forgo the class? In other words, is there anything in the class that cannot be accessed via the module?

[jrmoserbaltimore]

It's considered old API but I may be misreading. On a closer look it looks like they just moved the location, but it's the same function. I want to say I'm a politician, not a programmer, but it's more accurate to say I'm just a bad programmer 😄

That's fine, but I see another issue here, namely, what should happen, if the roots are complex? ulab doesn't support complex arrays for the simple reason that the overhead is really too much

Shrug, then don't. I'm mostly interested in being able to do stuff like:

try:
    from numpy.polynomial.polynomial import Polynomial
except ImportError:
    from ulab.numpy import Polynomial

or wherever in ulab it goes in the hierarchy. If there are limitations to ulab, I don't want to write a bunch of wrappers, or use facility A versus facility B with different syntax, throughout the code. Making certain considerations on certain limitations is usually okay, so long as the client code can function within those limitations—kind of like how a lot of code can import uzlib as zlib and avoid CPython-only calls.

is there anything in the class that cannot be accessed via the module?

Probably not. It'd be viable to include Polynomial written in Python as a wrapper around the methods, and not build it into ulab proper, as something that could be installed separately.

[v923z]

It's considered old API but I may be misreading. On a closer look it looks like they just moved the location, but it's the same function.

OK, so we could simply re-package the functions that we already have, and add the roots for n <= 3. n = 4 is already ugly, beyond that is possible only numerically (no closed formula). But if we do this, then you don't get the Polynomial class. Is that OK?

That's fine, but I see another issue here, namely, what should happen, if the roots are complex? ulab doesn't support complex arrays for the simple reason that the overhead is really too much

Shrug, then don't.

But you still can't avoid the difficulty of having complex roots of a real polynomial. Something like this:

x**2 + 1 = 0

polynomial.polyroots returns an array, which must then be able to hold complex entries.

I'm mostly interested in being able to do stuff like:

try:
    from numpy.polynomial.polynomial import Polynomial
except ImportError:
    from ulab.numpy import Polynomial

or wherever in ulab it goes in the hierarchy.

This last line would fail in ulab. I don't think you can do dot imports in micropython. If you know of an example, I would like to know about it. What you would do is

from ulab import numpy as np

np.polynomial....

If there are limitations to ulab, I don't want to write a bunch of wrappers, or use facility A versus facility B with different syntax, throughout the code. Making certain considerations on certain limitations is usually okay, so long as the client code can function within those limitations—kind of like how a lot of code can import uzlib as zlib and avoid CPython-only calls.

I wanted to bring this up: micropython implements a lot of stuff as python code, even if it is C in CPython. A good example is itertools. The reason for that is that you can then control the firmware size, and you include things via python, if you need it. This issue has also come up in relation to ulab, and we concluded that we could perhaps add a snippets folder with wrappers that you describe. An example is https://github.com/v923z/micropython-ulab/issues/278

[jrmoserbaltimore]

Fair enough. I didn't find the roots function in ulab, have to look further. Didn't realize it didn't do degree-2 polynomials either.

I'll think on this some more; lots of stuff I didn't consider or realize here.

[v923z]

Fair enough. I didn't find the roots function in ulab, have to look further.

No, it is not there, so don't look for it!

[v923z]

@jrmoserbaltimore Here is a new issue, you could comment on the specific questions there: https://github.com/v923z/micropython-ulab/issues/364

[v923z]

@jrmoserbaltimore I have just opened a draft under https://github.com/v923z/micropython-ulab/pull/366. With that, we should be able to support the roots function. I don't know, what exactly you want to do with the results, so it would be great, if you could comment on what exactly you need afterwards. Finding the roots is one thing, but I assume that you want to use the results somewhere.

Also, you should think about what you do, if you expect real output, but the results are slightly imaginary, i.e., you have something like

[1.0, 2.0, 3.0+0.00001j]

The 0.00001j is clearly a truncation error in this case.

In other words, do you need facilities to guard against such cases? If so, what should they be? Are you going to implement those in python, or do you need something in C? And so on.