[FEATURE REQUEST] adding scipy.integrate

[h-milz]

Earlier this year, I added a numeric integration module to Circuitpython's ulab.scipy as provided in my repo and written in the Adafruit Playground blog. The integration algorithms

• quad - Tanh-Sinh, Sinh-Sinh and Exp-Sinh quadrature
• romberg - Romberg quadrature
• simpson - Adaptive Simpson quadrature
• quadgk - Adaptive Gauss-Kronrod (G10,K21) quadrature

were taken from this paper are are used with the kind written permission of the author. His code is under the MIT license, and so is my port to micropython-ulab. The latest version of the patch is in lvgl_micropython.diff in my repo. As you can see, everything is already according to the respective coding conventions, and the compilation can be controlled by setting ULAB_SCIPY_HAS_INTEGRATE_MODULE to 0 or 1. Needless to say that numerical integration makes limited sense with fp32 math but switching on fp64 e.g. for ESP32 is easy (in fact, easier in Micropython than in Circuitpython). The mpconfigport.h and mpconfigboard.h code is there as well. The code was tested in Circuitpython with an Adafruit Feather M4 Express, UM Feather S2, and a Teensy 4.0 using Adafruit's Feather adapter, and on LVGL-Micropython on a generic ESP32-S3 module. Compiling for rp2 failed outside of my code so far, didn't dig into this yet. But since there are no hardware dependencies, I do not expect any problems on any platform that has enough flash.

If you are interested, let me know and I'll fire up a PR. Disregard the LVGL part for the moment - I'll separate the patches anyway.

[v923z]

It's never occurred to me that one might need more than simple sums on a microcontroller, but if you think that it might have its uses, then by all means, please, create a PR.

A couple of comments:

1. It seems to me that you implement four methods in the integrate module. Could you, please, make sure that each method can be switched off in https://github.com/v923z/micropython-ulab/blob/master/code/ulab.h?
2. Can you add some documentation to https://github.com/v923z/micropython-ulab/tree/master/docs? This would be a separate jupyter notebook. You can take e.g., https://github.com/v923z/micropython-ulab/blob/master/docs/scipy-signal.ipynb, and replace everything after END_OF_DEFS by your stuff. The restructured text files in https://github.com/v923z/micropython-ulab/tree/master/docs/manual/source are generated by https://github.com/v923z/micropython-ulab/blob/master/docs/ulab-convert.ipynb. Let me know if you have difficulties with that!
3. Change the version number in

   micropython-ulab/code/ulab.c

   Line 36 in 303e8d7

   #define ULAB_VERSION 6.6.1

   to 6.7.0!
4. Make sure that you know what you do with complex numbers! The simplest solution would be to just bail out for that case. Note that you would need code for that only, if ULAB_SUPPORTS_COMPLEX is defined as 1 in

   micropython-ulab/code/ulab.h

   Line 35 in 303e8d7

   #ifndef ULAB_SUPPORTS_COMPLEX

   What I'm trying to say is that you should put the exception that you throw inside pre-processor instructions, so that you first don't run into difficulties with the compilation, and second, you don't waste flash space, when not required.

[h-milz]

Hi Zoltan,

as for item 1, I think it would be useful to have more than one algorithm. CPython scipy.integrate also has more than one, mainly because different algos have different uses, strengths and weaknesses, and we're really talking about just a few hundred bytes on binary space. But granted, ulab users may want to constrain their image sizes so - ok, no problem. I'll make sure that I write (or refer to) in the docs what each algo is good for, to allow for informed choices. Given that many recent microcontroller boards come with plenty of flash it should not be such a big issue anyway.

4. in fact, these algos don't work with complex numbers and will throw a TypeError if you try. This being said, the integration algorithms should also be available if 'ULAB_SUPPORTS_COMPLEX' is set.

PR coming in a minute.

[v923z]

OK, thanks!

[v923z]

4. in fact, these algos don't work with complex numbers and will throw a TypeError if you try.  This being said, the integration algorithms should also be available if 'ULAB_SUPPORTS_COMPLEX' is set.

I don't see where that happens. On the other hand, for the case, when a compiler doesn't support complex numbers natively, we process the complex numbers as a sum of a real and an imaginary, e.g., here: https://github.com/v923z/micropython-ulab/blob/master/code/numpy/carray/carray.c. I don't know if this is necessary in signal.integrate.

[h-milz]

Ah, there I was in error. I added code like

   if (!mp_obj_is_type(args[1].u_obj, &mp_type_float)) {
        mp_raise_TypeError(MP_ERROR_TEXT("second argument must be a float"));
   }      

in all relevant places. Coming with the next PR.

In fact, if you don't do this, mp_obj_get_float() will only pick the real part, which may not be what the user expects.

[h-milz]

OK, this test

from math import *
from ulab import scipy
i = scipy.integrate
func = "x * sin(x) * exp(x)"
f = lambda x: x * sin(x) * exp(x)
a=1
b=2
print (f"integrating f = {func} from {a} to {b}")
try:
    print (i.quad(f, a, b))
except Exception as e:
    print (e)
try:
    print (i.quad(f, a+1j, b))
except Exception as e:
    print (e)
try:
    print (i.quad(f, a, b+1j))
except Exception as e:
    print (e)
try:
    print (i.quad(f, a, b, levels=6.1))
except Exception as e:
    print (e)
try:
    print (i.quad(f, a, b, levels=(6+1j)))
except Exception as e:
    print (e)

now says

integrating f = x * sin(x) * exp(x) from 1 to 2
(7.11263821415851, 5.438231077315757e-14)
second argument must be a float or int
third argument must be a float or int
can't convert float to int
can't convert complex to int

[v923z]

This test is not yet in the PR, right?

[h-milz]

Nah, I just created a new one where all is in I think.

[h-milz]

@v923z In the documentation section for scipy.integrate, I see "UsageError: Cell magic %%micropython not found." in the code snippets. What might be broken?

[v923z]

Implemented in #699

[v923z]

@v923z In the documentation section for scipy.integrate, I see "UsageError: Cell magic %%micropython not found." in the code snippets. What might be broken?

Did you run the cells that define the cell magic? I assume you're on linux, so the RAM disc works. Otherwise, on windows, you might have to install RAM disc, and then modify the pointer in the relevant section of the magic definition.