Bumbling at Writing a Coroutine

[Chris-24MAY60]

Hello, I am trying to understand how to write an ansynchronous coroutine. I have written a script, named uart.py, that I wish to run as my coroutine. The (larger) program runs beautifully! (but alas, only once) My Microdot server is still running, happily serving web pages and-so my suspicion is that my uart.py script still isn't a proper coroutine, although I have added the recommended await asyncio.sleep(0) line in attempt to ensure I've made it non-blocking. Sadly, no bueno. Below is my code.

import uasyncio as asyncio
import asyncio
import machine
from machine import UART
import time
from time import sleep_ms
import binascii
from binascii import hexlify
import network

async def process_uart():
#code to process the file
try:
uart2 = UART(2, tx=17, rx=18) # 18:RX 17:TX
uart2.init(baudrate=115200, bits=8, parity=None, stop=1)

    buf = bytearray(437)
    mv = memoryview(buf)
    idx = 0

    while idx < len(buf):
        if uart2.any():
            bytes_read = uart2.readinto(mv[idx:])
            #print('Got {} bytes of data'.format(bytes_read), hexlify(buf[idx:idx+bytes_read], b','))
            idx += bytes_read
            hex_data = hexlify(buf[3:435])
        #print('Final buffer =', hexlify(final_buf, b','))
    print(hex_data)#we need hex_data to be one long string of characters - but we have spaces between each byte

    rows = []
    for i in range(0, len(hex_data), 16):  # 16 characters represent 16 nibbles in hex
        fourth_nibble = int(hex_data[i+3:i+4], 16)
        eighth_nibble = int(hex_data[i+7:i+8], 16)
        tenth_nibble = int(hex_data[i+9:i+10], 16)

        # Calculate the value in the fourth column as per the provided formula
        value = (
            int(hex_data[i+14:i+15], 16) * 1048576 +
            int(hex_data[i+15:i+16], 16) * 65536 +
            int(hex_data[i+12:i+13], 16) * 4096 +
            int(hex_data[i+13:i+14], 16) * 256 +
            int(hex_data[i+10:i+11], 16) * 16 +
            int(hex_data[i+11:i+12], 16)
        )

        rows.append([fourth_nibble, eighth_nibble, tenth_nibble, str(value)])

    with open('another.csv', 'w') as file:
        for row in rows:
    # Convert each element in the row to a string before writing
            file.write(','.join(map(str, row)) + '\n')
            await asyncio.sleep(0)
except  Exception as e:
    print(f"An exception occurred: {e}")

[sophiedegran]

Hi
Have a look to github.com/peterhinch/micropython-async/blob/master/v3/docs/
They are lots of examples of async (look at auart.py)
Hope that help

[peterhinch]

Following on from @sophiedegran the way to asynchronously interface a UART is with StreamReader and StreamWriter instances, as in this sample. You can read about stream I/O in the tutorial.

[Chris-24MAY60]

Thanks, Peter, I imagine it's clear from reviewing my script, I've no experience with programming, much less the MicroPython language itself. The uart packet I am receiving is not of my own design, nor does it contain any newline characrers. I can say that a uart packet containing 438 bytes (always 438 bytes) represents a valid "occurrence", (I don't want to use the word event, because I think that has particular meaning to you in your parlance, and I don't understand, yet, what it means) ...while a packet of 6 bytes represents a "failed occurence". My users desrve to know if their attempt was unsuccessful. Can you tell me whether I must ask that these packet protocols need be revised, appending the newline character, so-that StreamReader can process them? And if the "line" is 438 hex bytes long, can StreamReader handle such length? Maybe before you reply, give me a day to see whether I can't answer these questions for myself. I do appreciate your help, and look forward to an educational day! Respectfully, Chris

…

On Sun, Jan 21, 2024, 5:48 AM Peter Hinch ***@***.***> wrote: Following on from @sophiedegran <https://github.com/sophiedegran> the way to asynchronously interface a UART is with StreamReader and StreamWriter instances, as in this sample <https://github.com/peterhinch/micropython-async/blob/master/v3/as_demos/auart.py>. You can read about stream I/O in the tutorial <https://github.com/peterhinch/micropython-async/blob/master/v3/docs/TUTORIAL.md#63-using-the-stream-mechanism> . — Reply to this email directly, view it on GitHub <#13479 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHS65MB3FHZ6QBVBNOJLYPT6BNAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DCOJXG43DS> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[Chris-24MAY60]

Hello Peter, After study, I still am not able to replace my blocking uart-read with a StreamReader version. import uasyncio as asyncio import asyncio import machine from machine import UART import time from time import sleep_ms import binascii from binascii import hexlify import network uart2 = UART(2, tx=17, rx=18) # 18:RX 17:TX uart2.init(baudrate=115200, bits=8, parity=None, stop=1) async def process_uart(): #code to process the file sreader = asyncio.StreamReader(uart2) while True: res = await sreader.readexactly(437) print("Received", res) buf = hexlify(res) print('Reporting', buf) hex_data = (buf[8:872]) print(hex_data) rows = [] for i in range(0, len(hex_data), 16): # 16 characters represent 16 nibbles in hex fourth_nibble = int(hex_data[i+3:i+4], 16) eighth_nibble = int(hex_data[i+7:i+8], 16) tenth_nibble = int(hex_data[i+9:i+10], 16) # Calculate the value in the fourth column as per the provided formula value = ( int(hex_data[i+14:i+15], 16) * 1048576 + int(hex_data[i+15:i+16], 16) * 65536 + int(hex_data[i+12:i+13], 16) * 4096 + int(hex_data[i+13:i+14], 16) * 256 + int(hex_data[i+10:i+11], 16) * 16 + int(hex_data[i+11:i+12], 16) ) rows.append([fourth_nibble, eighth_nibble, tenth_nibble, str(value)]) with open('putt_record.csv', 'w') as file: for row in rows: # Convert each element in the row to a string before writing file.write(','.join(map(str, row)) + '\n') await asyncio.sleep(0.025) async def main(): # Create tasks for asynchronous functions asyncio.create_task(process_uart()) while True: await asyncio.sleep(1) #Run the event loop asyncio.run(main()) I have a couple of questions. Firstly, even if I get this routine working -would this stand as a legitimate coroutine? I thought "with open(etc):" was a blocking operation. And secondly, where am I going wrong in deploying sreader.readexactly(437)? My other routine takes the uart packet and seemingly infallibly converts it to "another.csv" file. Every time. Clearly, the breakdown has to be my lack of understanding of how to deploy StreamReader. I've read the tutorial - and sadly - only gotten this far. Chris

…

On Sun, Jan 21, 2024 at 5:48 AM Peter Hinch ***@***.***> wrote: Following on from @sophiedegran <https://github.com/sophiedegran> the way to asynchronously interface a UART is with StreamReader and StreamWriter instances, as in this sample <https://github.com/peterhinch/micropython-async/blob/master/v3/as_demos/auart.py>. You can read about stream I/O in the tutorial <https://github.com/peterhinch/micropython-async/blob/master/v3/docs/TUTORIAL.md#63-using-the-stream-mechanism> . — Reply to this email directly, view it on GitHub <#13479 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHS65MB3FHZ6QBVBNOJLYPT6BNAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DCOJXG43DS> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[Chris-24MAY60]

I should have mentioned that the above works perfectly. Every once in a while. Occasionally, four or five "occurrences" in a row. But invariably, it wanders off into the weeds. My other routine never did that. On Sun, Jan 21, 2024 at 10:41 AM Christopher J Hahn ***@***.***> wrote:

…

Hello Peter, After study, I still am not able to replace my blocking uart-read with a StreamReader version. import uasyncio as asyncio import asyncio import machine from machine import UART import time from time import sleep_ms import binascii from binascii import hexlify import network uart2 = UART(2, tx=17, rx=18) # 18:RX 17:TX uart2.init(baudrate=115200, bits=8, parity=None, stop=1) async def process_uart(): #code to process the file sreader = asyncio.StreamReader(uart2) while True: res = await sreader.readexactly(437) print("Received", res) buf = hexlify(res) print('Reporting', buf) hex_data = (buf[8:872]) print(hex_data) rows = [] for i in range(0, len(hex_data), 16): # 16 characters represent 16 nibbles in hex fourth_nibble = int(hex_data[i+3:i+4], 16) eighth_nibble = int(hex_data[i+7:i+8], 16) tenth_nibble = int(hex_data[i+9:i+10], 16) # Calculate the value in the fourth column as per the provided formula value = ( int(hex_data[i+14:i+15], 16) * 1048576 + int(hex_data[i+15:i+16], 16) * 65536 + int(hex_data[i+12:i+13], 16) * 4096 + int(hex_data[i+13:i+14], 16) * 256 + int(hex_data[i+10:i+11], 16) * 16 + int(hex_data[i+11:i+12], 16) ) rows.append([fourth_nibble, eighth_nibble, tenth_nibble, str(value)]) with open('putt_record.csv', 'w') as file: for row in rows: # Convert each element in the row to a string before writing file.write(','.join(map(str, row)) + '\n') await asyncio.sleep(0.025) async def main(): # Create tasks for asynchronous functions asyncio.create_task(process_uart()) while True: await asyncio.sleep(1) #Run the event loop asyncio.run(main()) I have a couple of questions. Firstly, even if I get this routine working -would this stand as a legitimate coroutine? I thought "with open(etc):" was a blocking operation. And secondly, where am I going wrong in deploying sreader.readexactly(437)? My other routine takes the uart packet and seemingly infallibly converts it to "another.csv" file. Every time. Clearly, the breakdown has to be my lack of understanding of how to deploy StreamReader. I've read the tutorial - and sadly - only gotten this far. Chris On Sun, Jan 21, 2024 at 5:48 AM Peter Hinch ***@***.***> wrote: > Following on from @sophiedegran <https://github.com/sophiedegran> the > way to asynchronously interface a UART is with StreamReader and > StreamWriter instances, as in this sample > <https://github.com/peterhinch/micropython-async/blob/master/v3/as_demos/auart.py>. > You can read about stream I/O in the tutorial > <https://github.com/peterhinch/micropython-async/blob/master/v3/docs/TUTORIAL.md#63-using-the-stream-mechanism> > . > > — > Reply to this email directly, view it on GitHub > <#13479 (comment)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AF7CHS65MB3FHZ6QBVBNOJLYPT6BNAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DCOJXG43DS> > . > You are receiving this because you authored the thread.Message ID: > ***@***.*** > com> >

[Chris-24MAY60]

Don't help me. I'm making real headway.

…

On Sun, Jan 21, 2024 at 5:48 AM Peter Hinch ***@***.***> wrote: Following on from @sophiedegran <https://github.com/sophiedegran> the way to asynchronously interface a UART is with StreamReader and StreamWriter instances, as in this sample <https://github.com/peterhinch/micropython-async/blob/master/v3/as_demos/auart.py>. You can read about stream I/O in the tutorial <https://github.com/peterhinch/micropython-async/blob/master/v3/docs/TUTORIAL.md#63-using-the-stream-mechanism> . — Reply to this email directly, view it on GitHub <#13479 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHS65MB3FHZ6QBVBNOJLYPT6BNAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DCOJXG43DS> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[Chris-24MAY60]

Everything now working swimmingly - and (pending testing) seemingly - all the time. Appreciate your looking in on me. Thanks for guidance. Chris

…

On Sun, Jan 21, 2024 at 5:48 AM Peter Hinch ***@***.***> wrote: Following on from @sophiedegran <https://github.com/sophiedegran> the way to asynchronously interface a UART is with StreamReader and StreamWriter instances, as in this sample <https://github.com/peterhinch/micropython-async/blob/master/v3/as_demos/auart.py>. You can read about stream I/O in the tutorial <https://github.com/peterhinch/micropython-async/blob/master/v3/docs/TUTORIAL.md#63-using-the-stream-mechanism> . — Reply to this email directly, view it on GitHub <#13479 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHS65MB3FHZ6QBVBNOJLYPT6BNAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DCOJXG43DS> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[sophiedegran]

Nice Chris I was sure you could do it ;)

[Chris-24MAY60]

Next up, third-order and fourth-order polynomial curve-fit regressions on, say, 20 data-points. Jim Mussared has suggested pursuing this with the ulab guys. Exciting! Thanks for direction, Sophie. Chris

…

On Sun, Jan 21, 2024, 2:20 PM sophiedegran ***@***.***> wrote: Nice Chris I was sure you could do it ;) — Reply to this email directly, view it on GitHub <#13479 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHS3CEFRKSZRJXCHK2RLYPV2ALAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMBRGQ2TS> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[peterhinch]

Congratulations on getting it working.

It does seem as if a binary value is being converted to a hexadecimal string then being converted to integers. If I'm reading this correctly, there is some scope for cutting out the middleman...

[Chris-24MAY60]

Thanks for (continued) guidance. I have no doubt you're correct. https://www.poetryfoundation.org/poems/42891/stopping-by-woods-on-a-snowy-evening

…

On Mon, Jan 22, 2024, 8:27 AM Peter Hinch ***@***.***> wrote: Congratulations on getting it working. It does seem as if a binary value is being converted to a hexadecimal string then being converted to integers. If I'm reading this correctly, there is some scope for cutting out the middleman... — Reply to this email directly, view it on GitHub <#13479 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHSYQJECXWC632ETSA4DYPZZLLAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMBZGM2TM> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[peterhinch]

Hint:

>>> a = b'\x01\x02\x03\x04'
>>> m = memoryview(a)
>>> int.from_bytes(m[1:3], 'little')
770
>>> 0x302
770
>>> 

The use of a memoryview isn't essential:

>>> int.from_bytes(a[1:3], 'little')
770

but the memoryview saves needless allocation.

[Chris-24MAY60]

Peter, Thank you, again. My 432 bytes (after I strip 4 from the front and 2 from the back ...of my 438 byte read) comprises 54 records, each of 8 bytes. My "middleman", if I understand you correctly, heralds from a "little-endian" presentation of bytes 6, 7, and 8 ...which I wish were reversed, so I could simply concatenate them and convert from base 16 to base 10. I am very sorry, but am not following your hint. Compounded because, clumsy (not elegant) as my method is ...it works and I am driving to realize a rather (at least for me) complicated product - and I tend to accept "works" where "could use further refinement" is almost always true. Today I am studying how to brute force a 4th order polynomial regression of 13 observations. If you might connect me with someone who could guide me in those efforts, I would, again, be ever so grateful. And further, I might have as few as 4, or even 3 observations, reducing the order of the polynomial I can form. I am running on an ESP32S3 and have only a few different computations that need to be run over and again. So I don't want a big math library. I want a "only does what I need" library". Efficiency *does* matter to me. But in the name of continuing to make progress on my product, I might choose to delay pursuing your hint and I hope you will forgive me. Chris

…

On Tue, Jan 23, 2024, 7:42 AM Peter Hinch ***@***.***> wrote: Hint: >>> a = b'\x01\x02\x03\x04'>>> m = memoryview(a)>>> int.from_bytes(m[1:3], 'little')770>>> 0x302770>>> The use of a memoryview isn't essential: >>> int.from_bytes(a[1:3], 'little')770 but the memoryview saves needless allocation. — Reply to this email directly, view it on GitHub <#13479 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHSYW3KNMMLKTGFB6A7LYP642VAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMRQGY4DS> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[Chris-24MAY60]

Peter, I see now your hint prepares me, looking forward. I am notoriously slow, however, pride myself by compensating with low comprehension! I tend to proceed one task at a time - but am caught-up imagining I could avail myself to an immediate better result - if I were to lay-out the task, end-to-end for your consideration. Were I to do that, in general terms, could you share with me how you might go about it? It could possibly lend insight(s), for your forum's readers - to learn how you approach multi-step tasks. And selfishly, it would help me. Chris On Tue, Jan 23, 2024 at 8:44 AM Christopher J Hahn ***@***.***> wrote:

…

Peter, Thank you, again. My 432 bytes (after I strip 4 from the front and 2 from the back ...of my 438 byte read) comprises 54 records, each of 8 bytes. My "middleman", if I understand you correctly, heralds from a "little-endian" presentation of bytes 6, 7, and 8 ...which I wish were reversed, so I could simply concatenate them and convert from base 16 to base 10. I am very sorry, but am not following your hint. Compounded because, clumsy (not elegant) as my method is ...it works and I am driving to realize a rather (at least for me) complicated product - and I tend to accept "works" where "could use further refinement" is almost always true. Today I am studying how to brute force a 4th order polynomial regression of 13 observations. If you might connect me with someone who could guide me in those efforts, I would, again, be ever so grateful. And further, I might have as few as 4, or even 3 observations, reducing the order of the polynomial I can form. I am running on an ESP32S3 and have only a few different computations that need to be run over and again. So I don't want a big math library. I want a "only does what I need" library". Efficiency *does* matter to me. But in the name of continuing to make progress on my product, I might choose to delay pursuing your hint and I hope you will forgive me. Chris On Tue, Jan 23, 2024, 7:42 AM Peter Hinch ***@***.***> wrote: > Hint: > > >>> a = b'\x01\x02\x03\x04'>>> m = memoryview(a)>>> int.from_bytes(m[1:3], 'little')770>>> 0x302770>>> > > The use of a memoryview isn't essential: > > >>> int.from_bytes(a[1:3], 'little')770 > > but the memoryview saves needless allocation. > > — > Reply to this email directly, view it on GitHub > <#13479 (reply in thread)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AF7CHSYW3KNMMLKTGFB6A7LYP642VAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMRQGY4DS> > . > You are receiving this because you authored the thread.Message ID: > ***@***.*** > com> >

[peterhinch]

Feel free to describe the task. I may or may not be able to help, but there are plenty of others here, some of whom are very skilled at mathematics.

[Chris-24MAY60]

Thanks, again, Peter. I have just completed my own review of linear algebra. If you can tell me whether MicroPython offers: matrix_multiplication, matrix_transpose, and matrix_inverse, I think I can take it from there. Perhaps, if m-observations, nth-order polynomial regression doesn't already exist for MicroPython, I might be able to make a contribution! That would sure make me feel good. Chris

…

On Tue, Jan 23, 2024 at 12:09 PM Peter Hinch ***@***.***> wrote: Feel free to describe the task. I may or may not be able to help, but there are plenty of others here, some of whom are very skilled at mathematics. — Reply to this email directly, view it on GitHub <#13479 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHS5CVY6X7GW2C7AXBADYP74GVAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMRUGEYTM> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[Chris-24MAY60]

Peter! Never mind. (again) I found it! https://github.com/iyassou/umatrix This library provides my three asked-for matrix-operation capabilities. Replete with examples. What an incredible body of work that's been done with MicroPython. Chris On Tue, Jan 23, 2024 at 12:18 PM Christopher J Hahn ***@***.***> wrote:

…

Thanks, again, Peter. I have just completed my own review of linear algebra. If you can tell me whether MicroPython offers: matrix_multiplication, matrix_transpose, and matrix_inverse, I think I can take it from there. Perhaps, if m-observations, nth-order polynomial regression doesn't already exist for MicroPython, I might be able to make a contribution! That would sure make me feel good. Chris On Tue, Jan 23, 2024 at 12:09 PM Peter Hinch ***@***.***> wrote: > Feel free to describe the task. I may or may not be able to help, but > there are plenty of others here, some of whom are very skilled at > mathematics. > > — > Reply to this email directly, view it on GitHub > <#13479 (reply in thread)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AF7CHS5CVY6X7GW2C7AXBADYP74GVAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMRUGEYTM> > . > You are receiving this because you authored the thread.Message ID: > ***@***.*** > com> >

[rkompass]

After some reading into this it appears that the Wikipedia article is a nice description of the polynomial regression you are trying to code.
The last formula therein:
beta = (X^T X)⁻1 X^T y
for the polynomial coefficients beta is literally solved as it is written there both in ulab and in numpy (I found after checking the sources).
I.e. the matrix X^T X (the transpose of the Vandermonde matrix multiplied with itself, 5x5 in your case for 4-th order polynomial) is inverted.
But this gives away the advantage that this matrix is positive definite.
The preferred method should be to do a Cholesky decomposition of that matrix and solve X^T y = (X^T X) beta using that. Even better using the related LDL decomposition, which avoids computations of square roots.
In theory this should be both more precise and faster.
In Cholesky decomposition there is also the code for the cholesky decomposition. I can help to transform that into micropython code.
Btw. https://github.com/iyassou/umatrix has inverse only up to 4x4 matrix size, so it is not usable directly.

Would be interesting to see if this approach has advantages over the one in numpy and ulab.

[Chris-24MAY60]

Nope, I have to feed the CD a square matrix. Let me chew on it some more. Do you know of an online CD calculator that'll go to 5 x 5? I haven't found one. Thanks, Chris On Wed, Jan 24, 2024 at 9:37 AM Christopher J Hahn ***@***.***> wrote:

…

Hi Raul, So I've been pouring over my OLS (ordinary least squares) method - in an effort to understand the Cholesky decomposition. The first thing I did was multiply (X^T)(y), only to learn it yielded a five-row, one-column matrix. So next, I multiplied (X^T*X)(beta) and (no surprise to you, I imagine) got the same five-row, one-column matrix, assuring me that both sides of the Cholesky decomposition equation you gave me was, in fact true. Today, I'm going to attempt your suggested magic trick! Am I to understand that the Cholesky decomposition will operate on the "product-matrix", being only five rows by one column - until out falls two matrices, one of which so-happens to be the (5x5) (X^T*X)^-1 TOGETHER WITH the (5x1) Beta matrix! of polynomial coefficients, a4, a3, a2, a1, and a0?! I have to believe I have this somehow very wrong!! But I'm going to attempt firing-up Python this morning and running the Cholesky decomposition routine - to which you pointed me in you previous e-mail - and seeing what happens. If I AM way off base, please do drop me a quick line, doing what you can to correct my misconception(s), thereby saving me lost time. If I am understanding things - and this works - you'll have to tell me where to mail the check! ; ) Respectfully, Chris On Tue, Jan 23, 2024 at 5:36 PM Raul Kompaß ***@***.***> wrote: > You only need 5x5 matrices. > beta is the vector of polynomial coefficients. beta = (beta_0, beta_1, > .., beta_4)^T > with the polynomial being f(x)=beta_0 + beta_1 * x + beta_2 * x^2 + .. + > beta_4 * x^4 . > The 5x5 matrix (X^T X) should be computable in ridiculously simple > manner, if I'm not mistaken: > (X^T X)_ij = sum_k (x_k^(i+j)). Summation over all 13 x_k. > Cholesky decomposition code for python is here > <https://rosettacode.org/wiki/Cholesky_decomposition#Python>. Rosetta > code is our friend. > > Do you have 13 pairs of x and y values at hand? It's better to have > something to look at, and test with, while coding. > > — > Reply to this email directly, view it on GitHub > <#13479 (reply in thread)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AF7CHS3V3C7SN2WVPMTG4PTYQBCQNAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMRWGUZDC> > . > You are receiving this because you authored the thread.Message ID: > ***@***.*** > com> >

[rkompass]

Am I to understand that the
Cholesky decomposition will operate on the "product-matrix", being only
five rows by one column - until out falls two matrices, one of which
so-happens to be the (5x5) (X^T*X)^-1 TOGETHER WITH the (5x1) Beta matrix!

Almost right: the (X^T*X) is a 5x5 matrix.

from math import sqrt

def pprint(A):
    m, n = len(A), len(A[0])
    for i, r in enumerate(A):
        print(' [' if i else '[[', end='')
        for j, e in enumerate(r):
            print('{:6.2f}{}'.format(e, ',' if j < n-1 else ''), end='')           
        print('],' if i < m-1 else ']]')

def zeros(m, n):
    return [[0 for i in range(n)] for j in range(m)]

def matmul(A, B):  # https://rosettacode.org/wiki/Matrix_multiplication#Python
    return [[sum(x * B[i][col] for i,x in enumerate(row)) for col in range(len(B[0]))] for row in A]

def transpose(A):
    return list(zip(*A))

def cholesky(A):   # https://rosettacode.org/wiki/Cholesky_decomposition#Python
    L = zeros(len(A), len(A))
    for i, (Ai, Li) in enumerate(zip(A, L)):
        for j, Lj in enumerate(L[:i+1]):
            s = sum(Li[k]*Lj[k] for k in range(j))
            Li[j] = sqrt(Ai[i]-s) if (i==j) else (Ai[j]-s)/Lj[j]
    return L

x = [0, 0.065778, 0.121406, 0.165108, 0.2023745, 0.7878015, 0.8084525,
     0.830699, 0.8519465, 0.873236, 0.8948235, 0.917743, 0.9424455]
y = [0, -1.25, -2.5, -3.75, -5, -3.75, -2.5, -1.25, 0, 1.25, 2.5, 3.75, 5]

n = 4     # polynomial degree

def vandermonde(x, d):
    return [[e**j for j in range(d)] for e in x]

def vandermonde_xtx(x, d):   # d will be 5, it is the dimension of this matrix, i.e. n+1
    X = zeros(d, d)
    for j in range(d):
        X[0][j] = sum(el**j for el in x)
    for i in range(1, d):
        X[i][:-1] = X[i-1][1:]
        X[i][d-1] = sum(el**(i+d-1) for el in x)
    return X

v = vandermonde(x, n+1)
pprint(v); print()
vtv = matmul(transpose(v), v)
pprint(vtv); print()

xtx = vandermonde_xtx(x, n+1)
xty = [sum(x[i]**j * y[i] for i in range(len(x))) for j in range(n+1)]

pprint(xtx); print()
print(xty); print()
l = cholesky(xtx)
pprint(l); print()

returns here:

[[  1.00,  0.00,  0.00,  0.00,  0.00],
 [  1.00,  0.07,  0.00,  0.00,  0.00],
 [  1.00,  0.12,  0.01,  0.00,  0.00],
 [  1.00,  0.17,  0.03,  0.00,  0.00],
 [  1.00,  0.20,  0.04,  0.01,  0.00],
 [  1.00,  0.79,  0.62,  0.49,  0.39],
 [  1.00,  0.81,  0.65,  0.53,  0.43],
 [  1.00,  0.83,  0.69,  0.57,  0.48],
 [  1.00,  0.85,  0.73,  0.62,  0.53],
 [  1.00,  0.87,  0.76,  0.67,  0.58],
 [  1.00,  0.89,  0.80,  0.72,  0.64],
 [  1.00,  0.92,  0.84,  0.77,  0.71],
 [  1.00,  0.94,  0.89,  0.84,  0.79]]

[[ 13.00,  7.46,  6.07,  5.22,  4.54],
 [  7.46,  6.07,  5.22,  4.54,  3.97],
 [  6.07,  5.22,  4.54,  3.97,  3.48],
 [  5.22,  4.54,  3.97,  3.48,  3.07],
 [  4.54,  3.97,  3.48,  3.07,  2.71]]

[[ 13.00,  7.46,  6.07,  5.22,  4.54],
 [  7.46,  6.07,  5.22,  4.54,  3.97],
 [  6.07,  5.22,  4.54,  3.97,  3.48],
 [  5.22,  4.54,  3.97,  3.48,  3.07],
 [  4.54,  3.97,  3.48,  3.07,  2.71]]

[-7.5, 3.451842, 5.381219, 5.773444, 5.815036]

[[  3.61,  0.00,  0.00,  0.00,  0.00],
 [  2.07,  1.34,  0.00,  0.00,  0.00],
 [  1.68,  1.29,  0.16,  0.00,  0.00],
 [  1.45,  1.16,  0.23,  0.06,  0.00],
 [  1.26,  1.02,  0.26,  0.12,  0.01]]

Now what's still needed is forward substitution (and perhaps backward substitution, but probably not).
Which I have found already but have to check/improve yet.

Then the solution of xty = xtx * beta is there.

A lot of cross checking will be necessary too. I included computing the raw vandermonde matrix and multiplying it with it's transpose.
You see the result is the same.
Doing LU decomposition and solving the equation by that as an alternative - to be done yet.
But you see: There is not much to compute, once it is clear, how it should be done.

[Chris-24MAY60]

Hello, I've gotten my app working - but decided to preclude data from being taken (and sent) by my "measuring_node" ESP32S3 until the WiFi connection is made on the "reporting_node" ESP32S3 - the one running MicroPython. Immediately after the WiFi connection is made, I am trying to send a message, byte-wise, over the uart, to the "measuring_node" ESP32S3 that it can now spring into action, beginning to collect (and send) data. while not ap.isconnected(): await asyncio.sleep(1) # Add a short delay to avoid busy waiting swriter = asyncio.StreamWriter(uart2, {}) swriter.write(b'/x01/x00/x03/x62') await swriter.drain() # Transmission starts now. await asyncio.sleep(0.5) print('Connection successful') print(ap.ifconfig()) But my "measuring_node" ESP32S3 never comes out of initialization. I've reviewed Peter Hinch's tutorial as well as the MicroPython StreamReader documentation - but can't find an example of how to send a uart message byte-wise. Seems this should be easy. Any help would be appreciated. Chris

…

On Tue, Jan 23, 2024 at 12:09 PM Peter Hinch ***@***.***> wrote: Feel free to describe the task. I may or may not be able to help, but there are plenty of others here, some of whom are very skilled at mathematics. — Reply to this email directly, view it on GitHub <#13479 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHS5CVY6X7GW2C7AXBADYP74GVAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMRUGEYTM> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[Chris-24MAY60]

I'm sorry, never mind! My slashes were backwards and I was sending ascii and not bytes! Rookie mistake. On Wed, Jan 24, 2024 at 2:21 PM Christopher J Hahn ***@***.***> wrote:

…

Hello, I've gotten my app working - but decided to preclude data from being taken (and sent) by my "measuring_node" ESP32S3 until the WiFi connection is made on the "reporting_node" ESP32S3 - the one running MicroPython. Immediately after the WiFi connection is made, I am trying to send a message, byte-wise, over the uart, to the "measuring_node" ESP32S3 that it can now spring into action, beginning to collect (and send) data. while not ap.isconnected(): await asyncio.sleep(1) # Add a short delay to avoid busy waiting swriter = asyncio.StreamWriter(uart2, {}) swriter.write(b'/x01/x00/x03/x62') await swriter.drain() # Transmission starts now. await asyncio.sleep(0.5) print('Connection successful') print(ap.ifconfig()) But my "measuring_node" ESP32S3 never comes out of initialization. I've reviewed Peter Hinch's tutorial as well as the MicroPython StreamReader documentation - but can't find an example of how to send a uart message byte-wise. Seems this should be easy. Any help would be appreciated. Chris On Tue, Jan 23, 2024 at 12:09 PM Peter Hinch ***@***.***> wrote: > Feel free to describe the task. I may or may not be able to help, but > there are plenty of others here, some of whom are very skilled at > mathematics. > > — > Reply to this email directly, view it on GitHub > <#13479 (reply in thread)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/AF7CHS5CVY6X7GW2C7AXBADYP74GVAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMRUGEYTM> > . > You are receiving this because you authored the thread.Message ID: > ***@***.*** > com> >

[Chris-24MAY60]

Thanks, Raul, I got tied up today but now my plate is cleared. I will start with this tomorrow morning, sticking with it until I understand it. Thank you, very much, for putting this together. Chris

…

On Wed, Jan 24, 2024, 11:17 AM Raul Kompaß ***@***.***> wrote: Am I to understand that the Cholesky decomposition will operate on the "product-matrix", being only five rows by one column - until out falls two matrices, one of which so-happens to be the (5x5) (X^T*X)^-1 TOGETHER WITH the (5x1) Beta matrix! Almost right: the (X^T*X) is a 5x5 matrix. from math import sqrt def pprint(A): m, n = len(A), len(A[0]) for i, r in enumerate(A): print(' [' if i else '[[', end='') for j, e in enumerate(r): print('{:6.2f}{}'.format(e, ',' if j < n-1 else ''), end='') print('],' if i < m-1 else ']]') def zeros(m, n): return [[0 for i in range(n)] for j in range(m)] def matmul(A, B): # https://rosettacode.org/wiki/Matrix_multiplication#Python return [[sum(x * B[i][col] for i,x in enumerate(row)) for col in range(len(B[0]))] for row in A] def transpose(A): return list(zip(*A)) def cholesky(A): # https://rosettacode.org/wiki/Cholesky_decomposition#Python L = zeros(len(A), len(A)) for i, (Ai, Li) in enumerate(zip(A, L)): for j, Lj in enumerate(L[:i+1]): s = sum(Li[k]*Lj[k] for k in range(j)) Li[j] = sqrt(Ai[i]-s) if (i==j) else (Ai[j]-s)/Lj[j] return L x = [0, 0.065778, 0.121406, 0.165108, 0.2023745, 0.7878015, 0.8084525, 0.830699, 0.8519465, 0.873236, 0.8948235, 0.917743, 0.9424455]y = [0, -1.25, -2.5, -3.75, -5, -3.75, -2.5, -1.25, 0, 1.25, 2.5, 3.75, 5] n = 4 # polynomial degree def vandermonde(x, d): return [[e**j for j in range(d)] for e in x] def vandermonde_xtx(x, d): # d will be 5, it is the dimension of this matrix, i.e. n+1 X = zeros(d, d) for j in range(d): X[0][j] = sum(el**j for el in x) for i in range(1, d): X[i][:-1] = X[i-1][1:] X[i][d-1] = sum(el**(i+d-1) for el in x) return X v = vandermonde(x, n+1)pprint(v); print()vtv = matmul(transpose(v), v)pprint(vtv); print() xtx = vandermonde_xtx(x, n+1)xty = [sum(x[i]**j * y[i] for i in range(len(x))) for j in range(n+1)] pprint(xtx); print()print(xty); print()l = cholesky(xtx)pprint(l); print() returns here: [[ 1.00, 0.00, 0.00, 0.00, 0.00], [ 1.00, 0.07, 0.00, 0.00, 0.00], [ 1.00, 0.12, 0.01, 0.00, 0.00], [ 1.00, 0.17, 0.03, 0.00, 0.00], [ 1.00, 0.20, 0.04, 0.01, 0.00], [ 1.00, 0.79, 0.62, 0.49, 0.39], [ 1.00, 0.81, 0.65, 0.53, 0.43], [ 1.00, 0.83, 0.69, 0.57, 0.48], [ 1.00, 0.85, 0.73, 0.62, 0.53], [ 1.00, 0.87, 0.76, 0.67, 0.58], [ 1.00, 0.89, 0.80, 0.72, 0.64], [ 1.00, 0.92, 0.84, 0.77, 0.71], [ 1.00, 0.94, 0.89, 0.84, 0.79]] [[ 13.00, 7.46, 6.07, 5.22, 4.54], [ 7.46, 6.07, 5.22, 4.54, 3.97], [ 6.07, 5.22, 4.54, 3.97, 3.48], [ 5.22, 4.54, 3.97, 3.48, 3.07], [ 4.54, 3.97, 3.48, 3.07, 2.71]] [[ 13.00, 7.46, 6.07, 5.22, 4.54], [ 7.46, 6.07, 5.22, 4.54, 3.97], [ 6.07, 5.22, 4.54, 3.97, 3.48], [ 5.22, 4.54, 3.97, 3.48, 3.07], [ 4.54, 3.97, 3.48, 3.07, 2.71]] [-7.5, 3.451842, 5.381219, 5.773444, 5.815036] [[ 3.61, 0.00, 0.00, 0.00, 0.00], [ 2.07, 1.34, 0.00, 0.00, 0.00], [ 1.68, 1.29, 0.16, 0.00, 0.00], [ 1.45, 1.16, 0.23, 0.06, 0.00], [ 1.26, 1.02, 0.26, 0.12, 0.01]] Now what's still needed is forward substitution (and perhaps backward substitution, but probably not). Which I have found already but have to check/improve yet. Then the solution of xty = xtx * beta is there. A lot of cross checking will be necessary too. I included computing the raw vandermonde matrix and multiplying it with it's transpose. You see the result is the same. Doing LU decomposition and solving the equation by that as an alternative - to be done yet. But you see: There is not much to compute, once it is clear, how it should be done. — Reply to this email directly, view it on GitHub <#13479 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHSYSZZ6NY5SRRRVUAWTYQE6ZBAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMZVGY3TK> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[rkompass]

Strong improvements (~8..12) in interpolation mean square error up to polynomial degree 5.
The assumption was right:
Cholesky solve is better, but it crashes (tries to get the square root from negative value) for higher degree than 6.
But with ordinary (LU decomposition) solving the polynomial doesn't get better results anyway.

# poly degree, mean square error, cholesky solve
# 2 0.6066617
# 3 0.07040891
# 4 0.01398384
# 5 0.001057713
# 6 0.000780353
    
# poly degree, mean square error, ordinary solve
# 2 0.6066621
# 3 0.07040879
# 4 0.01398485
# 5 0.003148756
# 6 0.001345163
# 7 0.001312334
# 8 0.001207661
# 9 0.002023049

[Chris-24MAY60]

Hi Raul, Confess to being a bit lost. I don't understand how the last two matrices help us determine beta. I follow everything you've done except both the generation of the Cholesky (X^T*X) ...and its usefulness in helping us find beta. I used the OLS method and got these coefficients: a0 0.210950563440863 a1 -22.8906621119621 a2 -20.8906997505424 a3 67.6304434918384 a4 -14.1119736691275 Chris

…

On Wed, Jan 24, 2024 at 11:17 AM Raul Kompaß ***@***.***> wrote: Am I to understand that the Cholesky decomposition will operate on the "product-matrix", being only five rows by one column - until out falls two matrices, one of which so-happens to be the (5x5) (X^T*X)^-1 TOGETHER WITH the (5x1) Beta matrix! Almost right: the (X^T*X) is a 5x5 matrix. from math import sqrt def pprint(A): m, n = len(A), len(A[0]) for i, r in enumerate(A): print(' [' if i else '[[', end='') for j, e in enumerate(r): print('{:6.2f}{}'.format(e, ',' if j < n-1 else ''), end='') print('],' if i < m-1 else ']]') def zeros(m, n): return [[0 for i in range(n)] for j in range(m)] def matmul(A, B): # https://rosettacode.org/wiki/Matrix_multiplication#Python return [[sum(x * B[i][col] for i,x in enumerate(row)) for col in range(len(B[0]))] for row in A] def transpose(A): return list(zip(*A)) def cholesky(A): # https://rosettacode.org/wiki/Cholesky_decomposition#Python L = zeros(len(A), len(A)) for i, (Ai, Li) in enumerate(zip(A, L)): for j, Lj in enumerate(L[:i+1]): s = sum(Li[k]*Lj[k] for k in range(j)) Li[j] = sqrt(Ai[i]-s) if (i==j) else (Ai[j]-s)/Lj[j] return L x = [0, 0.065778, 0.121406, 0.165108, 0.2023745, 0.7878015, 0.8084525, 0.830699, 0.8519465, 0.873236, 0.8948235, 0.917743, 0.9424455]y = [0, -1.25, -2.5, -3.75, -5, -3.75, -2.5, -1.25, 0, 1.25, 2.5, 3.75, 5] n = 4 # polynomial degree def vandermonde(x, d): return [[e**j for j in range(d)] for e in x] def vandermonde_xtx(x, d): # d will be 5, it is the dimension of this matrix, i.e. n+1 X = zeros(d, d) for j in range(d): X[0][j] = sum(el**j for el in x) for i in range(1, d): X[i][:-1] = X[i-1][1:] X[i][d-1] = sum(el**(i+d-1) for el in x) return X v = vandermonde(x, n+1)pprint(v); print()vtv = matmul(transpose(v), v)pprint(vtv); print() xtx = vandermonde_xtx(x, n+1)xty = [sum(x[i]**j * y[i] for i in range(len(x))) for j in range(n+1)] pprint(xtx); print()print(xty); print()l = cholesky(xtx)pprint(l); print() returns here: [[ 1.00, 0.00, 0.00, 0.00, 0.00], [ 1.00, 0.07, 0.00, 0.00, 0.00], [ 1.00, 0.12, 0.01, 0.00, 0.00], [ 1.00, 0.17, 0.03, 0.00, 0.00], [ 1.00, 0.20, 0.04, 0.01, 0.00], [ 1.00, 0.79, 0.62, 0.49, 0.39], [ 1.00, 0.81, 0.65, 0.53, 0.43], [ 1.00, 0.83, 0.69, 0.57, 0.48], [ 1.00, 0.85, 0.73, 0.62, 0.53], [ 1.00, 0.87, 0.76, 0.67, 0.58], [ 1.00, 0.89, 0.80, 0.72, 0.64], [ 1.00, 0.92, 0.84, 0.77, 0.71], [ 1.00, 0.94, 0.89, 0.84, 0.79]] [[ 13.00, 7.46, 6.07, 5.22, 4.54], [ 7.46, 6.07, 5.22, 4.54, 3.97], [ 6.07, 5.22, 4.54, 3.97, 3.48], [ 5.22, 4.54, 3.97, 3.48, 3.07], [ 4.54, 3.97, 3.48, 3.07, 2.71]] [[ 13.00, 7.46, 6.07, 5.22, 4.54], [ 7.46, 6.07, 5.22, 4.54, 3.97], [ 6.07, 5.22, 4.54, 3.97, 3.48], [ 5.22, 4.54, 3.97, 3.48, 3.07], [ 4.54, 3.97, 3.48, 3.07, 2.71]] [-7.5, 3.451842, 5.381219, 5.773444, 5.815036] [[ 3.61, 0.00, 0.00, 0.00, 0.00], [ 2.07, 1.34, 0.00, 0.00, 0.00], [ 1.68, 1.29, 0.16, 0.00, 0.00], [ 1.45, 1.16, 0.23, 0.06, 0.00], [ 1.26, 1.02, 0.26, 0.12, 0.01]] Now what's still needed is forward substitution (and perhaps backward substitution, but probably not). Which I have found already but have to check/improve yet. Then the solution of xty = xtx * beta is there. A lot of cross checking will be necessary too. I included computing the raw vandermonde matrix and multiplying it with it's transpose. You see the result is the same. Doing LU decomposition and solving the equation by that as an alternative - to be done yet. But you see: There is not much to compute, once it is clear, how it should be done. — Reply to this email directly, view it on GitHub <#13479 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHSYSZZ6NY5SRRRVUAWTYQE6ZBAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMZVGY3TK> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[rkompass]

O.K.: In the Wikipdia article it is stated that the Polynomial coefficients, represented as a vector may be computed by the formula there (the last formula).

_beta = (VtV)⁻1 * Vty.

Where V is the Vandermonde determinant. And VtV the product of the transpose of V with V.
I followed the article to call it XtX, this may be confusing.

Now, to compute VtV it is not necessary to have V, you may shorten the calculation and don't need a matrix product.

def vandermonde(x, d):   # up to x^d  dimension will be (len(x), d+1) 
    return [[e**j for j in range(d+1)] for e in x]

def vandermonde_VtV(x, d):   # dimension will be (d+1, d+1)
    X = zeros(d+1, d+1)
    for j in range(d+1):
        X[0][j] = sum(el**j for el in x)
    for i in range(1, d+1):
        X[i][:-1] = X[i-1][1:]
        X[i][d] = sum(el**(i+d) for el in x)
    return X

The equality is demonstrated in the above codes results.

In order to compute (VtV)⁻1 * Vty (wich is 5x5 matrix times 5-vector) you need the inverse of VtV.
Instead of computing the inverse of VtV you may but shorten again and solve

VtV * _beta = Vty

This solves 1 vector, while computing the inverse has to solve for n vectors.
My code for this computation is here (hopefully complete).

def lu_decomp(A):  # https://rosettacode.org/wiki/LU_decomposition#Python
    """Decomposes a square matrix A by perm_mat(P) A = L U and returns L, U and P."""
    n = len(A)
    L, U, P = eye(n), zeros(n, n), [i for i in range(n)]
    for j in range(n):
        row = max(range(j, n), key=lambda i: abs(A[i][j]))
        if j != row:
            P[j], P[row] = P[row], P[j]
    A2 = [A[i] for i in P]             # A2 = matmul(perm_mat(P), A)
    for j in range(n):
        for i in range(j+1):
            U[i][j] = A2[i][j] - sum(U[k][j]*L[i][k] for k in range(i))
        for i in range(j, n):
            L[i][j] = (A2[i][j] - sum(U[k][j]*L[i][k] for k in range(j))) / U[j][j]
    return L, U, P

def forw_sub(L, b):  # Solve L*y = b. L is assumed to be lower triangular matrix
    d = [0]*len(b)
    d[0] = b[0]/L[0][0]
    for i in range(1, len(b)):
        d[i] = (b[i] - sum(d[j] * L[i][j] for j in range(i))) / L[i][i]
    return d

def backw_sub(U, y):  # Solve U*x = y. U is assumed to be upper triangular matrix
    n = len(y)
    x = [0]*n
    x[-1] = y[-1] / U[-1][-1]
    for i in range(n-2, -1, -1):
        x[i] = (y[i] - sum(x[j] * U[i][j] for j in range(i+1, n))) / U[i][i]
    return x

def lu_solve(L, U, b):  # solve L*U*x = b with lower and upper triangular matrices L and U
    y = forw_sub(L, b)
    return backw_sub(U, y)

def solve(A, b):        # solve A*x = b for x, A is arbitrary (hopefully well conditioned) square matrix
    L, U, P = lu_decomp(A)
    pb = [b[j] for j in P]
    return lu_solve(L, U, pb)

def solve_chol(A, b):  # Alternative to solve() if A is symmetric and positive definite
    L = cholesky(A)
    return lu_solve(L, transp(L), b)

solve and solve_chol differ in that instead of using the L and U matrices of an LU decomposition for the partial solutions using forward and backward substitution (which is quite equivalent to Gaussian algorithm) the Cholesky matrix and its transpose are used.

All this is Ordinary Least Squares, of course, but I'm not familiar now, how the Wikipedia formula derives from the OLS principle (may look that up).
But you might of course also do a gradient descent of the polynomial parameters. Did you do that? Should be also nice.

With the combined

def regression_poly(x, y, d):  # polynomial regression
    xtx = vandermonde_xtx(x, d)
    xty = [sum(x[i]**j * y[i] for i in range(len(x))) for j in range(d+1)]
    return solve_chol(xtx, xty)

I get
beta_chol: [-0.1224754, -5.102945, -150.4902, 296.6456, -132.5333]
which leads to mean_sq_err_chol: 0.01398384

With your parameters and

def horner(poly, x): # evaluates poly at x; returns poly[0]+poly[1]*x+..+poly[n-1]*x^(n-1)
    y = 0
    for c in reversed(poly):
        y = y * x + c
    return y

def mean_sq_diff(x, y):
    return sum((x[i]-y[i])**2 for i in range(len(x))) / len(x)

beta = [0.210950563440863, -22.8906621119621, -20.8906997505424, 67.6304434918384, -14.1119736691275]
print('poly:', beta, 'mse:', mean_sq_diff(y, [horner(poly, t) for t in x]))

it prints:
poly: [0.2109505, -22.89066, -20.8907, 67.63043, -14.11197] mse: 0.000780353

which is a better approximation. Impressing.

Now I'm wondering what you did. How works your OSL method in detail?

[Chris-24MAY60]

Hi Raul, Being in the very early stages, I've been using Libre Office Calc (Ubuntu) because I use the graphing of the generated polynomial to ensure that the (resulting) graphed curve "looks like reality". The motion I am studying always follows a particular path. I have learned that I absolutely cannot sacrifice my six decimal place precision! For example, when I use the full precision of 6 decimal places, I get a curve that beautifully describes the movement just made. The "fit" is smooth, with but one local minimum, again, matching the movement just made. But when I "throw away" precision, say, I use only two decimal places, the "fit" no longer matches the movement, presenting instead "solutions" with both a minimum *and* a maximum within the interval. Those "solutions" are of no value. Also, not only do the magnitudes of the coefficients change dramatically with changing precision, but so do which ones are positive and which ones are negative!! Additionally, when I maintain my six decimal places of precision, I see I can ignore half of the observations - as long as I am careful to choose the first two, the fifth and sixth, the tenth, and finally, the fourteenth. It turns out that the movement is so fluid that the other seven points are merely "redundant information" (i.e., they are almost exactly on the curve) and-so those other points do not further inform the solution. The difference between my six-observation solution and my thirteen observation ...is negligible if I will make my computations with all the precision of the observations. So cutting the number of observations roughly in half should speed up and simplify the computing of (good!) coefficients - but it takes a trained eye to determine which observations may safely be thrown out - and which are salient to the solution - and that implies a great deal of contorted logic running in my device, making unmonitored decisions! This is why I am looking for a method that always permits me to use all my observations - but still runs efficiently. (fast) I am not (yet) ruling-out what you provided to me yesterday. I haven't tried it yet using all the precision my observations deliver. You've done a lot for me and I so appreciate it. Give me a day or two to play around with what you've given me. Many times, very carefully-collected, high precision, observations "save the day". (i.e., no root-taking of negative numbers and so-on - although I lack your prowess in math and cannot tell you why) I will share with you what I discover - and if further work is needed, maybe my investigations can help shape that thinking. (then) Stay tuned! Chris

…

On Thu, Jan 25, 2024, 3:05 AM Raul Kompaß ***@***.***> wrote: O.K.: In the Wikipdia <https://en.wikipedia.org/wiki/Polynomial_regression> article it is stated that the Polynomial coefficients, represented as a vector may be computed by the formula there (the last formula). _beta = (VtV)⁻1 * Vty. Where V is the Vandermonde determinant. And VtV the product of the transpose of V with V. I followed the article to call it XtX, this may be confusing. Now, to compute VtV it is not necessary to have V, you may shorten the calculation and don't need a matrix product. def vandermonde(x, d): # up to x^d dimension will be (len(x), d+1) return [[e**j for j in range(d+1)] for e in x] def vandermonde_VtV(x, d): # dimension will be (d+1, d+1) X = zeros(d+1, d+1) for j in range(d+1): X[0][j] = sum(el**j for el in x) for i in range(1, d+1): X[i][:-1] = X[i-1][1:] X[i][d] = sum(el**(i+d) for el in x) return X The equality is demonstrated in the above codes results. In order to compute (VtV)⁻1 * Vty (wich is 5x5 matrix times 5-vector) you need the inverse of VtV. Instead of computing the inverse of VtV you may but shorten again and solve VtV * _beta = Vty This solves 1 vector, while computing the inverse has to solve for n vectors. My code for this computation is here (hopefully complete). def lu_decomp(A): # https://rosettacode.org/wiki/LU_decomposition#Python """Decomposes a square matrix A by perm_mat(P) A = L U and returns L, U and P.""" n = len(A) L, U, P = eye(n), zeros(n, n), [i for i in range(n)] for j in range(n): row = max(range(j, n), key=lambda i: abs(A[i][j])) if j != row: P[j], P[row] = P[row], P[j] A2 = [A[i] for i in P] # A2 = matmul(perm_mat(P), A) for j in range(n): for i in range(j+1): U[i][j] = A2[i][j] - sum(U[k][j]*L[i][k] for k in range(i)) for i in range(j, n): L[i][j] = (A2[i][j] - sum(U[k][j]*L[i][k] for k in range(j))) / U[j][j] return L, U, P def forw_sub(L, b): # Solve L*y = b. L is assumed to be lower triangular matrix d = [0]*len(b) d[0] = b[0]/L[0][0] for i in range(1, len(b)): d[i] = (b[i] - sum(d[j] * L[i][j] for j in range(i))) / L[i][i] return d def backw_sub(U, y): # Solve U*x = y. U is assumed to be upper triangular matrix n = len(y) x = [0]*n x[-1] = y[-1] / U[-1][-1] for i in range(n-2, -1, -1): x[i] = (y[i] - sum(x[j] * U[i][j] for j in range(i+1, n))) / U[i][i] return x def lu_solve(L, U, b): # solve L*U*x = b with lower and upper triangular matrices L and U y = forw_sub(L, b) return backw_sub(U, y) def solve(A, b): # solve A*x = b for x, A is arbitrary (hopefully well conditioned) square matrix L, U, P = lu_decomp(A) pb = [b[j] for j in P] return lu_solve(L, U, pb) def solve_chol(A, b): # Alternative to solve() if A is symmetric and positive definite L = cholesky(A) return lu_solve(L, transp(L), b) solve and solve_chol differ in that instead of using the L and U matrices of an LU decomposition for the partial solutions using forward and backward substitution (which is quite equivalent to Gaussian algorithm) the Cholesky matrix and its transpose are used. All this is Ordinary Least Squares, of course, but I'm not familiar now, how the Wikipedia formula derives from the OLS principle (may look that up). But you might of course also do a gradient descent of the polynomial parameters. Did you do that? Should be also nice. With the combined def regression_poly(x, y, d): # polynomial regression xtx = vandermonde_xtx(x, d) xty = [sum(x[i]**j * y[i] for i in range(len(x))) for j in range(d+1)] return solve_chol(xtx, xty) I get beta_chol: [-0.1224754, -5.102945, -150.4902, 296.6456, -132.5333] which leads to mean_sq_err_chol: 0.01398384 With your parameters and def horner(poly, x): # evaluates poly at x; returns poly[0]+poly[1]*x+..+poly[n-1]*x^(n-1) y = 0 for c in reversed(poly): y = y * x + c return y def mean_sq_diff(x, y): return sum((x[i]-y[i])**2 for i in range(len(x))) / len(x) beta = [0.210950563440863, -22.8906621119621, -20.8906997505424, 67.6304434918384, -14.1119736691275]print('poly:', beta, 'mse:', mean_sq_diff(y, [horner(poly, t) for t in x])) it prints: poly: [0.2109505, -22.89066, -20.8907, 67.63043, -14.11197] mse: 0.000780353 which is a better approximation. Impressing. Now I'm wondering what you did. How works your OSL method in detail? — Reply to this email directly, view it on GitHub <#13479 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHSYGXESM7NEPFOUSPTLYQIN7HAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DENBSGUZDE> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[rkompass]

I like improving things and like to do maths. So you are welcome, I had a lot of fun yesterday.
That your polynomial coefficients are so different and still better than those resulting from the computations I did (and presented above) is a riddle (solved below..).
On high-precision ordinary Python the polynomial coefficients obtained from both the LU and Cholesky variant are almost exactly the same. So the choice of Cholesky factorization is not essential. It seems to allow for a bit more precise computations in the 32bit FP scenario though. Also the experts tell it to be twice as fast.
Mean square error is (depending on polynomial order), everything computed with Python 3.11 double precision:

d:2 mse: 0.606662 poly:1.4465-44.8763*x+50.8236*x^2
d:3 mse: 0.070409 poly:0.2085-24.2705*x-8.1473*x^2+42.2442*x^3
d:4 mse: 0.013954 poly:-0.1148-5.5401*x-147.2678*x^2+290.9052*x^3-129.5501*x^4
d:5 mse: 0.001040 poly:-0.0167-14.0431*x-57.5786*x^2-15.6551*x^3+270.1908*x^4-175.7435*x^5
d:6 mse: 0.000302 poly:-0.0025-19.3595*x+41.0279*x^2-595.4462*x^3+1572.3122*x^4-1422.2920*x^5+431.4606*x^6
Chris' poly parameters:
d:4 mse: 0.075886 poly:0.2110-22.8907*x-20.8907*x^2+67.6304*x^3-14.1120*x^4

Update: I did a mistake:
So you have a 4-degree poly fit that is not as good as the full computation, as is expected from the strategy of selecting only partial data.
I had a typo-mistake in the above code horner(poly) should have been horner(beta) to take your parameters.
How did you do the fitting of coefficients? (What method?)

Anyway, from my surprise (at you apparently better fit):
The x values are very close to each other. Perhaps the VtV matrix is very ill-conditioned, which may make the results unreliable despite double precision.
I'm now interested in assembling a better set of matrix computation routines and the computation of the condition of a matrix is on the todo list.

[Chris-24MAY60]

Hi Raul, I made a mistake, too! AND, I learned that in order to get uniform results, I *did* need to include my (phony) (0, 0) starting-point, else the 4th order polynomial would not cross the y = 0 line. That, it turns out, was very important. Do you have an e-mail address to which I can send my workbook? Although done in Libre Office Calc, I have saved it in Excel format. Please let me know. Chris

…

On Thu, Jan 25, 2024 at 5:55 AM Raul Kompaß ***@***.***> wrote: I like improving things and like to do maths. So you are welcome, I had a lot of fun yesterday. That your polynomial coefficients are so different and still better than those resulting from the computations I did (and presented above) is a riddle (solved below..). On high-precision ordinary Python the polynomial coefficients obtained from both the LU and Cholesky variant are almost exactly the same. So the choice of Cholesky factorization is not essential. It seems to allow for a bit more precise computations in the 32bit FP scenario though. Also the experts tell it to be twice as fast. Mean square error is (depending on polynomial order), everything computed with Python 3.11 double precision: d:2 mse: 0.606662 poly:1.4465-44.8763*x+50.8236*x^2 d:3 mse: 0.070409 poly:0.2085-24.2705*x-8.1473*x^2+42.2442*x^3 d:4 mse: 0.013954 poly:-0.1148-5.5401*x-147.2678*x^2+290.9052*x^3-129.5501*x^4 d:5 mse: 0.001040 poly:-0.0167-14.0431*x-57.5786*x^2-15.6551*x^3+270.1908*x^4-175.7435*x^5 d:6 mse: 0.000302 poly:-0.0025-19.3595*x+41.0279*x^2-595.4462*x^3+1572.3122*x^4-1422.2920*x^5+431.4606*x^6 Chris' poly parameters: d:4 mse: 0.075886 poly:0.2110-22.8907*x-20.8907*x^2+67.6304*x^3-14.1120*x^4 Update: I did a mistake: So you have a 4-degree poly fit that is not as good as the full computation, as is expected from the strategy of selecting only partial data. I had a typo-mistake in the above code horner(poly) should have been horner(beta) to take your parameters. How did you do the fitting of coefficients? (What method?) Anyway, from my surprise (at you apparently better fit): The x values are very close to each other. Perhaps the VtV matrix is very ill-conditioned, which may make the results unreliable despite double precision. I'm now interested in assembling a better set of matrix computation routines and the computation of the condition of a matrix is on the todo list. — Reply to this email directly, view it on GitHub <#13479 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHS5S75SEIWJHV52NF5LYQJB4PAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DENBUGI2TI> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[rkompass]

You may zip your workbook it and append it to the discussion here (just drop it onto the message).

[rkompass]

With

import numpy as np

def pppoly(p):
    print('{:.4f}'.format(p[0]), end='')
    for i in range(1, len(p)):
        if i == 1:
            print('{:+.4f}*x'.format(p[1]), end='')
        else:
            print('{:+.4f}*x^{}'.format(p[i], i), end='')
    print()

x = np.array([0, 0.065778, 0.121406, 0.165108, 0.2023745, 0.7878015, 0.8084525,
     0.830699, 0.8519465, 0.873236, 0.8948235, 0.917743, 0.9424455])
y = np.array([0, -1.25, -2.5, -3.75, -5, -3.75, -2.5, -1.25, 0, 1.25, 2.5, 3.75, 5])

x_plt = np.arange(0., 1., 0.05)

import matplotlib.pyplot as plt
plt.plot(x, y, 'b*')
for d in range(1, 5):
    poly = np.polyfit(x,y,deg=d)
    y_reg = np.polyval(np.poly1d(poly), x)
    mse = np.mean((np.array(y)-y_reg)**2)
    print('d:{} mse: {:5f} poly:'.format(d, mse), end=''); pppoly(poly[::-1])
    plt.plot(x_plt, np.polyval(np.poly1d(poly), x_plt))

plt.legend(['Raw data', 'Poly d=1', 'Poly d=2', 'Poly d=3', 'Poly d=4'])
plt.show()

and the computed polynomials and fits are exactly like those obtained before.

You notice that important data points in the X range 0.2 .. 0.8 are missing and so trying to get the best fit doesn't make much sense. The fits differ considerably in their predictions for that range. Except in the situation that this x range is not relevant whatsoever. Which I cannot imagine.

[Chris-24MAY60]

Au contraire! I have collections of observations (of movements) with your mentioned 0.6 seconds characteristically gone "missing" and the prediction of the extent of the displacement of the movement my device makes is better than 3/16ths of an inch - for a 12" movement!! (~ +/- 0.8%) Validated, using a video camera... I have spent a year developing a "fits-in-your-pocket" device that very accurately describes an entire movement, using relatively few, very high resolution (position vs time) "observers". Neat, huh? Chris

…

On Thu, Jan 25, 2024, 9:52 AM Raul Kompaß ***@***.***> wrote: PolyFit.png (view on web) <https://github.com/micropython/micropython/assets/90282516/843e8897-fd29-4279-8f71-d9bcc45249da> With import numpy as np def pppoly(p): print('{:.4f}'.format(p[0]), end='') for i in range(1, len(p)): if i == 1: print('{:+.4f}*x'.format(p[1]), end='') else: print('{:+.4f}*x^{}'.format(p[i], i), end='') print() x = np.array([0, 0.065778, 0.121406, 0.165108, 0.2023745, 0.7878015, 0.8084525, 0.830699, 0.8519465, 0.873236, 0.8948235, 0.917743, 0.9424455])y = np.array([0, -1.25, -2.5, -3.75, -5, -3.75, -2.5, -1.25, 0, 1.25, 2.5, 3.75, 5]) x_plt = np.arange(0., 1., 0.05) import matplotlib.pyplot as pltplt.plot(x, y, 'b*')for d in range(1, 5): poly = np.polyfit(x,y,deg=d) y_reg = np.polyval(np.poly1d(poly), x) mse = np.mean((np.array(y)-y_reg)**2) print('d:{} mse: {:5f} poly:'.format(d, mse), end=''); pppoly(poly[::-1]) plt.plot(x_plt, np.polyval(np.poly1d(poly), x_plt)) plt.legend(['Raw data', 'Poly d=1', 'Poly d=2', 'Poly d=3', 'Poly d=4'])plt.show() and the computed polynomials and fits are exactly like those obtained before. You notice that important data points in the X range 0.2 .. 0.8 are missing and so trying to get the best fit doesn't make much sense. The fits differ considerably in their predictions for that range. Except in the situation that this x range is not relevant whatsoever. Which I cannot imagine. — Reply to this email directly, view it on GitHub <#13479 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHS6PMFXGUWHFVAX7LA3YQJ5UNAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DENBWG42DC> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[rkompass]

I don't understand what you are doing. How have the measurements with x between 0.2 and 0.8 gone "missing"?
If the predictions (of y??) are already better than needed, why do you need a fit?
Why do you need a fit with a formula (polynomial) that potentially is wrong in the middle range (0.2 .. 0.8 s) (see the different curves in the above graph) where you have exact data already?

Is x the time in seconds? Is y the displacement? Of what?
What are the observers? Is this a secret?

If you are interested: I can give an explanation now why the formula VtV * _beta = Vty returns the coefficients of the least-squares fit of the polynomial to the data.

[Chris-24MAY60]

What values did you get for the coeeficients? And can this be done using MicroPython on an ESP32S3?

…

On Thu, Jan 25, 2024, 9:52 AM Raul Kompaß ***@***.***> wrote: PolyFit.png (view on web) <https://github.com/micropython/micropython/assets/90282516/843e8897-fd29-4279-8f71-d9bcc45249da> With import numpy as np def pppoly(p): print('{:.4f}'.format(p[0]), end='') for i in range(1, len(p)): if i == 1: print('{:+.4f}*x'.format(p[1]), end='') else: print('{:+.4f}*x^{}'.format(p[i], i), end='') print() x = np.array([0, 0.065778, 0.121406, 0.165108, 0.2023745, 0.7878015, 0.8084525, 0.830699, 0.8519465, 0.873236, 0.8948235, 0.917743, 0.9424455])y = np.array([0, -1.25, -2.5, -3.75, -5, -3.75, -2.5, -1.25, 0, 1.25, 2.5, 3.75, 5]) x_plt = np.arange(0., 1., 0.05) import matplotlib.pyplot as pltplt.plot(x, y, 'b*')for d in range(1, 5): poly = np.polyfit(x,y,deg=d) y_reg = np.polyval(np.poly1d(poly), x) mse = np.mean((np.array(y)-y_reg)**2) print('d:{} mse: {:5f} poly:'.format(d, mse), end=''); pppoly(poly[::-1]) plt.plot(x_plt, np.polyval(np.poly1d(poly), x_plt)) plt.legend(['Raw data', 'Poly d=1', 'Poly d=2', 'Poly d=3', 'Poly d=4'])plt.show() and the computed polynomials and fits are exactly like those obtained before. You notice that important data points in the X range 0.2 .. 0.8 are missing and so trying to get the best fit doesn't make much sense. The fits differ considerably in their predictions for that range. Except in the situation that this x range is not relevant whatsoever. Which I cannot imagine. — Reply to this email directly, view it on GitHub <#13479 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHS6PMFXGUWHFVAX7LA3YQJ5UNAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DENBWG42DC> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[rkompass]

Yes. Absolutely and without going to the limits. All you need is there and the code is posted above already.
You have to copy together the needed functions. If I forgot one, tell me.

def regression_poly(x, y, d):  # polynomial regression
    xtx = vandermonde_xtx(x, d)
    xty = [sum(x[i]**j * y[i] for i in range(len(x))) for j in range(d+1)]
    return solve_chol(xtx, xty)

returns the polynomial coefficients, i.e. the _beta in a list. pppoly prints a nice representation of the polynomial like below.

for d in range(2, 7):
    poly = regression_poly(x, y, d)
    y_reg = [horner(poly, t) for t in x]
    mse = mean_sq_diff(y, y_reg)
    print('d:{} mse: {:5f} poly:'.format(d, mse), end=''); pppoly(poly)

gives the results for a range of polynomial degrees.
With micropython on ItsyBitsy nrf52840 board (32bit Floating Point) I got this result.

d:2 mse: 0.606662 poly:1.4464-44.8757*x+50.8230*x^2
d:3 mse: 0.070409 poly:0.2084-24.2686*x-8.1526*x^2+42.2479*x^3
d:4 mse: 0.013984 poly:-0.1225-5.1029*x-150.4902*x^2+296.6456*x^3-132.5333*x^4
d:5 mse: 0.001058 poly:-0.0113-14.4378*x-53.7463*x^2-27.6074*x^3+284.8140*x^4-181.8947*x^5
d:6 mse: 0.000780 poly:0.0073-17.3704*x-8.0583*x^2-270.1364*x^3+800.4865*x^4-657.6841*x^5+160.0000*x^6

which is precise for polynomial degree up to d=5. For d=6 the fit is worse than with ordinary Python, due to limited floating point resolution.
On ESP32S3 you should get the same results.

Still they are not robust without the x-values in the middle range (around 0.3-0.7 seconds), at least that is my impression.

P.S.: See the other remarks posted simultaneously.

[Chris-24MAY60]

Thanks, Raul, Please give me a few days to understand this stuff. I will share my result with you as soon as I have worked through it. Chris

…

On Thu, Jan 25, 2024, 1:15 PM Raul Kompaß ***@***.***> wrote: I don't understand what you are doing. How have the measurements with x between 0.2 and 0.8 gone "missing"? If the predictions (of y??) are already better than needed, why do you need a fit? Why do you need a fit with a formula (polynomial) that potentially is wrong in the middle range (0.2 .. 0.8 s) (see the different curves in the above graph) where you have exact data already? Is x the time in seconds? Is y the displacement? Of what? What are the observers? Is this a secret? If you are interested: I can give an explanation now why the formula VtV * _beta = Vty returns the coefficients of the least-squares fit of the polynomial to the data. — Reply to this email directly, view it on GitHub <#13479 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHS4Z4HU4QAZ7VVUATGTYQKVNHAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DENBZGA3TQ> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[Chris-24MAY60]

Hi Raul, Say, I just got to trying the code (above) and it appears I don't have the "math" module used in your 'from math import sqrt' statement. Can you please tell me where I can find the math library module - and where it needs be installed? I am running Ubuntu. Also, you'd asked whether my "observations" were a secret. They aren't - they're simply sensors arranged to capture the "fly-by" time of an object. If I can get this running, I'll take next-steps to deliver the fourth-order coefficients Thanks, Chris

…

On Wed, Jan 24, 2024 at 11:17 AM Raul Kompaß ***@***.***> wrote: Am I to understand that the Cholesky decomposition will operate on the "product-matrix", being only five rows by one column - until out falls two matrices, one of which so-happens to be the (5x5) (X^T*X)^-1 TOGETHER WITH the (5x1) Beta matrix! Almost right: the (X^T*X) is a 5x5 matrix. from math import sqrt def pprint(A): m, n = len(A), len(A[0]) for i, r in enumerate(A): print(' [' if i else '[[', end='') for j, e in enumerate(r): print('{:6.2f}{}'.format(e, ',' if j < n-1 else ''), end='') print('],' if i < m-1 else ']]') def zeros(m, n): return [[0 for i in range(n)] for j in range(m)] def matmul(A, B): # https://rosettacode.org/wiki/Matrix_multiplication#Python return [[sum(x * B[i][col] for i,x in enumerate(row)) for col in range(len(B[0]))] for row in A] def transpose(A): return list(zip(*A)) def cholesky(A): # https://rosettacode.org/wiki/Cholesky_decomposition#Python L = zeros(len(A), len(A)) for i, (Ai, Li) in enumerate(zip(A, L)): for j, Lj in enumerate(L[:i+1]): s = sum(Li[k]*Lj[k] for k in range(j)) Li[j] = sqrt(Ai[i]-s) if (i==j) else (Ai[j]-s)/Lj[j] return L x = [0, 0.065778, 0.121406, 0.165108, 0.2023745, 0.7878015, 0.8084525, 0.830699, 0.8519465, 0.873236, 0.8948235, 0.917743, 0.9424455]y = [0, -1.25, -2.5, -3.75, -5, -3.75, -2.5, -1.25, 0, 1.25, 2.5, 3.75, 5] n = 4 # polynomial degree def vandermonde(x, d): return [[e**j for j in range(d)] for e in x] def vandermonde_xtx(x, d): # d will be 5, it is the dimension of this matrix, i.e. n+1 X = zeros(d, d) for j in range(d): X[0][j] = sum(el**j for el in x) for i in range(1, d): X[i][:-1] = X[i-1][1:] X[i][d-1] = sum(el**(i+d-1) for el in x) return X v = vandermonde(x, n+1)pprint(v); print()vtv = matmul(transpose(v), v)pprint(vtv); print() xtx = vandermonde_xtx(x, n+1)xty = [sum(x[i]**j * y[i] for i in range(len(x))) for j in range(n+1)] pprint(xtx); print()print(xty); print()l = cholesky(xtx)pprint(l); print() returns here: [[ 1.00, 0.00, 0.00, 0.00, 0.00], [ 1.00, 0.07, 0.00, 0.00, 0.00], [ 1.00, 0.12, 0.01, 0.00, 0.00], [ 1.00, 0.17, 0.03, 0.00, 0.00], [ 1.00, 0.20, 0.04, 0.01, 0.00], [ 1.00, 0.79, 0.62, 0.49, 0.39], [ 1.00, 0.81, 0.65, 0.53, 0.43], [ 1.00, 0.83, 0.69, 0.57, 0.48], [ 1.00, 0.85, 0.73, 0.62, 0.53], [ 1.00, 0.87, 0.76, 0.67, 0.58], [ 1.00, 0.89, 0.80, 0.72, 0.64], [ 1.00, 0.92, 0.84, 0.77, 0.71], [ 1.00, 0.94, 0.89, 0.84, 0.79]] [[ 13.00, 7.46, 6.07, 5.22, 4.54], [ 7.46, 6.07, 5.22, 4.54, 3.97], [ 6.07, 5.22, 4.54, 3.97, 3.48], [ 5.22, 4.54, 3.97, 3.48, 3.07], [ 4.54, 3.97, 3.48, 3.07, 2.71]] [[ 13.00, 7.46, 6.07, 5.22, 4.54], [ 7.46, 6.07, 5.22, 4.54, 3.97], [ 6.07, 5.22, 4.54, 3.97, 3.48], [ 5.22, 4.54, 3.97, 3.48, 3.07], [ 4.54, 3.97, 3.48, 3.07, 2.71]] [-7.5, 3.451842, 5.381219, 5.773444, 5.815036] [[ 3.61, 0.00, 0.00, 0.00, 0.00], [ 2.07, 1.34, 0.00, 0.00, 0.00], [ 1.68, 1.29, 0.16, 0.00, 0.00], [ 1.45, 1.16, 0.23, 0.06, 0.00], [ 1.26, 1.02, 0.26, 0.12, 0.01]] Now what's still needed is forward substitution (and perhaps backward substitution, but probably not). Which I have found already but have to check/improve yet. Then the solution of xty = xtx * beta is there. A lot of cross checking will be necessary too. I included computing the raw vandermonde matrix and multiplying it with it's transpose. You see the result is the same. Doing LU decomposition and solving the equation by that as an alternative - to be done yet. But you see: There is not much to compute, once it is clear, how it should be done. — Reply to this email directly, view it on GitHub <#13479 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHSYSZZ6NY5SRRRVUAWTYQE6ZBAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMZVGY3TK> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[Chris-24MAY60]

Hello again, Raul, When I try running the script above, I get this: from: can't read /var/mail/math ./Cholesky_Decomposition.py: line 3: syntax error near unexpected token `(' ./Cholesky_Decomposition.py: line 3: `def pprint(A):' Chris

…

On Wed, Jan 24, 2024 at 11:17 AM Raul Kompaß ***@***.***> wrote: Am I to understand that the Cholesky decomposition will operate on the "product-matrix", being only five rows by one column - until out falls two matrices, one of which so-happens to be the (5x5) (X^T*X)^-1 TOGETHER WITH the (5x1) Beta matrix! Almost right: the (X^T*X) is a 5x5 matrix. from math import sqrt def pprint(A): m, n = len(A), len(A[0]) for i, r in enumerate(A): print(' [' if i else '[[', end='') for j, e in enumerate(r): print('{:6.2f}{}'.format(e, ',' if j < n-1 else ''), end='') print('],' if i < m-1 else ']]') def zeros(m, n): return [[0 for i in range(n)] for j in range(m)] def matmul(A, B): # https://rosettacode.org/wiki/Matrix_multiplication#Python return [[sum(x * B[i][col] for i,x in enumerate(row)) for col in range(len(B[0]))] for row in A] def transpose(A): return list(zip(*A)) def cholesky(A): # https://rosettacode.org/wiki/Cholesky_decomposition#Python L = zeros(len(A), len(A)) for i, (Ai, Li) in enumerate(zip(A, L)): for j, Lj in enumerate(L[:i+1]): s = sum(Li[k]*Lj[k] for k in range(j)) Li[j] = sqrt(Ai[i]-s) if (i==j) else (Ai[j]-s)/Lj[j] return L x = [0, 0.065778, 0.121406, 0.165108, 0.2023745, 0.7878015, 0.8084525, 0.830699, 0.8519465, 0.873236, 0.8948235, 0.917743, 0.9424455]y = [0, -1.25, -2.5, -3.75, -5, -3.75, -2.5, -1.25, 0, 1.25, 2.5, 3.75, 5] n = 4 # polynomial degree def vandermonde(x, d): return [[e**j for j in range(d)] for e in x] def vandermonde_xtx(x, d): # d will be 5, it is the dimension of this matrix, i.e. n+1 X = zeros(d, d) for j in range(d): X[0][j] = sum(el**j for el in x) for i in range(1, d): X[i][:-1] = X[i-1][1:] X[i][d-1] = sum(el**(i+d-1) for el in x) return X v = vandermonde(x, n+1)pprint(v); print()vtv = matmul(transpose(v), v)pprint(vtv); print() xtx = vandermonde_xtx(x, n+1)xty = [sum(x[i]**j * y[i] for i in range(len(x))) for j in range(n+1)] pprint(xtx); print()print(xty); print()l = cholesky(xtx)pprint(l); print() returns here: [[ 1.00, 0.00, 0.00, 0.00, 0.00], [ 1.00, 0.07, 0.00, 0.00, 0.00], [ 1.00, 0.12, 0.01, 0.00, 0.00], [ 1.00, 0.17, 0.03, 0.00, 0.00], [ 1.00, 0.20, 0.04, 0.01, 0.00], [ 1.00, 0.79, 0.62, 0.49, 0.39], [ 1.00, 0.81, 0.65, 0.53, 0.43], [ 1.00, 0.83, 0.69, 0.57, 0.48], [ 1.00, 0.85, 0.73, 0.62, 0.53], [ 1.00, 0.87, 0.76, 0.67, 0.58], [ 1.00, 0.89, 0.80, 0.72, 0.64], [ 1.00, 0.92, 0.84, 0.77, 0.71], [ 1.00, 0.94, 0.89, 0.84, 0.79]] [[ 13.00, 7.46, 6.07, 5.22, 4.54], [ 7.46, 6.07, 5.22, 4.54, 3.97], [ 6.07, 5.22, 4.54, 3.97, 3.48], [ 5.22, 4.54, 3.97, 3.48, 3.07], [ 4.54, 3.97, 3.48, 3.07, 2.71]] [[ 13.00, 7.46, 6.07, 5.22, 4.54], [ 7.46, 6.07, 5.22, 4.54, 3.97], [ 6.07, 5.22, 4.54, 3.97, 3.48], [ 5.22, 4.54, 3.97, 3.48, 3.07], [ 4.54, 3.97, 3.48, 3.07, 2.71]] [-7.5, 3.451842, 5.381219, 5.773444, 5.815036] [[ 3.61, 0.00, 0.00, 0.00, 0.00], [ 2.07, 1.34, 0.00, 0.00, 0.00], [ 1.68, 1.29, 0.16, 0.00, 0.00], [ 1.45, 1.16, 0.23, 0.06, 0.00], [ 1.26, 1.02, 0.26, 0.12, 0.01]] Now what's still needed is forward substitution (and perhaps backward substitution, but probably not). Which I have found already but have to check/improve yet. Then the solution of xty = xtx * beta is there. A lot of cross checking will be necessary too. I included computing the raw vandermonde matrix and multiplying it with it's transpose. You see the result is the same. Doing LU decomposition and solving the equation by that as an alternative - to be done yet. But you see: There is not much to compute, once it is clear, how it should be done. — Reply to this email directly, view it on GitHub <#13479 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHSYSZZ6NY5SRRRVUAWTYQE6ZBAVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMZVGY3TK> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[rkompass]

Mhm, how do you run simple python scripts? I feel like a beginner and use Thonny. But there are several possibilities.
Thonny allows me to see the script and click a green triangle to run it. I may switch from the attached board (ItsyBitsy at the moment) to locally installed Micropython to ordinary Python and back again.
The results are almost the same.

I am running Ubuntu.

I'm running Manjaro, but this shouldn't make a difference.
I put my script with all functions and the approximation in it here.
matrix.zip

Just unzip it and load it into Thonny.
Or run it from the terminal/shell with python matrix.py. You then see the results.
Disclaimer: It is not finished and contains errors (QR decomposition has flaws yes, for example).

[Chris-24MAY60]

Hi Raul, I've convinced myself that the ESP32S3 is not the place to perform my computations. That's fine, I just needed to know. Means that while my little sensor does an adequate job of wirelessly sending great results, with the degree of precision I need, I will need an app, running on a mobile device to produce the level of results interpretable to the user. Thank you, so much, for your help. I learned a lot from you. Every time I learn from smart people, I wish I'd learned more in their area. Too soon we grow old - and too late we grow smart. Chris Chris

…

On Fri, Jan 26, 2024, 4:28 PM Raul Kompaß ***@***.***> wrote: Mhm, how do you run simple python scripts? I feel like a beginner and use Thonny. But there are several possibilities. Thonny allows me to see the script and click a green triangle to run it. I may switch from the attached board (ItsyBitsy at the moment) to locally installed Micropython to ordinary Python and back again. The results are almost the same. I am running Ubuntu. I'm running Manjaro, but this shouldn't make a difference. I put my script with all functions and the approximation in it here. matrix.zip <https://github.com/micropython/micropython/files/14069811/matrix.zip> Just unzip it and load it into Thonny. Or run it from the terminal/shell with python matrix.py. You then see the results. Disclaimer: It is not finished and contains errors (QR decomposition has flaws yes, for example). — Reply to this email directly, view it on GitHub <#13479 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AF7CHS435X223JSM5ZVDKGTYQQUZ5AVCNFSM6AAAAABCCZ63IWVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DENRSGEZTO> . You are receiving this because you authored the thread.Message ID: ***@***.*** com>

[rkompass]

You are welcome, of course. Feel free to ask for an explanation on the algebra, once you are interested in it. I understood it once, forgot it many years ago and now I got again a grasp of it.
My feeling is: We older people can still learn a lot, but it takes much longer and we forget much faster...
On the other hand: When I saw my pupils then being a teacher: They forgot almost everything.